reportmeans=1;
reportstds=1;
runparammap=2;
reportbestfit=1;
doPCA1=2; % 7 parameters; TA and forcings
doPCA2=2; % 7 parameters; population and forcings
doPCA3=2; % 8 parameters; TA, population and forcings
getcontribseries=2; % to get series with the contributions of each individual forcing

% load all
if loaded==2
	load('temperature_recent')
	load('temperature_prehistoric')
	
	load('population')
	load('CO2_recent')
	load('CH4')
	load('N2O')
	
	load('CO2_prehistoric')
	load('solar')
	load('strataer')
	soi=load('soi.csv');

	loaded=1;
end

% range of modelling
max([min(population(:,1)),min(CO2_recent(:,1)),min(CH4(:,1)),min(N2O(:,3))]); % anthropogenic
min([max(population(:,1)),max(CO2_recent(:,1)),max(CH4(:,1)),max(N2O(:,3))]); % anthropogenic

ANModelmin=max([min(population(:,1)),min(CO2_recent(:,1)),min(CH4(:,1)),min(N2O(:,3)),min(solar(:,1)),min(strataer(:,1)),min(soi(:,1))]); % natural and anthro
ANModelmax=min([max(population(:,1)),max(CO2_recent(:,1)),max(CH4(:,1)),max(N2O(:,3)),max(solar(:,1)),max(strataer(:,1)),max(soi(:,1))]); % natural and anthro

maxrangeNModelmin=max([min(solar(:,1)),min(strataer(:,1)),min(soi(:,1))]); % natural
maxrangeNModelmax=min([max(solar(:,1)),max(strataer(:,1)),max(soi(:,1))]); % natural

% testing the base of the anomaly (1961-90)
mean(temperature_recent(find((temperature_recent(:,1)>1960).*(temperature_recent(:,1)<1991)),2));

% find bases of all anomalies
meanpop=mean(population(find((population(:,1)>1960).*(population(:,1)<1991)),2));
meanCO2=mean(CO2_recent(find((CO2_recent(:,1)>1960).*(CO2_recent(:,1)<1991)),2));
meanCH4=mean(CH4(find((CH4(:,1)>1960).*(CH4(:,1)<1991)),2));
meanN2O=mean(N2O(find((N2O(:,3)>1960).*(N2O(:,3)<1991)),4));
meansolar=mean(solar(find((solar(:,1)>1960).*(solar(:,1)<1991)),3));
meanstrataer=mean(strataer(find((strataer(:,1)>1960).*(strataer(:,1)<1991)),3));
meansoi=mean(soi(find((soi(:,1)>1960).*(soi(:,1)<1991)),4));

if reportmeans==1
	printf(['1961-90 means','\n',...
		'population: ',num2str(meanpop),'\n',...
		'CO2: ',num2str(meanCO2),'\n',...
		'CH4: ',num2str(meanCH4),'\n',...
		'N2O: ',num2str(meanN2O),'\n',...
		'solar irradiance: ',num2str(meansolar),'\n',...
		'stratospheric aerosols: ',num2str(meanstrataer),'\n',...
		'southern oscillation index: ',num2str(meansoi),'\n',...
		'\n'...
	])
end

% anomalies for 1961-90
popanom=population(find((population(:,1)>1960).*(population(:,1)<1991)),2)-meanpop;
CO2anom=CO2_recent(find((CO2_recent(:,1)>1960).*(CO2_recent(:,1)<1991)),2)-meanCO2;
CH4anom=CH4(find((CH4(:,1)>1960).*(CH4(:,1)<1991)),2)-meanCH4;
N2Oanom=N2O(find((N2O(:,3)>1960).*(N2O(:,3)<1991)),4)-meanN2O;
solaranom=solar(find((solar(:,1)>1960).*(solar(:,1)<1991)),3)-meansolar;
strataeranom=strataer(find((strataer(:,1)>1960).*(strataer(:,1)<1991)),3)-meanstrataer;
soianom=soi(find((soi(:,1)>1960).*(soi(:,1)<1991)),4)-meansoi;

