function [data_info] = data_artificial(mobj, mi, tv, varargin)
% use simbiology model object to generate artificial data. 
% 
% Dosing and measurement strategy are defined by the mcmc_info struct (mi),
% and the time points defined by the vector tv. Additional 
% name-value pair arguments can be specified. 
%
% REQUIRED INPUTS 
% mobj: simbiology model object. 
% 
% mi: mcmc_info struct. Type 'help mcmc_info' in the command line to read more 
% about this. 
% 
% tv: a vector of timepoints which the data in the output data array will 
% correspond to. These are the poitns at which the model will be simulated to 
% compute the data values. 
% 
% 
% OPTIONAL NAME VALUE PAIR ARGUMENTS - used to populate the data_info struct.
% 
% 'dataInfo': Human readable description of the data. Should be specified using 
% the format that makes it printable with the function fprintf. (so newline 
% characters \n at the ver least, for example.). If not specified, its value is 
% simply 'Artificial Data'
% 
% 'timeUnits': String specifying the units of the time axis. Default is 'seconds'
% 
% 'measuredNames': A 1 x number of measured species cell array of the strings 
% specifying which species are dosed. If not specified, the strings from the 
% mcmc_info struct will be used. When the measured species there is an aggregate 
% of a number of species in the model, a string [xxxxxx ' + ...'] is used, where 
% xxxx is the first string within the list of species whose concentrations get 
% added. 
% 
% 'dataUnits': A cell array of strings specifying the units of the measured 
% species. If not specified, the default is 'nM'. 
%
% 'dosedNames': A 1 x number of dosed species cell array of the strings specifying 
% which species are dosed. These strings should correspond to the strings in 
% the dosedNames field in the mcmc_info struct. If this is not specified, then 
% those very values from the mcmc_info struct are used. 
%
% 'doseUnits': A cell array of strings specifying the units of the 
% dosed species. If not specified, the default is 'nM'. 
%
% 'params', VALUE, where VALUE is a vector of nonnegative parameter values, 
% ordered according to mi.names_unord (the unordered array of parameter and 
% species names. )
% Note that the parameter values are NOT log transformed, ie, they all lie
% in the nonnegative orthant. Parameter 
% values are an optional argument, and if this argument is not specified, 
% then the geometric mean of the rows (each of which has 2 elements) of
% exp(mi.paramranges) will be as the parameter value. 
% 
% 'noise', VALUE; where VALUE is a vector of standard deviations of the gaussian
% noise added to the data, each element corresponding to one measured
% species. 
% 
% 'replicates', VALUE, where VALUE is a positive interger of the number
% of replicates. 
%
% OUTPUTS: This function returns a data_info struct with fields
% ------------------------------------------
%
% 'dataInfo': A human readable text description of the data. If not specified 
% as a name value pair input argument, the string 'Artificial Data' is used. 
% 
% 'timeVector': vector of timepoints, same as tv, a required positional input. 
% 
% 'timeUnits': units of the time vector. Allowed options are: 
% 'seconds', 'minutes', 'hours', 'days', 'weeks'. If no units are specified
% as a name value input, then the units are specified as 'seconds'. 
% 
% 'dataArray': An array contianing the raw data that is generated by simulating 
% the data according to the mcmc_info struct. Typically has dimensions 
% corresponding to timepoints x measured outputs x replicates x doses. 
% 
% 'measuredNames': A 1 x number of measured species cell array of the strings 
% specifying which species are dosed. These are not strings corresponding to
%  the species in the model. Takes from the corresponding name value pair input 
% argument. If not specified, the values are taken from the mcmc_info struct. 
% 
% 'dataUnits': A 1 x number measured species cell array of units corresponding to 
% the raw data in the dataArray. If no units are specified, then nM are used. 
%
% 'dimensionLabels': a 1 by length(size(data_info.dataArray)) cell array of 
% labels for the dimensions of the dataArray. 
%
% 'dosedNames': A 1 x number of dosed species cell array of the strings specifying 
% which species are dosed. These are not strings corresponding to the species 
% in the model. See mcmc_info constructor functions for that. 
% 
% 'dosedVals': A matrix of dose values of size
% # of dosed species by # of dose combinations
% 
% 'doseUnits': A 1 x number of dosed species cell array of strings specifying the 
% units of the dosed species. If no units are 
% 
% ------------------------------------------
%
%

