!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  
!       Program "SECATEURS" makes `pruning' of pedigrees easy !
!
!  'Pruning' reduces the number of levels to be fitted for genetic
!            effects, thus reducing the memory requirements
!
!  When fitting genetic effecs, we include animals which do not have
!  records themselves but are parents of animals represented in the data.
!  Using further pedigree information, we generally also include parents
!  of parents, parents of grand-parents, etc. 
!
!  Any parent without records connected to only one other animal in the
!  pedigree does not add any information, e.g. a parent with a single
!  offspring only and unknown parents him/herself. For the purpose of
!  analysis, this parent can be treated as unknown, and thus eliminated 
!  from the analysis. This is equivalent to 'absorbing' the equation for 
!  this animal, and referred to as `pruning'. For pedigrees spanning several
!  generations backwards, this should be done repeatedly. For direct genetic
!  effects, pruning is done downwards, i.e. scanning the pedigree from 
!  oldest to youngest animals.
!
!  Frequently, we fit parental genetic effects, most commonly maternal
!  effects. For these, only animals which have an offspring in the data
!  have a `record'. Hence animals in the data which do not have such
!  progeny are candidates for pruning, i.e. pruning needs to be carried out
!  upwards as well as downwards the pedigree. This can lead to substantial
!  numbers of parental animal effects which can be disregarded. If the
!  software used for the mixed model analysis allows separate pedigrees
!  or inverses of the numerator relationship matrix to be used for different
!  genetic efffects (e.g. ASREML allows for up to 6 NRM inverse files to
!  be specified), pedigrees for maternal of paternal effects should thus
!  be pruned specifically for these effects.
!
!  NB : Pruning carried out by "SECATEURS" assumes direct genetic and
!       paternal genetic effects are uncorrelated. 
!       Similarly, do not use "SECATEURS" if parental effects are fitted
!       and separate pedigrees/NRM inverses cannot be fitted.
!                                                 @c karin meyer, 2003

!  Open source version: Use algorithm of Meuwissen & Luo to calculate
!       inbreeding coefficients rather than Tier's procedure;
!       this version copes with > 23 generations              km 03/09
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

MODULE parameters
   integer, parameter :: maxan=1444555  ! max. no. of animal IDs
END MODULE

MODULE filenames
   character(len=30)  :: pedfile ='pedigree.dat                  '
   character(len=30)  :: datfile ='records.dat                   '
   character(len=30)  :: parfile ='prune.par                     '
   character(len=13), dimension(2) :: nrmfile =(/'nrminv.dat   ','nrminvmat.dat'/),&
 &                                    outfile =(/'pednew.dat   ','pednewmat.dat'/)
   character(len=14)  :: datnewfile='recordsnew.dat'
END MODULE

MODULE pedigreelist
   integer, save                              :: iadd, ndumm
   integer, dimension(:), allocatable, save   :: idvec,jdvec,kdvec, mdvec
   real(8), dimension(:,:), allocatable, save :: inbreed
END MODULE pedigreelist
  
MODULE pedlisttwo
   integer, dimension(:), allocatable, save   :: nannew
   integer, dimension(:,:), allocatable, save :: idwec,jdwec,kdwec
END MODULE pedlisttwo

MODULE positions
   integer, save                              :: ngeneff
   integer, dimension(2), save                :: ipos
END MODULE positions

MODULE options
   logical, save                              :: verbose, setup, karin, nrmout,&
 &                                               maternal, binary, recode,     &
 &                                               meuwissen 
   integer, save                              :: nfit, recl, maxgen
   integer, dimension(10), save               :: kfit
END MODULE options
                                                            ! end of modules

! ==================
  program SECATEURS
! ==================

  use pedigreelist; use pedlisttwo; use options; use filenames; use positions
  implicit none
  integer :: nanim, ndam, ii, jj, mm0,iun=99, mm
  logical :: lexist

!  get input options
   call cmdline

!  log file
   open(iun,file='PRUNE.log',status='unknown',form='formatted')
   write(iun,'(1x,15(1h*),a,18(1h*),2hkm,3(1h*)/)')                        &
 &                         '  Output from program "SECATEURS"  '
   call todayis (iun)

!  set up list of pedigrees
   call pedlist (nanim, iun)

!  declare arrays for ped info
   allocate(idwec(0:nanim,ngeneff), jdwec(nanim,ngeneff),   &
&           kdwec(nanim,ngeneff), nannew(ngeneff), stat=ii)
            if(ii/=0)stop 'alloc pedlst2'

! copy 'original' pedigree list
  do jj=1,ngeneff
     idwec(1:nanim,jj)=idvec(1:nanim)
     jdwec(1:nanim,jj)=jdvec(1:nanim)
     kdwec(1:nanim,jj)=kdvec(1:nanim)
  end do
  idwec(0,:)=0

! "prune" pedigrees
  do jj = 1, ngeneff
     call prunepeds (nanim, jj, ipos(jj), iun )  ! jj=effect no.
  end do
  deallocate( idvec, jdvec, kdvec, stat=ii); if(ii/=0)stop 'dealloc pedlst'

! add running codes to data file
  if(recode) then
    if(ngeneff==1)then
         mm=ipos(1)
    else
         mm=maxval(ipos)
    end if
    call addcodes (mm, recl)
  end if

  mm0=maxval(nannew);  allocate(inbreed(0:mm0,ngeneff),stat=ii)
  if(ii/=0)stop 'alloc inbreed'

! set up nrm-inverse
  if(nrmout)then
      if(karin) open(11,file='detl.dat')
      do jj = 1,ngeneff
         call inbreedng(jj,nannew(jj),iun); if( meuwissen ) cycle
         call a_inverse(nannew(jj),jj,iun)
      end do
  end if
  write(iun,'(1x,25(1h*),a,18(1h*),2hkm,3(1h*)/)')'END of "SECATEURS" '
  stop 'SECATEURS : all done !'

  end program secateurs

! ========================
  subroutine cmdline
! ========================

  use options; use positions; use filenames
  implicit none
  integer           :: iargc, i, k, l
  character(len=30) :: arg
  logical           :: lexist
  external getarg

! default values for options      
  verbose =.false.; nrmout=.false.; maternal=.false.; karin=.false.
  binary =.false.;  recode=.false.; meuwissen = .true. !.false.
  ngeneff=1; nfit=0; ipos=(/ 1, 3 /); maxgen = 0 !23
 
  inquire(file=parfile,exist=lexist); if (lexist) call readparfile

! no. of arguments
  k=iargc()
  if( k==0 .and. (.not.lexist))then
      write(*,'(a)')'Useage : secateurs [OPTIONS]'
      write(*,*)' OPTIONS available are :'
      write(*,*)'  -Ppedigree-file-name '
      write(*,*)'  -Ddata-file-name   '
      write(*,*)'  -Iposition-of-animal-ID'
      write(*,*)'  -MATposition-of-dam-ID '
      write(*,*)'  -N'
      write(*,*)'  -Gmax-no-generations  '
      stop 'Read the manual !'
  end if

  do i=1,k
     call getarg (i,arg); arg=adjustl(arg)
     call arg2opt
  end do

  return

  contains

! ======================
  subroutine readparfile
! ======================

  integer :: ii

  open(19,file=parfile,status='old',form='formatted',action='read')
  do 
     read(19,'(a)',iostat=ii)arg; if(ii/=0)exit
     call arg2opt
  end do
  close(19)

  return
  end subroutine readparfile

!     ==================
      subroutine arg2opt
!     ==================

!     write out NRM inverse
      if(arg(1:2).eq.'-n'.or.arg(1:2).eq.'-N')then
         nrmout=.true.

!     binary output for NRM inverse
      else if(arg(1:2).eq.'-b'.or.arg(1:2).eq.'-B')then
         binary = .true.
         nrmfile(1) = 'nrminv.bin   '
         nrmfile(2) = 'nrminvmat.bin'

!     maximum no. of generations
      else if(arg(1:2).eq.'-g'.or.arg(1:2).eq.'-G')then
         read(arg(3:),*) maxgen
         if( maxgen < 1 .or. maxgen > 31 ) then
             meuwissen = .true.
             write(*,*)'Using Meuwissen & Luo algorithm to calculate inbreeding'
         end if

!     data file
      else if(arg(1:2).eq.'-d'.or.arg(1:2).eq.'-D')then
         read(arg(3:),'(a)')datfile

!     pedigree file
      else if(arg(1:2).eq.'-p'.or.arg(1:2).eq.'-P')then
         read(arg(3:),'(a)')pedfile

!     prune maternal effects pedigree
      else if(arg(1:2).eq.'-m'.or.arg(1:2).eq.'-M')then
         maternal=.true.; ngeneff=2
         if(arg(3:4)=='at'.or.arg(3:4)=='aT'.or.arg(3:4)=='At'.or.           &
      &                                              arg(3:4)=='AT')then
            read(arg(5:),*)ipos(ngeneff)
         else
            read(arg(3:),*)ipos(ngeneff)
         end if

!     position of animal ID in data file
      else if(arg(1:2).eq.'-i'.or.arg(1:2).eq.'-I')then
         read(arg(3:),*)ipos(1)

!     order of polynomial fit in RR analysis
      else if(arg(1:2).eq.'-l'.or.arg(1:2).eq.'-L')then
         nfit=nfit+1
         if(nfit<11)read(arg(3:),*)kfit(nfit)

      else if(arg(1:2).eq.'-k'.or.arg(1:2).eq.'-K')then
         karin=.true.

!     recode data
      else if(arg(1:2).eq.'-r'.or.arg(1:2).eq.'-R')then
         recode=.true.
         read(arg(3:),*)recl

      else if(arg(1:2).eq.'-v'.or.arg(1:2).eq.'-V')then
         verbose=.true.

!     last argument is parameter file name
      else if(arg(1:1).ne.'-'.and.i==k)then
         read(arg,'(a)')parfile
         l=len_trim(parfile)
         if(l<26.and.parfile(l-3:l).ne.'.par')parfile(l+1:l+4)='.par'
         call readparfile

      else 
         write(*,*)'option not recognised',arg
      end if

      if(trim(pedfile)==trim(datfile))then
         write(*,*)'Data & pedigree files specified are the same - invalid !'
         stop 'erroneous options'
      end if

      return
      end subroutine arg2opt

      end subroutine cmdline

!     ==========================
      subroutine todayis (iun)
!     ==========================

      integer, intent(in)     :: iun
      integer                 :: ii=0, hostnm, getcwd
      character(len=40), save :: cwdir='not determined                '

      character(len=30),save  :: machine='Host name not determined     '
      character(len=8)        :: tday
      character(len=10)       :: now

!     fortran 95 intrinsic routine (should work on all systems)
      call date_and_time(tday,now)

      if(iun.eq.0)then
        write(*,'(10a)')' Today is ',tday(7:8),'/',tday(5:6),'/',          &
     &             tday(1:4),' -- Time is ',now(1:2),':',now(3:4)
      else
        write(iun,'(10a)')' Today is ',tday(7:8),'/',tday(5:6),'/',        &
     &             tday(1:4),' -- Time is ',now(1:2),':',now(3:4)
      end if

!     write out host name
!     UNIX specific 
      ii=hostnm(machine)
!     call getenv('HOST',machine)       ! Linux, Lahey LF95 specific routine

      if(ii.eq.0.and.iun>0)then
         write(iun,'(3a)')' Running on host : "',trim(machine),'" '
      else if(ii.eq.0)then
          write(*,'(3a)')' Running on host : "',trim(machine),'" '
      end if

!     UNIX specific - comment out for other systems
      ii=getcwd(cwdir)
!     call getenv('PWD',cwdir)          ! Linux, Lahey LF95 specific routine
      write(iun,*)'Working directory : ',trim(cwdir)
      write(iun,*)' '

      return
      end subroutine todayis

!==================================================
      subroutine pedlist (nanim, iun)
!===================================================

      use parameters, only: maxan;  use filenames, only: pedfile
      use options; use pedigreelist

      implicit none
      integer, intent(out) :: nanim
      integer, intent(in)  :: iun
      integer              :: ii, ipre, ian, imum, idad, jan, jdad, jmum
      integer              :: manim, nrec
      integer              :: ibig=2147483646
      integer, dimension(:), allocatable :: isex

!     input file : column separated list of animal, sire & dam IDs
!                  --> preferably sorted according to animal

      open(1, file=pedfile, status='old', form='formatted', action='read')

!     prelim : count no. of records in ped.file to get rough idea of max.no.
      call check_peds 
      
!     allocate arrays
      allocate (idvec(0:manim), jdvec(manim), kdvec(manim), isex(manim), stat=ii)
      if(ii>0)stop 'pedlist : alloc idvec'
      
! 1st pass : set up ordered list of IDs
  rewind(1); idvec(0) = ibig;  nanim=0
  do
     read(1,*,iostat=ii)ian, idad, imum;   if(ii<0)exit
     call add_list (ian)
     if(idad>0) call add_list (idad)
     if(imum>0) call add_list (imum)
  end do
      
! 2nd pass : set up list of recoded sire & dam IDs
  jdvec=0;  kdvec=0;  isex=0
  rewind(1); ii=0; nrec=0;  ipre=-99
  do 
     read(1,*,iostat=ii)ian, idad, imum ; if(ii/=0)exit
     nrec=nrec+1; if(verbose.and.mod(nrec,5000)==0)write(*,*)'  pedigree : at',nrec
     if(ian.eq.ipre)cycle ! check for & ignore doubles
     ipre=ian; call find_anim (ian,jan, nanim)
     if(idad>0)then
         call find_anim (idad,jdad, nanim);  jdvec(jan)=jdad
         if(isex(jdad)==0)then
            isex(jdad)=1
         else if(isex(jdad)==2)then
            write(*,*)'Inconsistent pedigree !'
            write(*,*)'Animal, sire, dam',ian,idad,imum
            write(*,*)'Sire ID has previously been found as *dam* ID'
!            stop 'correct pedigrees'
         end if
      end if
      if(imum>0)then
         call find_anim (imum,jmum, nanim); kdvec(jan)=jmum
         if(isex(jmum)==0)then
            isex(jmum)=2
         else if(isex(jmum)==1)then
            write(*,*)'Inconsistent pedigree !'
            write(*,*)'Animal, sire, dam',ian,idad,imum
            write(*,*)'dam ID has previously been found as *sire* ID'
            stop 'correct pedigrees'
         end if
      end if
  end do

  write(*,*)'Pedigree file                 =   ',pedfile
  write(*,*)'Total no. of animal IDs found =',nanim
  close(1);  deallocate(isex)
  return

  contains

!     ===============================
      subroutine add_list (ian)
!     ===============================

      integer, intent(inout) :: ian
      integer                :: ii, j

      if(idvec(nanim)>ian)then
         do ii=1,nanim
         if(idvec(ii).eq.ian)then
            return              !           animal already in list
         else if(idvec(ii)>ian)then
!           new animal, need to slot into correct spot in ordered list
            nanim=nanim+1;  if(nanim>maxan)stop 'too many animals'          
            do j=nanim,ii+1,-1
               idvec(j)=idvec(j-1)
            end do
            idvec(ii)=ian; return
         end if
         end do
      end if
!     new last animal
      nanim=nanim+1;  if(nanim>maxan)stop 'too many animals'          
      idvec(nanim)=ian

      return
      end subroutine add_list

! ======================
  subroutine check_peds
! ======================

   manim=0
   do
      read(1,*,iostat=ii)ian, idad, imum

      if(ii<0)then
         exit
      else if(ii>0)then
         write(*,*)'error occurred when reading ',pedfile
         stop '"PEDLIST"'
      end if

      if(ian>ibig.or.ian<1)then
          write(*,*)'invalid animal ID',ian
          stop 'error reading pedigree'
      end if 
      if(idad>ibig.or.idad<0)then
          write(*,*)'invalid sire ID',idad
          stop 'error reading pedigree'
      end if 
      if(imum>ibig.or.imum<0)then
          write(*,*)'invalid dam ID',imum
          stop 'error reading pedigree'
      end if 
      if(ian==idad.or.ian==imum.or.(idad>0.and.imum==idad))then
          write(*,*)'invalid pedigree',ian,idad,imum
          stop 'duplicate IDs'
      end if
      if(ian<idad.or.ian<imum)then
          write(*,*)'parental IDs must be numerically smaller than animal IDs !'
          write(*,*)'invalid pedigree', ian,idad,imum
          stop 'pedigree codes'
      end if
      manim=manim+1
      if(idad>0)manim=manim+1; if(imum>0)manim=manim+1
  end do

  return
  end subroutine check_peds

  end subroutine pedlist

! =======================================
  subroutine find_anim (ian, jan, nanim)
! =======================================

  use pedigreelist

  implicit none
  integer, intent(in)   :: ian, nanim
  integer, intent(out)  :: jan
  integer               :: k1, k2, kk

  jan=0
  if( ian<idvec(1) .or. ian>idvec(nanim) )then
      write(*,*)'error : animal/dam code not found !'
      write(*,*)'ID in data =', ian
      write(*,*)'pedigree list',0,idvec(1),idvec(nanim)
      write(*,*)'exceeds min/max of list',nanim; stop 'routine "find_anim"'
  end if
  k1=1; k2=nanim+1
  do 
      kk=k1+(k2-k1)/2                                                    
      if(ian<idvec(kk))then
         k2 = kk                    
      else if(ian.eq.idvec(kk))then
         jan = kk; return
      else if(ian>idvec(kk))then
         if(k1.eq.kk)then
            write(*,*)'error : animal/dam code not found ! '
            write(*,*)'ID in data =', ian
            write(*,*)'pedigree list',kk,idvec(kk-1:kk+1)
            stop 'routine "find_anim"'
         end if
         k1 = kk
      end if
  end do

  end subroutine find_anim

! ==============================================
  subroutine prunepeds (nanim, jj, jsubj, iun )
! ==============================================

  use options; use filenames; use pedigreelist; use pedlisttwo
  use positions, only: ngeneff

  implicit none
  integer, intent(in)                :: nanim, jj, jsubj, iun
  integer, dimension(:), allocatable :: noff, nnrec, mark
  real(8), dimension(:), allocatable :: dvec
  integer                            :: ibig=2147483646
  integer                            :: ipre,ian,jan,ii,idad,imum, &
&                                       nnew, nloop, npr, nrec, nprb, i
  logical                            :: lopen

! arrays for this subroutine
  allocate(noff(nanim),nnrec(0:nanim),mark(0:nanim),dvec(jsubj), &
&          stat=ii);  if(ii>0)stop 'alloc prune'

! count no. of  records per animal
  nnrec=0; ipre=-99
!     ...  data file
  inquire(file=datfile, opened=lopen)
  if(.not.lopen)then
      open(1,file=datfile,status='old',form='formatted',action='read')
  else
      rewind(1)
  end if
  nrec=0

  do 
     read(1,*,iostat=ii) dvec;  if(ii/=0)exit
     nrec=nrec+1
     if(verbose.and.mod(nrec,10000)==0)write(*,*)'  data file : at',nrec
     ian=dvec(jsubj)
     if( ian<1 .or. ian>ibig )then
         write(*,*)'invalid "animal" ID in data file : record no.',nrec
         write(*,*)'genetic effect no.',jj,'  position no. ',jsubj
         write(*,*)'code found ',ian
         if(jj==2.and.ian==0)                                          &
     &       write(*,*)'all dam IDs must be known to fit maternal effects !'
             stop 'check data file'
         end if
         if(ian /= ipre)then
            call find_anim (ian,jan, nanim); ipre=ian
         end if
         nnrec(jan)=nnrec(jan)+1
  end do
  write(iun,'(/a,t40,1h:,1x,a)')'Data file ',datfile
  write(iun,'(a,t40,1h=,i12)')'No. of records  ',nrec
  if(.not.recode) close(1)

! count no. of offspring per animal
  noff=0
  do i = 1,nanim
     idad = jdwec(i,jj); if(idad>0) noff(idad) = noff(idad)+1
     imum = kdwec(i,jj); if(imum>0) noff(imum) = noff(imum)+1
  end do
  if(verbose)then
     if(jj==1)open(2,file='OriginalNumbers',status='unknown', form='formatted')
     write(2,'(1x,t11,a,t19,a,t25,a,t32,a,t43,a)')'animal','nrec', &
&                                     'noff','sire ID','dam ID'
     do i = 1,nanim
        write(2,'(i8,i12,2i6,2i12)')i,idvec(i),nnrec(i),noff(i),jdvec(i),kdvec(i)
     end do
     if(jj==ngeneff)close(2)
  end if

! loop through pedigree list
  nloop=0
  do
     mark=0;  nloop=nloop+1
!    prune from 'top'
     do i=1,nanim
         if( idwec(i,jj)<1   &           ! this one has already been pruned
     &       .or. nnrec(i)>0  &          ! ....... has a record
     &       .or. noff(i)>1    &         ! ....... has more than one progeny
     &       .or. jdwec(i,jj)>0 .or. kdwec(i,jj)>0 &! ... has pedigree info
     &          ) cycle
!        ... this is an animal w/o records, less than 2 progeny & unknown
!            parents -> mark it for pruning !
         mark(i)=1
      end do
!     prune from 'bottom' (e.g. maternal effects)
      do i=1,nanim
         if(mark(i)>0                &! already marked for pruning
     &        .or. idwec(i,jj)<1     &! already pruned
     &        .or. nnrec(i)>0        &! has record(s)
     &        .or. noff(i)>0         &! has offspring
     &          ) cycle
!        ... animal w/o records and progeny -> prune
         mark(i)=2
      end do
      if(count(mask=(mark==1))==0)exit ! no parents pruned
      write(*,'(a,i4,2i7)')'Pruning : loop=',nloop,                       &
&                           count(mask=(mark==1)),count(mask=(mark==2))
      npr=0; nprb=0
      do i=1,nanim
!        ... prune parents
         if(jdwec(i,jj)>0)then
            if(mark(jdwec(i,jj))==1)jdwec(i,jj)=0
         end if
         if(kdwec(i,jj)>0)then
            if(mark(kdwec(i,jj))==1)kdwec(i,jj)=0
         end if
!        ... delete animal from list
         if(mark(i)==1)then
            idwec(i,jj)=0;  npr=npr+1
         else if(mark(i)==2)then
            idwec(i,jj)=0;  nprb=nprb+1
         end if
      end do         
  end do       !     pruning complete

! recode remaining animals from 1 to N
  mark=0; nnew=0
  do i = 1,nanim
     if(idwec(i,jj)==0)cycle
     nnew=nnew+1; mark(i)=nnew
     idwec(nnew,jj)=idwec(i,jj)
     nnrec(nnew)=nnrec(i)
     noff(nnew)=noff(i)
     if(jdwec(i,jj)>0)then
        jdwec(nnew,jj)=mark(jdwec(i,jj))
     else
        jdwec(nnew,jj)=0
     end if
     if(kdwec(i,jj)>0)then
        kdwec(nnew,jj)=mark(kdwec(i,jj))
     else
        kdwec(nnew,jj)=0
     end if
  end do
  nannew(jj)=nnew
  write(*,*)'After pruning',jj,' levels',nannew(jj),nanim

! write out new pedigree list
  if(.not.nrmout)then
      open(9,file=outfile(jj),action='write',status='unknown',form='formatted')
      do i=1,nnew
         write(9,'(3i12)')idwec(i,jj), idwec(jdwec(i,jj),jj), idwec(kdwec(i,jj),jj)
      end do
      close(9)
  end if

  if(verbose)then
     if(jj==1)open(3,file='PrunedNumbers',status='unknown',form='formatted' )
     write(3,'(1x,t11,a,t19,a,t25,a,t32,a,t43,a)')'animal','nrec',   &
     &                                     'noff','sire ID','dam ID' 
     do i=1,nnew
         write(3,'(i8,i12,2i6,2i12)')i,idwec(i,jj),nnrec(i),noff(i),     &
     &                              jdwec(i,jj), kdwec(i,jj)
     end do
     if(jj==ngeneff)close(3)
  end if
  call ped_counts (nfit, jj)

  deallocate (noff, nnrec, mark, dvec, stat=ii);if(ii/=0)stop 'deall prune'
  return 

  contains

!     ===================================
      subroutine ped_counts (nfit, irand)
!     ===================================

      integer, intent(in) :: nfit, irand

      character(len=21)   :: ff='PedigreeStructure.out'
      real(8)             :: xx, yy
      integer             :: n,i,n1,mm, ifit

      write(iun,'(a,i3/66(1h=))')'Pedigree Structure for random effect :  ',irand

      xx=nannew(jj)*100.d0/nanim
      write(iun,'(a,t45,1h=,i10)')'Original no. of animals ',nanim
      write(iun,'(a,t45,1h=,i10)')'No. of animals after pruning', &
     &                                                   nannew(jj)
      write(iun,'(a,t45,1h=,f10.1/)')' ... proportion (%) remaining',xx

      n=count(mask=(nnrec(1:nannew(jj))==0))
      write(iun,'(a,t45,1h=,i10)')'No. of levels w/out records',n
      n=count(mask=(nnrec(1:nannew(jj))>0))
      write(iun,'(a,t45,1h=,i10,f7.1,1h%)')'No. of levels with records', &
     &                                                             n,100.d0
      xx=100.d0/n
      yy=100.d0/nannew(jj)
      mm=maxval(nnrec(1:nannew(jj)))
      do i=1,min(mm,6)
      n=count(mask=(nnrec(1:nannew(jj))==i))
      write(iun,'(a,i1,a,t45,1h=,i10,f7.1,1h%)')' ... ',i,' record(s) ',n, &
     &                                                              n*xx
      end do
      n=count(mask=(nnrec(1:nannew(jj))>6.and.nnrec(1:nannew(jj))<11))
      if(n>0)write(iun,'(a,t45,1h=,i10,f7.1,1h%)')' ... 7-10 record(s) ',        &
     &                                                            n,n*xx
      n=count(mask=(nnrec(1:nannew(jj))>10.and.nnrec(1:nannew(jj))<21))
      if(n>0)write(iun,'(a,t45,1h=,i10,f7.1,1h%)')' ... 11-20 record(s) ',       &
     &                                                            n,n*xx
      n=count(mask=(nnrec(1:nannew(jj))>20))
      if(n>0)write(iun,'(a,t45,1h=,i10,f7.1,1h%/)')' ... >20 record(s) ',        &
     &                                                            n,n*xx
      if(nfit>0)then
         do ifit=1,nfit
         write(iun,'(/a,t45,1h=,i10)')'Minimum no. of records specified',       &
     &                                                    kfit(ifit)
         n=count( mask=(nnrec(1:nannew(jj))>=kfit(ifit)) )
         write(iun,'(a,i2,a,t45,1h=,i10,f7.1,1h%)')                              &
     &     'No. of animals with at least ',kfit(ifit),' records',n,n*xx

!        count progeny with min. no. of records
         mark=0
         do i=1,nannew(jj)
         if(nnrec(i)<kfit(ifit))cycle
         if(jdwec(i,jj)>0)mark(jdwec(i,jj))=mark(jdwec(i,jj))+1
         if(kdwec(i,jj)>0)mark(kdwec(i,jj))=mark(kdwec(i,jj))+1
         end do
         n=count(mask=(mark(:nannew(jj))>0))
         write(iun,'(a,i2,a)')'No. of parents which have  progeny which '
         write(iun,'(a,i2,a,t45,1h=,i10)')'   have at least ',kfit(ifit),      &
     &                                                ' records',n

         n=0
         do i=1,nannew(jj)
         if(nnrec(i)>=kfit(ifit) .or.                                          &
     &     ( jdwec(i,jj)>0 .and. nnrec(jdwec(i,jj))>=kfit(ifit) ) .or.         &
     &     (kdwec(i,jj)>0 .and. nnrec(kdwec(i,jj))>=kfit(ifit)) )n=n+1       
         end do
         write(iun,'(a,i2,a)')'No. of animals with at least ',kfit(ifit),' records' 
         write(iun,'(a,t45,1h=,i10)')'   themselves or on a parent ',n

         call count_all ( nannew(jj), kfit(ifit), -1, 0, n)
         write(iun,'(a,i2,a)')'No. of animals with at least ',kfit(ifit),' records' 
         write(iun,'(a,t45,1h=,i10,f7.1,1h%)')'   themselves, on a parent or sib(s)'&
     &,                                        n,n*yy
         call count_all ( nannew(jj), kfit(ifit), 1, 0, n)
         write(iun,'(a,t45,1h=,i10,f7.1,1h%)')'   ... in the data',n,n*xx
         end do
      end if

      n=count(mask=(noff(:nannew(jj))==0))
      write(iun,'(/a,t45,1h=,i10,f7.1,1h%)')'No. of animals w/out offspring',n,n*yy
      n=count(mask=(noff(:nannew(jj))>0))
      write(iun,'(a,t45,1h=,i10,f7.1,1h%)')'No. of animals with offspring',n,n*yy
      n=count(mask=(noff(:nannew(jj))>0.and.nnrec(1:nannew(jj))>0))
      write(iun,'(a,t45,1h=,i10,f7.1,1h%/)')' ... and records ',n,n*yy

      n=count(mask=(jdwec(:nannew(jj),jj)==0))
      write(iun,'(a,t45,1h=,i10)')'No. of animals with unknown sire',n
      n=count(mask=(kdwec(:nannew(jj),jj)==0))
      write(iun,'(a,t45,1h=,i10)')'No. of animals with unknown dam',n
      n=count(mask=(jdwec(:nannew(jj),jj)==0.and. &
     &                                   kdwec(:nannew(jj),jj)==0))
      write(iun,'(a,t45,1h=,i10)')'No. of animals with both parents unknown',n
      n=count(mask=(nnrec(1:nannew(jj))>0.and.jdwec(:nannew(jj),jj)==0))
      write(iun,'(a,t45,1h=,i10)')'No. of animals with records'
      write(iun,'(a,t45,1h=,i10)')' ... and unknown sire',n
      n=count(mask=(nnrec(1:nannew(jj))>0.and.kdwec(:nannew(jj),jj)==0))
      write(iun,'(a,t45,1h=,i10)')' ... and unknown dam',n
      n=count(mask=(nnrec(1:nannew(jj))>0.and.jdwec(:nannew(jj),jj)==0 &
     &                                  .and.kdwec(:nannew(jj),jj)==0))
      write(iun,'(a,t45,1h=,i10)')' ... and both parents unknown',n

      call count_codes ( nannew(jj),-1,-1,n,jdwec(:nannew(jj),jj) )
      write(iun,'(/a,t45,1h=,i10)')'No. of sires',n
      call count_codes ( nannew(jj),1,-1,n,jdwec(:nannew(jj),jj) )
      write(iun,'(a,t45,1h=,i10)')' ... with progeny in the data',n
      call count_codes ( nannew(jj),1,1,n,jdwec(:nannew(jj),jj) )
      write(iun,'(a,t45,1h=,i10)')' ... with records & progeny in data',n

      call count_codes ( nannew(jj),-1,-1,n,kdwec(:nannew(jj),jj) )
      write(iun,'(/a,t45,1h=,i10)')'No. of dams',n
      call count_codes ( nannew(jj),1,-1,n,kdwec(:nannew(jj),jj) )
      write(iun,'(a,t45,1h=,i10)')' ... with progeny in the data',n
      call count_codes ( nannew(jj),1,1,n,kdwec(:nannew(jj),jj) )
      write(iun,'(a,t45,1h=,i10)')' ... with records & progeny in data',n

      write(iun,'(/a)')'No. of animals with known/unpruned grand-parents'
      call count_granps (nannew(jj), n, jdwec(:nannew(jj),jj), &
     &                                          jdwec(:nannew(jj),jj) )
      write(iun,'(a,t45,1h=,i10)')' ... with paternal grandsire ',n
      call count_granps (nannew(jj), n, jdwec(:nannew(jj),jj), &
     &                                          kdwec(:nannew(jj),jj) )
      write(iun,'(a,t45,1h=,i10)')' ... with paternal granddam ',n
      call count_granps (nannew(jj), n, kdwec(:nannew(jj),jj), &
     &                                          jdwec(:nannew(jj),jj) )
      write(iun,'(a,t45,1h=,i10)')' ... with maternal grandsire ',n
      call count_granps (nannew(jj), n, kdwec(:nannew(jj),jj),  &
     &                                          kdwec(:nannew(jj),jj) )
      write(iun,'(a,t45,1h=,i10/)')' ... with maternal granddam ',n

      return
      end subroutine ped_counts

! ===================================================
  subroutine count_codes (nan, lim1, lim2, nn, idvec)
! ===================================================

  integer, intent(in)                 :: nan, lim1, lim2
  integer, intent(out)                :: nn
  integer, dimension(nan), intent(in) :: idvec
  integer                             :: i, j
  integer, dimension(nan)             :: id, nd

  nn = 0
  do i=1,nan
     if(nnrec(i)<lim1)cycle          ! not enough records for animal
     if(idvec(i)<1)cycle             ! parent unknown
     if(nnrec(idvec(i))<lim2)cycle   ! not enough records for parent
     do j=1,nn
         if(id(j)==idvec(i))then
            nd(j)=nd(j)+1; go to 10
         end if
     end do
     nn=nn+1;id(nn)=idvec(i); nd(nn)=1
10   continue
  end do

  return
  end subroutine count_codes

!     ==================================================
      subroutine count_granps (nan,  nn, idvec1, idvec2)
!     ==================================================

      integer, intent(in)                 :: nan
      integer, intent(out)                :: nn
      integer, dimension(nan), intent(in) :: idvec1, idvec2

      integer                             :: i, j, ii, jj

  nn=0
  do i=1,nan
      ii=idvec1(i)
      if(ii==0) cycle         ! parent is unknown
      jj=idvec2(ii)
      if(jj==0) cycle         ! grandparent is unknown      
      if(idvec(i)==idvec(jj))then
          write(*,*)'animal & grandparent identical',idvec(i)
          write(*,*)'  ... animal no. ',i
          write(*,*)'  ... parent no. ',ii
          stop 'Pedigree error ??'
      end if
      nn=nn+1
  end do

  return
  end subroutine count_granps

! ===============================================
  subroutine count_all ( nan, kk, lim1, lim2, nn)
! ===============================================

  integer, intent(in)   :: nan, kk, lim1, lim2
  integer, intent(out)  :: nn
  integer               :: i

  nn=0
  do i = 1,nan
      if(nnrec(i)<lim1) cycle
      if(nnrec(i)>=kk .or.                                 &
     &  (jdwec(i,jj)>0 .and. nnrec(jdwec(i,jj))>=kk ) .or. &
     &  (kdwec(i,jj)>0 .and. nnrec(kdwec(i,jj))>=kk ) .or. &
     &  (jdwec(i,jj)>0 .and. mark(jdwec(i,jj))>lim2)  .or. &
     &  (kdwec(i,jj)>0 .and. mark(kdwec(i,jj))>lim2 ) ) nn=nn+1
  end do

  return
  end subroutine count_all

 end subroutine prunepeds

!========================================
 subroutine inbreedng (jj, nanim, iun)
!========================================

use pedigreelist, only: inbreed; use filenames, only: outfile
use pedlisttwo; use options, only: meuwissen, maxgen

implicit none
integer, intent(in)           :: jj, nanim, iun
real(8), dimension(0:nanim)   :: fff
integer                       :: i, ii, n
real(8)                       :: fbar

if( .not. meuwissen ) then
    call bruce_code
else
    call meuwissen_code
end if

! write out list of new pedigrees together with inbreeding coefficients
open(9,file=outfile(jj),action='write',status='unknown',form='formatted')
do i = 1,nanim
   write(9,'(3i12,f10.4)')idwec(i,jj), idwec(jdwec(i,jj),jj),          &
&                         idwec(kdwec(i,jj),jj), (fff(i))*100.d0
end do
close(9)

fbar= 100.d0* sum(fff(1:nanim))/nanim
write(iun,'(/a,i4,a)')'Inbreeding coefficients for random effect',jj,' computed'
write(*,'(/a,i4,a)')'Inbreeding coefficients for random effect',jj,' computed'
write(iun,'(a,i8)')    'No. of animals                 =',nanim
n = count(mask=(fff>0))
write(iun,'(a,i8)')    'No. of inbred animals          =',n
write(iun,'(a,f8.4,a)')'Average inbreeding coefficient =',fbar,' (in %)'
write(*,'(a,i8)')    'No. of inbred animals          =',n
write(*,'(a,f8.4,a)')'Average inbreeding coefficient =',fbar,' (in %)'
fbar = 0.d0; ii = 0
do i = 1,nanim
   if( fff(i) > fbar)then
      fbar = fff(i);  ii=i
   end if
end do
write(iun,'(a,f8.4,a)')'Highest inbreeding coefficient =',fbar*100.0, &
&                                               ' (in %)'
write(iun,'(a,i12)')   '... for animal                 =',idwec(ii,jj)

return

contains

! =====================
  subroutine bruce_code
! =====================

integer, dimension(0:nanim)        :: kink, knik
integer, dimension(2,0:nanim)      :: ped
real(8), dimension(0:nanim)        :: lll, ddd
real, dimension(:), allocatable    :: depth
integer, dimension(:), allocatable :: linped
integer                            :: ia, is, id,  i, j, k, l, ka, kc, ke, kd, &
 &                                    ks, kk, kc2, ii, mglen
real(8)                            :: r, q

write(*,*)'Tier''s algorithm not available in this version'; stop

return
end subroutine bruce_code

! =========================
  subroutine meuwissen_code
! =========================

! inbreeding program from Meuwissen and Luo (1992), GSE 24: 305-313
  implicit none
  integer                       :: ian, jan, kan, is, id, ks, kd
  real(8)                       :: fi, r
  integer, dimension(nanim)     :: point
  integer, dimension(nanim,2:3) :: ped
  real(8), dimension(nanim)     :: l, d

! initialise
  point = 0; l = 0.d0; d = 0.d0; fff = 0.d0; fff(0) = -1.d0
  ped(:,2) = jdwec(1:nanim,jj); ped(:,3) = kdwec(1:nanim,jj)

! main loop
  do ian = 1, nanim
     is = ped(ian,2); id = ped(ian,3)
     ped(ian,2) = max(is,id); ped(ian,3) = min(is,id)
     d(ian) = 0.5d0 - 0.25d0 * (fff(is) + fff(id))
     if (is.eq.0. or .id.eq.0) then
         fff(ian) = 0.d0
     else if(ped(ian-1,2).eq.ped(ian,2).and.ped(ian-1,3).eq.ped(ian,3)) then
         fff(ian) = fff(ian-1)  ! full-sib, occuring in sequnce
     else
         fi = -1.d0; l(ian) = 1.d0; jan = ian
         do while (jan > 0 )
            kan = jan; r = 0.5d0 * l(kan)
            ks = ped(kan,2); kd = ped(kan,3)
            if (ks.gt.0) then
                do while (point(kan).gt.ks)
                   kan = point(kan)
                enddo
                l(ks) = l(ks) + r
                if (ks.ne.point(kan)) then
                     point(ks) = point(kan); point(kan) = ks
                endif
                if (kd.gt.0) then
                     do while (point(kan).gt.kd)
                        kan = point(kan)
                     enddo
                     l(kd) = l(kd) + r
                     if (kd.ne.point(kan)) then
                        point(kd) = point(kan); point(kan) = kd
                     endif
                endif
             endif
             fi = fi + l(jan) * l(jan) * d(jan)
             l(jan) = 0.0; kan = jan; jan = point(jan); point(kan) = 0
          enddo
          fff(ian) = fi
     endif
  end do ! ian
  return
  end subroutine meuwissen_code

 end subroutine inbreedng

! ===========================================
  subroutine a_inverse (nanim, jj, iun)
! ===========================================

  use filenames; use pedigreelist; use pedlisttwo; use options
  implicit none
  integer, intent(in)                  :: nanim, jj, iun
  integer, parameter                   :: ff=25
  real(8)                              :: detl, d
  integer                              :: mzhz, ll, ij, idad, imum, nrczhz,   &
&                                         irow, i1, i2, ian, ii
  integer, dimension(nanim)            :: ifirst, ivec
  integer, dimension(:), allocatable   :: inext, ivcol
  real(8), dimension(nanim)            :: rvec
  real(8), dimension(:), allocatable   :: zhz
  logical                              :: error

  if(binary)then
     open(44,file=nrmfile(jj),status='unknown',form='unformatted',action='write')
  else
     open(44,file=nrmfile(jj),status='unknown',form='formatted',action='write')
  end if
      
! set up spaces
  mzhz=nanim*ff
100  allocate(inext(0:mzhz), ivcol(0:mzhz), stat=ii)
  if( ii/=0 )then
      write(*,*)'cannot allocate arrays for NRM inverse'; stop 
  end if
  nrczhz = 0; ifirst = 0; ivcol(0) = 0; inext(0) = 0
! ... count how many elements we are going to generate 
  call nrczhz_count (error)
  deallocate(inext, ivcol,stat=ii); if(ii/=0)stop 'deall nrm'
  if(error)then
     mzhz=mzhz*1.5; go to 100 ! try again
  end if
  allocate(zhz(0:nrczhz),inext(0:nrczhz),ivcol(0:nrczhz),stat=ii)
  if(ii/=0)then
     write(*,*)'cannot allocate arrays for NRM inverse'
     write(*,*)'NRM inverse has',nrczhz,' non-zero elements'
     stop 
  end if
       
! set up NRM inverse
  detl = 0.d0; nrczhz = 0; ifirst = 0; inext(0) = 0; ivcol(0) = 0
  do ian = 1,nanim
      if(mod(ian,20000)==0)write(*,'(a,i10,a,i10)')'  NRM : at',ian,'  of',nanim
      idad = jdwec(ian,jj); imum = kdwec(ian,jj)
      if(idad>0 .or. imum>0 )then
         d = inbreed(ian,jj); detl=detl+dlog(d); d=1.d0/d
         i1 = min0(idad,imum); i2 = max0(idad,imum); ll = 0
         if( i1.gt.0 )then
            ll=ll+1; ivec(ll)=i1; rvec(ll)=0.25d0*d
            call lnkrow(i1,ll)
         end if
         if( i2.gt.0 )then
            ll=ll+1; ivec(ll)=i2; rvec(ll)=0.25d0*d
            call lnkrow(i2,ll)
         end if
         rvec(:ll)=-0.5d0*d; ll=ll+1; rvec(ll)=d; ivec(ll)=ian
         call lnkrow(ian,ll)
      else
         ivec(1)=ian; rvec(1)=1.d0
         call lnkrow(ian,1)
      end if
  end do

! write out NRM inverse
  do irow = 1, nanim
     ij = ifirst(irow)
     do while( ij>0 )
         if(binary)then
            write(44)irow,ivcol(ij),zhz(ij)
         else
            write(44,*)irow,ivcol(ij),zhz(ij)
         end if
         ij = inext(ij)
     end do
  end do ! irow

  close(44)
  write(*,*)' '
  write(*,*)'NRM inverse written out for random effect',jj
  write(*,*)'... no. of non-zero elements =',nrczhz
  write(*,*)' '
  write(iun,'(/a,i4)')'NRM inverse written out for random effect',jj
  write(iun,'(a,i12)')  '... no. of non-zero elements =',nrczhz
  write(iun,'(a,g16.8/)')'log determinant              =',detl
  if(karin) write(11,*)detl,nrczhz
  deallocate(ivcol, inext, zhz,stat=ii);  if(ii/=0)stop 'deall ainv'
  return

  contains

! ===========================
  subroutine lnkrow (irow,nn)
! ===========================

  implicit none
  integer, intent(in) :: irow, nn
  integer             :: kk, iplace, ipre, icol

  if(nn < 1) return
  ipre=0; iplace=ifirst(irow)
  do  kk = 1, nn
      icol=ivec(kk)
      do while( (iplace>0) .and. (ivcol(iplace).lt.icol) )
         ipre=iplace;  iplace=inext(iplace)
      end do
      if( (iplace > 0) .and. (icol == ivcol(iplace)) )then
         zhz(iplace)=zhz(iplace)+rvec(kk)
         ipre=iplace; iplace=inext(iplace)
      else
         nrczhz=nrczhz+1; zhz(nrczhz)=rvec(kk);  ivcol(nrczhz)=icol
         if(iplace.eq.0.and.ifirst(irow).eq.0)then
            ifirst(irow)=nrczhz; inext(nrczhz)=0
         else if(iplace.eq.0.and.ifirst(irow).ne.0)then
             inext(ipre)=nrczhz; inext(nrczhz)=0
         else if(iplace.ne.0.and.ipre.eq.0)then
             inext(nrczhz)=ifirst(irow); ifirst(irow)=nrczhz
         else
             inext(nrczhz)=inext(ipre); inext(ipre)=nrczhz
         end if
         ipre=nrczhz; iplace=inext(nrczhz)
      end if
  end do

  return
  end subroutine lnkrow

! =======================================
  subroutine pre_lnkrow (irow,nn, error)
! =======================================

  implicit none
  logical, intent(inout) :: error
  integer, intent(in)    :: irow, nn
  integer                :: kk, iplace, ipre, icol

  if( nn < 1 ) return
  ipre = 0;  iplace = ifirst(irow)
  do kk = 1, nn
     icol = ivec(kk)
      do while( (iplace>0) .and. (ivcol(iplace).lt.icol) )
         ipre = iplace; iplace = inext(iplace)
      end do
      if( (iplace.gt.0) .and. (icol.eq.ivcol(iplace)) )then
         ipre = iplace; iplace = inext(iplace)
      else
         nrczhz=nrczhz+1
         if(nrczhz.gt.mzhz)then
            print *,'"lnkrow" : dimensions exceeded !!'
            error=.true.; return
         end if
         ivcol(nrczhz) = icol
         if(iplace.eq.0.and.ifirst(irow).eq.0)then
            ifirst(irow) = nrczhz; inext(nrczhz) = 0
         else if(iplace.eq.0.and.ifirst(irow).ne.0)then
             inext(ipre) = nrczhz; inext(nrczhz) = 0
         else if(iplace.ne.0.and.ipre.eq.0)then
             inext(nrczhz) = ifirst(irow); ifirst(irow) = nrczhz
         else
             inext(nrczhz) = inext(ipre); inext(ipre) = nrczhz
         end if
         ipre = nrczhz; iplace = inext(nrczhz)
      end if
  end do

  return
  end subroutine pre_lnkrow

! ================================
  subroutine nrczhz_count (error)
! ================================

  logical, intent(out) :: error

  error=.false.
  do ian = 1,nanim
     idad = jdwec(ian,jj); imum = kdwec(ian,jj)
     i1=min0(idad,imum); i2=max0(idad,imum); ll=0
     if(i1.gt.0)then
         ll=ll+1; ivec(ll)=i1
     end if
     if(i2.gt.0)then
         call pre_lnkrow(i2,ll,error); if(error)return
         ll=ll+1; ivec(ll)=i2
     end if
     call pre_lnkrow(ian,ll,error); if(error)return
  end do
  nrczhz=nrczhz+nanim ! diagonal elements

  return
  end subroutine nrczhz_count

  end subroutine a_inverse

! ==================================
  subroutine addcodes (mm, nrl)
! ==================================

  use filenames; use pedlisttwo; use positions
  implicit none
  integer, intent(in)    :: nrl, mm
  character(len=nrl)     :: aa
  integer                :: ii, i, kk, k1, k2, ian
  logical                :: lopen
  integer, dimension(2)  :: ivec
  real(8), dimension(mm) :: dvec

 open(9,file=datnewfile,status='unknown',form='formatted',action='write')
 inquire(file = datfile, opened = lopen)
 if( .not. lopen) then
     open(1,file=datfile,status='old',form='formatted',action='read')
 else
     rewind(1)
 end if
 do
    read(1,'(a)',iostat=ii) aa; if(ii/=0)exit
    read(aa,*) dvec
    do i = 1,ngeneff
       ian = int(dvec(ipos(i)))
       k1 = 1; k2 = nannew(i)+1
       do 
          kk = k1+(k2-k1)/2                                                    
          if( ian < idwec(kk,i))then
              k2 = kk                    
          else if( ian == idwec(kk,i))then
              ivec(i) = kk; exit
          else if(ian > idwec(kk,i))then
              if( k1 == kk)stop 'err recode'
              k1 = kk
          end if
       end do
    end do
    write(9,'(a,2i8)') aa,ivec(:ngeneff)
  end do
  close(9);  close(1)

  return
  end subroutine addcodes