1 文本格式
using System;
namespace Legalsoft.Truffer
{
/// <summary>
/// General linear fit
/// </summary>
public class Fitlin
{
private int ndat { get; set; }
private int ma { get; set; }
private double[] x { get; set; }
private double[] y { get; set; }
private double[] sig { get; set; }
private bool[] ia { get; set; }
private double[] a { get; set; }
private double[,] covar { get; set; }
private double chisq { get; set; }
public UniVarRealMultiValueFun funcs;
public Fitlin(double[] xx, double[] yy, double[] ssig, UniVarRealMultiValueFun funks)
{
this.ndat = xx.Length;
this.x = xx;
this.y = yy;
this.sig = ssig;
this.funcs = funks;
ma = funcs.funk(x[0]).Length;
//a.resize(ma);
a = new double[ma];
//covar.resize(ma, ma);
covar = new double[ma, ma];
//ia.resize(ma);
ia = new bool[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 mfit = 0;
for (int j = 0; j < ma; j++)
{
if (ia[j])
{
mfit++;
}
}
if (mfit == 0)
{
throw new Exception("lfit: no parameters to be fitted");
}
double[,] temp = new double[mfit, mfit];
double[,] beta = new double[mfit, 1];
for (int i = 0; i < ndat; i++)
{
double[] afunc = funcs.funk(x[i]);
double ym = y[i];
if (mfit < ma)
{
for (int j = 0; j < ma; j++)
{
if (!ia[j])
{
ym -= a[j] * afunc[j];
}
}
}
double sig2i = 1.0 / Globals.SQR(sig[i]);
for (int j = 0, l = 0; l < ma; l++)
{
if (ia[l])
{
double wt = afunc[l] * sig2i;
for (int k = 0, m = 0; m <= l; m++)
{
if (ia[m])
{
temp[j, k++] += wt * afunc[m];
}
}
beta[j++, 0] += ym * wt;
}
}
}
for (int j = 1; j < mfit; j++)
{
for (int k = 0; k < j; k++)
{
temp[k, j] = temp[j, k];
}
}
GaussJordan.gaussj(temp, beta);
for (int j = 0, l = 0; l < ma; l++)
{
if (ia[l])
{
a[l] = beta[j++, 0];
}
}
chisq = 0.0;
for (int i = 0; i < ndat; i++)
{
double[] afunc = funcs.funk(x[i]);
double sum = 0.0;
for (int j = 0; j < ma; j++)
{
sum += a[j] * afunc[j];
}
chisq += Globals.SQR((y[i] - sum) / sig[i]);
}
for (int j = 0; j < mfit; j++)
{
for (int k = 0; k < mfit; k++)
{
covar[j, k] = temp[j, k];
}
}
for (int i = mfit; i < ma; i++)
{
for (int j = 0; j < i + 1; j++)
{
covar[i, j] = covar[j, i] = 0.0;
}
}
int kk = 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, kk], ref covar[i, j]);
}
for (int i = 0; i < ma; i++)
{
Globals.SWAP(ref covar[kk, i], ref covar[j, i]);
}
kk--;
}
}
}
}
}
2 代码格式
using System;
namespace Legalsoft.Truffer
{
/// <summary>
/// General linear fit
/// </summary>
public class Fitlin
{
private int ndat { get; set; }
private int ma { get; set; }
private double[] x { get; set; }
private double[] y { get; set; }
private double[] sig { get; set; }
private bool[] ia { get; set; }
private double[] a { get; set; }
private double[,] covar { get; set; }
private double chisq { get; set; }
public UniVarRealMultiValueFun funcs;
public Fitlin(double[] xx, double[] yy, double[] ssig, UniVarRealMultiValueFun funks)
{
this.ndat = xx.Length;
this.x = xx;
this.y = yy;
this.sig = ssig;
this.funcs = funks;
ma = funcs.funk(x[0]).Length;
//a.resize(ma);
a = new double[ma];
//covar.resize(ma, ma);
covar = new double[ma, ma];
//ia.resize(ma);
ia = new bool[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 mfit = 0;
for (int j = 0; j < ma; j++)
{
if (ia[j])
{
mfit++;
}
}
if (mfit == 0)
{
throw new Exception("lfit: no parameters to be fitted");
}
double[,] temp = new double[mfit, mfit];
double[,] beta = new double[mfit, 1];
for (int i = 0; i < ndat; i++)
{
double[] afunc = funcs.funk(x[i]);
double ym = y[i];
if (mfit < ma)
{
for (int j = 0; j < ma; j++)
{
if (!ia[j])
{
ym -= a[j] * afunc[j];
}
}
}
double sig2i = 1.0 / Globals.SQR(sig[i]);
for (int j = 0, l = 0; l < ma; l++)
{
if (ia[l])
{
double wt = afunc[l] * sig2i;
for (int k = 0, m = 0; m <= l; m++)
{
if (ia[m])
{
temp[j, k++] += wt * afunc[m];
}
}
beta[j++, 0] += ym * wt;
}
}
}
for (int j = 1; j < mfit; j++)
{
for (int k = 0; k < j; k++)
{
temp[k, j] = temp[j, k];
}
}
GaussJordan.gaussj(temp, beta);
for (int j = 0, l = 0; l < ma; l++)
{
if (ia[l])
{
a[l] = beta[j++, 0];
}
}
chisq = 0.0;
for (int i = 0; i < ndat; i++)
{
double[] afunc = funcs.funk(x[i]);
double sum = 0.0;
for (int j = 0; j < ma; j++)
{
sum += a[j] * afunc[j];
}
chisq += Globals.SQR((y[i] - sum) / sig[i]);
}
for (int j = 0; j < mfit; j++)
{
for (int k = 0; k < mfit; k++)
{
covar[j, k] = temp[j, k];
}
}
for (int i = mfit; i < ma; i++)
{
for (int j = 0; j < i + 1; j++)
{
covar[i, j] = covar[j, i] = 0.0;
}
}
int kk = 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, kk], ref covar[i, j]);
}
for (int i = 0; i < ma; i++)
{
Globals.SWAP(ref covar[kk, i], ref covar[j, i]);
}
kk--;
}
}
}
}
}
