//******************************************************
//	Programmer: Steve Martell
//	Project Name: Blue crab CMS model		
//	Date:
//	Version:  4
//	Comments:  Modified by Tom Miller to include N=N+R
//	
//******************************************************/
DATA_SECTION
	init_number m;

	init_number fm;			//Fishery month
	!!fm=fm/12;

	//Read in Catch data
	init_int scyr;
	init_int ecyr;
	init_vector ct(scyr,ecyr);
	int nobs;
	!!nobs=ecyr-scyr;
	//!!ct*=0;
	

	// read in estimated exploitation fractions
	init_int f_switch;			// f_switch=1 use adults, f_switch=2 use adutls and recruits
	init_int mu_switch;			// switch to include exploitation fraction in fvec (1=in, 0=out)
	init_int smuyr;
	init_int emuyr;
	init_vector mu(smuyr,emuyr);
	int nmu;				// number of obs of exploitaiton fraction
	
	!! nmu=(emuyr-smuyr);
	!!cout<<"f_switch "<<f_switch<<endl;
	!!cout<<"Mu_switch "<<mu_switch<<endl;

	//read in paired survey data
	init_int nsets;                   // Number of surveys
	init_vector relsel(1,nsets);
	init_ivector sSyr(1,nsets);
	init_ivector eSyr(1,nsets);
	init_vector sm(1,nsets);			
	
	init_matrix yt_r(1,nsets,sSyr,eSyr);
	init_matrix yt_a(1,nsets,sSyr,eSyr);
	!!sm=sm/12;
	//!!cout<<"Survey start year"<<endl<<sSyr<<endl;
	//!!cout<<"End start year"<<endl<<eSyr<<endl;

	//read in unpaired adult survey data
	init_int ny;
	init_ivector syt(1,ny);
	init_ivector eyt(1,ny);
	init_vector usm(1,ny);
	init_vector mode_pry(1,ny);
	init_vector cv_pry(1,ny);
	init_matrix yt(1,ny,syt,eyt);
	!!usm=usm/12;

	



PARAMETER_SECTION
	init_bounded_number log_no(0.0001,10,1);
	//init_bounded_vector sel(1,nsets,0.0001,2,2);
	init_bounded_vector kappa(1,nsets,0.0001,3,1);		//total variance
	init_bounded_vector log_rt(scyr,ecyr+1,0.0001,10,1);
	init_bounded_vector nu(scyr+1,ecyr,-15,15,1);		//
	init_bounded_number std_mu(0,100,1);
	!!std_mu=500;
	init_bounded_vector p(1,nsets,0.0001,1,1);		// proportion of total var that is measurement 
	!! p=0.5;
	
	
	!!kappa=0.065;
	
	!!log_no=log(max(ct)/0.007);
	
	!!log_rt=log(exp(log_no)/exp(-m)-exp(log_no));
	
	
	sdreport_number sd_no;
	
	objective_function_value f;
	
	number tau;
	
	vector na(1,nsets);
	vector nr(1,nsets);
	vector ny(1,nsets);
	
	vector q_a(1,nsets);		     // selectivity for adults in survey
	vector q_r(1,nsets);		     // selectivity for recruits in survey
	
	vector nt(scyr,ecyr+1);              // fully recruited abundance
	vector rt(scyr,ecyr+1);              // recruit abundance 
	vector ft(scyr,ecyr);                // exploitation fraction
	vector ft2(scyr,ecyr);		     // exploitation fraction based on n+r
	vector mu_dev(smuyr,emuyr);          // expoitation fraction deviations
	vector Nt(scyr,ecyr+1);              // total population abundance
	
	matrix epsilon(1,nsets,sSyr,eSyr);      //adult observation errors
	matrix delta(1,nsets,sSyr,eSyr);	//recruit observation errors
	matrix p_nu(1,nsets,sSyr+1,eSyr);	//predicted process errors
	
	
PROCEDURE_SECTION
	//--------MAIN--------
	
	pop_dynamics();
	//cout<<"OK after pop_dynamics"<<endl;
	exp_dev();
	//cout<<"OK after exploitation deviations"<<endl;
	observation_errors();
	//cout<<"OK after observation_errors"<<endl;
	process_errors();
	//cout<<"OK after process_errors"<<endl;
	calc_objective_function3();
	//cout<<"OK after calc_objective_function3"<<endl;
	
	//____________________
	
FUNCTION pop_dynamics
	int i;
	nt.initialize();
	nt(scyr)=exp(log_no);
	rt=exp(log_rt);
	
	//process error var.
	tau=min(elem_prod((1.-p),kappa));
	sd_no=exp(tau);
	//cout<<"tau = "<<tau<<endl;
	
	for(i=scyr;i<=ecyr;i++)
	{
		nt(i+1)=((nt(i)+rt(i))*exp(-m*fm)-ct(i))*exp(-m*(1.-fm));
		
		if(i>scyr) {
			nt(i+1)*=exp(sqrt(tau)*nu(i));
		}
		
		Nt(i)=nt(i)+rt(i);
	}
	switch(f_switch)
		{
	           case 1: 			// use just adults
	           		ft=elem_div(ct,nt(scyr,ecyr));
	           break;
		   case	2:			// use both adults and recruits
				ft=elem_div(ct,Nt(scyr,ecyr));
		   break;			
		}
//_________________________________________
FUNCTION exp_dev
	int i;
	mu_dev.initialize();
	
	
	
	for(i=smuyr;i<=emuyr;i++)
	{
		
		mu_dev(i)=log(mu(i)-ft(i));
	}
	
		
