!------------------------------------------------------------------------------
!                        Program GIN4IMPRINT
!
! Auxiliary program for WOMBAT to set up inverse of gametic relationship matrix
! Files written out allow imprinting effects to be fitted with a GIN covariance
! structure:
!   1. The .gin file -> default name: imprint.gin
!   2. The corresponding .codes file -> default name: imprint.codes
!      Gametes are coded as: animal ID * 10 + 1 (paternal) and 
!                            animal ID * 10 + 2 (maternal)
!                            => need to add corresponding column to data file!
! Input required:
!   Pedigree file (after pruning) with animal, sire and dam ID recoded in
!   running order and animals' inbreeding coefficients;
!   this can easily be generated by WOMBAT by running a set-up step for a
!   corresponding model without the imprinting effect(s) fitted
!
! Code for gametic relationship matrix written by Bruce Tier
!-------------------------------------------------------------km----02/2011-----
 
! =======================
  program Gin4Imprinting
! =======================
 
  implicit none
 
  integer, parameter :: ibig = 2147483647
! Output files for WOMBAT, edit names as desired
  character(len=13)  :: ginfile   = 'imprint.gin  ', &
&                       codesfile = 'imprint.codes'
! Pruned pedigree file with inbreeding coefficients written out by WOMBAT
  character(len=18)  :: ppedfile  = 'PrunedPedFile1.dat'
  integer            :: nanis, ii, ian, idad, imum, jan, jdad, jmum
  real(8)            :: xx
 
  write(*,*)'***********   Program "Gin4Imprinting"  **************'
 
! read pedigree file & count no. of animals
  open(31, file = ppedfile, status = 'old', form ='formatted', action ='read')
  do
     read(31, *, iostat=ii ) ian, idad, imum, jan, jdad, jmum, xx
     if( ii < 0 ) then
         exit ! end of file
     else if( ii > 0 ) then
         write(*,*) 'Some sort of error: reading file "',ppedfile,'"'
         write(*,*) 'Last running no. read correctly was', ian
         stop 'Check file!'
     end if
  end do
  write(*,*)'Highest animal ID found      =', ian, jan
 
! set up inverse of gametic relationship matrix
  nanis = ian;  rewind(31)
  call fbt
 
  write(*,*)'"Gametic inverse" complete'
  write(*,*)'Files written out: "',trim(ginfile),'"  and  "',trim(codesfile),'"'
  write(*,*)'******************************************************'
  stop
 
  contains
 
! ==============
  subroutine fbt
! ==============
 
  real(8), dimension(:), allocatable   :: eff, diag, pxp, pxg
  integer, dimension(:,:), allocatable :: ged, ped
  real(8)                              :: rlogl, ras, das, rad, ff, slogl
  integer                              :: npxp, i, ia, is, id, ka1, ka2, &
&                                         ks1, ks2, kd1, kd2, ngams
 
! memory requirement = nanis*14 words
  ngams = nanis*2
  allocate( pxg(0:ngams), ped(2,0:nanis), ged(0:2,0:ngams),  &
            eff(0:nanis), diag(0:ngams), pxp(0:ngams), stat=ii )
  ped = 0; ged = 0; eff = 0.d0; diag = 0.d0; pxp = 0.d0; pxg = 0.d0
 
! read in pedigree  & inbreeding and store elements
  eff(0)= - 1.d0; rlogl = 0.d0
 
! read in pedigree and store elements
  do i = 1, nanis
     read(31,*) ia, is, id, jan, jmum, jdad, ff
     xx = jan * 10 + 2
     if( xx > dble(ibig) ) then
         write(*,*)'Original ID too big for coding gametes as "ID*10+i"'
         write(*,*) jan
         stop
     end if
     ped(1,i) = is
     ped(2,i) = id
     eff(i) = ff
!    codes for gametes -> original ID*10+i; i=1: paternal, i=2: maternal
     ka2 = ia * 2; ka1 = ka2 - 1
     ged(0,ka1) = 10 * jan + 1  
     ged(0,ka2) = 10 * jan + 2
!    gametes contributions 
     ras = 2.d0 / (1.d0 - eff(is) )
     rad = 2.d0 / (1.d0 - eff(id) )
     rlogl = rlogl + log(ras) + log(rad)
     diag(ka1) = ras
     diag(ka2) = rad
     if( is > 0 ) then         ! sire gametes
        ks2 = is * 2
        ks1 = ks2 - 1
        ged(1,ka1) = ks1
        ged(2,ka1) = ks2
        ras = -0.5d0*ras
        pxg(ka1) = ras
        ras = -0.5d0*ras
        pxp(ks2) = pxp(ks2) + ras
        diag(ks1) = diag(ks1) + ras
        diag(ks2) = diag(ks2) + ras
    end if
    if( id > 0 ) then          ! dam gametes
       kd2 = id * 2
       kd1 = kd2 - 1
       ged(1,ka2) = kd1
       ged(2,ka2) = kd2
       rad = -0.5d0*rad
       pxg(ka2) = rad
       rad = -0.5d0*rad
       pxp(kd2) = pxp(kd2)+rad
       diag(kd1) = diag(kd1)+rad
       diag(kd2) = diag(kd2)+rad
    end if
    if( mod(ia,50000) == 0 ) write(*,'(2hAt,3(i12,2i8))') ia, ped(:,ia), ka1, &
&                                               ged(1:,ka1), ka2, ged(1:,ka2)
  end do
  close(31)
 
! write out files for WOMBAT
  open(408, file = codesfile, status = 'unknown', form = 'formatted', &
&                                                 action = 'write' )
  open(409, file = ginfile,  status = 'unknown', form = 'formatted',   &
&                                                 action = 'write' )
 
  write(409,*) -rlogl, -rlogl/log(10.0)   ! log determinant
  npxp = 0
  do i = 1, ngams
     write(408,'(i8,i15)') i, ged(0,i)
     if( pxg(i) /= 0.d0 ) then
         write(409,'(2i8,f16.8)') ged(1,i), i, pxg(i)
         write(409,'(2i8,f16.8)') ged(2,i), i, pxg(i)
         npxp = npxp + 2
     end if
     if( mod(i,2) == 0 .and. pxp(i) /= 0.d0) then
         write(409,'(2i8,f16.8)') i-1, i, pxp(i)
         npxp = npxp + 1
     end if
     write(409,'(2i8,f16.8)') i, i, diag(i)
  end do
  close(408); close(409) 
  write(*,*)'No. of gametes               =', ngams
  write(*,*)'No. of elements in inverse G =', ngams + npxp
 
  return
  end subroutine fbt
 
  end program Gin4Imprinting