program MYOWNDGD

       IMPLICIT NONE
       INTEGER MXGAM, DIMX
       PARAMETER (MXGAM=60) !option du type de matrice
       PARAMETER (DIMX=10000) !taille max de la matrice sous forme vecteur
       DOUBLE PRECISION X(DIMX),GAMMA(MXGAM),DET !X matrice var-cov, type matrice, determinant
       CHARACTER *72 BASENAME ! nom fichier
       INTEGER NV,NGAM,OPT,IFRM !NGAM=nb parametres NV=nb niveau var matrice OPT=option sortie .own et .gdg

       INTEGER I,J,K,LL,LENBN
       DOUBLE PRECISION C0,C1,C2

!modif (david,2013) pour faire passer modele anedependant
	Integer :: oante,ovar,NGAMMAANTE,col,rtest,jd,datatemps	
	Integer, allocatable ::  temps(:) !vecteur de stockage des valeurs de temps
	Integer, allocatable :: tante(:),NGAMADM(:) 
	DOUBLE PRECISION, allocatable :: GAMMAANTE(:)
	character :: tmodel*40,file*80,BASEFILE*80
	integer, allocatable :: typepara(:),typeorder(:)
integer,parameter::   io_p=40          ! unit number for parameter file




!      ASREML writes  BASENAME.own
!             runs    MYOWNGDG BASENAME
!
!             reads   BASENAME.gdg
! ------------------------------------------
!      First obtain BASENAME
! Unix & DIgital Fortran
       CALL GETARG(1,BASENAME)  !permet recup de l argument 1 ecrit apres executable pour creer les .gdg et savoir quel .own lire
! Lahey
!      CALL GETCL(BASENAME)
       
       LENBN=INDEX(BASENAME,' ')

! ------------------------------------------
!      NOW READ   BASENAME.own
      BASENAME(LENBN:)='.own'

      OPEN(21,file=BASENAME,status='old',err=920)
      READ(21,21,ERR=930,END=930) NV,NGAM,OPT !NV calcule tout seul par asreml, NGAM=chiffre apres OWN, OPT, donnee apres !F
21   FORMAT(4I5)

      IF(NGAM.GT.MXGAM) GOTO 930
 allocate(temps(NV))
 temps=0

BASEFILE='para'
if (OPT.ne.1) then
BASEFILE(5:)=str(OPT)
endif
 
!ma partie avec fonction speciale
!---------------------------------
!---------------------------------
 

call lectureparam
 print*,'lecture param OK'


if (tmodel.eq.'SAD') then
	
if (NGAMMAANTE.ne.NGAM) then
STOP 'ERROR NB PARAMETERS'
endif	

 	READ(21,22,ERR=930,END=930) (GAMMAANTE(I),I=1,NGAM)
 
 22   FORMAT(6D13.7)
      CLOSE(21)
 
!debut ecriture du .gdg

  IFRM=-1
	BASENAME(LENBN:)='.gdg'
      OPEN(21,file=BASENAME,status='unknown',err=950)
      WRITE(21,21) NV,NGAM,OPT,IFRM
 
      WRITE(21,22) (GAMMAANTE(I),I=1,NGAM)

 call makeinverseSAD
 print*,'makeinverse OK'
 call makeinversederivSAD
print*,'makederiv OK'
	close(21)

else if (tmodel.ne.'SAD') then
!partie de Gilmour
   print*,'A SAD model is not required, try another software'  
endif


!     This bit of code designed to make sure file is properly closed before returning

      OPEN(22,FILE=BASENAME,STATUS='OLD',ERR=950)
      READ(22,21,err=950) I
      CLOSE(22)

! ------------------------------------------
       STOP
910    STOP 'MYOWNGDG has no BASENAME argument on command line'
920    STOP 'MYOWNGDG could not open BASENAME.own'
930    STOP 'MYOWNGDG failed to read BASENAME.own'
940    STOP 'MYOWNGDG has invalid/inconsistent data in BASENAME.own'
950    STOP 'MYOWNGDG could not open BASENAME.gdg'
960    STOP 'MYOWNGDG has no BASENAME argument on command line'
 
 
 
       
contains



subroutine lectureparam 
 integer n, x(40),ef,stat,curr_eff,i,j,k,ii
 character xc(40)*40
integer, dimension(size(temps)) :: ttemps




