1 文本格式

using System;

namespace Legalsoft.Truffer
{
    /// <summary>
    /// Levenberg-Marquardt nonlinear fitting
    /// </summary>
    public class Fitmrq
    {
        private const int NDONE = 4;
        private const int ITMAX = 1000;
        private int ndat { get; set; }
        private int ma { get; set; }
        private int mfit { get; set; }
        private double[] x { get; set; }
        private double[] y { get; set; }
        private double[] sig { get; set; }
        private double tol { get; set; }
        private bool[] ia { get; set; }
        private double[] a { get; set; }
        private double[,] covar;
        private double[,] alpha;
        private double chisq { get; set; }
        public MultiFuncd funcs;

        public Fitmrq(double[] xx, double[] yy, double[] ssig, double[] aa, MultiFuncd funks, double TOL = 1.0e-3)
        {
            this.ndat = xx.Length;
            this.ma = aa.Length;
            this.x = xx;
            this.y = yy;
            this.sig = ssig;
            this.tol = TOL;
            this.funcs = funks;
            this.ia = new bool[ma];
            this.alpha = new double[ma, ma];
            this.a = aa;
            this.covar = new double[ma, ma];
            for (int i = 0; i < ma; i++)
            {
                ia[i] = true;
            }
        }

        public void hold(int i, double val)
        {
            ia[i] = false;
            a[i] = val;
        }

        public void free(int i)
        {
            ia[i] = true;
        }

        public void fit()
        {
            int done = 0;
            double alamda = .001;
            double[] atry = new double[ma];
            double[] beta = new double[ma];
            double[] da = new double[ma];
            mfit = 0;
            for (int j = 0; j < ma; j++)
            {
                if (ia[j])
                {
                    mfit++;
                }
            }
            double[,] oneda = new double[mfit, 1];
            double[,] temp = new double[mfit, mfit];
            mrqcof(a, alpha, beta);
            for (int j = 0; j < ma; j++)
            {
                atry[j] = a[j];
            }
            double ochisq = chisq;
            for (int iter = 0; iter < ITMAX; iter++)
            {
                if (done == NDONE)
                {
                    alamda = 0.0;
                }
                for (int j = 0; j < mfit; j++)
                {
                    for (int k = 0; k < mfit; k++)
                    {
                        covar[j, k] = alpha[j, k];
                    }
                    covar[j, j] = alpha[j, j] * (1.0 + alamda);
                    for (int k = 0; k < mfit; k++)
                    {
                        temp[j, k] = covar[j, k];
                    }
                    oneda[j, 0] = beta[j];
                }
                GaussJordan.gaussj( temp,  oneda);
                for (int j = 0; j < mfit; j++)
                {
                    for (int k = 0; k < mfit; k++)
                    {
                        covar[j, k] = temp[j, k];
                    }
                    da[j] = oneda[j, 0];
                }
                if (done == NDONE)
                {
                    covsrt( covar);
                    covsrt( alpha);
                    return;
                }
                for (int j = 0, l = 0; l < ma; l++)
                {
                    if (ia[l])
                    {
                        atry[l] = a[l] + da[j++];
                    }
                }
                mrqcof(atry, covar, da);
                if (Math.Abs(chisq - ochisq) < Math.Max(tol, tol * chisq))
                {
                    done++;
                }
                if (chisq < ochisq)
                {
                    alamda *= 0.1;
                    ochisq = chisq;
                    for (int j = 0; j < mfit; j++)
                    {
                        for (int k = 0; k < mfit; k++)
                        {
                            alpha[j, k] = covar[j, k];
                        }
                        beta[j] = da[j];
                    }
                    for (int l = 0; l < ma; l++)
                    {
                        a[l] = atry[l];
                    }
                }
                else
                {
                    alamda *= 10.0;
                    chisq = ochisq;
                }
            }
            throw new Exception("Fitmrq too many iterations");
        }

