function [stoichiometricMatrix,kineticParameters,initialConditions] = importSBML
% importSBML
% This function imports an SBML model and determines the data required
% to be used for simulation using the methods presented in the article:
% "Deterministic Modelling and Stochastic Simulation of Pathways
% using MATLAB".
%
% REQUIREMENTS:
% =============
% 1) The SBML Toolbox needs to be installed prior to using this function.
%    It can be found on www.sbml.org. The folder in which the function
%    'TranslateSBML' is located needs to be added to the MATLAB path.
% 2) It is assumed that all kinetic rate expressions are in the form of
%    mass action kinetics. NO(!) error checking with respect to the type of
%    the kinetic rate laws is done. 
% 3) Only SBML Level 2 Version 1 models can be imported.
% 4) If compartments are present, there size is required to be 1
% 5) The SBML model is not allowed to contain:
%       - rules
%       - events
%       - functions
%       - reversible reactions
%
% USAGE:
% ======
% [stoichiometricMatrix,kineticParameters,initialConditions] = importSBML
%
% The function opens a dialog that is asking the user to choose an SBML
% model file. The output arguments are the stoichiometric matrix, the
% kinetic parameters, and the initial conditions for the species.
%
% INFORMATION:
% ============
% This function is based on code included in the Systems Biology Toolbox
% for MATLAB, found at www.sbtoolbox.org. It should be seen as an example
% on how one can quickly convert SBML models in a more handy representation
% in terms of a MATLAB structure based on ODEs. The stoichiometric matrix
% and the kinetic parameters are subsequently determined from this more
% handy structure.
%
% Copyright (C) 2005 Henning Schmidt, henning@fcc.chalmers.se

% check if SBML Toolbox is present
if exist('TranslateSBML') ~= 3,
    error('Please make sure the SBML Toolbox is installed and included in the MATLAB path.');
end

% load an SBML model using the SBML Toolbox
SBMLmodel = TranslateSBML;

% Check things
if SBMLmodel.SBML_level ~= 2, error('Wrong SBML Level. Only Level 2 Version 1 models can be imported.'); end
if SBMLmodel.SBML_version ~= 1, error('Wrong SBML Version. Only Level 2 Version 1 models can be imported.'); end
if ~isempty(SBMLmodel.functionDefinition), error('SBML functions are not supported.'); end
if ~isempty(SBMLmodel.rule), error('SBML rules are not supported.'); end
if ~isempty(SBMLmodel.event), error('SBML events are not supported.'); end
if sum([SBMLmodel.reaction.reversible]) > 0, error('Reversible reactions are not allowed in the model.'); end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CREATE AN EMPTY SBSTRUCTURE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create empty SBstructure
% functions substructure
functionsStruct = struct('name',{},'arguments',{},'formula',{},'notes',{});
% states substructure
statesStruct = struct('name',{},'initialCondition',{},'ODE',{},'type',{},'compartment',{},'unittype',{},'notes',{});
% parameters substructure
parametersStruct = struct('name',{},'value',{},'type',{},'compartment',{},'unittype',{},'notes',{});
% variables substructure
variablesStruct = struct('name',{},'formula',{},'type',{},'compartment',{},'unittype',{},'notes',{});
% reactions substructure
reactionsStruct = struct('name',{},'formula',{},'notes',{},'reversible',{});
% event assignment substructure
eventassignmentStruct = struct('variable',{},'formula',{});
% event substructure
eventStruct = struct('name',{},'trigger',{},'assignment',eventassignmentStruct,'notes',{});
% Create SBstructure
SBstructure = struct('name','','notes','','functions',functionsStruct,'states',statesStruct,'parameters',parametersStruct,'variables',variablesStruct,'reactions',reactionsStruct,'events',eventStruct,'functionsMATLAB','');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NAME AND NOTES (INCLUDE IN STRUCTURE)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~strcmp(SBMLmodel.name,''),
    SBstructure.name = SBMLmodel.name;
elseif ~strcmp(SBMLmodel.id,''),
    SBstructure.name = SBMLmodel.id;
else
    SBstructure.name = 'Unknown model';
