Failure when first parameter has analytic derivative and next does not.

[rdnarayan]

It seems there may be a bookkeeping issue indexing the parameter derivatives when a parameter with user-supplied analytic derivatives is followed by one with numerical derivatives. This issue was observed after implementing the fix from the previous issue (#3), so there's a chance it could be due to that change. It seems the following change to line 1294 of MPFit.cs fixes the problem:

if (dside != null && dsidei == 3)
{
  ij += m; // bug fix
  continue;
}
//if (dside != null && dsidei == 3) continue; // replaced

The issue can be demonstrated by adding the following test to TestMPFit.cs:

   private static void TestGaussianWithDerivs()
   {
        double[] x = {-1.7237128E+00,1.8712276E+00,-9.6608055E-01,
    -2.8394297E-01,1.3416969E+00,1.3757038E+00,
    -1.3703436E+00,4.2581975E-02,-1.4970151E-01,
    8.2065094E-01};
        double[] y = {-4.4494256E-02,8.7324673E-01,7.4443483E-01,
    4.7631559E+00,1.7187297E-01,1.1639182E-01,
    1.5646480E+00,5.2322268E+00,4.2543168E+00,
    6.2792623E-01};
        double[] ey = new double[10];
        double[] p = { 0.0, 1.0, 1.0, 1.0 };       /\* Initial conditions _/
        double[] pactual = { 0.0, 4.70, 0.0, 0.5 };/_ Actual values used to make data_/
                                                   //double[] perror = new double[4];              /_ Returned parameter errors _/
        mp_par[] pars = new mp_par[4] /_ Parameter constraints */
                            {
                                new mp_par() { side = 3, deriv_debug = 0 },
                                new mp_par() { side = 1, deriv_debug = 0 },
                                new mp_par() { side = 3, deriv_debug = 0 },
                                new mp_par() { side = 1, deriv_debug = 0 }
                                //new mp_par() { deriv_debug = 1 },
                                //new mp_par() { deriv_debug = 1 },
                                //new mp_par() { deriv_debug = 1 },
                                //new mp_par() { deriv_debug = 1 }
                            };

    mp_result result = new mp_result(4);
    //result.xerror = perror;

    /* No constraints */

    for (uint i = 0; i < 10; i++) ey[i] = 0.5;

    CustomUserVariable v = new CustomUserVariable() { X = x, Y = y, Ey = ey };

    /* Call fitting function for 10 data points and 4 parameters (no
       parameters fixed) */

    var logger = new System.IO.StringWriter();
    int status = MPFit.Solve(GaussianFuncAndDerivs, 10, 4, p, pars, null, v, ref result, logger);

    Console.Write("*** TestGaussFitWithDerivs status = {0}\n", status);
    PrintResult(p, pactual, result);
    Console.WriteLine(logger.ToString());
}

private static int GaussianFuncAndDerivs(double[] p, double[] dy, IList<double>[] dvec, object vars)
{
    CustomUserVariable v = (CustomUserVariable)vars;
    double[] x, y, ey;
    double sig2;

    x = v.X;
    y = v.Y;
    ey = v.Ey;

    sig2 = p[3] * p[3];
    //var dvec1 = new double[dy.Length];
    //var residuals = new double[dy.Length];

    for (int i = 0; i < dy.Length; i++)
    {
        double xc = x[i] - p[2];
        double exp = Math.Exp(-0.5 * xc * xc / sig2);
        //residuals[i] = (y[i] - p[1] * exp - p[0]) / ey[i];
        dy[i] = (y[i] - p[1] * exp - p[0]) / ey[i];

        // NOTE: it would make sense to store the first 2 derivatives in vars since they don't change.
        if (dvec != null)
        {
            if (dvec[0] != null)
            {
                dvec[0][i] = -1.0 / ey[i];
            }
            if (dvec[1] != null)
            {
                //dvec1[i] = -exp / ey[i];
                dvec[1][i] = -exp / ey[i];
            }
            if (dvec[2] != null)
            {
                dvec[2][i] = -p[1] * xc * exp / (ey[i] * sig2);
            }
            if (dvec[3] != null)
            {
                dvec[3][i] = -p[1] * xc * xc * exp / (ey[i] * p[3] * sig2);
            }
        }
    }

    // Array assignment rather than element-wise causes failure.
    //dy = residuals;

    // Array mismatch exception
    //if (dvec != null)
    //{
    //  if (dvec[1] != null)
    //  {
    //      dvec[1] = dvec1;
    //  }
    //}
    return 0;
}