! pedsim.f90:          Creates a pedigree file for genomic predictions
!
!                      Fortran program written by Paul VanRaden 2007
!                      Animal Improvement Programs Lab, USDA
!                      Paul.VanRaden@ars.usda.gov
!
! OPTION DEFINITIONS:         (reads from pedsim.options file)
!                fullsib   = 0 for half sibs
!                          = 1 for full sibs
!                          = 2 for cross-classified parents
!                          =-1 for unrelated animals
!                ancgen    = 0 for sibs of animal
!                          = 1 for sibs of parents
!                ancstrs   =   number of ancestors (4, 8, etc.)
!                avgREL    =   average reliability of EBVs
!                avgREL2   =   average reliability of later EBVs
!                ancREL    = 0 does not output ancestor EBVs
!                          = X sets ancestor reliability of EBVs
!                maxped    =   maximum animals in pedigree file
!                maxan     =   maximum animals with data
!                maxdsdnt  =   max descendants without data
!                unkgrps   =   number of unknown parent groups
!                traits    =   number of traits in phenotype.data0
  use kinds
  implicit none
  real,       allocatable :: relpa(:), relpa2(:), parel(:)
  real,       allocatable :: REL(:) ,REL2(:), diag(:)
  real(8),    allocatable :: Fvec(:)
  integer,    allocatable :: sire(:), dam(:)
  integer,    allocatable :: ped(:,:)
  integer,    allocatable :: biryr(:), chip(:)
  character*1,allocatable :: sexv(:)
  character*3,allocatable :: oldyoung(:)
  integer maxped,maxan,maxdsdnt,traits
  integer fullsib,ancgen,idnum,pnt,pnt1,pnt2
  integer ancstrs,dsdnts,trt,unkgrps
  integer Ynum,start,stop,sibs,chips
  integer i,minyear/2000/
  real    avgREL,avgREL2,ancREL
  print *,'Output from pedsim.f90'
! -----------------------------------------------  Read option file
                    open ( 8,file='pedsim.options')

  read(8,*)           fullsib,ancgen,ancstrs,avgREL,avgREL2,ancREL
  print *,        '   fullsib,ancgen,ancstrs,avgREL,avgREL2,ancREL'
  print '(3i7,3f7.3)',fullsib,ancgen,ancstrs,avgREL,avgREL2,ancREL
  read(8,*)
  read(8,*)
  read(8,*)     maxped,maxan,maxdsdnt,chips,unkgrps,traits
  print *,  '   maxped,maxan,maxdsdnt,chips,unkgrps,traits'
  print '(6i7)',maxped,maxan,maxdsdnt,chips,unkgrps,traits
! ----------------------------------------------- Output data files
                    open (10,file='pedigree.file')
                    open (11,file='phenotype.data0')
                    open (12,file='phenotype.later0')
                    open (13,file='genotype.data0')
! ---------------------------------------------------
  allocate (relpa(maxped),sexv(maxped),diag(maxdsdnt))
  allocate (ped(2,maxped),oldyoung(maxped),chip(maxped))
  allocate (REL(0:maxped),parel(0:maxped),biryr(maxped))
  allocate (REL2(0:maxped),relpa2(0:maxped))
  allocate (sire(maxped), dam(maxped), Fvec(0:maxped))
  if (fullsib == -1) then
    print *,'Simulate information from unrelated animals'
  else
    print *,'Simulate pedigrees'
    end if
  if (fullsib == 1 .and. ancgen == 0) print *,'Fullsibs of animal'
  if (fullsib == 1 .and. ancgen == 1) print *,'Fullsibs of parents'
  if (fullsib == 0 .and. ancgen == 0) print *,'Halfsibs of animal'
  if (fullsib == 0 .and. ancgen == 1) print *,'Halfsibs of parents'
  if (fullsib == 2) print *,'Cross-classified parents'
  REL = .000001
  chip = 0
  Ynum = 0
! ------------------------------------------  Generate pedigree data
  dsdnts = maxdsdnt
  if (maxped - dsdnts < ancstrs + 1) then
    print *,'Decrease maxdsdnt or increase maxped in options'
    stop
    end if
  ped = 0
  do i=1,ancstrs
    biryr(i) = minyear
    if (i == (i/2)*2) then
      sexv(i) = 'F'
    else
      sexv(i) = 'M'
      chip(i) = 1
      end if
    REL(i) = ancREL
    REL2(i) = ancREL
    oldyoung(i) = 'old'
    end do ! i
  if (ancREL > 0.) then
    print *,ancstrs,' ancestors with phenotypes'
  else
    print *,ancstrs,' ancestors without phenotypes'
    end if
!                                 Sibs with data after ancestors w/o
  if (fullsib < 2) then
    do i=ancstrs+1,maxped-dsdnts
      if (i == (i/2)*2 .or. ancgen == 0) then
        if (fullsib > -1) ped(1,i) = 1
        if (fullsib == 1) ped(2,i) = 2
      else
        if (fullsib > -1) ped(1,i) = 3
        if (fullsib == 1) ped(2,i) = 4
        end if
      end do ! i
  else
!                                           Cross-classified parents
    start = ancstrs + 1
    sibs = (maxped-dsdnts) / ((ancstrs/2)**2)
    stop = ancstrs + sibs
    do pnt1 = 1,ancstrs,2
      do pnt2 = 2,ancstrs,2
        do i = start,stop
          ped(1,i) = pnt1
          ped(2,i) = pnt2
          end do ! i
        start = stop + 1
        stop = min(stop+sibs,maxped-dsdnts)
        end do ! pnt2
      end do ! pnt1
    end if
  do i=ancstrs+1,maxped-dsdnts
    Ynum = Ynum + 1
    REL(i) = avgREL
    REL2(i) = avgREL2
    sexv(i) = 'M'
    oldyoung(i) = 'old'
    biryr(i) = minyear + 5
    chip(i) = 1
    end do ! i
  print *,Ynum,' animals with genotypes and data (REL > RELpa)'
