diff --git a/analysis/randomSampling.m b/analysis/randomSampling.m index 7042370f..12f89184 100755 --- a/analysis/randomSampling.m +++ b/analysis/randomSampling.m @@ -1,4 +1,4 @@ -function [solutions, goodRxns]=randomSampling(model,varargin) +function [solutions, goodRxns, info]=randomSampling(model,varargin) % randomSampling Sample the flux solution space (entry point for all samplers). % % Dispatches to one of three sampling methods via the 'method' argument: @@ -75,6 +75,10 @@ % [randomObjective] vector of indexes of those reactions that are not % involved in loops or always carry zero flux and can be used as random % objective functions. Empty ([]) for the 'achr' and 'chrr' methods. +% info : struct +% [chrr] the diagnostics from sampleCHRR, including mveConverged (whether +% the maximum-volume ellipsoid rounding reached tolerance). Empty ([]) +% for the 'achr' and 'randomObjective' methods. % % Notes % ----- @@ -111,6 +115,7 @@ method=p.method; thinning=p.thinning; nBurnin=p.nBurnin; +info=[]; %only populated by the 'chrr' method if isempty(nSamples) nSamples=1000; end @@ -129,7 +134,7 @@ goodRxns = []; return; case 'chrr' - solutions = sampleCHRR(model, nSamples, thinning, nBurnin, seed); + [solutions, info] = sampleCHRR(model, nSamples, thinning, nBurnin, seed); goodRxns = []; return; case {'randomobjective', 'objective'} diff --git a/analysis/sampleCHRR.m b/analysis/sampleCHRR.m index 8c3d4853..dfa8f67b 100644 --- a/analysis/sampleCHRR.m +++ b/analysis/sampleCHRR.m @@ -118,6 +118,11 @@ [center, E, converged] = sampleMaxVolEllipse(A_full, b_full, x0); info.nDimensions = d; info.mveConverged = converged; +if ~converged + warning('RAVEN:warning', '%s', ['The maximum-volume ellipsoid rounding ' ... + 'did not converge; the samples may be poorly mixed. Inspect the ' ... + 'mveConverged field of the second output.']); +end % Rounded polytope {y : A_r*y <= b_r}, which contains the unit ball. A_r = A_full * E; diff --git a/testing/function_tests/tSampling.m b/testing/function_tests/tSampling.m index 2e984f20..c365addb 100644 --- a/testing/function_tests/tSampling.m +++ b/testing/function_tests/tSampling.m @@ -109,6 +109,14 @@ function randomSamplingCHRRDispatch(testCase) testCase.verifyLessThan(max(abs(resid(:))), 1e-6); end + function randomSamplingCHRRExposesInfo(testCase) + % The CHRR convergence diagnostics must be reachable through the + % documented randomSampling entry point, not only sampleCHRR. + evalc(['[sols, gr, sInfo] = randomSampling(testCase.model, 10, ' ... + '''method'', ''chrr'', ''thinning'', 5, ''nBurnin'', 20, ''seed'', 1);']); + testCase.verifyTrue(isfield(sInfo, 'mveConverged')); + end + function randomSamplingCHRRDeterministic(testCase) evalc(['s1 = randomSampling(testCase.model, 8, ''method'', ''chrr'', ' ... '''thinning'', 5, ''nBurnin'', 20, ''seed'', 42);']);