Merge pull request #51 from ICB-DCM/develop

Develop

@PestoOptions/PestoOptions.m

@@ -153,19 +153,18 @@
         %% get(Parameter|Property)Profiles options        
 
         % The following options are for getParameterProfiles only:
-        % Indices of the parameters for which the profile
-        % is calculated (default = 1:parameters.number).
+
+        % Indices of the parameters for which the profile is calculated. 
+        % Default: profile_optim_index will be set to 1:parameters.number
+        % if both indices are left empty
         parameter_index = [];
+        profile_optim_index = [];
+        profile_integ_index = [];
 
-        % Specifies the method for profile calculation
-        % 'optimization' (default), 'integration', 'mixed'
+        % How should profiles be computed? ('optimization', 'integration',
+        % 'mixed') Default: 'optimization'
         profile_method = 'optimization';
 
-        % Indices which specify which profile whould be optimized and
-        % which should be integrated, if profile_method is set to 'mixed' 
-        % (0: optimization, 1: integration)
-        parameter_method_index = [];
-
         % Indices of the properties for which the profile
         % is to be calculated (default = 1:properties.number).
         property_index = [];
@@ -210,18 +209,18 @@
             'update', 1.25);
 
         % options for Profile integration 
-        solver = struct('type', 'CVODE', ...
+        solver = struct('type', 'ode113', ...
             'algorithm', 'Adams', ...
             'nonlinSolver', 'Newton', ...
             'linSolver', 'Dense', ...
             'gamma', 0, ...
             'eps', 1e-8, ...
-            'minCond', 1e-12, ...
+            'minCond', 1e-10, ...
             'hessian', 'user-supplied', ... 
             'gradient', true, ...
             'MaxStep', 0.1, ...
-            'MinStep', 0.5e-4, ...
-            'MaxNumSteps', 1e7, ...
+            'MinStep', 1e-4, ...
+            'MaxNumSteps', 1e5, ...
             'RelTol', 1e-5, ...
             'AbsTol', 1e-8 ...
             );

@PestoPlottingOptions/PestoPlottingOptions.m

@@ -106,7 +106,7 @@
         % * .sp_ms: marker size for scatter plot (default = 5)
         % * .name: name of legend entry (default = 'S')
         S = struct('plot_type', 0, ...
-            'bins', 'conservative', ...
+            'bins', 'optimal', ...
             'scaling', [], ...
             'hist_col',  [0.7,0.7,0.7], ...
             'sp_col', [0.7,0.7,0.7], ...

README.md

@@ -10,7 +10,7 @@ PESTO features include:
 * Visualization routines
 * and more
 
-These function are demonstrated in several systems biology examples included in the `examples/` directory.
+These function are demonstrated in several systems biology examples included in the [`examples/`](examples/) directory.
 
 ## Installation
 
@@ -26,4 +26,4 @@ See [ ```LICENSE```](LICENSE) file in the PESTO source directory.
 
 ## Documentation
 
-Further PESTO documentation is available in the [ ```/doc```](doc) folder.
+PESTO usage is demonstrated in various [examples](examples/). Further documentation is available in [```doc/PESTO-doc.pdf```](doc/PESTO-doc.pdf).

doc/MatlabDocMaker.m

@@ -466,6 +466,8 @@ function create(varargin)
                         if s ~= 0
                             warn = [warn sprintf('LaTeX compiler output:\n') latexmsg];
                             latexerr = true;
+                        else
+                            copyfile(fullfile(latexdir,'refman.pdf'), fullfile(fileparts(mfilename('fullpath')),'PESTO-doc.pdf'))
                         end
                     else
                         fprintf('cannot compile LaTeX output: no refman.tex found...');

doc/config/Doxyfile.template

@@ -455,7 +455,7 @@ LOOKUP_CACHE_SIZE      = 0
 # normally produced when WARNINGS is set to YES.
 # The default value is: NO.
 
-EXTRACT_ALL            = YES
+EXTRACT_ALL            = NO
 
 # If the EXTRACT_PRIVATE tag is set to YES, all private members of a class will
 # be included in the documentation.
@@ -473,7 +473,7 @@ EXTRACT_PACKAGE        = NO
 # included in the documentation.
 # The default value is: NO.
 
-EXTRACT_STATIC         = YES
+EXTRACT_STATIC         = NO
 
 # If the EXTRACT_LOCAL_CLASSES tag is set to YES, classes (and structs) defined
 # locally in source files will be included in the documentation. If set to NO,
@@ -489,7 +489,7 @@ EXTRACT_LOCAL_CLASSES  = YES
 # included.
 # The default value is: NO.
 
-EXTRACT_LOCAL_METHODS  = YES
+EXTRACT_LOCAL_METHODS  = NO
 
 # If this flag is set to YES, the members of anonymous namespaces will be
 # extracted and appear in the documentation as a namespace called
@@ -862,6 +862,9 @@ RECURSIVE              = YES
 
 EXCLUDE                = doc \
                          DRAM \
+                         private \
+                         examples \
+                         tests \
                          README.md
 
 # The EXCLUDE_SYMLINKS tag can be used to select whether or not files or

examples/Pom1p_gradient_formation/mainPom1.m

@@ -79,7 +79,7 @@
 options.optionsMultistart = PestoOptions();
 options.optionsMultistart.obj_type = 'log-posterior';
 options.optionsMultistart.comp_type = 'sequential';
-options.optionsMultistart.fmincon = optimset(options.optionsMultistart.fmincon,...
+options.optionsMultistart.localOptimizerOptions = optimset(options.optionsMultistart.fmincon,...
    'Algorithm','interior-point',...
    'Display','off',...
    'GradObj','on',...

examples/conversion_reaction/mainConversionReaction.m

@@ -4,8 +4,10 @@
 % * getMultiStarts()
 % * getParameterProfiles()
 % * getParameterSamples()
-% * plotParameterUncertainty()
+% * getParameterConfidenceIntervals()
+% * getPropertyMultiStarts()
 % * getPropertyProfiles()
+% * getPropertySamples()
 % * getPropertyConfidenceIntervals()
 %
 % This example provides a model for the interconversion of two species 
@@ -24,8 +26,8 @@
 % 
 % Profile likelihood calculation is done using getParameterProfiles().
 %
-% Multi-chain Monte-Carlo sampling is performed by getParameterSamples() 
-% and plotted using plotParameterUncertainty().
+% Single-chain Monte-Carlo sampling is performed by getParameterSamples() 
+% and plotted.
 
 
 

examples/enzymatic_catalysis/mainEnzymaticCatalysis.m

@@ -9,7 +9,8 @@
 % Demonstrates furthermore:
 % * how to do sampling without multi-start local optimization beforehand
 % * the value of multi-start local optimization before sampling
-% * how to use Hessian information for optimization
+% * how to use the MEIGO toolbox for optimization
+% * how to compute profile likelihoods via ODE integration
 %
 % This example provides a model for the reaction of a species X_1 to a
 % species X_4, which is catalyzed by an enzyme X_2.
@@ -26,10 +27,9 @@
 % optimization based on these measurements, demonstrating the use of
 % getMultiStarts().
 %
-% The Profile likelihoods are calculated in two different ways: First, the
-% calculation is done by repeated reoptimization using 
-% getParameterProfiles(), then it is done by integrating the an ODE along
-% the profile path using integrateParameterProfiles().
+% The Profile likelihoods are calculated by integrating an ODE following
+% the profile path using getParameterProfiles with the option
+% optionsPesto.profile_method = 'integration'.
 
 
 
@@ -88,9 +88,10 @@
 
 %% Parameter Sampling
 % An adapted Metropolis-Hastings-Algorithm is used to explore the parameter
-% space. Without Multi-start local optimization, this is not extremly
-% effective, but for small problems, this is feasible and PESTO also allows
-% sampling without previous parameter optimization.
+% space. Without Multi-start local optimization, this is not recommended
+% (since it is in genereal ineffective), but for small problems, this is 
+% feasible and PESTO also allows sampling without previous parameter 
+% optimization.
 
 % Length of the chain
 % optionsPesto.MCMC.nsimu_warmup = 2e4;
@@ -124,17 +125,16 @@
 % uncomment:
 
 MeigoOptions = struct(...
-    'maxeval', 500, ...
+    'maxeval', 1000, ...
     'local', struct('solver', 'fmincon', ...
     'finish', 'fmincon', ...
     'iterprint', 0) ...
     );
 
-optionsMeigo = optionsPesto.copy();
-optionsMeigo.localOptimizer = 'meigo-ess';
-optionsMeigo.localOptimizerOptions = MeigoOptions;
-optionsMeigo.n_starts = 1;
-parameters = getMultiStarts(parameters, objectiveFunction, optionsMeigo);
+optionsPesto.localOptimizer = 'meigo-ess';
+optionsPesto.localOptimizerOptions = MeigoOptions;
+optionsPesto.n_starts = 1;
+parameters = getMultiStarts(parameters, objectiveFunction, optionsPesto);
 
 % Options for an alternative multi-start local optimization
 % 
@@ -154,23 +154,8 @@
 % The result of the sampling is compared with profile likelihoods.
 
 optionsPesto.profile_method = 'integration';
-optionsPesto.solver = struct('gamma', 10, ...
-    'type', 'ode113', ...
-    'algorithm', 'Adams', ...
-    'nonlinSolver', 'Newton', ...
-    'linSolver', 'Dense', ...
-    'eps', 1e-8, ...
-    'minCond', 1e-10, ...
-    'gradient', true, ...
-    'hessian', 'user-supplied', ... 
-    'MaxStep', 0.1, ...
-    'MinStep', 1e-5, ...
-    'MaxNumSteps', 1e4, ...
-    'RelTol', 1e-5, ...
-    'AbsTol', 1e-8, ...
-    'hessianStep', 1e-7, ...
-    'gradientStep', 1e-7 ...
-    );
+optionsPesto.solver.gamma = 10;
+optionsPesto.solver.hessian = 'user-supplied';
 
 parameters = getParameterProfiles(parameters, objectiveFunction, optionsPesto);
 

examples/erbb_signaling/mainErbBSignaling.m

@@ -68,8 +68,7 @@
 optionsMultistart.comp_type = 'sequential';
 optionsMultistart.mode      = 'text';
 optionsMultistart.rng       = 0;
-optionsMultistart.fmincon   = optimset(optionsMultistart.fmincon,...
-    'Algorithm','interior-point',...
+optionsMultistart.localOptimizerOptions = optimset('Algorithm','interior-point',...
     'GradObj', 'on',...
     'Display', 'iter', ...
     'MaxIter', 1000,...

examples/jakstat_signaling/mainJakstatSignaling.m

@@ -3,6 +3,9 @@
 % Demonstrates the use of:
 % * getMultiStarts()
 %
+% Demostrates furhtermore
+% * how to implement a user-supplied guess for intial parameters
+%
 % This example provides a model for the JakStat signaling pathway with an
 % time resolved input of the drug EPO. The model has been taken from the
 % papers "Identification of nucleocytoplasmic cycling as a remote sensor in 
@@ -81,21 +84,19 @@
 objectiveFunction = @(theta) logLikelihoodJakstat(theta, amiData);
 
 % PestoOptions
-optionsMultistart           = PestoOptions();
-optionsMultistart.n_starts  = 10;
-optionsMultistart.trace     = true;
-optionsMultistart.mode      = 'visual';
-optionsMultistart.fmincon   = optimset(...
+optionsMultistart          = PestoOptions();
+optionsMultistart.n_starts = 10;
+optionsMultistart.trace    = true;
+optionsMultistart.proposal = 'user-supplied';
+optionsMultistart.obj_type = 'log-posterior';
+optionsMultistart.mode     = 'visual';
+optionsMultistart.localOptimizer = 'fmincon';
+optionsMultistart.localOptimizerOptions = optimset(...
     'Algorithm','interior-point',...
     'GradObj', 'on',...
-    'Display', 'iter-detailed', ...
+    'Display', 'iter', ...
     'MaxIter', 800,...
-    'TolCon', 0,...
     'TolFun', 1e-10,...
-    'TolX', 1e-12,...
-    'MaxProjCGIter', 100,...
-    'TolPCG', 0.001,...
-    'PrecondBandWidth', Inf,...
     'MaxFunEvals', 1000*parameters.number);
 
 %% Perform Multistart optimization

examples/mRNA_transfection/mainTransfection.m

@@ -10,6 +10,11 @@
 % * getPropertySamples()
 % * getPropertyConfidenceIntervals()
 %
+% Demonstrates furthermore:
+% * how to carry out uncertainty analysis for local (non-global) optimum
+% * how to use the PSwarm toolbox for optimization (commented) and what
+%   problems may occur when no gradient based approach is used
+%
 % This example is a model for mRNA transfection, taken from the paper
 % "Single-cell mRNA transfection studies: Delivery, kinetics and statistics 
 % by numbers", by Leonhardt C et al., Nanomedicine: NBM, 2014, vol.10
@@ -22,7 +27,10 @@
 % these measurements, demonstrating the use of getMultiStarts(). The model 
 % fit is then visualized.
 % 
-% Profile likelihood calculation is done using getParameterProfiles().
+% Profile likelihood calculation is done by optimization and integration
+% using getParameterProfiles() with the option 
+% optionsMultistart.profile_method = 'mixed'to have comparison of both
+% methods.
 %
 % Multi-chain Monte-Carlo sampling is performed by getParameterSamples() 
 % and plotted using plotParameterUncertainty().
@@ -159,27 +167,25 @@
 
 % This section uses PSwarm, a particle swarm optimizer
 % (Install from http://www.norg.uminho.pt/aivaz/pswarm/ and uncomment)
+% 
+% optionsMultistartPSwarm = optionsMultistart.copy();
+% optionsMultistartPSwarm.localOptimizer = 'pswarm';
+% optionsMultistartPSwarm.localOptimizerOptions.MaxObj  = 25000;
+% optionsMultistartPSwarm.localOptimizerOptions.MaxIter = 1000;
+% optionsMultistartPSwarm.localOptimizerOptions.Size    = 100;
+% optionsMultistartPSwarm.localOptimizerOptions.Social  = 0.5;
+% optionsMultistartPSwarm.localOptimizerOptions.Cognitial = 0.9;
+% optionsMultistartPSwarm.localOptimizerOptions.IPrint  = -1;
+% 
+% parameters = getMultiStarts(parameters, objectiveFunction, optionsMultistartPSwarm);
 
-optionsMultistartPSwarm = optionsMultistart.copy();
-optionsMultistartPSwarm.localOptimizer = 'pswarm';
-optionsMultistartPSwarm.localOptimizerOptions.MaxObj  = 25000;
-optionsMultistartPSwarm.localOptimizerOptions.MaxIter = 1000;
-optionsMultistartPSwarm.localOptimizerOptions.Size    = 100;
-optionsMultistartPSwarm.localOptimizerOptions.Social  = 0.5;
-optionsMultistartPSwarm.localOptimizerOptions.Cognitial = 0.9;
-optionsMultistartPSwarm.localOptimizerOptions.IPrint  = -1;
-optionsMultistartPSwarm.n_starts = 20;
-
-parameters = getMultiStarts(parameters, objectiveFunction, optionsMultistartPSwarm);
 
-% This is an alternativ section which uses Multi-start local optimization
-% 
-% % Optimization
-% parameters = getMultiStarts(parameters, objectiveFunction, optionsMultistart);
+% This section uses multi-start local optimization
+parameters = getMultiStarts(parameters, objectiveFunction, optionsMultistart);
 
 %% Collection of results, check for bimodality
 
-% Check if a second mode was found
+% Check if a second optimum was found, which is good enough
 for iMode = 2 : optionsMultistart.n_starts
     if (parameters.MS.logPost(iMode) > 39.5)
         if (abs(parameters.MS.par(3,iMode) - parameters.MS.par(3,1)) > 0.1)
@@ -217,11 +223,12 @@
 % by using repeated reoptimization. The information about the profiles is
 % then written to the parameters struct.
 
+% Profiles are computed using optimization by getParameterProfiles()
 parameters = getParameterProfiles(parameters, objectiveFunction, optionsMultistart);
 
 % Computation for the second mode
 optionsMultistart.MAP_index = MAP_index2;
-optionsMultistart.parameter_index = [3, 4];
+optionsMultistart.fh = figure();
 parametersAlt = getParameterProfiles(parametersAlt, objectiveFunction, optionsMultistart);
 
 %% Single-chain Markov chain Monte-Carlo sampling -- Parameters