stdpopanom=std(popanom);
stdCO2anom=std(CO2anom);
stdCH4anom=std(CH4anom);
stdN2Oanom=std(N2Oanom);
stdsolaranom=std(solaranom);
stdstrataeranom=std(strataeranom);
stdsoianom=std(soianom);

if reportstds==1
	printf(['standard deviations','\n',...
		'population: ',num2str(stdpopanom),'\n',...
		'CO2: ',num2str(stdCO2anom),'\n',...
		'CH4: ',num2str(stdCH4anom),'\n',...
		'N2O: ',num2str(stdN2Oanom),'\n',...
		'solar irradiance: ',num2str(stdsolaranom),'\n',...
		'stratospheric aerosols: ',num2str(stdstrataeranom),'\n',...
		'southern oscillation index: ',num2str(stdsoianom),'\n'...
		'\n'...
	])
end

% anomalies for all
popanom=population(:,2)-meanpop;
CO2anom=CO2_recent(:,2)-meanCO2;
CH4anom=CH4(:,2)-meanCH4;
N2Oanom=N2O(:,4)-meanN2O;
solaranom=solar(:,3)-meansolar;
strataeranom=strataer(:,3)-meanstrataer;
soianom=soi(:,4)-meansoi;

% quantities over the model timespans
Nperiod=4; % 1 for 1876-2012 (maximal range), 2 for 1876-1960 (before massive GHGs), 3 for 1991-2012 (massive GHGs), 4 for 1984-2012 (maximal anthropogenic range)

if Nperiod==1
	NModelmin=maxrangeNModelmin;
	NModelmax=maxrangeNModelmax;
elseif Nperiod==2
	NModelmin=maxrangeNModelmin;
	NModelmax=1960;
elseif Nperiod==3
	NModelmin=1991;
	NModelmax=maxrangeNModelmax;
else % Nperiod==4
	NModelmin=ANModelmin;
	NModelmax=maxrangeNModelmax;
end

Ntemps=temperature_recent(find((temperature_recent(:,1)>(NModelmin-1)).*(temperature_recent(:,1)<(NModelmax+1))),2);
ANtemps=temperature_recent(find((temperature_recent(:,1)>(ANModelmin-1)).*(temperature_recent(:,1)<(ANModelmax+1))),2);

Nsolaranom=solar(find((solar(:,1)>(NModelmin-1)).*(solar(:,1)<(NModelmax+1))),3)-meansolar;
Nstrataeranom=strataer(find((strataer(:,1)>(NModelmin-1)).*(strataer(:,1)<(NModelmax+1))),3)-meanstrataer;
Nsoianom=soi(find((soi(:,1)>(NModelmin-1)).*(soi(:,1)<(NModelmax+1))),4)-meansoi;

ANpopanom=population(find((population(:,1)>(ANModelmin-1)).*(population(:,1)<(ANModelmax+1))),2)-meanpop;
ANCO2anom=CO2_recent(find((CO2_recent(:,1)>(ANModelmin-1)).*(CO2_recent(:,1)<(ANModelmax+1))),2)-meanCO2;
ANCH4anom=CH4(find((CH4(:,1)>(ANModelmin-1)).*(CH4(:,1)<(ANModelmax+1))),2)-meanCH4;
ANN2Oanom=N2O(find((N2O(:,3)>(ANModelmin-1)).*(N2O(:,3)<(ANModelmax+1))),4)-meanN2O;
ANsolaranom=solar(find((solar(:,1)>(ANModelmin-1)).*(solar(:,1)<(ANModelmax+1))),3)-meansolar;
ANstrataeranom=strataer(find((strataer(:,1)>(ANModelmin-1)).*(strataer(:,1)<(ANModelmax+1))),3)-meanstrataer;
ANsoianom=soi(find((soi(:,1)>(ANModelmin-1)).*(soi(:,1)<(ANModelmax+1))),4)-meansoi;

