*
* $Id$
*
* $Log$
* Revision 1.1  1996/01/16 17:08:04  mclareni
* Initial revision
*
*
#include "hbook/pilot.h"
*CMZ :  4.17/13 14/05/93  09.08.25  by  John Allison
*-- Author :
      SUBROUTINE HQSIG2 (NX0, NY0, VOLB0, IEDGE, IHQERR)
      INTEGER NX0, NY0, IEDGE, IHQERR
      REAL VOLB0
* Finds significant points which will become centres of the radial basis
*   functions for 2-D histograms.
*
* Histogram contents and variances are in Q (L2H +1,...) and Q (L2V + 1,...).
***** Note: this routine overwrites histogram contents and variances.
*
* NX0    =  number of x bins.
* NY0    =  number of y bins.
* VOLB0  =  volume of bins (normalised coordinate space).
* IEDGE  =  0, normal.
*        =  2, zero extremities (i.e., constrained fit recommended).
*        = <0, special case - no significant points found this pass.
* IHQERR =  0 if all's OK.
*
* We want significant points where the rate of change of slope,
*   or in N-dimensions, the N-dimensional Laplacian, is statistically
*   significant.
*
* The coordinates of significant points - stored in HCQCOM - are in a
*   "normalised" system in which the plot variables are mapped to the
*    range 0 -> 1.
*
* This routine must be called with NDIMC (the current no. of dimensions) = 2.
 
#include "hbook/hcqcom.inc"
#include "hbook/hcbook.inc"
 
* Miscellaneous local variables.
      LOGICAL ANY, THIS, EDGED
      INTEGER NX00, NBTOT1, NBINS1 (NDMAX), NX1, NY1, NSLRMS
      EQUIVALENCE (NX1, NBINS1 (1)), (NY1, NBINS1 (2))
      INTEGER L, IDIM, NFACTR
      INTEGER L1M, L1P, L2M, L2P, L2PP, IX, IY, LJ
      INTEGER NDIMW2 (NDMAX)
      INTEGER IXE, IYE, IEDX (2), IEDY (2)
      REAL DXN1, DYN1, VOLB1, FACTOR
      REAL AL, SL, V, V1M, V1P, V2M, V2P, VL, DELTA
      REAL AS, VS, SS
      REAL DEDGE, EDGEN (2)
 
      DATA EDGEN /0., 1./
 
* Check NDIMC
      IF (NDIMC .NE. 2) GO TO 230
 
* Clear parameters.
      IEDGE  = 0
      IHQERR = 0
 
* Set edge parameter - edge points this many bins from edge.
      IF (NDIMC .LT. NDIM) THEN
         DEDGE = 0.5
      ELSE
         DEDGE = 0.
      END IF
 
* DXN1, etc. and NBTOT1 are current values of bin sizes and total no. of bins.
      NX1    = NX0
      NY1    = NY0
      DXN1   = 1. / NX1
      DYN1   = 1. / NY1
      NBTOT1 = NX1 * NY1
      VOLB1  = VOLB0
 
* Search histogram to find significant points.  Look for
*   statistically significant values of Laplacian, AL.  When
*   all have been found at this resolution, rebin and look again.
* The first time any have been found, add edge points.
      ANY = .FALSE.
      EDGED = .FALSE.
*   (Loop back here after rebinning.)
   10 CONTINUE
 
* Determine cut so that probability of purely statistical fluctuation is
*   small.  The probability of an x standard deviation fluctuation is
*   approximately exp(-(x/1.2)**2) (a crude fit to data in PDB).  So to
*   make the probability of a stat. fluctn. less than 1 / (no. of bins)
*   (say) (SENSIT is a user factor that adjusts the sensitivity of the cut
*   - normally SENSIT = 1.).
*   Relax by factor 0.7 for 2-D.
      SLCUT = 0.7 * 1.2 * SQRT (ALOG (FLOAT (NBTOT1))) / SENSIT
 
* Record whether any significant points have been found THIS loop.
      THIS = .FALSE.
 
