*...:....1....:....2....:....3....:....4....:....5....:....6....:....7....:....8
*    
*     plot.f
*
*     Copyright (C) 1991.
*
*     Marie Davidian and A. Ronald Gallant 
*     P.O. Box 5513, Raleigh NC 27650-5513, USA   
*
*     Permission to use, copy, modify, and distribute this
*     software and its documentation for any purpose and
*     without fee is hereby granted, provided that the
*     above copyright notice appear in all copies and that
*     both that copyright notice and this permission notice
*     appear in supporting documentation.
*
*     This software is provided "as is" without any
*     expressed or implied warranty.
*
*...:....1....:....2....:....3....:....4....:....5....:....6....:....7....:....8
      subroutine plot()
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      save

      include "global.f"

      integer*4 i,j,ic

      character*1 digit
      character*12 fname


      dimension digit(0:9)

      data digit /'0','1','2','3','4','5','6','7','8','9'/

      integer*4 maxM

      parameter (maxM=99)

*     Compute points for plotting the joint density if jt=1.

      if (jt.eq.1) then
        call plotj()
        return
      end if

*     If M=1 the joint and the marginal are the same; thus, compute
*     the joint density if M=1.

      if (M.eq.1) then
        call plotj()
        return
      end if
       
*     Compute the marginals if jt=0.  Each marginal is put in a 
*     different file.  The names of these files are put in outfil.

      if (M.gt.maxM) then
        write(3,3901) 
        stop
      end if

      fname=outfil//'            '

      j=0
      do i=1,12
        if (fname(i:i).eq.'.') then
          j=i
        end if
        if ((fname(i:i).eq.' ').and.(j.eq.0)) then
          j=i
          fname(j+1:j+1)='p'
        end if
      end do

      if (j.eq.0) then
        fname(9:9)='.'
        j=9
      else 
        fname(j:j)='.'
      end if

      open(unit=19,file=outfil,status='unknown',form='formatted')

      do ic=1,M
        fname(j+2:j+3)=digit(10*(ic/10))//digit(ic-10*(ic/10))
        if (fname.eq.outfil) then
          fname(j+2:j+3)=digit(0)//digit(0)
        end if
        write(19,19001) ic,fname
        write(19,19002)
        write(19,19003)
        write(19,19004)
        write(19,19005)
        write(19,19006)
        write(19,19005)
        write(19,19006)
        write(19,19007)
        write(19,19007)
        write(19,19007)
        open(unit=9,file=fname,status='unknown',form='formatted')
        call plotm(ic)
        close(9)
      end do

      close(19)
      return
3901  format(1x,'error, plot, Can''t handle M that big.')
19001 format(1x,'Plot points for variable',i3,' written to file ',a12)
19002 format(1x,'   Ez    format(e30.18)')
19003 format(1x,'   Cz    format(e30.18)')
19004 format(1x,'   zinc  format(e30.18)')
19005 format(1x,'   z     format(e30.18)')
19006 format(1x,'   hplt  format(e30.18)')
19007 format(1x,'   .')
      end

      subroutine plotj()
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      save

      include "global.f"

      integer*4 i,j
      integer*4 lwork
      integer*4 izz

      real*8    zinc
      real*8    h,Eh,del,hplt
      real*8    work

      integer*4 mwork

      parameter (mwork=500)

      dimension zinc(mM)
      dimension del(mtheta)
      dimension work(mwork)
      dimension izz(mM)

      lwork=1+M+M*(M+1)/2+M+M
      if (lwork.gt.mwork) then
        write(3,3001) lwork
        stop
      end if

      lwork=2*ntermz
      if (lwork.gt.mwork) then
        write(3,3001) lwork
        stop
      end if

      open(unit=9,file=outfil,status='unknown',form='formatted')

      do i=1,M
        write(9,9001) Ez(i)
      end do

      do j=1,M
      do i=1,j
        write(9,9001) Cz(j*(j-1)/2+i)
      end do
      end do

      do i=1,M
        zinc(i)=sfac*dsqrt(Cz(i*(i+1)/2))/dfloat(ngrd)
      end do

      do i=1,M
        write(9,9001) zinc(i)
      end do

      do i=1,M
        izz(i)=-ngrd
      end do
      izz(1)=izz(1)-1

      do j=1,(2*ngrd+1)**M

        do i=1,M
          if (izz(i).lt.ngrd) then
            izz(i)=izz(i)+1
            go to 10
          else
            izz(i)=-ngrd
          end if
        end do
