*
* $Id$
*
* $Log$
* Revision 1.1  1996/01/16 17:08:04  mclareni
* Initial revision
*
*
#include "hbook/pilot.h"
*CMZ :  4.17/02 18/12/92  14.23.02  by  John Allison
*-- Author :
      SUBROUTINE HQPOIS (SIGCOV, DDERIV, DJN, DSIGA, DDSIGA, N, CHISQ,
     +ALOGLI, IHQERR)
      INTEGER N
#if defined(CERNLIB_DOUBLE)
      DOUBLE PRECISION SIGCOV (N, N)
      DOUBLE PRECISION DDERIV (N, N)
      DOUBLE PRECISION DJN (N), DSIGA (N), DDSIGA (N)
#endif
#if !defined(CERNLIB_DOUBLE)
      REAL SIGCOV (N, N)
      REAL DDERIV (N, N)
      REAL DJN (N), DSIGA (N), DDSIGA (N)
#endif
      INTEGER IHQERR
      REAL CHISQ, ALOGLI
* Perform unconstrained Poisson (binned) likelihood fit to bin contents.
* IHQERR = 0 if all's OK.
* CHISQ is calculated by comparing the maximum likelihood solution to the trial
*   histogram data; hence it is not the minimum chi-squared.
* ALOGLI is log likelihood.
 
#include "hbook/hcqcom.inc"
#include "hbook/hcqcor.inc"
#include "hbook/hcbook.inc"
#include "hbook/hcunit.inc"
 
      CHARACTER*80 CHQMES
      LOGICAL PPRINT
      LOGICAL EVALDD, DESPRT
      INTEGER J, K, L, LB, LB0, NWW
      INTEGER IX, IY, IZ
      INTEGER NSECS, ITER, ITERDD, ITRMX0, ITERMX, IDDINT
      REAL OPS, UFACT
      REAL V (NDMAX)
      REAL X, Y, Z, FF, FFMIN, STEP, STEPMI
      EQUIVALENCE (X, V (1)), (Y, V(2)), (Z, V(3))
#if defined(CERNLIB_DOUBLE)
      PARAMETER (NWW = 2)
      DOUBLE PRECISION DCHISQ, DLOGLI, DLOGL0, DCHEK1, DCHEK2, DN
      DOUBLE PRECISION VEXPDN
      DOUBLE PRECISION HQDN, HQDJN
      DOUBLE PRECISION SS
#endif
#if !defined(CERNLIB_DOUBLE)
      PARAMETER (NWW = 1)
      REAL DCHISQ, DLOGLI, DLOGL0, DCHEK1, DCHEK2, DN
      REAL VEXPDN
      REAL HQDN, HQDJN
      REAL SS
#endif
 
      PARAMETER (ITRMX0 = 50, IDDINT = 3)
      PARAMETER (UFACT = 110000.)
      PARAMETER (DCHEK1 = 1.D-3)
      PARAMETER (DCHEK2 = 10. * DCHEK1)
      PARAMETER (STEPMI = 1. / 32.)
 
 
 
      CHISQ  = 0.
      ALOGLI = 0.
      IHQERR = 0
 
* Set function flag for exponential function.
      IMQFUN = 2
 
* Calculate constant factor.
      VCONST = VOLBIN
 
      OPS = FLOAT (NSIG ** 2) * FLOAT (NBTOT)
      PPRINT = VPRINT
      IF (PPRINT) THEN
         WRITE (LOUT, 10000) IDMQ
         IF (NDIM .EQ. 1) THEN
            WRITE (LOUT, 10100) NSIG, NBTOT
         ELSE IF (NDIM .EQ. 2) THEN
            WRITE (LOUT, 10200) NSIG, NBTOT, NX, NY
         ELSE IF (NDIM .EQ. 3) THEN
            WRITE (LOUT, 10300) NSIG, NBTOT, NX, NY, NZ
         ELSE
            GO TO 430
         END IF
         NSECS = OPS / UFACT
         WRITE (LOUT, 10400) NSECS, 10 * NSECS
         WRITE (LOUT, 10500) ITRMX0
      END IF
 
