MODULE nrtype
  IMPLICIT NONE
	
  Integer , parameter      :: I4B = SELECTED_INT_KIND(9)
  Integer , parameter      :: dp  = SELECTED_REAL_KIND(15,307)
end MODULE nrtype 
!------------------
MODULE amvars
  use nrtype 
  IMPLICIT NONE


  Integer(I4B) , parameter :: NOBS = 32000  !number of individuals
  Integer(I4B) , parameter :: NPERS =    1  !number of periods/individual
  Integer(I4B) , parameter :: NRECS = NOBS*NPERS   
  
  Integer(I4B) , parameter :: MAX_MCR = 300 !# of Monte Carlo iterations
  
  Integer(I4B) , parameter :: MSTATE = 150  !number of states
  Integer(I4B) , parameter :: CAP_J  = 4    !number of choices
  Integer(I4B) , parameter :: NPARM  = 4    !number of parameters
  Real(dp)     , parameter :: BETA   = 0.90 !discount factor 

  Integer(I4B) , parameter :: NPARME = 4    ! parms to estimate Debug
  Integer(I4B) , parameter, dimension(NPARM) :: ss = (/1,1,1,1/)


  TYPE mcrt 
     Integer(I4B) :: nplit                !# PI-K iterations
     Real(dp)     :: nplfcn               !NPL loglike
     Real(dp) , dimension(NPARM) :: nplq  !NPL parm est
     Real(dp)     :: nplt                 !Seconds to compute NPL est
     Integer(I4B) :: nfxpit               !# NFXP iterations
     Real(dp)     :: nfxpfcn              !NFXP loglike
     Real(dp) , dimension(NPARM) :: nfxpq !NFXP parm est
     Real(dp)     :: nfxpt                !Secnds to compute NFXP est
  END TYPE mcrt

  TYPE(mcrt)   :: mcr

  Real(dp)     , parameter :: TOL1   = 1.0e-6 !used in amoeba
  Real(dp)     , parameter :: TOL2   = 5.0e-6 !used to define NPL Fixed Point
  Real(dp)     , parameter :: EVFTOL = 1.0e-8 !used to define EVF fixed point
  Integer(I4B) , parameter :: MAXITER = 900   !maximum number of iters (EVF)

  Real(dp)     , parameter :: TINYP   = 1.0e-20  !trap low probs
  Real(dp)     , parameter :: SUPDP   = 100.0_dp 

  Real(dp) , parameter                      ::  econst = 0.0_dp
  Real(dp) , dimension(MSTATE,CAP_J)        ::  u    !formed in Util
  Real(dp) , dimension(MSTATE,CAP_J)	    ::  v    !formed in EVF 
  Real(dp) , dimension(MSTATE)              ::  EV   !formed in EVF
  Real(dp) , dimension(MSTATE,MSTATE,CAP_J) ::  ff   !formed in Starts
  Real(dp) , dimension(MSTATE,NPARM,CAP_J)  ::  HH   !formed in Starts
  
  Real(dp) , dimension(:), allocatable		:: f0    !used in NEWTON     

!	amin, and amax define the choices available by state
  Integer(I4B) , dimension(MSTATE)  :: amin, amax 
  Integer(I4B) , dimension(NRECS,5) :: dat !data vector for debugging
end MODULE amvars
!-----------------------------------
MODULE nrutil 
use nrtype

Interface imaxloc
   Module procedure imaxloc_r, imaxloc_i
end interface

Interface outerprod
   Module procedure outerprod_d
end interface

Interface swap
   Module procedure swap_i, swap_r, swap_rv, &
        masked_swap_rs, masked_swap_rv,masked_swap_rm
end interface
	
Interface assert_eq
   Module procedure assert_eq2, assert_eq3, assert_eq4, assert_eqn
end interface

CONTAINS
  Function assert_eq2(n1,n2,string)
    !Report and die if integers not all equal (used for size checking)
    Character(LEN=*), intent(in) :: string
    Integer , intent(in) :: n1, n2
    Integer :: assert_eq2
    
    if (n1==n2) then 
       assert_eq2 = n1
    else
       write(*,*) 'nerror: an assert_eq failed with this tag: ',string
       STOP 'program terminated by assert_eq2'
    end if
  end Function assert_eq2
  
  Function assert_eq3(n1,n2,n3,string)
    Character(LEN=*), intent(in) :: string
    Integer , intent(in) :: n1, n2, n3
    Integer :: assert_eq3
    
    if (n1 == n2 .AND. n2 == n3) then
       assert_eq3 = n1
    else
       write(*,*) 'nerror: an assert_eq failed with this tag: ',string
       STOP 'program terminated by assert_eq3'
    end if
  end Function assert_eq3
  
  Function assert_eq4(n1,n2,n3,n4,string)	
    Character(LEN=*), intent(in) :: string
    Integer , intent(in) :: n1, n2, n3, n4
    Integer :: assert_eq4
    
    if (n1 == n2 .AND. n2 == n3 .AND. n3 == n4) then
       assert_eq4 = n1
    else
       write(*,*) 'nerror: an assert_eq failed with this tag: ',string
       STOP 'program terminated by assert_eq4'
    end if
  end Function assert_eq4
  
  Function assert_eqn(nn,string)
    Character(LEN=*), intent(in) :: string
    Integer , intent(in), dimension(:) :: nn
    Integer :: assert_eqn
    
    if ( all(nn(2:) == nn(1)) ) then
       assert_eqn = nn(1)
    else
       write(*,*) 'nerror: an assert_eq failed with this tag: ',string
       STOP 'program terminated by assert_eqn'
    end if
  end Function assert_eqn
  
  Function ifirstloc(mask)
    !index of first occurrence of .true. in a logical vector
    Logical , dimension(:), intent(in) :: mask
    Integer , dimension(1) :: loc
    Integer :: ifirstloc
    loc=maxloc(merge(1,0,mask))
    ifirstloc = loc(1)
    if (.not. mask(ifirstloc)) ifirstloc = size(mask)+1
  end Function ifirstloc

  Function imaxloc_r(arr)
    !index of maxloc on an array
    Real(dp), dimension(:), intent(in) :: arr
    Integer :: imaxloc_r
    Integer, dimension(1) :: imax
    imax=maxloc(arr(:))
    imaxloc_r = imax(1)
  end Function imaxloc_r
  
  Function imaxloc_i(iarr)
    Integer , dimension(:), intent(in) :: iarr
    Integer :: imaxloc_i
    Integer , dimension(1) :: imax
    imax = maxloc(iarr(:))
    imaxloc_i = imax(1)
  end Function imaxloc_i
  
  Function iminloc(arr)
    !index of minloc of an array
    Real(dp), dimension(:), intent(in) :: arr
    Integer , dimension(1) :: imin
    Integer :: iminloc
    imin = minloc(arr(:))
    iminloc = imin(1)
  end Function iminloc
  
  Subroutine nerror(string)
    Character(LEN=*), intent(in) :: string
    
    write(*,*) 'nerror:  ',string
    STOP 'Program terminated by nerror.'
  end SUBROUTINE nerror
  
  Function outerprod_d(a,b)
    Real(dp), dimension(:), intent(in) :: a,b
    Real(dp), dimension(size(a),size(b)) :: outerprod_d
    outerprod_d = spread(a,dim=2,ncopies=size(b)) * &
         spread(b,dim=1,ncopies=size(a))
  end Function outerprod_d
  
  Subroutine swap_i(a,b)
    Integer , intent(inout) :: a,b
    Integer  :: dum
    dum=a
    a=b
    b=dum
  end SUBROUTINE swap_i
  
  SUBROUTINE swap_r(a,b)
    Real(dp) , intent(inout) :: a,b
    Real(dp) :: dum
    dum = a
    a = b
    b = dum
  end SUBROUTINE swap_r
  
  SUBROUTINE swap_rv(a,b)
    Real(dp), dimension(:), intent(inout) :: a,b
    Real(dp), dimension(Size(a)) :: dum
    dum=a
    a = b
    b = dum
  end SUBROUTINE swap_rv
  
  SUBROUTINE masked_swap_rs(a,b,mask)
    Real(dp), intent(inout) :: a,b
    Logical, intent(in)     :: mask
    Real(dp) :: swp
    if (mask) then
       swp = a
       a = b
       b = swp
    end if
  end SUBROUTINE masked_swap_rs
  
  SUBROUTINE masked_swap_rv(a,b,mask)
    Real(dp), dimension(:), intent(inout) :: a,b
    Logical,  dimension(:), intent(in)    :: mask
    Real(dp) , dimension(Size(a))         :: swp
    where (mask)
       swp = a
       a = b
       b = swp
    end where
  end SUBROUTINE masked_swap_rv
  
  SUBROUTINE masked_swap_rm(a,b,mask)
    Real(dp), dimension(:,:), intent(inout)  :: a,b
    Logical , dimension(:,:), intent(in)     :: mask
    Real(dp), dimension(Size(a,1),Size(a,2)) :: swp
    where (mask)
       swp = a
       a = b
       b = swp
    end where
  end SUBROUTINE masked_swap_rm
