%% CALLED BY NCfitImCFparforFailGbl2.m %%
function [resMatStd, resMat, selTimesStd, selIntensStd, FiltTimesELr, NormIntensELr] =...
    NCscurImCF_3parfor(dataMatrix, AUCfinalTime, currSpotArea, sols, bl, minTime)
%Major revision for Early-Late data cuts to improve accuracof 'r'. Removed legacy iterative method.
%Significant Modification for Parfor  
%***************************************************************

%##########################################################################
 %*******************************************  New Stage 1***************************************************************
%Preallocate
resMatStd= zeros(1,27);
resMat= zeros(1,27);
%Set internal variables sent to matlab fit function 
me=200;
meL=750;
mi=25;  %50
miL=250;
%***********************************
 rmsStg1=0;
 rmsStg1I(1)= 0;
 slps=1;

  filterTimes=[];
  normIntens=[];
  nn=1;
  numFitTpts=0;
%Build filterTimes and normIntens from spot dataMatrix selection codes produced in filter section  
for n=1:size(dataMatrix,2)
    if (((dataMatrix(3,n)==1))||(dataMatrix(3,n)==3)||(dataMatrix(3,n)==2)...
            ||(dataMatrix(3,n)==0))
        filterTimes(nn)= dataMatrix(1,n);
        normIntens(nn)=dataMatrix(4,n);  
        nn=nn+1;
    end
end

%------------------------------------------------------------------

%++++++++++++++++++++++++++++++++++++++++

filterTimes=filterTimes';  
selTimesStd=filterTimes; 
normIntens=normIntens';  
selIntensStd=normIntens;

%normIntens  %debugging parfor gbl  200330 good values
%afgj

  lastTptUsed= 1;
  lastIntensUsed= 1;
  thresGT2TmStd=0;  
 try
    
    lastTptUsed=max(filterTimes);
    lastIntensUsed=normIntens(length(normIntens));
    lastIntensUsedStd= lastIntensUsed;
    
    lastTptUsedStd= lastTptUsed;
    Tpt1Std= filterTimes(1);    
    numFitTptsStd=nnz((normIntens(:)>=0)==1);
    thresGT2 = find(((normIntens(:)>2)==1), 1);
    if isempty(thresGT2)
     thresGT2TmStd=0;
    else
    thresGT2TmStd = filterTimes(thresGT2);
    end
    numTptsGT2Std = nnz((normIntens(:)>=2)==1); %nnz(filterTimes(find(filterTimes>=thresGT2Tm)));
    K_Guess = max(normIntens);
    numTimePts = length(filterTimes);
    opts =  fitoptions('Method','Nonlinear','Robust','On',...
        'DiffMinChange',1.0E-11,'DiffMaxChange',0.001,...
        'MaxFunEvals',me, 'MaxIter', mi, 'TolFun', 1.0E-12, 'TolX', 1.0E-10, 'Lower', [K_Guess*0.5,0,0], 'StartPoint', [K_Guess,filterTimes(floor(numTimePts/2)),0.30], 'Upper', [K_Guess*2.0,max(filterTimes),1.0],'Display','off');
    ftype = fittype('K / (1 +  exp(-r* (t - l )))','independent','t','dependent',['K','r','l'],'options',opts);

    % carry out the curve fitting process
    [fitObject, errObj] = fit(filterTimes,normIntens,ftype);
    coeffsArray = coeffvalues(fitObject);
    rmsStg1=errObj.rsquare;
    rmsStg1I(slps)= errObj.rsquare;
    sDat(slps,1)=errObj.rsquare; 
    K = coeffsArray(1); sDat(slps,2)= coeffsArray(1);  % Carrying Capacity
    l = coeffsArray(2); sDat(slps,3)= coeffsArray(2); %lag time
    r = coeffsArray(3); sDat(slps,4)= coeffsArray(3); % rateS

    % integrate (from first to last time point)
    numVals = size(filterTimes);
    numVals = numVals(1);
    t_begin = 0;
    t_end = AUCfinalTime;
    AUC = (K/r*log(1+exp(-r*(t_end-l)))-K/r*log(exp(-r*(t_end-l)))) - (K/r*log(1+exp(-r*(t_begin-l)))-K/r*log(exp(-r*(t_begin-l))));
    MSR = r;

    rsquare = errObj.rsquare;  
    confObj = confint(fitObject,0.9); % get the 90% confidence 
    
    NANcond=0; stdNANcond=0;  %stdNANcond added to relay not to attempt ELr as there is no curve to find critical point
    confObj_filtered=confObj; 
    Klow = confObj(1,1); sDat(slps,5)=confObj(1,1);
    Kup = confObj(2,1); sDat(slps,6)=confObj(2,1);
    llow = confObj(1,2); sDat(slps,7)=confObj(1,2);
    lup = confObj(2,2); sDat(slps,8)=confObj(2,2);
    rlow = confObj(1,3); sDat(slps,9)=confObj(1,3);
    rup = confObj(2,3); sDat(slps,10)=confObj(2,3);
    if(isnan(Klow)||isnan(Kup)||isnan(llow)||isnan(lup)||isnan(rlow)||isnan(rup))
    NANcond=1; stdNANcond=1; %stdNANcond added to relay not to attempt ELr as there is no curve to find critical point
    end  
    
    %rup  %debugging parfor gbl 200330
    %Klow
    %asdfjj114
