! markersim.f90:       Simulates markers and QTLs for genomic studies
!
!                      Fortran subroutine written by Paul VanRaden 2007
!                      Animal Improvement Programs Lab, USDA
!                      Paul.VanRaden@ars.usda.gov
!
! OPTION DEFINITIONS:         (reads from markersim.options file)
!                loci      =   total markers plus QTL loci
!                qtls      =   loci with effects (or set qtls=loci)
!                curvparm  =   parameter for heavy-tailed QTL distribution
!                seednum   =   starting value for simulation (integer)
!                chromes   =   number of chromosome pairs
!                minfreq   =   minimum marker allele frequency
!                maxfreq   =   maximum marker allele frequency
!                QTLfreq   =   0 gives uniform density to QTL frequencies
!                          =   1 uses sqrt to give quadratic density
!                          =   2 uses 4th root to give U-shaped density
!                Morgans   =   length of chromosomes
  USE IFPORT
  implicit none
  real*8, allocatable :: trans(:,:),Gmtx(:,:),Gcorr(:)
  real,   allocatable :: effect(:,:)
  real,   allocatable :: freq(:),bigqtl(:),effnorm(:)
  integer,allocatable :: chro(:,:),chiplist(:,:)
  integer,allocatable :: locate(:),location(:),reduce(:,:),offset(:,:)
  integer,allocatable :: marker(:),maxqtl(:),snps(:)
  integer,allocatable :: chromnum(:),nloci(:),genotype(:)
  real*8, allocatable :: qtlvar(:)
  integer loci,qtls,maxloc,sd,ind
  integer seednum,chromes,prevchrome,seed1(1)
  integer correlated
  integer prevmin,locus,maxloci/0/,chips,chip,QTLfreq/0/
  integer mloci,loc,next,i,j,k,n,idnum,trt,maxtrt,traits,ntcp
  integer storage
  parameter (storage=10000000)
  real    minfreq,maxfreq,thisqtl,curvparm,freqblock/2000./,Morgans/1.0/
  real*4  random_normal,uniform,corruni,var
  real*8  sum,sumvar,MAF
  integer*8    kbase/10/
  character*15 markername
  character*3  oldyoung
  character*1  sex
  print *,'Output from markersim.f90'
! -----------------------------------------------  Read option file
  open ( 9,file='markersim.options')
  read(9,*)                loci,qtls,curvparm,correlated,traits
  print *,           '     loci,qtls,curvparm,correlated,traits'
  print '(i9,i7,f7.3,2i7)',loci,qtls,curvparm,correlated,traits
  read(9,*)
  read(9,*)
  read(9,*)                seednum,chromes,minfreq,maxfreq,QTLfreq
  print *,             '   seednum,chromes,minfreq,maxfreq,QTLfreq'
  print '(i9,i7,2f9.3,i5)',seednum,chromes,minfreq,maxfreq,QTLfreq
!
! -----------------------------------------------  Input data files
  if (correlated == 1) open (10,file='genetic.cor')
                       open (12,file='chips.txt')
! ............................................... Output data files
                       open (15,file='true.effects')
                       open (16,file='true.frequency')
                       open (18,file='chromosome.data')
                       open (25,file='simseed')
! ----------------------------------------------------------------
  write (25,*) ' seednum,    loci,    qtls,   traits,  correlated'
  write (25,'(5i9)') seednum, loci, qtls, traits, correlated
  maxloc = 3 * loci / chromes
  print *
  print *,freqblock,' length of blocks with common frequency'
  print '(f10.7,a)', 1. - loci/(loci + chromes*Morgans*freqblock), &
         ' prob of new frequency block between SNPs'
  read (12,*)
! -------------------------------------------  Read number of chips
  do
    read (12,*,end=5) chips
    end do
 5 continue
  print *
  print *,chips,' chip definitions read from chips.txt'
  rewind 12
!
  allocate (locate(chromes*storage),location(loci))
  allocate (freq(loci),marker(loci),genotype(loci))
  allocate (chromnum(0:loci),nloci(chromes))
  allocate (chro(maxloc,chromes))
  allocate (snps(chips),chiplist(loci,chips))
  allocate (reduce(chips,2),offset(chips,2))
