/* FUNCTION mrqmin */
/* Marquart-Levenworth Nonlinear Curve Fit adapted to PL/SQL         */

CREATE OR REPLACE FUNCTION mrqmin(pX IN IC50_ENGINE.NumTabType,
                                  pY IN IC50_ENGINE.NumTabType, 
                                  pSig IN IC50_ENGINE.NumTabType,
                                  pNData IN NUMBER,
                                  pA IN OUT IC50_ENGINE.NumTabType,
                                  pIa IN IC50_ENGINE.NumTabType,
                                  pMa IN NUMBER,
                                  pCovar IN OUT IC50_ENGINE.MultiDNumTabType,
                                  pAlfa IN OUT IC50_ENGINE.MultiDNumTabType,
                                  pChiSq IN OUT NUMBER,
                                  pALambda IN OUT NUMBER )
RETURN NUMBER

AS 
/* minimizes the chi square fit between a set of pNData data points (x,y),
   with standard deviations pSig of y at x, and a nonlinear function.
   the function has pMa parameters (= 4 for IC50) given as pA().
   pIa(i) indicates if the parameter is fitted (= 1) or fixed (= 0).
   pCovar is the covariance matrix.
   pAlfa is the linearized fitting matrix.
   pChiSq is the chi square.
   pALambda is the retry parameter.
   when pAlambda = 0, pCovar is the covariance matrix output, and
   pAlfa is the curvature matrix.
   fixed parameters have 0 covariance. */
 j      NUMBER;
 k      NUMBER;
 l      NUMBER;
 pMfit  NUMBER;
 vError NUMBER;
/* ochisq, atry, beta, da, oneda - are defined in PACK_IC50_ENGINE.SQL as
   gMrqOldChiSq NUMBER;
   gMrqAtry NumTabType;
   gMrqBeta NumTabType;
   gMrqDa NumTabType;
   gMrqOneda MultiDNumTabType; */
 
BEGIN

/* place this before IF (not like in book) since sets pMfit each time
   unlike C, which can keep it. */
  pMfit := 0;
  FOR j IN 1..pMa LOOP
    IF pIa(j) <> 0 THEN
      pMfit := pMfit + 1;
    END IF;
  END LOOP;
 
  IF pALambda < 0.0 THEN
    pALambda := 0.001;
    vError := mrqcof(pX, pY, pSig, pNData, pA, pIa, pMa,
              pAlfa,IC50_Engine.gMrqBeta,pChiSq);
    IC50_Engine.gMrqOldChiSq := pChiSq;
    FOR j IN 1..pMa LOOP
      IC50_Engine.gMrqAtry(j) := pA(j);
    END LOOP;
  END IF;
  
  FOR j IN 1..pMfit LOOP
    FOR k IN 1..pMfit LOOP
/* pCovar(j,k) = pAlfa(j,k) */
      pCovar((j - 1) * 4 + k) := pAlfa((j - 1) * 4 + k); 
    END LOOP;
/* pCovar(j,j) = pAlfa(j,j) + ... */
    pCovar((j - 1) * 4 + j) :=
      pAlfa((j - 1) * 4 + j) * (1.0 + pALambda);
/* oneda(j,1) = ... */
    IC50_Engine.gMrqOneDa((j - 1) * 4 + 1) := IC50_Engine.gMrqBeta(j);
  END LOOP;
  
  vError := gaussj(pCovar, pMfit, IC50_Engine.gMrqOneDa, 1);
  IF vError <> 0 THEN 
    RETURN IC50_ENGINE.sFail + 2;
  END IF;
  
  FOR j IN 1..pMfit LOOP
/* ... = oneda(j,1) */
    IC50_Engine.gMrqDa(j) := IC50_Engine.gMrqOneDa((j - 1) * 4 + 1);
  END LOOP;
  
  IF pALambda = 0.0 THEN 
    covsrt(pCovar, pMa, pIa, pMfit);
    RETURN 0;
  END IF;

  j := 0;
  FOR l IN 1..pMa LOOP
    IF pIa(l) <> 0 THEN
      j := j + 1;
      IC50_Engine.gMrqAtry(l) := pA(l) + IC50_Engine.gMrqDa(j);
    END IF;
  END LOOP;

/*  IF (IC50_Engine.gMrqAtry(3) > 0 AND                            */
/*      IC50_Engine.gMrqAtry(4) >= IC50_Engine.gMrqAtry(1) AND     */
/*      IC50_Engine.gMrqAtry(2) < 10 AND                           */
/*      IC50_Engine.gMrqAtry(2) > -10 AND                          */
/*      (IC50_Engine.gMrqAtry(2) < -0.001 OR                       */
/*       IC50_Engine.gMrqAtry(2) > 0.001)) THEN                    */

  vError := mrqcof(pX, pY, pSig, pNData, IC50_Engine.gMrqAtry, pIa, pMa,
    pCovar, IC50_Engine.gMrqDa, pChiSq);
  
  IF pChiSq < IC50_Engine.gMrqOldChiSq THEN
    pALambda := pALambda * 0.1;
    IC50_Engine.gMrqOldChiSq := pChiSq;
    FOR j IN 1..pMfit LOOP
      FOR k IN 1..pMfit LOOP
/* pAlfa(j,k) = pCovar(j,k) */
        pAlfa((j - 1) * 4 + k) := pCovar((j - 1) * 4 + k);
      END LOOP;
      IC50_Engine.gMrqBeta(j) := IC50_Engine.gMrqDa(j);
    END LOOP;
    FOR l IN 1..pMa LOOP
      pA(l) := IC50_Engine.gMrqAtry(l);
    END LOOP;
  ELSE
    pALambda := pALambda * 10.0;
    pChiSq := IC50_Engine.gMrqOldChiSq;
  END IF;

/*  ELSE                                             */
/*    RETURN 1;                                      */
/*  END IF;                                          */

 RETURN 0;
END;
/
SHOW ERRORS FUNCTION mrqmin