! phenosim.f90:        Simulates phenotype.data for genomic predictions
!                      Writes phenotype.data and phenotype.later
!
!                      Fortran program written by Paul VanRaden  2006
!                      Animal Improvement Programs Lab, USDA
!                      Paul.VanRaden@ars.usda.gov
!
! OPTION DEFINITIONS:         (reads from phenosim.options file)
!                maxped    =   maximum animals in pedigree file
!                maxan     =   maximum animals with data
!                maxdsdnt  =   max descendants without data
!                sameREL   = 0 uses reliabilities in phenotype.data0
!                          = n sets all REL the same as trait n
!                herit     =   assumed heritability of traits
  implicit none
  real*8, dimension(:,:),  allocatable :: BV,EBV,EBVanc,Yvec,Yvec1,Yvec2
  real*8, dimension(:),    allocatable :: BVvec
  real,   dimension(:,:),  allocatable :: PTA,PA,relpta,parel
  real,   dimension(:,:),  allocatable :: PTA2,PA2,relpta2,parel2
  real,   dimension(:,:),  allocatable :: REL,REL2,errvar,errvar2
  real,   dimension(:,:),  allocatable :: relpa,relpa2
  real,   dimension(:),    allocatable :: ptavec,relvec,pavec,parelvec
  integer,dimension(:,:),  allocatable :: ped
  integer,dimension(:),    allocatable :: Ylist,listY
  integer,dimension(:),    allocatable :: eqnlist,ancstr,ids,eqns
  integer,dimension(:),    allocatable :: listdsd,dsdlist,dsdnt
  integer,dimension(:),    allocatable :: idvec,genotyped
  character*1,dimension(:),allocatable :: sexv
  character*3,dimension(:),allocatable :: oldyoung,oldyoung2
  character*1 sex
  character*3 oldyng
  character*8 birth
  integer*8   nmiss
  integer loci,qtls,correlated
  integer maxped,maxan,maxdsdnt,maxped2
  integer seed,idnum,maxn,seed1(1)
  integer dsdnts,nhash,nbv,step,steps
  integer traits,mintrt,sameREL
  integer Ynum,prevmin,biryr
  integer n13,n16,n18,n26,trt,bigtrt/1/
  integer nped,RELdif/0/
  integer ngen,ngen2,nphen,nphen2
  integer i,j,nn,pnt,pnt1,pnt2,locus,anim,n12,RELdrop/0/
  real    herit/.25/
  real    DEprog,regr,regr2,random_normal,rand
  real    var,err1,err2,uniform
  real*8  clength
! -----------------------------------------------  Read option file
  print *,'Output from phenosim.f90'

  open ( 9,file='phenosim.options')
  read(9,*)     maxped,maxan,maxdsdnt,sameREL
  print *,  '   maxped,maxan,maxdsdnt,sameREL'
  print '(7i7)',maxped,maxan,maxdsdnt,sameREL
!
! -----------------------------------------------  Input data files
! sort pedfile by ascending birth year !
                    open (10,file='pedigree.file')
                    open (11,file='phenotype.data0')
                    open (12,file='phenotype.later0')
                    open (13,file='hap.filled')
                    open (14,file='genomic.bv')
                    open (25,file='simseed')
! ...................................................  Output files
                    open (17,file='phenotypic.dyd')
                    open (18,file='phenotype.later')
                    open (19,file='phenotype.data')
                    open (21,file='trait.names')
! -------------------------------------------------  Declare arrays
  read (25,*)
  read (25,*) seed, loci, qtls, trt, correlated
  maxped2 = 2 * maxped
  allocate (eqns(0:maxped2))
  allocate (Ylist(maxan),listY(maxped),ped(2,maxped),ids(maxped2))
  allocate (eqnlist(maxped),ancstr(maxan))
  allocate (listdsd(maxped),dsdlist(maxdsdnt),dsdnt(maxdsdnt))
  allocate (sexv(maxped))
  allocate (idvec(maxped),genotyped(maxped))
  allocate (oldyoung(maxped),oldyoung2(maxped))
  oldyoung  = 'yng'
  oldyoung2 = 'yng'