end MODULE nrutil
!------------------------------------------------------------
!------------------------------------------------------------
MODULE nr
  use nrtype
  IMPLICIT NONE
  
CONTAINS
  
  SUBROUTINE amoeba(p,y,ftol,func,iter)
    use nrtype; use nrutil, ONLY : assert_eq, imaxloc, iminloc, nerror, swap
    IMPLICIT NONE
    Integer , intent(out) :: iter
    Real(dp), intent(in)  :: ftol
    Real(dp), dimension(:), intent(inout) :: y
    Real(dp), dimension(:,:), intent(inout) :: p
    
    Interface
       Function func(x)
         use nrtype
         IMPLICIT NONE
         Real(dp), dimension(:), intent(in) :: x
         Real(dp) :: func
       end Function func
    end Interface
    
    Integer, parameter :: ITMAX = 5000
    !	Mininization of the function, func, in N dimensions by the downhill 
    !	simplex method of Nelder and Mead.  The (N+1) x N maxtrix p is input. Its
    !	N+1 rows are N-dimensional vectors that are the vertices of the starting
    !	simplex.  Also input is the vector y of the length N+1, whose components
    !	must be preinitialized to the values of func evaluated to at the
    !	N+1 vertices (rows) of p; and ftol the fractional convergence tolerance
    !	to be achieved in the FUNCTION VALUE.  On output, p and y will have
    !	been reset to N+1 new points all within ftol of a minimum function value,
    !	and iter gives the number of function evaluations taken.
    !	Parameter: The maximum allowed number of function evaluations.
    
    Integer :: ihi, ndim								!Global variables
    Real(dp), dimension(size(p,2))  :: psum
    CALL amoeba_private
    
  Contains
    
    Subroutine amoeba_private
      IMPLICIT NONE
      Integer  :: i, ilo, inhi
      Real(dp) :: rtol, ysave, ytry, ytmp
      
      ndim=assert_eq(size(p,2),size(p,1)-1,size(y)-1,'amoeba')
      iter = 0
      psum(:) = sum(p(:,:),dim=1)
      do
         ilo    = iminloc(y(:))  !determine which point is highest (worst)
         ihi    = imaxloc(y(:))  !next highest and lowest (best)
         ytmp   = y(ihi)
         y(ihi) = y(ilo)
         inhi   = imaxloc(y(:))
         y(ihi) = ytmp
         !		print *,"Iteration:   ",iter
         !		print '(A,E24.12)',"Minimum [y(ilo)] ",y(ilo)
         !		print '(A,E24.12)',"Maximum [y(ihi)] ",y(ihi)
         !		print '(A)',"   Parameters at low and high points."
         !		print '(3(3F24.12,/))',p(ilo,:)
         !		print *," "
         !     print '(3(3F24.12,/))',p(ihi,:)
!     print *," "
         rtol   = 2.0_dp*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo)))
         !compute fraction range from highest to lowest and return if satisfactory
         if (rtol < ftol) then	!if returning put best point and value in slot 1
            call swap(y(1),y(ilo))
            call swap(p(1,:),p(ilo,:))
            RETURN
         end if
         if (iter >= ITMAX) call nerror('ITMAX exceeded in amoeba')
         !begin new iteration.  First extrapolate by a factor -1 through the
         !face of the simplex across from the highest point, i.e., reflect
         !the simplex from the high point.
         ytry=amotry(-1.0_dp)
         if (ytry <= y(ilo)) then	!gives result better than the best point
            ytry = amotry(2.0_dp)	!so try an additional extrapolation by
            iter = iter + 1			!a factor of 2
         else if (ytry >= y(inhi)) then !the reflected point is worse than the
            ysave=y(ihi)					 !2nd highest, so look for an intermediate
            ytry = amotry(0.5_dp)		 !lower pt, i.e., do a one-dimensional 
            iter = iter+1					 !contraction
            
            if (ytry >= ysave) then	!cant seem to get rid of that high point. 
               !Better contract around the lowest (best)
               !point.
               p(:,:) = 0.5_dp*(p(:,:)+spread(p(ilo,:),1,size(p,1)))
               do i=1,ndim+1
                  if (i /= ilo) y(i) = func(p(i,:))
               end do
               iter = iter+ndim	!keep track of function evaluations
               psum(:) = sum(p(:,:),dim=1)
            end if
         end if
      end do     !go back for the test of doneness and the next iteration
    end SUBROUTINE amoeba_private
    
    Function amotry(fac)
      IMPLICIT NONE
      Real(dp), intent(in) :: fac
      Real(dp) :: amotry
      !Extrapolates by a factor, fac, through the face of the simplex across from
      !the highest point, tries it, and replaces the high point if the new point
      !is better.
      Real(dp) :: fac1, fac2, ytry
      Real(dp) , dimension(size(p,2)) :: ptry
      
      fac1 = (1.0_dp-fac)/ndim
      fac2 = fac1-fac
      ptry(:) = psum(:)*fac1-p(ihi,:)*fac2
      ytry=func(ptry)		      !evaluate function at trial point
      if (ytry < y(ihi)) then    !If it's better than the highest, then replace
         y(ihi) = ytry				!the highest.
         psum(:) = psum(:) - p(ihi,:) + ptry(:)
         p(ihi,:)= ptry(:)
      end if
      amotry = ytry
    end Function amotry
  end SUBROUTINE amoeba
!--------------------------------------------------------------
SUBROUTINE ludcmp(a,indx,d)
  USE nrtype; USE nrutil, ONLY : assert_eq,imaxloc,nerror,outerprod,swap
  IMPLICIT NONE
  REAL(dp), DIMENSION(:,:), INTENT(INOUT) :: a
  INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: indx
  REAL(dp), INTENT(OUT) :: d
  REAL(dp), DIMENSION(size(a,1)) :: vv
  REAL(dp), PARAMETER :: TINY=1.0e-20_dp
  INTEGER(I4B) :: j,n,imax
  
  n=assert_eq(size(a,1),size(a,2),size(indx),'ludcmp')
  !	print *,"ludcmp:: size of matrix = ",n
  d=1.0
  vv=maxval(abs(a),dim=2)
  if (any(vv == 0.0)) call nerror('singular matrix in ludcmp')
  vv=1.0_dp/vv
  do j=1,n
     imax=(j-1)+imaxloc(vv(j:n)*abs(a(j:n,j)))
     if (j /= imax) then
        call swap(a(imax,:),a(j,:))
        d=-d
        vv(imax)=vv(j)
     end if
     indx(j)=imax
     if (a(j,j) == 0.0) a(j,j)=TINY
     a(j+1:n,j)=a(j+1:n,j)/a(j,j)
     a(j+1:n,j+1:n)=a(j+1:n,j+1:n)-outerprod(a(j+1:n,j),a(j,j+1:n))
  end do
END SUBROUTINE ludcmp
!---------------------------------------------------------
SUBROUTINE lubksb(a,indx,b)
  USE nrtype; USE nrutil, ONLY : assert_eq
  IMPLICIT NONE
  REAL(dp), DIMENSION(:,:), INTENT(IN) :: a
  INTEGER(I4B), DIMENSION(:), INTENT(IN) :: indx
  REAL(dp), DIMENSION(:), INTENT(INOUT) :: b
  INTEGER(I4B) :: i,n,ii,ll
  REAL(dp) :: summ
  n=assert_eq(size(a,1),size(a,2),size(indx),'lubksb')
  ii=0
  do i=1,n
     ll=indx(i)
     summ=b(ll)
     b(ll)=b(i)
     if (ii /= 0) then
        summ=summ-dot_product(a(i,ii:i-1),b(ii:i-1))
     else if (summ /= 0.0) then
        ii=i
     end if
     b(i)=summ
  end do
  do i=n,1,-1
     b(i) = (b(i)-dot_product(a(i,i+1:n),b(i+1:n)))/a(i,i)
  end do
END SUBROUTINE lubksb


SUBROUTINE GIH(qe,Grd,VCM,fval,Func)
  use amvars
  IMPLICIT NONE
  
  Interface
     Function Func(x)
       use amvars
       IMPLICIT NONE
       Real(dp), intent(in), dimension(:) :: x
       Real(dp) :: Func
     end Function Func
  end Interface
  
  
  Real(dp), intent(in), dimension(:) :: qe !estimated parameters 
  Real(dp), intent(out), dimension(:)   :: Grd  !Gradient
  Real(dp), intent(out), dimension(:,:) :: VCM !Inverse Hessian
  Real(dp), intent(inout)               :: fval
  
  ! Calculates numerical approximation of gradient g and outer product
  ! approximation to estimate Hessian
  ! Returns Gradient Grd, Inverse Hessian, and function value
  ! Dennis and Schnabel (1996) recommend setting the stepsize, h, for the
