clear;clc;close all

%%% Matlab Climate Data Toolbox
%%% (https://www.mathworks.com/matlabcentral/fileexchange/70338-climate-data-toolbox-for-matlab)
%%% is prerequisite to determine land-sea mask

lon=0:2:358;;lat=-90:2:90;
NIOlon=find(lon>=40 & lon<=100);NIOlat=find(lat>=0 & lat<=30);
WNPlon=find(lon>=100 & lon<=180);WNPlat=find(lat>=0 & lat<=50);
ENPlon=find(lon>=180 & lon<=260);ENPlat=find(lat>=0 & lat<=50);
NATlon=find(lon>=260 & lon<=360);NATlat=find(lat>=0 & lat<=50);
SIOlon=find(lon>=30 & lon<=135);SIOlat=find(lat>=-50 & lat<=0);
SPOlon=find(lon>=135 & lon<=295);SPOlat=find(lat>=-50 & lat<=0);
SATlon=find(lon>=295 & lon<=360);SATlat=find(lat>=-50 & lat<=0);
[basin,X,Y]=basin_CNN(lon,lat);xland=island(Y,X);
% basinID: 1:NIO; 2:WNP; 3:ENP; 4:NAT; 5:SIO; 6: SPO; 7:SAT
basinID=4;

%%% Load large-scale CESM data for 1877-2000 and preprocess
load('demo.mat')
yrlen=size(SSTcesm,3);
%%% Crop to selected basin, set land data to 0
for i=1:yrlen
    i
    dummy=(SSTcesm(:,:,i)-SSTmean)/SSTstd; dummy(basin~=basinID)=0;dummy(xland)=0; XTest(:,:,1,i)=dummy(NATlon,NATlat); % SST
    dummy=(SDcesm(:,:,i)-SDmean)/SDstd; dummy(basin~=basinID)=0;dummy(xland)=0; XTest(:,:,2,i)=dummy(NATlon,NATlat); % Q
    dummy=(VORTcesm(:,:,i)-VORTmean)/VORTstd; dummy(basin~=basinID)=0;dummy(xland)=0; XTest(:,:,3,i)=dummy(NATlon,NATlat); % VORT
    dummy=(VWScesm(:,:,i)-VWSmean)/VWSstd; dummy(basin~=basinID)=0;dummy(xland)=0; XTest(:,:,4,i)=dummy(NATlon,NATlat); % VWS
end

%%% Load transfer-learned CNN
load('NAT_TC_ensembleCNNs.mat')

%%% Predict seasonal TC activity
for i=1:length(Transfernet)
disp(['Ensemble # ',num2str(i),' is testing'])
net = Transfernet{i};
YTest(:,i) = predict(net,XTest);
end

%%% Note that, regression output layer can lead to negative values, which is meaningless for seasonal TC activity 
YTest=max(YTest,0);

figure;plot(year,YTest,'-k');
hold on;plot(year,nanmean(YTest,2),'-r','linewidth',3)