10      continue

        do i=1,M
          z(i)=Ez(i)+dfloat(izz(i))*zinc(i)
        end do

        do i=1,M
          write(9,9001) z(i)
        end do

        call hyx(z,x,h,del,
     &    theta(iaA),theta(iaB0),theta(iaB),theta(iaR),
     &    M,L,ntermz,ntermx,Kz,Kx,
     &    work)

        call ehyx(z,x,Eh,del,
     &    theta(iaA),theta(iaB0),theta(iaB),theta(iaR),
     &    M,L,ntermz,ntermx,Kz,Kx,
     &    work)

        hplt=h/Eh
        if (hplt.le.1.d-60) hplt=0.d0

        write(9,9001) hplt

      end do

      close(9)
      return
3001  format(1x,'error, plot, increase mwork to at least',i5)
9001  format(e30.18)
      end

      subroutine plotm(ic)
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      save

      include "global.f"

*     Local data

      real*8 b0,b,V
      real*8 zinc,zic
      real*8 hplt
      real*8 work
      real*8 T,W,u,Tjj
      real*8 endpts
      real*8 sum
      real*8 fctr
      real*8 P,PwA,Pwz,Eh,Ehdel
      real*8 R0,tQt,uQu,detR,detV
      real*8 weight

      integer*4 ic
      integer*4 igrid,ier
      integer*4 lwork
      integer*4 npoint

      integer*4 i,j,ii,jj
      integer*4 izz

*     Local parameters

      real*8 eps
      real*8 rootpi
      real*8 root2

      integer*4 mwork
      integer*4 mT

      parameter (rootpi=1.7724538509055160272d0)
      parameter (root2 =1.4142135623730950488d0)
      parameter (mT=30)
      parameter (eps=1.d-13)
      parameter (mwork=500)

      dimension b0(mM),b(mM),V(mM*mM)
      dimension T(mT),W(mT),u(mM),work(mwork)
      dimension endpts(2)
      dimension izz(mM)
      dimension PwA(mtheta),Pwz(mM),Ehdel(mtheta)
      dimension R0(mM*mM)

      npoint=min0(npt,mT)

*     Error checks

      lwork=1+M+(M*(M+1))/2+M+M
      if (lwork.gt.mwork) then
        write(3,3901) lwork
        stop
      end if

      lwork=2*ntermz
      if (lwork.gt.mwork) then
        write(3,3901) lwork
        stop
      end if

      lwork=2*(M*M+M+2)
      if (lwork.gt.mwork) then
        write(3,3901) lwork
        stop
      end if

      lwork=((M*(M-1))/2)**2
      if (lwork.gt.mwork) then
        write(3,3901) lwork
        stop
      end if

      lwork=npoint
      if (lwork.gt.mwork) then
        write(3,3901) lwork
        stop
      end if

      lwork=2*ntermz
      if (lwork.gt.mwork) then
        write(3,3901) lwork
        stop
      end if

*     The algorithm is an importance function method with a conditional
*     normal as the importance function using quadrature to integrate 
*     rather than using Monte Carlo.

      write(9,9001) Ez(ic)
      write(9,9001) Cz((ic*(ic+1))/2)
      zinc=sfac*dsqrt(Cz((ic*(ic+1))/2))/dfloat(ngrd)
      write(9,9001) zinc

*     Get the parameters b0, b, and V of the conditional normal.

      call dmcvt(Cz,R0,M,1,0)
      call cndist(Ez,R0,M,ic,b0,b,V,work,mwork)
      call dmcvt(V,V,M-1,0,1)
      call factor(V,work,M-1,eps,ier)
      if (ier.ne.0) write(3,3902) ier

      detV=1.d0
      do i=1,(M-1)
        detV=detV*V((i*(i+1))/2)
      end do 
      detV=dabs(detV)