!                                  Descendants with genotype, no data
  do i=maxped-dsdnts+1,maxped
    if (i == (i/2)*2) then
      sexv(i) = 'F'
    else
      sexv(i) = 'M'
      end if
    if (fullsib > -1) then
      if (ancgen == 0) then
        ped(1,i) = 1
        ped(2,i) = 2
      else
        ped(1,i) = ancstrs + 1
        ped(2,i) = ancstrs + 2
        end if
      if (fullsib == 2) ped(2,i) = ancstrs + (i - maxped + dsdnts + 3) / 3
!     if (fullsib == 2) ped(2,i) = maxped - dsdnts
      end if
    REL(i) = 0.
    REL2(i) = avgREL2
    oldyoung(i) = 'yng'
    biryr(i) = minyear + 10
    chip(i) = 1 + i - i / chips * chips
    end do ! i
!                                                  Compute inbreeding
  do i = 1,maxped
    sire(i) = max(ped(1,i), 0)
    dam(i)  = max(ped(2,i), 0)
    end do ! i
  Fvec = 0
  call meuw(sire,dam,Fvec,maxped)
!                                                 Write pedigree file
  do i = 1,maxped
    idnum = i + 10000
    if (ped(1,i) == 0) then
      pnt1 = 0
      if (unkgrps > 0) pnt1 = -1 - i + i/unkgrps*unkgrps
    else
      pnt1 = ped(1,i) + 10000
      end if
    if (ped(2,i) == 0) then
      pnt2 = 0
      if (unkgrps > 0) pnt2 = -i + (i-1)/unkgrps*unkgrps
    else
      pnt2 = ped(2,i) + 10000
      end if
    write(10,'(a,i12,2i9,i5,a4,f5.1,i5)') &
      sexv(i),idnum,pnt1,pnt2,biryr(i),'0101',Fvec(i)*100.,0
    end do ! i
  print *,dsdnts,' descendants with genotypes but no data'
  write(13,'(a)') 'idnum  chip'
  write(13,'(a)')  
!                                                 Write genotyped list
  do i = 1,maxped
    if (sexv(i) == ' ') sexv(i) = 'F'
    if (REL(i) > .99) REL(i) = .99
    idnum = i + 10000
!                                                  Only sire genotypes
!   if ((i > ancstrs) .or. &
!       (fullsib == 2 .and. i .ne. i/2*2)) then
!                                                     Or all genotypes
!   if (i > ancstrs .or. fullsib == 2) then
    if (i > ancstrs .or. ancREL > 0.) then
      write(11,'(i10,1x,a,1x,a,i5,99(2f10.2,2f4.0))') &
         idnum,sexv(i),oldyoung(i),traits,(0.,0. &
        ,REL(i)*100,relpa(i),trt=1,traits)
!                                            Write later phenotype REL      
      write(12,'(i10,1x,a,1x,a,i5,99(2f10.2,2f4.0))') &
         idnum,sexv(i),oldyoung(i),traits,(0.,0. &
        ,REL2(i)*100,relpa2(i),trt=1,traits)
!
      write(13,'(i10,i5)') idnum,chip(i)
      end if
    end do ! i
  print *
  print *,'Program pedsim completed normally'
  stop
  end
!--------------------------------------------------------------------
  subroutine meuw(ss,dd,f,n)
!           Ignacio Aguilar and Ignacy Misztal, University of Georgia
!                                               01/29/09 - - 07/06/10
!                     Meuwissen & Luo algorithm to compute inbreeding
   use kinds
   implicit none
   integer :: ss(n),dd(n)
   real(8) :: f(0:n)
!  integer :: ss(:),dd(:)
!  real(8) :: f(0:)
   integer, allocatable:: point(:)
   real(8), allocatable:: l(:),d(:)
   integer :: is,id,i,j,k,n,ks,kd
   real(8) :: r,fi
   real :: t1
      
!   n=size(ss,dim=1)
   allocate(point(n),l(n),d(n))
   point=0;l=0;d=0
   f(0)=-1
   do i=1,n
      is=ss(i); id=dd(i)
      ss(i)=max(is,id)
      dd(i)=min(is,id)
      d(i)=.5-.25*(f(is)+f(id))
      if (is==0 .or. id==0) then
         f(i)=0
         cycle
      endif
      fi=-1
      l(i)=1.0
      j=i
      do while (j/=0)
         k=j
         r=0.5*l(k)
         ks=ss(k)
         kd=dd(k)
         if (ks>0) then
            do while(point(k)>ks)
               k=point(k)
            enddo
            l(ks)=l(ks)+r
            if (ks/=point(k))then
               point(ks)=point(k)
               point(k)=ks
            endif
            if (kd>0)then
               do while (point(k)>kd)
                  k=point(k)
               enddo
               l(kd)=l(kd)+r
               if (kd/=point(k)) then
                  point(kd)=point(k)
                  point(k)=kd
               endif
            endif
         endif

         fi=fi+l(j)*l(j)*d(j)
         l(j)=0
         k=j
         j=point(j)
         point(k)=0
      enddo
      f(i)=fi
   enddo
   f(0)=0
end subroutine
! -------------------------------------------------------------------
