% Function that takes care of running the model and producting the model
% prediction that will be compared with the data. 
% params - a structure with the model parameters to be used
% returns retVal - a vector with the model prediction. Make sure that the
% returned vector has the same dimensions as the data

function retVal = model_driver( params) 
 
    opt = odeset('MaxStep',1);
    Tspan = params.Tspan;
    %% run with diestrus parameters
    X0 = model_IC(params,1);
    [T, Y] = ode23s(@(t,y) modelPopulation_v103(t,y,params), Tspan , X0, opt);
    X0 = Y(end,:);
    %calculate period and duty cycle
    ft = (round(size(Y, 1)/2)):size(Y, 1);
    timecourse_tmp = mean(Y(ft,3:3:end),2) ;
    ac = xcorr(timecourse_tmp,timecourse_tmp);  
    [~,loc] = findpeaks(ac./max(ac), 'MinPeakHeight',0.1);
    DT = Tspan(2)-Tspan(1);

    if isempty(max(diff(loc)*DT))        
                retVal(1, 1) = 0 ;
                retVal(2,1) = 1;
           
    else
                retVal(1,1) = max(diff(loc)*DT) ; 
                retVal(2,1) = mean(timecourse_tmp > ((max(timecourse_tmp) -min(timecourse_tmp))/2));
           
    end
    
    traj(1).T = T;
    traj(1).Y = Y;
    
    %% run with estrus parameters
    paramstmp = params;
    paramstmp.fitparams(1:4) = paramstmp.fitparams(1:4) + paramstmp.fitparams(5:8);
     X0tmp = model_IC(params,1);
    [T, Y] = ode23s(@(t,y) modelPopulation_v103(t,y,paramstmp), Tspan , X0tmp, opt);
    X0tmp = Y(end,:);
    
    ft = (round(size(Y, 1)/2)):size(Y, 1);
    timecourse_tmp = mean(Y(ft,3:3:end),2) ;
    ac = xcorr(timecourse_tmp,timecourse_tmp) ;

   
   [~,loc] = findpeaks(ac./max(ac), 'MinPeakHeight',0.1);
    
    if isempty(max(diff(loc)*DT))        
                retVal(1,2) = 0 ;
                retVal(2,2) = 1;
      
    else
                retVal(1,2) = max(diff(loc)*DT) ; 
                retVal(2,2) = mean(timecourse_tmp > ((max(timecourse_tmp) -min(timecourse_tmp))/2));
      
    end
    
    
    traj(2).T = T;
    traj(2).Y = Y;
    
    %% run diestrus parameters + 5Hz extenal stimulation
    params.p3_external = 5;
    [T, Y] = ode23s(@(t,y) modelPopulation_v103(t,y,params), Tspan , X0, opt);

    %calculate period and duty cycle
    ft = (round(size(Y, 1)/2)):size(Y, 1);
    timecourse_tmp = mean(Y(ft,3:3:end),2) ;
    ac = xcorr(timecourse_tmp,timecourse_tmp);  
    [~,loc] = findpeaks(ac./max(ac), 'MinPeakHeight',0.1);
    DT = Tspan(2)-Tspan(1);

    if isempty(max(diff(loc)*DT))        
                retVal(1,3) = 0 ;
                retVal(2,3) = 1;
         
    else
                retVal(1,3) = max(diff(loc)*DT) ; 
                retVal(2,3) = mean(timecourse_tmp > ((max(timecourse_tmp) -min(timecourse_tmp))/2));
        
    end
    
    traj(3).T = T;
    traj(3).Y = Y;
    
    %% run estrus parameters + 5Hz external stimulations
    paramstmp.p3_external = 5;
    [T, Y] = ode23s(@(t,y) modelPopulation_v103(t,y,paramstmp), Tspan , X0tmp, opt);

    %calculate period and duty cycle
    ft = (round(size(Y, 1)/2)):size(Y, 1);
    timecourse_tmp = mean(Y(ft,3:3:end),2) ;
    ac = xcorr(timecourse_tmp,timecourse_tmp);  
    [~,loc] = findpeaks(ac./max(ac), 'MinPeakHeight',0.1);
    DT = Tspan(2)-Tspan(1);

    if isempty(max(diff(loc)*DT))        
                retVal(1,4) = 0 ;
                retVal(2,4) = 1;
      
    else
                retVal(1,4) = max(diff(loc)*DT) ; 
                retVal(2,4) = mean(timecourse_tmp > ((max(timecourse_tmp) -min(timecourse_tmp))/2));   
    end
 
  
    traj(4).T = T;
    traj(4).Y = Y;
end

%ydata = [1/19.257];% 0.149]