C
       subroutine gibbs(n,nc,ifood,itemp,itime,iscore,indexftptm,
     1      z,gamma,amu,theta,sqdelta,sig2,indexfd,ntfood,iseed)
c      purpose: to generate a gibbs iteration
       implicit real*8 (a-h,o-z)
       real*8 z(n),gamma(0:9)
       real*8 amu(1:12,1:4,1:9),theta(1:12,1:4),sqdelta(1:12)
       real*8 sig2(1:12)
       integer ifood(n),itemp(n),itime(n),iscore(n)
       integer nc(4),indexftptm(1:12,1:4,1:9)
       integer indexfd(12),ntfood
       integer iseed
c......generate z(n) given amu and gamma
       call Gz(n,ifood,itemp,itime,iscore,amu,gamma,
     1         sig2,indexfd,z,iseed)
       call Gsqdelta(amu,theta,nc,indexfd,indexftptm,
     1               ntfood,sqdelta,iseed)
c......generate sig2
       call Gsig2(n,nc,z,ifood,itemp,itime,indexfd,
     1         indexftptm,amu,sig2,iseed)
c......generate theta given amu and delta
       call Gtheta(amu,nc,sqdelta,indexfd,indexftptm,theta,iseed)
c......generate amu given z,theta,sqdelta
       call Gamu(n,nc,ifood,itemp,itime,iscore,z,indexfd,
     1            indexftptm,theta,sqdelta,sig2,amu,iseed)
       return
       end

       subroutine Gz(n,ifood,itemp,itime,iscore,amu,gamma,
     1               sig2,indexfd,z,iseed)
c      generate z(n) given amu and gamma
       implicit real*8 (a-h,o-z)
       real*8 z(n),gamma(0:9),sig2(1:12)
       real*8 amu(1:12,1:4,1:9)
       real*8 aleft,aright,amean,tn,sig
       integer ifood(n),itemp(n),itime(n),iscore(n)
       integer iseed,indexfd(12)
       logical la,lb
       la=.false.
       lb=.false.
       do 300 i=1,n
         if (indexfd(ifood(i)) .eq. 0) goto 300
         amean=amu(ifood(i),itemp(i),itime(i))
         sig=dsqrt(sig2(ifood(i)))
         if (iscore(i) .eq. 1) then
          la=.true.
          aright=-amean/sig
          tn=YTUVN(aleft,aright,la,lb,iseed)
          z(i)=amean+tn*sig
          la=.false.
         endif
         if ((iscore(i) .ge. 2) .and. (iscore(i) .le. 8)) then
          aleft =(gamma(iscore(i)-1)-amean)/sig
          aright=(gamma(iscore(i))-amean)/sig
          tn=YTUVN(aleft,aright,la,lb,iseed)
          z(i)=amean+tn*sig
         endif
         if (iscore(i) .eq. 9) then
          lb=.true.
          aleft=(gamma(8)-amean)/sig
          tn=YTUVN(aleft,aright,la,lb,iseed)
          z(i)=amean+tn*sig
          lb=.false.
         endif
 300   continue
       return
       end



       subroutine Gsqdelta(amu,theta,nc,indexfd,indexftptm,
     1            ntfood,sqdelta,iseed)
c      generate sqdelta given amu and theta
       implicit real*8 (a-h,o-z)
       real*8 amu(1:12,1:4,1:9),theta(1:12,1:4),sqdelta(1:12)
       real*8 sumin
       real*8  scaler,shape,r(1)
       integer nc(4),indexfd(12),ntfood
       integer indexftptm(1:12,1:4,1:9)
       integer iseed
       external DRNGAM
       do 340 i=1,12
         if (indexfd(i) .eq. 0) goto 340
         sumin=(amu(i,1,1)-theta(i,1))**2
         do 350 j=1,4
           do 355 l=2,nc(j)
             if (indexftptm(i,j,l) .eq. 1) goto 355
             sumin=sumin+(amu(i,j,l)-theta(i,j))**2
 355       continue
 350     continue
cc        shape=23.0d0/2.0d0-1.0d0+3.0000000001d0
cc        scaler=sum/2.0d0+0.047699d0
c        shape=23.0d0/2.0d0-1.0d0+15.767646d0
c        scaler=sumin/2.0d0+0.6567029d0
c        shape=23.0d0/2.0d0-1.0d0+3.2571845d0 !using entrees
c        scaler=sumin/2.0d0+0.4161029d0 !using entrees
         shape=23.0d0/2.0d0-1.0d0+2.298d0 !fit using pastries data
         scaler=sumin/2.0d0+0.0049d0 !fit using pastries data (10*variance)
