* Fisher macros;

%macro FirstLast(SeasonStart,SeasonEnd,nSeasons,nv,nd,status);

* Finds first and last year visited;
  &SeasonStart = .;
  &SeasonEnd = .;
  do _i=1 to &nSeasons;
     if &SeasonStart = . and &nv{_i} > 0 then &SeasonStart = _i;
	 if &SeasonEnd = . and &nv{&nSeasons-_i+1} > 0 
        then &SeasonEnd = &nSeasons - _i + 1;
	 if &nv{_i} = 0 then &status{_i} = .;
	 else if &nd{_i} = 0 then &status{_i} = 0;
	 else &status{_i} = 1;
  end;
  if &SeasonStart = . then abort;

%mend FirstLast;

%macro Zeros(nMissing,missingYears,SeasonStart,SeasonEnd,nd);

* Finds years with zero detections or no visits;
  &nMissing = 0;
  do _i=&SeasonStart to &SeasonEnd;
     if &nd{_i} = 0 then do;
        &nMissing = &nMissing + 1;
        &missingYears{&nMissing} = _i;
     end;
  end;

%mend Zeros;

%macro ProcessHistory(prob,nSeasons,st,status,nMissing,missingYears,
  SeasonStart,SeasonEnd,Psi,pHistory,epsilon,gamma);

* Determines likelihood for a particular site history;

* Process history;
  if &nMissing > 0 then do;
	 * Copy status to st;
     do ii=1 to &nSeasons; &st{ii} = &status{ii}; end;
     k = 2**&nMissing - 1;
     &prob = 0;
	 do j=0 to k;
        do m = 1 to &nMissing;
		   &st{&missingYears{m}} = int(j/2**(m-1))-2*int(j/2**m);
		end;
        _pr = 0;
        do i=&SeasonStart to &SeasonEnd;
		   if i=&SeasonStart then do;
              if &st{i} = 1 
                 then _pr = &Psi * &pHistory;
			     else _pr = 1 - Ψ
		   end;
		   else do;
             if &st{i-1} = 0 and &st{i} = 0 then _pr = _pr * (1-&gamma);
		     if &st{i-1} = 0 and &st{i} = 1 then _pr = _pr * &gamma*&pHistory;
			 if &st{i-1} = 1 and &st{i} = 0 then _pr = _pr * ε
             if &st{i-1} = 1 and &st{i} = 1 then _pr = _pr * (1-&epsilon)*&pHistory;
		   end;
        end;
        &prob = &prob + _pr;
	 end;	 
  end;
  else do;
     do i=&SeasonStart to &SeasonEnd;
	    if i=&SeasonStart then &prob = &Psi * &pHistory;
		                  else &prob = &prob * (1-&epsilon) * &pHistory; 
	 end;
  end;

%mend ProcessHistory;

%macro p2FromP1eg(p2,p1,e,g);

* Computes p2 from p1, e (extinction probability), and gamma
  (colonization probability);

* p2 must be a variable name but p1, e, and g can be variable
  names, expressions, or numbers;

  &p2 = (&p1)*(1-(&e)) + (1-(&p1))*(&g);

* Check for numerical slippage past boundaries;
  if &p2 < 0 then &p2 = 0;
  if &p2 > 1 then &p2 = 1;

%mend p2FromP1eg;

%macro phiFromp1p2e(phi,p1,p2,epsilon);

* Estimates phi from p1, p2, and g;

* phi must be a variable name but p1, p2, and e can be variable
  names, expressions or numbers;

  *if (&p2) > (1-(&epsilon)) then do;
  *   if ((&p1)+(&p2)) > 1 then &phi = (&p1)*(1-(&epsilon)-(&p2))/((&p1)*(&p2)-(&p1)-(&p2)+1);
  *                    else &phi = (1-(&epsilon)-(&p2))/(&p2);
  *end;
  *else do;
  *   if (&p1) <= (&p2)
  *      then &phi = (1-(&epsilon)-(&p2))/(1-(&p2));
  *      else &phi = (&p1)*(1-(&epsilon)-(&p2))/((&p2)-(&p1)*(&p2));
  *end;
  if (&p2) > (1-(&epsilon)) then do;
     if ((&p1)+(&p2)) > 1 
        then &phi = (&p1)*(1-(&epsilon)-(&p2))/((&p1)*(&p2)-(&p1)-(&p2)+1);
        else &phi = (1-(&epsilon)-(&p2))/(&p2);
  end;
  else do;
     if (&p1) <= (&p2)
        then &phi = (1-(&epsilon)-(&p2))/(1-(&p2));
        else &phi = (&p1)*(1-(&epsilon)-(&p2))/((&p2)-(&p1)*(&p2));
  end;
  if &phi < -1 then &phi = -1;
  if &phi >  1 then &phi = 1;