! function evaluation as sqrt(macheps)*typx(x), where typx(x) is the "typical"
! value of the parameter x.  macheps is the smallest number such that 
! 1+macheps > 1. This formula assumes that the function has optimized to 
! precision macheps. Taking the sqrt of macheps gives precision of the left-
! half of the digits.
! The instrinsic function EPSILON() returns macheps, about 1.e-16 and yields
! very small (probably too small) step sizes.

! 7/22/01 include trap for coefficients at zero. Assumes scaling of variables
! so estimated coefficients will be within the unit interval (takes midpt).
!
  Integer(I4B) 		:: i, k
  Real(dp)    :: h, tsize
  Real(dp) , parameter :: macheps = EPSILON(qe(1)) !est parameters same type
  Real(dp) , parameter :: TINY = 1.0e-20
  Real(dp) , allocatable, dimension(:,:) :: qdir, G, J, HESS
  Real(dp) , allocatable, dimension(:)   :: pt, b
  Real(dp) , dimension(0:NPARME)         :: y      !function values
  Integer(I4B) , dimension(:), allocatable :: indx
  
  if ( (size(qe) .NE. size(Grd)) .OR. (size(qe) .NE. size(VCM,dim=1)) ) then
     print *,"Error:  Matrices in GIH do not conform."
     print *,"   Should equal the value of NPARME, which is ",NPARME
     print *,"   size of qe=",size(qe)
     print *,"   size of Grd=",size(Grd)
     print *,"   size of VDM=",size(VCM,dim=1)
     print *,"Kill job."
     stop
  end if
  
  allocate(qdir(NPARME,NPARME))  !needed to build basis for amoeba
  allocate(G(NRECS,0:NPARME))    !hold function evaluations
  allocate(J(NRECS,NPARME))      !matrix NRECS x NPARME of ind obs gradient
  if (.not. allocated(f0)) allocate(f0(NRECS))
  allocate(pt(1:NPARME)) 
  
  !	set up direction matrix 
  qdir = 0.0_dp
  forall(i=1:NPARME) qdir(i,i) = 1.0_dp
  
  !	calculate baseline loglikelihood value
  !	print *,"Debug (GIH):: In GIH."
  if (fval == 0.0) then  !do only when calculating standard errors
     f0(:) = 0.0_dp !f0 is a global variable contain ind contr to loglike
     !		print *,"Calculate baseline function evaluation: "
     !		print *,"  parameters (qe): ",qe
     !		print *,"Now invoke Func(qe)"
     y(0) = -Func(qe)  !note y(0) is loglike
     G(:,0) = f0(:) !individual observations ll values
     !		print *,"Debug:: GIH: baseline loglikelihood value=",y(0)
  else 
     y(0) = fval    !must be loglikelihood (not -loglikelihood value)
     G(:,0) = f0(:) !available from last call to Func during optimization
  end if
  
  !	print *," " 
  !	print *,"Debug:: GIH:  macheps=",macheps
  !	print *," " 
  do i =1,NPARME
     if (abs(qe(i)) < TINY) then !trap coef = 0
        tsize = 0.5_dp !assumed scaling to within unit interval
     else 
        tsize = abs(qe(i))
     end if
     h = SQRT(macheps)*tsize  !stepsize
     pt(:) = qe + qdir(i,:)*h
     !		print *," " 
     !		print *," parameter: ",i," pt of evaluation: "
     !		print '(3F12.8)',pt(:)
     !		print *," stepsize: ",h
     y(i) = -Func( pt(:) )
     !		print *," (GIH) loglikelihood value= ",y(i)
     G(:,i) = f0(:) !f0 evaluated in Func
     J(:,i) =	(G(:,i) - G(:,0)) / h 
     Grd(i) = sum(J(:,i))
     !		print *,"Parameter            Sum Gradient                Step Size"
     !		print '(I8,1X,E24.12,1X,E24.12)', i,Grd(i),h
  end do !i
  
  print *," "
  !	calculate outer product approximation of Hessian
  deallocate(G)
  allocate(HESS(NPARME,NPARME))
  !	allocate(VCM(NPARME,NPARME))
  allocate(indx(NPARME))
  allocate(b(NPARME))
  HESS = matmul(transpose(J),J)
  Call ludcmp(HESS,indx,h) !on output HESS contains LU decomposition
  do k=1,NPARME
     b = 0.0_dp
     b(k) = 1.0_dp
     Call lubksb(HESS,indx,b)
     VCM(:,k) = b !build inverse column by column
  end do
  !	print *," "
  !	print *,"Inverse Hessian (Variance-Covariance Matrix): "
  !	print '(5(E14.8,1X))',VCM
  forall(k=1:NPARME) b(k) = VCM(k,k)
  !	print *," "
  !	print *,"Standard Errors: "
  !	print '(5(E14.8,1X))',sqrt(b)
end Subroutine GIH
!--------------------------------------------------------------
SUBROUTINE Newton(qe,fval,Func)
  use amvars 
  IMPLICIT NONE

  Real(dp) , intent(inout) , dimension(:)  :: qe   !estimated parameters
  Real(dp) , intent(inout)                 :: fval !new function value

  Interface
     Function Func(qe)
       use amvars
       IMPLICIT NONE
       Real(dp) , intent(in), dimension(:) :: qe
       Real(dp) :: Func
     end FUNCTION Func
  end Interface

  Real(dp)                                 :: oldfval !old function value

  Integer :: k
  Real(dp) , dimension(:), allocatable :: dqe, qe1, ser, Grd
  Real(dp) , dimension(:,:), allocatable :: VCM 

!	print *,"In Newton..."
!	print *,"Size of qe=",size(qe)
  allocate(dqe(Size(qe)))
  allocate(qe1(Size(qe)))
  allocate(ser(Size(qe)))
  allocate(Grd(Size(qe)))
  allocate(VCM(Size(qe),Size(qe)))
  
  Grd(:) = 0.0_dp
  VCM(:,:) = 0.0_dp

!	print *," "
!	print *,"Calculate gradient, inverse Hessian (VCM) and next value of loglike."
  Call GIH(qe,Grd,VCM,fval,Func)  !calculate the Gradient and Inverse Hessian
!	print *," "
!	print *," "
  print *,"Parm       Estimate              SE            Grad            Step"
  print *," ------------------------------------------------------------------"
  dqe(:)  = -matmul(VCM,Grd)
  forall(k=1:NPARME) ser(k) = sqrt(VCM(k,k))
  do k=1,NPARME
     print '(I4,4(1X,E15.8))',k,qe(k),ser(k),Grd(k),dqe(k)
  end do
  print *," ------------------------------------------------------------------"
  !	print *, "maximum step ", maxval(abs(dqe))

  qe = qe - dqe 

!	fval should be bigger than oldfval (both are negative numbers)
!	oldfval = fval
!	do k = 1,15
!		print *," "
!		print *,"updated parameters: "
!		print '(3E24.12,/,3E24.12,/,3E24.12)',qe
!		fval   = -Func( qe ) !new function value
!		if (fval > 1.0*oldfval) then
!			exit ! do loop
!		else
!			qe = qe + (0.5**k)*dqe
!			print *, "notweN", k
!			print *,"Function change: from ", oldfval, " to ",   fval
!			if (k == 15) then
!				stop
!			end if
!		end if
!	end do
!
!	print *," "
!	print *,"Function evaluation (at next step): ",fval
!	print *," "
  deallocate(Grd)
  deallocate(VCM)
end SUBROUTINE Newton
!-------------------------------------------------
end MODULE nr
!-------------------------------------------------
!-------------------------------------------------
MODULE NFXPP  
  use amvars
  use nr
  IMPLICIT NONE

  Real(dp), dimension(NPARM) :: q  !fixed and est parms NFXP

CONTAINS
!------------------
  FUNCTION Func(qe)  !For NFXP
    use nrtype; use amvars
    IMPLICIT NONE

    interface 
       Subroutine UTIL(x,pflag)
         use nrtype; use amvars
         IMPLICIT NONE
         Real(dp), intent(in), dimension(NPARM) :: x
         Integer(I4B) , intent(in)              :: pflag
       end subroutine UTIL
end interface

interface
   Subroutine EVF(x)
     use amvars
     IMPLICIT NONE
     Real(dp), intent(in), dimension(NPARM) :: x
   end Subroutine EVF
end interface

Real(dp) , intent(in), dimension(:) :: qe

Real(dp) , parameter :: TINY = 1.0e-20

Real(dp) :: Func

Integer(I4B) ::  i,j,nn,cx,a,nx
Real(dp) :: ll, tval

!  update (scope within NFXP) q parameter vector for estimated elements