open(io_p,file=BASEFILE)

 call readline(io_p,x,xc,n)
   if (xc(1) /= 'MODEL') then
      write(201,*)'MODEL not found'
      stop
  endif 
      
  call readline(io_p,x,xc,n) 
  tmodel=xc(1)
  
  
 if (tmodel.eq.'SAD') then
 
 call readline(io_p,x,xc,n)
   if (xc(1) /= 'ANTEDEP_ORDER') then
      write(201,*)'ANTEDEP_ORDER not found'
      stop
  endif 
  
      
  call readline(io_p,x,xc,n) 
  oante=x(1)
  
  
allocate(tante(oante))
  
  
   call readline(io_p,x,xc,n)
   if (xc(1) /= 'MODEL_ANTEDEP') then
       write(201,*)'MODEL_ANTEDEP not found'
       stop
   endif
   
   
	call readline(io_p,x,xc,n)
	tante(1:oante)=x(1:oante)
	
call readline(io_p,x,xc,n)
   if (xc(1) /= 'MODEL_INNOVATION_VARIANCE') then
       write(201,*)'MODEL_INNOVATION_VARIANCE not found'
       stop
   endif	
   call readline(io_p,x,xc,n)
	ovar=x(1)
	
	
   NGAMMAANTE=oante+sum(tante(1:oante))+1+ovar
   allocate(GAMMAANTE(oante+sum(tante(1:oante))+1+ovar),typepara(NGAM),typeorder(NGAM))
   GAMMAANTE=0.0d0
   typepara=0
   typeorder=0
do i=1,oante
do j=0,tante(i)
typepara(sum(tante(1:i-1))+i+j)=i
typeorder(sum(tante(1:i-1))+i+j)=j
enddo
enddo   
do i=0,ovar
typeorder(sum(tante(1:oante))+oante+i+1)=i
enddo


call readline(io_p,x,xc,n)
   if (xc(1) /= 'DATA_FILE') then
       write(201,*)'DATA_FILE not found'
       stop
   endif	
   call readline(io_p,x,xc,n)
	file=xc(1)

call readline(io_p,x,xc,n)
   if (xc(1) /= 'COLUMN_FACTOR') then
       write(201,*)'COLUMN_FACTOR not found'
       stop
   endif	
   call readline(io_p,x,xc,n)
	col=x(1)
	
	
open(12,file=file)
jd=1

do
call readline(12,x,xc,n)
  if (n == -1) exit
 
 datatemps=x(col)
 call test
 if(jd.eq.NV+1) exit
enddo   
close (io_p)

! ttemps=QuickSort(temps) 
 !  temps=ttemps
 call tri(temps) 
    endif
end subroutine lectureparam


subroutine test
integer :: k
rtest=0
do k=1,jd

if (temps(k).ne.datatemps) then
rtest=rtest+1
endif
enddo
if (rtest.eq.jd) then
temps(jd)=datatemps
jd=jd+1
endif
end subroutine test

subroutine tri(t)
implicit none
integer, dimension (:), intent(inout) :: t
integer :: temp
integer :: i,j
do i=1,size(t)-1
	do j=i+1,size(t)
	if (t(i).gt.t(j)) then
	temp=t(i)
	t(i)=t(j)
	t(j)=temp
	endif
	enddo
enddo
end subroutine tri

recursive function quicksort(InList) result(OutList)
    
    INTEGER,DIMENSION(:) :: InList
    INTEGER,DIMENSION(size(InList)) :: OutList
    INTEGER,DIMENSION(size(InList)) :: SupList, OrderedSupList, InfList, OrderedInfList
    INTEGER :: pivot
    INTEGER :: i,j, InfListSize, SupListSize

