server-scripts/qhtcp-workflow/apps/matlab/easy/NCscurImCF_3parfor.m

%% CALLED BY NCfitImCFparforFailGbl2.m %%
function [resMatStd, resMat, selTimesStd, selIntensStd, FiltTimesELr, NormIntensELr] =...
  NCscurImCF_3parfor(dataMatrix, AUCfinalTime, currSpotArea, sols, bl, minTime)

  % Preallocate
  resMatStd=zeros(1,27);
  resMat=zeros(1,27);

  % Set internal variables sent to matlab fit function 
  me=200;
  meL=750;
  mi=25;
  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;
  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  
  catch
    % 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

  if (exist('K','var')&& exist('r','var') && exist('l','var'))
    t=(0:1:200);   
    Growth=K ./ (1 +  exp(-r.* (t - l )));
    fitblStd=min(Growth);
  end
   
  cutTm(1:4)=1000;  %-1 means cuts not used or NA
  
  % 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));
  maxRateTime=0; %[];  %Std shows []; ELr shows 0;  %parfor
  resMatStd(15)=0; %maxRateTimestdMeth;
  resMatStd(16)=lastTptUsedStd;
  if isempty(Tpt1Std)
    Tpt1Std=777;
  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; % yhe Standard method has no cuts .:.no cutTm  
  resMatStd(25)=999;
  resMatStd(26)=999;
  resMatStd(27)=999;

  % 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 symmetric 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
      % 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

  % 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
  catch
    % 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
   
  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 % jh diag
    if isempty(thresGT2Tm)
      thresGT2Tm=0
    end
  catch
    thresGT2Tm=0;
    numTptsGT2=0;
  end

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