c        shape=23.0d0/2.0d0-1.0d0+2.0298d0 !fit using pastries data
c        scaler=sumin/2.0d0+0.0039d0 !fit using pastries data (100*variance)
         call DRNGAM(1,shape,r)
         sqdelta(i)=scaler/r(1)
 340   continue
       return
       end

       subroutine Gsig2(n,nc,z,ifood,itemp,itime,indexfd,
     1         indexftptm,amu,sig2,iseed)
       implicit real*8 (a-h,o-z)
       real*8  amu(1:12,1:4,1:9),sig2(1:12)
       real*8  z(n)
       real*8  scaler,shape,r(1)
       real*8  sum(12)
       integer nc(4),indexfd(12)
       integer indexftptm(1:12,1:4,1:9)
       integer iseed
c      integer icount(12)
       integer ifood(n),itemp(n),itime(n)
       external DRNGAM
c      shape=23.0d0*36.0d0/2.0d0+4.4198d0  !fit from entree data
       shape=23.0d0*36.0d0/2.0d0+8.44d0  !fit using pastries data x 10
c      shape=23.0d0*36.0d0/2.0d0+2.64d0  !fit using pastries data x 100
       do 100 i=1,12
         sum(i)=0.0d0
c        icount(i)=0
 100   continue
       do 200 i=1,n
        sum(ifood(i))=sum(ifood(i))+(z(i)-
     1   amu(ifood(i),itemp(i),itime(i)))**2
c       icount(ifood(i))=icount(ifood(i))+1
 200   continue
       do 300 k=1,12
        if (indexfd(k) .eq. 0) goto 300
c       scaler=sum(k)/(2.0d0)+17.875216   !fit from entrees data
        scaler=sum(k)/(2.0d0)+0.42505d0   !fit using pastries data x 10
c       scaler=sum(k)/(2.0d0)+0.09392d0   !fit using pastries data x 100
        call DRNGAM(1,shape,r)
        sig2(k)=scaler/r(1)
c       write(*,*) 'k=',k,'scaler=',scaler,'(828)count=',icount(k)
 300   continue
       return
       end


       subroutine Gtheta(amu,nc,sqdelta,indexfd,indexftptm,theta,iseed)
c......generate theta given amu and delta
       implicit real*8 (a-h,o-z)
       real*8 amu(1:12,1:4,1:9),theta(1:12,1:4),sqdelta(1:12)
       real*8 sum(4),sigma,amean
       integer nc(4),indexfd(12)
       integer indexftptm(1:12,1:4,1:9)
       integer iseed
       real*8 aleft,aright,tn
c      real*8 avgmu,alambda
       integer num(4)
       logical la,lb
c......initialize
       la=.false.
       lb=.false.
       do 360 i=1,12
         if (indexfd(i) .eq. 0) goto 360
         sum(1)=amu(i,1,1)
         num(1)=1
         do 362 i1=2,4
          sum(i1)=0.0d0
          num(i1)=0
 362     continue
         do 365 j=1,4
            do 370 l=2,nc(j)
              if (indexfd(i) .eq. 0) goto 370
              if (indexftptm(i,j,l) .eq. 1) goto 370
              sum(j)=sum(j)+amu(i,j,l)
              num(j)=num(j)+1
 370        continue
 365     continue
c      avgmu=sum(1)/dfloat(num(1))
c      alambda=.0072912d0*10.0d0/(.0072912d0*10.0d0+
c    1           sqdelta(i)/dfloat(num(1))) !using pastries
c      amean=alambda*avgmu+(1.0d0-alambda)*.53138d0 !using pastries
c      sigma=dsqrt( (1.0d0-alambda)*.0072912d0*10.0d0 ) !using pastries
       amean=sum(1)/dfloat(num(1))
       sigma=dsqrt(sqdelta(i)/dfloat(num(1)))
       aleft=(theta(i,2)-amean)/sigma
       lb=.true.
       tn=YTUVN(aleft,aright,la,lb,iseed)
       lb=.false.
       theta(i,1)=amean+sigma*tn
       do 367 j=2,4