% annualise data for N2O, strataer, soi - doing this after the filtering by year so there are no initial-few-lines complications
ANN2Oanom=transpose(mean(reshape(ANN2Oanom,[12 numel(ANN2Oanom)/12])));
Nsoianom=transpose(mean(reshape(Nsoianom,[12 numel(Nsoianom)/12])));
ANsoianom=transpose(mean(reshape(ANsoianom,[12 numel(ANsoianom)/12])));

% but the strataer data has some months missing, so I can't reshape in the same way - but all of the missing data is in the last year, 2012
if or(Nperiod==1,Nperiod==3,Nperiod==4)
	Nstrataeranom=[transpose(mean(reshape(Nstrataeranom(1:end-9),[12 numel(Nstrataeranom(1:end-9))/12])));mean(Nstrataeranom(end-8:end))];
else
	Nstrataeranom=transpose(mean(reshape(Nstrataeranom,[12 numel(Nstrataeranom)/12])));
end
ANstrataeranom=[transpose(mean(reshape(ANstrataeranom(1:end-9),[12 numel(ANstrataeranom(1:end-9))/12])));mean(ANstrataeranom(end-8:end))];

% control structure parameters
% sampling parameters
% for each run: set sampleres, whether to record, which models to test
sampleres=2; % 1 for low (0.1), 2 for high (0.01)
record=2; % 1 to record
testNModel=1; % 1 to test
testANModel=2; % 1 to test
% check that the forcing vectors are filtered for the desired periods
% next set each of the sampling vectors below

impactsample=[-0:0.1:1]; % this has no influence unless it is used in the next few lines

if sampleres==1
	sampleinterval=0.1;
else
	sampleinterval=0.01;
end

SolImpactsample=[0:sampleinterval:0.1];
StrataerImpactsample=[0.1:sampleinterval:0.3];
ENSOImpactsample=[0:sampleinterval:0.2];
CO2Impactsample=[0:sampleinterval:1.0];
CH4Impactsample=[0:sampleinterval:1.0];
N2OImpactsample=[0:sampleinterval:1.0];

% initialise
minNModelRMSD=Inf;
minNModelRMSDparams=[NaN,NaN,NaN];

minANModelRMSD=Inf;
minANModelRMSDparams=[NaN,NaN,NaN,NaN,NaN,NaN];

if record==1
	results=[];
end