* Estimate starting values.
* First a least squares fit to linear sum of multiquadrics to get function
*   values (actually only needed for empty bins).
      IF (PPRINT) WRITE (LOUT,10600)
      CALL HQLSQU (SIGCOV, DJN, N, CHISQ, IHQERR)
      IF (IHQERR .NE. 0) GO TO 470
      IF (PPRINT) WRITE (LOUT, 10700) CHISQ
* Rebin.
      CALL HQBIN (NBINS, 1, IHQERR)
      IF (IHQERR .NE. 0) GO TO 470
* Unconstrained, linear least squares fit to log of (contents/VCONST) for
*   starting values.  This is an ad hoc method of getting starting values.
*   Some of the parameters were tuned by experience, and do not have much
*   foundation.  E.g., FFMIN is the minimum value substituted for empty bins
*   (to avoid log (0) = -INF) - otherwise a value from the previous step is
*   used.  Also I tried various algorithms for the variance of empty bins but
*   leaving it at 1 seems OK.
      IF (PPRINT) WRITE (LOUT, 10800)
      FFMIN  = 1. / NMQFUL
      IF (NDIM .EQ. 1) THEN
         DO 10    IX = 1, NX
            X = (IX - 0.5) / NX
            L = IX
            IF (Q (L1H + L) .EQ. 0.) THEN
               FF = HQDN (V)
               FF = MAX (FFMIN, FF)
               Q (L1H + L) = LOG (FF / VCONST)
*              Q (L1V + L) is already 1., as required.
            ELSE
               Q (L1V + L) = Q (L1V + L) / (Q (L1H + L)) ** 2
               Q (L1H + L) = LOG (Q (L1H + L) / VCONST)
            END IF
   10    CONTINUE
      ELSE IF (NDIM .EQ. 2) THEN
         DO 30    IX = 1, NX
            X = (IX - 0.5) / NX
            DO 20    IY = 1, NY
               Y = (IY - 0.5) / NY
               L = (IY - 1) * NX + IX
               IF (Q (L2H + L) .EQ. 0.) THEN
                  FF = HQDN (V)
                  FF = MAX (FFMIN, FF)
                  FF = MIN (1., FF)
                  Q (L2H + L) = LOG (FF / VCONST)
*                 Q (L2V + L) is already 1., as required.
               ELSE
                  Q (L2V + L) = Q (L2V + L) / (Q (L2H + L)) ** 2
                  Q (L2H + L) = LOG (Q (L2H + L) / VCONST)
               END IF
   20       CONTINUE
   30    CONTINUE
      ELSE IF (NDIM .EQ. 3) THEN
         DO 60    IX = 1, NX
            X = (IX - 0.5) / NX
            DO 50    IY = 1, NY
               Y = (IY - 0.5) / NY
               DO 40    IZ = 1, NZ
                  Z = (IZ - 0.5) / NZ
                  L = (IZ - 1) * NX * NY + (IY - 1) * NX + IX
                  IF (Q (L3H + L) .EQ. 0.) THEN
                     FF = HQDN (V)
                     FF = MAX (FFMIN, FF)
                     Q (L3H + L) = LOG (FF / VCONST)
*                    Q (L3V + L) is already 1., as required.
                  ELSE
                     Q (L3V + L) = Q (L3V + L) / (Q (L3H + L)) ** 2
                     Q (L3H + L) = LOG (Q (L3H + L) / VCONST)
                  END IF
   40          CONTINUE
   50       CONTINUE
   60    CONTINUE
      ELSE
         GO TO 430
      END IF
      CALL HQLSQU (SIGCOV, DJN, N, CHISQ, IHQERR)
      IF (IHQERR .NE. 0) GO TO 470
      IF (PPRINT) WRITE (LOUT, 10900) CHISQ
* Rebin.
      CALL HQBIN (NBINS, 1, IHQERR)
      IF (IHQERR .NE. 0) GO TO 470
 
* OK, let's go.
      IF (PPRINT) WRITE (LOUT,11000)
 
      IF (LHQDJN .NE. 0) THEN
