! Program geno2seq.f90              Convert genotypes.raw to sequences.readdepth
!                                   Simulate read depth (Poisson) for each allele
!
!                                   Written Nov 2013 by Chuanyu Sun, NAAB
!                                                 and Paul VanRaden, AIPL

	implicit none
	integer random_Poisson
	
	character*50 markername
	integer nchr,loci,totloci,nloci,id,chip,maxloc
	integer chipn, reduce1, offset1, reduce2, offset2
	integer seed0(1),ngen,location,chips,chip2
	integer,allocatable :: chrloci(:), chiplist(:,:)
	real,   allocatable :: chiperr(:), chipdepth(:)
	character*1, allocatable :: genorow(:),seqrow(:)
	integer code(36)
	
	integer nA1,nA2
    real(kind=8) uniform,gentest
	real read_depth,reduction
    integer(kind=4),allocatable :: rpoi(:)
    integer :: seed,npoi
	integer rnum,rsum
	integer i,j,k,pos,tmp1,tmp2,tmp3
    logical first/.TRUE./
	
! -----------------------------------------------  Input data files
    open ( 9,file='geno2seq.options')
    open (10,file='simseed')
    open (20,file='chromosome.data')
    open (21,file='chips.txt')
    open (30,file='genotypes.raw')
! ----------------------------------------------- Output data files
    open (60,file='sequences.readdepth')
! ----------------------------------------------- Read options file
	read(9,*)  chip2
	if (chip2 > 0) print *,'Chip2 animals assumed to have read depth 35'
	if (chip2 == 2) print *,'Sequenced animals also have chip2 with read depth 35'
	print *
! ----------------------------------------------- Read seed file
	CALL RANDOM_SEED(size=k)
	read(10,*) 
	read(10,*)  seed0(1)
	print *,'Simulation seed = ',seed0(1)
! ----------------------------------------------- Initialise values
	call random_seed (put = seed0)

! ----------------------------------------------- Create code number
	do i=1,10
		code(i)=i-1+48
	end do
	do i=65,90
		code(i-54)=i
	end do
  print '(a,36a1)','Read depth codes = ',(code(i),i=1,36)
  print *
! ----------------------------------------------  Read chromosome.data
	read(20,*) 
	do
		read(20,*,end=5) markername, nchr, loci, totloci, location, chips
    end do
	5 continue
	print *, nchr, ' chromosomes read from chromosome.data'
	print *, totloci, ' total loci read from chromosome.data'
	print *, chips, ' chips defined from chromosome.data'
	rewind(20)
		
	allocate(chrloci(nchr), genorow(totloci), chiplist(totloci,chips))
	allocate(chipdepth(chips), chiperr(chips))
	read(20,*) 
	do i=1,totloci
		read(20,*)       markername,nchr,loci, nloci, location, chips, &
		                (chiplist(i,chip), chip=1,chips)
		chrloci(nchr)=loci
    end do
	print *,       '   chromosome numbers and total markers',nchr,totloci
	print *
! ---------------------------------------------------  Read chips.txt
	print *
	print *,'   Chip  Read_Depth  Error_Rate'
	read(21,*)
	do
		read(21,*,end=21) chipn, reduce1, offset1, reduce2, offset2, &
                   chipdepth(chipn), chiperr(chipn)
		if (chipdepth(chipn) > 35) chipdepth(chipn) = 35
		print '(i7,f10.2,f10.6)', chipn, chipdepth(chipn), chiperr(chipn)
		end do
	21 continue
	print *, chipn, ' chip depth and error rates read from chips.txt'
	ngen = 0
	npoi=totloci
	!maxloc=2*MAXVAL(chrloci)
	maxloc=2*totloci
	allocate(seqrow(maxloc))
	allocate(rpoi(npoi))
	print *,'begin to read genotypes.raw'
! ----------------------------------------------  Read genotypes.raw
	do
		read(30,'(3i10,1x,2000000A1)',end=6)   id,chip,loci,(genorow(i),i=1,loci)
		ngen = ngen + 1
		read_depth = chipdepth(chip)