j=1
do i=1,NPARM
   if (ss(i) == 1) then
      q(i) = qe(j)         !q global
      if (j==NPARME) exit  !skip out of loop
      j = j+1
   end if
end do

!	if (.not. allocated(f0)) allocate(f0(NRECS))
!	f0(:) = 0.0_dp !for use in Newton

!	print *," In (NFXP) Func: ......"

!	Calculate utility at new parameter values q 
!	print '(A,/,3(3E24.12,/))'," CALL UTIL q=",q

CALL UTIL(q,0) !calculates utility for this set of home locations

!	print *,"Utility done."
!	print *," CALL EVF"

CALL EVF(q)

!	print *,"Debug(NFXP):: EVF done."
!	print *,"Debug(NFXP):: EVF finished, now calculate likelihood function."

ll = 0.0_dp

!	print *," Calculate loglikelihood."
!	print *," Number of observations: ",NRECS
!	print *," # of Rows of dat(:,:) =",SIZE(dat,1)
!	print *," # of Cols of dat(:,:) =",SIZE(dat,2)
do nn =1,NRECS 
   cx = dat(nn,3) !current state
   a  = dat(nn,4) !decision
   nx = dat(nn,5) !next period's state                
   
!	calculate weighted value of EVF function for revised conditional
!	choice probability; assume populations scaled to unit interval.  
!	See EVF Function.

!		if (nn <= 5) print '(A,4I5)',"nn,cx,a,nx ",nn,cx,a,nx 

   tval  =  v(cx,a) - EV(cx) 
!		f0(nn) = tval 
   ll = ll + tval  
!		if (nn < 10) then
!				print '(A,I3,F12.6,2I4,2F12.6)',"nn,tval,cw,a,v, EV= ",&
!					nn,tval,cx,a,v(cx,a),EV(cx)
!		end if
end do !nn		
!	print *,"   (func)  function value is (-ll): ",-ll

func = -ll    !returns -loglike  to min{-loglike}
end FUNCTION func  !NFXP
!-----
end MODULE NFXPP
!------------------------------------------------------------
MODULE NPLP
  use nrtype
  use amvars 
  IMPLICIT NONE
!define functions and subroutines used by NPL algorithm
!variables defined with NPLP scope

  Real(dp) , dimension(MSTATE,NPARM,CAP_J) :: ht
  Real(dp) , dimension(MSTATE,CAP_J)       :: et
  Real(dp) , dimension(MSTATE,CAP_J)       :: lp
  Real(dp) , dimension(MSTATE,MSTATE)		  :: IFUPI
  Real(dp) , dimension(NPARM)              :: q

CONTAINS
  Subroutine EP(x,CCP)
    use amvars
    IMPLICIT NONE
!calculates logit probabilities for a given parameter vector qe, and 
!CCP as embedded in ht and et arrays. 
!	note ht and et are variables within NPL and depend on Choice probs 

    Real(dp) , dimension(NPARM), intent(in) :: x      !vector of parms
    Real(dp) , dimension(MSTATE,CAP_J), intent(out) :: CCP 

    Integer :: i, a

    Real(dp) , dimension(MSTATE,CAP_J) :: numer, tnum
    Real(dp) , dimension(NPARM)  :: z
    Real(dp)  :: nscale, denom

!	print *,"In EP."
!	print *,"Size of x: ",size(x)
!	print '(A,4F10.5)',"x (passed): ",x
!	print '(A,4F10.5)',"q (global): ",q

    numer(:,:) = 0.0
    CCP(:,:)   = 0.0
    z=x
    if (abs(x(3)) .GT. SUPDP) then
       print *,"(EP) trapped large x(3)=",x(3)
       print *,"z(3) set to TINYP."
       z(3)=TINYP  
    else 
       z(3) = 1.0_dp/(1.0_dp+exp(-x(3)))
    end if
!	print *,"Account for logit structure of x(3)."
!	print '(A,4F10.5)',"z: ",z
!	print *,"Calculation of PHI(P) using ht and et"
    do i =1,MSTATE
!		print *,"i=",i
       do a = amin(i),amax(i) 
!			print *,"ht(i,:,",a,"): "
!			print '(4F10.5)',ht(i,:,a)
!			print '(A,F10.5)',"et(i,a)=",et(i,a)
          tnum(i,a) = dot_product(ht(i,:,a),z) + et(i,a)
       end do !a
    end do !i
    nscale = maxval(tnum)
    do i = 1,MSTATE
       denom = 0.0_dp
       do a=amin(i),amax(i)
          denom = denom + exp((tnum(i,a) - nscale))
       end do
       do a=amin(i),amax(i)
          CCP(i,a) = exp(tnum(i,a) - nscale)/denom
       end do !a
    end do !i
    
!	print *,"In subroutine EP."
!	print *,"CCP as calculated from ht and et"
!	print '(4I10)',(i,i=1,CAP_J)
!	do i=1,MSTATE
!		print '(4F10.5)',CCP(i,:)
!	end do !a
end Subroutine EP
!--------------------------------------
Function Func(qe)   !NPL 
  use amvars 
  IMPLICIT NONE
  Real(dp) , intent(in), dimension(:)   :: qe
  Real(dp) :: Func

  Real(dp) , dimension(MSTATE,CAP_J) :: CCP   !working version

  Real(dp) :: ll
  Integer(I4B) :: n, xi, ai, i, j 

  CCP(:,:) = 0.0_dp  !initialize
	
! convert estimates parameters in full parameter vector
!	print *,"Debug(Func-NPL):: parameters passed qe."
!	print *,"Size of qe=",size(qe)
!	print '(3F12.6)',qe(:)
!	print *," "
!	print *,"Global vars q: "
!	print '(4F12.6)',q(:)
  j=1
  do i=1,NPARM
     if (ss(i) == 1) then
        q(i) = qe(j)            !q local
        if (j==NPARME) exit
        j = j+1
     end if
  end do
 
!	print *,"Updated q:"
!	print '(4F12.6)',q(:)

  CALL EP(q,CCP)   !returns CCP for q parms
	
  ll = 0.0_dp
  f0(:) = 0.0_dp
  do n=1,NRECS
     xi = dat(n,3)
     ai = dat(n,4)
!		if (n <= 5) then
!			print '(A,3I5,F10.5)',"in func: n",n,xi,ai,CCP(xi,ai)
!		end if
     if (CCP(xi,ai) < TINYP) then
        print *,"(NPL(func)  CCP =",CCP(xi,ai)," which is less than TINYP, reset to TINYP."
        CCP(xi,ai) = TINYP
     end if
     f0(n) = log(CCP(xi,ai))
     ll = ll + f0(n)
  end do
  Func = -ll       !returns -loglike to min(-loglike) 
end FUNCTION Func
!---------------------------------
Function Had1(x,y)
  use nrtype
  IMPLICIT NONE
  Real(dp) , dimension(:), intent(in) :: x,y
  Real(dp) , dimension(size(x,1)) :: Had1, g
  Integer :: l
  
  do l=1,size(x,1)
     g(l) = x(l)*y(l)
  end do
  Had1 = g
end Function Had1
!---------------------------
Function Had2(x,y)
  use nrtype
  IMPLICIT NONE
  Real(dp) , dimension(:),   intent(in) :: x 
  Real(dp) , dimension(:,:), intent(in) :: y
  Real(dp) , dimension(size(x,1),size(y,2)) :: Had2,g
  Integer :: l
  
!	print *,"In Function Hadamard."
!	print *,"Number of rows (from x): ",size(x,1)
!	print *,"Number of columns (from y): ",size(y,2)
  do l=1,size(x,1)
     g(l,:) = x(l) * y(l,:)
  end do
  Had2 = g
!	print *,"In Hadamard...direct product is"
!	print *,G  
end Function Had2
!---------------------------
Subroutine VF(q,CCP0,CCP1)
  use amvars
  IMPLICIT NONE

  Real(dp), intent(in),  dimension(NPARM)      :: q  !util fcn parms
  Real(dp), intent(in),  dimension(MSTATE,CAP_J)  :: CCP0 !initial CCP
  Real(dp), intent(out), dimension(MSTATE,CAP_J)  :: CCP1 !updated

Interface
   Subroutine UTIL(q,pflag)
     use amvars
     IMPLICIT NONE
     Real(dp), intent(in), dimension(NPARM) :: q
     Integer(I4B), intent(in)               :: pflag
   end Subroutine Util
end Interface


Real(dp), dimension(MSTATE)			:: VP, tval, denom
Real(dp) :: maxv
Integer(I4B) :: a,i


!	print *,"Debug in VF, parameters q: "
!	print '(4F12.6)',q
!	print *," "

!	form utility function (x,a)
CALL UTIL(q,0)

!	print *,"Utility flows u(x,a) in VF" 
!	print '(4I10)',(a,a=1,CAP_J)
!	do i=1,MSTATE
!		print '(4F10.5)',u(i,:)
!	end do