! ------------------------------ Read # traits from phenotype.data0
  read(11,*,end=990) &
      idnum,sexv(1),oldyoung(1),traits
  rewind 11
  print *
  print *,traits,' traits in phenotype.data0'
  mintrt = min(traits,10)
  if (sameREL > 0 .and. traits > 1) print *, &
    'All traits simulated with same Reliability as trait',sameREL
  if (sameREL > traits) go to 988
!                                          Write trait.names file
  do trt = 1,traits
    if (trt < 10) then
      write(21,'(i2,1x,a,i1,12x,f5.2)') trt,'Trait',trt,herit
    else
      write(21,'(i2,1x,a,i2,11x,f5.2)') trt,'Trait',trt,herit
      end if
    end do ! trt
  allocate (PTA(maxped,traits),PA(maxped,traits))
  allocate (PTA2(maxped,traits),PA2(maxped,traits))
  allocate (EBV(maxdsdnt,traits),EBVanc(maxan,traits))
  allocate (REL(0:maxped,traits),REL2(0:maxped,traits))
  allocate (Yvec(maxped,traits))
  allocate (relpta(maxped,traits),parel(maxped,traits))
  allocate (relpta2(maxped,traits),parel2(maxped,traits))
  allocate (relpa(maxped,traits),relpa2(maxped,traits))
  allocate (errvar(maxped,traits),errvar2(maxped,traits))
  allocate (BV(traits,maxped))
  allocate (Yvec2(maxped,traits))
  allocate (ptavec(traits),relvec(traits))
  allocate (pavec(traits),parelvec(traits))
  allocate (BVvec(traits))
!                                                          Set seed
  seed1(1) = seed
  if (seed > 0) call random_seed (put = seed1)
!                                                  Absoft statement
! if (seed > 0) uniform = rand(seed)
  REL = .000001
  listdsd = 0
  dsdlist = 0
  dsdnt = 0
  Ynum = 0
  eqns(0) = 0
  nn = nhash(0,-1,ids,maxped2,maxn)
! ----------------------------------------------  Read pedigree data
  nped = 0
  eqnlist = 0
  do
    read(10,*,end=3) sex,idnum,pnt1,pnt2,birth
    read(birth,'(i4)') biryr
    nped = nped + 1
    sexv(nped) = sex
    ped(1,nped) = pnt1
    ped(2,nped) = pnt2
    if(pnt1 <= 0) ped(1,nped) = 0
    if(pnt2 <= 0) ped(2,nped) = 0
    idvec(nped) = idnum
!                                                    Hash individual
    nn = nhash(idnum,1,ids,maxped2,maxn)
    if (maxn .ne. nped) then
      print *,'Duplicate IDs in pedigree.file !'
      stop
      end if
    eqns(nn) = nped
    end do
 3 print *
   print *,nped,' records read from pedigree file'
!                                                       Hash parents
  do i = 1,nped
    do pnt = 1,2
      ped(pnt,i) = eqns(nhash(ped(pnt,i),0,ids,maxped2,maxn))
      if (ped(pnt,i) >= i) print *,'progeny ',i &
        ,' must precede parent ',ped(pnt,i),' in pedigree.file'
      if (ped(pnt,i) >= i) ped(pnt,i) = 0
      end do ! pnt
    end do ! i
! -------------------------------------------------  Read genomic.bv
  nbv = 0
  do
    read (14,*,end=5) idnum,traits,(bvvec(trt),trt=1,traits)
    nbv = nbv + 1
    nn = nhash(idnum,0,ids,maxped2,maxn)
    if (nn == 0) then
      print *,idnum,' ID in genomic.bv ' &
             ,'not found in pedigree.file for nbv ',nbv
      cycle
      end if
    do trt = 1,traits
      bv(trt,eqns(nn)) = bvvec(trt)
      end do
    end do
 5 print *
  print *,nbv,' animals read from genomic.bv'
