function [rscoredata]=rscore_mn(pfgnmrdata,ncomp,Specrange,Options)
%   [rscoredata]=rscore_mn(pfgnmrdata,ncomp,Specrange,Options,Nug,Fixed,Dinit)
%   R-SCORE (Speedy Component Resolution) fitting of
%   PFG-NMR diffusion data (aka DOSY data)
%
%   -------------------------------INPUT--------------------------------------
%   pfgnmrdata      Data structure containing the PFGNMR experiment.
%                   containing the following members:
%                       filename: the name and path of the original file
%                       np: number of (real) data points in each spectrum
%                       wp: width of spectrum (ppm)
%                       sp: start of spectrum (ppm)
%                       dosyconstant: gamma.^2*delts^2*DELTAprime
%                       Gzlvl: vector of gradient amplitudes (T/m)
%                       ngrad: number of gradient levels
%                       Ppmscale: scale for plotting the spectrum
%                       SPECTRA: spectra (processed)
%
%   ncomp           Number of componets (spectra) to fit.
%
%
%   ---------------------------INPUT - OPTIONAL-------------------------------
%
%   Specrange       The spectral range (in ppm) in which the score fitting
%                   will be performed. if set to [0] defaults will be used.
%                   DEFAULT is [sp wp+sp];
%   Options         Optional parameters%
%                   Options(1): Method to obtain the initial guesses
%                       0 (DEFAULT) obtained by monoexponential fitting of
%                         the summed amplitudes. Centered on the fitted
%                         value and differing by a factor of 1.5
%                       1 random values between 0 and 10 (seconds).
%                    
%                   Options(2): Constraints on the spectra
%                       0 No constraint (DEFAULT)
%                       1 Non negativity
%                       2 Non negativity - fast algorithm. Requires the
%                         N-way Toolbox by R. Bro (http://www.models.kvl.dk)
%                   Options(3): Shape of the decay
%                       0 T2 (2 paramteter exponential decay)
%                       1 T1 (3 paramater inversion recovery
%               
%              
%  
%
%   --------------------------------------OUTPUT------------------------------
%   rscoredata      Structure containg the data obtained after score with
%                  the follwing elements:
%
%                  COMPONENTS: The fitted spectra. In the optimal case
%                              representing the true spectra of the
%                              indiviual componets in the mixture
%                  DECAYS:     The decys of the fitted spectra
%                  fit_time:   The time it took to run the fitting
%                  d2: vector of echo times (seconds)
%                  wp: width of spectrum (ppm)
%                  sp: start of spectrum (ppm)
%                  Ppmscale: scale for plotting the spectrum
%                  filename: the name and path of the original file
%                  Rval: fitted relaxation coeffcients (seconds))
%                  ncomp: number of fitted components
%                  relp: relavtive proportion of the fitted components
%                  Options: options used for the fitting
%
%   Example:
%
%   See also: dosy_mn, score_mn, decra_mn, mcr_mn, varianimport,
%             brukerimport, jeolimport, peakpick_mn, dosyplot_mn,
%             dosyresidual, dosyplot_gui, scoreplot_mn, decraplot_mn,
%             mcrplot_mn
%
%   References:
%   Speedy Component Resolution: An Improved Tool for Processing
%   Diffusion-Ordered SPectroscopy Data. Anal. Chem., 2008 (In Press)
%
%   This is a part of the GNAT
%   Copyright  2017  <Mathias Nilsson>%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License along
%   with this program; if not, write to the Free Software Foundation, Inc.,
%   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
%   Dr. Mathias Nilsson
%   School of Chemistry, University of Manchester,
%   Oxford Road, Manchester M13 9PL, UK
%   Telephone: +44 (0) 161 306 4465
%   Fax: +44 (0)161 275 4598
%   mathias.nilsson@manchester.ac.uk
%
% SCORE processing of
% PFG-NMR diffusion data (aka DOSY data)
%
% see also 'core_mn,decra_mn,dosy_mn and mcr_mn'
%


tic;
disp('RSCORE')

Opts=Options;




rscore_params=[];
rscore_params.ncomp=ncomp;




if Specrange==0
    %do nothing
else
    if length(Specrange)~=3 && length(Specrange)~=2
        error('RSCORE: Specrange should have 2 or 3 elements')
    end
    if Specrange(1)<pfgnmrdata.sp
        disp('RSCORE: Specrange(1) is too low. The minumum will be used')
        Specrange(1)=pfgnmrdata.sp;
    end
    if Specrange(2)>(pfgnmrdata.wp+pfgnmrdata.sp)
        disp('RSCORE: Specrange(2) is too high. The maximum will be used')
        Specrange(2)=pfgnmrdata.wp+pfgnmrdata.sp;
    end
    for k=1:length(pfgnmrdata.Ppmscale)
        if (pfgnmrdata.Ppmscale(k)>Specrange(1))
            begin=k-1;
            break;
        end
    end
    
    for k=begin:length(pfgnmrdata.Ppmscale)
        if (pfgnmrdata.Ppmscale(k)>=Specrange(2))
            endrange=k;
            break;
        end
    end
    %make a new stucture
    pfgnmrdata.sp=pfgnmrdata.Ppmscale(begin);
    pfgnmrdata.wp=pfgnmrdata.Ppmscale(endrange)-pfgnmrdata.Ppmscale(begin);
    pfgnmrdata.Ppmscale=pfgnmrdata.Ppmscale(begin:endrange);
    pfgnmrdata.SPECTRA=pfgnmrdata.SPECTRA(begin:endrange,:);
    pfgnmrdata.np=length(pfgnmrdata.Ppmscale);
end

rscore_params.pfgnmrdata=pfgnmrdata;

switch Opts(3)
    
    case 0
        disp('RSCORE: T2 fit (2 parameter exponential decay)')
        rscore_params.model=@t2;
        rscore_params.model2=@t2_monofit;
        
    case 1
        disp('RSCORE: T1 fit (3 parameter inversion recovery)')
        rscore_params.model=@t1;
        rscore_params.model2=@t1_monofit;
    otherwise
        error('RSCORE: Illegal Options(3)')
end







switch Opts(1)
    
    case 0
        disp('RSCORE: Using monoexponential fitting for initialisation of difussion coefficients')
        opts_fmin=optimset('fminsearch');
        opts_fmin=optimset(opts_fmin,'Display','Final');
        opts_fmin=optimset(opts_fmin,'MaxFunEvals',10000);
        
        tmp=sum(pfgnmrdata.SPECTRA);
        
        fitparm(1)=tmp(1);
        fitparm(2)= 4;
        switch Opts(3)
            case 1
                %T1 fit (3 parameter inversion recovery)')
                fitparm(3)= -fitparm(1);
        end
        fitparm=real(fitparm);
        [fitd]= fminsearch(@rscore_params.model2,fitparm,opts_fmin,pfgnmrdata.d2,sum(pfgnmrdata.SPECTRA));
        
        rcenter=fitd(2);
        if mod(ncomp,2)==0 %even number
            rval=linspace(rcenter/(ncomp/2*1.5),rcenter*(ncomp/2*1.5),ncomp);
        elseif ncomp > 2
            rval=linspace(rcenter/((ncomp-1)/2*1.5),rcenter*((ncomp-1)/2*1.5),ncomp);
        else
            rval=rcenter;
        end
        rval(ncomp+1)=-1  ;
        nguesses=1;
        
        
    case 1
        disp('RSCORE: Using random values for initialisation of difussion coefficients')
        %rng('shuffle');
        %rval=rand(1,ncomp)*10; %random starting values for T1/T2
        %if T1 inversion recovery
       % rval(ncomp+1)=-1;
       
        opts_fmin=optimset('fminsearch');
        opts_fmin=optimset(opts_fmin,'Display','Final');
        opts_fmin=optimset(opts_fmin,'MaxFunEvals',10000);
        
        tmp=sum(pfgnmrdata.SPECTRA);
        
        fitparm(1)=tmp(1);
        fitparm(2)= 4;
        switch Opts(3)
            case 1
                %T1 fit (3 parameter inversion recovery)')
                fitparm(3)= -fitparm(1);
        end
        fitparm=real(fitparm);
        [fitd]= fminsearch(@rscore_params.model2,fitparm,opts_fmin,pfgnmrdata.d2,sum(pfgnmrdata.SPECTRA));
        
        nguesses=100;
        rcenter=fitd(2);
        rval=abs(normrnd(rcenter,0.05,[nguesses ncomp]));
        rval(:,ncomp+1)=0;
        
        
    otherwise
        error('RSCORE: Illegal Options(1)')
end

switch Opts(2)
    
    case 0
        disp('RSCORE: No constraints')
        nnegflag=0;
    case 1
        disp('RSCORE: Non negativity constraint')
        nnegflag=1;
    otherwise
        error('RSCORE: Illegal Options(2)')
end

switch Opts(6)
    
    case 0
        disp('RSCORE: Using SCORE (residuals minimisation) inner loop')
        OUTSCORE=0;
    case 1
        disp('RSCORE: Using OUTSCORE (cross-talk minimisation) inner loop')
        OUTSCORE=1;
    otherwise
        error('RSCORE: Illegal Options(6)')
end

%lb=zeros(1,ncomp);
lb=[];
ub=[];
% ub=inf(1,ncomp);
% ub(ncomp+1)=0;
rscore_params.pfgnmrdata.d2;
numel(rscore_params.pfgnmrdata.d2);

opts1=optimset('fmincon');
opts1=optimset(opts1,'Display','Iter');
for p=1:nguesses
[relaxcoef(:,p) , resid(:,p)]= fmincon(@rscore_inner,rval(p,:),[],[],[],[],lb,ub,[],opts1,rscore_params);
end 
[resid, index]=min(resid);                                                % Here choosing only the minimum residual value out of all other residuals obtained by the different number of guesses
relaxcoef=relaxcoef(:,index);

C=zeros(numel(pfgnmrdata.d2),rscore_params.ncomp);
S=zeros(rscore_params.ncomp,rscore_params.pfgnmrdata.np);
for k=1:rscore_params.ncomp %#ok<*FXUP>
    inval(1)=1;
    inval(2)=relaxcoef(k);
    C(:,k)=rscore_params.model(inval,rscore_params.pfgnmrdata.d2);
end


X=pfgnmrdata.SPECTRA;
if nnegflag==1
    switch Opts(3)
        case 1
            %T1 fit (3 parameter inversion recovery)')
            for k=1:pfgnmrdata.np
                S(:,k)=lsqnonneg(C,(-X(k,:)'));
            end
            
        otherwise
            for k=1:pfgnmrdata.np
                S(:,k)=lsqnonneg(C,X(k,:)');
            end
    end
    
    
else
    switch Opts(3)
        case 1
            %T1 fit (3 parameter inversion recovery)')
            S=C\(-X');
            
        otherwise
            S=C\(X');
    end
    
end


% figure
% plot(S')

endt=toc;
h=fix(endt/3600);
m=fix((endt-3600*h)/60);
s=endt-3600*h - 60*m;
fit_time=[h m s];
disp(['RSCORE: Fitting time was: ' num2str(fit_time(1)) ' h  ' num2str(fit_time(2)) ' m and ' num2str(fit_time(3),2) ' s']);


rscoredata.RESIDUALS=X-S'*C';
rscoredata.filename=pfgnmrdata.filename;
rscoredata.Gzlvl=pfgnmrdata.Gzlvl;
rscoredata.COMPONENTS=S;


switch Opts(3)
    case 1
        %T1 fit (3 parameter inversion recovery)')
        rscoredata.DECAYS=-C;
        
    otherwise
        rscoredata.DECAYS=C;
end
rscoredata.fit_time=fit_time;
rscoredata.Ppmscale=pfgnmrdata.Ppmscale;
rscoredata.ncomp=ncomp;
rscoredata.Options=Opts;
rscoredata.Rval=relaxcoef;
rscoredata.resid=resid;
%find the relative contributions
relp=sum(abs(rscoredata.COMPONENTS),2);
rscoredata.relp=relp;
rscoredata.np=pfgnmrdata.np;
rscoredata.sp=pfgnmrdata.sp;
rscoredata.wp=pfgnmrdata.wp;
rscoredata.SPECTRA=pfgnmrdata.SPECTRA;
rscoredata.d2=pfgnmrdata.d2;


rscoreplot(rscoredata);



%rscoredata=[];
return
    function [resid]=rscore_inner(Tval,params)
        
        C=zeros(numel(params.pfgnmrdata.d2),params.ncomp);
        S=zeros(params.ncomp,params.pfgnmrdata.np);
        X=params.pfgnmrdata.SPECTRA;
        
        for k=1:(params.ncomp)
            inval(1)=1;
            inval(2)=Tval(k);
            inval(3)=Tval(params.ncomp+1);
            
            C(:,k)=params.model(inval,params.pfgnmrdata.d2);
        end
        
        if nnegflag==1
            for k=1:params.pfgnmrdata.np
                S(:,k)=lsqnonneg(C,X(k,:)');
            end
        elseif nnegflag==2
            for k=1:params.pfgnmrdata.np
                S(:,k)=fastnnls(C,X(k,:)');
            end
        else
            S=C\X';
        end
        
        %% Minimisation Criteria
        if params.ncomp<=1  % Can't OUTSCORE one component, revert to SCORE
            OUTSCORE=0;
        end
        
        
        switch OUTSCORE
            case 1 % area normalisation (OUTSCORE)
                
                for m=1:ncomp
                    S(m,:)=abs(S(m,:))./sum(abs(S(m,:)));
                end
                CrossTalk=ones(1,pfgnmrdata.np);
                
                CompPairs=nchoosek(1:ncomp,2);
                [rows , ~]=size(CompPairs);
                for n=1:rows % the sum of pairwise .* of spectra
                    CrossTalk=CrossTalk+(S(CompPairs(n,1),:).*S(CompPairs(n,2),:));
                end
                resid = sum(CrossTalk);
            case 0 % Residuals Minimisation (SCORE)
                M=C*S;
                resid=sum(sum(((M'-X).^2)));
        end
        
    end
    function sse = t2_monofit(a,xdata,ydata)
        % for fitting a  pure t2 exponential
        xdata=xdata';
        y =  a(1)*exp(-xdata/a(2));
        
        sse=ydata-y;
        sse=sum(sse.^2,2);
    end
    function [y, J] = t2(a,xdata)
        %T2:  M(t) = (M(0) - M(inf))*exp(-t/T2) + M(inf)
        %where M(0) is the magnetization at time zero (i.e., the full
        %magnetization excited by the observe pulse) and M(inf) is the
        %xy-magnetization at infinite time (zero unless the peak is sitting on an
        %offset baseline). So using:
        %T2: M(t) = M(0)*exp(-t/T2)
        xdata=xdata';
        y =a(1)*exp(-xdata/a(2));
        
        if nargin>1
            J=zeros(length(xdata),length(a));
            J(:,1)=exp(-xdata./a(2));
            J(:,2)=(a(1)*exp(-xdata/a(2)).*xdata)./a(2).^2;
            
        end
    end
    function sse = t1_monofit(a,xdata,ydata)
        % for fitting a  pure t1 exponential
        xdata=xdata';
        y =  (a(1)-a(3))*exp(-xdata/a(2)) +a(3);
        
        sse=ydata-y;
        sse=sum(sse.^2,2);
    end
    function [y, J] = t1(a,xdata)
        % for fitting a number of pure exponentials
        % this function works for both T1 and T2 but the initial guesses
        % should be different. However, I will use a 2 parameter for T2 as bot signal
        %and noise should decay to zero.
        %
        %T1: M(t)=(M(0) - M0)*exp(-t/T1)+M0
        %where M0 is the equilibrium Z magnetization and M(0) is the magnetization
        %at time zero (e.g., immediately after the 180° pulse for an inversion
        %recovery T1 experiment). Notice that this equation will fit inversion
        %recovery data (for which M(0) is approximately equal to -M0) or saturation
        %recovery data (for which M(0) is 0).
        %
        %T2:  M(t) = (M(0) - M(inf))*exp(-t/T2) + M(inf)
        %where M(0) is the magnetization at time zero (i.e., the full
        %magnetization excited by the observe pulse) and M(inf) is the
        %xy-magnetization at infinite time (zero unless the peak is sitting on an
        %offset baseline).        
        
        xdata=xdata';
        %y=zeros(1,length(xdata));
        y =(a(1) - a(3))*exp(-xdata/a(2)) + a(3);
        if nargin>1
            J=zeros(length(xdata),length(a));
            J(:,1)=exp(-xdata./a(2));
            J(:,2)=((a(1) - a(3))*exp(-xdata/a(2)).*xdata)./a(2).^2;
            J(:,3)=1- exp(-xdata./a(2));
            
        end
    end

end