end
% Indicate in notes that the model has been created by import from SBML.
SBstructure.notes = SBMLmodel.notes;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SOME COUNTER DEFINITIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
numberStates = 1;
numberParameters = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SPECIES (INCLUDE IN STRUCTURE)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
speciesODElist = [];
speciesODElistIndex = 1;
% Cycle through all the species and determine how they are to be included
% in the model structure
for k1 = 1:length(SBMLmodel.species),
    speciesName = SBMLmodel.species(k1).id;
    speciesCompartmentName = SBMLmodel.species(k1).compartment;
    % determine the initial value for species
    initialCondition = [];
    if SBMLmodel.species(k1).isSetInitialAmount,
        initialCondition = SBMLmodel.species(k1).initialAmount;
        speciesUnits = 'amount';
    elseif SBMLmodel.species(k1).isSetInitialConcentration,
        % convert concentration to amount
        initialCondition = SBMLmodel.species(k1).initialConcentration;
        speciesUnits = 'concentration';
    else
        % if nothing is chosen, assume concentration units
        speciesUnits = 'concentration';
    end
    % Process all the species that DO NOT have a boundary condition
    % In case the species is non-constant it is included as for Level 1
    % In case the species is constant it is included as parameter
    if SBMLmodel.species(k1).constant == 0,
        if SBMLmodel.species(k1).boundaryCondition == 0,  
            % add species to the ODE list (for use below in ODE
            % construction)
            speciesODElist(speciesODElistIndex).name = speciesName;
            speciesODElist(speciesODElistIndex).stateIndex = numberStates;
            speciesODElist(speciesODElistIndex).units = speciesUnits;
            speciesODElistIndex = speciesODElistIndex + 1;
            % Include species as state
            SBstructure.states(numberStates).name = speciesName;
            % Check if an initial value is set for a species
            if length(initialCondition)==0,
                error(sprintf('No initial value defined for species ''%s''\n',speciesName));
            end
            SBstructure.states(numberStates).initialCondition = initialCondition;
            % Note that this state is a 'species'
            % Notes are not used by the toolbox for other things as
            % documentation
            SBstructure.states(numberStates).notes = '';
            SBstructure.states(numberStates).type = 'isSpecie';
            SBstructure.states(numberStates).compartment = SBMLmodel.species(k1).compartment;
            SBstructure.states(numberStates).unittype = speciesUnits;
            % Initialize the ODE field for this state
            SBstructure.states(numberStates).ODE = '';
            % Increment
            numberStates = numberStates+1;
            % The case where numberRulesSpecies > 1 does not have to be treated
            % here since already an error message is set above. This error
            % message will lead to an empty model structure at the end. Even
            % if the error is set the rest of the model will be processed.
            % This has the advantage that all errors can be
            % detected at once and makes the control structure of this script
            % simpler.
        else
            % Include boundary species as a parameter
            SBstructure.parameters(numberParameters).name = speciesName;
            % Check if an initial amount is set for a species
            if length(initialCondition)==0,
                error(sprintf('No initial value defined for species ''%s''\n',speciesName));
            end
            % Include initial value
            SBstructure.parameters(numberParameters).value = initialCondition;
            % Set note that type of the parameter is 'boundaey species'
            SBstructure.parameters(numberParameters).notes = '';
            % initialize type, compartment, and unittype fields
            SBstructure.parameters(numberParameters).type = 'isParameter';
            SBstructure.parameters(numberParameters).compartment = '';
            SBstructure.parameters(numberParameters).unittype = '';
            % Increment
            numberParameters = numberParameters+1;
        end
    else
        % Species defined as being constant!
        % Thus include it as parameter
        SBstructure.parameters(numberParameters).name = speciesName;
        % Check if an initial amount is set for species
        if length(initialCondition)==0,
            error(sprintf('No initial value defined for species ''%s''\n',speciesName));
        end
        % Include Initial value in SBstructure
        SBstructure.parameters(numberParameters).value = initialCondition;
        % Set note that type of the parameter is 'constant species'
        SBstructure.parameters(numberParameters).notes = '';
        % initialize type, compartment, and unittype fields
        SBstructure.parameters(numberParameters).type = 'isParameter';
        SBstructure.parameters(numberParameters).compartment = '';
        SBstructure.parameters(numberParameters).unittype = '';
        % Increment
        numberParameters = numberParameters+1;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PARAMETERS (INCLUDE IN STRUCTURE)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start with global parameters