! ----------------------------------  Read list of genotyped animals
  genotyped = 0
  ngen = 0
  read(13,*,end=6) 
  read(13,*,end=6)
  do
    read(13,*,end=6) idnum
    ngen = ngen + 1
    nn = nhash(idnum,0,ids,maxped2,maxn)
    if (nn == 0) then
      print *,idnum,' ID in hap.filled ' &
             ,'not found in pedigree.file for ngen ',ngen
      stop
      end if
    eqnlist(ngen) = eqns(nn)
    genotyped(eqns(nn)) = ngen
    end do
 6 print *
  print *,ngen,' animals read from hap.filled'
! --------------------------------------------  Read phenotype.data0
  dsdnts = 0
  nphen = 0
  print *,'First ten and last animal in phenotype.data0'
  print *,'  N       ID       PTA     PA     REL      RELpa     RELprog'
  do
    nphen = nphen + 1
    read(11,*,end=7) &
      idnum,sex,oldyng,traits &
     ,(ptavec(trt),pavec(trt),relvec(trt),parelvec(trt) &
     ,trt=1,traits)
    nn = nhash(idnum,0,ids,maxped2,maxn)
    if (nn == 0) then
      print *,idnum,' ID in phenotype.data0 ' &
             ,'not found in pedigree.file for nphen ',nphen
      cycle
      end if
    sexv(eqns(nn)) = sex
    oldyoung(eqns(nn)) = oldyng
    do trt = 1,traits
!                                 Give all traits same REL (optional)
      if (sameREL > 0) then
         relvec(trt) = relvec(sameREL)
         parelvec(trt) = parelvec(sameREL)
         end if
!                                         Check REL of only 1st trait
      if (relvec(trt) .ne. relvec(1)) RELdif = 1
!                                          Adjust REL to remove RELpa
!                                           b = (DE-DEpa)/(DE-DEpa+k)
      DEprog = relvec(trt) / (100. - relvec(trt)) &
             - parelvec(trt) / (100. - parelvec(trt))
      regr = DEprog / (DEprog + 1.)
      if (regr <= 0.) then
        if (regr < 0.) print *,'Negative regression ! = ',regr
        regr = .000001
        end if
!                                          Error ratio for diagonal
      errvar(eqns(nn),trt) = (1./regr) - 1.
      REL(eqns(nn),trt) = regr
      relpa(eqns(nn),trt) = parelvec(trt) / 100.
!                                                      Store PTA, REL
      pta(eqns(nn),trt) = ptavec(trt)
      relpta(eqns(nn),trt) = relvec(trt)
      pa(eqns(nn),trt) = pavec(trt)
      parel(eqns(nn),trt) = parelvec(trt)
      end do ! trt
!                                             List animals to predict
    if (oldyng == 'yng') then
      dsdnts = dsdnts + 1
      regr = 0.
    else
!                                              List animals with data
      Ynum = Ynum + 1
      end if
    if(nphen <= 10) print '(2i10,2f10.2,3f6.2,i9)' &
      ,nphen,idnum,ptavec(1),pavec(1)  &
      ,relvec(1),parelvec(1),regr*100
    end do ! read genotyped
 7  nphen = nphen - 1
  if (nphen > 10) print '(2i10,2f10.2,3f6.2,i9)'   &
     ,nphen,idnum &
     ,ptavec(1),pavec(1),relvec(1),parelvec(1),regr*100.
  print *
  print *,nphen,' records read from phenotype.data0'
  if (RELdif == 1) print *,'WARNING: trait 1 reliability used for all'
!
! ----------------------------  Read later reliabilies of same animals
  n12 = 0
  nphen2 = 0
  print *
  print *,'First ten and last animal in phenotype.later0'
  print *,'  N       ID       PTA     PA     REL      RELpa     RELprog'
  do
    nphen2 = nphen2 + 1
    read(12,*,end=8) &
      idnum,sex,oldyng,traits &
     ,(ptavec(trt),pavec(trt),relvec(trt),parelvec(trt) &
     ,trt=1,traits)
    n12 = n12 + 1
    nn = nhash(idnum,0,ids,maxped2,maxn)
    if (nn == 0) then
      print *,idnum,' ID in phenotype.later0 ' &
             ,'not found in pedigree.file for nphen2 ',nphen2
      cycle
      end if
    oldyoung2(eqns(nn)) = oldyng
    do trt = 1,traits
