Merge pull request #84 from ICB-DCM/develop

Small fixes

@PestoOptions/PestoOptions.m

@@ -172,7 +172,7 @@
         % How should profiles be computed (if no more precise options are 
         % set like profile_optim_index or profile_integ_index)? 
         % Possibilities: {'optimization' (default), 'integration', 'mixed'}
-        profile_method = 'optimization';
+        profile_method = 'default';
 
 
 

examples/enzymatic_catalysis/mainEnzymaticCatalysis.m

@@ -9,10 +9,12 @@
 % Demonstrates furthermore:
 % * how to do sampling without multi-start local optimization beforehand
 % * how to use the MEIGO toolbox for optimization (commented)
-% * how to compute profile likelihoods via ODE integration
+% * how to compute profile likelihoods via the integration method
 % * how to use plotting functions outside the get... routines
 % * the reliability of sampling and profiling in the case of
 %   non-identifiabilites
+% * how to use diagnosis tool (e.g. plotMCMCDiagnosis and
+%   plotMultiStartDiagnosis)
 %
 % 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.
@@ -29,6 +31,11 @@
 % optimization based on these measurements, demonstrating the use of
 % getMultiStarts().
 %
+% Parameter sampling is done first without prior information from
+% optimization, then with information from optimization.
+%
+% Parameter optimization is done using multi-start local optimization.
+%
 % The Profile likelihoods are calculated by integrating an ODE following
 % the profile path using getParameterProfiles with the option
 % optionsPesto.profile_method = 'integration'.
@@ -99,8 +106,8 @@
 
 % PT (with only 1 chain -> AM) specific options:
 optionsSampling.samplingAlgorithm = 'PT';
-optionsSampling.PT.nTemps         = 5;
-optionsSampling.PT.exponentT      = 4;    
+optionsSampling.PT.nTemps         = 6;
+optionsSampling.PT.exponentT      = 6;    
 optionsSampling.PT.regFactor      = 1e-8;
 optionsSampling.PT.temperatureAdaptionScheme = 'Lacki15'; %'Vousden16';
 
@@ -113,6 +120,8 @@
 % Run the sampling
 parameters = getParameterSamples(parameters, objectiveFunction, optionsSampling);
 
+% Use a diagnosis tool to see, how plotting worked (see burn-in etc.)
+plotMCMCdiagnosis(parameters, 'parameters');
 
 %% Calculate Confidence Intervals
 % Confidence Intervals for the Parameters are inferred from the local 
@@ -151,6 +160,9 @@
 optionsPesto.n_starts = 5;
 parameters = getMultiStarts(parameters, objectiveFunction, optionsPesto);
 
+% Use a diagnosis tool to see, how optimization worked
+plotMultiStartDiagnosis(parameters);
+
 
 %% Calculate Confidence Intervals
 % Confidence Intervals for the Parameters are inferred from the local 

getParameterConfidenceIntervals.m

@@ -108,7 +108,7 @@
 
         % Confidence intervals computed using sample
         if isfield(parameters,'S')
-            parameters.CI.S(i,:,k) = prctile(parameters.S.par(i,:),50 + 100*[-alpha(k)/2, alpha(k)/2]);
+            parameters.CI.S(i,:,k) = prctile(parameters.S.par(i,:,1),50 + 100*[-alpha(k)/2, alpha(k)/2]);
         end
     end
 end

getParameterProfiles.m

@@ -90,18 +90,33 @@
 end
 
 % Process, which profiles should be computed in which manner
+if strcmp(options.profile_method, 'default')
+    if (isempty(options.profile_optim_index) && isempty(options.profile_integ_index))
+        options.profile_method = 'optimization';
+    elseif (~isempty(options.profile_optim_index) && isempty(options.profile_integ_index))
+        options.profile_method = 'optimization';
+    elseif (isempty(options.profile_optim_index) && ~isempty(options.profile_integ_index))
+        options.profile_method = 'integration';
+    elseif (~isempty(options.profile_optim_index) && ~isempty(options.profile_integ_index))
+        options.profile_method = 'mixed';
+    end
+end    
 if isempty(options.parameter_index)
     switch options.profile_method
         case 'optimization'