* Fill working space with phi_j for each bin.
         IF (PPRINT) WRITE (LOUT, 11100)
         IF (NDIM .EQ. 1) THEN
            DO 80    IX = 1, NX
               X = (IX - 0.5) / NX
               L = IX
               LB0 = NSIG * (L - 1)
               DO 70    K = 1, NSIG
                  LB = NWW * (LB0 + K - 1) + 1
                  SS = HQDJN (V, K)
                  CALL UCOPY (SS, Q (LHQDJN + LB), NWW)
   70          CONTINUE
   80       CONTINUE
         ELSE IF (NDIM .EQ. 2) THEN
            DO 110   IX = 1, NX
               X = (IX - 0.5) / NX
               DO 100   IY = 1, NY
                  Y = (IY - 0.5) / NY
                  L = (IY - 1) * NX + IX
                  LB0 = NSIG * (L - 1)
                  DO 90    K = 1, NSIG
                     LB = NWW * (LB0 + K - 1) + 1
                     SS = HQDJN (V, K)
                     CALL UCOPY (SS, Q (LHQDJN + LB), NWW)
   90             CONTINUE
  100          CONTINUE
  110       CONTINUE
         ELSE IF (NDIM .EQ. 3) THEN
            DO 150   IX = 1, NX
               X = (IX - 0.5) / NX
               DO 140   IY = 1, NY
                  Y = (IY - 0.5) / NY
                  DO 130   IZ = 1, NZ
                     Z = (IZ - 0.5) / NZ
                     L = (IZ - 1) * NX * NY + (IY - 1) * NX + IX
                     LB0 = NSIG * (L - 1)
                     DO 120   K = 1, NSIG
                        LB = NWW * (LB0 + K - 1) + 1
                        SS = HQDJN (V, K)
                        CALL UCOPY (SS, Q (LHQDJN + LB), NWW)
  120                CONTINUE
  130             CONTINUE
  140          CONTINUE
  150       CONTINUE
         ELSE
            GO TO 430
         END IF
      END IF
 
* Begin iterating.
      DESPRT = .FALSE.
      DLOGL0 = -1.D30
      STEP   = 1.
      ITER   = 0
      ITERMX = ITRMX0
      ITERDD = IDDINT
      GO TO 340
 
* Start of iteration loop
  160 CONTINUE
      ITER   = ITER   + 1
      ITERDD = ITERDD + 1
      IF (ITER .GT. ITERMX) THEN
         WRITE (LOUT, 11200)
         GO TO 420
      END IF
  170 CONTINUE
      EVALDD = ITERDD .GE. IDDINT .OR. STEP .LT. 1.
      IF (EVALDD) ITERDD = 0
 
* Evaluate the matrices A and B and solve A * SIGA = B.
      CALL UZERO (DDSIGA, 1, NWW * NSIG)
      IF (EVALDD) CALL UZERO (DDERIV, 1, NWW * NSIG ** 2)
* Calculate DDSIGA and DDERIV lower triangle.
      IF (NDIM .EQ. 1) THEN
         DO 200   IX = 1, NX
            X = (IX - 0.5) / NX
            L = IX
            VEXPDN = VCONST * EXP (HQDN (V))
            LB0 = NSIG * (L - 1)
            DO 190   K = 1, NSIG
               IF (LHQDJN .NE. 0) THEN
                  LB = NWW * (LB0 + K - 1) + 1
                  CALL UCOPY (Q (LHQDJN + LB), SS, NWW)
               ELSE
                  SS = HQDJN (V, K)
               END IF
               DDSIGA (K) = DDSIGA (K) + (Q (L1H + L) - VEXPDN) * SS
               IF (EVALDD) THEN
                  DJN (K) = SS
                  DO 180   J = 1, K
                     DDERIV (J, K) = DDERIV (J, K) + VEXPDN * DJN (K) *
     +               DJN (J)
  180             CONTINUE
               END IF
  190       CONTINUE
  200    CONTINUE
      ELSE IF (NDIM .EQ. 2) THEN
         DO 240   IX = 1, NX
            X = (IX - 0.5) / NX
            DO 230   IY = 1, NY
               Y = (IY - 0.5) / NY
               L = (IY - 1) * NX + IX
               VEXPDN = VCONST * EXP (HQDN (V))
               LB0 = NSIG * (L - 1)
               DO 220   K = 1, NSIG
                  IF (LHQDJN .NE. 0) THEN
                     LB = NWW * (LB0 + K - 1) + 1
                     CALL UCOPY (Q (LHQDJN + LB), SS, NWW)
                  ELSE
                     SS = HQDJN (V, K)
                  END IF
                  DDSIGA (K) = DDSIGA (K) + (Q (L2H + L) - VEXPDN) * SS
                  IF (EVALDD) THEN
                     DJN (K) = SS
                     DO 210   J = 1, K
                        DDERIV (J, K) = DDERIV (J, K) + VEXPDN * DJN
     +                  (K) * DJN (J)
  210                CONTINUE
                  END IF
  220          CONTINUE
  230       CONTINUE
  240    CONTINUE
      ELSE IF (NDIM .EQ. 3) THEN
         DO 290   IX = 1, NX
            X = (IX - 0.5) / NX
            DO 280   IY = 1, NY
               Y = (IY - 0.5) / NY
               DO 270   IZ = 1, NZ
                  Z = (IZ - 0.5) / NZ
                  L = (IZ - 1) * NX * NY + (IY - 1) * NX + IX
                  VEXPDN = VCONST * EXP (HQDN (V))
                  LB0 = NSIG * (L - 1)
                  DO 260   K = 1, NSIG
                     IF (LHQDJN .NE. 0) THEN
                        LB = NWW * (LB0 + K - 1) + 1
                        CALL UCOPY (Q (LHQDJN + LB), SS, NWW)
                     ELSE
                        SS = HQDJN (V, K)
                     END IF
                     DDSIGA (K) = DDSIGA (K) + (Q (L3H + L) - VEXPDN) *
     +               SS
                     IF (EVALDD) THEN
                        DJN (K) = SS
                        DO 250   J = 1, K
                           DDERIV (J, K) = DDERIV (J, K) + VEXPDN *
     +                     DJN (K) * DJN (J)
  250                   CONTINUE
                     END IF
  260             CONTINUE
  270          CONTINUE
  280       CONTINUE
  290    CONTINUE
      ELSE
         GO TO 430
      END IF