!	print *,"Debug(VF):: Call EVF."
!	CALL EVF(q) 
!	print *,"Choice Specific Utility Flows from EVF"
!	print '(4I10)',(i,i=1,CAP_J)
!	do i=1,MSTATE
!		print '(4F10.5)',v(i,:)
!	end do

!	print *,"Debug:: CCP as calculated from EVF."
!	CALL Eprob(CCP1)  !uses v(:,:)

!	print *,"lp values"
!	do i=1,MSTATE
!		print '(4F10.5)',lp(i,:)
!	end do
!	print *," "
!	print *,"u(:,a) + lp(:,a)"
!	do i=1,MSTATE
!		print '(4F10.5)',u(i,:)+lp(i,:)
!	end do
	tval = 0.0_dp  !MSTATEx1 vector
!	print *,"tval by a (maybe) should be 5 x 1 vector"
 do a=1,CAP_J
!		print '(A,I3,5F10.5)',"a= ",a,Had1(CCP0(:,a),u(:,a)+lp(:,a))
    tval = tval + Had1(CCP0(:,a),u(:,a)+lp(:,a))
end do
		
VP(:) = matmul(IFUPI,tval)

!	print *,"In VF"
!	print *," "
!	print *,"Value function for parameters q, and CCP0"
!	print '(5(F10.5,2X))',VP
!	print *," "


! calculate choice specific value value v
v(:,:) = 0.0_dp
do a=1,CAP_J
   v(:,a) = u(:,a) + BETA*matmul(ff(:,:,a),VP)
end do !a

!	print *,"Choice Specific utility, v"
!	print '(4I10)',(i,i=1,CAP_J)
!	do i=1,MSTATE
!		print '(4F10.5)',v(i,:)
!	end do

!	print *,"Debug:: CCP as calculated from VF -- v(cx,a)."
CALL Eprob(CCP1)  !uses v(:,:)
!	print *," "

end Subroutine VF
!-----
end MODULE NPLP
!============================================================
program am7
  use amvars
!Version 7b: increase MSTATE to 100 change gam probabilities
!Version 7a: turn on NPL estimator 
!Version 7 : production version for monte carlo runs
!   copied from Version 6. July 7, 2003
!   writes out to file run,iter-NPL,q(NPL),iter-NFXP,q(NFXP)
!
!Version 6: recode utility function AND gets both NPL and NFXP to run 
!	 u(x,a) = 0.0 if a=1 (no sexual activity)
!   u(x,a) = q(1)*log(sa(cx)) - 1[a=2]*q(2) - 1[a=3]*log(sa(cx)), a=2,3
!	 u(x,a) = q(4)
!   sa(:) = (0,5,15,25,preg) assumed level of sexual activity to set utility
!   values 
!Date: July 6, 2003
!Copied from Version5 
!

! Program codes Nested Pseudo Likelihood Estimator of Aguirregabiria and
! Mira (Econometrica, 2002) and the NFXP Estimator of Rust for a simple
! model of pregnancy and contraception.
!
!
!
  Interface
     Subroutine Eprob(CCP) !calculates CCP from choice-specific util v(x,a)
       use amvars
       IMPLICIT NONE
       Real(dp), intent(out), dimension(MSTATE,CAP_J) :: CCP
     end Subroutine Eprob
  end Interface
  
  Interface 
     Subroutine gendat
       use amvars
       IMPLICIT NONE
     end Subroutine gendat
  end Interface
  
  Interface
     Subroutine NFXP(x)
       use amvars; use NFXPP
       use nr
       IMPLICIT NONE
       Real(dp) , intent(inout), dimension(NPARM) :: x
     end Subroutine NFXP
  end interface
  
  Interface
     Subroutine NPL(KITER,CCP,x)
       use amvars; use nr
!       use IMSLF90 
       Integer(I4B), intent(in) :: KITER
       Real(dp), intent(in), dimension(MSTATE,CAP_J) :: CCP
       Real(dp), intent(inout), dimension(NPARM)       :: x
     end Subroutine NPL
  end interface
  
  Interface 
	Subroutine sumstat(CCP) !initial estimate of CCP   
   use amvars
   Real(dp), intent(out), dimension(MSTATE,CAP_J) :: CCP
 end Subroutine sumstat
end Interface

! ================================================================
! ================================================================	

! main program

Real(dp),parameter,dimension(NPARM) :: qstar= (/1.,1.,2.1972,1.0/)
Real(dp) , dimension(NPARM)				:: q !parameters
Real(dp) , dimension(MSTATE,CAP_J) 		:: CCP

Real(dp)     :: t1,t2     !timing variables
Integer(I4B) :: i,KK, r

!		Decision variables:
!	a=1 -> not sexually active
!	a=2 -> sexually active, high level of contraception
!	a=3 -> sexually active, low level of contraception
!	a=4 -> dummy variable for absorbing state
!		State variables
!	cx=1 -> not sexually active
!	cx=2 -> low level of sexual activity
!	cx=3 -> middle level of sexual activity
!	cx=4 -> high level of sexual activity
!	cx=5 -> pregnant (absorbing state)


!	print *,"initial random number seed."

CALL RANDOM_SEED()  !initialize random number generater

!	print *,"Call Starts."
print *," "

!	open(unit=8,FILE="am7b3.dat")

!  Define Choice SET
amin(1) = 1;  amin(2:MSTATE-1) = 2; amin(MSTATE) = CAP_J
amax(1:MSTATE-1) = 3;  amax(MSTATE) = CAP_J 

print *,"MSTATE  MAX_MCR      NOBS              qstar"
print '(I6,3X,I7,I10,4F8.4)',MSTATE,MAX_MCR,NOBS,qstar(:)
		
do r = 1,MAX_MCR

!	initialize mcr
   mcr%nplit   = 0
   mcr%nplfcn  = 0.0_dp
   mcr%nplq(:) = 0.0_dp
   mcr%nplt    = 0.0_dp
   
   mcr%nfxpit   = 0
   mcr%nfxpfcn  = 0.0_dp
   mcr%nfxpq(:) = 0.0_dp
   mcr%nfxpt    = 0.0_dp

   u(:,:) = 0.0_dp
   v(:,:) = 0.0_dp
   CCP(:,:) = 0.0_dp
   EV(:)  = 0.0_dp
   ff(:,:,:) = 0.0_dp
   HH(:,:,:) = 0.0_dp
   dat(:,:)  = 0       !integer
		
!		print *,"Monte Carlo Iteration = ",r
!		print '(A,4F10.5)',"qstar: ",qstar(:)
!		print *," "

   CALL STARTS  !initialize matrices Pi, gam, HH, ff

   CALL UTIL(qstar,0) !fill utility matrix u(x,a)
   CALL EVF(qstar)  !solve for fixed point for q0
   CALL Eprob(CCP)
   CALL CPU_TIME(t1) !Fortran 95 instrinsic function
   CALL gendat	     !generated simulated histories (fills dat(NRECS,5))
   CALL CPU_TIME(t2) !Fortran 95 instrinsic function
   !		print *,"Time(seconds) required to generate a dataset= ",t2-t1

!		print *,"generated data for iteration ",r
   CALL sumstat(CCP) !initializes CCP0
!		print *,"Choice Probabilities."
!		do i=1,MSTATE
!			print '(4F10.6)',CCP(i,:)
!		end do

   KK =  25
   
   q = qstar     !initialize parameters 
   do i=1,NPARM  !initialize parameters for NPL
      if (ss(i) == 1) then
         q(i) = 0.1 ! 9*qstar(i)   !sv
      end if
   end do

!		print *,"CALL NPL."
   Call CPU_TIME(t1)  !Fortran 95 instrinic function
   CALL NPL(KK,CCP,q) !assign fields of mcr within NPL
   CALL CPU_TIME(t2)  
   mcr%nplt = t2-t1
   !		print *,"Seconds required for NPL iteration= ",t2-t1
!		print *," "
   
   q = qstar  ! set global parm vector
   do i=1,NPARM
      if (ss(i) == 1) then
         q(i) = 0.9*qstar(i)   !sv
      end if
   end do

!		print *,"CALL NFXP."
   Call CPU_TIME(t1)  !Fortran 95 instrinic function
!		print '(A,4F10.5)',"parameters: ",q(:)
   CALL NFXP(q) !assign fields of mcr within NFXP
   CALL CPU_TIME(t2)  
   mcr%nfxpt = t2-t1
!		print *,"Seconds required for NFXP iteration= ",t2-t1
   write(unit=*,FMT='(I5,2(I5,F18.8,4F10.5,F12.4))') r,mcr
end do
!	close(unit=8)
end program am7
!-----------------------------------------------------------------
Subroutine Eprob(CCP) !uses v(x,a) calculated in EVF to set CCP
use amvars
Real(dp) , intent(out), dimension(MSTATE,CAP_J) :: CCP