! ------------------------------------------  Read chip definitions
  read (12,*,end=990)
  do
    read (12,*,end=10) chip, reduce(chip,1), offset(chip,1), &
                             reduce(chip,2), offset(chip,2)
    end do
 10 continue
  maxtrt = min(traits,10)
  ntcp = traits*(traits+1)/2
  allocate (effect(traits,loci),qtlvar(traits))
  allocate (bigqtl(traits),maxqtl(traits))
  allocate (effnorm(traits),Gcorr(ntcp))
  allocate (trans(traits,traits),Gmtx(traits,traits))
  effect = 0.
  if (curvparm <= 1.) then
    print *,'Normal distribution used for QTL effects'
    curvparm = 1.0
  else
    print *,'Heavy-tailed QTL distribution used'
    print *,' Normal SD   Heavy SD'
    do sd = 1,7
      print '(i7,f9.3)',sd,sd*curvparm**(sd - 2.)
      end do
    end if
  if (correlated == 1) then
    read(10,*,end=999)
    read(10,*,end=999)
    end if
  Gmtx = 0.
  ind=0
! ------------------------------------   Read genetic correlations
  do i=1,traits
    Gmtx(i,i) = 1.0
    if (correlated == 1) &
      read(10,*,end=999) (Gmtx(i,j),j=i,traits)
    do j=i,traits
      Gmtx(j,i) = Gmtx(i,j)
      ind=ind + 1
      Gcorr(ind)=Gmtx(i,j)
      end do ! j
!   print *, (Gmtx(i,j),j=1,traits)
    end do ! i
  if (correlated == 1) then
    print *
    print *,'Trait transformation matrix'
    end if
!                          CHOLESKY DECOMPOSITION Gcorr=Trans*Trans'
  call triang(Gcorr,NTCP,traits)
  ind=0
  do j=1,traits
    do i=j,traits
      ind=ind + 1
      trans(j,i)=0.D0
      trans(i,j)=Gcorr(ind)
      end do ! i
    if (correlated == 1) &
      print '(99f8.4)',(trans(j,i),i=1,traits)
    end do ! j
  if (correlated == 1) then
    print *
    print *,'Trait correlation matrix'
    end if
!                                      MULTIPLY TO GET  Trans*Trans'
  do j=1,traits
    do i=1,traits
      sum=0.d0
      do k=1,traits
        sum=sum + trans(i,k)*trans(j,k)
        end do ! k
      Gmtx(i,j)=sum
      end do ! i
    end do ! j
  do i=1,traits
    if (correlated == 1) &
      print '(99f8.4)',(Gmtx(i,j),j=1,traits)
    end do ! i
  print *
!                                                          Set seed
! if (seednum > 0) uniform = rand(seednum)
  seed1(1) = seednum
  if (seednum > 0) call random_seed (put = seed1)
  locate = 0
! ------------------------------  Generate marker and QTL locations
  do locus = 1,loci
!                                       Random: find empty location
    do i = 1,1000
!     uniform = rand(0)
!     call random_number(uniform)
      next = drand(0)*chromes*storage + 1
      if (locate(next) == 0) then
        locate(next) = locus
        go to 3
        end if
      end do ! i
    print *,'QTL location table too full!'
    stop
 3  continue
    end do ! locus
! ----------------------------------------   Condense QTL locations
  next = 0
  do i = 1,chromes*storage
    if (locate(i) > 0) then
      next = next + 1
      location(next) = i
      end if
    end do ! i
  if (next .ne. loci) print *,'loci mismatch ',next,loci
  qtlvar = 0.d0
  bigqtl = 0.
  print *,'  marker, chrome,  position,  freq,    effects on traits'
  prevmin = 0
  chromnum(0) = 0
  sumvar = 0.d0
  MAF = 0.d0
  call random_number(corruni)
! -------------------------------  Renumber loci within chromosomes
  do locus = 1,loci
    chromnum(locus) = (location(locus) - 1)/storage + 1
    if (chromnum(locus) .ne. chromnum(locus - 1)) loc = 0
    loc = loc + 1
    nloci(chromnum(locus)) = loc
    maxloci = max(maxloci,loc)
    if (maxloci > maxloc) then
      print *,maxloc,' maxloc too small, increase !'
      stop
      end if
    chro(loc,chromnum(locus)) = locus
!   ------------------------------   Allele frequencies and effects
!   uniform = rand(0)
!                                  Frequencies for haplotype blocks
    call random_number(uniform)
    if (uniform > loci / (loci + chromes*Morgans*freqblock)) then
      call random_number(corruni)
      call random_number(uniform)
!     print *,'locus, new value of corruni = ', locus, corruni
    else
!                                         Slightly change frequency
!                                              within 2% increments
      call random_number(uniform)