!                                                 Call poisson
		rsum=0
		do i=1, npoi
			rpoi(i) = random_Poisson(read_depth,first)
			rsum=rsum+rpoi(i)
			!print *,i,npoi,read_depth,rpoi(i)
		end do
		!print *,ngen,read_depth,npoi,rsum,1.0*rsum/npoi
		
		tmp1=0
		do i=1,npoi !limit to 35
!                                  Optional: Limit max read depth = 35
			if(rpoi(i) > 35)  rpoi(i)=35
!                                  Full read depth for less dense chips
			if(loci < totloci .and. chip2 > 0)  rpoi(i)=35
!
!                                  Full read depth for chip2 + sequence
			if(loci == totloci .and. chiplist(i,2) > 0 .and. chip2 == 2)  rpoi(i)=35
			if(rpoi(i)==0) tmp1=tmp1+1
		end do
		! print *, 'animalID:: the number of unread locus:: ',id,'  ',tmp1
		
		tmp1=0
		do i=1,loci
!                                      Missing genotypes have 0 read depth
			if(genorow(i)=='5') rpoi(i) = 0
			if(rpoi(i)==0) then
				write(seqrow(2*i-1),'(A1)') '0'
				write(seqrow(2*i),'(A1)')  '0'
			else
!                                            Generate sequence read errors
				if(genorow(i)=='0') then
					gentest = chiperr(chip)
					if(loci < totloci .and. chip2 > 0) gentest = 0.
					if(chip == 1 .and. chiplist(i,2) > 0 .and. chip2 == 2) gentest = 0.
				end if
				if(genorow(i)=='2') then
					gentest = 1. - chiperr(chip)
					if(loci < totloci .and. chip2 > 0) gentest = 1.
					if(chip == 1 .and. chiplist(i,2) > 0 .and. chip2 == 2) gentest = 1.
				end if
!                                                            Heterozygotes
				if(genorow(i)=='1') gentest = .5
				nA1=0
				nA2=0
!                                                  Simulate sequence reads
				do j=1,rpoi(i)
					call random_number(uniform)
					if(uniform < gentest) then
						nA1=nA1+1
					else
						nA2=nA2+1
					end if
				end do
!                                  Limit to maximum hex code = 35
            if (nA1 > 35 .or. nA2 > 35) then
               reduction = 35. / max(nA1,nA2)
!              print *,'Before = ',nA1,nA2,reduction
               nA1 = nint(nA1 * reduction)
               nA2 = nint(nA2 * reduction)
!              print *,'After  = ',nA1,nA2
               endif
				write(seqrow(2*i-1),'(A1)') code(nA1+1)
				write(seqrow(2*i),'(A1)')  code(nA2+1)
			end if
		end do
		!------------------------      output
		write(60,'(3i10,1x,2000000A1)') id, chip, loci, (seqrow(j),j=1,loci*2)
    end do ! geno file

	6 print *,ngen,' records read from genotypes.raw'
  print *
  print *,'Program geno2seq finished normally'
  stop
  end
  
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!MODULE random
! A module for random number generation from the following distributions:
!
!     Distribution                    Function/subroutine name
!
!     Normal (Gaussian)               random_normal
!     Exponential                     random_exponential
!     Poisson                         random_Poisson


INTEGER FUNCTION random_Poisson(mu, first)
!**********************************************************************
!     Translated to Fortran 90 by Alan Miller from:
!                                    RANLIB
!
!           Library of Fortran Routines for Random Number Generation
!
!                           Compiled and Written by:
!
!                                Barry W. Brown
!                                 James Lovato
!
!                    Department of Biomathematics, Box 237
!                    The University of Texas, M.D. Anderson Cancer Center
!                    1515 Holcombe Boulevard
!                    Houston, TX      77030
!
! This work was supported by grant CA-16672 from the National Cancer Institute.

!                    GENerate POIsson random deviate

!                              Function

!     Generates a single random deviate from a Poisson
!     distribution with mean mu.

