How can I run the DREAM_Suite on mutliple data files without modifying the main code

I have this code that runs DREAM_Suite (attached) for data fitting and optimization using a Markov Chain algorithm. I can run this code on a single data file at a time, but I'd like to be able to run this optimization on multiple data files of the same kind, one after the other. Is it possible to do this through a simple modification in this "driver" script and withought having to modify/recompile the main code and associated functions in the DREAM package. (One has to download/install the DREAM package to run the attched)...

I'd also like to silence the postprocess output which is horendous as designed (muliple output figures and files, etc). I only need the file with the fit parameters (called therein 'DREAM_ZS_output.txt') as output from the individual data file fitting.

I am running Matlab 2025a version.

[Stephen23]

Did you try calling DREAM_Suite in a FOR-loop?

https://www.mathworks.com/help/matlab/getting-started-with-matlab.html

[dpb]

You uploaded either two copies of the driver m file or the same file twice with the same name? Not easy to do anything with that...please straighten out the confusion in the proper manner; either delete the extra or rename one or both to be distinct.

[dpb]

type driver_example

% clear all
% close all
% fclose('all')
% 
% addpath 'C:\Program Files\MATLAB\R2025a\DREAM-Suite-main\MATLAB'
% Setup the Dream_Suite optimization/fittin, on single data file 

tab1=readtable('datafile1.csv');

lipC =tab1.Cl; % lipid concetration in mM
pepC =tab1.Cpt; % lipid concetration in mM
normI =tab1.normInt;


%% DEFINE YOUR MEASUREMENT DATA, [X,Y, ERR]


CL0=lipC*1e3; 
CPT0=pepC*1e3; 
normI0=abs(normI); 

data0= [normI0 CL0 CPT0]; 

%% Sort data in the ascending order by 1st colunm (here x = normI)
%% X become Y for fitting purposes
data1=sortrows(data0,-1, 'ascend'); 
x=data1(:,1); 
y=data1(:,2);
dy=0.1;

%% The setting "DREAMPar.lik = 12" below says gaussian residuals
%% Problem settings defined by user
DREAMPar.d = 3;                         % Dimension of the problem
DREAMPar.lik = 12;                      % Gaussian likelihood with measurement error

%% Provide information parameter space and initial sampling
Par_info.initial = 'latin';             % Latin hypercube sampling
Par_info.boundhandling = 'reflect';     % Explicit boundary handling

% %Bounds for n, kD and Iinf
Par_info.min = [1     1e-3       round(max(normI),1)-1];                    % If 'latin', min values 
Par_info.max = [100   1e+3       round(max(normI),1)+1];                    % If 'latin', max values


%% Define name of function (.m file) for posterior exploration
Func_name = 'model_example';

%% Load the data. X is passed to model as plugin.X

plugin.X = x;
Meas_info.Y = y;
Meas_info.Sigma = dy;

%% Optional settings
options.modout = 'yes';                % Return model (function) simulations of samples (yes/no)?

%% Define method to use {'dream','dream_zs','dream_d','dream_dzs','mtdream_zs'}
method = 'dream_zs';

switch method
    case {'dream','dream_d'}
        DREAMPar.N = 10;                            % # Markov chains
        DREAMPar.T = 20000;                          % # generations
    case {'dream_zs','dream_dzs','mtdream_zs'}
        DREAMPar.N = 3;                             % # Markov chains
        DREAMPar.T = 10000;                          % # generations
        %DREAMPar.T = 150;  % fast check that program runs to completion
end

if strcmp(method,'dream_d') || strcmp(method,'dream_dzs')
    Par_info.steps = 1000*ones(1,DREAMPar.d);       % # discrete steps
end

%% Run DREAM-Suite - Paul
 [chain,output,FX,Z,logL] = DREAM_Suite(method,Func_name,DREAMPar,Par_info,Meas_info,options,[],plugin);

[dpb]

The second one

type driver_example

clear all
close all
fclose('all')

addpath 'C:\Program Files\MATLAB\R2025a\DREAM-Suite-main\MATLAB'

% Setup the Dream_Suite optimization/fittin, on single data file 
tab1=readtable('datafile2.csv');

lipC =tab1.Cl;
pepC =tab1.Cpt; 
normI =tab1.normInt;


%% DEFINE YOUR MEASUREMENT DATA, [X,Y, ERR]

CL0=lipC*1e3; 
CPT0=pepC*1e3; 
normI0=abs(normI); 

data0= [normI0 CL0 CPT0]; 

