!!------------------------------------------------------------------------
! $Id: american_mc.f90,v 1.29 1997/11/12 02:03:48 svmcguir Exp svmcguir $
!!Written by Scott A. Klasky
!!Copyright Scott A. Klasky, Geoffrey Fox, Mike Boughton, Scott McGuire
!! and NPAC, Syracuse University
!!Program to calculate American options
!!
!!First attempt: 10,19,1997-
!!Modified 10/21/1997 by Scott V. McGuire
!!   NOTES: Limitations on this program include the following:
!!          Only made to work with 1 underlying (will change data structures
!!          later).
!!------------------------------------------------------------------------

program american
  use PATHS
  use MOD_INT
  use input_params
  implicit none

  type (PATH),dimension(:,:),save,allocatable    :: eta,sorteta

  real*8,dimension(:,:),save,allocatable         :: payoff, stockprice
  real*8,dimension(:,:),save,allocatable         :: sortpayoff 
  integer,dimension(:,:),save,allocatable        :: sort_order

  integer                                        :: max_nn
  integer                                        :: max_iter
  integer                                        :: ntime, npath, nn
  integer                                        :: tl, local_path, iter, otl
  integer,parameter                              :: max_params = 30
  real*8,dimension(max_params)                   :: params
  character*10,parameter                         :: fname="path.stock"
  real*8                                         :: total_energy
  real*8                                         :: energy, oldenergy, newenergy
  real*8                                         ::  oldpe, oldke, newpe
  real*8                                         :: newke,deltae
  real*8                                         :: spin_new 
  real*8                                         :: r1
  real*8,parameter                               :: HALF = 0.5d0
  integer,parameter                              :: SEED = 100
  real*8,dimension(:),save,allocatable           :: min_energy 
  real*8                                         :: temp
  integer                                        :: anealing
  integer                                        :: tnum
  integer                                        :: flips
  integer                                        :: annealing
  real*8                                         :: r
  integer                                        :: i,j
  integer                                        :: jiter
  integer                                        :: accept,reject 
  integer                                        :: oldaccept
  real*8                                         :: avgpayoff,c1,c2,mdp,dp