!                              Arguments

!     mu --> The mean of the Poisson distribution from which
!            a random deviate is to be generated.
!                              REAL mu

!                              Method

!     For details see:

!               Ahrens, J.H. and Dieter, U.
!               Computer Generation of Poisson Deviates
!               From Modified Normal Distributions.
!               ACM Trans. Math. Software, 8, 2
!               (June 1982),163-179

IMPLICIT NONE

!     TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
!     COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL

!     SEPARATION OF CASES A AND B

!     .. Scalar Arguments ..
REAL, INTENT(IN)    :: mu
LOGICAL, INTENT(IN) :: first
!     ..
!     .. Local Scalars ..
REAL a0, a1, a2, a3, a4, a5, a6, a7, b1, b2, c, c0, c1, c2, c3, d, del,      &
     difmuk, e, fk, fx, fy, g, omega, p, p0, px, py, q, s, t, u, v, x, xx
INTEGER j, k, kflag, l, m
real    random_exponential, random_normal
LOGICAL :: full_init
!     ..
!     .. Local Arrays ..
REAL fact(10), pp(35)

REAL    :: zero = 0.0, half = 0.5, one = 1.0, two = 2.0,   &
           vsmall = TINY(1.0), vlarge = HUGE(1.0)
!       :: integral
!     ..
!     .. Intrinsic Functions ..
INTRINSIC ABS, LOG, EXP, FLOAT, IFIX, MAX, MIN, SIGN, SQRT
!     ..
!     .. Data statements ..
DATA a0, a1, a2, a3, a4, a5, a6, a7/-.5, .3333333, -.2500068, .2000118, &
                                    -.1661269, .1421878, -.1384794, .1250060/
DATA fact/1., 1., 2., 6., 24., 120., 720., 5040., 40320., 362880./

SAVE pp, s, d, l, m, p, q, p0, full_init

!     ..
!     .. Executable Statements ..
IF (mu .GT. 10.0) THEN
!     C A S E  A. (RECALCULATION OF S, D, L IF MU HAS CHANGED)

  IF (first) THEN
    s = SQRT(mu)
    d = 6.0*mu*mu

!             THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
!             PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484)
!             IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .

    l = IFIX(mu-1.1484)
    full_init = .false.
  END IF


!     STEP N. NORMAL SAMPLE - random_normal() FOR STANDARD NORMAL DEVIATE

  g = mu + s*random_normal()
  IF (g.GT.0.0) THEN
    random_Poisson = IFIX(g)

!     STEP I. IMMEDIATE ACCEPTANCE IF random_Poisson IS LARGE ENOUGH

    IF (random_Poisson.GE.l) RETURN

!     STEP S. SQUEEZE ACCEPTANCE - SAMPLE U

    fk = FLOAT(random_Poisson)
    difmuk = mu - fk
    CALL RANDOM_NUMBER(u)
    IF (d*u .GE. difmuk*difmuk*difmuk) RETURN
  END IF

!     STEP P. PREPARATIONS FOR STEPS Q AND H.
!             (RECALCULATIONS OF PARAMETERS IF NECESSARY)
!             .3989423=(2*PI)**(-.5)  .416667E-1=1./24.  .1428571=1./7.
!             THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
!             APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
!             C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.

  IF (.NOT. full_init) THEN
    omega = .3989423/s
    b1 = .4166667E-1/mu
    b2 = .3*b1*b1
    c3 = .1428571*b1*b2
    c2 = b2 - 15.*c3
    c1 = b1 - 6.*b2 + 45.*c3
    c0 = 1. - b1 + 3.*b2 - 15.*c3
    c = .1069/mu
    full_init = .true.
  END IF

  IF (g.LT.0.0) GO TO 50

!             'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)

  kflag = 0
  GO TO 70

!     STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)

  40 IF (fy-u*fy .LE. py*EXP(px-fx)) RETURN