* A convention here is that the pointers for the 4 points surrounding the point
*   of interest (whose pointer is L) are:
*   L1M, L1P for 1st dimension Minus and Plus one bin.
*   L2M, L2P for 2nd dimension Minus and Plus one bin.
      DELTA = MIN (DXN1, DYN1)
      SLRMS = 0.
      NSLRMS = 0
      DO 30    IY = 2, NY1 - 1
         YBIN = (IY - 0.5) * DYN1
         DO 20    IX = 2, NX1 - 1
            XBIN = (IX - 0.5) * DXN1
            L   = (IY - 1) * NX1 + IX
            L1M = L - 1
            L1P = L + 1
            L2M = L - NX1
            L2P = L + NX1
            AL = Q (L2H + L1M) + Q (L2H + L1P) +
     +      Q (L2H + L2M) + Q (L2H + L2P) - 4. * Q (L2H + L)
            V   = Q (L2V + L)
            V1M = Q (L2V + L1M)
            V1P = Q (L2V + L1P)
            V2M = Q (L2V + L2M)
            V2P = Q (L2V + L2P)
            IF (V   .EQ. 0.) V   = VMEAN2
            IF (V1M .EQ. 0.) V1M = VMEAN2
            IF (V1P .EQ. 0.) V1P = VMEAN2
            IF (V2M .EQ. 0.) V2M = VMEAN2
            IF (V2P .EQ. 0.) V2P = VMEAN2
            VL = V1M + V1P + V2M + V2P + 16. * V
            SL = ABS (AL / SQRT (VL))
            IF (SL .GT. 1.E-6) THEN
               SLRMS = SLRMS + SL ** 2
               NSLRMS = NSLRMS + 1
            END IF
            IF (SL .GE. SLCUT) THEN
               CALL HQSIGA (BINV, Q (L2H + L), VOLB1, DELTA, SL, IHQERR)
               IF (IHQERR .EQ. 0) THEN
                  ANY = .TRUE.
                  THIS = .TRUE.
               ELSE IF (IHQERR .LT. 0) THEN
                  IHQERR = 0
               ELSE
                  GO TO 240
               END IF
            END IF
   20    CONTINUE
   30 CONTINUE
      SLRMS = SQRT (SLRMS / NSLRMS)
 
* If this is the first time any significant points have been found...
      IF ((ANY .AND. .NOT. EDGED) .OR. (NDIMC .NE. NDIM)) THEN
         EDGED = .TRUE.
* Record bin size.
         IF (NDIM .EQ. NDIMC) THEN
            NXN = NX1
            NYN = NY1
         END IF
*************************************************************
* New strategy - add points if slope is significant.
* Low x...
         IX = 1
         XBIN = 0.
         DO 40 IY = 1, NY1
            YBIN = (IY - 0.5) * DYN1
            L   = (IY - 1) * NX1 + IX
            L1P = L + 1
            AS  = Q (L2H + L1P) - Q (L2H + L)
            V   = Q (L2V + L)
            V1P = Q (L2V + L1P)
            IF (V   .EQ. 0.) V   = VMEAN2
            IF (V1P .EQ. 0.) V1P = VMEAN2
            VS = V1P + V
            SS = ABS (AS / SQRT (VS))
            IF (SS .GE. SLCUT) THEN
               CALL HQSIGA (BINV, Q (L2H + L), VOLB1, DELTA,
     +         SS, IHQERR)
               IF (IHQERR .EQ. 0) THEN
                  ANY = .TRUE.
                  THIS = .TRUE.
               ELSE IF (IHQERR .LT. 0) THEN
                  IHQERR = 0
               ELSE
                  GO TO 240
               END IF
            END IF
   40    CONTINUE
* High x...
         IX = NX1
         XBIN = 1.
         DO 50 IY = 1, NY1
            YBIN = (IY - 0.5) * DYN1
            L   = (IY - 1) * NX1 + IX
            L1M = L - 1
            AS  = Q (L2H + L1M) - Q (L2H + L)
            V   = Q (L2V + L)
            V1M = Q (L2V + L1M)
            IF (V   .EQ. 0.) V   = VMEAN2
            IF (V1M .EQ. 0.) V1M = VMEAN2
            VS = V1M + V
            SS = ABS (AS / SQRT (VS))
            IF (SS .GE. SLCUT) THEN
               CALL HQSIGA (BINV, Q (L2H + L), VOLB1, DELTA,
     +         SS, IHQERR)
               IF (IHQERR .EQ. 0) THEN
                  ANY = .TRUE.
                  THIS = .TRUE.
               ELSE IF (IHQERR .LT. 0) THEN
                  IHQERR = 0
               ELSE
                  GO TO 240
               END IF
            END IF
   50    CONTINUE
