rates=readtable('D://bin_cc_rates2/procrates.txt');
rates.Properties.VariableNames={'M-3','M-2','M-1','M0','M1','M2','M3','M4','M5','M6','M7','M8','M9'};

moments=readtable('D://bin_cc_rates2/bin.txt');
names=strcat({'c'},{'M-3','M-2','M-1','M0','M1','M2','M3','M4','M5','M6','M7','M8','M9'});
names2=strcat({'r'},{'M-3','M-2','M-1','M0','M1','M2','M3','M4','M5','M6','M7','M8','M9'});
moments.Properties.VariableNames={names{:},names2{:}};

%Columns correspond to moment rates for the -3rd to 9th moments
r.tmsc=rates(1:6:end,:); %self-collection cloud
r.tmsr=rates(2:6:end,:); %self-collection rain
r.tmauc=rates(3:6:end,:); %autoconversion cloud
r.tmaur=rates(4:6:end,:); %autoconversion rain
r.tmacc=rates(5:6:end,:); %accretion cloud
r.tmacr=rates(6:6:end,:); %accretion rain
r.tmf=array2table(-table2array(r.tmsc)-table2array(r.tmsr)-table2array(r.tmauc)...
    +table2array(r.tmaur)-table2array(r.tmacc)+table2array(r.tmacr)); %total budget
r.tmc=array2table(table2array(r.tmsc)+table2array(r.tmauc)+table2array(r.tmacc)); %budget cloud
r.tmr=array2table(-table2array(r.tmsr)+table2array(r.tmaur)+table2array(r.tmacr)); %budget rain
r.tmf.Properties.VariableNames={'M-3','M-2','M-1','M0','M1','M2','M3','M4','M5','M6','M7','M8','M9'};
r.tmc.Properties.VariableNames={'M-3','M-2','M-1','M0','M1','M2','M3','M4','M5','M6','M7','M8','M9'};
r.tmr.Properties.VariableNames={'M-3','M-2','M-1','M0','M1','M2','M3','M4','M5','M6','M7','M8','M9'};

%% 
procs={'tmsc','tmsr','tmauc','tmaur','tmacc','tmacr','tmf','tmf','tmf','tmf','tmc','tmr'};
procns={'SelfColl','SelfColl','Autocon','Autocon','Accret','Accret','All','All Cloud Dom.','All Cloud Rain Mix','All Rain Dom.','All','All'};

nmom=13;

%This code written to be fairly flexible in analysis. Only this settings
%is used in the paper:
mtyp='f'; %c=cloud, r=rain, f=full

if strcmp(mtyp,'c')
    jj=[5,6]; %Just the processes that depend on both cloud and rain
elseif strcmp(mtyp,'r')
    jj=[5,6];
else
    jj=7;
end
pms=-3:9;
colors=get(0,'defaultaxescolororder');

