! Multivariate Universal Maximum Likelihood Estimation (MUMLE)
! Developed by 
!                      Lianhong Gu
!                      Environmental Sciences Division
!                      Oak Ridge National Lab
!                      lianhong-gu@ornl.gov
!
! Version 1.0, July, 2006

      subroutine coapproxyfunc(nkrigpoints,nxdim,xsampkrig,
     &    ysampkrig,nvarscp,indexvarscp,iupdatelikelihood,
     &    iwhattodo,ithetap0,ixnewpos,xnew,
     &    ypred,grad_xnew,rootmeanerror)
      implicit none
      include '../maxlikelihood/colikelihood.h'

!!!!!!!!!!!!!!!!!!!! Arguments !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-------------------Inputs--------------------------------------
!nkrigpoints: the number of available points for kriging.
!nxdim: the dimension in the design parameter vector x 
!xsampkrig: the coordinates of the sampled points for kriging. Note that these
!           points may or may not be the same as those points used for estimating
!           the parameters in the likelihood function
!
!  ****************************************************************************
!  *Note different variables are stacked in the vector (xsampkrig,ysampkrig)  *
!  *in an orderly fashion. The number of different variables is given by      *
!  *nvarscp. The number of points in each variable is indicated by indexvarscp*
!  ****************************************************************************
!
!ysampkrig: the value at xsampkrig
!nvarscp:   the number of variables for cokriging
!indexvarscp: the number of samples in each variable for cokriging
!iupdatelikelihood: whether to update the likelihood function
!           =0: the first entering, so initilize the parameters to be optimized
!               and do likelihood function optimization
!           =1: update the likelihood parameters using the saved values from
!               the last optimization as the initial guesses
!           =2: do not do likelihood optimization
!iwhattodo: =1, only make a predition at xnew
!           =2, make a prediction and calculate the derivative at xnew
!           =3, only calculate the derivative at xnew
!ithetap0: =ithetap: =1 set pdist to 2.0, =2: pdist to be optimized   
!ixnewpos: which variable does xnew belong to.
!xnew:      the coordinates of the point to be kriged
!
!---------------------- Outputs --------------------------------
!ypred: the predicted value at xnew
!rootmeanerror: the root mean square error of ypred
!grad_xnew: the derivative at xnew
      
      integer nkrigpoints,nxdim,iupdatelikelihood,iwhattodo,
     &    nvarscp,indexvarscp(nvarscp),ithetap0,ixnewpos
      double precision xsampkrig(nkrigpoints,nxdim),
     &     ysampkrig(nkrigpoints),xnew(nxdim),ypred,
     &     rootmeanerror,grad_xnew(nxdim)
!
!!!!!!!!!!!!!!!!!!! Local variables !!!!!!!!!!!!!!!!!!!!!!!!!
!
      double precision ftol
      parameter(ftol=1.0d-8)
      external df1dimd_like,gradlikefunc,
     &         likelihoodfunc,f1dim_like,colikelihoodgrad
      integer i,j,nplike,info,nii,index,irowinc,istart,
     &        icompactpostri,icompactposoff
      integer nbd(maxnvars*(maxnvars+1)*
     &            (2*maxdesignparam+1)/2)
      double precision plike(maxnvars*(maxnvars+1)*
     &            (2*maxdesignparam+1)/2),
     &    x(nxdim),flikemin,xsamp(nxdim),
     &    theta((maxnvars*(maxnvars+1))/2,maxdesignparam),
     &    pdist((maxnvars*(maxnvars+1))/2,maxdesignparam),
     &    theta1(nxdim),pdist1(nxdim),
     &    covars(maxnvars),correls(maxnvars*(maxnvars-1)/2),
     &    rx(nkrigpoints),fb(maxkbasvar),corr_like,
     &    xmean(maxdesignparam),
     &    gradient(maxkbasvar,nxdim)
     
      double precision gradlike(maxnvars*(maxnvars+1)*
     &                        (2*maxdesignparam+1)/2)
      save nplike
      save plike,theta,pdist,xmean,covars,correls

      if(iupdatelikelihood.le.1)then