        public void mrqcof(double[] a, double[,] alpha, double[] beta)
        {
            double ymod = 0.0;
            double[] dyda = new double[ma];
            for (int j = 0; j < mfit; j++)
            {
                for (int k = 0; k <= j; k++)
                {
                    alpha[j, k] = 0.0;
                }
                beta[j] = 0.0;
            }
            chisq = 0.0;
            for (int i = 0; i < ndat; i++)
            {
                funcs.funk(x[i], a, ref ymod,  dyda);
                double sig2i = 1.0 / (sig[i] * sig[i]);
                double dy = y[i] - ymod;
                for (int j = 0, l = 0; l < ma; l++)
                {
                    if (ia[l])
                    {
                        double wt = dyda[l] * sig2i;
                        for (int k = 0, m = 0; m < l + 1; m++)
                        {
                            if (ia[m])
                            {
                                alpha[j, k++] += wt * dyda[m];
                            }
                        }
                        beta[j++] += dy * wt;
                    }
                }
                chisq += dy * dy * sig2i;
            }
            for (int j = 1; j < mfit; j++)
            {
                for (int k = 0; k < j; k++)
                {
                    alpha[k, j] = alpha[j, k];
                }
            }
        }

        public void covsrt(double[,] covar)
        {
            for (int i = mfit; i < ma; i++)
            {
                for (int j = 0; j < i + 1; j++)
                {
                    covar[i, j] = covar[j, i] = 0.0;
                }
            }
            int k = mfit - 1;
            for (int j = ma - 1; j >= 0; j--)
            {
                if (ia[j])
                {
                    for (int i = 0; i < ma; i++)
                    {
                        Globals.SWAP(ref covar[i, k], ref covar[i, j]);
                    }
                    for (int i = 0; i < ma; i++)
                    {
                        Globals.SWAP(ref covar[k, i], ref covar[j, i]);
                    }
                    k--;
                }
            }
        }
    }
}
 

2 代码格式

using System;

namespace Legalsoft.Truffer
{
    /// <summary>
    /// Levenberg-Marquardt nonlinear fitting
    /// </summary>
    public class Fitmrq
    {
        private const int NDONE = 4;
        private const int ITMAX = 1000;
        private int ndat { get; set; }
        private int ma { get; set; }
        private int mfit { get; set; }
        private double[] x { get; set; }
        private double[] y { get; set; }
        private double[] sig { get; set; }
        private double tol { get; set; }
        private bool[] ia { get; set; }
        private double[] a { get; set; }
        private double[,] covar;
        private double[,] alpha;
        private double chisq { get; set; }
        public MultiFuncd funcs;

        public Fitmrq(double[] xx, double[] yy, double[] ssig, double[] aa, MultiFuncd funks, double TOL = 1.0e-3)
        {
            this.ndat = xx.Length;
            this.ma = aa.Length;
            this.x = xx;
            this.y = yy;
            this.sig = ssig;
            this.tol = TOL;
            this.funcs = funks;
            this.ia = new bool[ma];
            this.alpha = new double[ma, ma];
            this.a = aa;
            this.covar = new double[ma, ma];
            for (int i = 0; i < ma; i++)
            {
                ia[i] = true;
            }
        }

        public void hold(int i, double val)
        {
            ia[i] = false;
            a[i] = val;
        }

        public void free(int i)
        {
            ia[i] = true;
        }