c        avgmu=sum(j)/dfloat(num(j))
c        alambda=.0072912d0*10.0d0/(.0072912d0*10.0d0+
c    1           sqdelta(i)/dfloat(num(j))) !using pastries
c        amean=alambda*avgmu+(1.0d0-alambda)*.53138d0 !using pastries
c        sigma=dsqrt( (1.0d0-alambda)*.0072912d0*10.0d0 ) !using pastries
         amean=sum(j)/dfloat( num(j) )
         sigma=dsqrt(sqdelta(i)/dfloat(num(j)))
         if (j .lt. 4) then
          aleft=(theta(i,j+1)-amean)/sigma
          aright=(theta(i,j-1)-amean)/sigma
         else
          aright=(theta(i,j-1)-amean)/sigma
          la=.true.
         endif
         tn=YTUVN(aleft,aright,la,lb,iseed)
         theta(i,j)=amean+sigma*tn
         la=.false.
 367   continue
 360   continue
       return
       end

       subroutine Gamu(n,nc,ifood,itemp,itime,iscore,z,indexfd,
     1          indexftptm,theta,sqdelta,sig2,amu,iseed)
c......generate amu givenz,theta,sqdelta
       implicit real*8 (a-h,o-z)
       real*8 amu(1:12,1:4,1:9),theta(1:12,1:4),sqdelta(1:12)
       real*8 z(n)
       real*8 sumz(1:12,1:4,1:9),sig2(1:12)
       real*8 aleft,aright,tn
       real*8 amean,sigma
       integer ifood(n),itemp(n),itime(n),iscore(n)
       integer nc(4),indexftptm(1:12,1:4,1:9),indexfd(12)
       integer iseed
       logical la,lb
c......initialize
       la=.false.
       lb=.false.
       do 380 i=1,12
          if (indexfd(i) .eq. 0) goto 380
          do 385 j=1,4
             do 390 l=1,nc(j)
                sumz(i,j,l)=0.0d0
 390         continue
 385      continue
 380   continue
       do 400 i=1,n
         if (indexfd(ifood(i)) .eq. 0) goto 400
         sumz(ifood(i),itemp(i),itime(i))=
     1        sumz(ifood(i),itemp(i),itime(i))+z(i)
 400   continue
       do 410 i=1,12
         if (indexfd(i) .eq. 0) goto 410
         lb=.true.
         amean=(sumz(i,1,1)*sqdelta(i)+theta(i,1)*sig2(i))
     1         / (36.0d0*sqdelta(i)+sig2(i))
         sigma=dsqrt(sqdelta(i)*sig2(i)/(36.0d0*sqdelta(i)+sig2(i)))
         aleft=amu(i,1,3)
         aleft=max(aleft,amu(i,2,3))
         do 415 j=3,4
           aleft=max(aleft,amu(i,j,2))
 415     continue
         aleft=(aleft-amean)/sigma
         tn=YTUVN(aleft,aright,la,lb,iseed)
         lb=.false.
c........get amu(i,1,1)
         amu(i,1,1)=amean+sigma*tn
         do 420 j=1,4
            do 430 l=2,nc(j)
c.............get two different means and sigmas
              if (indexftptm(i,j,l) .eq. 1) goto 430
              amean=(sumz(i,j,l)*sqdelta(i)+sig2(i)*theta(i,j))/
     1              (36.0d0*sqdelta(i)+sig2(i))
c.............get the bounds of amu's               
c.............begin first row
              if ((j .eq. 1) .and. (l .eq. 3)) then
                aright=(amu(i,1,1)-amean)/sigma 
                aleft=(amu(i,2,3)-amean)/sigma
                aleft=max((amu(i,1,6)-amean)/sigma,aleft)
              endif
              if ((j .eq. 1) .and. (l .eq. 6)) then
                aright=(amu(i,1,3)-amean)/sigma
                aleft=(amu(i,2,6)-amean)/sigma
                aleft=max((amu(i,1,7)-amean)/sigma,aleft)
              endif
              if ( (j .eq. 1) .and. (l .lt. nc(1)) 
     1             .and. (l .gt. 6) ) then
                aright=(amu(i,j,l-1)-amean)/sigma
                aleft=(amu(i,j,l+1)-amean)/sigma
                aleft=max(aleft,(amu(i,j+1,l)-amean)/sigma)
              endif
              if ( (j .eq. 1) .and. (l .eq. nc(1)) ) then
                aright=(amu(i,j,l-1)-amean)/sigma
                aleft=(amu(i,j+1,l)-amean)/sigma
              endif