Integer(I4B) :: xi, a
Real(dp)     :: denom,vscale

CCP = 0.0_dp
vscale = maxval(v)
do xi=1,MSTATE
   denom = 0.0_dp
   do a=amin(xi),amax(xi)
      denom = denom + exp(v(xi,a)-vscale)
   end do
   do a=amin(xi),amax(xi)
      CCP(xi,a) = exp(v(xi,a) - vscale)/denom
   end do
end do !xi
end Subroutine Eprob
!-----------------------------------------------------------------
Subroutine EVF(x)
!	Calculate the Expected value function for parameter vector x.
  use amvars
  IMPLICIT NONE
  Real(dp), intent(in), dimension(NPARM) :: x
  
  Real(dp), parameter :: TINY = 1.0e-20
  
  Real(dp), dimension(MSTATE) :: EV0 !work
  Real(dp), dimension(MSTATE) :: p		!tmp value for cond transition prob
  
  Integer(I4B) :: iter,cx,a
  Real(dp) :: vscale,tev2,tval,maxdiff
  
  !	Initialize EV0, EV

  EV = 0.0_dp  	!global
  v  = 0.0_dp		!global


!	print '(A,4F10.5)',"Debug(EVF):: parameters: ",x
	
  do iter=1,MAXITER
     EV0 = EV
     vscale = maxval(v(:,:)) !use same scale factor for all states
!		print *,"vscale=",vscale
     do cx=1,MSTATE   	!cx = current state, loop over states
        tval = 0.0_dp
        do a = amin(cx),amax(cx)  !state specific choice set 
!				print *,"Debug(EVF):: state action: ",cx,a
           p(:) = ff(cx,:,a)  	!conditional (on a) transition matix
!				print '(A,5F7.3)',"  transition prob: ",p(:)
           tev2 = vscale + dot_product(p(:),(EV0(:) - vscale)) 	
!				print *,"tev2 = <p(:),EV0-vscale> ",tev2
           v(cx,a) = u(cx,a) + BETA*tev2
!				print *,"  u =",u(cx,a),"   v(cx,a)=",v(cx,a)
           tval = tval + exp(v(cx,a)-vscale)
        end do !a
		
!			print *,"--------- cx = ",cx,"  ---------  tval = ",tval

        if (tval <= 0.0_dp) then
           print *,"ERROR:  NUMERICAL PROBLEMS in EVF."
           print *,"   tval is non-positive."
           print '(A,F10.5)',"tval= ",tval
           print '(A,/,7E20.12)'," parameter vector: ",x(:)
           print '(A,/,4F10.5)',"exp(v(cx,:) =",exp(v(cx,:))
           print *," "
           print *,"Job Stopped."
           stop
        end if
        if (tval .gt. 0.0_dp .and. tval .lt. TINYP) then
           print *,"tval < TINYP, tval= ",tval," tval reset to TINYP."
           tval = TINYP
        end if
        EV(cx) = vscale + log(tval) 
	
     end do !cx
		
     maxdiff = maxval( abs(EV - EV0) )

     if (maxdiff < EVFTOL) exit !iter do loop

  end do !iter

!	if (iter .GE. MAXITER) then
!		print *,"Expected value function not optimized."
!		print *,"Reached Maximum number of iterations."
!		print *,"Maxdiff= ",maxdiff
!		print *,"Job stopped."
!		stop
!	else
!		print *," "
!		print *,"Fix Point Found!! - Expected Value Function calculated."
!		print *,"Interation = ",iter,"  Maximum Deviation=",maxdiff
!		print *,"EVF::   EV(:)"
!		print '(5F10.5)',EV(:) 
!		print *,"v(cx,a)"
!		print *,"cx  v(cx,:)"
!		do cx=1,MSTATE
!			print '(I2,3X,4F10.5)',cx,(v(cx,a),a=1,CAP_J) 
!		end do 
!	end if !iter
end subroutine EVF
!------------------------------------------------------------
Subroutine gendat !generates data vector dat(NRECS,5)
  use amvars
!	Uses choice-specific maximum values calculated by EVF, v(x,a).  v(x,a) global var

  Real(dp) , dimension(CAP_J)   :: rv, e1  !type 1 extreme value
  Real(dp) , dimension(MSTATE)  :: p       !conditional probability vector
  
  Real(dp)     :: z, maxv, cdf
  Integer(I4B) :: n, t, nt, xi, astar, xj, j, a

!	CALL RANDOM_SEED()  !initialize random number generater
  nt = 0
  do n=1,NOBS 
     do t=1,NPERS 
        nt=nt+1
        dat(nt,1) = n
        dat(nt,2) = t
        ! generate initial state xi
        if (t==1) then
           CALL Random_number(z)       !a scalar
           xi = 1+floor(MSTATE*z)  	!initial state, generates xi
           if (xi > MSTATE) xi=MSTATE		
        else
           xi = dat(nt-1,5)
        end if
        ! generate taste shocks
        CALL Random_number(rv)     !rv: CAP_J-1 x 1 vector 
        e1(:) = -log(-log(rv(:)))	!type one Extreme Value
        ! select action a
        maxv = -9.9E+10
        astar = 0
        do a=amin(xi),amax(xi)
           if (v(xi,a) + e1(a) .GT. maxv) then
              maxv = v(xi,a)+e1(a)
              astar = a
           end if
           if (astar == 0) then
              print *,"ERROR(gendat):: Error in simulation."
              print *,"    Optimal choice not defined."
              print *,"    current state, xi: ",xi
              print *,"Job Killed."
              stop
           end if
        end do !a
        ! define transition probability
        if (astar == 1 .AND. xi > 1) then
           print *,"Debug(gendat):: ERROR:: astar == 1 and xi=",xi
           print *,"Job Killed."
           stop
        end if
        p(:) = ff(xi,:,astar)       !MSTATE x 1 vector
        if (xi == MSTATE) then
           xj = MSTATE
        else
           CALL Random_number(z)	!a scalar
           !				print '(A,5F8.5)',"p(:) ",p(:)
           cdf = 0.0_dp
           xj = 0
           do j=1,MSTATE
              cdf = cdf + p(j)
              if (z <= cdf) then
                 xj = j
                 exit !do
              end if
           end do
           if (xj==0) then
              print *,"ERROR: destination state is unassigned."
              print *,"Job Killed."
              stop
           end if
        end if
        
        dat(nt,3) = xi
        dat(nt,4) = astar
        dat(nt,5) = xj
     end do !t
  end do !n
  !	print *,"debug(gendat):: last call to random number generator: ",z
end SUBROUTINE gendat
!------------------------------------------------------------
!------------------------------------------------------------
Subroutine NFXP(x)
  use amvars; use NFXPP !procedures for NFXPP
  use nr
  IMPLICIT NONE
  
  Real(dp) , intent(inout), dimension(NPARM) :: x 
  
  Real(dp) , parameter                    :: SF = 0.10
  Real(dp) , dimension(NPARM)             :: ass  = SF      !amoeba step size
  
  Real(dp) , dimension(:),   allocatable  :: qe
  Real(dp) , dimension(:,:), allocatable  :: qdir !direction vector
  Real(dp) , dimension(:),   allocatable  :: y
  Real(dp) , dimension(:,:), allocatable  :: pt
  
  Integer(I4B) :: i, j, iters
  Real(dp)     :: fval
  
  
  iters = 0
  
  if (.not. allocated(qe))   allocate(qe(NPARME))
  if (.not. allocated(qdir)) allocate(qdir(NPARME,NPARME))
  if (.not. allocated(pt))   allocate(pt(NPARME+1,NPARME))
  if (.not. allocated(y))    allocate(y(NPARME+1))
  
  q = x  !initialize parameters, q scope is NFXPP
  !	print '(A,4F10.5)',"Debug(NFXP) parms:  ",q(:)
  j = 1
  do i=1,NPARM
     if (ss(i) == 1) then
        qe(j) = q(i)
        j=j+1
     end if
  end do
  
  if (sum(ss) /= NPARME) then
     print *,"Error:: Mismatch between number of parameters to estimate &
          & and selection vector."
     print *,"NPARME = ",NPARME
     print *,"Selection vector= ",ss
     print *,"Job Stopped."
     stop
  end if
  
  !	set initial direction vector
  qdir(1:NPARME,1:NPARME) = 0.0_dp
  !	print *,"qdir: shd be all zeros."
  !	print '(4F5.1)',(qdir(i,:),i=1,4)
  
  forall (i=1:NPARME) qdir(i,i) = 1.0_dp
  
  !	build basis and set up  y
  !	print *,"Debug(NFXP): build basis."
  !	print *,"qdir: "
  !	do i=1,NPARME
  !		print '(4F5.1)',qdir(i,:)
  !	end do
  
  !	print *,"ass: "
  !	print '(4F4.1)',ass(:)
  !	print *," "
  pt(1,:) = qe  !vector containing parameters
  y(1) = Func( pt(1,:) )
  j=2
  do i=1,NPARM
     if (ss(i) == 1) then
        if (ass(i) == 0) then
           print *,"ERROR:: ass(",i,") is zero, must be nonzero."
           print *,"Change control deck and restart."
           print *,"Job Stopped."
           stop
        end if
        !			print '(A,4F10.5)',"qdir(j-1,:): ",qdir(j-1,:)
        !			print '(A,4F10.5)',"ass          ",ass(:)
        pt(j,:) = pt(1,:) + qdir(j-1,:)*ass(i)  !ass amoeba step size
        !			print '(A,4(F10.5,1X))',"Debug(NFXP):: pt( 1,:): ",pt(1,:)
        !			print '(A,I2,A,4(F10.5,1X))',"Debug(NFXP):: pt(",j,",:): ",pt(j,:)
        y(j) = Func(pt(j,:)) ! to be consistent with Amoeba, j=-loglike
        j=j+1
     end if
  end do
  
  CALL amoeba(pt,y,TOL1,Func,iters)
  
  !	iters = Number of Amoeba function evaluations
  mcr%nfxpit = iters
  
  !  loglikelihood value is in -y(1)
  mcr%nfxpfcn = -y(1)
  !	fval = y(1)
  fval = 0.0_dp
  
  ! transfer estimated parameters to q vector
  qe = pt(1,:)
  j=1
  do i=1,NPARM
     if (ss(i) == 1) then
        q(i) = qe(j)
        if (j == NPARME) exit !do loop
        j = j+1
     end if
  end do
  mcr%nfxpq(:) = q 
  x = q  !define parameter vector to send out 
  
  !  no need to calculate standard errors
  !	CALL Newton(qe,fval,Func)
  !	deallocate(f0)
  deallocate(y)
  deallocate(qe)
  deallocate(qdir)
  deallocate(pt)