! update parameters and correlation functions in the likelihood function
!       
! nbd is an INTEGER array of dimension n that must be set by the
! user to the type of bounds imposed on the variables:
! nbd(i)=0 if x(i) is unbounded,
!        1 if x(i) has only a lower bound,
!        2 if x(i) has both lower and upper bounds, 
!        3 if x(i) has only an upper bound.              

        nlikepoints=nkrigpoints
        nvars=nvarscp
        ithetap=ithetap0
        do i=1,nvars
          indexvars(i)=indexvarscp(i)
        enddo
        ndesignparam=nxdim
        do i=1,nlikepoints
          ysamplike(i)=ysampkrig(i)       
        enddo
      
! Centralize the coordinates to stablize the linear system
        do j=1,nxdim
          xmean(j)=0.0d0
          do i=1,nlikepoints
            xmean(j)=xmean(j)+xsampkrig(i,j)
          enddo
          xmean(j)=xmean(j)/dble(nlikepoints)
          do i=1,nlikepoints
            xsamplike(i,j)=xsampkrig(i,j)-xmean(j)
          enddo
        enddo
        
        call setinit_bounds(nvars,ndesignparam,ithetap,
     &  iupdatelikelihood,plikemin,plike,plikemax,nbd,nplike)

        isitlastcall=.false.

! to use Lbfgsb_2_4, no need to call the cost function first
        call Lbfgsb_2_4(nplike,plike,flikemin,plikemin,plikemax,
     &      nbd,colikelihoodgrad,ftol,info)
                
! to use frprmn, the cost function needs to be called first
!        call likelihoodfunc(nplike,plike,flikemin)
!        call frprmn(plike,nplike,ftol,flikemin,plikemin,
!     &       plikemax,gradlikefunc,f1dim_like,df1dimd_like)
     
!        call dfpmin(plike,nplike,ftol,flikemin,likelihoodfunc,
!     &     plikemin,plikemax,gradlikefunc,f1dim_like,df1dimd_like)

        write(*,*)flikemin,plike(1),plike(2),plike(3),info
        call gradlikefunc(nplike,plike,flikemin,gradlike)
        write(*,*)gradlike(1),gradlike(2),gradlike(3)
       
        call nongradopt(nplike,likelihoodfunc,f1dim_like,plike,
     &      plikemin,plikemax,ftol,flikemin)
        
        write(*,*)flikemin,plike(1),plike(2),plike(3)
        call gradlikefunc(nplike,plike,flikemin,gradlike)
        write(*,*)gradlike(1),gradlike(2),gradlike(3)
                
        isitlastcall=.true.          
        call likelihoodfunc(nplike,plike,flikemin)
        call colikeliparampart(nplike,plike,ithetap,
     &   ndesignparam,nvars,
     &   theta(1:(nvars*(nvars+1))/2,1:ndesignparam),
     &   pdist(1:(nvars*(nvars+1))/2,1:ndesignparam),
     &   covars,correls)
     
      endif
!      
      do j=1,ndesignparam
        x(j)=xnew(j)-xmean(j)
      enddo
      istart=0
      do i=1,ixnewpos-1
        istart=istart+kbasvar(i)
      enddo
!
      if(iwhattodo.eq.1.or.iwhattodo.eq.2)then
! make a prediction at x
        call cofuncbasis(ixnewpos,index,ndesignparam,
     &           x,fb)
        irowinc=0
        do i=1,nvars
          index=icompactpostri(ixnewpos,i)
          do j=1,ndesignparam
            theta1(j)=theta(index,j)
            pdist1(j)=pdist(index,j)
          enddo
          if(ixnewpos.ne.i)then
            index=icompactposoff(ixnewpos,i)
          endif
          do nii=1,indexvars(i)
            irowinc=irowinc+1
            do j=1,ndesignparam
              xsamp(j)=xsamplike(irowinc,j)
            enddo
            if(ixnewpos.eq.i)then
              rx(irowinc)=covars(i)*corr_like(x,xsamp,
     &                 theta1,pdist1,ndesignparam)
            else
              rx(irowinc)=dsqrt(covars(i)*covars(ixnewpos))*
     &              correls(index)*corr_like(x,xsamp,
     &                 theta1,pdist1,ndesignparam)
            endif
          enddo
        enddo
        
        ypred=0.0d0
        do i=1,kbasvar(ixnewpos)
          ypred=ypred+fb(i)*coeffbas(i+istart)
        enddo
        do i=1,nlikepoints
          ypred=ypred+delvector(i)*rx(i)
        enddo

        rootmeanerror=covars(ixnewpos)
        do i=1,kbasis
          call symmatindex(kbasis,i,kpvt)
          flikemin=0.0d0
          do j=1,kbasis
            if(j.gt.istart.and.j.le.(istart+kbasvar(ixnewpos)))then
              flikemin=flikemin+FtRinvF_inv(kpvt(j))*fb(j-istart)
            endif
          enddo
          if(i.gt.istart.and.i.le.(istart+kbasvar(ixnewpos)))then
            rootmeanerror=rootmeanerror+fb(i-istart)*flikemin
          endif
        enddo

        do i=1,kbasis
          flikemin=0.0d0
          do j=1,nlikepoints
            flikemin=flikemin+fbassamp(j,i)*rx(j)
          enddo
          if(i.gt.istart.and.i.le.(istart+kbasvar(ixnewpos)))then
            rootmeanerror=rootmeanerror-2.0d0*fb(i)*flikemin
          endif
        enddo

        do i=1,nlikepoints
          call symmatindex(nlikepoints,i,kpvt)
          flikemin=0.0d0
          do j=1,nlikepoints
            flikemin=flikemin+work(kpvt(j))*rx(j)
          enddo
          rootmeanerror=rootmeanerror+rx(i)*flikemin
        enddo
        if(rootmeanerror.lt.0.0d0)then