% fitting
if runparammap==1
	for SolImpact=SolImpactsample
		printf(['value ',num2str(find(SolImpactsample==SolImpact)),' of ',num2str(numel(SolImpactsample)),'\n']) % progress report
		if find(SolImpactsample==SolImpact)==numel(SolImpactsample)
			printf('\n')
		end
		for StrataerImpact=StrataerImpactsample
			for ENSOImpact=ENSOImpactsample
				if testNModel==1
					NModel=0 ...
						+SolImpact*Nsolaranom/stdsolaranom ...
						-StrataerImpact*Nstrataeranom/stdstrataeranom ...
						-ENSOImpact*Nsoianom/stdsoianom;
					NModelDev=Ntemps-NModel;
					
					%AvgNModelDev=sum(NModelDev)/length(NModelDev);
					NModelRMSD=sqrt((sum(NModelDev.^2))/length(NModelDev));
					%NModelMAE=(sum(abs(NModelDev)))/length(NModelDev);
		
					if NModelRMSD<minNModelRMSD
						minNModelRMSD=NModelRMSD;
						minNModelRMSDparams=[SolImpact,StrataerImpact,ENSOImpact];
					end
				end
				if testANModel==1
					for CO2Impact=CO2Impactsample
						for CH4Impact=CH4Impactsample
							for N2OImpact=N2OImpactsample
								ANModel=0 ...
									+SolImpact*ANsolaranom/stdsolaranom ...
									-StrataerImpact*ANstrataeranom/stdstrataeranom ...
									-ENSOImpact*ANsoianom/stdsoianom ...
									+CO2Impact*ANCO2anom/stdCO2anom ...
									+CH4Impact*ANCH4anom/stdCH4anom ...
									+N2OImpact*ANN2Oanom/stdN2Oanom;
								ANModelDev=ANtemps-ANModel;
								
								%AvgANModelDev=sum(ANModelDev)/length(ANModelDev);
								ANModelRMSD=sqrt((sum(ANModelDev.^2))/length(ANModelDev));
								%ANModelMAE=(sum(abs(ANModelDev)))/length(ANModelDev);
		
								if ANModelRMSD<minANModelRMSD
									minANModelRMSD=ANModelRMSD;
									minANModelRMSDparams=[SolImpact,StrataerImpact,ENSOImpact,CO2Impact,CH4Impact,N2OImpact];
								end
								if record==1
									results=[results;[CO2Impact,SolImpact,ANModelRMSD]];
								end
							end
						end
					end
				end
			end
		end
	end
	
	if reportbestfit==1
		reportstring=[''];
		if testNModel==1
			Nreportstring=['best natural model has:','\n',...
				'RMSD: ',num2str(minNModelRMSD),'\n',...
				'SolImpact: ',num2str(minNModelRMSDparams(1)),'\n',...
				'StrataerImpact: ',num2str(minNModelRMSDparams(2)),'\n',...
				'ENSOImpact: ',num2str(minNModelRMSDparams(3)),'\n',...
				'\n'...
			];

			reportstring=[reportstring,Nreportstring];
		end
		if testANModel==1
			ANreportstring=['best anthropogenic+natural model has:','\n',...
				'RMSD: ',num2str(minANModelRMSD),'\n',...
				'SolImpact: ',num2str(minANModelRMSDparams(1)),'\n',...
				'StrataerImpact: ',num2str(minANModelRMSDparams(2)),'\n',...
				'ENSOImpact: ',num2str(minANModelRMSDparams(3)),'\n',...
				'CO2Impact: ',num2str(minANModelRMSDparams(4)),'\n',...
				'CH4Impact: ',num2str(minANModelRMSDparams(5)),'\n',...
				'N2OImpact: ',num2str(minANModelRMSDparams(6)),'\n'...
			];

			reportstring=[reportstring,ANreportstring];
		end
		printf(reportstring)
	end
	if record==1
		save('results','results')
	end
end

% principal component analysis
% it doesn't matter whether we use the anomaly or original data, but the anomaly data is already filtered for the period 1984-2012
% I would say 'assuming that Nperiod=4 above', but it doesn't matter because we're only interested in the ANModel vectors

pkg load statistics % this requires octave-forge statistics package

% it's surprising that octave doesn't do this automatically, but we need to normalise so that each has unit variance
prepcaANpopanom=sqrt(1/var(ANpopanom))*ANpopanom;
prepcaANtemps=sqrt(1/var(ANtemps))*ANtemps;
prepcaANCO2anom=sqrt(1/var(ANCO2anom))*ANCO2anom;
prepcaANCH4anom=sqrt(1/var(ANCH4anom))*ANCH4anom;
prepcaANN2Oanom=sqrt(1/var(ANN2Oanom))*ANN2Oanom;
prepcaANsolaranom=sqrt(1/var(ANsolaranom))*ANsolaranom;
prepcaANstrataeranom=sqrt(1/var(ANstrataeranom))*ANstrataeranom;
prepcaANsoianom=sqrt(1/var(ANsoianom))*ANsoianom;