!                                 Give all traits same REL (optional)
      if (sameREL > 0) then
         relvec(trt) = relvec(sameREL)
         parelvec(trt) = parelvec(sameREL)
         end if
!                                         Check REL of only 1st trait
      if (relpta(nphen2,trt) .ne. relpta(nphen2,1)) RELdif = 1
!                                          Adjust REL to remove RELpa
!                                           b = (DE-DEpa)/(DE-DEpa+k)
      DEprog = relvec(trt) / (100. - relvec(trt)) &
             - parelvec(trt) / (100. - parelvec(trt))
      regr = DEprog / (DEprog + 1.)
      if (regr <= 0.) then
        if (regr < 0.) print *,'Negative regression ! = ',regr
        regr = .000001
        end if
!                                          Error ratio for diagonal
      errvar2(eqns(nn),trt) = (1./regr) - 1.
      REL2(eqns(nn),trt) = regr
      relpa2(eqns(nn),trt) = parelvec(trt) / 100.
!                                                Store later PTA, REL
      pta2(eqns(nn),trt) = ptavec(trt)
      relpta2(eqns(nn),trt) = relvec(trt)
      pa2(eqns(nn),trt) = pavec(trt)
      parel2(eqns(nn),trt) = parelvec(trt)
      end do ! trt
    if(nphen2 <= 10) print '(2i10,2f10.2,3f6.2,i9)'     &
     ,nphen2,idnum &
     ,ptavec(1),pavec(1),relvec(1),parelvec(1),regr*100.
    end do ! read ped
 8  nphen2 = nphen2 - 1
  if (nphen2 > 10) print '(2i10,2f10.2,3f6.2,i9)'        &
     ,nphen2,idnum &
     ,ptavec(1),pavec(1),relvec(1),parelvec(1),regr*100.
  print *
  print *,nphen2,' records read from phenotype.later0'
! -----------------------------------------   Compute error from REL
  do i = 1,nped
    do trt = 1,traits
!                                               Prevent REL decreases
      if (parel2(i,trt) < parel(i,trt)) &
        parel2(i,trt) = parel(i,trt) + .000001
      if (relpta2(i,trt) < relpta(i,trt)) then
        if (sameREL == 0 .or. trt == sameREL) then
          RELdrop = RELdrop + 1
          if (RELdrop <= 10) print * &
          ,'Unexpected decrease in REL for trait ',trt,' from ' &
          ,relpta(i,trt),' to ',relpta2(i,trt),' for ID ',idnum
          end if
        relpta2(i,trt) = relpta(i,trt) + .000001
        end if
      end do ! trt
    do trt = 1,traits
!                                             List animals to predict
      if (relpta2(i,trt) - parel2(i,trt) < .1) then
        regr = 0.
      else
!                                          Adjust REL to remove RELpa
!                                           b = (DE-DEpa)/(DE-DEpa+k)
        DEprog = relpta2(i,trt) / (100. - relpta2(i,trt)) &
               - parel2(i,trt) / (100. - parel2(i,trt))
        regr = DEprog / (DEprog + 1.)
        if (regr <= 0.) regr = .000001
!                                          Error ratio for diagonal
        errvar2(i,trt) = (1./regr) - 1.
        REL2(i,trt) = regr
!       print *,'DEprog regr =',DEprog,regr
        if (regr < 0.) print *,'Negative regression ! = ',regr
        end if
      relpa2(i,trt) = parel2(i,trt) / 100.
      end do ! trt
    if (RELdrop > 0) print *,RELdrop,' unexpected decreases in REL'
    end do 
  print *
  n13 = 0
  n16 = 0
! ---------------------------------------------------------------------
  print *,ynum,' animals had REL - RELpa > .001'