% {   
catch %ME
    % if no data is given, return zeros
   AUC = 0;MSR = 0;K = 0;r = 0;l = 0;rsquare = 0;Klow = 0;Kup = 0;
   rlow = 0;rup = 0;lup = 0;llow = 0; 
   NANcond=1; stdNANcond=1;  %stdNANcond added to relay not to attempt ELr as there is no curve to find critical point
end   %end Try
%}
    if (exist('K','var')&& exist('r','var') && exist('l','var'))
    t = (0:1:200);   
    Growth = K ./ (1 +  exp(-r.* (t - l )));
    fitblStd= min(Growth);   %jh diag
    end
   
    cutTm(1:4)= 1000;  %-1 means cuts not used or NA
%{    
l  %debugging parfor
K
r
Klow
k
%}    
%***Preserve for ResultsStd+++++++++++++++++++++++++++++

resMatStd(1)= AUC;
resMatStd(2)= MSR;
resMatStd(3)= K;
resMatStd(4)= r;
resMatStd(5)= l;
resMatStd(6)= rsquare;
resMatStd(7)= Klow;
resMatStd(8)= Kup;
resMatStd(9)= rlow;
resMatStd(10)= rup;
resMatStd(11)= llow;
resMatStd(12)= lup;
resMatStd(13)= currSpotArea;
resMatStd(14)= lastIntensUsedStd;  %filtNormIntens(length(filtNormIntens));

%spline fit unneccessary and removed;therefore No maxRateTime assoc'd w/spline fit
%try
%resMatStd(15)= maxRateTime;
%catch
    maxRateTime= 0; %[];  %Std shows []; ELr shows 0;  %parfor
    resMatStd(15)= 0; %maxRateTimestdMeth;
%end
 

resMatStd(16)= lastTptUsedStd;   %filterTimes(length(filterTimes));
     if isempty(Tpt1Std)
        Tpt1Std= 777;  %0.000002;  %
     end
resMatStd(17)= Tpt1Std;
resMatStd(18)= bl; %perform in the filter section of NCfitImCFparfor
resMatStd(19)= fitblStd; %Taken from NCfil... and not affected by NCscur...changes
resMatStd(20)= minTime; %Not affected by changes made in NCscur...for refined 'r'
resMatStd(21)= thresGT2TmStd;
resMatStd(22)= numFitTptsStd;
resMatStd(23)= numTptsGT2Std;

resMatStd(24)= 999; %The Standard method has no cuts .:.no cutTm  
resMatStd(25)= 999;
resMatStd(26)= 999;
resMatStd(27)= 999;


%if SwitchEvsEL==3  %Remove 'SwitchEvsEL==...' temporary SWITCH when Hartman decides what he wants
%*********************************************************************************
%ELr         New Experimental data through L+deltaS Logistic fit for 'Improved r'  Fitting
%*********************************************************************************
FiltTimesELr= [];  %{ii}= filterTimes;
NormIntensELr= [];  %{ii}= normIntens;  
normIntens= selIntensStd;
filterTimes= selTimesStd;
stdIntens= selIntensStd;
tmpIntens= selIntensStd;
stdTimes= selTimesStd;

if stdNANcond==0
%Determine critical points and offsets for selecting Core Data based on 
%Standard curve fit run.  Put diff into NImStartupImCF02.m calling source
%to reduce repeated execution since it doesn't change.
 %fd4=diff(sym('K / (1 +  exp(-r* (t - l )))'),4);
 %sols=solve(fd4);
tc1= eval(sols(2));
tc2= eval(sols(3));
LL= l;  %eval(sols(1));  %exactly the same as 'l' from std. fit method-Save time
rsTmStd= LL-tc1;   %%riseTime (first critical point to L)
deltS= rsTmStd/2;
    tc1Early= tc1-deltS;    %AKA- tc1AdjTm %2*tc1 -LL
    L_Late= LL+deltS;
    
           
    tc1EdatPt= find(filterTimes>(tc1Early),1,'first');
        cutTm(1)= filterTimes(2);
        cutDatNum(1)= 2;
        cutTm(2)= tc1Early;
        cutDatNum(2)= tc1EdatPt-1;
    
    L_LDatPt= find(filterTimes< L_Late,1,'last');
    tc2LdatPt= find(filterTimes< tc2+rsTmStd,1,'last');
    
        cutTm(3)= L_Late;
        cutDatNum(3)= L_LDatPt;
    