%mend phiFromp1p2e;

%macro egFromp1p2phi(e,g,p1,p2,phi);

* Determines e (extinction probability) and g (colonization
  probability) from p1,p2, and phi;

* e and g must be variable names but p1, p2, and phi can be
  variable names, expressions or numbers;

  if (&phi) < 0 then do;
	 if ((&p1)+(&p2)) <= 1 
        then &e = 1-(&p2)*(1+(&phi));
	    else &e = (1-(&p2))*((&p1)-(1-(&p1))*(&phi))/(&p1);
  end;
  else do;
     if (&p2) < (&p1)
        then &e = 1 - (&p2) - (&phi)*((&p2)-(&p2)*(&p1))/(&p1); 
        else &e = (1-(&p2))*(1-(&phi));
  end;

* Check for numerical accuracy slipping past boundaries;
  if &e < 0 then &e = 0;
  if &e > 1 then &e = 1;
  &g = ((&p2)-(&p1)*(1-(&e)))/(1-(&p1));
  if &g < 0 then &g = 0;
  if &g > 1 then &g = 1;

%mend egFromp1p2phi;

%macro Analysis(data,nSeasons,psiType,phiType,pType,psi=0.3,phi=0,p=0.3,a=0.5,b=0,pA=0.3,pB=0.7,palpha=0.5,gconv=1E-12);

* General models;

* Set model ID;
  %let id = &psiType.&phiType.&pType;

* Initialize output datasets;
  data pe&id;
  data fs&id;
  %if (pType = TwoGroup) %then %do; data ae&id; %end;
  