* Low y...
         IY = 1
         YBIN = 0.
         DO 60 IX = 1, NX1
            XBIN = (IX - 0.5) * DXN1
            L   = (IY - 1) * NX1 + IX
            L2P = L + NX1
            AS  = Q (L2H + L2P) - Q (L2H + L)
            V   = Q (L2V + L)
            V2P = Q (L2V + L2P)
            IF (V   .EQ. 0.) V   = VMEAN2
            IF (V2P .EQ. 0.) V2P = VMEAN2
            VS = V2P + V
            SS = ABS (AS / SQRT (VS))
            IF (SS .GE. SLCUT) THEN
               CALL HQSIGA (BINV, Q (L2H + L), VOLB1, DELTA,
     +         SS, IHQERR)
               IF (IHQERR .EQ. 0) THEN
                  ANY = .TRUE.
                  THIS = .TRUE.
               ELSE IF (IHQERR .LT. 0) THEN
                  IHQERR = 0
               ELSE
                  GO TO 240
               END IF
            END IF
   60    CONTINUE
* High y...
         IY = NY1
         YBIN = 1.
         DO 70 IX = 1, NX1
            XBIN = (IX - 0.5) * DXN1
            L   = (IY - 1) * NX1 + IX
            L2M = L - NX1
            AS  = Q (L2H + L2M) - Q (L2H + L)
            V   = Q (L2V + L)
            V2M = Q (L2V + L2M)
            IF (V   .EQ. 0.) V   = VMEAN2
            IF (V2M .EQ. 0.) V2M = VMEAN2
            VS = V2M + V
            SS = ABS (AS / SQRT (VS))
            IF (SS .GE. SLCUT) THEN
               CALL HQSIGA (BINV, Q (L2H + L), VOLB1, DELTA,
     +         SS, IHQERR)
               IF (IHQERR .EQ. 0) THEN
                  ANY = .TRUE.
                  THIS = .TRUE.
               ELSE IF (IHQERR .LT. 0) THEN
                  IHQERR = 0
               ELSE
                  GO TO 240
               END IF
            END IF
   70    CONTINUE
* End of new strategy.
*************************************************************************
******************************************************************
* Note: the following has been modified by commenting out some statements
*   to suppress calls to other HQSIGn routines.  This also has the effect
*   of to forcing the adding of corner points.
******************************************************************
* Now we are going to present HQSIG1 with the faces of this 2-D histogram.
*   HQSIG1 will think it's a 1-D histogram.  When it comes to adding a point,
*   however, HQSIGA must know:
*     (a) the original dimensionality, NDIM,
*     (b) the current dimensionality, NDIMC,
*     (c) the coordinate(s) being treated, (NDIMWH (I), I = 1, NDIMC) and
*     (d) the value of the remaining coordinates,
*           (BINV (I), I = NDIMC + 1, NDIM), if any, destined for
*           (SIGV (NSIG, (NDIMWH(I))), I = NDIMC + 1, NDIM).
* Preserve current values of NDIMWH (NDIMC is assumed to be 2 at this point).
         DO 80    IDIM = 1, NDIMC
            NDIMW2 (IDIM) = NDIMWH (IDIM)
   80    CONTINUE
* Find sig. points along the 4 edges...
         NDIMC = 1
* x edges...
         DO 100   IX = 1, NX1, NX1 - 1
            IXE = IX / NX1 + 1
            DO 90    IY = 1, NY1
               L = (IY - 1) * NX1 + IX
               Q (L1H + IY) = Q (L2H + L)
               Q (L1V + IY) = Q (L2V + L)
   90       CONTINUE
            NDIMWH (1) = NDIMW2 (2)
            NDIMWH (2) = NDIMW2 (1)
            BINV (2) = EDGEN (IXE) -
     +      (2. * EDGEN (IXE) - 1.) * DEDGE * DXN1
            IEDX (IXE) = 0
**************            CALL HQSIG1 (NY1, VOLB1, IEDX (IXE), IHQERR)
            IF (IHQERR .NE. 0) GO TO 240
  100    CONTINUE
* y edges...
         DO 120   IY = 1, NY1, NY1 - 1
            IYE = IY / NY1 + 1
            DO 110   IX = 1, NX1
               L = (IY - 1) * NX1 + IX
               Q (L1H + IX) = Q (L2H + L)
               Q (L1V + IX) = Q (L2V + L)
  110       CONTINUE
            NDIMWH (1) = NDIMW2 (1)
            NDIMWH (2) = NDIMW2 (2)
            BINV (2) = EDGEN (IYE) -
     +      (2. * EDGEN (IYE) - 1.) * DEDGE * DYN1
            IEDY (IYE) = 0
***************            CALL HQSIG1 (NX1, VOLB1, IEDY (IYE), IHQERR)
            IF (IHQERR .NE. 0) GO TO 240
  120    CONTINUE
* Restore current values of NDIMWH.
         NDIMC = 2
         DO 130   IDIM = 1, NDIMC
            NDIMWH (IDIM) = NDIMW2 (IDIM)
  130    CONTINUE
