!!----------------------------------------------------------------------------
!! $Id: find_nn_list.f90,v 1.16 1997/11/25 15:56:21 scott Exp $ 
!! subroutine find_nearest_neighbor_list 
!!    This subroutine will find the nearest neighbor lists
!!    Author: Scott A. Klasky
!!    Date  : 10/20/97
!!    The metric function that we will use by default will help compute
!!    the nearest neighbor list
!!----------------------------------------------------------------------------
subroutine find_nearest_neighbor_list(max_params,params,payoff1,payoff2,eta)
  use PATHS
  use input_params
  implicit none
  integer,intent(IN)                             :: max_params
  real*8,dimension(max_params),intent(IN)        :: params
  real*8,dimension(:,:),intent(IN)               :: payoff1,payoff2
  type (PATH),dimension(:,:),intent(INOUT)       :: eta

  integer                                        :: metric_toggle

  interface
     subroutine standard_metric(max_params,params,payoff1,payoff2,eta)
       use PATHS
       use input_params
       implicit none
       integer,intent(IN)                             :: max_params
       real*8,dimension(max_params),intent(IN)        :: params
       type (PATH),dimension(:,:),intent(INOUT)       :: eta
       real*8,dimension(:,:),intent(IN)               :: payoff1,payoff2

     end subroutine standard_metric
  end interface


  metric_toggle = int(params(P_METRIC))

  select case(metric_toggle)

  case(1)
     call standard_metric(max_params,params,payoff1,payoff2,eta)
     write(6,*) 'nnnnn',eta(1,1)%nnn

  case default
     write(*,*)'Invalid value for metric type' 
     stop

  end select


  return
end subroutine find_nearest_neighbor_list

!!----------------------------------------------------------------------------
!! subroutine   standard_metric
!!    This subroutine will read in parameters used for the MC iterations
!!    Author: Scott A. Klasky
!!    Date  : 10/20/97
!!    for this two dimensional case,
!!----------------------------------------------------------------------------
subroutine standard_metric(max_params,params,payoff1,payoff2,eta)
  use PATHS
  use input_params
  implicit none
  integer,intent(IN)                             :: max_params   
  real*8,dimension(max_params),intent(IN)        :: params
  type (PATH),dimension(:,:),intent(INOUT)       :: eta
  real*8,dimension(:,:),intent(IN)               :: payoff1,payoff2

  integer                                        :: nn
  integer                                        :: i,j       !! eta(i,j)
  integer                                        :: inn,jnn   !! for nn calcs 
  real*8                                         :: r, time, stock1,stock2
  integer,parameter                              :: maxn=1000
  real*8,dimension(maxn)                         :: g
  integer,dimension(maxn)                        :: is,js,index
  integer                                        :: ipoints, l, icount
  real*8,parameter                               :: CONST=0.5d0

  open(22,file='nnlist.out')

  nn= int(params(P_NNN))
  eta(:,:)%nnn = 0

  do j = 1, int(params(P_NPATH))
     do i = 1, int(params(P_NTIME)) 
        stock = payoff(i,j)
        time  = 1.0d0 * i

        ipoints = 0
        do jnn = max(1,j - nn), min(int(params(P_NPATH)),j + nn)
           do inn = i,i ! Only equal time neighbors
           !do inn = max(1,i - nn), min(int(params(P_NTIME)),i + nn)
              if ( (inn.eq.i).and.(jnn.eq.j)) then
              else if (inn.eq.i) then
                 ipoints = ipoints + 1

                 g(ipoints) =                                               &
                      1.0d0*(time - 1.0d0*inn)**2 +          &
!                      params(P_METCO)*(stock- payoff(inn,jnn))**2 
                      params(P_METCO)*(j- jnn)**2 

                 is(ipoints) = inn
                 js(ipoints) = jnn
                 index(ipoints) = ipoints
              end if
           end do !inn
        end do !jnn
        ! Now I must sort this 
        call batcher(g,index,ipoints)
        ! Now I put the first `nn` back in the nearest neighbor list
        icount = 0
        do l = 1, ipoints
           if ( icount.lt.nn ) then
              icount = icount + 1
              eta(i,j)%nni(l) = is(index(l))
              eta(i,j)%nnj(l) = js(index(l))
              eta(i,j)%nns(l) = g(index(l))
              eta(i,j)%nng(l) = 1.0d0 
           else
              eta(i,j)%nnn = icount
              exit
           end if
        end do
        if (eta(i,j)%nnn.eq.0) then !just take 2 closest points
           eta(i,j)%nnj(1) = max(1,j-1)
           eta(i,j)%nnj(2) = min(int(params(P_NPATH)),j+1)
           eta(i,j)%nni(1) = i 
           eta(i,j)%nni(2) = i 
           eta(i,j)%nns(1) = 1.0d0
           eta(i,j)%nns(2) = 1.0d0
           eta(i,j)%nng(1) = 1.0d0
           eta(i,j)%nng(2) = 1.0d0
           eta(i,j)%nni(1) = max(1,i-1)
           eta(i,j)%nni(2) = min(int(params(P_NTIME)),i+1)
           eta(i,j)%nnj(1) = j
           eta(i,j)%nnj(2) = j 
           eta(i,j)%nns(1) = 1.0d0
           eta(i,j)%nns(2) = 1.0d0
           eta(i,j)%nng(1) = 1.0d0
           eta(i,j)%nng(2) = 1.0d0
           eta(i,j)%nnn    = 4
        end if
     end do !i
  end do !j

  close(22)
  close(66);close(67)
  return
end subroutine standard_metric