% PCA for forcings and TA
if doPCA1==1
	X=[prepcaANtemps, prepcaANCO2anom, prepcaANCH4anom, prepcaANN2Oanom, prepcaANsolaranom, prepcaANstrataeranom, prepcaANsoianom];
	%u % empirical mean
	
	[COEFF,SCORE,latent] = princomp(X)
		% coeff: the principal components (columns) built out of the original basis (rows)
		% score: the matrix X in PC space
		% eigenvalues of the covariance matrix X . . . actually I think the covariace matrix is cov(X)
	
	% project forcings onto temperature
	reducedCOEFF=COEFF(:,1:2); % most (74%) of the action, according to the latent vector, is in PC1 and PC2
	
	CO2ontemp=dot(reducedCOEFF(1,:),reducedCOEFF(2,:))/norm(reducedCOEFF(1,:))/norm(reducedCOEFF(1,:))
	CH4ontemp=dot(reducedCOEFF(1,:),reducedCOEFF(3,:))/norm(reducedCOEFF(1,:))/norm(reducedCOEFF(1,:))
	N2Oontemp=dot(reducedCOEFF(1,:),reducedCOEFF(4,:))/norm(reducedCOEFF(1,:))/norm(reducedCOEFF(1,:))
	solarontemp=dot(reducedCOEFF(1,:),reducedCOEFF(5,:))/norm(reducedCOEFF(1,:))/norm(reducedCOEFF(1,:))
	strataerontemp=dot(reducedCOEFF(1,:),reducedCOEFF(6,:))/norm(reducedCOEFF(1,:))/norm(reducedCOEFF(1,:))
	soiontemp=dot(reducedCOEFF(1,:),reducedCOEFF(7,:))/norm(reducedCOEFF(1,:))/norm(reducedCOEFF(1,:))
end

% PCA for population and forcings
if doPCA2==1
	X=[prepcaANpopanom, prepcaANCO2anom, prepcaANCH4anom, prepcaANN2Oanom, prepcaANsolaranom, prepcaANstrataeranom, prepcaANsoianom];
	
	[COEFF,SCORE,latent] = princomp(X)
		% coeff: the principal components (columns) built out of the original basis (rows)
		% score: the matrix X in PC space
		% eigenvalues of the covariance matrix X . . . actually I think the covariace matrix is cov(X)
	
	% project population onto forcings
	reducedCOEFF=COEFF(:,1:2); % most (77%) of the action, according to the latent vector, is in PC1 and PC2
	
	% regarding the direction in which these projections are done: they are an indication of covariance, and the direction doesn't say anything about causation. projecting onto the larger vector, which the population vector tends to be here,, is better because some of the others are small, leading to curiously-large results even when two vectors are almost perpendicular
	CO2onpop=dot(reducedCOEFF(1,:),reducedCOEFF(2,:))/norm(reducedCOEFF(1,:))/norm(reducedCOEFF(1,:))
	CH4onpop=dot(reducedCOEFF(1,:),reducedCOEFF(3,:))/norm(reducedCOEFF(1,:))/norm(reducedCOEFF(1,:))
	N2Oonpop=dot(reducedCOEFF(1,:),reducedCOEFF(4,:))/norm(reducedCOEFF(1,:))/norm(reducedCOEFF(1,:))
	solaronpop=dot(reducedCOEFF(1,:),reducedCOEFF(5,:))/norm(reducedCOEFF(1,:))/norm(reducedCOEFF(1,:))
	strataeronpop=dot(reducedCOEFF(1,:),reducedCOEFF(6,:))/norm(reducedCOEFF(1,:))/norm(reducedCOEFF(1,:))
	soionpop=dot(reducedCOEFF(1,:),reducedCOEFF(7,:))/norm(reducedCOEFF(1,:))/norm(reducedCOEFF(1,:))
end

