%clear env
clear
%% set up model parameters
params.N = 1; 
params.d1 = 0.25;
params.d2 = 0.25;
params.d3 = 10; 
params.p1 =4;
params.p2 = 320;
params.p3 = 1;
params.maxRate = 30000;
params.p3Basal =  0.08; 
params.p1_basal = 0;
params.p2_basal = 0;
params.KD = 0.3; 
params.KN = 32;
params.Kr1 = 1200;
params.Kr2 = 1200;
params.KI = 2;
params.p3_external = 0;
params.Tmax =500;
params.Dt =0.1;
params.Tspan = 0:params.Dt:params.Tmax;

%% data
data.X = [30 0 0 24] ;

%% set up SMC params
SMCparams.epsilons =  (54:-2:1)'; %set tolerances 
SMCparams.particles = 1000; % set number of particles
SMCparams.Knn = 50; % set number of neighbours for adjusting perturbation kernel
SMCparams.paramN = 9; % set number of parameters to be inferred
SMCparams.minV = [-3 ; -3   ; -3  ; -3; -2  ;  -2  ; -3  ; -3; -5]; %lower bound for inferred parameters
SMCparams.maxV = [ 3 ;  3   ;  3  ;  0;  2  ;   2  ;  3  ;  3;  2]; %upper bound for inferred parameters
SMCparams.Sigma = diag([0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05]); % initial perturbation
SMCparams.prefix = 'results'; % string for naming ouput files
SMCparams.modelparams = params; % append model parameters to the SMC parameter structure



%% run ABC SMC
[particles, W] = my_abc_SMC(...
    @model_propoasal, ...
    SMCparams, ...
    @model_driver,...
    @model_distance, ...
    data ) ;