!                                                   Generate phenotypes
  do i = 1,nped
    do trt = 1,traits
      if (REL(i,trt) > REL2(i,trt) .and. REL2(i,trt) > 0.) then
        print *,'Unexpected decrease in REL from ' &
          ,REL(i,trt),' to ',REL2(i,trt),' for ID ',idvec(i)
!       stop
        end if
      if (REL2(i,trt) > .99) REL2(i,trt) = .99
      if (REL2(i,trt) > .000001) then
        err2 = random_normal() * sqrt(1./REL2(i,trt) - 1.)
        Yvec2(i,trt) = BV(trt,i) + err2
      else
        err2 = 0.
        Yvec2(i,trt) = 0.
        REL2(i,trt) = .000001
        end if
      if (REL(i,trt) > .99) REL(i,trt) = .99
      if (REL(i,trt) < REL2(i,trt)) then
        err1 = random_normal() * sqrt(1./REL(i,trt) - 1./REL2(i,trt))
        Yvec(i,trt) = BV(trt,i) + err1 + err2
        if (REL(i,trt) <= .0) Yvec(i,trt) = Yvec2(i,trt)
      else
        Yvec(i,trt) = Yvec2(i,trt)
        end if
      end do ! trt
    end do ! i
! -----------------------------------------------  Write phenotypic.dyd
  do i = 1,nped
    write(17,'(i9,i4,99f10.3)') idvec(i) &
      ,traits,(Yvec(i,trt),Yvec2(i,trt),trt=1,traits)
    n13 = n13 + 1
    end do ! i
! -------------------------------------------  Compute parent averages
  pta = 0.
  pta2 = 0.
  pa = 0.
  pa2 = 0.
  parel = 0.
  parel2 = 0.
  do i = 1,nped
    do trt = 1,traits
      do pnt = 1,2
        if (ped(pnt,i) > 0) then
          pa(i,trt) = pa(i,trt) + .5 * pta(ped(pnt,i),trt)
          pa2(i,trt) = pa2(i,trt) + .5 * pta2(ped(pnt,i),trt)
          parel(i,trt) = parel(i,trt) + .25 * relpta(ped(pnt,i),trt)
          parel2(i,trt) = parel2(i,trt) + .25 * relpta2(ped(pnt,i),trt)
           end if
        end do ! pnt
!                                                Compute PTA from Yvec
      DEprog = relpta(i,trt) / (100. - relpta(i,trt)) &
             - parel(i,trt) / (100. - parel(i,trt))
      regr = DEprog / (DEprog + 1.)
      DEprog = relpta2(i,trt) / (100. - relpta2(i,trt)) &
             - parel2(i,trt) / (100. - parel2(i,trt))
      regr2 = DEprog / (DEprog + 1.)
      if (oldyoung(i) == 'old') then
        pta(i,trt) = regr * (Yvec(i,trt) &
                   - pa(i,trt)) + pa(i,trt)
      else
        pta(i,trt) = pa(i,trt)
        relpta(i,trt) = parel(i,trt)
        end if
      if (oldyoung2(i) == 'old') then
         pta2(i,trt) = regr2 * (Yvec2(i,trt) &
                     - pa2(i,trt)) + pa2(i,trt)
      else
        pta2(i,trt) = pa2(i,trt)
        relpta2(i,trt) = parel2(i,trt)
        end if
      end do ! trt
    end do ! i
!                                                Write phenotype.data
  n18 = 0
  do i = 1,nped
!   if (relpta(i,1) > 0. .or. oldyoung(i) == 'yng') then
    if (genotyped(i) > 0) then
      write(19,'(i10,1x,a,1x,a,i5,99(2f10.2,2f4.0))') &
         idvec(i),sexv(i),oldyoung(i),traits &
        ,(pta(i,trt),pa(i,trt),relpta(i,trt),parel(i,trt) &
        ,trt=1,traits)