end Subroutine NFXP
!------------------------------------------------------------
Subroutine NPL(MAX_KK,CCP,x)
!  use IMSLF90
  use amvars
  use NPLP    !routines for NPL
  use nr
  IMPLICIT NONE
  
  Integer  , intent(in) :: MAX_KK
  Real(dp) , intent(in), dimension(MSTATE,CAP_J) :: CCP
  Real(dp) , intent(inout), dimension(NPARM)    :: x
  
  Interface
     Subroutine UTIL(x,pflag)
       use amvars
       IMPLICIT NONE
       Real(dp) , intent(in), dimension(NPARM) :: x
       Integer(I4B), intent(in)                :: pflag
     end Subroutine UTIL
  end interface
  
  Real(dp) , parameter                    :: SF = 0.10
!  Real(dp) , parameter                    :: macheps = EPSILON(SF)
!  Real(dp) , parameter                    :: TINY = 1.0e-10
  Real(dp) , dimension(MSTATE,CAP_J)  :: CCP0, CCP1 !work
  Real(dp) , dimension(MSTATE,MSTATE) :: FUP
  Real(dp) , dimension(MSTATE,MSTATE) :: f 
  Real(dp) , dimension(MSTATE,MSTATE) :: tFUP, I_M, MTM
  Real(dp) , dimension(MSTATE,NPARM)  :: H, HP
  Real(dp) , dimension(MSTATE)        :: p, tet,b
  

  Real(dp) , dimension(NPARM)             :: ass  = SF      !amoeba step size
  Real(dp) , dimension(:),   allocatable  :: qe
  Real(dp) , dimension(:,:), allocatable  :: qdir !direction vector
  Real(dp) , dimension(:),   allocatable  :: y
  Real(dp) , dimension(:,:), allocatable  :: pt
  
  Integer(I4B) , dimension(MSTATE)        :: indx
  Integer(I4B) , dimension(1:2)           :: tti  !location of max diff
  
  Integer(I4B) :: kk, i, a, j, iters
  Real(dp) :: supdif, fval,d,tsize
  
  
  !	allocate vectors for estimation
  !	Estimate parameters using initial estimate of CCP, ht, et etc.
  if (.not. allocated(qe))   allocate(qe(NPARME))
  if (.not. allocated(qdir)) allocate(qdir(NPARME,NPARME))
  if (.not. allocated(pt))   allocate(pt(NPARME+1,NPARME))
  if (.not. allocated(y))    allocate(y(NPARME+1))
  if (.not. allocated(f0))   allocate(f0(NRECS))
  !  All vectors are initialized, now start iterations
  
  CCP0 = CCP  !input choice probs
  
  q = x 		!assign parameters


!  if (abs(q(i)) < TINY) then !trap coef = 0
!     tsize = 0.5_dp !assumed scaling to within unit interval
!  else 
!     tsize = abs(q(i))
!  end if
 
 
  indx(:) = 0 
  
  estimate: do kk=1,MAX_KK
     tFUP = 0.0
     FUP  = 0.0
     IFUPI = 0.0
     do a=1,CAP_J
        p = CCP0(:,a)  !1xMSTATE (from sumstat)
        f = ff(:,:,a) !MxM
        tFUP = Had2(p,f)  !set for debugging purposes
        FUP = FUP + tFUP
     end do
     
     I_M(:,:)  = 0.0
     forall (i=1:MSTATE) I_M(i,i) = 1.0 !identity matrix
     MTM = I_M - BETA*FUP

!     CALL DLINRG (MSTATE,MTM,MSTATE,IFUPI,MSTATE)  !DLINRG double precision
!  DLINRG is the IMSL routine to find the inverse of a general maxtrix 
!  find inverse of MTM using Numerical Recipes

     call ludcmp(MTM,indx,d)
     do a=1,MSTATE
        b(:) = 0.0
        b(a) = 1.0_dp
        call lubksb(MTM,indx,b) !build inverse column by column
        IFUPI(:,a) = b
     end do
        
! check accuracy of inverse
!     I_M = MATMUL(IFUPI,MTM)
!     print *,"Product of MTM and inverse shd be identity matrix."
!    do a=1,MSTATE
!        write(unit=*,fmt='(5F10.5)') (I_M(a,j),j=1,MSTATE)
!     end do
     HP = 0.0
     do a=1,CAP_J
        p = CCP0(:,a)   !1xM
        H = HH(:,:,a) !MSTATE x NPARM
        HP = HP + Had2(p,H) 
     end do !a
     
     ! calculate direct product of P and (econst - log(p))
     tet = 0.0         !1xM
     do a=1,CAP_J
        p = CCP0(:,a)  !1xM
        do j=1,MSTATE
           if (p(j) > 0.0_dp) then  !check PP > 0
              if (p(j) < TINYP) then
                 print *,"(NPL)  p(j) < TINYP, p(j)=",p(j)," reset to TINYP"
                 p(j) = TINYP
              end if
              lp(j,a) = econst - log(p(j))  !natural log
           else
              lp(j,a) = 0.0
           end if
        end do
        tet = tet + Had1(p,lp(:,a))
     end do !a
     
     ! redefine MTM
     ! calculate h-tilde vector (ht) ht : M x NPARM * CAP_J 
     do a=1,CAP_J
        f = ff(:,:,a)
        H = HH(:,:,a)
        ! redefine MTM
        MTM = 0.0 
        MTM = BETA*matmul(f,IFUPI)
        ht(:,:,a) = H + matmul(MTM,HP)
        ! calculate e-titde vector 
        et(:,a) = matmul(MTM,tet) 
     end do !a
     
     j=1
     do i=1,NPARM
        if (ss(i) == 1) then
           qe(j) = q(i)
           j=j+1
        end if
     end do
     
     !	set initial direction vector
     qdir(1:NPARME,1:NPARME) = 0.0_dp
     
     !		print *,"Debug(NPL):: qdir, shd be all zeros."
     !		print '(4F5.1)',(qdir(i,:),i=1,NPARME)
     
     forall (i=1:NPARME) qdir(i,i) = 1.0_dp
     
     !		print *,"Debug(NPL):: qdir, shd be I matrix."
     !		print '(4F5.1)',(qdir(i,:),i=1,NPARME)
     
     !	build basis and set up  y
     pt(1,:) = qe  !vector containing parameters
     !		print '(A,4F10.5)',"Debug(NPL):: pt( 1,:)=",pt(1,:)
     !	to be consistent with Amoeba, define y as -loglike
     y(1) = Func( pt(1,:) )
     !		print *,"(NPL): Loglikelihood evaluate: ",-y(1)	
     j=2
     do i=1,NPARME
        pt(j,:) = pt(1,:) + qdir(j-1,:)*ass(i)  !ass amoeba step size
        !			print '(A,I2,A,4F10.5)',"Debug(NPL):: pt(",j,",:)= ",pt(j,:)
        y(j) = Func(pt(j,:)) ! to be consistent with Amoeba, j=-loglike
        j=j+1
     end do
     
     CALL amoeba(pt,y,TOL1,Func,iters)
     
     qe=pt(1,:)
     j=1
     do i=1,NPARM
        if (ss(i) == 1) then
           q(i) = qe(j)
           if (j==NPARME) exit !do loop
           j=j+1
        end if
     end do
     
     !	Calculate Value Function associated with parameters q, transition matrix ff
     !  and utility as represented by HH.
     
     CCP1 = 0.0_dp !initialize for call to VF
     
     !calculates value function at new parms generate new CCP1
     CALL VF(q,CCP0,CCP1) 
     
     supdif = maxval(abs(CCP0-CCP1))
     tti(:) = maxloc(abs(CCP0-CCP1))
     if ( supdif < TOL2) then
        CCP0 = CCP1
        exit estimate
     end if
     CCP0 = CCP1
  end do estimate
  
  !	kk = = Number of PI-K iterations
  
  mcr%nplit = kk
  mcr%nplfcn = -y(1)
  mcr%nplq(:) = q
  
  x=q !return estimated parameters to main
  fval = 0.0_dp 
  !	CALL Newton(qe,fval,Func)
  deallocate(f0)
  deallocate(qdir)
  deallocate(qe)
  deallocate(pt)
  deallocate(y)