%Select Core Data Set  (Remove Early data before critical point)         
ints=[];
ints(1:L_LDatPt-tc1EdatPt+2)= (tmpIntens(L_LDatPt));
ints(2:end)= tmpIntens(tc1EdatPt:L_LDatPt);
ints(1)= tmpIntens(1);
tms=[];
tms(1:L_LDatPt-tc1EdatPt+2)= (stdTimes(L_LDatPt));
tms(2:end)= stdTimes(tc1EdatPt:L_LDatPt);
tms(1)= stdTimes(1);
%-----------------------------------------------
%Include/Keep late data that define K *********
if length(tmpIntens(tc2LdatPt:end))> 4
   KlastInts= tmpIntens(tc2LdatPt:end);
   KlastTms= stdTimes(tc2LdatPt:end);
   lengthKlast= length(tmpIntens(tc2LdatPt:end));
   ints(end:(end+ lengthKlast-1))= KlastInts;
   tms(end:(end+ lengthKlast-1 ))= KlastTms;
   
   cutTm(4)= tc2+rsTmStd;
   cutDatNum(4)= tc2LdatPt-1;
   
   
else
  lengthKlast= length(tmpIntens(tc2LdatPt-1:end)); 
  if lengthKlast>1
     KlastInts= tmpIntens(end-(lengthKlast-1):end);
     KlastTms= stdTimes(end-(lengthKlast-1):end);
     ints(end:(end+ lengthKlast-1 ))= KlastInts;
     tms(end:(end+ lengthKlast-1 ))= KlastTms;
  end

cutTm(4)= stdTimes(tc2LdatPt-1);
cutDatNum(4)= tc2LdatPt-2;   %length(stdTimes(end-(lengthKlast-1):end));
   
end
    

%************************************************
Ints=[];
Tms=[];
Ints= ints';
Tms= tms';

try
  filterTimes= Tms; filterTimes4= Tms; 
  normIntens= Ints;   normIntens4=Ints;
  
%----------------------------------------------------------------------------------------------
    
   %classic symetric logistic curve fit setup restated as COMMENTS for reference convenience----------------------------
    %opts= fitoptions is the same as for Std and so is redundant     
    %opts =  fitoptions('Method','Nonlinear','Robust','On',...
     %   'DiffMinChange',1.0E-11,'DiffMaxChange',0.001,...
      %  'MaxFunEvals',me, 'MaxIter', mi, 'TolFun', 1.0E-12, 'TolX', 1.0E-10, 'Lower', [K_Guess*0.5,0,0], 'StartPoint', [K_Guess,filterTimes(floor(numTimePts/2)),0.30], 'Upper', [K_Guess*2.0,max(filterTimes),1.0]);
    
    ftype = fittype('K / (1 +  exp(-r* (t - l )))','independent','t','dependent',['K','l','r'],'options',opts);
    fitObject=[]; errObj=[];
       % carry out the curve fitting process
    [fitObject, errObj] = fit(Tms,Ints,ftype);
      coeffsArray = coeffvalues(fitObject);
    r3 = coeffsArray(3); %sDat(slps,4)= coeffsArray(3); % rateS
    
   if (exist('K','var')&& exist('r','var') && exist('l','var'))
    t = (0:1:200);   
    GrowthELr = K ./ (1 +  exp(-r.* (t - l )));
    fitblELr= min(GrowthELr);   %jh diag
   end
   
   catch ME
    % if no data is given, return zeros
   AUC = 0;MSR = 0;K = 0;r = 0;l = 0;rsquare = 0;Klow = 0;Kup = 0;
   rlow = 0;rup = 0;lup = 0;llow = 0; %normIntens=[];
  
end   %end Try

end  %if stdNANcond=0   
 %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
 %Update values if r is better(higher) with removal of early data
