expressionMCS(parameters,expression,n,varargin)

function varargout = expressionMCS(parameters,expression,n,varargin)

% Monte Carlo Simulation for generic Matlab commands.
% 
% This practical function provides a tool for very fast development
% of Monte Carlo Simulation frameworks.
% 
% syntax: [  E(expression1) , E(expression2) , ...  ] = expressionMCS( initialConditions , command , N , expression1, expression2, ... ) 
% 
% example 1
% the following line determined the expected value of a random number uniformly distributed
% over the [0,1] interval 
%     
% Ey = expressionMCS({},'y = rand',1000,'y')
% 
% The expressionMCS function executes the command 'y = rand' for 1000 times and determines 
% the expected value of the expression 'y'
% 
% 
% example 2
% the following code determines the first and second order moments of a random number with N(mu,sigma)
% distribution
%  
% mu = 1; sigma = 2; 
% [y,y2] = expressionMCS({'mu' mu ; 'sigma' sigma},'y = normrnd(mu,sigma)',10000,'y','y^2'); 
% s = (y2 - y^2)^.5% 
% 
% The main difference from example 1, is that some variable (mu and sigma) are passed as input.
% If the user wishes to pass inputs to the expressionMCS function, such as
% 
% var1name = var1value;
% var2name = var2value;
% var3name = var3value;
% ...
%     
% the inputs must be wrapped within a cell matrix with following format
% 
% initialConditions = {
% 'var1name'  var1value
% 'var2name'  var2value
% 'var3name'  var3value
% ...
% }
% 
% 
% example 3
% the following code deals with a random variable 'y' with exponential distribution
% the first sampled expression 'y>= 0 && y < 1' is equal to one if 'y' belongs to [0,1) and zero otherwise. The expectation 
% is the probability for 'y' to belong to [0,1), and it is tored within the output variable 'p1'
% Similarily, 'p2' is the probability for 'y' to belong to [1,2), and so on
% 
% [p1,p2,p3,p4,Expectation] = expressionMCS({},'y = exprnd(1)',50000,'y>= 0 && y < 1','y>= 1 && y < 2','y>= 2 && y < 3','y>= 3','y');
% bar([.5 1.5 2.5 3.5],[p1,p2,p3,p4]);
% hold on;
% plot([0:.1:4],exppdf([0:.1:4],1),'r');
% hold off;
%
%
% 
% example 4
% Consider the network
%              o
%             /|\
%            / | \
%        X1 /  |  \ X4
%          /   |   \
%         /    |    \
%      A o   X3|     o B
%         \    |    /
%          \   |   /
%        X2 \  |  / X5
%            \ | /
%             \|/
%              o
%              
% where distances X1,X2,X3,X4,X5 are independent random variable exponentially distributed 
% with mean u = [u1,u2,u3,u4,u5].
% We want to calculate the probability p for the shortest path from node A
% to node B, to be greater than 1/2
% 
% 
% u = [0.25 0.4 0.1 0.3 0.2];
% parameters = {'u' u};
% expression = [];
% expression = [expression 'X = exprnd(u); S = min([X(1) + X(4) ; X(2) + X(5) ; X(1) + X(3) + X(5) ;  X(2) + X(3) + X(4) ]);'];
% expression = [expression 'I = S>=1/2;'];
% p = expressionMCS(parameters,expression,10000,'I')

progress = 0;
prog_bar = [0:10:100];
t = cputime;

for i=1:size(parameters,1)
    variable =  parameters{i,2};
    eval([parameters{i,1} ' =  variable;']);
end


eval([expression ';']);
for i =1:nargin - 3
    eval(['varargout{i} = ' varargin{i} ';']);
end



for i = 2:n
    
    
    this_prog = prog_bar(find(abs(prog_bar - 100*i/n) == min(abs(prog_bar - 100*i/n))));
    if this_prog ~= progress
        progress = this_prog;
        time = (cputime - t)*(100/progress - 1)/60;
        fprintf('\n%i percent of the process done;  time left approximately %1.1f minutes',progress,time);
    end
    
    eval([expression ';']);
    for i =1:nargin - 3
%         eval(['varargout{i} = ( varargout{i}.*(n-1) + (' varargin{i} '))./n;']);
        eval(['varargout{i} =  varargout{i} + (' varargin{i} ');']);
    end
end

    for i =1:nargin - 3
        varargout{i}  =  varargout{i}/n;
    end