!     print *,'locus new old =', locus, corruni, &
!       (int(corruni * 50) + uniform) / 50.
      corruni = (int(corruni * 50) + uniform) / 50.
      end if
!                                                              QTLs
    if (locate(location(locus)) <= qtls) then
      freq(locus) = uniform
!                                       Make frequencies correlated
      call random_number(uniform)
      if (uniform > .1) freq(locus) = corruni
!                                 Convert uniform to quadratic freq
!                                 2nd pass converts to U-shape freq
      do j = 1, QTLfreq
        if (freq(locus) > .5) then
          freq(locus) = .5 + .5 * sqrt(abs(freq(locus) - .5)) / sqrt(.5)
        else
          freq(locus) = .5 - .5 * sqrt(abs(freq(locus) - .5)) / sqrt(.5)
          end if
        end do
!     print *,'new old freq = ', freq(locus), uniform
      sumvar = sumvar + 2. * freq(locus) * (1. - freq(locus))
      MAF = MAF + .5 - abs(freq(locus) - .5)
    else
!                                                           Markers
      freq(locus) = minfreq + uniform * (maxfreq - minfreq)
!                                       Make frequencies correlated
      call random_number(uniform)
      if (uniform > .1) then
        freq(locus) = minfreq + corruni * (maxfreq - minfreq)
        end if
!                                                     Sequence data  !!! Should put QTLs in markers
      chip = 2
      if ((locus+offset(chip,1)) / reduce(chip,1) * reduce(chip,1) /= locus+offset(chip,1) .and. &
          (locus+offset(chip,2)) / reduce(chip,2) * reduce(chip,2) /= locus+offset(chip,2)) then
        do j = 1, QTLfreq
          if (freq(locus) > .5) then
            freq(locus) = .5 + .5 * sqrt(abs(freq(locus) - .5)) / sqrt(.5)
          else
            freq(locus) = .5 - .5 * sqrt(abs(freq(locus) - .5)) / sqrt(.5)
            end if
          end do
        end if
      end if
    end do ! locus
  print *
  print *,sumvar,' sum of QTL variances (2pq)'
  print *,MAF/qtls,' average MAF of QTLs'
  print *
! -------------------------------------   Genetic variance from QTL
  do locus = 1,loci
    if (locate(location(locus)) <= qtls) then
      var = 2. * freq(locus) * (1. - freq(locus))
!                                               Independent traits
      do trt=1,traits
        effnorm(trt) = random_normal()
        end do ! trt
!                                                Correlated traits
      do trt=1,traits
        effect(trt,locus) = 0.
        do j = 1,traits
          effect(trt,locus) = effect(trt,locus) &
                            + trans(trt,j)*effnorm(j)
          end do ! j
!                                                      Heavy tailed
        effect(trt,locus) = effect(trt,locus) * &
          curvparm**(abs(effect(trt,locus)) - 2.)
        effect(trt,locus) = effect(trt,locus) * sqrt(1./(sumvar))
        thisqtl = effect(trt,locus)**2 * var
        qtlvar(trt) = qtlvar(trt) + thisqtl
        if (thisqtl > bigqtl(trt)) then
          bigqtl(trt) = thisqtl
          maxqtl(trt) = locus
          end if
        end do ! trt
      end if
    if (locus <= 10 .or. locus == loci) &
      print '(2i9,i10,f9.5,2x,99f10.5)',locus,chromnum(locus) &
       ,location(locus),freq(locus) &
       ,(effect(trt,locus),trt=1,maxtrt)
!                                Write combined marker and QTL data
    write (15,'(f8.4,i10,99f10.5)') &
        freq(locus),location(locus) &
        ,(effect(trt,locus),trt=1,traits)
    end do ! locus
  if (loci > 10) print *,'Only first 10 QTLs listed above'
  mloci = loci
  print *
  if (qtls < loci) then
    mloci = loci - qtls
    print *,qtls,' of the ',loci,' loci had effects' 
    print *,'     QTL,  chrome,   position,   freq,' &
           ,'     effects on traits'
    next = 0
    do locus = 1,loci
      if (locate(location(locus)) <= qtls) then
        next = next + 1
        if (next <= 10) print '(2i9,i10,f9.5,2x,99f10.5)',locus &
          ,chromnum(locus),location(locus),freq(locus) &
          ,(effect(trt,locus),trt=1,maxtrt)
      else
        marker(locus - next) = locus
        end if
      end do ! locus
  else
    print *,qtls,'QTL and marker locations assumed to be identical'
    do locus = 1,loci
      marker(locus) = locus
      end do ! locus
    end if
  print *
  print *,'Trait  Genetic var   Largest QTL var    Freq     Effect   Locus '
  do trt = 1,traits
    print '(i3,f12.5,f15.5,a,2f10.5,i10)' &
      ,trt,qtlvar(trt),bigqtl(trt)*100./qtlvar(trt),' %  ' &
      ,freq(maxqtl(trt)),effect(trt,maxqtl(trt)),maxqtl(trt)
    end do
  print *
  do i = 1,chromes
    print *,'chromosome ',i,' contains ',nloci(i),' loci'
    end do
  write (18,'(a,10x,2a)') 'SNPname','chrome  within  overall  location' &
    ,'  n_chips   chip1   chip2   etc'
  prevchrome = 0
  n = 0
  snps = 0
  chiplist = 0