!  computational error
          rootmeanerror=0.0d0
        endif
        rootmeanerror=dsqrt(rootmeanerror)
      endif
!
      if(iwhattodo.eq.2.or.iwhattodo.eq.3)then
! calculate the derivative of the approximated function with respect to x
        call cogradbasx(ixnewpos,kbasvar(ixnewpos),ndesignparam,
     &    x,gradient(1:kbasvar(ixnewpos),1:ndesignparam))
        do i=1,ndesignparam
          grad_xnew(i)=0.0d0
          do j=1,kbasvar(ixnewpos)
            grad_xnew(i)=grad_xnew(i)+gradient(j,i)*coeffbas(j+istart)
          enddo
        enddo
        
        irowinc=0
        do i=1,nvars
          if(i.eq.ixnewpos)then
            flikemin=covars(ixnewpos)
          else
            index=icompactposoff(i,ixnewpos)
            flikemin=dsqrt(covars(i)*covars(ixnewpos))*
     &               correls(index)
          endif
          index=icompactpostri(i,ixnewpos)
          do j=1,ndesignparam
            theta1(j)=theta(index,j)
            pdist1(j)=pdist(index,j)
          enddo
          do j=1,indexvars(i)
            irowinc=irowinc+1
            do index=1,ndesignparam
              xsamp(index)=xsamplike(irowinc,index)
            enddo
            call gradcorratx(x,xsamp,ndesignparam,theta1,
     &                       pdist1,rx) 
            do index=1,ndesignparam
              grad_xnew(index)=grad_xnew(index)+flikemin*
     &                         rx(index)*delvector(irowinc)
            enddo
          enddo
        enddo
      endif
      return
      end
     
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&           
      subroutine colikelihoodgrad(nplike,plike,dogradient,
     &      fvalue,gradlike,ierr)
      implicit none
      include '../maxlikelihood/colikelihood.h'
!
! calculate the transformed likelihood function value and derivatives
! with respect to input parameters if specified
!
!----------- Arguments ------------------------------------
      integer nplike,ierr
      logical dogradient
      double precision plike(nplike),fvalue,gradlike(nplike)
!
!=> nplike: the number of likelihood parameters to be estimated
!=> plike:  the likelihood parameters to be estimated
!=> dogradient: whether to do derivative calculations
!<= fvalue: the transformed likelihood function value
!<= gradlike: derivatives at plike
!<= ierr: error messages
!    =0: every thing ok
!    <0: parameter out of bounds. |ierr| denotes the out-of-bound parameter
!    =11: the correlation matrix is not postitive definite (the determinant is negative)

!----------- End of list of arguments ---------------------       
!------------- Local variables ----------------------------
      integer i,j,n,t,inert(3),index,nii,
     &   irowinc,getwhichvar,ivar,jvar,icount,k,
     &   icompactposoff,icompactpostri,iloop,istart,jstart
      double precision lnabsdet,signunit,delta,corr_like,
     &    covars(nvars),theta(nvars*(nvars+1)/2,ndesignparam),
     &    pdist(nvars*(nvars+1)/2,ndesignparam),
     &    correls(nvars*(nvars-1)/2),xi(ndesignparam),
     &    xj(ndesignparam),theta1(ndesignparam),
     &    pdist1(ndesignparam),derR(nlikepoints*(nlikepoints+1)/2),
     &    sumi,sumt,gradtheta(ndesignparam),gradp(ndesignparam)
      
      do i=1,nplike
        if(plike(i).lt.plikemin(i).or.
     &     plike(i).gt.plikemax(i))then
