! Generates (co)variances for random regression models given shapes of
! variance, correlations, and variances avergaed over lactations.
! Ignacy Misztal, 10/22/1999 - 10/26/1999

use kinds; use denseop
implicit none
integer,parameter::ntrait=2   !number of traits in G0
integer,parameter::nreg=5     !number of random regressions
integer,parameter::day1=1,daylast=305
real (8)::g0(ntrait,ntrait)   ! G averaged over lactation
real (8)::gdidj(ntrait,ntrait,305,305)   ! g0 for each combination of test day
real (8)::gr(ntrait*nreg,ntrait*nreg)    ! to be estimated
integer,parameter::neq=ntrait**2 * nreg**2
real (r8)::x(neq),xx(neq,neq),xy(neq),sol(neq)   ! least square equations
integer::i,j,k,l,m,n,t1,t2

! set up g0
g0=reshape( (/1.0, 2., 2., 4.0/),(/2,2/))

! check f and corr
print*,' INITIAL VARIANCES AND CORRELATIONS'
print*,'dim    var(dim)  corr(1,dim) corr(150,dim) corr(300,dim)'
do i=1,350,20
   print'(i3,4g13.3)',i,f(i),corr(1,i),corr(150,i),corr(300,i)
enddo   

! calculate g0(dim1,dim2)
 do i=1,305
    do j=1,305
       do k=1,ntrait
          do l=1,ntrait
             gdidj(k,l,i,j)=g0(k,l)* sqrt(f(i)*f(j)) * corr(i,j)
          enddo
       enddo    
    enddo
 enddo
 
 Print*,' calculated Gdidj'
 
 ! set up equations
 xx=0; xy=0; sol=0
 do i=1,305,5
    do j=1,305,5
       do k=1,ntrait
          do l=1,ntrait
           ! First set up x
             x=0
             do m=1,nreg
                do n=1,nreg
                   x(address(k,l,m,n))=i**(m-1) /(305.**(m-1)) * &
		                       j**(n-1) /(305.**(n-1))
                enddo  
             enddo
             ! set up xx and xy
             do m=1,neq
                do n=1,neq
                   xx(m,n)=xx(m,n)+x(m)*x(n)
                enddo
                xy(m)=xy(m)+x(m)*gdidj(k,l,i,j)
             enddo  
         enddo
       enddo    
    enddo
 enddo
 
 print*,' Equations set up'
	     
! Weights can be added above if variances of different traits are markedly
! different
! Solve
call solve_s(xx,xy,sol,rank=i)
print*,size(xy),' equations, rank=',i
        
! Create variance-covariance matrices
 do k=1,ntrait
    do l=1,ntrait
       do m=1,nreg
          do n=1,nreg
            gr(m+(k-1)*nreg,n+(l-1)*nreg)=sol(address(k,l,m,n))
          enddo  
       enddo
    enddo
 enddo    

   
call printmat(gr,'gr','(8f10.3)')
! print variances for trait 1 and correlations between day 150 and given day

! combinations of trait 1 only: 
t1=1; t2=1
print*,' ESTIMATED VARIANCES AND COVARIANCES'
print*,' DIM   var(dim) corr(1,dim)  corr(150,dim) corr(300,dim)'
do i=1,350,20
   print '(i3,4g13.3)',i,covarrr(i,i,t1,t2,gr),&
           covarrr(1,i,t1,t2,gr)/ &
	         sqrt(covarrr(1,1,t1,t2,gr)* covarrr(i,i,t1,t2,gr)),&
           covarrr(150,i,t1,t2,gr)/ &
	         sqrt(covarrr(150,150,t1,t2,gr)* covarrr(i,i,t1,t2,gr)),&
           covarrr(300,i,t1,t2,gr)/ &
	         sqrt(covarrr(300,300,t1,t2,gr)* covarrr(i,i,t1,t2,gr))
enddo

CONTAINS

function covarrr(i,j,t1,t2,gr)
real (8)::covarrr,gr(:,:),pi(size(gr,1)),pj(size(gr,1))
integer::i,j,k,t1,t2
! covarr=pi'Gr pj, where pi is coefficients of random regression at
! day i
!
pi=0; pj=0

do k=1,nreg
   pi((t1-1)*nreg+K)=I**(K-1) / (305.**(k-1))
   PJ((t2-1)*nreg+K)=J**(K-1) / (305.**(k-1))
enddo
covarrr=dot_product(matmul(pi,gr),pj)
end function

 function address(nt1,nt2,nr1,nr2)
 ! returns address of equation containing the given combination of trait
 ! and random regressions; sorted within traits
integer::address,nt1,nt2,nr1,nr2
address=(nt1+(nr1-1)*ntrait-1)*ntrait*nreg+ nt2+(nr2-1)*ntrait
end function

function f(dim) result (var)
! calculates variances over days in milk standarizing function f1
integer::dim,i
real::var
real,save::ave=0
! standarize if necessary
if (ave == 0) then 
   do i=1,305
      ave=ave+f1(i)
   enddo
   ave=ave/304
   print*,'ave=',ave!
endif
var=f1(dim)/ave
end function
      


function f1(dim) result (var)
! calculates variances over days in milk
integer::dim,i
real (8)::var
real::points(5)=(/1,50,200,250, 305/)
! Function given by points
real::values(5)=(/1., .6, .5, .8, .6/)
!Extrapolate
if (dim < 1 .or. dim > points(size(points))) then
print*,'dim too large: ',dim
endif
do i=1,size(points)-1
 if (dim >= points(i) .and. dim <= points(i+1)) then
    var=values(i)+(values(i+1)-values(i))*(dim-points(i))/(points(i+1)-points(i))
    exit
  endif
enddo
end function



function corr(dim1,dim2)
! calculates correlations between two different days in milk
! using the autocorrelative approach; the autocorrelation coeeficient is set 
! to such a number so that corr(1,305)=corr1_305
integer::dim1,dim2
real::corr,rho
real,parameter::corr1_305=.4
rho=exp(log(corr1_305)/304.)
corr=rho**(abs(dim1-dim2))
end function

end