!     STEP E. EXPONENTIAL SAMPLE - random_exponential() FOR STANDARD EXPONENTIAL
!             DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
!             (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)

  50 e = random_exponential()
  CALL RANDOM_NUMBER(u)
  u = u + u - one
  t = 1.8 + SIGN(e, u)
  IF (t.LE. (-.6744)) GO TO 50
  random_Poisson = IFIX(mu+s*t)
  fk = FLOAT(random_Poisson)
  difmuk = mu - fk

!             'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)

  kflag = 1
  GO TO 70

!     STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)

  60 IF (c*ABS(u) .GT. py*EXP(px+e) - fy*EXP(fx+e)) GO TO 50
  RETURN

!     STEP F. 'SUBROUTINE' F. CALCULATION OF PX, PY, FX, FY.
!             CASE random_Poisson .LT. 10 USES FACTORIALS FROM TABLE FACT

  70 IF (random_Poisson.GE.10) GO TO 80
  px = -mu
  py = mu**random_Poisson/fact(random_Poisson+1)
  GO TO 110

!             CASE random_Poisson .GE. 10 USES POLYNOMIAL APPROXIMATION
!             A0-A7 FOR ACCURACY WHEN ADVISABLE
!             .8333333E-1=1./12.  .3989423=(2*PI)**(-.5)

  80 del = .8333333E-1/fk
  del = del - 4.8*del*del*del
  v = difmuk/fk
  IF (ABS(v).GT.0.25) THEN
    px = fk*LOG(one + v) - difmuk - del
  ELSE
    px = fk*v*v* (((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0) - del
  END IF
  py = .3989423/SQRT(fk)
  110 x = (half - difmuk)/s
  xx = x*x
  fx = -half*xx
  fy = omega* (((c3*xx + c2)*xx + c1)*xx + c0)
  IF (kflag) 40, 40, 60

!---------------------------------------------------------------------------
!     C A S E  B.    mu < 10
!     START NEW TABLE AND CALCULATE P0 IF NECESSARY

ELSE
  IF (first) THEN
    m = MAX(1, IFIX(mu))
    l = 0
    p = EXP(-mu)
    q = p
    p0 = p
  END IF

!     STEP U. UNIFORM SAMPLE FOR INVERSION METHOD

  DO
    CALL RANDOM_NUMBER(u)
    random_Poisson = 0
    IF (u.LE.p0) RETURN

!     STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
!             PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
!             (0.458=PP(9) FOR MU=10)

    IF (l.EQ.0) GO TO 150
    j = 1
    IF (u.GT.0.458) j = MIN(l, m)
    DO k = j, l
      IF (u.LE.pp(k)) GO TO 180
    END DO
    IF (l.EQ.35) CYCLE

!     STEP C. CREATION OF NEW POISSON PROBABILITIES P
!             AND THEIR CUMULATIVES Q=PP(K)

    150 l = l + 1
    DO k = l, 35
      p = p*mu/FLOAT(k)
      q = q + p
      pp(k) = q
      IF (u.LE.q) GO TO 170
    END DO
    l = 35
  END DO

  170 l = k
  180 random_Poisson = k
  RETURN
END IF

END FUNCTION random_Poisson
! ----------------------------------------------------------------

REAL FUNCTION random_exponential( )

! Adapted from Fortran 77 code from the book:
!     Dagpunar, J. 'Principles of random variate generation'
!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM
! A NEGATIVE EXPONENTIAL DlSTRIBUTION WlTH DENSITY PROPORTIONAL
! TO EXP(-random_exponential), USING INVERSION.

IMPLICIT NONE

!     Local variable
REAL  :: r
REAL  :: zero = 0.0

DO
  CALL RANDOM_NUMBER(r)
  IF (r .GT. zero) EXIT
END DO

random_exponential = -LOG(r)
RETURN

END FUNCTION random_exponential
! -----------------------------------------------------------------------
REAL FUNCTION random_normal()

! Function from NAG
! The module uses the intrinsic RANDOM_NUMBER function supplied with the
! compiler.   It has only been tested using the Lahey LF90 compiler.

! Alan Miller

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

! 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
! --------------------------------------------------------------------