end Subroutine NPL
!-----------------------------------
Subroutine STARTS
  use amvars
  IMPLICIT NONE
  
  Real(dp), dimension(MSTATE,CAP_J) 	 :: Pi
  Real(dp), dimension(MSTATE,MSTATE)   :: gam
  Real(dp), dimension(2:3)             :: rho 
  Integer(I4B) , parameter             :: ssa = 20  !span of sex activity
  Integer(I4B) :: i, cx, a, nx, j
  
  rho(2) = 0.005_dp;  rho(3) = 0.01_dp
  Pi(:,:) = 0.0_dp  !set default at 0.0
  do cx=1,MSTATE
     do a = 2,3
        Pi(cx,a) = (1.0_dp - rho(a))**(cx-1)
     end do
  end do
  a = 4          !dummy argument
  Pi(MSTATE,a) = 1.0_dp
  Pi(1:MSTATE-1,a) = 0.0_dp
  
  !	print *,"Debug(Starts):: Pi matrix."
  !	print '(4F7.4)',(Pi(cx,:),cx=1,MSTATE)
  
  gam(:,:) = 0.0_dp !set default at 0.0, prob distribution on sexual activity
  !	gam(cx,nx) = for current state cx, probability of sex level nx 
  do cx=1,MSTATE-1
     nx = min(MSTATE-1,cx+ssa)
     if (cx > 1) then
        gam(cx,cx:nx) = 1.0_dp/real(nx-cx+1)
     else !cx == 1
        gam(cx,cx+1:nx) = 1.0_dp/real(nx-cx)
     end if
  end do
  
  !	print *,"Debug(Starts):: gam matrix."
  !	print '(5F7.4)',(gam(cx,:),cx=1,MSTATE)
  
  ! build transition matrix for each decision variable a
  !	initialize ff(x,x',a)
  ff(:,:,:) = 0.0_dp
  !  pregnant
  ff(MSTATE,MSTATE,4) = 1.0_dp  !absorbing state for all actions a
  
  !	a=1 no sex
  ff(1,1,1) = 1.0_dp
  !	a=2 have sex use High effective BC
  !	a=3 have sex use Low effective BC
  do a=2,3
     do i=1,MSTATE-1
        do j=2,MSTATE-1
           if (i==1) then
              ff(i,j,a) = gam(i,j)
           else
              ff(i,j,a) = gam(i,j)*Pi(i,a)
           end if
        end do !j
        if (i .NE. 1) then
           ff(i,MSTATE,a) = 1.0_dp - sum(ff(i,:,a))
        end if
     end do !i
  end do !a
  !
  
  !	print *,"Debug(Starts):: ff(nx,cx,a) matrix."
  do a=1,CAP_J
     !		print *,"decision rule a=",a
     do cx=1,MSTATE
        !			print '(5F7.4)',(ff(cx,j,a),j=1,MSTATE)
     end do !cx
     !		print *," "
  end do !a
  
  !	print *," "
  
  !	Form HH matrix: HH(MSTATE,NPARM,CAP_J)
  !	Initialize and then set nonzero elements
  HH = 0.0_dp   	
  !	a=1 HH(:,:,1) remains as zero
  HH(MSTATE,4,4) = -1.0_dp   !a=4 pregnancy state
  do cx=2,MSTATE-1
     HH(cx,1,2:3) = log((real(cx)))
     HH(cx,2,2)   = -1.0_dp   !a=2
     HH(cx,3,3)   = -log((real(cx)))  !a=3
  end do
  !	print *," "
  !	print *,"Debug(Starts):: HH matrix: HH(MSTATE,NPARM,a)."
  do a=1,CAP_J
     !		print *,"Decision rule= ",a
     do i=1,MSTATE
        !			print '(4F7.3)',HH(i,:,a)
     end do
  end do
  !	print *," "
end Subroutine Starts
!----------------------------------------------------
Subroutine SUMSTAT(CCP)
  use amvars 
  IMPLICIT NONE
  Real(dp)     , intent(out), dimension(MSTATE,CAP_J)	:: CCP
  Integer(I4B) , dimension(MSTATE,MSTATE,CAP_J) 		   :: cnt
  Integer(I4B) , dimension(MSTATE,MSTATE)       		   :: scnt
  Integer(I4B) , dimension(MSTATE)					:: nn
  Integer(I4B) :: n, i, a, j
  
  cnt(:,:,:)  = 0
  scnt(:,:)   = 0
  nn(:)       = 0
  
  CCP(:,:) = 0.0_dp
  
  !	assumed structure of dat
  !		field	1	: individual number
  !		field 2	: period (w/in individual) number
  !		field 3	: current state (cx)
  !		field 4	: action (a)
  !		field	5	: next state    (nx)
  
  do n=1,NRECS
     i=dat(n,3)
     a=dat(n,4)
     j=dat(n,5)
     cnt(i,j,a) = cnt(i,j,a) + 1
     nn(i) = nn(i) + 1
  end do !n
  
  scnt(:,:) = sum(cnt(:,:,:),dim=3)
  
  do i=1,MSTATE
     do a=1,CAP_J
        CCP(i,a) = sum(cnt(i,:,a))/real(nn(i))
     end do !a
  end do !i
  
  !	print *,"Debug(sumstat):: empirical choice probs."
  !	do i=1,MSTATE
  !		print '(4F10.5)',CCP(i,:)
  !	end do
  !	print *," "
  !	print '(A,I10)',"cnt(MSTATE,MSTATE,CAP_J)= ",cnt(MSTATE,MSTATE,CAP_J)
end Subroutine sumstat
!---------------------------------------------------------------------------
SUBROUTINE UTIL(q,pflag)
  use amvars
  IMPLICIT NONE
  Real(dp), intent(in), dimension(NPARM) :: q !parameters
  Integer(I4B), intent(in)	:: pflag !controls printing

!	a component for the mbe sequence of programs
!	Fill u(cx,a) - a global variable. 
!		cx is the current state  
!		a  is the choice selected 

  Real(dp) , parameter :: TINY = 1.0e-20
  Real(dp) , dimension(NPARM) :: z
  Integer(I4B) :: a, i
  
  if (pflag == 1) then
     print *,"Debug(Util) "
     print '(A)',"parameters passed: "
     print '(3E24.12)',q
     print *," "
     print *," Logit form of q(3): ",1.0_dp/(1.0_dp+exp(-q(3)))
  end if
  
  
  !	q(1) = utility  of sex
  !	q(2) = psychic cost of high contraceptive effort (fixed cost)
  !	q(3) = psychic cost of low contraceptive effort (barrier method)
  !  q(4) = utility of pregnancy (pregnancy an absorbing state)
  !
  
  u(:,:) = 0.0_dp
  
  z = q
  if (q(3) .gt. SUPDP) then
     print *,"(Util) q(3) > SUPDP, q(3)= ",q(3)," z(3) set to TINYP."
     z(3) = TINYP  
  else
     z(3) = 1.0_dp/(1.0_dp+exp(-q(3))) !logit transform of q(3)
  end if
  
  !	form utility function (x,a)
  do a=1,CAP_J
     do i=1,MSTATE
        u(i,a) = dot_product(HH(i,:,a),z)
     end do !i
  end do !a 
  
  if (pflag == 1) then
     print *,"Parameter estimates in u(x,a)."
     print '(4I10)',(i,i=1,CAP_J)
     do i=1,MSTATE
        print '(10F7.4)',u(i,:)
     end do
     print *," "
  end if !pflag
end Subroutine UTIL
!-----------------------------------------------------