% Copyright (c) 2018, Vipul Singhal, Caltech
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the Software is
% furnished to do so, subject to the following conditions:

% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.

% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
% SOFTWARE.



pdefaults = geomean(exp(mi.paramranges), 2); % !TODO check this
ms = mi.measuredSpecies;
noisedefaults = zeros(length(ms), 1);
dimensionLabels = {'time points', 'measured species', 'replicates', 'doses'};
p = inputParser;

addParameter(p, 'replicates', 1, @isnumeric);  % 
addParameter(p, 'params', pdefaults, @isnumeric);
addParameter(p, 'noise', noisedefaults, @isnumeric);  % the default is 0
addParameter(p, 'timeUnits', 'seconds', @ischar);  
addParameter(p, 'dataUnits', 'nM');  
addParameter(p, 'dosedVals', mi.dosedVals);  
addParameter(p, 'dosedNames', mi.dosedNames);  
addParameter(p, 'doseUnits', 'nM');  
addParameter(p, 'params', pdefaults);
addParameter(p, 'dataInfo', ['Artificial Data']);
addParameter(p, 'measuredNames', mi.measuredSpecies); % default: empty strings
addParameter(p, 'dimensionLabels', dimensionLabels);  

parse(p, varargin{:});
p = p.Results;

configsetObj = getconfigset(mobj, 'active'); 
set(configsetObj, 'StopTime', tv(end));

noisevec = p.noise;
% initialize the data array
da = zeros(length(tv), length(mi.measuredSpecies),...
 p.replicates, size(mi.dosedVals, 2));
ms = mi.measuredSpecies;
% for each dose, simulate the model

dv = mi.dosedVals;

% set parameters and species initial concentrations in the model
paramnames = mi.names_unord;
paramvals = p.params;
for i = 1:length(paramnames)
    p1 = sbioselect(mobj.parameters, 'Name',  paramnames{i});
    if ~isempty(p1)
        set(p1, 'Value', paramvals(i);
    end
end

for i = 1:length(paramnames)
    s1 = sbioselect(mobj.species, 'Name', paramnames{i});
    if ~isempty(s1)
        set(s1, 'InitialValue', paramvals(i))
    end
end

% set dose values, simulate model, and populate output data array. 
for dID = 1:size(dv, 2)
	% set the dose value using the mcmc_info struct
    for i = 1:length(p.dosedNames)
        s1 = sbioselect(mobj.species, 'Name', p.dosedNames{i});
        if ~isempty(s1)
            set(s1, 'InitialValue', p.dosedVals(i))
        end
    end

    % simulate the model. 
	sd = sbiosimulate(mobj);
    sd = resample(sd, tv); 
        for msID = 1:length(ms)
            measuredSpecies = ms{msID};
            spSD = selectbyname(sd, measuredSpecies);
            summed_trajectories = sum(spSD.Data, 2);
   			 %  add noise if needed. 
   			for rID = 1:p.replicates
            	da(:, msID, rID, dID) = ...
                    summed_trajectories + noisevec(msID)*randn(length(tv), 1);
            end
        end
end

% make the data_info struct
di = struct('dataInfo', {p.dataInfo}, ...
            'timeVector' {tv}, ...
            'timeUnits', {p.timeUnits}
            'dataArray' {da},...
            'measuredNames', {p.measuredNames},...
            'dataUnits', {p.dataUnits}
            'dimensionLabels', {p.dimensionLabels}, ...
            'dosedNames' {p.dosedNames},...
            'dosedVals', {p.dosedVals}, ...
            'doseUnits', {p.doseUnits})