for k = 1:length(SBMLmodel.parameter),
    parameterName = SBMLmodel.parameter(k).id;
    parameterValue = SBMLmodel.parameter(k).value;
    % Include parameter as parameter
    SBstructure.parameters(numberParameters).name = parameterName;
    % Check if a value is set for the parameter
    if length(parameterValue)==0,
        error(sprintf('No value defined defined for parameter ''%s''\n',parameterName));
    end
    % Include value in SBstructure
    SBstructure.parameters(numberParameters).value = parameterValue;
    % Set note that type of the parameter is 'global'
    SBstructure.parameters(numberParameters).notes = '';
    % initialize type, compartment, and unittype fields
    SBstructure.parameters(numberParameters).type = 'isParameter';
    SBstructure.parameters(numberParameters).compartment = '';
    SBstructure.parameters(numberParameters).unittype = '';
    % Increment
    numberParameters = numberParameters+1;
end
% Then include the local parameters within the different reactions
% Name collisions are avoided by adding the reaction name in front of each
% local parameter.
for k1 = 1:length(SBMLmodel.reaction),
    % get the current reaction name
    reactionName = SBMLmodel.reaction(k1).id;
    if ~isempty(SBMLmodel.reaction(k1).kineticLaw),
        for k2 = 1:length(SBMLmodel.reaction(k1).kineticLaw.parameter),
            % get the current local parameter in reaction reactionName
            parameterName = SBMLmodel.reaction(k1).kineticLaw.parameter(k2).id;
            parameterValue = SBMLmodel.reaction(k1).kineticLaw.parameter(k2).value;
            % check if a value is set for the parameter
            if length(parameterValue)==0,
                error(sprintf('No value defined defined for parameter ''%s'' in reaction ''%s''\n',parameterName,reactionName));
            end
            % add the reaction name to the parameter and update the kinetic rate law
            parameterName = char([double(reactionName) double('_') double(parameterName)]);  % fastest strcat ;)
            % Include parameter in the structure
            SBstructure.parameters(numberParameters).name = parameterName;
            SBstructure.parameters(numberParameters).value = parameterValue;
            SBstructure.parameters(numberParameters).notes = '';
            % initialize type, compartment, and unittype fields
            SBstructure.parameters(numberParameters).type = 'isParameter';
            SBstructure.parameters(numberParameters).compartment = '';
            SBstructure.parameters(numberParameters).unittype = '';
            numberParameters = numberParameters+1;
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% COMPARTMENTS (INCLUDE IN STRUCTURE)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 1:length(SBMLmodel.compartment),
    compartmentValue = SBMLmodel.compartment(k).size;
    % check compartment size and return an error in case it is "NaN"
    if compartmentValue ~= 1,
        error(sprintf('All compartment sizes are required to be 1.'));
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% REACTIONS (INCLUDE IN STRUCTURE)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read in reaction information into the model structure. Replace colliding 
% parameter names with their prefixed versions (reaction name). Reactions 
% are included into the model structure in the 'reactions' field. This leads 
% to shorter ODE definitions for the specie states, which then contain the 
% reaction names with the correct stoichiometries. The 'reactions' field
% contains three fields: 'name', 'formula', 'notes'
%
% We differentiate the case where:
%   - single compartment model with constant compartment size = 1
%     (nothing is done additionally)
%   - multi compartment model or single compartment model with compartment 
%     size different from one. In this case all species names in the
%     reactions are exchanged agains speciesname/speciesCompartmentSize
%
% Cycle through all reactions
for k1 = 1:length(SBMLmodel.reaction),
    reactionName = SBMLmodel.reaction(k1).id;
    kineticLaw = SBMLmodel.reaction(k1).kineticLaw;
    % check if a kineticLaw is given
    if length(kineticLaw)==0,
        error(sprintf('No kinetic law is given for reaction ''%s''\n',reactionName));
    end
    % Process the kineticLaw formula by replacing the local parameter names
    % by reactionname_parameternames and replace the mathml expressions
    if ~isempty(kineticLaw),
        if length(kineticLaw.formula) >= length(kineticLaw.math),
            formula = kineticLaw.formula;
        else
            formula = kineticLaw.math;
        end
        if length(kineticLaw.parameter) > 0,
            % SBML 2 can have formula in formula or math. Check which string is
            % longer and use this one (quick and dirty fix - since it has been
            % observer that TranslateSBML can return strange stuff in the math
            % field sometimes)
            for k2 = 1:length(kineticLaw.parameter),
                parameterName = kineticLaw.parameter(k2).id;
                newParameterName = char([double(reactionName) double('_') double(parameterName)]); % fastest strcat
                formula = exchangeStringInString(formula,parameterName,newParameterName);
            end
        end
    else
        formula = '';
    end
    % Include reaction information in the structure
    SBstructure.reactions(k1).name = reactionName;
    SBstructure.reactions(k1).formula = formula;
    SBstructure.reactions(k1).notes = '';
    SBstructure.reactions(k1).reversible = SBMLmodel.reaction(k1).reversible;    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PROCESS REACTION ODES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Add the reaction variable names with correct stoichiometries to the 'ODE' field
