%
restoredefaultpath

clear
close all
clc

addpath('./../../tools')
addpath('./../../soda-ssp-mo')
soda -docinstall


timeObsVec = [96:1:216];
%timeObsVec = [0:1:96];
timeSimVec = timeObsVec;
 
nRows = 15;
nCols = 7;
nLays = 5;


load('./../20110930a/worker001/soildepth.mat')
load('./../20110930a/worker001/initial-pressure-head.mat')
%load('./../20110930a/worker001/preshead-0-at-lower-bound.mat')
load('./../20110930a/worker001/artificial-obs.mat')
load('./../20110930a/worker001/water-balance-components.mat')


rowSel = ismember(artifObs.times,timeObsVec);


presHead52 = presHeadIni(:,:,2:5);
presHead1 = presHeadIni(:,:,1);
presHead1Old = presHead1;

artifObs_h = artifObs.h(rowSel,:);
artifObs_Q = artifObs.Q(rowSel,1);
artifObs_times = artifObs.times(rowSel,1);

mConstants = {'artifObs_h',artifObs_h;...
              'artifObs_Q',artifObs_Q;...
              'artifObs_times',artifObs_times};

mStatesNOKFInit(1,1:nRows*nCols) = presHead1Old(:)';
mStatesNOKFInit(1,nRows*nCols+(1:nRows*nCols*(nLays-1))) = presHead52(:)';
mStatesNOKFInit(1,nRows*nCols*nLays+[1:6]) = 0; % initial 'discharge','vInit','vPrec','vSink','vSeep','vEnd'

upperstatebound = soilDepth(:);

sodaPar.parNames = {'LOG10rSink','LOG10Ks'};
sodaPar.parNamesTex = {'log_{10}(r_{sink})','log_{10}(K_s)'};
sodaPar.parSpaceLoBound = log10([1.0e-3, 0.01]);
sodaPar.parSpaceHiBound = log10([0.7   , 2]);
sodaPar.stateNamesKF = expandVar('presHead1',repmat(true,[nRows,nCols,1]));
sodaPar.stateNamesNOKF = [expandVar('presHead1Old',repmat(true,[nRows,nCols,1])),...
                          expandVar('presHead52',repmat(true,[nRows,nCols,nLays-1])),...
                          'discharge','vInit','vPrec','vSink','vSeep','vEnd'];
sodaPar.stateSpaceLoBound = -2.50*ones(1,nRows*nCols);
sodaPar.stateSpaceHiBound = upperstatebound'; %(upperstatebound = soildepth)
sodaPar.nMeasurements = 30;
sodaPar.modelName = 'swmsRunFun';
sodaPar.obsPoints = {'timeSimVec',timeObsVec};
sodaPar.simPoints = {'timeSimVec',timeSimVec};
sodaPar.measErr = 1e-20;
sodaPar.measNames = {'artifObs_h'};
sodaPar.nModelEvalsMax = 3280;
sodaPar.plotEnsemble = false;
sodaPar.optMethod = 'likelihood';
sodaPar.nCompl = 4;
sodaPar.nSamples = 40;
sodaPar.nEnsembleMembers = 1;%ceil(1.2*numel(sodaPar.stateNamesKF))+1;
sodaPar.verboseOutput = true;
sodaPar.initState = 'reference';
sodaPar.initStateRef = presHead1(:);
sodaPar.initStateErr = 0;
sodaPar.nOutputs = nRows*nCols;
sodaPar.errModelCallStr = 'stochForce = stochForceWhiteNoise(sodaPar,iIter,stochForce);';
sodaPar.stochForcePars = {'parPhi',0;'covArray',eye(numel(sodaPar.stateNamesKF))*1e-16};
%sodaPar.statesViz = [40,50,94];
sodaPar.forceRecalcFraction = 0;
sodaPar.vizWindow = timeSimVec([1,end]);
sodaPar.badScore = 1e12;
sodaPar.crashAndBurn = false;
sodaPar.drawInterval = 5;
sodaPar.ssp.nWorkers = 4;

sodaParBasic = sodaPar;

for k=1%[1,2]

    sodaPar = sodaParBasic;
    
    sodaPar.optEndTime = datenum(now+3/24);
    switch k
        case 1
            sodaPar.startFromUniform = false;
%             sodaPar.objCallStrs = {['gridOfErrColorTwoSpace(sodaPar,ensemblePred-mConstants{1,2},parVec,',...
%                 char(39),'modeStr',char(39),',',char(39),'error',char(39),',',...
%                 char(39),'yUnitStr',char(39),',',char(39),'m',char(39),');objSco',...
%                         're = objectivefun1(sodaPar,mStatesNOKF,vInnTrue,vSinkTrue);'];...
%                         'objScore = objectivefun2(sodaPar,mStatesNOKF,vSeepTrue);';...
%                         'objScore = objectivefun3(sodaPar,mStatesNOKF,mConstants);'};
            sodaPar.objCallStrs = {['objScore = objectivefun1(sodaPar,mStatesNOKF,vInnTrue,vSinkTrue);'];...
                        'objScore = objectivefun2(sodaPar,mStatesNOKF,vSeepTrue);';...
                        'objScore = objectivefun3(sodaPar,mStatesNOKF,mConstants);'};
            sodaPar.modeStr = 'scemua';                    
        case 2
            sodaPar.startFromUniform = false;
            sodaPar.objCallStrs = {['gridOfInnColorTwoSpace(sodaPar,mStatesInn,parVec,',...
                char(39),'modeStr',char(39),',',char(39),'innovations',char(39),',',...
                char(39),'yUnitStr',char(39),',',char(39),'m',char(39),');objSco',...
                        're = objectivefun1(sodaPar,mStatesNOKF,vInnTrue,vSinkTrue);'];...
                        'objScore = objectivefun2(sodaPar,mStatesNOKF,vSeepTrue);'};