c.............end first row
c.............begin second row
              if ((j .eq. 2) .and. (l .eq. 3)) then
                aright=(amu(i,1,1)-amean)/sigma 
                aright=min(aright,(amu(i,1,3)-amean)/sigma)
                aleft=(amu(i,2,4)-amean)/sigma
                aleft=max((amu(i,3,3)-amean)/sigma,aleft)
              endif
              if ((j .eq. 2) .and. (l .eq. 4)) then
                aright=(amu(i,2,3)-amean)/sigma
                aleft=(amu(i,2,5)-amean)/sigma
                aleft=max((amu(i,3,4)-amean)/sigma,aleft)
              endif
              if ((j .eq. 2) .and. (l .eq. 5)) then
                aright=(amu(i,2,4)-amean)/sigma
                aleft=(amu(i,2,6)-amean)/sigma
                aleft=max((amu(i,3,5)-amean)/sigma,aleft)
              endif
              if ( (j .eq. 2) .and. (l .le. nc(3)) 
     1             .and. (l .gt. 5) ) then
                aright=(amu(i,j,l-1)-amean)/sigma
                aright=min(aright, (amu(i,j-1,l)-amean)/sigma)
                aleft=(amu(i,j,l+1)-amean)/sigma
                aleft=max(aleft,(amu(i,j+1,l)-amean)/sigma)
              endif
              if ( (j .eq. 2) .and. (l .lt. nc(2)) 
     1             .and. (l .gt. nc(3)) ) then
                aright=(amu(i,j,l-1)-amean)/sigma
                aright=min(aright, (amu(i,j-1,l)-amean)/sigma)
                aleft=(amu(i,j,l+1)-amean)/sigma
              endif
              if ( (j .eq. 2) .and. (l .eq. nc(2)) ) then
                aright=(amu(i,j,l-1)-amean)/sigma
                aright=min(aright, (amu(i,j-1,l)-amean)/sigma)
                la=.true.
              endif
c.............end second row
c.............begin third row
              if ((j .eq. 3) .and. (l .eq. 2)) then
                aright=(amu(i,1,1)-amean)/sigma
                aleft=(amu(i,3,3)-amean)/sigma
                aleft=max((amu(i,4,2)-amean)/sigma,aleft)
              endif
              if ( (j .eq. 3) .and. (l .le. nc(4))
     1             .and. (l .gt. 2) ) then
                aright=(amu(i,j,l-1)-amean)/sigma
                aright=min(aright, (amu(i,j-1,l)-amean)/sigma)
                aleft=(amu(i,j,l+1)-amean)/sigma
                aleft=max(aleft,(amu(i,j+1,l)-amean)/sigma)
              endif
              if ( (j .eq. 3) .and. (l .lt. nc(3))
     1             .and. (l .gt. nc(4)) ) then
                aright=(amu(i,j,l-1)-amean)/sigma
                aright=min(aright, (amu(i,j-1,l)-amean)/sigma)
                aleft=(amu(i,j,l+1)-amean)/sigma
              endif
              if ( (j .eq. 3) .and. (l .eq. nc(3)) ) then
                aright=(amu(i,j,l-1)-amean)/sigma
                aright=min(aright, (amu(i,j-1,l)-amean)/sigma)
                la=.true.
              endif
c.............end third row
c.............begin forth row
              if ((j .eq. 4) .and. (l .eq. 2)) then
                aright=(amu(i,1,1)-amean)/sigma
                aright=min(aright,(amu(i,j-1,2)-amean)/sigma)
                aleft=(amu(i,4,3)-amean)/sigma
              endif
              if ( (j .eq. 4) .and. (l .lt. nc(4))
     1             .and. (l .gt. 2) ) then
                aright=(amu(i,j,l-1)-amean)/sigma
                aright=min(aright, (amu(i,j-1,l)-amean)/sigma)
                aleft=(amu(i,j,l+1)-amean)/sigma
              endif
              if ( (j .eq. 4) .and. (l .eq. nc(4)) ) then
                aright=(amu(i,j,l-1)-amean)/sigma
                aright=min(aright, (amu(i,j-1,l)-amean)/sigma)
                la=.true.
              endif
c.............end forth row
c........get amu(i,j,l)
              tn=YTUVN(aleft,aright,la,lb,iseed)
              la=.false.
              lb=.false.
              amu(i,j,l)=amean+sigma*tn
 430        continue
 420     continue
 410   continue
       return
       end