print*,'tailletemps',size(InList,1)
if(size(InList,1) < 2) then
       OutList(1) = Inlist(1)
    else
       ! Le pivot sera le premier élément de la liste
       pivot = InList(1)
 
       ! On trie la liste
       InfListSize = 0
       SupListSize = 0
       do i = 2, size(InList,1)
          if(InList(i) < Pivot) then
             InfListSize = InfListSize + 1
             InfList(InfListSize) = InList(i)
          elseif(InList(i) >= Pivot) then
             SupListSize = SupListSize + 1
             SupList(SupListSize) = InList(i)
          end if
       enddo
 print*,'trie OK'
 print*,'la liste',suplist(:)
       ! On recompose la liste
       if(InfListSize < 1) then
          OrderedSupList = QuickSort(SupList(1:SupListSize))
          OutList = (/ Pivot, (OrderedSupList(j), j=1,SupListSize) /)
       elseif(SupListSize < 1) then
          OrderedInfList = QuickSort(InfList(1:InfListSize))
          OutList = (/ (OrderedInfList(j), j=1,InfListSize), Pivot /)
       else
          OrderedInfList = QuickSort(InfList(1:InfListSize))
          OrderedSupList = QuickSort(SupList(1:SupListSize))
 
          OutList = (/ (OrderedInfList(j), j=1,InfListSize), Pivot, (OrderedSupList(j), j=1,SupListSize) /)
       endif
    end if
print*,'recomp OK'
end function quicksort

subroutine makeinverseSAD
!subroutine qui calcul G-1=L'D-1L, selon formule G-1(ij)=sum(l=1 to NV)(sum(m=1 to NV)(L(mi)*L(lj)*D-1(ml)))
integer :: i,j,k,l,m
real :: val,theta,thetap,var
do i=1,NV
do j=max(i-oante,1),i
val=0.0d0

	do l=1,NV
	theta=0.0d0
	if (l.eq.j) then
	theta=1.0d0
	elseif (l.gt.j.and.l.le.j+oante) then
	do k=0,tante(l-j)
	theta=theta-gammaante(sum(tante(1:l-j-1))+l-j+k)*temps(l)**(k)
	enddo
	endif
	thetap=0.0d0
	if (l.eq.i) then
	thetap=1.0d0
	elseif (l.gt.i.and.l.le.i+oante) then
	do k=0,tante(l-i)
	thetap=thetap-gammaante(sum(tante(1:l-i-1))+l-i+k)*temps(l)**(k)
	enddo
	endif
		var=1.0d0
		do k=0,ovar
		var=var*exp(gammaante(sum(tante(1:oante))+oante+k+1)*temps(l)**(k))
		enddo
		!print*,'variance au temps',temps(l),':',var
		!print*,'para utilise',gammaante(sum(tante(1:oante))+oante+1),gammaante(sum(tante(1:oante))+oante+2)
	val=val+(theta*thetap/var)
	enddo

		if (val.ne.0.0d0) then
		WRITE(21,'(I5,I5,G18.10)') i,j,val
		!print*,'i',i,'j',j,'val',val
		endif
enddo
enddo


end subroutine makeinverseSAD



subroutine makeinversederivSAD
!subroutine qui calcul G-1=L'D-1L, selon formule G-1(ij)=sum(l=1 to NV)(sum(m=1 to NV)(L(mi)*L(lj)*D-1(ml)))
integer :: i,j,k,l,m,ii
real :: val,theta,thetap,var