% for the correct species. The 'correct' species are only non-boundary and 
% non-constant species for which no rules exist. These are listed in
% the 'speciesODElist'
%
for k1 = 1:length(speciesODElist),
    % get name of species and the index of the corresponding state
    speciesName = speciesODElist(k1).name;
    stateIndex = speciesODElist(k1).stateIndex;
    % cycle through all the reactions and look in which reactions the
    % species is changed
    for k2 = 1:length(SBstructure.reactions),
        % get reaction name for inclusion in ODEstring
        reactionName = SBstructure.reactions(k2).name;
        % cycle through the reactants of the current reaction
        % (ordering of reactions in SBstructure and SBMLmodel are equal)
        for k3 = 1:length(SBMLmodel.reaction(k2).reactant),
            reactantName = SBMLmodel.reaction(k2).reactant(k3).species;
            if length(SBMLmodel.reaction(k2).reactant(k3).stoichiometryMath) > 0,
                % stoichiometry given as formula. this is not supported =>
                % error
                error(sprintf('Stoichiometry for species ''%s'' in reaction ''%s''\ngiven as MATH expression. This is not supported.\n',speciesName,reactionName));
                reactantStoichiometry = double(1);  % just a value - its an error anyway!
            else
                reactantStoichiometry = abs(double(SBMLmodel.reaction(k2).reactant(k3).stoichiometry) / double(SBMLmodel.reaction(k2).reactant(k3).denominator));
            end
            % If species name and reactant name are equal then add reaction
            % to species ODE
            if strcmp(reactantName,speciesName),
                % construct the string to add to the ODE of the current reactant
                if reactantStoichiometry ~= 1,
                    ODEstringAdd = char([double('-') double(num2str(reactantStoichiometry)) double('*') double(reactionName)]); % fastest strcat
                else
                    ODEstringAdd = char([double('-') double(reactionName)]); % fastest strcat
                end
                % add the ODEstringAdd to the 'ODE' field of this species
                SBstructure.states(stateIndex).ODE = char([double(SBstructure.states(stateIndex).ODE) double(ODEstringAdd)]); % fastest strcat
                break;
            end
        end
        % cycle through the products of the current reaction
        % (ordering of reactions in SBstructure and SBMLmodel are equal)
        for k3 = 1:length(SBMLmodel.reaction(k2).product),
            productName = SBMLmodel.reaction(k2).product(k3).species;
            if length(SBMLmodel.reaction(k2).product(k3).stoichiometryMath) > 0,
                % stoichiometry given as formula. this is not supported =>
                % error
                error(sprintf('Stoichiometry for species ''%s'' in reaction ''%s''\ngiven as MATH expression. This is not supported.\n',speciesName,reactionName));
                productStoichiometry = double(1);  % just a value - its an error anyway!
            else
                productStoichiometry = abs(double(SBMLmodel.reaction(k2).product(k3).stoichiometry) / double(SBMLmodel.reaction(k2).product(k3).denominator));
            end
            % If species name and product name are equal then add reaction
            % to species ODE
            if strcmp(productName,speciesName),
                % construct the string to add to the ODE of the current product
                if productStoichiometry ~= 1,
                    ODEstringAdd = char([double('+') double(num2str(productStoichiometry)) double('*') double(reactionName)]); % fastest strcat
                else
                    ODEstringAdd = char([double('+') double(reactionName)]); % fastest strcat
                end
                % add the ODEstringAdd to the 'ODE' field of this species
                SBstructure.states(stateIndex).ODE = char([double(SBstructure.states(stateIndex).ODE) double(ODEstringAdd)]); % fastest strcat
                break;
            end
        end
    end
