!!------------------------------------------------------------------------
! $Id: american_mc.long.f90,v 1.1 1997/11/12 02:03:48 svmcguir Exp $
!!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).
!!------------------------------------------------------------------------

!! Want to define the Data structure before we begin, because this will 
!! make it much easier to code with the proper data structure.

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, ol, 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, oldpe, oldke, newpe, newke,deltae
real*8                                         :: spin_new 
real*8                                         :: r1
real*8,parameter                               :: HALF = 0.5d0
integer,parameter                              :: SEED = 100

!!SVM 10/21/1997
!real*8,dimension(:,:,:),save,allocatable       :: eta_save  !! Too much memory
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, ool
integer                                        :: accept,reject 
real*8                                         :: avgpayoff
real*8,dimension(:,:),save,allocatable         :: covm

!! 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))

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

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

   allocate (min_energy(int(params(P_ANNEALS))))

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

!   covm(0,0) = params(P_SIGMA00) !! This should probably come in 
!   covm(0,1) = params(P_SIGMA01) !! with the paths.
!   covm(1,0) = params(P_SIGMA10)
!   covm(1,1) = params(P_SIGMA11)

   write(20,*) 'input parameters:'
   do tl = 1, max_params
      write(20,*) tl,params(tl)
   end do

   if (int(params(P_SORT)) <> 0) then
      call sort_stock(max_params,params,stockprice,sort_order,payoff)
   else
      write(6,*) 'no sorting'
      write(20,*) 'no sorting'
      payoff = stockprice 
   end if
   call calculate_payoff(max_params,params,payoff,sortpayoff)

   call find_nearest_neighbor_list(max_params,params,sortpayoff,covm,eta)
	write(6,*) 'nnn',eta(1,1)%nnn
   call initialize_spins(max_params,params,sortpayoff,eta)

   if (int(params(P_PE_PATH)) == 0) then
      oldenergy =  compute_energy(max_params,params,       &
                                  sortpayoff,eta)
   else
      oldenergy =  compute_energy_path(max_params,params,  &
                                       sortpayoff,eta,sort_order)
   end if

   write(6,*) 0,oldenergy
   write(20,*) 0,oldenergy

!! Now we are ready to start a Metropolis algorithm
!! Metropolis works as follows:
!! We will loop for a max # of iterations (defined in a file)
!!    Loop over time (from end time to 0 time)
!!       Loop over stocks: which are sorted
!!          Compute the energy from the nearest neighbor list
!!          Change the \eta_{ij}
!!          Compute the new energy from the nearest neighbor list
!!          if the energy goes down then
!!            accept the change 
!!          else 
!!            pick a random # from 0 to 1
!!            if the random # is less than 0.5 then
!!               accept the change
!!            else
!!               reject the change
!!            end if
!!          end if
!!       end loop over stock prices
!!   end loop over time