* Add corner points if adjacent edges have points.
         DO 150   IX = 1, NX1, NX1 - 1
            IXE = IX / NX1 + 1
            DO 140   IY = 1, NY1, NY1 - 1
               IYE = IY / NY1 + 1
               IF (IEDX (IXE) .GE. 0 .AND. IEDY (IYE) .GE. 0) THEN
                  XBIN = EDGEN (IXE) -
     +            (2. * EDGEN (IXE) - 1.) * DEDGE * DXN1
                  YBIN = EDGEN (IYE) -
     +            (2. * EDGEN (IYE) - 1.) * DEDGE * DYN1
                  L = (IY - 1) * NX1 + IX
                  CALL HQSIGA (BINV, Q (L2H + L), VOLB1, SMALLD,
     +            SLCUT, IHQERR)
                  IF (IHQERR .EQ. 0) THEN
                     ANY = .TRUE.
                     THIS = .TRUE.
                  ELSE IF (IHQERR .LT. 0) THEN
                     IHQERR = 0
                  ELSE
                     GO TO 240
                  END IF
               END IF
  140       CONTINUE
  150    CONTINUE
* Are the extremities zero?
         DO 170   IX = 1, NX1, NX1 - 1
            DO 160   IY = 1, NY1, NY1 - 1
               L = (IY - 1) * NX1 + IX
               IF (Q (L2H + L) .NE. 0) GO TO 200
  160       CONTINUE
  170    CONTINUE
         DO 180   IXE = 1, 2
            IF (IEDX (IXE) .GE. 0) GO TO 200
  180    CONTINUE
         DO 190   IYE = 1, 2
            IF (IEDY (IYE) .GE. 0) GO TO 200
  190    CONTINUE
* If so, mark for zero asyptotic constrained fit.
         IEDGE = 2
  200    CONTINUE
      END IF
 
* If no points have been found, or points have been found this time round (i.e.,
*   cease searching if points have been found but not this time round).
**** (Not sure about the wisdom of this (why did I introduce it?) so disable for
****   now.)
*     IF (.NOT. ANY .OR. THIS) THEN
      IF (NTUPLE) THEN
* Rebin - the data is from an ntuple.  This simplifies things immensely.
* Double significance by increasing bin size by cube root of 2.
         FACTOR = 1.5
         NFACTR = 10. * FACTOR
         IF (NX1 .GE. NFACTR .AND. NY1 .GE. NFACTR) THEN
            NX1 = NX1 / FACTOR
            NY1 = NY1 / FACTOR
            NBTOT1 = NX1 * NY1
            DXN1 = 1. / NX1
            DYN1 = 1. / NY1
            VOLB1 = DXN1 * DYN1
* Get contents and errors.
            CALL HQBIN (NBINS1, 0, IHQERR)
            IF (IHQERR. NE. 0) GO TO 240
            GO TO 10
         END IF
      ELSE
* Rebin - double up.
***** NOTE: bin doubling leaves the odd bin over.  Not knowing what to do with
*****   it I simply leave it - possibly a significant point might be missed but
*****   the above code gave it at least one chance.  When it comes to choosing
*****   bins for finding significant points for N-tuples, choose 2**N bins.
         IF (NX1 .GE. 20 .AND. NY1 .GE. 20) THEN
            NX00 = NX1
            NX1 = NX1 / 2
            NY1 = NY1 / 2
            DXN1 = 2. * DXN1
            DYN1 = 2. * DYN1
            VOLB1 = 4. * VOLB1
            NBTOT1 = NX1 * NY1
            DO 220   IY = 1, NY1
               DO 210   IX = 1, NX1
                  L = 2 * (IY - 1) * NX00 + 2 * (IX - 1) + 1
                  L1P = L + 1
                  L2P = L + NX00
                  L2PP = L + NX00 + 1
                  LJ = (IY - 1) * NX1 + IX
                  Q (L2H + LJ) = Q (L2H + L) + Q (L2H + L1P) + Q (L2H +
     +            L2P) + Q (L2H + L2PP)
                  Q (L2V + LJ) = Q (L2V + L) + Q (L2V + L1P) + Q (L2V +
     +            L2P) + Q (L2V + L2PP)
  210          CONTINUE
  220       CONTINUE
            GO TO 10
         END IF
      END IF
*     END IF
 
* If no points have been found this pass, use IEDGE as a flag.
      IF (.NOT. ANY) IEDGE = -1
 
      GO TO 240
 
* Errors.
  230 CONTINUE
      CALL HBUG ('NDIMC not 2.', 'HQSIG2', IDMQ)
      IHQERR = 10
 
  240 CONTINUE
 
      END