il=0;
il2=0;
for j=jj
    RMSEx=zeros(nmom+1,nmom+1);

    if strcmp(mtyp,'c')
        if procs{j}(end)=='c' %then it is a cloud moment process
            inds=moments.cM3>0 & abs(r.(procs{j}).M0(:)./r.tmc.M0(:))>0.01;
        else % then it is a rain moment process
            inds=moments.cM3>0 & abs(r.(procs{j-1}).M0(:)./r.tmc.M0(:))>0.01;
        end
    elseif strcmp(mtyp,'r')
        if procs{j}(end)=='c'
            inds=moments.cM3>0 & abs(r.(procs{j}).M0(:)./r.tmc.M0(:))>0.01;
        elseif strcmp(procs{j},'tmsr')
            inds=moments.rM3>0 & abs(r.(procs{j}).M0(:)./r.tmr.M0(:))>0.01;
        else
            inds=moments.cM3>0 & abs(r.(procs{j-1}).M0(:)./r.tmc.M0(:))>0.01;
        end
        
    elseif strcmp(mtyp,'f')
        if j==7
            inds=(moments.cM3+moments.rM3)*pi/6*1e6>1;
        elseif j==8
            inds=moments.cM3./(moments.cM3+moments.rM3)>.99;
        elseif j==9
            inds=moments.cM3./(moments.cM3+moments.rM3)>.1 & moments.cM3./(moments.cM3+moments.rM3)<=.99;
        else
            inds=moments.cM3./(moments.cM3+moments.rM3)<=.1;
        end
    end
    
    if strcmp(mtyp,'c')
        M0= moments.cM0(inds);
        M3= moments.cM3(inds);
        M0b=moments.rM0(inds);
        M3b=moments.rM3(inds);
        yl='Cloud';
    elseif strcmp(mtyp,'r')
        M0= moments.rM0(inds);
        M3= moments.rM3(inds);
        M0b=moments.cM0(inds);
        M3b=moments.cM3(inds);
        yl='Rain';
    elseif strcmp(mtyp,'f')
        M0=moments.cM0(inds)+moments.rM0(inds);
        M3=moments.cM3(inds)+moments.rM3(inds);
        yl='Full';
    else
        display('invalid moment type')
        break
    end
    
    nM3=M3./M0;
    if strcmp(mtyp,'c')||strcmp(mtyp,'r');nM3b=M3b./M0b;end
    
    pvn=r.(procs{j}).Properties.VariableNames;
    for imr=1:13 %loop over moment rates
        
        clear inputarray
        mypvn=pvn{imr};
        
        if strcmp(mypvn,'M3') && any(strcmp(procs{j},{'tmsc','tmsr','tmaur','tmacr','tmf'}))
            R2x(imr)=NaN;
            continue
        end
        if strcmp(mypvn,'M0') && strcmp(procs{j},'tmacr')
            R2x(imr)=NaN;
            continue
        end
        
        y=log(abs(r.(procs{j}).(mypvn)(inds)));        
        
        clear RMSE R2 RMSEm6 R2m6
        RMSE=nan(nmom+1,nmom+1);
        RMSEm6=nan(nmom+1,nmom+1);
        for i=1:13 %loop over 1st extra predictor
            if strcmp(mtyp,'f');kst=i;else;kst=1;end
            for k=kst:13 %loop over 2nd extra predictor
                clear inputarray

                if ~strcmp(mypvn,'M0') && ~strcmp(mypvn,'M6')
                    continue
                end
                
                inputarray(:,2)=log(M0);
                inputarray(:,3)=log(nM3);
                inputarray(:,1)=1;
                iia=4;
                if any(strcmp(procs{j},{'tmacc','tmacr'}))
                    inputarray(:,4)=log(M0b);
                    inputarray(:,5)=log(nM3b);
                    iia=iia+2;
                end
                
                vnx=moments.Properties.VariableNames{i};
                vny=moments.Properties.VariableNames{k};
                vnx2=moments.Properties.VariableNames{i+13};
                vny2=moments.Properties.VariableNames{k+13};
                if strcmp(mtyp,'c')
                    Mx=moments.(vnx)(inds);
                    My=moments.(vny2)(inds);
                elseif strcmp(mtyp,'r')
                    vnx2=moments.Properties.VariableNames{i+13};
                    vny2=moments.Properties.VariableNames{k+13};
                    Mx=moments.(vnx)(inds);
                    My=moments.(vny2)(inds);
                else
                    vnx2=moments.Properties.VariableNames{i+13};
                    vny2=moments.Properties.VariableNames{k+13};
                    Mx=moments.(vnx)(inds)+moments.(vnx2)(inds);
                    My=moments.(vny)(inds)+moments.(vny2)(inds);
                end
                
                z=i-4;
                if z==0 || z==3
                    R2(i,k)=NaN;
                    RMSE(i,k)=NaN;
                    continue
                end
                z2=k-4;
                if z2==0 || z2==3
                    R2(i,k)=NaN;
                    RMSE(i,k)=NaN;
                    continue
                end
                
                nMx=Mx./(M0.^(1-z/3).*M3.^(z/3));
                if strcmp(mtyp,'f')
                    nMy=My./M0./(nMx.*(M3./M0).^(z2/3));
                else
                    nMy=My./(M0.^(1-z2/3).*M3.^(z2/3));
                end
                
                inputarray(:,iia)=log(nMx);
                
                if strcmp(mtyp,'f')
                    if i~=k
                        inputarray(:,iia+1)=log(nMy);
                    end
                else
                    inputarray(:,iia+1)=log(nMy);
                end
                
                [~,~,~,~,stats]=regress(y,inputarray);
                
                if strcmp(mypvn,'M0') || strcmp(mypvn,'M3')|| strcmp(mypvn,'M6')
                    R2(i,k)=stats(1);
                    RMSE(i,k)=stats(4);
                else
                    R2x(i,k)=R2x(i,k)+stats(1);
                    if i==k
                        RMSEx(i,k)=RMSEx(i,k)+stats(4);
                    else
                        RMSEx(i,k)=RMSEx(i,k)+stats(4)/2;
                    end
                end
            end        
        end
        
        if strcmp(mypvn,'M0') || strcmp(mypvn,'M3') || strcmp(mypvn,'M6')
            if strcmp(mypvn,'M0');figure('position',[1 1 300 300]);subplot(2,2,1);
            else; subplot(2,2,2); end
            
            RMSE(RMSE==0)=NaN;
            RMSE(1:3,:)=NaN;
            RMSE(:,1:3)=NaN;
            for iii=1:length(RMSE)
                RMSE(iii,iii)=NaN;
            end
            pcolor(-3.5:9.5,-3.5:9.5,RMSE)
            xlabel(['Normalized ',yl,' Moment'])
            ylabel(['Normalized ',yl,' Moment'])
            colorbar
            if strcmp(mypvn,'M0')
                title(['(a) RMSE ',mypvn,' Tend.'])
                RMSE0=RMSE;
            else
                title(['(b) RMSE ',mypvn,' Tend.'])  
                RMSE6=RMSE;
            end
            set(gca,'xlim',[0.5 9.5],'ylim',[0.5 9.5])
        end
    end    
end