%% Sort data in the ascending order by 1st colunm (here x = normI)
%% X become Y for fitting purposes on this complicated model
data1=sortrows(data0,-1, 'ascend'); 
x=data1(:,1); 
y=data1(:,2);
dy=0.0001+0.1*y;


%% The setting "DREAMPar.lik = 12" below says gaussian residuals
%% Problem settings defined by user
DREAMPar.d = 3;                         % Dimension of the problem
DREAMPar.lik = 12;                      % Gaussian likelihood with measurement error

%% Provide information parameter space and initial sampling
Par_info.initial = 'latin';             % Latin hypercube sampling
Par_info.boundhandling = 'reflect';     % Explicit boundary handling

% % Bounds for n, kD and Iinf
Par_info.min = [1     1e-3       round(max(normI),1)-1];                    % If 'latin', min values 
Par_info.max = [100   1e+3       round(max(normI),1)+1];                    % If 'latin', max values


%% Define name of function (.m file) for posterior exploration
Func_name = 'model_example';

%% Load the data. X is passed to model as plugin.X

plugin.X = x;
Meas_info.Y = y;
Meas_info.Sigma = dy;

%% Optional settings
options.modout = 'yes';                % Return model (function) simulations of samples (yes/no)?

%% Define method to use {'dream','dream_zs','dream_d','dream_dzs','mtdream_zs'}
method = 'dream_zs';

switch method
    case {'dream','dream_d'}
        DREAMPar.N = 10;                            % # Markov chains
        DREAMPar.T = 20000;                          % # generations
    case {'dream_zs','dream_dzs','mtdream_zs'}
        DREAMPar.N = 3;                             % # Markov chains
        DREAMPar.T = 10000;                          % # generations
        %DREAMPar.T = 150;  % fast check that program runs to completion
end

if strcmp(method,'dream_d') || strcmp(method,'dream_dzs')
    Par_info.steps = 1000*ones(1,DREAMPar.d);       % # discrete steps
end

%% Run DREAM-Suite - Paul
 [chain,output,FX,Z,logL] = DREAM_Suite(method,Func_name,DREAMPar,Par_info,Meas_info,options,[],plugin);

[dpb]

type model_example

function [Y] = model_MCMC_FL(pa, plugin)

CPT=10; % constant

n = pa(1); kD = pa(2); Iinf=pa(3);
X = plugin.X;

Y = CPT*(X-1).*(n./(Iinf-1)) + kD*(X-1).*(n./(Iinf-X));

There's nothing of interest here other than the functional being used except the function is missing the END statement that well-written code should use even though can get away without it in one-function m-files.

Answer 1

[dpb]

From a quick perusal of the two driver scripts, it looks as though if you were to convert those to a function it would be straightforward as @Stephen23 suggests to use a loop...

Since you stated "...I only need the file with the fit parameters (called therein 'DREAM_ZS_output.txt') as output", the following example function adds a sequence number to the base file name and returns it...the top level routine can save that as a cellstr array for later access.

If you need to be able to modify the other input parameters, you'll have to add those variables into the argument list of passed variables and build an input file from which to read the needed data for each case.