*     Get Gauss-Hermite quadrature points.

      call gaussq(4,npoint,0.d0,0.d0,0,endpts,work,T,W)

      do i=1,npoint
        T(i)=root2*T(i)
      end do

      fctr=1.d0/(root2*rootpi**M)

      do i=1,lR
        R0(i)=theta(ibR+i)
      end do

      call drinv(R0,M)

      detR=1.d0
      do i=1,M
        detR=detR*theta(ibR+(i*(i+1))/2)
      end do
      detR=dabs(detR)

      do igrid=-ngrd,ngrd

        zic=Ez(ic)+dfloat(igrid)*zinc

        do i=1,(M-1)
          izz(i)=1
        end do
        izz(1)=izz(1)-1

        sum=0.d0

        do j=1,npoint**(M-1)

          do jj=1,(M-1)
            if (izz(jj).lt.npoint) then
              izz(jj)=izz(jj)+1
              go to 10
            else
              izz(jj)=1
            end if
          end do
 10       continue

          do jj=1,(M-1)
            u(jj)=0.d0
          end do

          tQt=0.d0

          do jj=1,(M-1)
            Tjj=T(izz(jj))
            tQt=tQt+Tjj**2
            do ii=1,jj
              u(ii)=u(ii)+V(jj*(jj-1)/2+ii)*Tjj
            end do
          end do

          tQt=-tQt/2.d0

          ii=0
          do jj=1,M
            if (jj.eq.ic) then
              z(jj)=zic
            else 
              ii=ii+1
              z(jj)=b0(ii)+b(ii)*zic+u(ii)
            end if 
          end do

          do ii=1,M
            u(ii)=0.d0
          end do

          do jj=1,M
          do ii=1,jj
           u(ii)=u(ii)+R0(jj*(jj-1)/2+ii)*z(jj)
          end do
          end do

          uQu=0.d0
          do ii=1,M
            uQu=uQu+u(ii)**2
          end do
          uQu=-uQu/2.d0

          call pzx(u,x,P,PwA,Pwz,theta(iaA),
     &      M,L,ntermz,ntermx,work)

          call ehyx(u,x,Eh,Ehdel,
     &      theta(iaA),theta(iaB0),theta(iaB),theta(iaR),
     &      M,L,ntermz,ntermx,Kz,Kx,work)

          weight=1.d0
          do jj=1,(M-1)
            weight=weight*W(izz(jj))
          end do

          sum=sum+weight*(P/Eh)*fctr*(detV/detR)*exp(uQu-tQt)

        end do

        hplt=sum
        if (hplt.le.1.d-60) hplt=0.d0

        write(9,9001) zic
        write(9,9001) hplt

      end do

      return
3901  format(1x,'error, plotm, increase mwork to at least',i5)
3902  format(1x,'error, plotm, error return from factor, ier=',i5)
9001  format(e30.18)
      end
      
*     Given mu and sigma of an M variate multivariate normal X and an
*     index ic, returns b0, b, and V such that the M-1 variate normal
*     with mean b0 + b*X2 and variance V is the conditional distribution
*     of X1 = [X(1),...,X(ic-1),X(ic+1),...,X(M)]' given X2 = X(ic).

*     The formulas for the conditional normal are from p. 37 of
*     T. W. Anderson, An Introduction to Multivariate Statistical
*     Analysis, Second Edition, Wiley, 1984.

      subroutine cndist(mu,sigma,M,ic,b0,b,V,work,lwork)
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      save

      real*8 mu,sigma
      real*8 b0,b,V
      real*8 work

      integer*4 M,ic,lwork
      integer*4 i,j,ii,jj
      integer*4 ibS11,ibS12,ibS22,ibS22inv
      integer*4 ibU1,ibU2
      integer*4 ibW
      integer*4 length
      integer*4 iaS11,iaS12,iaS22,iaS22inv
      integer*4 iaU1,iaU2
      integer*4 iaW

      dimension work(lwork)
      dimension mu(M),sigma(M,M)
      dimension b0(M-1),b(M-1),V((M-1)**2)

      ibS11    = 0
      ibS12    = ibS11    + (M-1)**2
      ibS22    = ibS12    + (M-1)
      ibS22inv = ibS22    + 1
      ibU1     = ibS22inv + 1
      ibU2     = ibU1     + (M-1)
      ibW      = ibU2     + 1
      length   = ibW      + (M-1)**2

      iaS11    = ibS11    + 1
      iaS12    = ibS12    + 1
      iaS22    = ibS22    + 1
      iaS22inv = ibS22inv + 1
      iaU1     = ibU1     + 1
      iaU2     = ibU2     + 1
      iaW      = ibW      + 1

      if (length.gt.lwork) then
        write(3,3901) length
        stop
      end if

      work(iaU2)=mu(ic)

      work(iaS22)=sigma(ic,ic)

      work(iaS22inv)=1.d0/work(iaS22)

      ii=0
      do i=1,M
        if (i.ne.ic) then
          ii=ii+1
          work(ibU1+ii)=mu(i)
          work(ibS12+ii)=sigma(i,ic)
        end if
      jj=0
      do j=1,M
        if ((i.ne.ic).and.(j.ne.ic)) then
          jj=jj+1
          work(ibS11+(M-1)*(jj-1)+ii)=sigma(i,j)
        end if
      end do
      end do

      call dgmprd(work(iaS12),work(iaS22inv),b,M-1,1,1)
      call dgmprd(b,work(iaU2),work(iaW),M-1,1,1)
      call dgmsub(work(iaU1),work(iaW),b0,M-1,1)

      call dgmprd(b,work(iaS12),work(iaW),M-1,1,M-1)
      call dgmsub(work(iaS11),work(iaW),V,M-1,M-1)

      return
3901  format(1x,'error, cndist, increase lwork to at least',i5)
      end