% PCA for temp, population and forcings all together (realised that I don't have eniough space on the dashboard to present twoo separate PCAs)
if doPCA3==1
	X=[prepcaANtemps, prepcaANpopanom, prepcaANCO2anom, prepcaANCH4anom, prepcaANN2Oanom, prepcaANsolaranom, prepcaANstrataeranom, prepcaANsoianom];
	
	[COEFF,SCORE,latent] = princomp(X)
		% coeff: the principal components (columns) built out of the original basis (rows)
		% score: the matrix X in PC space
		% eigenvalues of the covariance matrix X . . . actually I think the covariace matrix is cov(X)
	
	% project forcings onto temperature
	reducedCOEFF=COEFF(:,1:2); % most (77%) of the action, according to the latent vector, is in PC1 and PC2
	
	popontemp=dot(reducedCOEFF(1,:),reducedCOEFF(2,:))/norm(reducedCOEFF(1,:))/norm(reducedCOEFF(1,:))
	CO2ontemp=dot(reducedCOEFF(1,:),reducedCOEFF(3,:))/norm(reducedCOEFF(1,:))/norm(reducedCOEFF(1,:))
	CH4ontemp=dot(reducedCOEFF(1,:),reducedCOEFF(4,:))/norm(reducedCOEFF(1,:))/norm(reducedCOEFF(1,:))
	N2Oontemp=dot(reducedCOEFF(1,:),reducedCOEFF(5,:))/norm(reducedCOEFF(1,:))/norm(reducedCOEFF(1,:))
	solarontemp=dot(reducedCOEFF(1,:),reducedCOEFF(6,:))/norm(reducedCOEFF(1,:))/norm(reducedCOEFF(1,:))
	strataerontemp=dot(reducedCOEFF(1,:),reducedCOEFF(7,:))/norm(reducedCOEFF(1,:))/norm(reducedCOEFF(1,:))
	soiontemp=dot(reducedCOEFF(1,:),reducedCOEFF(8,:))/norm(reducedCOEFF(1,:))/norm(reducedCOEFF(1,:))
	
	% regarding the direction in which these projections are done: they are an indication of covariance, and the direction doesn't say anything about causation. projecting onto the larger vector, which the population vector tends to be here,, is better because some of the others are small, leading to curiously-large results even when two vectors are almost perpendicular
	temponpop=dot(reducedCOEFF(1,:),reducedCOEFF(1,:))/norm(reducedCOEFF(2,:))/norm(reducedCOEFF(2,:))
	CO2onpop=dot(reducedCOEFF(1,:),reducedCOEFF(3,:))/norm(reducedCOEFF(2,:))/norm(reducedCOEFF(2,:))
	CH4onpop=dot(reducedCOEFF(1,:),reducedCOEFF(4,:))/norm(reducedCOEFF(2,:))/norm(reducedCOEFF(2,:))
	N2Oonpop=dot(reducedCOEFF(1,:),reducedCOEFF(5,:))/norm(reducedCOEFF(2,:))/norm(reducedCOEFF(2,:))
	solaronpop=dot(reducedCOEFF(1,:),reducedCOEFF(6,:))/norm(reducedCOEFF(2,:))/norm(reducedCOEFF(2,:))
	strataeronpop=dot(reducedCOEFF(1,:),reducedCOEFF(7,:))/norm(reducedCOEFF(2,:))/norm(reducedCOEFF(2,:))
	soionpop=dot(reducedCOEFF(1,:),reducedCOEFF(8,:))/norm(reducedCOEFF(2,:))/norm(reducedCOEFF(2,:))
end

% contributions of individual forcings

if getcontribseries==1
	% we do ANModel 1984-2012
	SolImpact=0.02;
	StrataerImpact=0.05;
	ENSOImpact=0.08;
	CO2Impact=0.06;
	CH4Impact=0.03;
	N2OImpact=0.01;

	numyears=maxrangeNModelmax-ANModelmin+1;

	contribdata=[...
		transpose([ANModelmin:maxrangeNModelmax]),CO2Impact*ANCO2anom/stdCO2anom,repmat(1,[numyears,1]);...
		transpose([ANModelmin:maxrangeNModelmax]),CH4Impact*ANCH4anom/stdCH4anom,repmat(2,[numyears,1]);...
		transpose([ANModelmin:maxrangeNModelmax]),N2OImpact*ANN2Oanom/stdN2Oanom,repmat(3,[numyears,1]);...
		transpose([ANModelmin:maxrangeNModelmax]),SolImpact*ANsolaranom/stdsolaranom,repmat(4,[numyears,1]);...
		transpose([ANModelmin:maxrangeNModelmax]),-StrataerImpact*ANstrataeranom/stdstrataeranom,repmat(5,[numyears,1]);...
		transpose([ANModelmin:maxrangeNModelmax]),-ENSOImpact*ANsoianom/stdsoianom,repmat(6,[numyears,1])...
	];
end