do ii=1,NGAM
 if (typepara(ii).gt.0) then

   do i=1,NV
   do j=max(i-oante,1),i
   val=0.0d0
	!dL'/dpara*G-1L, theta qui change
	do l=1,NV
	
	theta=0.0d0
	if (l-j.eq.typepara(ii)) then
	theta=-temps(l)**typeorder(ii)
	endif
	thetap=0.0d0
	if (l.eq.i) then
	thetap=1.0d0
	elseif (l.gt.i.and.l.le.i+oante) then
	do k=0,tante(l-i)
	thetap=thetap-gammaante(sum(tante(1:l-i-1))+l-i+k)*temps(l)**(k)
	enddo
	endif
		var=1.0d0
		do k=0,ovar
		var=var*exp(gammaante(sum(tante(1:oante))+oante+k+1)*temps(l)**(k))
		enddo
	val=val+(theta*thetap/var)
	enddo



	!L'*G-1dL/dpara, thetap qui change
	do l=1,NV
	
	theta=0.0d0
	if (l.eq.j) then
	theta=1.0d0
	elseif (l.gt.j.and.l.le.j+oante) then
	do k=0,tante(l-j)
	theta=theta-gammaante(sum(tante(1:l-j-1))+l-j+k)*temps(l)**(k)
	enddo
	endif
	
	thetap=0.0d0
	if(l-i.eq.typepara(ii)) then
	thetap=-temps(l)**typeorder(ii)
	endif
	
	
	
		var=1.0d0
		do k=0,ovar
		var=var*exp(gammaante(sum(tante(1:oante))+oante+k+1)*temps(l)**(k))
		enddo
		
	val=val+(theta*thetap/var)
	enddo
	
   !if (val.ne.0.0d0) then
   WRITE(21,'(I5,I5,G18.10)') i,j,val
   !print*,'i',i,'j',j,'val',val
   !endif
   enddo
   enddo


 elseif (typepara(ii).eq.0) then !que L'dG-1/dparaL
 do i=1,NV
   do j=max(i-oante,1),i
   val=0.0d0
 
	do l=1,NV
	theta=0.0d0
	if (l.eq.j) then
	theta=1.0d0
	elseif (l.gt.j.and.l.le.j+oante) then
	do k=0,tante(l-j)
	theta=theta-gammaante(sum(tante(1:l-j-1))+l-j+k)*temps(l)**(k)
	enddo
	endif
	thetap=0.0d0
	if (l.eq.i) then
	thetap=1.0d0
	elseif (l.gt.i.and.l.le.i+oante) then
	do k=0,tante(l-i)
	thetap=thetap-gammaante(sum(tante(1:l-i-1))+l-i+k)*temps(l)**(k)
	enddo
	endif
		var=1.0d0
		do k=0,ovar
		var=var*exp(-gammaante(sum(tante(1:oante))+oante+k+1)*temps(l)**(k))
		enddo
		if (typeorder(ii).eq.0) then
		var=var*(-1)
		else
		var=var*(-temps(l)**(typeorder(ii)))
		endif
	
	val=val+(theta*thetap*var)
	enddo

		!if (val.ne.0.0d0) then
		WRITE(21,'(I5,I5,G18.10)') i,j,val
		!print*,'i',i,'j',j,'val',val
		!endif
 enddo
 enddo
 endif
enddo
end subroutine makeinversederivSAD

subroutine readline(unit,x,xc,n,nohash)
!subroutine of Ignacy Misztal http://nce.ads.uga.edu/~ignacy/newprograms.html
! Reads a line from unit and decomposes it into numeric fields in x and
! alphanumeric fields in xc; the number of fields is stored in n.
! If optional variable nohash is absent, 
!  all characters after character # are ignored.

  integer unit,x(40),n,stat,i
  integer,optional::nohash
  character xc(40)*(*), a*200
 
  do
     read(unit,'(a)',iostat=stat)a
     if (stat /= 0) then
         n=-1
         return
     endif
     a=adjustl(a)               ! ignore leading spaces
     
     if (present(nohash)) then 
        exit
       else	!ignore characters after #
        i=scan(a,'#')
        if (i == 0) exit
        if (i>1) then
            a=a(1:i-1)
            exit
        endif
     endif  	
  enddo
  
  call nums2(a,n,x,xc)
  return
end subroutine


subroutine nums2(a,n,x,xc)
!subroutine of Ignacy Misztal http://nce.ads.uga.edu/~ignacy/newprograms.html
! separates array a into items delimited by blanks. character elements are
! put into optional character vector xc, decoded numeric values 
! into optional real vector x, and n contains the number of items. The 
! dimension of x and xc can be lower than n.
! A modification of nums() from f77 to f90
! Now accepts real numbers
! 2/23/2000

character (*)::a
character (*),optional::xc(:)
integer,optional::x(:)
integer::n,curr,first,last,lena,stat,i

curr=1;lena=len(a);n=0

do 
  ! search for first nonspace
  first=0
  do i=curr,lena
     if (a(i:i) /= ' ') then
        first=i
	exit
     endif
  enddo
  if (first == 0) exit
  
  
  ! search for first space
  curr=first+1
  last=0
  do i=curr,lena
      if (a(i:i) == ' ') then
        last=i
        exit
      endif
  enddo
  
  if (last == 0) last=lena
  
  n=n+1
  if (present(xc)) then
     if (size(xc) >= n) xc(n)=a(first:last)
  endif   
  if (present(x)) then
     if (size(x) >= n) then
        read(a(first:last),'(i12)',iostat=stat)x(n)
        if (stat /=0) x(n)=0
     endif
  endif   
	
  curr=last+1
enddo
end subroutine

character(len=20) function str(k)
!   "Convert an integer to string."
    integer, intent(in) :: k
    write (str, *) k
    str = adjustl(str)
end function str

end program MYOWNDGD