        public void fit()
        {
            int done = 0;
            double alamda = .001;
            double[] atry = new double[ma];
            double[] beta = new double[ma];
            double[] da = new double[ma];
            mfit = 0;
            for (int j = 0; j < ma; j++)
            {
                if (ia[j])
                {
                    mfit++;
                }
            }
            double[,] oneda = new double[mfit, 1];
            double[,] temp = new double[mfit, mfit];
            mrqcof(a, alpha, beta);
            for (int j = 0; j < ma; j++)
            {
                atry[j] = a[j];
            }
            double ochisq = chisq;
            for (int iter = 0; iter < ITMAX; iter++)
            {
                if (done == NDONE)
                {
                    alamda = 0.0;
                }
                for (int j = 0; j < mfit; j++)
                {
                    for (int k = 0; k < mfit; k++)
                    {
                        covar[j, k] = alpha[j, k];
                    }
                    covar[j, j] = alpha[j, j] * (1.0 + alamda);
                    for (int k = 0; k < mfit; k++)
                    {
                        temp[j, k] = covar[j, k];
                    }
                    oneda[j, 0] = beta[j];
                }
                GaussJordan.gaussj( temp,  oneda);
                for (int j = 0; j < mfit; j++)
                {
                    for (int k = 0; k < mfit; k++)
                    {
                        covar[j, k] = temp[j, k];
                    }
                    da[j] = oneda[j, 0];
                }
                if (done == NDONE)
                {
                    covsrt( covar);
                    covsrt( alpha);
                    return;
                }
                for (int j = 0, l = 0; l < ma; l++)
                {
                    if (ia[l])
                    {
                        atry[l] = a[l] + da[j++];
                    }
                }
                mrqcof(atry, covar, da);
                if (Math.Abs(chisq - ochisq) < Math.Max(tol, tol * chisq))
                {
                    done++;
                }
                if (chisq < ochisq)
                {
                    alamda *= 0.1;
                    ochisq = chisq;
                    for (int j = 0; j < mfit; j++)
                    {
                        for (int k = 0; k < mfit; k++)
                        {
                            alpha[j, k] = covar[j, k];
                        }
                        beta[j] = da[j];
                    }
                    for (int l = 0; l < ma; l++)
                    {
                        a[l] = atry[l];
                    }
                }
                else
                {
                    alamda *= 10.0;
                    chisq = ochisq;
                }
            }
            throw new Exception("Fitmrq too many iterations");
        }

        public void mrqcof(double[] a, double[,] alpha, double[] beta)
        {
            double ymod = 0.0;
            double[] dyda = new double[ma];
            for (int j = 0; j < mfit; j++)
            {
                for (int k = 0; k <= j; k++)
                {
                    alpha[j, k] = 0.0;
                }
                beta[j] = 0.0;
            }
            chisq = 0.0;
            for (int i = 0; i < ndat; i++)
            {
                funcs.funk(x[i], a, ref ymod,  dyda);
                double sig2i = 1.0 / (sig[i] * sig[i]);
                double dy = y[i] - ymod;
                for (int j = 0, l = 0; l < ma; l++)
                {
                    if (ia[l])
                    {
                        double wt = dyda[l] * sig2i;
                        for (int k = 0, m = 0; m < l + 1; m++)
                        {
                            if (ia[m])
                            {
                                alpha[j, k++] += wt * dyda[m];
                            }
                        }
                        beta[j++] += dy * wt;
                    }
                }
                chisq += dy * dy * sig2i;
            }
            for (int j = 1; j < mfit; j++)
            {
                for (int k = 0; k < j; k++)
                {
                    alpha[k, j] = alpha[j, k];
                }
            }
        }

        public void covsrt(double[,] covar)
        {
            for (int i = mfit; i < ma; i++)
            {
                for (int j = 0; j < i + 1; j++)
                {
                    covar[i, j] = covar[j, i] = 0.0;
                }
            }
            int k = mfit - 1;
            for (int j = ma - 1; j >= 0; j--)
            {
                if (ia[j])
                {
                    for (int i = 0; i < ma; i++)
                    {
                        Globals.SWAP(ref covar[i, k], ref covar[i, j]);
                    }
                    for (int i = 0; i < ma; i++)
                    {
                        Globals.SWAP(ref covar[k, i], ref covar[j, i]);
                    }
                    k--;
                }
            }
        }
    }
}