end
% Now all ODEs should be defined - check if everything is correct!
% Finally cycle through all the states and check if all ODEs have been defined
for k = 1:length(SBstructure.states),
    stateName = SBstructure.states(k).name;
    stateODE = SBstructure.states(k).ODE;
    if length(stateODE) == 0,
        % Issue an error
        error(sprintf('No ODE defined for state ''%s''',stateName));
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DETERMINE STOICHIOMETRIC MATRIX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
componentNames = {SBstructure.states.name};
reactionNames = {SBstructure.reactions.name};
N = [];
for k = 1:length(SBstructure.states),
    ODE = SBstructure.states(k).ODE;
    % delete all white spaces from the ODE
    ODE = regexprep(ODE,'\s*','');
    ODE = char([double(ODE) double('+')]);  % fastest strcat
    % first explode the ODE in additive terms
    terms = [];
    termIndex = 1;
    % check the sign of the first term (first character in string)
    if ODE(1) == '-',
        signCurrent = -1;
        lastIndex = 2;
    elseif ODE(1) == '+',
        signCurrent = +1;
        lastIndex = 2;
    else
        signCurrent = +1;
        lastIndex = 1;
    end
    % explode in terms, check if term has the right format and if the
    % second term features a reactionname. then construct the row of the
    % stoichiometric matrix
    Nrow = zeros(1,length(reactionNames));
    startk = lastIndex;
    for k = startk:length(ODE),
        % check for positive term
        if ODE(k) == '+' || ODE(k) == '-',
            element = ODE(lastIndex:k-1);
            % check the element if composed of term in the right format
            multIndex = find(element == '*');
            if length(multIndex) == 0,
                stoichiometry = signCurrent*1;
                reactionterm = element;
            elseif length(multIndex) == 1,
                absStoichiometry = str2double(element(1:multIndex-1));
                if isempty(absStoichiometry),
                    errorFlag = 1;
                    break;
                end
                stoichiometry = signCurrent*absStoichiometry;
                reactionterm = element(multIndex+1:end);
            else
                % to many multiplication signs (only one allowed)
                errorFlag = 1;
                break;
            end
            % find the index of the reaction name and add the
            % stoichiometric information to Nrow
            indexReaction = strmatch(reactionterm,reactionNames,'exact');
            if isempty(indexReaction),
                errorFlag = 1;
                break;
            end
            Nrow(indexReaction) = stoichiometry;
            % increment
            termIndex = termIndex + 1;
            lastIndex = k+1;
            if ODE(k) == '+',
                signCurrent = +1;
            else
                signCurrent = -1;
            end
        end
    end
    N = [N; Nrow];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GET PARAMETER VALUES BY A VERY SIMPLE APPROACH
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
kineticParameters = [];
regexpressions_components = {};
for k = 1:length(componentNames),
    regexpressions_components{k} = strcat('\<',strtrim(componentNames{k}),'\>');
end
regexpressions_parameters = {};
parameterValuesStrings = {};
parameterNames = {SBstructure.parameters.name};
for k = 1:length(parameterNames),
    regexpressions_parameters{k} = strcat('\<',strtrim(parameterNames{k}),'\>');
    parameterValuesStrings{k} = num2str(SBstructure.parameters(k).value);
end
for k = 1:length(SBstructure.reactions),
    evalString = SBstructure.reactions(k).formula;
    evalString = regexprep(evalString,regexpressions_components,'1');
    evalString = regexprep(evalString,regexpressions_parameters,parameterValuesStrings);
    kineticParameters(k) = eval(evalString);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GET INITIAL CONDITIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
initialConditions = [SBstructure.states.initialCondition];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OUTPUT VALUES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
stoichiometricMatrix = N;
kineticParameters = kineticParameters(:);
initialConditions = initialConditions(:);
return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CHECK IF A CERTAIN SPECIES APPEARS AS PRODUCT AND/OR REACTANT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% In the case that a non-boundary species is defined by a rule
% we need to check if it is also present as product or reactant in any
% reaction. In this case the SBML model is not correctly defined and
% reactionSpeciesPresent = 1
function [reactionSpeciesPresent] = checkReactionsForSpecies(speciesName,SBMLmodel)
% Initialize return value
reactionSpeciesPresent = 0;
reactionComponents = {};
for k = 1:length(SBMLmodel.reaction),
    reactionComponents = {reactionComponents{:},SBMLmodel.reaction(k).reactant.species,SBMLmodel.reaction(k).product.species};
end
if ~isempty(strmatch(speciesName,reactionComponents,'exact')),
    reactionSpeciesPresent = 1;
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EXCHANGE STRING IN STRING - WITH CHECK OF ALLOWED CHARACTERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find occurrences of stringOld in fullString and exchange it with
% stringNew. 
function [processedString] = exchangeStringInString(fullString,stringOld,stringNew);
% really simple!!!
exprString = char([double('\<') double(stringOld) double('\>')]);
processedString = regexprep(fullString, exprString, stringNew);
return