! penalize out-of-bounds parameters
          fvalue=1.0d+100
          ierr=-i
          return
        endif
      enddo

      call colikeliparampart(nplike,plike,ithetap,
     &   ndesignparam,nvars,
     &   theta(1:(nvars*(nvars+1))/2,1:ndesignparam),
     &   pdist(1:(nvars*(nvars+1))/2,1:ndesignparam),     
     &   covars,correls)
     
      t=0
      do j=1,nlikepoints
        jvar=getwhichvar(j,nvars,indexvars)
        do n=1,ndesignparam
          xi(n)=xsamplike(j,n)
        enddo
        do i=1,j
          t=t+1
          if(i.eq.j)then
            work(t)=covars(jvar)
          else
            ivar=getwhichvar(i,nvars,indexvars)
            k=icompactpostri(ivar,jvar)
            do n=1,ndesignparam
              xj(n)=xsamplike(i,n)
              theta1(n)=theta(k,n)
              pdist1(n)=pdist(k,n)
            enddo
            if(ivar.eq.jvar)then
              work(t)=covars(jvar)*corr_like(xj,xi,
     &               theta1,pdist1,ndesignparam)
            else
              k=icompactposoff(ivar,jvar)
              work(t)=dsqrt(covars(jvar)*covars(ivar))*
     &          correls(k)*corr_like(xj,xi,
     &               theta1,pdist1,ndesignparam)
            endif
          endif
        enddo
      enddo

      call dspfa(work,nlikepoints,kpvt,i)    
      j=11
      call dspdi(work,nlikepoints,kpvt,derR,inert,delvector,j)
! work now stores the upper triangle of the inverse covariance matrix
! determinant=derR(1)*10.0d0**derR(2)
      signunit=dsign(1.0d0,derR(1))
      lnabsdet=dlog(dabs(derR(1)))+derR(2)*dlog(10.0d0)
      
      if(signunit.lt.0.0d0)then
!        write(*,*)'Wrong! the covariance matrix has a negative'
!        write(*,*)'determinant. Program stops in likelihoodfunc'        
        write(*,*)signunit,lnabsdet
        ierr=11
        return
      endif
      
      t=0
      kbasis=0
      do i=1,nvars
        do j=1,indexvars(i)
          t=t+1
          do n=1,ndesignparam
            xi(n)=xsamplike(t,n)
          enddo
          call cofuncbasis(i,kbasvar(i),ndesignparam,
     &           xi,delvector)
          do n=1,kbasvar(i)
            fbassamp(t,kbasis+n)=delvector(n)
          enddo
        enddo
        if(i.gt.1)then
          do j=t-indexvars(i)+1,t
            do n=1,kbasis
              fbassamp(j,n)=0.0d0
            enddo
          enddo
          do j=1,t
            do n=kbasis+1,kbasis+kbasvar(i)
              fbassamp(j,n)=0.0d0
            enddo
          enddo
        endif
        kbasis=kbasis+kbasvar(i)
      enddo
      call solvebassystem(delta)

      fvalue=delta+lnabsdet
!
! The actual likelihood function =exp(-(fvalue+nlikepoints*ln(2*pi))/2)

      if(dogradient.eqv..false..or.dogradient.eqv..FALSE.)then
        ierr=0
        return
      endif
      if(isitlastcall.eqv..true..or.isitlastcall.eqv..TRUE.)then
        ierr=0
        return
      endif
