% ccc
distcomp.feature( 'LocalUseMpiexec', false);
format long g
parpool(24)
Topt=3; % This should be same with the OED part. You can take it from its saved mat file. Normally it is 1 seconds
dt=1e-3;
Nopt=Topt/dt+1;
ta=0:dt:Topt;
Nita=100;
% The true parameters

be=1/0.020;
bi=1/0.040;
we=1;
wi=0.7;
wee=1.2;
wei=2.0;
wie=0.70;
wii=0.400;

% Auxiliary Parameters:

Fmaxe=100;
Fmaxi=50;
vce=70;
vci=35;
ae=0.04;
ai=0.04;
% qx=[Fmaxe Fmaxi vce vci ae ai];
tht=[be bi we wi wee wei wie wii Fmaxe Fmaxi vce vci ae ai]';


Amax=100;

fbase=15/Topt;
Nu=5;
wbase=2*pi*fbase;
WE=wbase*(1:Nu)';
AE=Amax.*ones(Nu,1);
lbp=zeros(length(tht),1);
ubp=2*tht;
for is=1:10
for oi=1:Nita
    PE=unifrnd(-pi,pi,Nu,1);
    ve=0;
    vi=0;
    for i=1:Nopt
        t=ta(i);
        u=sum(AE.*cos(WE.*t+PE));
        ge=Fmaxe/(exp(ae*(vce - ve)) + 1);
        gev=(Fmaxe*ae*exp(ae*(vce - ve)))/(exp(ae*(vce - ve)) + 1)^2;
        gevv=(2*Fmaxe*ae^2*exp(2*ae*(vce - ve)))/(exp(ae*(vce - ve)) + 1)^3 - (Fmaxe*ae^2*exp(ae*(vce - ve)))/(exp(ae*(vce - ve)) + 1)^2;
        gi=Fmaxi/(exp(ai*(vci - vi)) + 1);
        giv=(Fmaxi*ai*exp(ai*(vci - vi)))/(exp(ai*(vci - vi)) + 1)^2;
        givv=(2*Fmaxi*ai^2*exp(2*ai*(vci - vi)))/(exp(ai*(vci - vi)) + 1)^3 - (Fmaxi*ai^2*exp(ai*(vci - vi)))/(exp(ai*(vci - vi)) + 1)^2;
        
        vedot=be*(-ve+ge*wee-gi*wei+we*u);
        vidot=bi*(-vi-gi*wii+ge*wie+wi*u);
        re=ge;
        ve=ve+vedot*dt;
        vi=vi+vidot*dt;
        rea(i)=re;
    end
    rta=rea*dt;
    uaar=unifrnd(0,1,1,Nopt);
    shi=rta>uaar;
    shinz=find(shi);
    Tspks=ta(shinz);
    uspks(oi).spks=Tspks';
    uspks(oi).A=AE;
    uspks(oi).P=PE;
    uspks(oi).W=WE;    
end
Nin=16;
Np=length(tht);
fs1='%d %d %d %f %f ';
fs2=strjoin(repelem({'%f'},Np));
fs3='\n';
fs=strcat(fs1,fs2,fs3);
fsd=strcat(fs2,fs3);
sa=num2str(Amax);
sw=num2str(wbase);
sn=num2str(Nita);
snu=num2str(Nu);
fname=strcat('results_','nit_',sn,'_amax_',sa,'_wbase_',sw,'_nu_',snu,'.txt');
fnameb=strcat('results_','nit_',sn,'_amax_',sa,'_wbase_',sw,'_nu_',snu,'.dat');
parfor k=1:Nin;
    th0=unifrnd(lbp,ubp);
    foptions=optimset('Algorithm','interior-point', 'Diagnostics','on', 'Display','iter', 'GradObj', 'off', 'Hessian','off','DerivativeCheck','off','TolFun',1e-6,'TolX',eps,'TolCon',1e-6,'MaxFunEvals',30000,'FunValCheck','off','MaxIter',500);
    [theta,fvalp,exitflagp,outputp,lambdap,gradp,hessianp] = fmincon(@(th)jhujed_obj(th,uspks,Topt,dt,Nita),th0,[],[],[],[],lbp,ubp,[],foptions);
    fid=fopen(fname,'a+');
    fidb=fopen(fnameb,'a+');
    % THT(k,:)=theta';
    % stht(is).data=THT;
    fprintf(fid,fs,is,k,exitflagp,fvalp,norm(gradp),theta');
    fprintf(fidb,fsd,theta');
    fclose(fid);
    fclose(fidb);
end
end
delete(gcp)