' ******************************************************************
' Woolf test
' ******************************************************************

' ------------------------------------------------------------------
' Error codes
' ------------------------------------------------------------------

#DEFINE FOk 0
' No error

#DEFINE FSing -2
' Singularity

#DEFINE MatErrDim -3
' Non-compatible dimensions

' ------------------------------------------------------------------
' Global variable
' ------------------------------------------------------------------

COMMON SHARED ErrCode AS INTEGER
' Error code from the latest function evaluation

' ******************************************************************

SUB Woolf_Conform(      Obs()   AS INTEGER, _
                        Calc()  AS DOUBLE,  _
                        N_estim AS INTEGER, _
                  BYREF G       AS DOUBLE,  _
                  BYREF DoF     AS INTEGER)

  DIM AS INTEGER Lb, Ub, I

  Lb = LBOUND(Obs) : Ub = UBOUND(Obs)

  IF LBOUND(Calc) <> Lb OR UBOUND(Calc) <> Ub THEN
    ErrCode = MatErrDim
    EXIT SUB
  END IF

  FOR I = Lb TO Ub
    IF Obs(I) <= 0 OR Calc(I) <= 0 THEN
      ErrCode = FSing
      EXIT SUB
    END IF
  NEXT I

  G = 0
  ErrCode = FOk

  FOR I = Lb TO Ub
    G = G + Obs(I) * LOG(Obs(I) / Calc(I))
  NEXT I

  G = 2 * G
  DoF = Ub - Lb - N_estim
END SUB

SUB Woolf_Indep(      Obs() AS INTEGER, _
                BYREF G     AS DOUBLE,  _
                BYREF DoF   AS INTEGER)

  DIM AS INTEGER Lb1, Ub1, Lb2, Ub2, I, J

  Lb1 = LBOUND(Obs, 1) : Ub1 = UBOUND(Obs, 1)
  Lb2 = LBOUND(Obs, 2) : Ub2 = UBOUND(Obs, 2)

  FOR I = Lb1 TO Ub1
    FOR J = Lb2 TO Ub2
      IF Obs(I, J) <= 0 THEN
        ErrCode = FSing
        EXIT SUB
      END IF
    NEXT J
  NEXT I

  DIM AS INTEGER SumLin(Lb1 TO Ub1), SumCol(Lb2 TO Ub2), Sum

  DIM AS DOUBLE  Prob, Calc

  ErrCode = FOk

  FOR I = Lb1 TO Ub1
    SumLin(I) = 0
    FOR J = Lb2 TO Ub2
      SumLin(I) = SumLin(I) + Obs(I, J)
    NEXT J
    Sum = Sum + SumLin(I)
  NEXT I

  FOR J = Lb2 TO Ub2
    SumCol(J) = 0
    FOR I = Lb1 TO Ub1
      SumCol(J) = SumCol(J) + Obs(I, J)
    NEXT I
  NEXT J

  G = 0
  FOR I = Lb1 TO Ub1
    Prob = SumLin(I) / Sum
    FOR J = Lb2 TO Ub2
      Calc = SumCol(J) * Prob
      G = G + Obs(I, J) * LOG(Obs(I, J) / Calc)
    NEXT J
  NEXT I

  G = 2 * G
  DoF = (Ub2 - Lb2) * (Ub1 - Lb1)
END SUB