! 
! section for derivatives starts
!      
! delvector=R^(-1)(Y-FB) was computed in solvebassystem and will not be altered from now
!----Section for calculating derivatives with respect to theta and pdist-------      
      icount=0
      do iloop=1,ithetap
        do k=1,nvars*(nvars+1)/2
          call ivarjvartri(k,nvars,indexvars,
     &                  ivar,jvar,istart,jstart)
          if(ivar.ne.jvar)then
            nii=icompactposoff(ivar,jvar)
          endif
          do n=1,ndesignparam
            theta1(n)=theta(k,n)
            pdist1(n)=pdist(k,n)
          enddo
          ierr=0
          do j=1,nlikepoints
            do i=1,j
              ierr=ierr+1
              derR(ierr)=0.0d0
            enddo
          enddo
          do t=1,ndesignparam
            icount=icount+1
            do j=jstart,jstart+indexvars(jvar)-1
              do n=1,ndesignparam
                xj(n)=xsamplike(j,n)
              enddo
              if(ivar.eq.jvar)then
                kpvt(1)=j-1
              else
                kpvt(1)=istart+indexvars(ivar)-1
              endif
              do i=istart,kpvt(1)
                do n=1,ndesignparam
                  xi(n)=xsamplike(i,n)
                enddo
                ierr=icompactpostri(i,j)
                if(ivar.eq.jvar)then
                  if(iloop.eq.1)then
                    call gradcorr_liketheta(xi,xj,
     &                theta1,pdist1,ndesignparam,gradtheta)
                    derR(ierr)=covars(ivar)*gradtheta(t)
                  else
                    call gradcorr_likep(xi,xj,
     &                theta1,pdist1,ndesignparam,gradp)
                    derR(ierr)=covars(ivar)*gradp(t)
                  endif
                else            
                  if(iloop.eq.1)then
                    call gradcorr_liketheta(xi,xj,
     &                theta1,pdist1,ndesignparam,gradtheta)
                    derR(ierr)=dsqrt(covars(ivar)*covars(jvar))*
     &              correls(nii)*gradtheta(t)  
                  else
                    call gradcorr_likep(xi,xj,
     &                theta1,pdist1,ndesignparam,gradp)
                    derR(ierr)=dsqrt(covars(ivar)*covars(jvar))*
     &              correls(nii)*gradp(t)
                  endif
                endif
              enddo
            enddo
            gradlike(icount)=0.0d0
            do i=1,nlikepoints
              call symmatindex(nlikepoints,i,kpvt)
              sumi=0.0d0
              do j=1,nlikepoints
                sumi=sumi+derR(kpvt(j))*delvector(j)
              enddo
              gradlike(icount)=
     &              gradlike(icount)+delvector(i)*sumi
            enddo
            gradlike(icount)=-gradlike(icount)
!
!Now, the trace part        
            sumi=0.0d0
            do i=1,nlikepoints
              call symmatindex(nlikepoints,i,kpvt)
              sumt=0.0d0
              do j=1,nlikepoints
                sumt=sumt+work(kpvt(j))*derR(kpvt(j))        
              enddo
              sumi=sumi+sumt
            enddo
! sumi is the trace
            gradlike(icount)=gradlike(icount)+sumi
          enddo
        enddo
      enddo
!-------------------------------------------------------------------
! ------- Section for calculating derivatives with respect to variance
      do k=1,nvars
        icount=icount+1
        ierr=0
        do j=1,nlikepoints
          do i=1,j
            ierr=ierr+1
            derR(ierr)=0.0d0
          enddo
        enddo
! first, vertical elements
        do i=1,k-1
          ierr=icompactpostri(i,k)
          do t=1,ndesignparam
            theta1(t)=theta(ierr,t)
            pdist1(t)=pdist(ierr,t)
          enddo
          ierr=icompactposoff(i,k)
          call istartjstart(i,k,nvars,indexvars,istart,jstart)
          do n=istart,istart+indexvars(i)-1
            do t=1,ndesignparam
              xi(t)=xsamplike(n,t)
            enddo
            do j=jstart,jstart+indexvars(k)-1
              do t=1,ndesignparam
                xj(t)=xsamplike(j,t)
              enddo
              ivar=icompactpostri(n,j)
              derR(ivar)=0.5d0*dsqrt(covars(i)/covars(k))*
     &               correls(ierr)*corr_like(xi,xj,
     &               theta1,pdist1,ndesignparam)
            enddo
          enddo
        enddo
        ierr=icompactpostri(k,k)
        do t=1,ndesignparam
          theta1(t)=theta(ierr,t)
          pdist1(t)=pdist(ierr,t)
        enddo
        call istartjstart(k,k,nvars,indexvars,istart,jstart)
        do j=jstart,jstart+indexvars(k)-1
          do t=1,ndesignparam
            xi(t)=xsamplike(j,t)
          enddo
          do n=istart,j
            ivar=icompactpostri(n,j)
            if(n.eq.j)then
              derR(ivar)=1.0d0
            else
              do t=1,ndesignparam
                xj(t)=xsamplike(n,t)
              enddo
              derR(ivar)=corr_like(xi,xj,
     &               theta1,pdist1,ndesignparam)
            endif
          enddo
        enddo