!  integer::temps=1000000,anneals=10
  integer                                        :: temps,anneals
  real*8::tbase,minenergy,exponent,prob=0.50
  real*8,dimension(:),allocatable::minspin
  integer                                        :: ipot,jpot,kpot,JI,JJ,nnl

  !! read in the parameters (Stock Price, Strike price, etc.
  !! read in the data from the european MC code
  !! Next sort the list in the Stock Price Dimension
  !! Next find the nearest neighbor list
  !! Set up the initial spins
  !! Find the initial energy of the system

  call random_seed(SEED)

  call read_params(max_params,params)

  !! allocate the memory
  ntime = int(params(P_NTIME))
  npath = int(params(P_NPATH))
  nn    = int(params(P_NNN))
  exponent = params(P_EXP)

  allocate ( eta(ntime,npath),sorteta(ntime,npath) ) 
  allocate ( payoff(ntime,npath) ) 
  allocate ( stockprice(ntime,npath) ) 
  allocate ( sort_order(ntime,npath) )
  allocate ( sortpayoff(ntime,npath) )

  do  tl = 1 , ntime 
     do local_path = 1 , npath 
        allocate(eta(tl,local_path)%nns(nn)) 
        allocate(eta(tl,local_path)%nng(nn)) 
        allocate(eta(tl,local_path)%nnt(nn)) 
        allocate(eta(tl,local_path)%nni(nn)) 
        allocate(eta(tl,local_path)%nnj(nn)) 
        allocate(sorteta(tl,local_path)%nns(nn)) 
        allocate(sorteta(tl,local_path)%nng(nn)) 
        allocate(sorteta(tl,local_path)%nnt(nn)) 
        allocate(sorteta(tl,local_path)%nni(nn)) 
        allocate(sorteta(tl,local_path)%nnj(nn)) 
     end do !local_path
  end do !tl

  allocate(minspin(int(params(P_NPATH))))

  call read_data(fname,max_params,params,stockprice)
  open(20,file='am.out')

  write(20,*) 'input parameters:'
  do tl = 1, max_params
     write(20,*) tl,params(tl)
  end do
  temps = int(params(P_TEMPS))
  anneals = int(params(P_ANNEALS))

  call sort_stock(max_params,params,stockprice,sort_order,payoff)
  call calculate_payoff(max_params,params,payoff,sortpayoff)
  call find_nearest_neighbor_list(max_params,params,sortpayoff,eta)
  call initialize_spins(max_params,params,sortpayoff,eta)




  jiter = 0
  c1 =1.0d0
  call random_seed(SEED+annealing) ! start off with a different random #
  call initialize_spins(max_params,params,sortpayoff,eta)
  jiter = 0 

  do otl = int(params(P_NTIME))-1 , 2,-1 
     ! Find the largest potential energy change from a spin flip
     mdp = 0
     do j=1,npath
        dp = compute_lpotential_energy_path(otl,j,max_params,params, &
             sortpayoff,eta,sort_order)
        eta(otl,j)%spin = 1-eta(otl,j)%spin
        dp = dp - compute_lpotential_energy_path(otl,j,max_params,params, &
             sortpayoff,eta,sort_order)
        eta(otl,j)%spin = 1-eta(otl,j)%spin
        dp = abs(dp)
        if(dp > mdp) then
           mdp = dp
        end if
     end do

     ! Make each kinetic energy change from a spin flip
     ! dominate potential changes
     c2 = npath*mdp/0.5

     eta(otl,:)%spin = eta(otl+1,:)%spin
     minspin = eta(otl,:)%spin
     minenergy = compute_energy_path(max_params,params,sortpayoff,eta,  &
                                     sort_order)
     do annealing = 1, anneals
        tbase = -(c2+mdp)/log(prob)
        accept = 0;reject = 0; oldaccept=100
        do tnum = 1 , temps
           temp = exp(log(tbase) - exponent*log(1.0d0*tnum))
           call random_number(r)
           local_path = 2 + int(r* (params(P_NPATH)-2) )

!           oldpe = -compute_lpotential_energy_path(otl,local_path, &
!                max_params,params,sortpayoff,eta,sort_order)

! replaced with---------------------------------------------
            oldpe = 0.0d0
          
            do jpot = 1 , int(params(P_NPATH))
               if (sort_order(otl,jpot) == local_path) exit 
               !this is the jpot that we want
            end do !jpot
!
            do kpot =1,int(params(P_NTIME))
              if(eta(kpot,sort_order(kpot,jpot))%spin == 1) then
                oldpe = oldpe +  sortpayoff(kpot,sort_order(kpot,jpot))
                exit !kpot
              end if
            end do
					  oldpe = - oldpe
!!-----------------------------------------------------
!
!           oldke = -compute_lkinetic_energy(otl,local_path,max_params, &
!                params,sortpayoff,eta)


!replaced with ------------------------------------------------


            oldke = 0.0d0

!! For local energy changes only loop over the NN list of (t,p) termi

            do nnl = 1 , eta(otl,local_path)%nnn
               jpot = eta(otl,local_path)%nni(nnl)
               kpot = eta(otl,local_path)%nnj(nnl)
               do ipot = 1, eta(jpot,kpot)%nnn
                  JI = eta(jpot,kpot)%nni(ipot)
                  JJ = eta(jpot,kpot)%nnj(ipot)
                  oldke=oldke +                  &
                       (eta(jpot,kpot)%spin-0.5d0)*(eta(JI,JJ)%spin-0.5d0)
               end do !ipot
            end do !nnl
! now we must change the energy of the t,p point too... since this is not in the! nn list
            jpot = otl 
            kpot = local_path 
            do ipot = 1, eta(jpot,kpot)%nnn
               JI = eta(jpot,kpot)%nni(ipot)
               JJ = eta(jpot,kpot)%nnj(ipot)
               oldke=oldke +                  &
                    (eta(jpot,kpot)%spin-0.5d0)*(eta(JI,JJ)%spin-0.5d0)
            end do !i
         
            oldke=-oldke / eta(jpot,kpot)%nnn
!---------------------------------------------------------
!

           spin_new = 1.0d0 - eta(otl,local_path)%spin
           eta(otl,local_path)%spin = spin_new 
!           newpe = -compute_lpotential_energy_path(otl,local_path,max_params, &
!                params,sortpayoff,eta,sort_order)

! replaced with---------------------------------------------
            newpe = 0.0d0
          
            do jpot = 1 , int(params(P_NPATH))
               if (sort_order(otl,jpot) == local_path) exit 
               !this is the jpot that we want
            end do !jpot

            do kpot =1,int(params(P_NTIME))
              if(eta(kpot,sort_order(kpot,jpot))%spin == 1) then
                newpe = newpe +  sortpayoff(kpot,sort_order(kpot,jpot))
                exit !kpot
              end if
            end do
						newpe = -newpe
!-----------------------------------------------------

!           newke = -compute_lkinetic_energy(otl,local_path,max_params, &
!                params,sortpayoff,eta)

! replaced with---------------------------------------------
            newke = 0.0d0

!! For local energy changes only loop over the NN list of (t,p) termi

            do nnl = 1 , eta(otl,local_path)%nnn
               jpot = eta(otl,local_path)%nni(nnl)
               kpot = eta(otl,local_path)%nnj(nnl)
               do ipot = 1, eta(jpot,kpot)%nnn
                  JI = eta(jpot,kpot)%nni(ipot)
                  JJ = eta(jpot,kpot)%nnj(ipot)
                  newke=newke +                  &
                       (eta(jpot,kpot)%spin-0.5d0)*(eta(JI,JJ)%spin-0.5d0)
               end do !ipot
            end do !nnl
!! now we must change the energy of the t,p point too... since this is not in the! nn list
             jpot = otl
            kpot = local_path 
            do ipot = 1, eta(jpot,kpot)%nnn
               JI = eta(jpot,kpot)%nni(ipot)
               JJ = eta(jpot,kpot)%nnj(ipot)
               newke=newke +                  &
                    (eta(jpot,kpot)%spin-0.5d0)*(eta(JI,JJ)%spin-0.5d0)
            end do !i
         
            newke=-newke / eta(jpot,kpot)%nnn
!-----------------------------------------------------

           deltae = ((c1*newpe + c2*newke) - (c1*oldpe + c2*oldke))/temp
           if( (deltae <= 0)) then
              !accept by doing nothing
					accept = accept + 1
           else
              call random_number(r)
              if((r >= exp(-deltae))) then 
                 !reject by changing it back
                 eta(otl,local_path)%spin = 1.0d0 - spin_new
                 reject = reject +1
              else
                 accept = accept + 1
              end if
           end if
           if (mod(tnum,int(params(P_DUMP))).eq.0) then
              oldenergy =  compute_energy_path(max_params,params, &
                   sortpayoff,eta,sort_order)
              write(6,133) accept,reject,temp,tnum,oldenergy,c2
              write(20,133) accept,reject,temp,tnum,oldenergy,c2
              if ( (oldaccept.eq.0).and.(accept.eq.0)) then
								write(6,*) 'going to exit'
								exit
							end if
              oldaccept = accept
							accept=0;reject=0
           end if
        end do !tnum
        jiter = jiter + 1
        call dump_dat1d("spins",jiter,jiter*1.0d0,max_params, &
                         params,sortpayoff,eta)
        avgpayoff = sorted_average_payoff(max_params,params,sortpayoff, &
                                          eta,sort_order)
        oldenergy = compute_energy_path(max_params,params,sortpayoff,eta, &
                     sort_order)
        if(oldenergy < minenergy) then
           minenergy = oldenergy
           minspin(:) = eta(otl,:)%spin
        end if
        jiter = jiter + 1
        write(6,134)  annealing,tbase,temp,oldenergy,minenergy,avgpayoff
        write(20,134) annealing,tbase,temp,oldenergy,minenergy,avgpayoff
134     format(i4,1P,5e14.5)
     end do !anneals
     eta(otl,:)%spin = minspin(:)
  end do !otl
  avgpayoff = sorted_average_payoff(max_params,params,sortpayoff, &
                                    eta,sort_order)
  oldenergy = compute_energy_path(max_params,params,sortpayoff,eta, &
                     sort_order)

  write(6,*) 'average payoff/energy=',avgpayoff,oldenergy
  write(20,*) 'average payoff/energy=',avgpayoff,oldenergy

  jiter = jiter + 1
  call dump_dat1d("spins",jiter,jiter*1.0d0,max_params, &
                  params,sortpayoff,eta)

133   format(2i8,1P,e12.2,i15,2e12.2)
  close(120)

end program american