FUNCTION dvariable get_q_MLE(dvector& y, dvar_vector& n)
	{
		//compute conditional MLE for q_a
		int i;
		dvariable q;
		dvariable xy;
		dvariable x2;
		
		xy=0.;
		x2=0.;
		
		for(i=y.indexmin();i<=y.indexmax();i++)
		{
			if(y(i)>0){
				xy+=y(i)*n(i);                // maintain running totals
				x2+=n(i)*n(i);
			}
		}
		q=xy/x2;
		return(q);
	}

FUNCTION observation_errors
	int i,j;
	//Modify to accomodate missing data  DONE
	
	epsilon.initialize();
	delta.initialize();
	q_a.initialize();
	q_r.initialize();
	
	na.initialize();
	nr.initialize();
	
 	for(j=1;j<=nsets;j++)
	{
		
		//MLE for q_a
		q_a(j)=get_q_MLE(yt_a(j),nt(sSyr(j),eSyr(j)));
		q_r(j)=q_a(j)*relsel(j);
		
		//q_r(j)=q_a(j)*relsel(j);
		//q_r(j)=q_a(j)*sel(j);
		
		
		
		
		
		//do adult indices first 
		
		
		for(i=sSyr(j);i<=eSyr(j);i++)
		
		{
			if(yt_a(j,i)>0){
				epsilon(j,i)=log(yt_a(j,i))-log(q_a(j))-log(nt(i)*exp(-m*sm(j)));
				na(j)+=1.;	                                                     // sample size
			}
		}
		
		//now do juvenile indices
		
		
		for(i=sSyr(j);i<=eSyr(j);i++)
		
		{
			if(yt_r(j,i)>0){
				delta(j,i)=log(yt_r(j,i))-log(q_r(j))-log(rt(i)*exp(-m*sm(j)));
				nr(j)+=1.;                                                           // sample size
			}
			
			//if(yt_a(j,i)>0)epsilon(j,i)-=sum(epsilon(j))/na(j);
		}
		
	}
	
	//cout<<"e = "<<sum(epsilon)<<endl;
	//cout<<"d = "<<sum(delta)<<endl;
	
	
FUNCTION process_errors
	//For each index, calculate the expected process errors 
	//based on relative abundance indices and pop dy model.
	
	//Modify to accomodate missing data  DONE
	
	int i,j;
	p_nu.initialize();
	ny.initialize();
	
	
	
	for(j=1;j<=nsets;j++)
	
	{
		for(i=sSyr(j);i<eSyr(j);i++)
		{
			if(yt_a(j,i+1)>0&&yt_a(j,i)>0&&yt_r(j,i)>0){
				p_nu(j,i+1)=log(yt_a(j,i+1)/q_a(j));
				p_nu(j,i+1)-=log(((nt(i)+rt(i))*exp(-m*fm)-ct(i))*exp(-m*(1.-fm)));
				ny(j)+=1.;
			}
		}
		//cout<<p_nu(j)<<endl;
		//cout<<"p_nu = "<<sum(p_nu)<<endl;
	}
	//cout<<"p_nu = "<<sum(p_nu)<<endl;
	
FUNCTION calc_objective_function3	
	//Modify to accomodate missing data DONE(changes in sample size na)

	int j;
	dvar_vector fvec(1,nsets+1);
	fvec.initialize();
	
	
	
	for(j=1;j<=nsets;j++)
	
	{
			fvec[j]=0.5*(na(j)+nr(j)+ny(j))*log(kappa(j))
						+0.5*(na(j)+nr(j))*log(p(j))
						+0.5*ny(j)*log((1-p(j)))
						+1./(2.*kappa(j))
						*(norm2(epsilon(j))/p(j)+norm2(delta(j))/p(j)+norm2(p_nu(j))/(1.-p(j)));
	}
	
	if(mu_switch==1) {
				fvec[nsets+1] =(nmu/2)*norm2(mu_dev);
			}	
	cout<<"fvec"<<fvec<<endl;	
	f=sum(fvec);
		
	
REPORT_SECTION
	int i;
	report<<"Nt"<<endl<<nt<<endl;
	report<<"Rt"<<endl<<rt<<endl;
	report<<"nsets"<<endl<<nsets<<endl;
	report<<"sSyr & eSyr"<<endl<<sSyr<<endl<<eSyr<<endl;
	report<<"Adult Survey"<<endl;
	for(i=1;i<=nsets;i++){
		report<<yt_a(i)/q_a(i)<<endl;
	}
	report<<"Recruit Survey"<<endl;
	for(i=1;i<=nsets;i++){
			report<<yt_r(i)/q_r(i)<<endl;
	}
	
	//fishing mortality rate
	report<<"Ft"<<endl<<ft<<endl;
	report<<"Ct"<<endl<<ct<<endl;
	report<<"Ft2"<<endl<<ft2<<endl;
	report<<"Exp. Deviations"<<endl<<mu_dev<<endl;

	//residuals
	report<<"delta"<<endl<<(delta)<<endl;
	report<<"epsilon"<<endl<<(epsilon)<<endl;
	report<<"p_nu"<<endl<<(p_nu)<<endl;
	report<<"nu"<<endl<<nu<<endl;
	report<<"predicted ratio"<<p<<endl;
	report<<"q_a"<<endl<<q_a<<endl;
	report<<"tau = "<<tau<<endl;
	report<<"std_mu"<<std_mu<<endl;