try
   if r3>r && stdNANcond==0 
     r = r3; sDat(slps,4)= sDat(slps,4); % rateS 
     K = coeffsArray(1); sDat(slps,2)= coeffsArray(1);  % Carrying Capacity
     l = coeffsArray(2); sDat(slps,3)= coeffsArray(2); %lag time
     
    coeffsArray = coeffvalues(fitObject);
    rmsStg1=errObj.rsquare;
    rmsStg1I(slps)= errObj.rsquare;
    sDat(slps,1)=errObj.rsquare; 
    
     %jh diagnostics*************    
    numFitTpts=nnz((normIntens(:)>=0)==1);
    thresGT2 = find(((normIntens(:)>2)==1), 1);
    thresGT2Tm = filterTimes(thresGT2); 
    numTptsGT2 = nnz((normIntens(:)>=2)==1);
    numTimePts = length(filterTimes);
    
    AUC = (K/r*log(1+exp(-r*(t_end-l)))-K/r*log(exp(-r*(t_end-l)))) - (K/r*log(1+exp(-r*(t_begin-l)))-K/r*log(exp(-r*(t_begin-l))));
    MSR = r3;
    %***************************
    
    rsquare = errObj.rsquare;  
    confObj = confint(fitObject,0.9); % get the 90% confidence 
    NANcond=0;

    confObj_filtered=confObj; 
    Klow = confObj(1,1); sDat(slps,5)=confObj(1,1);
    Kup = confObj(2,1); sDat(slps,6)=confObj(2,1);
    llow = confObj(1,2); sDat(slps,7)=confObj(1,2);
    lup = confObj(2,2); sDat(slps,8)=confObj(2,2);
    rlow = confObj(1,3); sDat(slps,9)=confObj(1,3);
    rup = confObj(2,3); sDat(slps,10)=confObj(2,3);
    if(isnan(Klow)||isnan(Kup)||isnan(llow)||isnan(lup)||isnan(rlow)||isnan(rup))
    NANcond=1;    
    end
    
    filterTimes= Tms; 
    normIntens= Ints;   
     resMat(17)= .00002;  
     resMat(18)= bl;  
     resMat(19)= fitblELr; 
     resMat(20)= minTime; 
    else    % r is better than r3 so use the Std data in the ELr result sheet
     filterTimes=selTimesStd; 
     normIntens=selIntensStd;
     
     lastTptUsed= lastTptUsedStd; %Reinstall Std values for jh diags
     
     Tpt1=filterTimes(1);
     try
     if isempty(Tpt1)
        Tpt1= 0.00002;  %777;
     end
     catch
         Tpt1= 0.00002;  %777;
     end
     
resMat(17)= Tpt1;
     numFitTpts= numFitTptsStd;
     numTptsGT2= numTptsGT2Std;
     thresGT2Tm = thresGT2TmStd;
     
     cutTm(1:4)= 1000;  %-1 means cuts not used or NA
     
     resMat(18)= bl;  %only applicable to Std curve Fit; ELr superceeds and makes meaningless
     resMat(19)= fitblStd; %only applicable to Std curve Fit; ELr superceeds and makes meaningless
     resMat(20)= minTime; %only applicable to Std curve Fit; ELr superceeds and makes meaningless
    end  %if r3>r1
    
   %rup
   %asdf352
% {   
catch ME
    % if no data is given, return zeros
   AUC = 0;MSR = 0;K = 0;r = 0;l = 0;rsquare = 0;Klow = 0;Kup = 0;
   rlow = 0;rup = 0;lup = 0;llow = 0; %normIntens=[];
  
end   %end Try
%}    
resMat(1)= AUC;
resMat(2)= MSR;
resMat(3)= K;
resMat(4)= r;
resMat(5)= l;
resMat(6)= rsquare;
resMat(7)= Klow;
resMat(8)= Kup;
resMat(9)= rlow;
resMat(10)= rup;
resMat(11)= llow;
resMat(12)= lup;
resMat(13)= currSpotArea;
resMat(14)= lastIntensUsed;  %filtNormIntens(length(filtNormIntens));

%spline fit unneccessary and removed therfor no max spline rate time->set 0
    maxRateTime=0;  %ELr will show 0; Std will show [] 
    resMat(15)= maxRateTime;

resMat(16)= lastTptUsed;   %filterTimes(length(filterTimes));


try                   %if Std fit used no cuts .:. no cutTm 
resMat(24)= cutTm(1);
resMat(25)= cutTm(2);
resMat(26)= cutTm(3);
resMat(27)= cutTm(4);
catch
resMat(24)= 999;     %if Std fit used no cuts .:. no cutTm  
resMat(25)= 999;
resMat(26)= 999;
resMat(27)= 999;
end

FiltTimesELr= filterTimes;
NormIntensELr= normIntens;

%**********************************************************************************
%########################################################################## 

  lastTptUsed=max(filterTimes);
  lastIntensUsed=normIntens(length(normIntens));  
  
    if (exist('K','var')&& exist('r','var') && exist('l','var'))
    t = (0:1:200);   
    Growth = K ./ (1 +  exp(-r.* (t - l )));
    fitbl= min(Growth);   %jh diag
    end
  %} 
  try
 if isempty(thresGT2Tm),thresGT2Tm=0;end  %jh diag
  catch
      thresGT2Tm= 0;
      numTptsGT2= 0;
  end

resMat(21)= thresGT2Tm;
resMat(22)= numFitTpts;
resMat(23)= numTptsGT2;
 
 
end  %function end