* Fill DDERIV upper triangle.
      IF (EVALDD) THEN
         DO 310   K = 1, NSIG
            DO 300   J = 1, K - 1
               DDERIV (K, J) = DDERIV (J, K)
  300       CONTINUE
  310    CONTINUE
      END IF
* Copy to DSIGA and SIGCOV ready for inversion.
      CALL UCOPY (DDSIGA, DSIGA, NWW * NSIG)
      CALL UCOPY (DDERIV, SIGCOV, NWW * NSIG ** 2)
 
* Call CERN library routine DEQINV (F010) for simultaneous equations to
*   solve A * DSIGA = B.
#if defined(CERNLIB_DOUBLE)
      CALL DEQINV (NSIG, SIGCOV, NSIG, DJN, IHQERR, 1, DSIGA)
#endif
#if !defined(CERNLIB_DOUBLE)
      CALL REQINV (NSIG, SIGCOV, NSIG, DJN, IHQERR, 1, DSIGA)
#endif
      IF (IHQERR .NE. 0) GO TO 440
 
* Update SIGA.
  320 CONTINUE
      DO 330   J = 1, NSIG
         SIGA (J) = SIGA (J) + STEP * DSIGA (J)
  330 CONTINUE
 
* Evaluate log likelihood for convergence testing, and chi-squared.
  340 CONTINUE
      DLOGLI = 0.
      DCHISQ = 0.
      IF (NDIM .EQ. 1) THEN
         DO 350   IX = 1, NX
            X = (IX - 0.5) / NX
            L = IX
            DN = HQDN (V)
            VEXPDN = VCONST * EXP (DN)
            DLOGLI = DLOGLI + Q (L1H + L) * DN - VEXPDN
            DCHISQ = DCHISQ + (Q (L1H + L) - VEXPDN) ** 2 /
     +      Q (L1V + L)
  350    CONTINUE
      ELSE IF (NDIM .EQ. 2) THEN
         DO 370   IX = 1, NX
            X = (IX - 0.5) / NX
            DO 360   IY = 1, NY
               Y = (IY - 0.5) / NY
               L = (IY - 1) * NX + IX
               DN = HQDN (V)
               VEXPDN = VCONST * EXP (DN)
               DLOGLI = DLOGLI + Q (L2H + L) * DN - VEXPDN
               DCHISQ = DCHISQ + (Q (L2H + L) - VEXPDN) ** 2 /
     +         Q (L2V + L)
  360       CONTINUE
  370    CONTINUE
      ELSE IF (NDIM .EQ. 3) THEN
         DO 400   IX = 1, NX
            X = (IX - 0.5) / NX
            DO 390   IY = 1, NY
               Y = (IY - 0.5) / NY
               DO 380   IZ = 1, NZ
                  Z = (IZ - 0.5) / NZ
                  L = (IZ - 1) * NX * NY + (IY - 1) * NX + IX
                  DN = HQDN (V)
                  VEXPDN = VCONST * EXP (DN)
                  DLOGLI = DLOGLI + Q (L3H + L) * DN - VEXPDN
                  DCHISQ = DCHISQ + (Q (L3H + L) - VEXPDN) ** 2 /
     +            Q (L3V + L)
  380          CONTINUE
  390       CONTINUE
  400    CONTINUE
      ELSE
         GO TO 430
      END IF
      IF (PPRINT) THEN
         IF (ITER .EQ. 0) WRITE (LOUT, 11300)
         WRITE (LOUT, 11400) ITER, DLOGLI, DCHISQ
      END IF
 