function outfile=dreamDriver(datafile,Func_name)
    % Setup the Dream_Suite optimization/fitting, on single data file supplied
    % as function argument.  Optionally pass in the evaluation function m-file
    % name, default is 'model_example'
    % Returns modified sequnteially numbered output file name as char string
    
    % set the default function name if not provided...
    if nargin==1, Func_name='model_example'; end    % set default name as desired
    % read the input data file
    tab1=readtable(datafile);
    
    lipC =tab1.Cl;
    pepC =tab1.Cpt; 
    normI =tab1.normInt;
    %% DEFINE YOUR MEASUREMENT DATA, [X,Y, ERR]
    CL0=lipC*1e3; 
    CPT0=pepC*1e3; 
    normI0=abs(normI); 
    data0= [normI0 CL0 CPT0]; 
    
    %% Sort data in the ascending order by 1st colunm (here x = normI)
    %% X become Y for fitting purposes on this complicated model
    data1=sortrows(data0,-1, 'ascend'); 
    x=data1(:,1); 
    y=data1(:,2);
    dy=0.0001+0.1*y;
    %% The setting "DREAMPar.lik = 12" below says gaussian residuals
    %% Problem settings defined by user
    DREAMPar.d = 3;                         % Dimension of the problem
    DREAMPar.lik = 12;                      % Gaussian likelihood with measurement error
    %% Provide information parameter space and initial sampling
    Par_info.initial = 'latin';             % Latin hypercube sampling
    Par_info.boundhandling = 'reflect';     % Explicit boundary handling
    
    % % Bounds for n, kD and Iinf
    Par_info.min = [1     1e-3       round(max(normI),1)-1];                    % If 'latin', min values 
    Par_info.max = [100   1e+3       round(max(normI),1)+1];                    % If 'latin', max values
    %% Define name of function (.m file) for posterior exploration
    % Use the dummy argument value passed...or default if not passed - 
    %Func_name = 'model_example';
    
    %% Load the data. X is passed to model as plugin.X
    plugin.X = x;
    Meas_info.Y = y;
    Meas_info.Sigma = dy;
    
    %% Optional settings
    options.modout = 'yes';                % Return model (function) simulations of samples (yes/no)?
    
    %% Define method to use {'dream','dream_zs','dream_d','dream_dzs','mtdream_zs'}
    method = 'dream_zs';
    switch method
        case {'dream','dream_d'}
            DREAMPar.N = 10;                            % # Markov chains
            DREAMPar.T = 20000;                          % # generations
        case {'dream_zs','dream_dzs','mtdream_zs'}
            DREAMPar.N = 3;                             % # Markov chains
            DREAMPar.T = 10000;                          % # generations
            %DREAMPar.T = 150;  % fast check that program runs to completion
    end
    if strcmp(method,'dream_d') || strcmp(method,'dream_dzs')
        Par_info.steps = 1000*ones(1,DREAMPar.d);       % # discrete steps
    end
    %% Run DREAM-Suite - Paul
    [chain,output,FX,Z,logL] = DREAM_Suite(method,Func_name,DREAMPar,Par_info,Meas_info,options,[],plugin);
  % move default file to a sequentially-nunmbered file of the same base output file name  
  outfile=sprintf('DREAM_ZS_output_%02d.txt',i); % create sequential output file name
  movefile('DREAM_ZS_output.txt',outfile)        % move the default file to the sequentially named one
end

Then can simply call it as

d=dir('datafile*.csv');         % get the list of desired data files with appropriate wildcard spec
for i=1:numel(d)                                    % iterate over the list
  % run the app, save the list of output file names....  
  resultFiles{i}=dreamDriver(fullfile(d(i).folder,d(i).name));
end

The above will run all the cases; it appears that the only option in this driver script (now function) to change anything about output is the one parameter

    options.modout = 'yes';                % Return model (function) simulations of samples (yes/no)?

You can try setting that to 'no' and see if that helps or not, otherwise you'll have to see if there are any other options to control how much and what ouput is generated from the documentation. Lacking those, you'll need to go into the guts of the program and find the outputs that you don't want and either comment out those sections of code or add some additional granularity or options to the options if the authors did not incorporate any other choices.

The only alternative lacking that facility I would see without code modifications from what is available to inspect here would be to do similar as to what is done here in renaming the default output file to something unique--find the standard names of the other files not wanted and DELETE() them after each pass through the loop and CLOSE all the figures. You could try setting the default top level figure visibility to 'off' and see if the app respects that and if so, while they would still have been generated, they wouldn't be visibile. Some apps like this tend to consider themselves as the only kid in town and all important and will force visibility so that may not work. Those details of how the DREAM package works are out of the scope of the MATLAB Answers forum, you'll be better served to ask those details in its support group.

For the above workaround/solution, you could alternatively use <fileparts> to get the base input data file name and modify it to reflect output and then have a set of resultant files that are related to the input files. Of course, if you don't name a separate series of input files with something that identifies them uniquely but use the same naming convention every time, you'll still not be able to identify them or will overwrite the prior ones of the same name. Note that the array of filenames variable is independent of the names of the files, themselves; the file names will still be the base output name the app uses with just a sequence number inserted -- the sequence number will still be from 1:NumberInputFiles so a new run will replace that many existing from a prior run unless they're named something unique to avoid doing that as mentioned above.

[dpb]

ADDENDUM

From the SWITCH statement on method and the example output file name, one guesses that the code generates the output file name as

outputFilename=sprintf('%s_output.txt',upper(method));

so a generic solution would need to look for and modify the specific name created if it is intended to run alternate methods. That would also entail passing the method to the function to be able to switch it as well, of course.

Note the capitalization mismatch -- if running under Windows, it won't matter; file names are case insensitive although case-preserving. If under the other OS variety, then they are also case sensitive so it would matter.