diff --git a/.gitignore b/.gitignore index b956acc..9664f43 100644 --- a/.gitignore +++ b/.gitignore @@ -48,4 +48,8 @@ WorkSpace/OldProablyTrash/*/* WorkSpace/OldProablyTrash/*/*/* */*/*/tmpPropensityFunctions/* */*/tmpPropensityFunctions/* -*/tmpPropensityFunctions/* \ No newline at end of file +*/tmpPropensityFunctions/* +Examples/Iterative*.mat +WorkSpace/Josh/Fake*.csv +Workspace/*/* +Workspace/* diff --git a/WorkSpace/EricModel/TMPEricMH.mat b/WorkSpace/EricModel/TMPEricMH.mat deleted file mode 100644 index 98d8e6a..0000000 Binary files a/WorkSpace/EricModel/TMPEricMH.mat and /dev/null differ diff --git a/WorkSpace/EricModel/TMPEricMHDusp1.mat b/WorkSpace/EricModel/TMPEricMHDusp1.mat deleted file mode 100644 index 7f7e5f4..0000000 Binary files a/WorkSpace/EricModel/TMPEricMHDusp1.mat and /dev/null differ diff --git a/WorkSpace/EricModel/TMPEricMHGR.mat b/WorkSpace/EricModel/TMPEricMHGR.mat deleted file mode 100644 index 6f2199f..0000000 Binary files a/WorkSpace/EricModel/TMPEricMHGR.mat and /dev/null differ diff --git a/src/+ssit/+fsp/adaptiveFspSolve.m b/src/+ssit/+fsp/adaptiveFspSolve.m index 288b233..139e1a9 100644 --- a/src/+ssit/+fsp/adaptiveFspSolve.m +++ b/src/+ssit/+fsp/adaptiveFspSolve.m @@ -520,7 +520,7 @@ if propensities{i}.isFactorizable y(i) = wt(i).*propensities{i}.stateDependentFactor(x,parameters); else - y(i) = propensities{i}.hybridJointFactor(t,[]); + y(i) = propensities{i}.hybridJointFactor(t,x,parameters,v); end end end diff --git a/src/+ssit/@FspMatrix/FspMatrix.m b/src/+ssit/@FspMatrix/FspMatrix.m index 8aa2a8e..0a3037a 100644 --- a/src/+ssit/@FspMatrix/FspMatrix.m +++ b/src/+ssit/@FspMatrix/FspMatrix.m @@ -166,8 +166,15 @@ w = wt(1)*[obj.terms{1}.matrix*vJ1;... obj.terms{1}.propensity.ODEstoichVector]; for i = 2:length(obj.terms) - w = w + wt(i)*[obj.terms{i}.matrix*vJ1;... - obj.terms{i}.propensity.ODEstoichVector]; + if obj.terms{i}.isFactorizable + w = w + wt(i)*[obj.terms{i}.matrix*vJ1;... + obj.terms{i}.propensity.ODEstoichVector]; + else + tmp = ssit.FspMatrixTerm.generateHybridMatrixTerm(t, obj.terms{i}.propensity, obj.terms{i}.matrix, parameters, ... + obj.terms{i}.numConstraints, v(end-length(upstreamODEs)+1:end)'); + w = w + [tmp*vJ1;... + obj.terms{i}.propensity.ODEstoichVector]; + end end end @@ -310,7 +317,7 @@ % A = A + obj.terms{i}.propensity.hybridFactor(t,v2)*obj.terms{i}.matrix; A = A + wt(i)*obj.terms{i}.matrix; else - A = A + ssit.FspMatrixTerm.generateHybridgMatrixTerm(t, obj.terms{i}.propensity, obj.terms{i}.matrix, parameters, obj.terms{i}.numConstraints, v2); + A = A + ssit.FspMatrixTerm.generateHybridMatrixTerm(t, obj.terms{i}.propensity, obj.terms{i}.matrix, parameters, obj.terms{i}.numConstraints, v2); end end diff --git a/src/+ssit/@Propensity/Propensity.m b/src/+ssit/@Propensity/Propensity.m index 0219cb1..5aec391 100644 --- a/src/+ssit/@Propensity/Propensity.m +++ b/src/+ssit/@Propensity/Propensity.m @@ -112,8 +112,8 @@ else if ~isempty(obj.jointDependentFactor) y = obj.jointDependentFactor(t, x, parameters); - elseif ~isempty(obj.hybridFactor) - y = obj.hybridFactor(t, x, parameters, upstreamSpecies); + elseif ~isempty(obj.hybridJointFactor) + y = obj.hybridJointFactor(t, x, parameters, upstreamSpecies); end end end @@ -244,7 +244,7 @@ oneSym = str2sym('1'); % change to parfor? - parfor iRxn = 1:n_reactions + for iRxn = 1:n_reactions prop_vars = symvar(symbolicExpression{iRxn}); hybridFactor =[]; prefixNameLocal = [prefixName,'_',num2str(iRxn)]; @@ -332,7 +332,7 @@ % are created as much slower anonymous functions. this % will require keeping track since the variables need % to be sent individually to the anonymous functions. - if ~isempty(logicTerms{iRxn})&&(isfield(logicTerms{iRxn},'logT')||isfield(logicTerms{iRxn},'logE')) + if ~isempty(logicTerms{iRxn})&&(isfield(logicTerms{iRxn},'logT')||isfield(logicTerms{iRxn},'logJ')) obj{iRxn}.anonymousT = true; hybridFactor = sym2propfun(expr_t, true, false, nonXTpars(:,1), speciesStoch, varODEs, logicTerms(iRxn)); anyLogical(iRxn) = true; @@ -346,7 +346,7 @@ % hybridFactor = sym2mFun(expr_t, true, false, nonXTpars(:,1), speciesStoch, varODEs); end - if ~isempty(logicTerms{iRxn})&&(isfield(logicTerms{iRxn},'logX')||isfield(logicTerms{iRxn},'logE')) + if ~isempty(logicTerms{iRxn})&&(isfield(logicTerms{iRxn},'logX')||isfield(logicTerms{iRxn},'logJ')) obj{iRxn}.anonymousX = true; anyLogical(iRxn) = true; % elseif sum(contains(string(symvar(expr_x)),'logT')) @@ -420,16 +420,21 @@ for i2=1:length(upstreamODEs) expr_tx = subs(expr_tx,upstreamODEs{i2},varODEs(i2)); end - obj{iRxn}.(jntFactorName) = sym2propfun(expr_tx, true, true, nonXTpars(:,1), speciesStoch, varODEs, logicTerms(iRxn)); + [obj{iRxn}.(jntFactorName),expr_dt_vec_dodei] = ... + sym2propfun(expr_tx, true, true, nonXTpars(:,1), speciesStoch, varODEs, logicTerms(iRxn), true); obj{iRxn}.isTimeDependent = true; expr_t_vec(iRxn) = sym(NaN); expr_x_vec(iRxn) = sym(NaN); - - % if computeSens - % for iPar = 1:n_pars - % expr_tx_vec_sens(iRxn,iPar) = signHybridFactor*diff(expr_tx,nonXTpars{iPar,1}); - % end - % end + + if ~isempty(expr_dt_vec_dodei) + expr_dt_vec_dode(iRxn,:) = expr_dt_vec_dodei; + end + + if computeSens + for iPar = 1:n_pars + expr_tx_vec_sens(iRxn,iPar) = signHybridFactor*diff(expr_tx,nonXTpars{iPar,1}); + end + end end end @@ -560,13 +565,13 @@ n = [0,0,0]; stNew = st; for i=1:3 - J = strfind(st,logTypes{i}); + J = strfind(stNew,logTypes{i}); for j = 1:length(J) - K = strfind(st,'('); + K = strfind(stNew,'('); k1 = max(K(KJ(j))); - logE = st(k1:k2); + logE = stNew(k1:k2); if contains(logE,'t')&&max(contains(logE,species)) n(1)=n(1)+1; logicTerms.logJ{n(1),1} = logE; @@ -704,14 +709,6 @@ end end -for i = 1:length(nonXTpars) - exprStr = strrep(exprStr, nonXTpars{i}, ['Parameters(',num2str(i),')']); -end - -for i = length(nonXTpars):-1:1 - exprStr = strrep(exprStr, ['parameters',num2str(i)], ['Parameters(',num2str(i),')']); -end - for i=1:length(logicTerms) if isfield(logicTerms{i},'logT') for j=1:size(logicTerms{i}.logT,1) @@ -724,15 +721,24 @@ exprStr=strrep(exprStr,logicTerms{i}.logX{j,2},logicTerms{i}.logX{j,1}); end end - if isfield(logicTerms{i},'logE') + if isfield(logicTerms{i},'logJ') % time_dep = true; % state_dep = true; for j=1:size(logicTerms{i}.logX,1) - exprStr=strrep(exprStr,logicTerms.logX{j,2},logicTerms.logX{j,1}); + exprStr=strrep(exprStr,logicTerms{i}.logX{j,2},logicTerms{i}.logX{j,1}); end end end +for i = 1:length(nonXTpars) + exprStr = strrep(exprStr, nonXTpars{i}, ['Parameters(',num2str(i),')']); +end + +for i = length(nonXTpars):-1:1 + exprStr = strrep(exprStr, ['parameters',num2str(i)], ['Parameters(',num2str(i),')']); +end + + if (time_dep && state_dep) fhandle_var = ['@(t, x',parStr,')']; for i = 1:length(varNames) diff --git a/src/CommandLine/SSIT.m b/src/CommandLine/SSIT.m index cc42c20..2780d4a 100644 --- a/src/CommandLine/SSIT.m +++ b/src/CommandLine/SSIT.m @@ -1049,7 +1049,7 @@ function sampleDataFromFSP(obj,fspSoln,saveFile) end function [fimTotal,mleCovEstimate,fimMetrics] = evaluateExperiment(obj,... - fimResults,cellCounts) + fimResults,cellCounts,priorCoVariance) % This function evaluates the provided experiment design (in % "cellCounts" and produces an array of FIMs (one for each % parameter set. @@ -1057,15 +1057,23 @@ function sampleDataFromFSP(obj,fspSoln,saveFile) obj fimResults cellCounts + priorCoVariance = [] end + Ns = size(fimResults,2); Nt = size(fimResults,1); Np = size(fimResults{1,1},1); fimTotal = cell(1,Ns); mleCovEstimate = cell(1,Ns); + if isempty(priorCoVariance) + PriorFIM = zeros(Np); + else + PriorFIM = inv(priorCoVariance); + end + for is=1:Ns - fimTotal{is} = 0*fimResults{1,is}; + fimTotal{is} = PriorFIM; %0*fimResults{1,is}; for it=1:Nt fimTotal{is} = fimTotal{is} + cellCounts(it)*fimResults{it,is}; @@ -1225,10 +1233,13 @@ function sampleDataFromFSP(obj,fspSoln,saveFile) obj.dataSet.linkedSpecies = linkedSpecies; Q = contains(obj.dataSet.dataNames,{'time','Time','TIME'}); - if sum(Q)==1 - obj.dataSet.app.ParEstFitTimesList.Items = {}; + if sum(Q)>=1 obj.dataSet.app.ParEstFitTimesList.Value = {}; - col_time = find(Q); + obj.dataSet.app.ParEstFitTimesList.Items = {}; + if sum(Q)>1 + warning('Provided data more than one entry with keyword "time"') + end + col_time = find(Q,1); obj.dataSet.app.DataLoadingAndFittingTabOutputs.fittingOptions.fit_time_index = col_time; obj.dataSet.app.DataLoadingAndFittingTabOutputs.fittingOptions.fit_times = sort(unique(cell2mat(obj.dataSet.DATA(:,col_time)))); for i=1:length(obj.dataSet.app.DataLoadingAndFittingTabOutputs.fittingOptions.fit_times) @@ -1236,6 +1247,7 @@ function sampleDataFromFSP(obj,fspSoln,saveFile) obj.dataSet.app.ParEstFitTimesList.Value{i} = num2str(obj.dataSet.app.DataLoadingAndFittingTabOutputs.fittingOptions.fit_times(i)); end % We need to make sure that the fitting times are included in the solution times. + else error('Provided data set does not have required column named "time"') end @@ -1717,6 +1729,13 @@ function sampleDataFromFSP(obj,fspSoln,saveFile) fitAlgorithm = 'fminsearch'; end + % parse fitting options + allFitOptions.suppressFSPExpansion = true; + fNames = fieldnames(fitOptions); + for i=1:length(fNames) + allFitOptions.(fNames{i}) = fitOptions.(fNames{i}); + end + if isempty(obj.propensitiesGeneral) obj = formPropensitiesGeneral(obj); end @@ -1735,12 +1754,11 @@ function sampleDataFromFSP(obj,fspSoln,saveFile) if strcmp(obj.solutionScheme,'FSP') % Set solution scheme to FSP. [FSPsoln,bounds] = obj.solve; % Solve the FSP analysis obj.fspOptions.bounds = bounds;% Save bound for faster analyses - objFun = @(x)-obj.computeLikelihood(exp(x),FSPsoln.stateSpace); % We want to MAXIMIZE the likelihood. - - if isfield(fitOptions,'suppressFSPExpansion')&&fitOptions.suppressFSPExpansion + if allFitOptions.suppressFSPExpansion tmpFSPtol = obj.fspOptions.fspTol; obj.fspOptions.fspTol = inf; end + objFun = @(x)-obj.computeLikelihood(exp(x),FSPsoln.stateSpace); % We want to MAXIMIZE the likelihood. elseif strcmp(obj.solutionScheme,'ode') % Set solution scheme to ode. objFun = @(x)-obj.computeLikelihoodODE(exp(x)); % We want to MAXIMIZE the likelihood. end @@ -1749,19 +1767,13 @@ function sampleDataFromFSP(obj,fspSoln,saveFile) switch fitAlgorithm case 'fminsearch' - allFitOptions.suppressFSPExpansion = true; - fNames = fieldnames(fitOptions); - for i=1:length(fNames) - allFitOptions.(fNames{i}) = fitOptions.(fNames{i}); - end - [x0,likelihood,~,otherResults] = fminsearch(objFun,x0,allFitOptions); case 'fminunc' obj.fspOptions.fspTol = inf; objFun = @obj.minusLogL; % We want to MAXIMIZE the likelihood. x0 = log(parGuess); - [x0,likelihood] = fminunc(objFun,x0,fitOptions,FSPsoln.stateSpace,true); + [x0,likelihood] = fminunc(objFun,x0,allFitOptions,FSPsoln.stateSpace,true); case 'particleSwarm' obj.fspOptions.fspTol=inf; @@ -1772,25 +1784,27 @@ function sampleDataFromFSP(obj,fspSoln,saveFile) initSwarm = repmat(x0',fitOptions.SwarmSize-1,1); initSwarm = [x0';initSwarm.*(1+0.1*randn(size(initSwarm)))]; fitOptions.InitialSwarmMatrix = initSwarm; - [x0,likelihood] = particleswarm(OBJps,length(x0),LB,UB,fitOptions); + [x0,likelihood] = particleswarm(OBJps,length(x0),LB,UB,allFitOptions); case 'mlSearch' % Not yet working efficiently. - allFitOptions.maxIter=1000; - allFitOptions.burnIn=30; - allFitOptions.updateRate=10; - allFitOptions.guessRate=1000; - allFitOptions.proposalDistribution=@(x)x+0.01*randn(size(x)); - allFitOptions.useFIMforSearch = false; - allFitOptions.CovFIMscale = 0.6; - allFitOptions.suppressFSPExpansion = true; - allFitOptions.logForm = true; - allFitOptions.plotFunVals = false; - allFitOptions.proposalDistributionWide=@(x)x+randn(size(x)); - - fNames = fieldnames(fitOptions); + defaultFitOptions.maxIter=1000; + defaultFitOptions.burnIn=30; + defaultFitOptions.updateRate=10; + defaultFitOptions.guessRate=1000; + defaultFitOptions.proposalDistribution=@(x)x+0.01*randn(size(x)); + defaultFitOptions.useFIMforSearch = false; + defaultFitOptions.CovFIMscale = 0.6; + defaultFitOptions.suppressFSPExpansion = true; + defaultFitOptions.logForm = true; + defaultFitOptions.plotFunVals = false; + defaultFitOptions.proposalDistributionWide=@(x)x+randn(size(x)); + + fNames = fieldnames(defaultFitOptions); for i=1:length(fNames) - allFitOptions.(fNames{i}) = fitOptions.(fNames{i}); + if ~isfield(allFitOptions,fNames{i}) + allFitOptions.(fNames{i}) = defaultFitOptions.(fNames{i}); + end end if allFitOptions.logForm @@ -1805,26 +1819,29 @@ function sampleDataFromFSP(obj,fspSoln,saveFile) case 'MetropolisHastings' - allFitOptions.isPropDistSymmetric=true; - allFitOptions.thin=1; - allFitOptions.numberOfSamples=1000; - allFitOptions.burnIn=100; - allFitOptions.progress=true; - allFitOptions.proposalDistribution=@(x)x+0.1*randn(size(x)); - allFitOptions.numChains = 1; - allFitOptions.useFIMforMetHast = false; - allFitOptions.CovFIMscale = 0.6; - allFitOptions.suppressFSPExpansion = true; - allFitOptions.logForm = true; + defaultFitOptions.isPropDistSymmetric=true; + defaultFitOptions.thin=1; + defaultFitOptions.numberOfSamples=1000; + defaultFitOptions.burnIn=100; + defaultFitOptions.progress=true; + defaultFitOptions.proposalDistribution=@(x)x+0.1*randn(size(x)); + defaultFitOptions.numChains = 1; + defaultFitOptions.useFIMforMetHast = false; + defaultFitOptions.CovFIMscale = 0.6; + defaultFitOptions.suppressFSPExpansion = true; + defaultFitOptions.logForm = true; j=1; while exist(['TMPmh_',num2str(j),'.mat'],'file') j=j+1; end - allFitOptions.saveFile = ['TMPmh_',num2str(j),'.mat']; - fNames = fieldnames(fitOptions); + defaultFitOptions.saveFile = ['TMPmh_',num2str(j),'.mat']; + + fNames = fieldnames(defaultFitOptions); for i=1:length(fNames) - allFitOptions.(fNames{i}) = fitOptions.(fNames{i}); + if ~isfield(allFitOptions,fNames{i}) + allFitOptions.(fNames{i}) = defaultFitOptions.(fNames{i}); + end end if allFitOptions.logForm @@ -1909,7 +1926,7 @@ function sampleDataFromFSP(obj,fspSoln,saveFile) pars = exp(x0); - if strcmp(obj.solutionScheme,'FSP')&&isfield(fitOptions,'suppressFSPExpansion') + if strcmp(obj.solutionScheme,'FSP')&&allFitOptions.suppressFSPExpansion obj.fspOptions.fspTol = tmpFSPtol; end diff --git a/src/CommandLine/SSITMultiModel.m b/src/CommandLine/SSITMultiModel.m index 200a5c2..298a747 100644 --- a/src/CommandLine/SSITMultiModel.m +++ b/src/CommandLine/SSITMultiModel.m @@ -141,8 +141,13 @@ % all models for the provided parameter combination. Nmods = length(SMM.logLikelihoodFunctions); logLs = zeros(1,Nmods); + + G = cell(1,Nmods); + for i = 1:Nmods + G{i} = SMM.logLikelihoodFunctions{i}; + end for i = 1:Nmods - logLs(i) = SMM.logLikelihoodFunctions{i}(parameterGuess); + logLs(i) = G{i}(parameterGuess); end totalLogLikelihood = sum(logLs); @@ -267,7 +272,7 @@ if allFitOptions.logForm&&min(eig(FIMlog))<1 disp('Warning -- FIM has one or more small eigenvalues. Reducing proposal width to 10x in those directions. MH Convergence may be slow.') - FIMlog = FIMlog + 1*eye(length(FIMlog)); + FIMlog = FIMlog + 0.1*eye(length(FIMlog)); end covLog = FIMlog^-1; diff --git a/src/gui/data_loading/makeSeparatePlotOfData.m b/src/gui/data_loading/makeSeparatePlotOfData.m index de2a0b0..ad02fd1 100644 --- a/src/gui/data_loading/makeSeparatePlotOfData.m +++ b/src/gui/data_loading/makeSeparatePlotOfData.m @@ -7,193 +7,242 @@ function makeSeparatePlotOfData(app,smoothWindow,fignums,usePanels) end %% This function creates a histogram from the loaded data. -NT = length(app.ParEstFitTimesList.Value); -% NdMod = max(1,length(size(app.DataLoadingAndFittingTabOutputs.fitResults.current))-1); -NdMod = size(app.NameTable.Data,1); +numTimes = length(app.ParEstFitTimesList.Value); +NdModFit = max(1,length(size(app.DataLoadingAndFittingTabOutputs.fitResults.current))-1); NdDat = length(app.SpeciesForFitPlot.Value); +if NdModFit~=NdDat + error('Model and Data size do not match in fitting routine') +end -Plts_to_make = zeros(1,NdMod); +% Plts_to_make = zeros(1,NdModFit); +% Density and Cumulative Distributions plots for DistType = 0:1 - if isempty(fignums) - figure - else - figure(fignums(DistType+1)) + % Choose which figures to plot in + if isempty(fignums); figHandle = figure; + else; figHandle = figure(fignums(DistType+1)); end + if DistType==0; set(figHandle,'Name','Cumulative Distributions versus Time') + elseif DistType==1; set(figHandle,'Name','Probability Distributions versus Time') end - if NT<=4 - subPlotSize = [1,NT]; + + % initialize histogram plots + if numTimes<=4 + subPlotSize = [1,numTimes]; else - subPlotSize(2) = ceil(sqrt(NT)); - subPlotSize(1) = ceil(NT/subPlotSize(2)); + subPlotSize(2) = ceil(sqrt(numTimes)); + subPlotSize(1) = ceil(numTimes/subPlotSize(2)); end + modelHistTime = cell(numTimes,NdDat); + dataHistTime = cell(numTimes,NdDat); - for iTime = 1:NT - it = find(strcmp(app.ParEstFitTimesList.Value,app.ParEstFitTimesList.Value{iTime})); + % loop over number of time points + for iTime = 1:numTimes + % Find index of time in model solution times set + % iTime = find(strcmp(app.ParEstFitTimesList.Value,app.ParEstFitTimesList.Value{iTime})); if usePanels subplot(subPlotSize(1),subPlotSize(2),iTime) end - FNHists = gca; -% hold off + FigHandleHists = gca; - L = {}; + LegNames = {}; xm = 0; ym = 0; - % Switches plotted data, (x1, x2, x3,...), and legend depending on if the + % Switches plotted data, (e.g., x1, x2, x3,...), and legend depending on if the % species was selected for plotting. for icb=1:NdDat - cb = contains(app.SpeciesForFitPlot.Items,app.SpeciesForFitPlot.Value{icb}); - Plts_to_make(cb) = 1; - - soMod = setdiff([1:NdMod],find(cb~=0)); - soDat = setdiff([1:NdDat],icb); - le = app.NameTable.Data{icb,2}; - -% if cb - matTensor = double(app.DataLoadingAndFittingTabOutputs.dataTensor); + cb = max(contains(app.SpeciesForFitPlot.Items,app.SpeciesForFitPlot.Value{icb})); + % Plts_to_make(cb) = 1; + + indsToSumOver = setdiff([1:NdModFit],icb); + % soDat = setdiff([1:NdDat],icb); + speciesName = app.NameTable.Data{icb,2}; + + %% Histograms for the data. + if cb % If this species selected for plotting + dataTensor = double(app.DataLoadingAndFittingTabOutputs.dataTensor); if length(size(app.DataLoadingAndFittingTabOutputs.dataTensor))==1 - H1 = matTensor; + H1 = dataTensor; else - try - H1 = squeeze(matTensor(it,:,:,:,:,:,:,:,:,:)); - catch - whos + try H1 = squeeze(dataTensor(iTime,:,:,:,:,:,:,:,:,:)); + catch; whos; end end - if ~isempty(soDat) - H1 = sum(H1,soDat); + if ~isempty(indsToSumOver) + H1 = sum(H1,indsToSumOver); end - H1 = H1/sum(H1(:)); % Normalize data before plotting. - + H1 = H1(:)/sum(H1(:)); % Normalize data before plotting. + dataHistTime{iTime,icb} = H1; + + % make data histogram plots if DistType -% stairs(FNHists,[0:length(H1)-1],smoothBins(H1,smoothWindow)); hold on - stairs(FNHists,[0:length(H1)-1],smoothBins(H1,smoothWindow),'linewidth',3); - hold(FNHists,'on') + stairs(FigHandleHists,[0:length(H1)-1],smoothBins(H1,smoothWindow),'linewidth',3); + hold(FigHandleHists,'on') ym = max(ym,max(H1)); else - stairs(FNHists,[0:length(H1)-1],cumsum(H1),'linewidth',3); - hold(FNHists,'on') + stairs(FigHandleHists,[0:length(H1)-1],cumsum(H1),'linewidth',3); + hold(FigHandleHists,'on') ym=1; end - L{end+1} = [le,'-data']; + + % add name to legend + LegNames{end+1} = [speciesName,'-data']; xm = max(xm,length(H1)); - - % if ~isempty(soMod) - % M = squeeze(app.DataLoadingAndFittingTabOutputs.fitResults.current(it,:,:,:,:,:,:,:,:,:,:,:,:)); - % H1 = squeeze(sum(M,soMod)); - % else - H1 = squeeze(app.DataLoadingAndFittingTabOutputs.fitResults.current(it,:,:,:,:,:,:,:,:,:,:,:,:)); - % end - H1(end+1)=0; + + %% Histograms for the model + if ~isempty(indsToSumOver) + modelTensor = squeeze(app.DataLoadingAndFittingTabOutputs.fitResults.current(iTime,:,:,:,:,:,:,:,:,:,:,:,:)); + H1 = squeeze(sum(modelTensor,indsToSumOver)); + else + H1 = squeeze(app.DataLoadingAndFittingTabOutputs.fitResults.current(iTime,:,:,:,:,:,:,:,:,:,:,:,:)); + end + modelHistTime{iTime,icb} = H1; + + % make model histogram plots if DistType - stairs(FNHists,[0:length(H1)-1],smoothBins(H1,smoothWindow),'linewidth',3); - hold(FNHists,'on') + stairs(FigHandleHists,[0:length(H1)-1],smoothBins(H1,smoothWindow),'linewidth',3); + hold(FigHandleHists,'on') ym = max(ym,max(H1)); else - stairs(FNHists,[0:length(H1)-1],cumsum(H1),'linewidth',3); - hold(FNHists,'on') + stairs(FigHandleHists,[0:length(H1)-1],cumsum(H1),'linewidth',3); + hold(FigHandleHists,'on') ym=1; end - L{end+1} = [le,'-mod']; -% end - end - title(['t = ',app.ParEstFitTimesList.Value{iTime}]) -% if mod(iTime,subPlotSize(2))==1 - ylabel('Probability') -% ylim = get(gca,'ylim'); -% else -% set(gca,'yticklabels',{}) -% set(gca,'ylim',ylim); -% end -% if iTime>NT-subPlotSize(2) - xlabel('Mol. Count') -% else -% set(gca,'xticklabels',{}) -% end -set(gca,'fontsize',15) - end -end + % add name to legend + LegNames{end+1} = [speciesName,'-mod']; + end -%% Make a trajectory plot for model. -NdMod = size(app.NameTable.Data,1); -if isempty(fignums) - figure -else - figure(fignums(3)) -end -FNTraj = gca; -T_array = eval(app.FspPrintTimesField.Value); -for it = 1:length(T_array) - if ~isempty(app.FspTabOutputs.solutions{it}) - for i=1:NdMod - INDS = setdiff([1:NdMod],i); - - % Add effect of PDO. - px = app.FspTabOutputs.solutions{it}.p; - if ~isempty(app.FIMTabOutputs.distortionOperator) - px = app.FIMTabOutputs.distortionOperator.computeObservationDist(px); + %% Adjust histograms plot style + title(['t = ',app.ParEstFitTimesList.Value{iTime}]) + if mod(iTime,subPlotSize(2))==1 + ylabel('Probability') + % ylim = get(gca,'ylim'); + % else + % set(gca,'yticklabels',{}) + % set(gca,'ylim',ylim); end - if ~isempty(INDS) - mdist{i} = double(px.sumOver(INDS).data); - else - mdist{i} = double(px.data); + if iTime>numTimes-subPlotSize(2) + xlabel('Num. Mol.') + % else + % set(gca,'xticklabels',{}) end + set(gca,'fontsize',15) end - for j=1:NdMod - mns(it,j) = [0:length(mdist{j})-1]*mdist{j}; - mns2(it,j) = [0:length(mdist{j})-1].^2*mdist{j}; - end + end + legend(LegNames) +end + +%% Make trajectory plots for model. +% NdModAll = size(app.NameTable.Data,1); +if isempty(fignums); figHandle = figure; +else; figHandle = figure(fignums(3)); +end +set(figHandle,'Name','Mean +/- STD vs. time') + +meanVarTrajAxis = gca; + +tArrayModel = eval(app.FspPrintTimesField.Value); +% if length(tArrayModel)==numTimes +% tArrayPlot = tArrayModel; +% else +% tArrayPlot = app.DataLoadingAndFittingTabOutputs.fittingOptions.dataTimes; +% end + +for iTime = 1:length(tArrayModel) + for j=1:NdDat + soInds = find(~strcmp(app.SpeciesForFitPlot.Items,app.SpeciesForFitPlot.Value{j})); + Z = double(app.FspTabOutputs.solutions{iTime}.p.sumOver(soInds).data); + mnsMod(iTime,j) = [0:length(Z)-1]*Z; + mns2Mod(iTime,j) = [0:length(Z)-1].^2*Z; end end -vars = mns2-mns.^2; + + % for iTimeModel = 1:length(T_array) + % if ~isempty(app.FspTabOutputs.solutions{iTimeModel}) + % for i=1:NdModAll + % INDS = setdiff([1:NdDat],i); + % + % % Add effect of PDO. + % px = app.FspTabOutputs.solutions{iTimeModel}.p; + % if ~isempty(app.FIMTabOutputs.distortionOperator) + % px = app.FIMTabOutputs.distortionOperator.computeObservationDist(px); + % end + % if ~isempty(INDS) + % mdist{i} = double(px.sumOver(INDS).data); + % else + % mdist{i} = double(px.data); + % end + % end + % for j=1:NdModAll + % mns(iTimeModel,j) = [0:length(mdist{j})-1]*mdist{j}; + % mns2(iTimeModel,j) = [0:length(mdist{j})-1].^2*mdist{j}; + % end + % end + % end +% end +vars = mns2Mod-mnsMod.^2; cols = ['b','r','g','m','c','k']; cols2 = [.90 .90 1.00; 1.00 .90 .90; .90 1.00 .90; .60 .60 1.00; 1.00 .60 .60; .60 1.00 .60]; LG = {}; -for iplt=1:NdMod - if Plts_to_make(iplt) - BD = [mns(:,iplt)'+sqrt(vars(:,iplt)'),mns(end:-1:1,iplt)'-sqrt(vars(end:-1:1,iplt)')]; - TT = [T_array(1:end),T_array(end:-1:1)]; - fill(FNTraj,TT,BD,cols2(iplt,:)); - hold(FNTraj,'on'); - plot(FNTraj,T_array,mns(:,iplt),cols(iplt),'linewidth',2); +for iplt=1:NdDat + % if Plts_to_make(iplt) + BD = [mnsMod(:,iplt)'+sqrt(vars(:,iplt)'),mnsMod(end:-1:1,iplt)'-sqrt(vars(end:-1:1,iplt)')]; + TT = [tArrayModel(1:end),tArrayModel(end:-1:1)]; + fill(meanVarTrajAxis,TT,BD,cols2(iplt,:)); + hold(meanVarTrajAxis,'on'); + plot(meanVarTrajAxis,tArrayModel,mnsMod(:,iplt),cols(iplt),'linewidth',2); LG{end+1} = [char(app.NameTable.Data(iplt,2)),' Model Mean \pm std']; LG{end+1} = [char(app.NameTable.Data(iplt,2)),' Model Mean']; - end + % end end title('Trajectory of means and standard deviations') %% Add data to trajectory plot. -Ned = length(size(app.DataLoadingAndFittingTabOutputs.dataTensor))-1; -for iplt = 1:Ned -% if Plts_to_make(iplt) - so = setdiff([1:Ned],iplt)+1; - matTensor = double(app.DataLoadingAndFittingTabOutputs.dataTensor); - if ~isempty(so) - matTensor = squeeze(sum(matTensor,so)); - end - P = matTensor./sum(matTensor,2); - mn = P*[0:size(matTensor,2)-1]'; - mn2 = P*([0:size(matTensor,2)-1]').^2; - var = mn2-mn.^2; -% plot(FNTraj,app.DataLoadingAndFittingTabOutputs.fittingOptions.dataTimes,mn,[cols(iplt),... -% 'o'],'MarkerSize',12,'MarkerFaceColor',cols(iplt)) - errorbar(FNTraj,app.DataLoadingAndFittingTabOutputs.fittingOptions.dataTimes,... - mn,sqrt(var),[cols(iplt),... - 'o'],'MarkerSize',12,'MarkerFaceColor',cols(iplt),'linewidth',3) - LG{end+1} = [char(app.NameTable.Data(iplt,2)),' Data mean \pm std']; -% end +for iTime = 1:numTimes + for j=1:NdDat + mnsDat(iTime,j) = [0:length(dataHistTime{iTime,j})-1]*dataHistTime{iTime,j}; + mns2Dat(iTime,j) = [0:length(dataHistTime{iTime,j})-1].^2*dataHistTime{iTime,j}; + end +end +varDat = mns2Dat-mnsDat.^2; +T_array = app.DataLoadingAndFittingTabOutputs.fittingOptions.dataTimes; +for j=1:NdDat + errorbar(meanVarTrajAxis,T_array,... + mnsDat(:,j),sqrt(varDat(:,j)),[cols(j),... + 'o'],'MarkerSize',12,'MarkerFaceColor',cols(j),'linewidth',3) + LG{end+1} = [char(app.NameTable.Data(j,2)),' Data mean \pm std']; end -legend(FNTraj,LG) + +% Ned = length(size(app.DataLoadingAndFittingTabOutputs.dataTensor))-1; +% for iplt = 1:Ned +% % if Plts_to_make(iplt) +% so = setdiff([1:Ned],iplt)+1; +% dataTensor = double(app.DataLoadingAndFittingTabOutputs.dataTensor); +% if ~isempty(so) +% dataTensor = squeeze(sum(dataTensor,so)); +% end +% P = dataTensor./sum(dataTensor,2); +% mn = P*[0:size(dataTensor,2)-1]'; +% mn2 = P*([0:size(dataTensor,2)-1]').^2; +% var = mn2-mn.^2; +% errorbar(meanVarTrajAxis,app.DataLoadingAndFittingTabOutputs.fittingOptions.dataTimes,... +% mn,sqrt(var),[cols(iplt),... +% 'o'],'MarkerSize',12,'MarkerFaceColor',cols(iplt),'linewidth',3) +% LG{end+1} = [char(app.NameTable.Data(iplt,2)),' Data mean \pm std']; +% % end +% end +legend(meanVarTrajAxis,LG) xlabel('Time') ylabel('Response') -set(FNTraj,'fontsize',15) +set(meanVarTrajAxis,'fontsize',15) if max(T_array)>0 - set(FNTraj,'xlim',[0,max(T_array)]) + set(meanVarTrajAxis,'xlim',[0,max(T_array)]) end -%% + +%% Make plots of Likelihood functions versus time. if isempty(fignums) figure else @@ -223,6 +272,29 @@ function makeSeparatePlotOfData(app,smoothWindow,fignums,usePanels) set(gca,'fontsize',15); end +%% Make Joint Density Plots if 2 or more species in data set +% if NdDat==2 +% ModelTensor = app.DataLoadingAndFittingTabOutputs.fitResults.current; +% DataTensor = app.DataLoadingAndFittingTabOutputs.fitResults.currentData; +% figHandle = figure; +% for iTime = 1:size(ModelTensor,1) +% Z = log10(squeeze(ModelTensor(iTime,:,:))+1e-6); +% subplot(subPlotSize(1),subPlotSize(2),iTime) +% contourf(Z) +% end +% figHandle = figure; +% for iTime = 1:size(ModelTensor,1) +% Z = log10(squeeze(DataTensor(iTime,:,:))+1e-6); +% subplot(subPlotSize(1),subPlotSize(2),iTime) +% contourf(Z) +% end +% end + + + + + + end diff --git a/tests/poissonTest.m b/tests/poissonTest.m index 408133b..8e14452 100644 --- a/tests/poissonTest.m +++ b/tests/poissonTest.m @@ -275,6 +275,25 @@ function TestFIM(testCase) testCase.verifyEqual(diff<0.001, true, ... 'FIM Calculation is not within 1e-4% Tolerance'); + end + + function TestFIMwPrior(testCase) + % In this test case, we check that the FIM calculation using + % the FSP matches to the analytical expresion for the Poisson + % model when there is a appropriate Prior: + Model = testCase.Poiss; + Model.solutionScheme = 'fspSens'; + Model.fspOptions.fspTol = 1e-6; + SensSoln = Model.solve; + fspFIM = Model.computeFIM(SensSoln.sens); + t = Model.tSpan; + SIGprior = diag(rand(1,2)); + FIMwPrior = Model.evaluateExperiment(fspFIM, [1:length(t)],... + SIGprior); + FIMwoPrior = Model.evaluateExperiment(fspFIM,[1:length(t)]); + diff = max(abs(FIMwPrior{1}-FIMwoPrior{1}-inv(SIGprior)),[],"all"); + testCase.verifyEqual(diff<1e-5, true, ... + 'FIM Prior Calculation is not within 1e-4% Tolerance'); end