-            options.profile_optim_index = 1:parameters.number;
+            if isempty(options.profile_optim_index)
+                options.profile_optim_index = 1:parameters.number;
+            end
             if ~isempty(options.profile_integ_index)
-                options.profile_optim_index(options.profile_integ_index) = [];
+                error('Some profiles seem to be computed twice. Please redefine consistent options!');
             end
 
         case 'integration'
-            options.profile_integ_index = 1:parameters.number;
+            if isempty(options.profile_integ_index)
+                options.profile_integ_index = 1:parameters.number;
+            end
             if ~isempty(options.profile_optim_index)
-                options.profile_integ_index(options.profile_optim_index) = [];
+                error('Some profiles seem to be computed twice. Please redefine consistent options!');
             end
 
         case 'mixed'
@@ -113,27 +128,34 @@
             if length(unique([options.profile_optim_index options.profile_integ_index])) < length([options.profile_optim_index options.profile_integ_index])
                 error('Some profiles seem to be computed twice. Please redefine consistent options!');
             end
+
+        otherwise
+            error('Unknown profile computationg method');
     end
     options.parameter_index = sort(unique([options.profile_optim_index options.profile_integ_index]));
 else
     switch options.profile_method
         case 'optimization'
             options.profile_optim_index = options.parameter_index;
             if ~isempty(options.profile_integ_index)
-                options.profile_optim_index(options.profile_integ_index) = [];
+                error('Some profiles seem to be computed twice. Please redefine consistent options!');
             end
 
         case 'integration'
             options.profile_integ_index = options.parameter_index;
             if ~isempty(options.profile_optim_index)
-                options.profile_integ_index(options.profile_optim_index) = [];
+                error('Some profiles seem to be computed twice. Please redefine consistent options!');
             end
 
         case 'mixed'
-            if (length(unique([options.profile_optim_index, options.profile_integ_index])) < length(unique([options.parameter_index, options.profile_optim_index, options.profile_integ_index])))
+            if (length(unique([options.profile_optim_index, options.profile_integ_index])) ~= length([options.profile_optim_index, options.profile_integ_index])) ...
+                    || (length([options.profile_optim_index, options.profile_integ_index]) ~= length(options.parameter_index))
                 error('Inconsistent settings for indices in profile calculation.');
             end
             options.parameter_index = sort(unique([options.profile_optim_index options.profile_integ_index]));
+
+        otherwise
+            error('Unknown profile computationg method');
     end
 end
 

plotMCMCdiagnosis.m

@@ -35,7 +35,7 @@
 if nargin >= 2 && ~isempty(varargin{1})
     type = varargin{1};
     if ~max(strcmp({'parameters','log-posterior'},type))
-        error('The ''type'' of plot is unknown. ''type'' can only be ''parameter'' or ''log-posterior''.')
+        error('The ''type'' of plot is unknown. ''type'' can only be ''parameters'' or ''log-posterior''.')
     end
 end
 

private/performPT.m

@@ -77,7 +77,7 @@
    res.temperatures = nan(nIter, nTemps);
 
    beta = linspace(1,1/nTemps,nTemps).^exponentT;
-   if strcmp(temperatureAdaptionScheme,'Vousden16')
+   if strcmp(temperatureAdaptionScheme,'Vousden16') && nTemps > 1
       beta(end) = 0;
    end
    T = ones(1,nTemps);

tests/ConversionReactionTest.m

@@ -115,6 +115,7 @@ function testSamplingPT(testCase)
 
             optionsSampling = PestoSamplingOptions();
             optionsSampling.nIterations = 100;
+            optionsSampling.mode = 'silent';
 
             optionsSampling.samplingAlgorithm   = 'PT';
             optionsSampling.PT.nTemps           = 3;