! then, horizontal elements 
        do i=k+1,nvars
          ierr=icompactpostri(k,i)
          do t=1,ndesignparam
            theta1(t)=theta(ierr,t)
            pdist1(t)=pdist(ierr,t)
          enddo
          ierr=icompactposoff(k,i)
          call istartjstart(k,i,nvars,indexvars,istart,jstart)
          do n=istart,istart+indexvars(i)-1
            do t=1,ndesignparam
              xi(t)=xsamplike(n,t)
            enddo
            do j=jstart,jstart+indexvars(k)-1
              do t=1,ndesignparam
                xj(t)=xsamplike(j,t)
              enddo
              ivar=icompactpostri(n,j)
              derR(ivar)=0.5d0*dsqrt(covars(i)/covars(k))*
     &               correls(ierr)*corr_like(xi,xj,
     &               theta1,pdist1,ndesignparam)
            enddo
          enddo
        enddo
 
        gradlike(icount)=0.0d0
        do i=1,nlikepoints
          call symmatindex(nlikepoints,i,kpvt)
          sumi=0.0d0
          do j=1,nlikepoints
            sumi=sumi+derR(kpvt(j))*delvector(j)
          enddo
          gradlike(icount)=
     &        gradlike(icount)+delvector(i)*sumi
        enddo
        gradlike(icount)=-gradlike(icount)
!
!Now, the trace part        
        sumi=0.0d0
        do i=1,nlikepoints
          call symmatindex(nlikepoints,i,kpvt)
          sumt=0.0d0
          do j=1,nlikepoints
            sumt=sumt+work(kpvt(j))*derR(kpvt(j))        
          enddo
          sumi=sumi+sumt
        enddo

! sumi is the trace
        gradlike(icount)=gradlike(icount)+sumi
      enddo
!--------------------------------------------------------------------    
! ------- Section for calculating derivatives with respect to correls
      do k=1,nvars*(nvars-1)/2
        icount=icount+1
        ierr=0
        do j=1,nlikepoints
          do i=1,j
            ierr=ierr+1
            derR(ierr)=0.0d0
          enddo
        enddo
        call ivarjvaroff(k,nvars,indexvars,
     &                  ivar,jvar,istart,jstart)
        ierr=icompactpostri(ivar,jvar)
        do t=1,ndesignparam
          theta1(t)=theta(ierr,t)
          pdist1(t)=pdist(ierr,t)
        enddo
        do i=istart,istart+indexvars(ivar)-1
          do t=1,ndesignparam
            xi(t)=xsamplike(i,t)
          enddo
          do j=jstart,jstart+indexvars(jvar)-1
            do t=1,ndesignparam
              xj(t)=xsamplike(j,t)
            enddo
            ierr=icompactpostri(i,j)
            derR(ierr)=dsqrt(covars(ivar)*covars(jvar))*
     &         corr_like(xi,xj,theta1,pdist1,ndesignparam)
          enddo
        enddo

        gradlike(icount)=0.0d0
        do i=1,nlikepoints
          call symmatindex(nlikepoints,i,kpvt)
          sumi=0.0d0
          do j=1,nlikepoints
            sumi=sumi+derR(kpvt(j))*delvector(j)
          enddo
          gradlike(icount)=
     &        gradlike(icount)+delvector(i)*sumi
        enddo
        gradlike(icount)=-gradlike(icount)

!Now, the trace part        
        sumi=0.0d0
        do i=1,nlikepoints
          call symmatindex(nlikepoints,i,kpvt)
          sumt=0.0d0
          do j=1,nlikepoints
            sumt=sumt+work(kpvt(j))*derR(kpvt(j))        
          enddo
          sumi=sumi+sumt
        enddo
! sumi is the trace
        gradlike(icount)=gradlike(icount)+sumi
      enddo
!--------------------------------------------------------------
      ierr=0
      return
      end
      
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      subroutine solvebassystem(delta)
      implicit none
      include '../maxlikelihood/colikelihood.h'
