Index: /trunk/archive/scripts/src/arraytest.c =================================================================== --- /trunk/archive/scripts/src/arraytest.c (revision 3943) +++ /trunk/archive/scripts/src/arraytest.c (revision 3943) @@ -0,0 +1,59 @@ +# include + +// test psMinimizeLMChi2 +int main (int argc, char **argv) { + + int i; + psArray *x; + psVector *y; + static char *data = "test line"; + + // create data array + x = psArrayAlloc (100); + x->n = 0; + for (i = 0; i < 500; i++) { + psArrayAdd (x, 100, psStringCopy (data)); + } + + // create data vector + y = psVectorAlloc (100); + y->n = 0; + for (i = 0; i < y->n; i++) { + n = y->n; + y->data.F32[n] = 2*i; + psVectorExtend (y, 100, 1); + // increments n by 1, extends length if needed by 100 + } + exit (0); +} + +static psHash *timers = NULL; + +bool psTimerStart (char *name) { + + psTime *start; + + if (timers == NULL) timers = psHashAlloc (16); + + start = psTimeGetTime (PS_TIME_UTC); + + psHashAdd (timers, name, start); + + return (TRUE); +} + +psF64 psTimerMark (char *name) { + + psTime *start; + psTime *mark; + psF64 delta; + + if (timers == NULL) return (0); + + start = psHashLookup (timers, name); + if (start == NULL) return (0); + + mark = psTimeGetTime (PS_TIME_UTC); + delta = psTimeDelta (mark, start); + return (delta); +} Index: /trunk/archive/scripts/src/modeltest.c =================================================================== --- /trunk/archive/scripts/src/modeltest.c (revision 3943) +++ /trunk/archive/scripts/src/modeltest.c (revision 3943) @@ -0,0 +1,119 @@ +# include +# include + +/* pars: x, sx, sky, I, */ +psF64 fgaussOD (psVector *dpar, psVector *par, psVector *x) { + + double X, S, Z, R, f; + + X = x->data.F32[0] - par->data.F32[0]; + S = 1.0 / (par->data.F32[1]*par->data.F32[1]); + Z = -0.5*X*X*S; + R = exp (Z); + f = par->data.F32[2]*R + par->data.F32[3]; + + if (dpar != NULL) { + dpar->data.F32[0] = par->data.F32[2]*R*X*S; + dpar->data.F32[1] = dpar->data.F32[0]*X/par->data.F32[1]; + dpar->data.F32[2] = R; + dpar->data.F32[3] = 1; + } + return (f); +} + +/* pars: x, sx, */ +psF64 fgaussODs (psVector *dpar, psVector *par, psVector *x) { + + double X, S, Z, R, f; + + X = x->data.F32[0] - par->data.F32[0]; + S = 1.0 / (par->data.F32[1]*par->data.F32[1]); + Z = -0.5*X*X*S; + f = exp (Z); + + if (dpar != NULL) { + dpar->data.F32[0] = f*X*S; + dpar->data.F32[1] = dpar->data.F32[0]*X/par->data.F32[1]; + } + return (f); +} + +// test psMinimizeLMChi2 +int main_test (int argc, char **argv) { + + int i; + bool rc; + psF32 dx; + psArray *x; + psVector *y; + psVector *v; + psVector *p; + + psTraceSetLevel (".", 3); + + psMinimization *myMin = psMinimizationAlloc (20, 0.001); // (iterations, tolerance) + + // create data vector + x = psArrayAlloc (100); + y = psVectorAlloc (100, PS_TYPE_F32); + for (i = 0; i < x[0].n; i++) { + v = psVectorAlloc (1, PS_TYPE_F32); + v[0].data.F32[0] = i; + x[0].data[i] = (psPtr) v; + dx = i - 30; + y[0].data.F32[i] = exp (-(dx*dx)/(2*3*3)); + // xo = 30.0, sigma = 3.0 + } + p = psVectorAlloc (2, PS_TYPE_F32); + p[0].data.F32[0] = 25.0; + p[0].data.F32[1] = 5.0; + + rc = psMinimizeLMChi2(myMin, NULL, p, NULL, x, y, NULL, fgaussODs); + + for (i = 0; i < p[0].n; i++) { + fprintf (stderr, "par %d: %f\n", i, p[0].data.F32[i]); + } + + exit (0); +} + +// test psMinimizeLMChi2 +int main (int argc, char **argv) { + + int i; + bool rc; + psF32 dx; + psArray *x; + psVector *y; + psVector *v; + psVector *p; + + psTraceSetLevel (".", 3); + + psMinimization *myMin = psMinimizationAlloc (20, 0.001); // (iterations, tolerance) + + // create data vector + x = psArrayAlloc (100); + y = psVectorAlloc (100, PS_TYPE_F32); + for (i = 0; i < x[0].n; i++) { + v = psVectorAlloc (1, PS_TYPE_F32); + v[0].data.F32[0] = i; + x[0].data[i] = (psPtr) v; + dx = i - 30; + y[0].data.F32[i] = 5.0*exp (-(dx*dx)/(2*3*3)) + 4.0; + // sigma = 3.0, xo = 30.0, A = 5.0, B = 4.0 + } + p = psVectorAlloc (4, PS_TYPE_F32); + p[0].data.F32[0] = 25.0; + p[0].data.F32[1] = 5.0; + p[0].data.F32[2] = 8.0; + p[0].data.F32[3] = 2.0; + + rc = psMinimizeLMChi2(myMin, NULL, p, NULL, x, y, NULL, fgaussOD); + + for (i = 0; i < p[0].n; i++) { + fprintf (stderr, "par %d: %f\n", i, p[0].data.F32[i]); + } + + exit (0); +} Index: /trunk/archive/scripts/src/poly2dtest.c =================================================================== --- /trunk/archive/scripts/src/poly2dtest.c (revision 3943) +++ /trunk/archive/scripts/src/poly2dtest.c (revision 3943) @@ -0,0 +1,67 @@ +# include "psphot.h" +# include +double drand48(void); + +int main (int argc, char **argv) { + + int i, j, N; + psPolynomial2D *poly; + + srand48(0); + + // *** test 2D fit without errors + { + psVector *x = psVectorAlloc (1000, PS_TYPE_F64); + psVector *y = psVectorAlloc (1000, PS_TYPE_F64); + psVector *z = psVectorAlloc (1000, PS_TYPE_F64); + psVector *dz = NULL; + + N = 0; + for (i = 0; i < 50; i++) { + for (j = 0; j < 20; j++, N++) { + x->data.F64[N] = i; + y->data.F64[N] = j; + z->data.F64[N] = 3 + 5*i - 2*j; + } + } + + // use the code in psUtils.c + poly = Polynomial2DAlloc (1, 1, PS_POLYNOMIAL_ORD); + poly = VectorFitPolynomial2DOrd_EAM (poly, x, y, z, dz); + + for (i = 0; i < poly->nX + 1; i++) { + for (j = 0; j < poly->nY + 1; j++) { + fprintf (stderr, "x^%d y^%d : %6.3f +/- %6.3f\n", i, j, poly->coeff[i][j], poly->coeffErr[i][j]); + } + } + } + + + // *** test 2D fit with errors + { + psVector *x = psVectorAlloc (1000, PS_TYPE_F64); + psVector *y = psVectorAlloc (1000, PS_TYPE_F64); + psVector *z = psVectorAlloc (1000, PS_TYPE_F64); + psVector *dz = psVectorAlloc (1000, PS_TYPE_F64); + + N = 0; + for (i = 0; i < 50; i++) { + for (j = 0; j < 20; j++, N++) { + x->data.F64[N] = i; + y->data.F64[N] = j; + dz->data.F64[N] = 1*(drand48() - 0.5); + z->data.F64[N] = 3 + 5*i - 2*j + dz->data.F64[N]; + } + } + + // use the code in psUtils.c + poly = Polynomial2DAlloc (2, 2, PS_POLYNOMIAL_ORD); + poly = VectorFitPolynomial2DOrd_EAM (poly, x, y, z, dz); + + for (i = 0; i < poly->nX; i++) { + for (j = 0; j < poly->nY; j++) { + fprintf (stderr, "x^%d y^%d : %f +/- %f\n", i, j, poly->coeff[i][j], poly->coeffErr[i][j]); + } + } + } +} Index: /trunk/archive/scripts/src/polytest.c =================================================================== --- /trunk/archive/scripts/src/polytest.c (revision 3943) +++ /trunk/archive/scripts/src/polytest.c (revision 3943) @@ -0,0 +1,73 @@ +# include "psphot.h" +# include +double drand48(void); + +int main (int argc, char **argv) { + + int i; + psPolynomial1D *poly; + + // *** test 1D fit without errors + { + psVector *x = psVectorAlloc (1000, PS_TYPE_F64); + psVector *y = psVectorAlloc (1000, PS_TYPE_F64); + psVector *dy = NULL; + psVector *dy2 = NULL; + + srand48(0); + + for (i = 0; i < 1000; i++) { + x->data.F64[i] = i; + y->data.F64[i] = 3 + 5*i - 2*i*i +0.4*i*i*i; + } + + // use the code in psUtils.c + poly = psPolynomial1DAlloc (4, PS_POLYNOMIAL_ORD); + poly = VectorFitPolynomial1DOrd_EAM (poly, x, y, dy); + + for (i = 0; i < poly->n; i++) { + fprintf (stderr, "x^%d : %f +/- %f\n", i, poly->coeff[i], poly->coeffErr[i]); + } + + // use the code in psMinimize.c + poly = psPolynomial1DAlloc (4, PS_POLYNOMIAL_ORD); + poly = VectorFitPolynomial1DOrd (poly, x, y, dy2); + + for (i = 0; i < poly->n; i++) { + fprintf (stderr, "x^%d : %f +/- %f\n", i, poly->coeff[i], poly->coeffErr[i]); + } + } + + // *** test 1D fit with errors + { + psVector *x = psVectorAlloc (1000, PS_TYPE_F64); + psVector *y = psVectorAlloc (1000, PS_TYPE_F64); + psVector *dy = psVectorAlloc (1000, PS_TYPE_F64); + psVector *dy2 = psVectorAlloc (1000, PS_TYPE_F64); + + srand48(0); + + for (i = 0; i < 1000; i++) { + x->data.F64[i] = i; + dy->data.F64[i] = 10*(drand48() - 0.5); + dy2->data.F64[i] = PS_SQR(dy->data.F64[i]); + y->data.F64[i] = 3 + 5*i - 2*i*i +0.4*i*i*i + dy->data.F64[i]; + } + + // use the code in psUtils.c + poly = psPolynomial1DAlloc (4, PS_POLYNOMIAL_ORD); + poly = VectorFitPolynomial1DOrd_EAM (poly, x, y, dy); + + for (i = 0; i < poly->n; i++) { + fprintf (stderr, "x^%d : %f +/- %f\n", i, poly->coeff[i], poly->coeffErr[i]); + } + + // use the code in psMinimize.c + poly = psPolynomial1DAlloc (4, PS_POLYNOMIAL_ORD); + poly = VectorFitPolynomial1DOrd (poly, x, y, dy2); + + for (i = 0; i < poly->n; i++) { + fprintf (stderr, "x^%d : %f +/- %f\n", i, poly->coeff[i], poly->coeffErr[i]); + } + } +}