%%
%%STEP 1: Parameter factors generation (crea una matriz aleatoria de
%%parámetros que siguen una distribución normal)
%                  1     2      3     4        5       6       7      
parameter_names={'INa','INaL','Ito','ICaL', 'ICaNa', 'ICaK', 'IKr', ...
                'IKs','IK1','INCX','INaK','INab', 'IKb', 'IpCa', 'ICab'};
%                 8     9     10     11     12     13      14     15   
%                 

n_parameters = length(parameter_names);

trials = 1000; % number of trials (rows); parameters (columns)
mu=1;
sigma = 0.2; % standar deviation
rng('shuffle');
p=normrnd(mu, sigma, [ trials,n_parameters]); 
%% Simulación de la population sin fármaco
CL=1000; %pacing cycle length in ms
beats=500; %number of beats in the simulation
%%%%%%%%%% Las siguientes lineas sirven para inicializar los 'workers' y hacer
% simulaciones en paralelo (1 por worker) y que terminen antes
tic
myCluster=parcluster('local');
poolobj = gcp('nocreate'); % If no pool, do not create new one.
if isempty(poolobj)
    poolobj = parpool(myCluster.NumWorkers);
end
poolobj.IdleTimeout= 120;
disp(datetime)
disp(['Usando ' num2str(poolobj.NumWorkers) 'workers'])
%%%%%%%%%%%%%%
% Lanza la simulaciones
parfor k=1:trials
    p_SA = p(k,:); % 15 parameters
    [time{k},X{k},IsJs{k}] = runORdmD(p_SA, 3)
end
%%% guarda la simulación y borra los elementos no necesarios
delete(poolobj);
toc
disp('Esperando a guardar');
s.name='PopulationORdmD_raw';
s.params=p;
s.result=[(time'),(X'),(IsJs')];
save(['PopulationORdmD_raw.mat'], 's')
clearvars -except p trials
%% Simulación de la population con un bloqueo del 90% de IKs (necesario para el calibrado)
CL=1000; %pacing cycle length in ms
beats=500; %number of beats in the simulation
tic
myCluster=parcluster('local');
% myCluster.NumWorkers=64;
poolobj = gcp('nocreate'); % If no pool, do not create new one.
if isempty(poolobj)
    poolobj = parpool(myCluster.NumWorkers);
end
poolobj.IdleTimeout= 120;
disp(datetime)
disp(['Usando ' num2str(poolobj.NumWorkers)])
p_IKs90=p;
p_IKs90(:,8)=p_IKs90(:,8)*0.1; % IKs es el 8 parámetro
parfor k=1:trials
    p_SA = p_IKs90(k,:); % 15 parameters
    [time{k},X{k},IsJs{k}] = runORdmD(p_SA, 3)
end
delete(poolobj);
tocclc
disp('Esperando a guardar');
s.name='PopulationORdmD_raw_IKs90';
s.params=p_IKs90;
s.result=[(time'),(X'),(IsJs')];
save(['PopulationORdmD_raw_IKs90.mat'], 's')
clearvars -except p trials
%% Simulación de la population con un bloqueo del 70% de IKr (necesario para el calibrado)
CL=1000; %pacing cycle length in ms
beats=500; %number of beats in the simulation
tic
myCluster=parcluster('local');
poolobj = gcp('nocreate'); % If no pool, do not create new one.
if isempty(poolobj)
    poolobj = parpool(myCluster.NumWorkers);
end
poolobj.IdleTimeout= 120;
disp(datetime)
disp(['Usando ' num2str(poolobj.NumWorkers)])
p_IKr70=p;
p_IKr70(:,7)=p_IKr70(:,7)*0.3;
parfor k=1:trials
    p_SA = p_IKr70(k,:); % 15 parameters
    [time{k},X{k},IsJs{k}] = runORdmD(p_SA, 3)
end
delete(poolobj);
toc
disp('Esperando a guardar');
s.name='PopulationORdmD_raw_IKr70';
s.params=p_IKr70;
s.result=[(time'),(X'),(IsJs')];
save(['PopulationORdmD_raw_IKr70.mat'], 's')
clearvars -except p trials
%% Simulación de la population con un bloqueo del 50% de IK1 (necesario para el calibrado)
CL=1000; %pacing cycle length in ms
beats=500; %number of beats in the simulation
tic
myCluster=parcluster('local');
% myCluster.NumWorkers=64;
poolobj = gcp('nocreate'); % If no pool, do not create new one.
if isempty(poolobj)
    poolobj = parpool(myCluster.NumWorkers);
end
poolobj.IdleTimeout= 120;
disp(datetime)
disp(['Usando ' num2str(poolobj.NumWorkers)])
p_IK1_50=p;
p_IK1_50(:,9)=p_IK1_50(:,9)*0.5;
parfor k=1:trials
    p_SA = p_IK1_50(k,:); % 15 parameters
    [time{k},X{k},IsJs{k}] = runORdmD(p_SA, 3)
end
delete(poolobj);
toc
disp('Esperando a guardar');
s.name='PopulationORdmD_raw_IK1_50';
s.params=p_IK1_50;
s.result=[(time'),(X'),(IsJs')];
save(['PopulationORdmD_raw_IK1_50.mat'], 's')
clearvars -except p