* Check convergence.
      IF (DLOGLI - DLOGL0 .GT. DCHEK1) THEN
         DLOGL0 = DLOGLI
         IF (STEP .LT. 1.) STEP = 2. * STEP
         GO TO 160
      ELSE IF (DLOGLI - DLOGL0 .LT. -DCHEK2) THEN
* Decreasing log liklelihood - back up...
         DO 410   J = 1, NSIG
            SIGA (J) = SIGA (J) - STEP * DSIGA (J)
  410    CONTINUE
         IF (STEP .GE. STEPMI) THEN
* ...and try smaller step, re-evaluating 2nd derivative matrix next time...
            STEP = 0.5 * STEP
            GO TO 320
         ELSE
            IF (DESPRT) THEN
               GO TO 450
            ELSE
* ...or, indesperation, re-evaluate 2nd derivative now.
               DESPRT = .TRUE.
               GO TO 170
            END IF
         END IF
      END IF
 
* Set output variables.
  420 CONTINUE
      CHISQ  = DCHISQ
      ALOGLI = DLOGLI
 
      GO TO 470
 
  430 CONTINUE
      WRITE (CHQMES, '(I3, '' dimensions not programmed yet.'')') NDIM
      IHQERR = 20
      GO TO 460
  440 CONTINUE
      WRITE (CHQMES, '(''Error'', I5, '' in DEQINV'')') IHQERR
      GO TO 460
  450 CONTINUE
      WRITE (CHQMES, '(''Likelihood decreasing in spite of rescue'',
     +'' operations!'')')
      IHQERR = 40
      GO TO 460
 
  460 CONTINUE
      CALL HBUG (CHQMES, 'HQPOIS', IDMQ)
 
  470 CONTINUE
 
10000 FORMAT (1X, 'Multiquadric smoothing of histogram/ntuple', I7)
10100 FORMAT (I6, ' parameters for', I7, ' bins.')
10200 FORMAT (I6, ' parameters for', I7, ' bins (', I3, ' x', I3, ').')
10300 FORMAT (I6, ' parameters for', I7, ' bins (', I3, ' x', I3, ' x',
     +I3, ').')
10400 FORMAT (3X, 'This will take anything from', I10, ' to', I10,
     +' CERN-unit-secs.')
10500 FORMAT (3X, 'We allow, initially,', I4, ' iterations.')
10600 FORMAT (1X, 'First a least squares fit to a sum of multiquadrics'/
     +3X, '(to get function values for next step).')
10700 FORMAT (3X, 'Chi-squared was:', G12.5)
10800 FORMAT (
     +1X, 'Next a least squares fit of log of contents to a sum of '/
     +3X, 'multiquadrics (with help of function values derived above).')
10900 FORMAT (3X, 'Chi-squared was:', G12.5)
11000 FORMAT (
     +1X, 'Now an unconstrained Poisson likelihood fit of binned'/
     +3X, 'data to exponential of sum of multiquadrics.')
11100 FORMAT (3X, '(Storing function values to speed subsequent ',
     +'calculation.)')
11200 FORMAT (1X, 'Not converged - result accepted nevertheless.')
11300 FORMAT (1X,
     +'HQPOIS: Iteration  log likelihood         chi-squared')
11400 FORMAT (I16, 2X, 2G20.13)
 
      END