! -------------------------------------   Write positions of markers
  do locus = 1,mloci
    write (markername(1:6),'(a6)') 'Marker'
    write (markername(7:15),'(i9)') locus
    do i = 7,15
      if (markername(i:i) == ' ') markername(i:i) = '0'
      end do
    if (chromnum(marker(locus)) > prevchrome) then
      prevchrome = chromnum(marker(locus))
      n = 1
      end if
    do chip = 1,chips
      if ((locus+offset(chip,1)) / reduce(chip,1) * reduce(chip,1) == locus+offset(chip,1) .or. &
          (locus+offset(chip,2)) / reduce(chip,2) * reduce(chip,2) == locus+offset(chip,2)) then
        snps(chip) = snps(chip) + 1
        chiplist(locus,chip) = snps(chip)
        end if
      end do ! chip
!    if (chips > 1 .and. locus / reduce * reduce == locus) then
!      snps(2) = snps(2) + 1
!      chiplist(locus,2) = snps(2)
!      end if
!    if (chips > 2 .and. locus / reduce**2 * reduce**2 == locus) then
!      snps(3) = snps(3) + 1
!      chiplist(locus,3) = snps(3)
!      end if
!    if (chips > 3 .and. locus / reduce**3 * reduce**3 == locus) then
!      snps(4) = snps(4) + 1
!      chiplist(locus,4) = snps(4)
!      end if
    write (18,'(a,i5,2i9,i12,i5,99i9)')  &
      markername,chromnum(marker(locus)),n,locus &
     ,(location(marker(locus)) - (chromnum(marker(locus))-1)*storage)*kbase &
     ,chips,(chiplist(locus,chip),chip=1,chips)
!   ----------------------------------------------   Write true.frequency
    write (16,'(a,i10,f9.5)') markername,locus,freq(marker(locus))
    n = n + 1
    end do ! locus
  print *
  print *,'Program markersim completed normally'
  stop
 990 print *,'chips.txt file is missing'
  stop
 999 print *,'phenotype.data0 file is missing'
  stop
  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 rand(0)
! 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,   &
                             vsmall = TINY(1.0), vlarge = HUGE(1.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
!--------------------------------------------------------------------
!            DOUBLE PRECISION SUBROUTINE TO PERFORM CHOLESKY
!            DECOMPOSITION       STEVE P SMITH   JUNE 1982
!   IF A IS POSITIVE DEFINITE THERE EXISTS A UNIQUE LOWER
!   TRIANGULAR MATRIX L SUCH THAT A=LL'.  A MUST CONTAIN THE
!   UPPER TRIANGULAR ELEMENTS IN ROW ORDER.  THE LOWER TRIANGULAR
!   ELEMENTS OF L IN COLUMN ORDER WILL BE RETURNED IN THE SAME
!   SPACE (A IS OVERWRITTEN BY L, BOTH HALF-STORED)
!
      SUBROUTINE TRIANG(A,LENGTH,M)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      DIMENSION A(LENGTH)
      N=1
      K=M
      DO 25 L=1,M
        X=DSQRT(A(N))
        A(N)=X
        IF (L.EQ.M) GO TO 30
        N=N+1
        DO 10 I=N,K
10        A(I)=A(I)/X
        KK=K
        DO 20 J=N,K
          X=A(J)
          IF(X .EQ. 0.D0) GO TO 20
          IDIF=KK - J + 1
          DO 15 I=J,K
15          A(I+IDIF)=A(I+IDIF) - X*A(I)
20        KK=KK + K - J + 1
        N=K+1
25      K=M-L+K
30    RETURN
      END