!
! Calculating the best linear estimates of the coefficients of the basis functions,
! delta (the residual square error), and matrices needed for prediction and estimation
! of prediction error

      integer kcopy(kbasis),inert(3)
      double precision FR(kbasis,nlikepoints),
     &    a(kbasis*(kbasis+1)/2),rhs(kbasis),sum,
     &    acopy(max0(kbasis*(kbasis+1)/2,nlikepoints)),
     &    amat(kbasis*(kbasis+1)/2),
     &    rhscopy(kbasis),delta
      
      integer n,j,t,i,mark,irowinc,index,nii
      
      do i=1,nlikepoints
        call symmatindex(nlikepoints,i,kpvt)
        do j=1,kbasis
          FR(j,i)=0.0d0
          do t=1,nlikepoints
            FR(j,i)=FR(j,i)+work(kpvt(t))*fbassamp(t,j)
          enddo
        enddo
      enddo  
      
      do i=1,kbasis
! right hand side of the linear system
        rhs(i)=0.0d0
        do j=1,nlikepoints
          rhs(i)=rhs(i)+FR(i,j)*ysamplike(j)
        enddo
! make a copy
        rhscopy(i)=rhs(i)
      enddo

! coefficient matrix, symmetric      
      t=0
      do j=1,kbasis
        do i=1,j
          t=t+1
          a(t)=0.0d0
          do n=1,nlikepoints    
            a(t)=a(t)+FR(i,n)*fbassamp(n,j)         
          enddo
! make a copy
          amat(t)=a(t)
        enddo
      enddo
      
! solve the basis function linear system
      call dspfa(a,kbasis,kpvt,mark)
      do i=1,kbasis*(kbasis+1)/2
        acopy(i)=a(i)
      enddo
      do i=1,kbasis
        kcopy(i)=kpvt(i)
      enddo
      if(isitlastcall.eqv..true..or.isitlastcall.eqv..TRUE.)then
! calculate the inverse only
        do j=1,kbasis*(kbasis+1)/2
          FtRinvF_inv(j)=a(j)
        enddo
        j=1
        call dspdi(FtRinvF_inv,kbasis,kpvt,coeffbas,
     &             inert,delvector,j)
        do i=1,kbasis
          kpvt(i)=kcopy(i)
        enddo
      endif
      call dspsl(a,kbasis,kpvt,rhs)
      do i=1,kbasis
        coeffbas(i)=rhs(i)
      enddo
      
! one step improvement
      do i=1,kbasis
        call symmatindex(kbasis,i,kpvt)
        sum=0.0d0
        do t=1,kbasis
          sum=sum+amat(kpvt(t))*coeffbas(t)
        enddo
        rhs(i)=(sum-rhscopy(i))*1.0d+9
      enddo
      call dspsl(acopy,kbasis,kcopy,rhs)
      do i=1,kbasis
        coeffbas(i)=coeffbas(i)-rhs(i)*1.0d-9
      enddo
      
      do i=1,nlikepoints
        acopy(i)=0.0d0
        do j=1,kbasis
          acopy(i)=acopy(i)+coeffbas(j)*fbassamp(i,j)
        enddo
        acopy(i)=ysamplike(i)-acopy(i)
      enddo      

      do t=1,nlikepoints
        call symmatindex(nlikepoints,t,kpvt)
        delvector(t)=0.0d0
        do i=1,nlikepoints
          delvector(t)=delvector(t)+work(kpvt(i))*acopy(i)
        enddo
      enddo
      delta=0.0d0
      do t=1,nlikepoints
        delta=delta+delvector(t)*acopy(t)
      enddo
      
      if(isitlastcall.eqv..true..or.isitlastcall.eqv..TRUE.)then
! All parameters in the likelihood function have been optimized. Prepare
! matrices for MSE calculations
        do i=1,kbasis
          call symmatindex(kbasis,i,kpvt)
          do j=1,nlikepoints
            fbassamp(j,i)=0.0d0
            do t=1,kbasis
              fbassamp(j,i)=fbassamp(j,i)+FtRinvF_inv(kpvt(t))*
     &                      FR(t,j)
            enddo
          enddo
        enddo
        t=0
        do j=1,nlikepoints
          do i=1,j
            t=t+1
            work(t)=-work(t)
            do n=1,kbasis
              work(t)=work(t)+fbassamp(i,n)*FR(n,j)
            enddo
          enddo
        enddo 
      endif
      return
      end
!######################################################

      DOUBLE PRECISION FUNCTION df1dimd_like(x)
      INTEGER NMAX
      double precision x
      PARAMETER (NMAX=50)
CU    USES gradlikefunc
      INTEGER j,ncom
      double precision df(NMAX),pcom(NMAX),xicom(NMAX),
     &    xt(NMAX)
,dummy
      COMMON /f1com/ pcom,xicom,ncom
      do 11 j=1,ncom
        xt(j)=pcom(j)+x*xicom(j)
