diff --git a/src/Numerics/Optimization/DelimitedArray.cs b/src/Numerics/Optimization/DelimitedArray.cs index 8de0cc30..ee2d5a66 100644 --- a/src/Numerics/Optimization/DelimitedArray.cs +++ b/src/Numerics/Optimization/DelimitedArray.cs @@ -61,7 +61,7 @@ namespace MathNet.Numerics.Optimization { int _offset; int _count; - T[] _array; + readonly T[] _array; public DelimitedArray(T[] array, int offset, int count) { diff --git a/src/Numerics/Optimization/MPFit.cs b/src/Numerics/Optimization/MPFit.cs index 024b6af0..7c3c9492 100644 --- a/src/Numerics/Optimization/MPFit.cs +++ b/src/Numerics/Optimization/MPFit.cs @@ -47,45 +47,41 @@ namespace MathNet.Numerics.Optimization { public static class MpFit { - public const string MPFIT_VERSION = "1.1"; - /* Error codes */ - public const int MP_ERR_INPUT = 0; /* General input parameter error */ - public const int MP_ERR_NAN = -16; /* User function produced non-finite values */ - public const int MP_ERR_FUNC = -17; /* No user function was supplied */ - public const int MP_ERR_NPOINTS = -18; /* No user data points were supplied */ - public const int MP_ERR_NFREE = -19; /* No free parameters */ - public const int MP_ERR_MEMORY = -20; /* Memory allocation error */ - public const int MP_ERR_INITBOUNDS = -21; /* Initial values inconsistent w constraints*/ - public const int MP_ERR_BOUNDS = -22; /* Initial constraints inconsistent */ - public const int MP_ERR_PARAM = -23; /* General input parameter error */ - public const int MP_ERR_DOF = -24; /* Not enough degrees of freedom */ + public const int MpErrNan = -16; /* User function produced non-finite values */ + public const int MpErrFunc = -17; /* No user function was supplied */ + public const int MpErrNpoints = -18; /* No user data points were supplied */ + public const int MpErrNfree = -19; /* No free parameters */ + public const int MpErrInitbounds = -21; /* Initial values inconsistent w constraints*/ + public const int MpErrBounds = -22; /* Initial constraints inconsistent */ + public const int MpErrParam = -23; /* General input parameter error */ + public const int MpErrDof = -24; /* Not enough degrees of freedom */ /* Potential success status codes */ - public const int MP_OK_CHI = 1; /* Convergence in chi-square value */ - public const int MP_OK_PAR = 2; /* Convergence in parameter value */ - public const int MP_OK_BOTH = 3; /* Both MP_OK_PAR and MP_OK_CHI hold */ - public const int MP_OK_DIR = 4; /* Convergence in orthogonality */ - public const int MP_MAXITER = 5; /* Maximum number of iterations reached */ - public const int MP_FTOL = 6; /* ftol is too small; no further improvement*/ - public const int MP_XTOL = 7; /* xtol is too small; no further improvement*/ - public const int MP_GTOL = 8; /* gtol is too small; no further improvement*/ + public const int MpConvergedChiSquared = 1; /* Convergence in chi-square value */ + public const int MpConvergedParameter = 2; /* Convergence in parameter value */ + public const int MpConvergedBoth = 3; /* Both MP_OK_PAR and MP_OK_CHI hold */ + public const int MpConvergedOrthogonality = 4; /* Convergence in orthogonality */ + public const int MpMaxIterations = 5; /* Maximum number of iterations reached */ + public const int MpFtol = 6; /* ftol is too small; no further improvement*/ + public const int MpXtol = 7; /* xtol is too small; no further improvement*/ + public const int MpGtol = 8; /* gtol is too small; no further improvement*/ #if FLOAT_PRECISION /* Float precision */ - public const float MP_MACHEP0 =1.19209e-07; - public const float MP_DWARF = 1.17549e-38; - public const float MP_GIANT = 3.40282e+38; - public const float MP_RDWARF = 1.3278686946331594e-018; - public const float MP_RGIANT = 1844673472786071600; + const float MP_MACHEP0 =1.19209e-07; + const float MP_DWARF = 1.17549e-38; + const float MP_GIANT = 3.40282e+38; + const float MP_RDWARF = 1.3278686946331594e-018; + const float MP_RGIANT = 1844673472786071600; #else /* Double precision numeric constants */ - public const double MP_MACHEP0 = 2.2204460e-16; - public const double MP_DWARF = 2.2250739e-308; - public const double MP_GIANT = 1.7976931e+308; - public const double MP_RDWARF = 1.8269129289596699331800430554921e-153; - public const double MP_RGIANT = 1.3407807799935081109978164571307e+153; + const double MP_MACHEP0 = 2.2204460e-16; + const double MP_DWARF = 2.2250739e-308; + const double MP_GIANT = 1.7976931e+308; + const double MP_RDWARF = 1.8269129289596699331800430554921e-153; + const double MP_RGIANT = 1.3407807799935081109978164571307e+153; #endif /* Expand for full description of Solve and lmdif functions @@ -378,11 +374,11 @@ namespace MathNet.Numerics.Optimization conf.stepfactor = 100.0; conf.nprint = 1; conf.epsfcn = MP_MACHEP0; - conf.maxiter = 200; - conf.douserscale = 0; - conf.maxfev = 0; + conf.MaxIterations = 200; + conf.DoUserScale = 0; + conf.MaxEvaluations = 0; conf.covtol = 1e-14; - conf.nofinitecheck = 0; + conf.NoFiniteCheck = 0; if (config != null) { @@ -393,11 +389,11 @@ namespace MathNet.Numerics.Optimization if (config.stepfactor > 0) conf.stepfactor = config.stepfactor; if (config.nprint >= 0) conf.nprint = config.nprint; if (config.epsfcn > 0) conf.epsfcn = config.epsfcn; - if (config.maxiter > 0) conf.maxiter = config.maxiter; - if (config.douserscale != 0) conf.douserscale = config.douserscale; + if (config.MaxIterations > 0) conf.MaxIterations = config.MaxIterations; + if (config.DoUserScale != 0) conf.DoUserScale = config.DoUserScale; if (config.covtol > 0) conf.covtol = config.covtol; - if (config.nofinitecheck > 0) conf.nofinitecheck = config.nofinitecheck; - conf.maxfev = config.maxfev; + if (config.NoFiniteCheck > 0) conf.NoFiniteCheck = config.NoFiniteCheck; + conf.MaxEvaluations = config.MaxEvaluations; } info = 0; @@ -407,17 +403,17 @@ namespace MathNet.Numerics.Optimization if (funct == null) { - return MP_ERR_FUNC; + return MpErrFunc; } if ((m <= 0) || (xall == null)) { - return MP_ERR_NPOINTS; + return MpErrNpoints; } if (npar <= 0) { - return MP_ERR_NFREE; + return MpErrNfree; } fnorm = -1.0; @@ -469,7 +465,7 @@ namespace MathNet.Numerics.Optimization } if (nfree == 0) { - info = MP_ERR_NFREE; + info = MpErrNfree; return info; } @@ -481,7 +477,7 @@ namespace MathNet.Numerics.Optimization (pars[i].limited[0] != 0 && (xall[i] < pars[i].limits[0])) || (pars[i].limited[1] != 0 && (xall[i] > pars[i].limits[1]))) { - info = MP_ERR_INITBOUNDS; + info = MpErrInitbounds; return info; } if ((pars[i].isFixed != 0) && @@ -489,7 +485,7 @@ namespace MathNet.Numerics.Optimization (pars[i].limited[1] != 0) && (pars[i].limits[0] >= pars[i].limits[1])) { - info = MP_ERR_BOUNDS; + info = MpErrBounds; return info; } } @@ -514,17 +510,17 @@ namespace MathNet.Numerics.Optimization /* Sanity checking on input configuration */ if ((npar <= 0) || (conf.ftol <= 0) || (conf.xtol <= 0) || - (conf.gtol <= 0) || (conf.maxiter < 0) || + (conf.gtol <= 0) || (conf.MaxIterations < 0) || (conf.stepfactor <= 0)) { - info = MP_ERR_PARAM; + info = MpErrParam; return info; } /* Ensure there are some degrees of freedom */ if (m < nfree) { - info = MP_ERR_DOF; + info = MpErrDof; return info; } @@ -635,7 +631,7 @@ namespace MathNet.Numerics.Optimization */ if (iter == 1) { - if (conf.douserscale == 0) + if (conf.DoUserScale == 0) { for (j = 0; j < nfree; j++) { @@ -698,7 +694,7 @@ namespace MathNet.Numerics.Optimization /* ( From this point on, only the square matrix, consisting of the triangle of R, is needed.) */ - if (conf.nofinitecheck != 0) + if (conf.NoFiniteCheck != 0) { /* Check for overflow. This should be a cheap test here since FJAC has been reduced to a (small) square matrix, and the test is @@ -719,7 +715,7 @@ namespace MathNet.Numerics.Optimization if (nonfinite != 0) { - info = MP_ERR_NAN; + info = MpErrNan; return info; } } @@ -752,14 +748,14 @@ namespace MathNet.Numerics.Optimization /* * test for convergence of the gradient norm. */ - if (gnorm <= conf.gtol) info = MP_OK_DIR; + if (gnorm <= conf.gtol) info = MpConvergedOrthogonality; if (info != 0) goto L300; - if (conf.maxiter == 0) goto L300; + if (conf.MaxIterations == 0) goto L300; /* * rescale if necessary. */ - if (conf.douserscale == 0) + if (conf.DoUserScale == 0) { for (j = 0; j < nfree; j++) { @@ -984,16 +980,16 @@ namespace MathNet.Numerics.Optimization if ((Math.Abs(actred) <= conf.ftol) && (prered <= conf.ftol) && (p5*ratio <= one)) { - info = MP_OK_CHI; + info = MpConvergedChiSquared; } if (delta <= conf.xtol*xnorm) { - info = MP_OK_PAR; + info = MpConvergedParameter; } if ((Math.Abs(actred) <= conf.ftol) && (prered <= conf.ftol) && (p5*ratio <= one) && (info == 2)) { - info = MP_OK_BOTH; + info = MpConvergedBoth; } if (info != 0) { @@ -1003,27 +999,27 @@ namespace MathNet.Numerics.Optimization /* * tests for termination and stringent tolerances. */ - if ((conf.maxfev > 0) && (nfev >= conf.maxfev)) + if ((conf.MaxEvaluations > 0) && (nfev >= conf.MaxEvaluations)) { /* Too many function evaluations */ - info = MP_MAXITER; + info = MpMaxIterations; } - if (iter >= conf.maxiter) + if (iter >= conf.MaxIterations) { /* Too many iterations */ - info = MP_MAXITER; + info = MpMaxIterations; } if ((Math.Abs(actred) <= MP_MACHEP0) && (prered <= MP_MACHEP0) && (p5*ratio <= one)) { - info = MP_FTOL; + info = MpFtol; } if (delta <= MP_MACHEP0*xnorm) { - info = MP_XTOL; + info = MpXtol; } if (gnorm <= MP_MACHEP0) { - info = MP_GTOL; + info = MpGtol; } if (info != 0) { @@ -1074,35 +1070,35 @@ namespace MathNet.Numerics.Optimization } /* Compute and return the covariance matrix and/or parameter errors */ - if (result != null && (result.covar != null || result.xerror != null)) + if (result != null && (result.FinalParameterCovarianceMatrix != null || result.FinalparameterUncertainties != null)) { mp_covar(nfree, fjac, ldfjac, ipvt, conf.covtol, wa2); - if (result.covar != null) + if (result.FinalParameterCovarianceMatrix != null) { /* Zero the destination covariance array */ - for (j = 0; j < (npar*npar); j++) result.covar[j] = 0; + for (j = 0; j < (npar*npar); j++) result.FinalParameterCovarianceMatrix[j] = 0; /* Transfer the covariance array */ for (j = 0; j < nfree; j++) { for (i = 0; i < nfree; i++) { - result.covar[ifree[j]*npar + ifree[i]] = fjac[j*ldfjac + i]; + result.FinalParameterCovarianceMatrix[ifree[j]*npar + ifree[i]] = fjac[j*ldfjac + i]; } } } - if (result.xerror != null) + if (result.FinalparameterUncertainties != null) { - for (j = 0; j < npar; j++) result.xerror[j] = 0; + for (j = 0; j < npar; j++) result.FinalparameterUncertainties[j] = 0; for (j = 0; j < nfree; j++) { double cc = fjac[j*ldfjac + j]; if (cc > 0) { - result.xerror[ifree[j]] = Math.Sqrt(cc); + result.FinalparameterUncertainties[ifree[j]] = Math.Sqrt(cc); } } } @@ -1110,22 +1106,21 @@ namespace MathNet.Numerics.Optimization if (result != null) { - result.version = MPFIT_VERSION; - result.bestnorm = mp_dmax1(fnorm, fnorm1); - result.bestnorm *= result.bestnorm; - result.orignorm = orignorm; - result.status = info; - result.niter = iter; - result.nfev = nfev; - result.npar = npar; - result.nfree = nfree; - result.npegged = npegged; - result.nfunc = m; + result.BestNorm = mp_dmax1(fnorm, fnorm1); + result.BestNorm *= result.BestNorm; + result.OriginalNorm = orignorm; + result.Status = info; + result.Iterations = iter; + result.Evaluations = nfev; + result.ParameterCount = npar; + result.FreeParameterCount = nfree; + result.PeggedParameterCount = npegged; + result.ResidualCount = m; /* Copy residuals if requested */ - if (result.resid != null) + if (result.FinalResiduals != null) { - for (j = 0; j < m; j++) result.resid[j] = fvec[j]; + for (j = 0; j < m; j++) result.FinalResiduals[j] = fvec[j]; } } @@ -2428,15 +2423,6 @@ namespace MathNet.Numerics.Optimization * form the inverse of r in the full upper triangle of r. */ -#if IF0 - for (j=0; jInitial step bound public double stepfactor; - /// Range tolerance for covariance calcu + /// Range tolerance for covariance public double covtol; /// @@ -68,10 +68,10 @@ namespace MathNet.Numerics.Optimization /// errors/covariances are estimated based on input /// parameter values, but no fitting iterations are done. /// - public int maxiter; + public int MaxIterations; /// Maximum number of function evaluations - public int maxfev; + public int MaxEvaluations; /// public int nprint; @@ -81,16 +81,13 @@ namespace MathNet.Numerics.Optimization /// 1 = yes, user scale values in diag; /// 0 = no, variables scaled internally /// - public int douserscale; + public int DoUserScale; /// /// Disable check for infinite quantities from user? /// 0 = do not perform check /// 1 = perform check /// - public int nofinitecheck; - - // /// Placeholder pointer - must set to 0 - //mp_iterproc iterproc; + public int NoFiniteCheck; } } diff --git a/src/Numerics/Optimization/MpResult.cs b/src/Numerics/Optimization/MpResult.cs index 1a9c7fd9..03e51908 100644 --- a/src/Numerics/Optimization/MpResult.cs +++ b/src/Numerics/Optimization/MpResult.cs @@ -48,47 +48,44 @@ namespace MathNet.Numerics.Optimization public class MpResult { /// Final chi^2 - public double bestnorm; + public double BestNorm; /// Starting value of chi^2 - public double orignorm; + public double OriginalNorm; /// Number of iterations - public int niter; + public int Iterations; /// Number of function evaluations - public int nfev; + public int Evaluations; /// Fitting status code - public int status; + public int Status; /// Total number of parameters - public int npar; + public int ParameterCount; /// Number of free parameters - public int nfree; + public int FreeParameterCount; /// Number of pegged parameters - public int npegged; + public int PeggedParameterCount; /// Number of residuals (= num. of data points) - public int nfunc; + public int ResidualCount; /// Final residuals nfunc-vector, or 0 if not desired - public double[] resid; + public double[] FinalResiduals; /// Final parameter uncertainties (1-sigma) npar-vector, or 0 if not desired - public double[] xerror; + public double[] FinalparameterUncertainties; /// Final parameter covariance matrix npar x npar array, or 0 if not desired - public double[] covar; - - /// MPFIT version string - public string version; + public double[] FinalParameterCovarianceMatrix; public MpResult(int numParameters) { - xerror = new double[numParameters]; + FinalparameterUncertainties = new double[numParameters]; } } } diff --git a/src/UnitTests/OptimizationTests/TestMPFit.cs b/src/UnitTests/OptimizationTests/TestMPFit.cs index 7890ad07..d583ed0e 100644 --- a/src/UnitTests/OptimizationTests/TestMPFit.cs +++ b/src/UnitTests/OptimizationTests/TestMPFit.cs @@ -331,27 +331,27 @@ namespace MathNet.Numerics.UnitTests.OptimizationTests if (x == null) return; Console.Write(" CHI-SQUARE = {0} ({1} DOF)\n", - result.bestnorm, result.nfunc - result.nfree); - Console.Write(" NPAR = {0}\n", result.npar); - Console.Write(" NFREE = {0}\n", result.nfree); - Console.Write(" NPEGGED = {0}\n", result.npegged); - Console.Write(" NITER = {0}\n", result.niter); - Console.Write(" NFEV = {0}\n", result.nfev); + result.BestNorm, result.ResidualCount - result.FreeParameterCount); + Console.Write(" NPAR = {0}\n", result.ParameterCount); + Console.Write(" NFREE = {0}\n", result.FreeParameterCount); + Console.Write(" NPEGGED = {0}\n", result.PeggedParameterCount); + Console.Write(" NITER = {0}\n", result.Iterations); + Console.Write(" NFEV = {0}\n", result.Evaluations); Console.Write("\n"); if (xact != null) { - for (i = 0; i < result.npar; i++) + for (i = 0; i < result.ParameterCount; i++) { Console.Write(" P[{0}] = {1} +/- {2} (ACTUAL {3})\n", - i, x[i], result.xerror[i], xact[i]); + i, x[i], result.FinalparameterUncertainties[i], xact[i]); } } else { - for (i = 0; i < result.npar; i++) + for (i = 0; i < result.ParameterCount; i++) { Console.Write(" P[{0}] = {1} +/- {2}\n", - i, x[i], result.xerror[i]); + i, x[i], result.FinalparameterUncertainties[i]); } } }