%             sodaPar.objCallStrs = {['objScore = objectivefun1(sodaPar,mStatesNOKF,vInnTrue,vSinkTrue);'];...
%                         'objScore = objectivefun2(sodaPar,mStatesNOKF,vSeepTrue);'};
            sodaPar.modeStr = 'reset';
    end
    
    [evalResults,critGelRub,sequences,metropolisRejects,sodaPar] = soda(sodaPar);

    switch sodaPar.modeStr
        case 'soda'
            save(['soda-results-',num2str(size(evalResults,1)),'-',...
            num2str(sodaPar.nCompl),'-',...
            num2str(sodaPar.nSamples),'-',...
            num2str(sodaPar.nEnsembleMembers),'-',...
            num2str(sodaPar.randSeed),'-',...
            sodaPar.initState,...
            '.mat']);
        case 'scemua'
            save(['scemua-results-',num2str(size(evalResults,1)),'-',...
            num2str(sodaPar.nCompl),'-',...
            num2str(sodaPar.nSamples),'-',...
            num2str(sodaPar.randSeed),...
            '.mat']);
        case 'reset'
            save(['reset-results-',num2str(size(evalResults,1)),'-',...
            num2str(sodaPar.nCompl),'-',...
            num2str(sodaPar.nSamples),'-',...
            num2str(sodaPar.randSeed),...
            '.mat']);
    end

    subplotScreen(1,1,1)
    sodaPlotSeq(sodaPar,sequences,metropolisRejects,...
        'plotMode','subaxes','checkered',true)
    set(gcf,'inverthardcopy','off','paperpositionmode','auto')
    print([sodaPar.modeStr,'-sodaPlotSeq.eps'],'-depsc2','-r300','-loose')
    print([sodaPar.modeStr,'-sodaPlotSeq.png'],'-dpng','-r300')

    subplotScreen(1,1,1)
    sodaMargHist(sodaPar,evalResults)
    set(gcf,'inverthardcopy','off','paperpositionmode','auto')
    print([sodaPar.modeStr,'-sodaMargHist.eps'],'-depsc2','-r300','-loose')
    print([sodaPar.modeStr,'-sodaMargHist.png'],'-dpng','-r300')
    
end



sodaPar.objTexNames = {'RMSE_{total sink volume}','RMSE_{discharge}','RMSE_{pressure head}'};

subplotScreen(1,1,1,'rect',[100    100   900   500])
handles = sodaMatrixOfScatter(sodaPar,'eval-par',sequences,metropolisRejects,...
    'nHistory','all','yTrue',{log10([0.01;0.30]),log10(0.35)});
set(handles(1),'ylim',[sodaPar.parSpaceLoBound(1),sodaPar.parSpaceHiBound(1)])
set(handles(2),'ylim',[sodaPar.parSpaceLoBound(2),sodaPar.parSpaceHiBound(2)])
set(gcf,'inverthardcopy','on','paperpositionmode','auto')
print([sodaPar.modeStr,'MatrixOfScatterEvalPar.eps'],'-depsc2','-r300','-loose')
print([sodaPar.modeStr,'MatrixOfScatterEvalPar.png'],'-dpng','-r300')



subplotScreen(1,1,1,'rect',[100    100   700   600])
paretoOptimal = evalResults(:,sodaPar.paretoCol)>-2;
TMP1 = evalResults;
TMP1(~paretoOptimal,:)=NaN;
for iCompl = 1:sodaPar.nCompl
    TMP2(:,:,iCompl) = TMP1(iCompl:sodaPar.nCompl:end,:);
end

TMP3 = evalResults;
TMP3(paretoOptimal,:)=NaN;
for iCompl = 1:sodaPar.nCompl
    TMP4(:,:,iCompl) = TMP3(iCompl:sodaPar.nCompl:end,:);
end
handles = sodaMatrixOfScatter(sodaPar,'obj-obj',TMP2,TMP4,...
    'nHistory','all','excludeDiag',true,'excludeTriup',true,...
    'colorList',zeros(sodaPar.nCompl,3));
%set(handles(3,2),'xtick',[0.01:0.005:0.03])
set(handles(2,1),'xlim',[0,1],'ylim',[0,0.025]);
set(handles(3,1),'xlim',[0,1],'ylim',[0,0.032]);
set(handles(3,2),'xlim',[0,0.025],'ylim',[0,0.032]);
% set(handles(2,1),'yticklabel',num2str(get(handles(2,1),'ytick')','%.3f'))
% set(handles(3,1),'yticklabel',num2str(get(handles(3,1),'ytick')','%.3f'))


set(gcf,'inverthardcopy','on','paperpositionmode','auto')
print([sodaPar.modeStr,'MatrixOfScatterObjObj.eps'],'-depsc2','-r300','-loose')
print([sodaPar.modeStr,'MatrixOfScatterObjObj.png'],'-dpng','-r300')