* Find maximum likelihood estimates;
proc nlmixed data=&data technique=nrridg gconv=&gconv maxiter=5000;
  by zone;
  parms 
    %if(&psiType=dot) %then _psi1 = &psi;
    %if(&psiType=Year) %then _psi1-_psi&nSeasons = &psi;
    %if(&PsiType=logit) %then _a = &a _b = &b;

    %if(&phiType=dot) %then _phi1 = &phi;
    %if(&phiType=Year) %then _phi1-_phi%eval(&nSeasons-1) = 0;

    %if(&pType=dot) %then _p1 = &p;
    %if(&pTYpe=Year) %then _p1-_p&nSeasons = &p;
    %if(&pType=TwoGroup) %then %do; _pA1 = &pA _pB1 = &pB _alpha1 = &palpha  %end;
	%if(&pTYpe=YearTwoGroup) %then %do; _pA1-_pA&nSeasons = &pA _pB1-_pB&nSeasons = &pB _alpha1-_alpha&nSeasons = &palpha %end;
	;

  bounds 0 <= 
    %if(&psiType=dot) %then _Psi1;
    %if(&psiType=Year) %then _Psi1-_Psi&nSeasons;

    %if(&pType=dot) %then _p1;
    %if(&pTYpe=Year) %then _p1-_p&nSeasons;
    %if(&pType=TwoGroup) %then %do; _pA1 _pB1 _alpha1 %end;
	%if(&pTYpe=YearTwoGroup) %then %do; _pA1-_pA&nSeasons _pB1-_pB&nSeasons _alpha1-_alpha&nSeasons %end; <= 1

    %if(&phiType=dot) %then , -1 <= _phi1 <= 1;
    %if(&phiType=Year) %then , -1 <= _phi1-_phi%eval(&nSeasons-1) <= 1;;

  * Count arrays;
  array nv{&nSeasons} nv1-nv&nSeasons;
  array nd{&nSeasons} nd1-nd&nSeasons;

  * Status arrays;
  array _status{&nSeasons} _status1-_status&nSeasons;
  array _missingSeasons{&nSeasons} _missingSeasons1-_missingSeasons&nSeasons;
  array _st{&nSeasons} _st1-_st&nSeasons;

  * Parameter arrays;
  array _epsilon{%eval(&nSeasons-1)} _epsilon1-_epsilon%eval(&nSeasons-1);
  array _gamma{%eval(&nSeasons-1)} _gamma1-_gamma%eval(&nSeasons-1);
  array _Psi{&nSeasons} _Psi1-_Psi&nSeasons;
  array _pHistory{&nSeasons} _pHistory1-_pHistory&nSeasons;
  array _p{&nSeasons} _p1-_p&nSeasons;
  array _phi{%eval(&nSeasons-1)} _phi1-_phi%eval(&nSeasons-1);
  array _pA{&nSeasons} _pA1-_pA&nSeasons;
  array _pB{&nSeasons} _pB1-_pB&nSeasons;
  array _alpha{&nSeasons} _alpha1-_alpha&nSeasons;

  * Initialize values that are not parameters to keep SAS from complaining;
  _pHistory{1} = 0;
  do _i=2 to &nSeasons;
     _epsilon{_i-1} = 0;
	 _gamma{_i-1} = 0;
	 _pHistory{_i} = 0;
  end;

  %if (&pType=dot) %then %do;
      _pA{1} = 0;
      _pB{1} = 0;
      _alpha{1} = 0;
      do _i=2 to &nSeasons;
         _p{_i} = _p{1};
		 _pA{_i} = 0;
		 _pB{_i} = 0;
		 _alpha{_i} = 0;
      end;
  %end;

  if ("&pType"="year") then do;
      do _i=1 to &nSeasons;
	     _pA{_i} = 0;
		 _pB{_i} = 0;
		 _alpha{_i} = 0;
	  end;
  end;

 %if (&pType=TwoGroup) %then %do;
      _p{1} = 0;
      do _i=2 to &nSeasons;
	     _pA{_i} = _pA{1};
		 _pB{_i} = _pB{1};
		 _alpha{_i} = _alpha{1};
	     _p{_i} = 0;
      end;
  %end;
 
  %if (&pType=YearTwoGroup) %then %do;
      do _i=1 to &nSeasons;
	     _p{_i} = 0;
      end;
  %end;

  %if (&psiType=dot) %then %do;
      do _i=2 to &nSeasons;
	     _psi{_i} = _psi{1};
	  end;
  %end;

  %if (&phiType=dot) %then %do;
      do _i=2 to &nSeasons-1;
	     _phi{_i} = _phi{1};
	  end;
  %end;

  %if (&phiType=0) %then %do;
      do _i=1 to &nSeasons-1;
	     _phi{_i} = 0;
	  end;
  %end;

  %if (&psiType=logit) %then %do;
      do _i=1 to &nSeasons;
	     _psi{_i} = 1 - 1/(1+exp(_a+_b*_i));
	  end;
  %end;

  * Determine epsilon and gamma values from psi and phi values;
  do _i = 2 to &nSeasons;
     %egFromp1p2phi(_epsilon{_i-1},_gamma{_i-1},_psi{_i-1},_psi{_i},_phi{_i-1});
  end;

  * Create variable that will always exist (again so SAS will not complain);
  _x = 0;

  * Find first and last year visited;
  %FirstLast(_SeasonStart,_SeasonEnd,&nSeasons,nv,nd,_status)

  * Find years with zero detections or no visits;
  %Zeros(_nMissing,_missingSeasons,_SeasonStart,_SeasonEnd,nd);
 
  * Determine detection probabilities for each history;
  %if (&pType = dot or &pType = Year) %then %do;
  do _i=_SeasonStart to _SeasonEnd;
     _pHistory{_i} = 1;
     if nv{_i} > 0 then do;
	    if nd{_i} = 0 then _pHistory{_i} = (1-_p{_i})**nv{_i};
		else if nd{_i} = nv{_i} then _pHistory{_i} = _p{_i}**nv{_i};
		else _pHistory{_i} = _p{_i}**nd{_i}*(1-_p{_i})**(nv{_i}-nd{_i});
	 end;
  end;
  %end;

  %if (&pType = TwoGroup or &pType = YearTwoGroup) %then %do;
  * Determine detection probabilities for each history;
  do _i=_SeasonStart to _SeasonEnd;
     _pHistory{_i} = 1;
     if nv{_i} > 0 then do;
        if nd{_i} = 0 then 
           _pHistory{_i} = _alpha{_i}*(1-_pA{_i})**nv{_i} +
                         (1-_alpha{_i})*(1-_pB{_i})**nv{_i};
		else if nd{_i} = nv{_i} then 
           _pHistory{_i} = _alpha{_i}*_pA{_i}**nv{_i} + (1-_alpha{_i})*_pB{_i}**nv{_i};
		else _pHistory{_i} = _alpha{_i}*_pA{_i}**nd{_i}*(1-_pA{_i})**(nv{_i}-nd{_i}) +
           (1-_alpha{_i})*_pB{_i}**nd{_i}*(1-_pB{_i})**(nv{_i}-nd{_i});
	 end;
  end;
  %end;

  * Process history;
  %ProcessHistory(_prob,&nSeasons,_st,_status,_nMissing,_missingSeasons,
    _SeasonStart,_SeasonEnd,_Psi{i},_pHistory{i},_epsilon{i-1},_gamma{i-1});

  _logL = log(_prob);
  model _x ~ general(_logL);

  %if (&pType = TwoGroup) %then %do;
  estimate 'p 1 visit' _pA1*_alpha1 + _pB1*(1-_alpha1);
  estimate 'p 5 visits' _alpha1*(1-(1-_pA1)**5) + (1-_alpha1)*(1-(1-_pB1)**5);
  %end;

  %if (&psiType = logit) %then %do;
  estimate 'Average annual change' 
    100*((1 - 1/(1+exp(_a+_b*8))) - (1 - 1/(1+exp(_a+_b*1))))/7;
  estimate 'psi(2002)' 1 - 1/(1+exp(_a+_b*1));
  estimate 'psi(2003)' 1 - 1/(1+exp(_a+_b*2));
  estimate 'psi(2004)' 1 - 1/(1+exp(_a+_b*3));
  estimate 'psi(2005)' 1 - 1/(1+exp(_a+_b*4));
  estimate 'psi(2006)' 1 - 1/(1+exp(_a+_b*5));
  estimate 'psi(2007)' 1 - 1/(1+exp(_a+_b*6));
  estimate 'psi(2008)' 1 - 1/(1+exp(_a+_b*7));
  estimate 'psi(2009)' 1 - 1/(1+exp(_a+_b*8));
  %end;

  %do i=1 %to %eval(&nSeasons-1);
  estimate "extinction probability &i" _epsilon{&i};
  estimate "colonization probability &i" _gamma{&i};
  %end;

  ods output ParameterEstimates=PE&id FitStatistics=FS&id 
    AdditionalEstimates=AE&id;