!                                               Write phenotype.later
      write(18,'(i10,1x,a,1x,a,i5,99(2f10.2,2f4.0))') &
         idvec(i),sexv(i),oldyoung2(i),traits &
        ,(pta2(i,trt),pa2(i,trt),relpta2(i,trt),parel2(i,trt) &
        ,trt=1,traits)
      n18 = n18 + 1
      end if
    end do ! i
  print *,n18,' phenotype.data  records written to unit 19'
  print *,n18,' phenotype.later records written to unit 18'
  print *
  print *,'Program phenosim completed normally'
    stop
 988 print *,'Option sameREL exceeds number of traits in data'
    stop
 990 print *,'phenotype.data0 file is missing'
    stop
 996 print *,'Unexpected end of file, unit 14 (genomic.bv)'
    stop
 998 print *,'Increase maxan or decrease anmls'
    stop
 999 print *,'Unexpected end of file, unit 10 (pedigree.file)'
  end
! -----------------------------------------------------------------------
! Function from NAG
! This file contains:
! 1. The module RANDOM for generating random deviates from a number of
!    continuous and discrete distributions.   The module also includes a
!    simple quicksort routine.
! 2. A program to demonstrate its use.   The program includes code for the
!    log-gamma function.
! The module uses the intrinsic RANDOM_NUMBER function supplied with the
! compiler.   It has only been tested using the Lahey LF90 compiler.

! The two parts are separated by a line of exclamation marks.

! Alan Miller

! The compilers own random number generator, SUBROUTINE RANDOM_NUMBER(r),
! is used to provide a source of uniformly distributed random numbers.

REAL FUNCTION random_normal()

! Adapted from the following Fortran 77 code
!      ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
!      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
!      VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.

!  The function random_normal() returns a normally distributed pseudo-random
!  number with zero mean and unit variance.

!  The algorithm uses the ratio of uniforms method of A.J. Kinderman
!  and J.F. Monahan augmented with quadratic bounding curves.

IMPLICIT NONE
REAL                      :: zero = 0.0, half = 0.5, one = 1.0, two = 2.0

!     Local variables
REAL     :: s = 0.449871, t = -0.386595, a = 0.19600, b = 0.25472,           &
            r1 = 0.27597, r2 = 0.27846, u, v, x, y, q, rand

!     Generate P = (u,v) uniform in rectangle enclosing acceptance region

DO
  u = rand(0)
  v = rand(0)
  call random_number (u)
  call random_number (v)
  v = 1.7156 * (v - half)

!     Evaluate the quadratic form
  x = u - s
  y = ABS(v) - t
  q = x**2 + y*(a*y - b*x)

!     Accept P if inside inner ellipse
  IF (q .LT. r1) EXIT
!     Reject P if outside outer ellipse
  IF (q .GT. r2) CYCLE
!     Reject P if outside acceptance region
  IF (v**2 .LT. -4.0*LOG(u)*u**2) EXIT
END DO

!     Return ratio of P's coordinates as the normal deviate
random_normal = v/u
RETURN

END FUNCTION random_normal
! --------------------------------------------------------------------
integer function nhash (id,mode,idtab,isize,maxn)
  integer mode,nprime/433/,ie/16/,i,j,isize,maxn &
    ,idtab(isize),id
!                                   instructions for hashing:
!                           mode=-1 initiates parameters, zeroes idtab
!                               = 0 looks up id only
!                               = 1 puts new id's into idtab
  if (mode == -1) then
    maxn = 0
    do i=1,isize
      idtab(i) = 0
    enddo
    nhash = -1
    return
  endif
  if (id == 0) then
    nhash = 0
    return
  endif
  nhash = abs(mod(nprime*id,isize)) + 1
  do j = 1,isize
    if (idtab(nhash) == 0) then
      if (mode == 1) then
        idtab(nhash) = id
        maxn = maxn + 1
        if (maxn .gt. .9*isize) then
          print *,'hash table size',isize &
           ,'  is too small, increase it'
          stop 261
        endif
      else
        nhash = 0
      endif
      return
    endif
    if (idtab(nhash) == id) return
    nhash = mod(nhash+ie,isize) + 1
  enddo
  print *,'Function NHASH dropped out of the loop'
  stop 41
end
