diff --git a/workflow/templates/easy/NCfitImCFparforFailGbl2.m b/workflow/templates/easy/NCfitImCFparforFailGbl2.m index ee2ee97e..262ce053 100755 --- a/workflow/templates/easy/NCfitImCFparforFailGbl2.m +++ b/workflow/templates/easy/NCfitImCFparforFailGbl2.m @@ -218,11 +218,8 @@ function [par4scanselIntensStd,par4scanselTimesStd,par4scanTimesELr,par4scanInte outCstd(ii,:)=resMatStd; %{ii, par4resMatStd}; end - %*********************************19_1001*********************************** - %To accomodate parfor copy par4scan thru global p4 functions inside of - %parfor loop --then outside to par4Gbl_Main8b.m - - %************************************************************************** + % To accomodate parfor + % Copy par4scan thru global p4 functions inside of parfor loop --then outside to par4Gbl_Main8b.m fileExt='.txt'; filePrefix='FitResultsComplete_'; fileNamePlate=[filePrefix fileSuffix fileExt]; diff --git a/workflow/templates/easy/NCscurImCF_3parfor.m b/workflow/templates/easy/NCscurImCF_3parfor.m index ce9b241d..608cb46e 100755 --- a/workflow/templates/easy/NCscurImCF_3parfor.m +++ b/workflow/templates/easy/NCscurImCF_3parfor.m @@ -1,437 +1,356 @@ %% 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; + 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) + % 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; + ||(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; + end + filterTimes=filterTimes'; + selTimesStd=filterTimes; + normIntens=normIntens'; + selIntensStd=normIntens; + lastTptUsed=1; + lastIntensUsed=1; thresGT2TmStd=0; - try - + try lastTptUsed=max(filterTimes); lastIntensUsed=normIntens(length(normIntens)); - lastIntensUsedStd= lastIntensUsed; - - lastTptUsedStd= lastTptUsed; - Tpt1Std= filterTimes(1); + lastIntensUsedStd=lastIntensUsed; + lastTptUsedStd=lastTptUsed; + Tpt1Std=filterTimes(1); numFitTptsStd=nnz((normIntens(:)>=0)==1); - thresGT2 = find(((normIntens(:)>2)==1), 1); + thresGT2=find(((normIntens(:)>2)==1), 1); if isempty(thresGT2) thresGT2TmStd=0; else - thresGT2TmStd = filterTimes(thresGT2); + 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); + 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; + 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 + 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 + % 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); + 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 + 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 + 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; + 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 -cutTm(4)= stdTimes(tc2LdatPt-1); -cutDatNum(4)= tc2LdatPt-2; %length(stdTimes(end-(lengthKlast-1):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 -end - - -%************************************************ -Ints=[]; -Tms=[]; -Ints= ints'; -Tms= tms'; - -try - filterTimes= Tms; filterTimes4= Tms; - normIntens= Ints; normIntens4=Ints; + cutTm(1:4)=1000; %-1 means cuts not used or NA -%---------------------------------------------------------------------------------------------- - - %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 + % 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)); -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; + maxRateTime=0; %[]; %Std shows []; ELr shows 0; %parfor + resMatStd(15)=0; %maxRateTimestdMeth; - 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; + 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; %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 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]); - 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 + 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 - %rup - %asdf352 -% { -catch ME + 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 + + % 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 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; - -%********************************************************************************** -%########################################################################## + 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 + 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 - %} - try - if isempty(thresGT2Tm),thresGT2Tm=0;end %jh diag catch - thresGT2Tm= 0; - numTptsGT2= 0; + thresGT2Tm=0; + numTptsGT2=0; end -resMat(21)= thresGT2Tm; -resMat(22)= numFitTpts; -resMat(23)= numTptsGT2; - - -end %function end + resMat(21)=thresGT2Tm; + resMat(22)=numFitTpts; + resMat(23)=numTptsGT2; +end \ No newline at end of file