!! Open questions are how we loop, acceptance criteria, nearest neighbor
!! calculations (how many defines a nearest neighbor, how do we calculate
!! the metric (for stock price dimension)

!! for right now, calculate the entire energy, later only calculate the
!! energy terms that are affected by a change

   if(int(params(P_SORT)) == 2) then
      avgpayoff = sorted_average_payoff(max_params,params,sortpayoff,eta,sort_order)
   else
      avgpayoff = average_payoff(max_params,params,sortpayoff,eta,sort_order)
   end if
   write(6,*) 'average payoff=',avgpayoff
   write(20,*) 'average payoff=',avgpayoff

   jiter = 0
	call dump_data1d("svst",0.0d0,max_params,params,stockprice)
	call dump_data1d("pvst",0.0d0,max_params,params,sortpayoff)
!    call dump_data("eta",0.0d0,max_params,params,eta)
!   call dump_datar("payoff",0.0d0,max_params,params,sortpayoff)

   do annealing = 1, int(params(P_ANNEALS))

      write(6,*) 'annealing=',annealing
      write(20,*) 'annealing=',annealing
      call random_seed(SEED+annealing) ! start off with a different random #
      call initialize_spins(max_params,params,sortpayoff,eta)
      do tnum = 1 , int(params(P_TEMPS)) 
         write(6,*) '=============================temp loop',tnum
         write(20,*) '=============================temp loop',tnum
         accept = 0
         reject = 0
         if (tnum == 1) then
             temp = params(P_INIT_TEMP)
         else
              temp = params(P_INIT_TEMP)/log(1.0d0*tnum)
         end if
         do ool = 1, int(params(P_OUTL))
            write(6,*) 'ool=',ool
            write(20,*) 'ool=',ool

            do otl =  int(params(P_NTIME))-1, 1 , -1
               write(6,*) 'otl=',otl
               do iter = 1 , int(params(P_MCITERS))
                  call random_number(r)
                  tl = min(int(r*(1+params(P_WINDOW)))+otl,                  &
                           int(params(P_NTIME))-1) 
                  call random_number(r)
                  ol = int(r*params(P_NPATH)) + 1 
                  if (int(params(P_PE_PATH)) == 0) then

                     oldpe = compute_lpotential_energy(tl,ol,max_params,      &
                                                      params,sortpayoff,eta)
                     oldke = compute_lkinetic_energy(tl,ol,max_params,        &
                                                    params,sortpayoff,eta)
                  else
                      oldpe = compute_lpotential_energy_path(tl,ol,max_params, &
                                                       params,sortpayoff,eta, &
                                                       sort_order)
                      oldke = compute_lkinetic_energy(tl,ol,max_params, &
                                                       params,sortpayoff,eta) &
                  end if
                  spin_new = 1.0d0 - eta(tl,ol)%spin
                  eta(tl,ol)%spin = spin_new 

                  if (int(params(P_PE_PATH)) == 0) then

                     newpe = compute_lpotential_energy(tl,ol,max_params,      &
                                                      params,sortpayoff,eta)
                     newke = compute_lkinetic_energy(tl,ol,max_params,        &
                                                    params,sortpayoff,eta)
                  else
                      newpe = compute_lpotential_energy_path(tl,ol,max_params, &
                                                       params,sortpayoff,eta, &
                                                       sort_order)
                      newke = compute_lkinetic_energy(tl,ol,max_params, &
                                                       params,sortpayoff,eta) &
                  end if
                  deltae = (-params(P_C1)*newpe -params(P_C2)*newke/temp) - (-params(P_C1)*oldpe - params(P_C2)*oldke/temp)
                  if( (deltae <= 0)) then
                     !accept by doing nothing
                  else
                     call random_number(r)
                     if((r >= exp(deltae))) then 
                        !reject by changing it back
                        eta(tl,ol)%spin = 1.0d0 - spin_new
                        reject = reject +1
                     else
                        accept = accept + 1
                     end if
                  end if
                  if (mod(iter,10000)==0) then
                     if (int(params(P_PE_PATH)) == 0) then
                        oldenergy =  compute_energy(max_params,params,       &
                                                    sortpayoff,eta)
                     else

                        oldenergy =  compute_energy_path(max_params,params,  &
                                                    sortpayoff,eta,sort_order)
                     end if
                     write(6,*)  annealing,tnum,ool,iter,oldenergy
                     write(20,*) annealing,tnum,ool,iter,oldenergy
                  end if
               end do !iter
               if (int(params(P_PE_PATH)) == 0) then
                  oldenergy =  compute_energy(max_params,params,       &
                                              sortpayoff,eta)
               else
                  oldenergy =  compute_energy_path(max_params,params,  &
                                                   sortpayoff,eta,sort_order)
               end if
               write(6,*)  annealing,tnum,ool,iter,oldenergy
               write(20,*) annealing,tnum,ool,iter,oldenergy

               if(int(params(P_SORT)) == 2) then
                  avgpayoff = sorted_average_payoff(max_params,params,       &
                                                    sortpayoff,eta,sort_order)
               else
                  avgpayoff = average_payoff(max_params,params,sortpayoff,   &
                                             eta,sort_order)
               end if
               write(6,*) 'average payoff=',avgpayoff
               write(20,*) 'average payoff=',avgpayoff
               write(6,*) 'accept,reject=',accept,reject
               write(20,*) 'accept,reject=',accept,reject
            end do !otl
         end do !ool
      end do !tnum
      jiter = jiter + 1
      call dump_data("eta",1.0d0*jiter,max_params,params,eta)

      if (int(params(P_PE_PATH)) == 0) then
         min_energy(annealing) = compute_energy(max_params,params,           &
                                                sortpayoff,eta)
      else
         min_energy(annealing) = compute_energy_path(max_params,params,      &
                                                    sortpayoff,eta,sort_order)
      end if
      write(6,*) annealing,min_energy(annealing)
      write(20,*) annealing,min_energy(annealing)
   end do !annealing

    !Now we should pick out the result with the lowest energy and do
    !something with it
   close(20)

end program american