11    continue
      call gradlikefunc(ncom,xt,dummy,df)
      df1dimd_like=0.0d0
      do 12 j=1,ncom
        df1dimd_like=df1dimd_like+df(j)*xicom(j)
12    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.

      subroutine gradlikefunc(nplike,plike,fvalue,gradlike)
      implicit none
      integer nplike,ierr
      double precision plike(nplike),fvalue,gradlike(nplike)
      logical dogradient
      
      dogradient=.true.     
      call colikelihoodgrad(nplike,plike,dogradient,
     &      fvalue,gradlike,ierr)
      return
      end
      
      subroutine likelihoodfunc(nplike,plike,fvalue)
      implicit none
      integer nplike,ierr
      double precision plike(nplike),fvalue,gradlike(nplike)
      logical dogradient
      
      dogradient=.false.     
      call colikelihoodgrad(nplike,plike,dogradient,
     &      fvalue,gradlike,ierr)
      return
      end
      
      double precision function f1dim_like(x)
      INTEGER NMAX
      double precision x
      PARAMETER (NMAX=50)
CU    USES likelihood
      INTEGER j,ncom
      double precision pcom(NMAX),xicom(NMAX),xt(NMAX)
      COMMON /f1com/ pcom,xicom,ncom

      do 11 j=1,ncom
        xt(j)=pcom(j)+x*xicom(j)
11    continue
      call likelihoodfunc(ncom,xt,f1dim_like)
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.

!#################################################################
      subroutine gradfunkmin(ndim,x,fatx,fjac)
      implicit none
      
      integer ndim,j
      double precision x(ndim),fatx,fjac(ndim),h,eps,
     &  zero,temp,fxh
      
      call likelihoodfunc(ndim,x,fatx)
      
      eps = 1.0d-5
      zero=0.0d0
      
      do j = 1, ndim
        temp = x(j)
        h = eps*dabs(temp)
        if (h .eq. zero) h = eps
        x(j) = temp + h
        call likelihoodfunc(ndim,x,fxh)
        x(j) = temp
        fjac(j) = (fxh - fatx)/h
      enddo

      return
      end subroutine gradfunkmin

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&     
      subroutine randpermut(npoints,ndim,x)
      implicit none
!
! conduct random permutation
      integer npoints,ndim,i,j,index
      double precision x(npoints,ndim),xtemp(npoints),
     &                ran2
      
      do i=1,ndim
        do j=1,npoints
          xtemp(j)=x(j,i)
        enddo
        do j=1,npoints
          index=int(dble(npoints-j+1)*ran2()+1.0d0)
          x(j,i)=xtemp(index)
          xtemp(index)=xtemp(npoints-j+1)
        enddo
      enddo    
      return
      end

c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
c
c####################################################################
c random number generator
c
      double precision function ran2()
      implicit none
      integer im1,im2,imm1,ia1,ia2,iq1,iq2,ir1,ir2,ntab,ndiv
      double precision am,eps,rnmx
      parameter(im1=2147483563,im2=2147483399,am=1.0d0/dble(im1),
     &imm1=im1-1,ia1=40014,ia2=40692,iq1=53668,iq2=52774,ir1=
     &12211,ir2=3791,ntab=32,ndiv=1+imm1/ntab,eps=1.2d-7,
     &rnmx=1.0d0-eps)
      integer idum2,j,k,iv(ntab),iy,idum
      save iv,iy,idum2,idum
      data idum2/123456789/,iv/ntab*0/,iy/0/,idum/-1/
      
      if(idum.le.0) then
        idum=max0(-idum,1)
        idum2=idum
        do 11 j=ntab+8,1,-1
          k=idum/iq1
          idum=ia1*(idum-k*iq1)-k*ir1
          if(idum.lt.0) idum=idum+im1
          if(j.le.ntab) iv(j)=idum
11      continue
        iy=iv(1)
      end if
      k=idum/iq1
      idum=ia1*(idum-k*iq1)-k*ir1
      if(idum.lt.0)idum=idum+im1
      k=idum2/iq2
      idum2=ia2*(idum2-k*iq2)-k*ir2
      if(idum2.lt.0) idum2=idum2+im2
      j=1+iy/ndiv
      iy=iv(j)-idum2
      iv(j)=idum
      if(iy.lt.1)iy=iy+imm1
      ran2=dmin1(am*dble(iy),rnmx)
      return
      end
!######################################################