sampleCbModel

PURPOSE ^

sampleCbModel Sample the solution-space of a constraint-based model

SYNOPSIS ^

function [modelSampling,samples] = sampleCbModel(model,sampleFile,samplerName,options)

DESCRIPTION ^

sampleCbModel Sample the solution-space of a constraint-based model

 [modelSampling,samples] = sampleCbModel(model,sampleFile,samplerName,options)

INPUTS
 model         COBRA model structure
 sampleFile    File names for sampling output files

OPTIONAL INPUTS
 samplerName   Name of the sampler to be used to sample the solution
 (default 'ACHR')
 options       Options for sampling and pre/postprocessing (default values
               in parenthesis)
   nWarmupPoints           Number of warmup points (5000)
   nFiles                  Number of output files (10)
   nPointsPerFile          Number of points per file (1000)
   nStepsPerPoint          Number of sampler steps per point saved (200)
   nPointsReturned         Number of points loaded for analysis (2000)
   nFilesSkipped           Number of output files skipped when loading
                           points to avoid potentially biased initial samples (2)
   removeLoopsFlag         Attempt to remove loops from the model (false)
   removeLoopSamplesFlag   Remove sampling data for reactions involved in
                           loops (true)

OUTPUTS
 modelSampling Cleaned up model used in sampling
 samples       Uniform random samples of the solution space

 Examples of use:

 1)    Sample a model called 'superModel' using default settings and save the
       results in files with the common beginning 'superModelSamples'

           [modelSampling,samples] = sampleCbModel(superModel,'superModelSamples');

 2)    Sample a model called 'hyperModel' using default settings except with a total of 50 sample files
       saved and with 5000 sample points returned. 
       Save the results in files with the common beginning 'hyperModelSamples'

           options.nFiles = 50;
           options.nPointsReturned = 5000;
           [modelSampling,samples] = sampleCbModel(hyperModel,'hyperModelSamples');

 Markus Herrgard 8/14/06

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [modelSampling,samples] = sampleCbModel(model,sampleFile,samplerName,options)
0002 %sampleCbModel Sample the solution-space of a constraint-based model
0003 %
0004 % [modelSampling,samples] = sampleCbModel(model,sampleFile,samplerName,options)
0005 %
0006 %INPUTS
0007 % model         COBRA model structure
0008 % sampleFile    File names for sampling output files
0009 %
0010 %OPTIONAL INPUTS
0011 % samplerName   Name of the sampler to be used to sample the solution
0012 % (default 'ACHR')
0013 % options       Options for sampling and pre/postprocessing (default values
0014 %               in parenthesis)
0015 %   nWarmupPoints           Number of warmup points (5000)
0016 %   nFiles                  Number of output files (10)
0017 %   nPointsPerFile          Number of points per file (1000)
0018 %   nStepsPerPoint          Number of sampler steps per point saved (200)
0019 %   nPointsReturned         Number of points loaded for analysis (2000)
0020 %   nFilesSkipped           Number of output files skipped when loading
0021 %                           points to avoid potentially biased initial samples (2)
0022 %   removeLoopsFlag         Attempt to remove loops from the model (false)
0023 %   removeLoopSamplesFlag   Remove sampling data for reactions involved in
0024 %                           loops (true)
0025 %
0026 %OUTPUTS
0027 % modelSampling Cleaned up model used in sampling
0028 % samples       Uniform random samples of the solution space
0029 %
0030 % Examples of use:
0031 %
0032 % 1)    Sample a model called 'superModel' using default settings and save the
0033 %       results in files with the common beginning 'superModelSamples'
0034 %
0035 %           [modelSampling,samples] = sampleCbModel(superModel,'superModelSamples');
0036 %
0037 % 2)    Sample a model called 'hyperModel' using default settings except with a total of 50 sample files
0038 %       saved and with 5000 sample points returned.
0039 %       Save the results in files with the common beginning 'hyperModelSamples'
0040 %
0041 %           options.nFiles = 50;
0042 %           options.nPointsReturned = 5000;
0043 %           [modelSampling,samples] = sampleCbModel(hyperModel,'hyperModelSamples');
0044 %
0045 % Markus Herrgard 8/14/06
0046 
0047 % Default options
0048 nWarmupPoints = 5000;
0049 nFiles = 10;
0050 nPointsPerFile = 1000;
0051 nStepsPerPoint = 200;
0052 nPointsReturned = 2000;
0053 nFilesSkipped = 2;
0054 removeLoopsFlag = false;
0055 removeLoopSamplesFlag = true;
0056 
0057 if (nargin < 3)
0058     samplerName = 'ACHR';
0059 end
0060 
0061 % Handle options
0062 if (nargin > 3)
0063     if (isfield(options,'nWarmupPoints'))
0064         nWarmupPoints = options.nWarmupPoints;
0065     end
0066     if (isfield(options,'nFiles'))
0067         nFiles = options.nFiles;
0068     end
0069     if (isfield(options,'nPointsPerFile'))
0070         nPointsPerFile = options.nPointsPerFile;
0071     end
0072     if (isfield(options,'nStepsPerPoint'))
0073         nStepsPerPoint = options.nStepsPerPoint;
0074     end
0075     if (isfield(options,'nPointsReturned'))
0076         nPointsReturned = options.nPointsReturned;
0077     end
0078     if (isfield(options,'nFilesSkipped'))
0079         nFilesSkipped = options.nFilesSkipped;
0080     end
0081     if (isfield(options,'removeLoopsFlag'))
0082         removeLoopsFlag = options.removeLoopsFlag;
0083     end
0084     if (isfield(options,'removeLoopSamplesFlag'))
0085         removeLoopSamplesFlag = options.removeLoopSamplesFlag;
0086     end
0087 end
0088 
0089 
0090 fprintf('Prepare model for sampling\n');
0091 % Prepare model for sampling by reducing bounds
0092 [nMet,nRxn] = size(model.S);
0093 fprintf('Original model: %d rxns %d metabolites\n',nRxn,nMet);
0094 
0095 % Reduce model
0096 fprintf('Reduce model\n');
0097 modelRed = reduceModel(model, 1e-6, false,false);
0098 [nMet,nRxn] = size(modelRed.S);
0099 fprintf('Reduced model: %d rxns %d metabolites\n',nRxn,nMet);
0100 save modelRedTmp modelRed;
0101 
0102 modelSampling = modelRed;
0103 
0104 switch samplerName
0105     case 'ACHR'
0106         % Use Artificial Centering Hit-and-run
0107         
0108         fprintf('Create warmup points\n');
0109         % Create warmup points for sampler
0110         warmupPts= createHRWarmup(modelSampling,nWarmupPoints);
0111 
0112         save sampleCbModelTmp modelSampling warmupPts
0113 
0114         fprintf('Run sampler for a total of %d steps\n',nFiles*nPointsPerFile*nStepsPerPoint);
0115         % Sample model
0116         ACHRSampler(modelSampling,warmupPts,sampleFile,nFiles,nPointsPerFile,nStepsPerPoint);
0117 
0118     otherwise
0119         error(['Unknown sampler: ' samplerName]);
0120 end
0121 
0122 fprintf('Load samples\n');
0123 % Load samples
0124 nPointsPerFileLoaded = ceil(nPointsReturned/(nFiles-nFilesSkipped));
0125 if (nPointsPerFileLoaded > nPointsPerFile)
0126    error('Attempted to return more points than were saved'); 
0127 end
0128 samples = loadSamples(sampleFile,nFiles,nPointsPerFileLoaded,nFilesSkipped);
0129 samples = samples(:,round(linspace(1,size(samples,2),nPointsReturned)));
0130 % Fix reaction directions
0131 [modelSampling,samples] = convRevSamples(modelSampling,samples);

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003