run;

* Add in model information;
  data ae&id;
    length label $ 255;
	length model $ 255;
	length psiModel $ 255;
    length phiModel $ 255;
    length pModel $ 255;
	set ae&id;
    drop psiModel phiModel pModel;
    if "&psiType" = "dot" then psiModel = "psi(.)";
	if "&psiType" = "Year" then psiModel = "psi(Year)";
	if "&psiType" = "logit" then psiModel = "psi(logit)";
	if "&phiType" = "dot" then phiModel = ",phi(.)";
	if "&phiType" = "Year" then phiModel = ",phi(Year)";
	if "&phiType" = "0" then phiModel = ",phi=0";
	if "&pType" = "dot" then pModel = ",p(.)";
	if "&pType" = "Year" then pModel = ",p(Year)";
	if "&pType" = "TwoGroup" then pModel = ",p(2 groups)";
	if "&pType" = "YearTwoGroup" then pModel = ",p(Year*2 groups)";
    model = trim(psiModel) || trim(phiModel) || trim(pModel);
  run;
  data pe&id;
    length parameter $ 255;
	length model $ 255;
	length psiModel $ 255;
    length phiModel $ 255;
    length pModel $ 255;
	set pe&id;
    drop psiModel phiModel pModel;
    if "&psiType" = "dot" then psiModel = "psi(.)";
	if "&psiType" = "Year" then psiModel = "psi(Year)";
	if "&psiType" = "logit" then psiModel = "psi(logit)";
	if "&phiType" = "dot" then phiModel = ",phi(.)";
	if "&phiType" = "Year" then phiModel = ",phi(Year)";
	if "&phiType" = "0" then phiModel = ",phi=0";
	if "&pType" = "dot" then pModel = ",p(.)";
	if "&pType" = "Year" then pModel = ",p(Year)";
	if "&pType" = "TwoGroup" then pModel = ",p(2 groups)";
	if "&pType" = "YearTwoGroup" then pModel = ",p(Year*2 groups)";
    model = trim(psiModel) || trim(phiModel) || trim(pModel);
  run;
  data fs&id;
    length model $ 255;
	length psiModel $ 255;
    length phiModel $ 255;
    length pModel $ 255;
    set fs&id;
	keep zone model AIC;
    drop psiModel phiModel pModel;
    if descr = "AIC (smaller is better)";
	AIC = value;
    if "&psiType" = "dot" then psiModel = "psi(.)";
	if "&psiType" = "Year" then psiModel = "psi(Year)";
	if "&psiType" = "logit" then psiModel = "psi(logit)";
	if "&phiType" = "dot" then phiModel = ",phi(.)";
	if "&phiType" = "Year" then phiModel = ",phi(Year)";
	if "&phiType" = "0" then phiModel = ",phi=0";
	if "&pType" = "dot" then pModel = ",p(.)";
	if "&pType" = "Year" then pModel = ",p(Year)";
	if "&pType" = "TwoGroup" then pModel = ",p(2 groups)";
	if "&pType" = "YearTwoGroup" then pModel = ",p(Year*2 groups)";
    model = trim(psiModel) || trim(phiModel) || trim(pModel);
  run;
%mend Analysis;