// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* // ** Copyright UCAR (c) 1992 - 2014 // ** University Corporation for Atmospheric Research(UCAR) // ** National Center for Atmospheric Research(NCAR) // ** Research Applications Laboratory(RAL) // ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA // ** See LICENCE.TXT if applicable for licence details // ** 2014/05/26 09:53:32 // *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* /////////////////////////////////////////////////////////////// // uTest.cc // // TTest object for umath // // Mike Dixon, RAP, NCAR, P.O.Box 3000, Boulder, CO, 80307-3000, USA // // May 2014 // /////////////////////////////////////////////////////////////// #include #include #include #include using namespace std; static void testExponentialFit(double a0, double a1, double a2, double noiseMean, double noiseSdev); static void testPolynomialOrder3(double a0, double a1, double a2, double a3, double noiseMean, double noiseSdev); static void testLinearFit(double mm, double cc, double noiseMean, double noiseSdev); static void testNewtRaph(); static void cubic(double x, double *val, double *deriv); // main int main(int argc, char **argv) { ////////////////////// // Checking uCurveFit testExponentialFit(0.1, 5.0, -5.0, 0.0, 0.1); testExponentialFit(1.0, 0.3, -0.3, 0.0, 0.1); testExponentialFit(0.0, 0.25, -0.25, 0.0, 0.02); testPolynomialOrder3(13.0, 5.0, 0.25, -1.25, 0.0, 0.2); testPolynomialOrder3(10.0, 2.0, -1.5, 1.20, 0.0, 0.2); testLinearFit(0.33333, 150, 0.0, 5.0); testLinearFit(4.25, 9.99, 0.0, 2.0); testLinearFit(-2.50, -8.88, 0.0, 0.1); testNewtRaph(); return 0; } ///////////////////////////////////////////////////////// // test exponential fit static void testExponentialFit(double a0, double a1, double a2, double noiseMean, double noiseSdev) { double aa[3]; aa[0] = a0; aa[1] = a1; aa[2] = a2; double xx[100]; double yy[100]; for (int ii = 0; ii < 100; ii++) { double noise = STATS_normal_gen(noiseMean, noiseSdev); double val = ii; xx[ii] = val; yy[ii] = aa[0] + aa[1] * exp(val * aa[2]) + noise; } cerr << "====================================" << endl; cerr << "Exponential fit" << endl; cerr << "Input aa[0]: " << aa[0] << endl; cerr << "Input aa[1]: " << aa[1] << endl; cerr << "Input aa[2]: " << aa[2] << endl; cerr << "Noise mean: " << noiseMean << endl; cerr << "Noise sdev: " << noiseSdev << endl; double stdErr, rSquared; uExponentialFit(100, xx, yy, aa, &stdErr, &rSquared); cerr << "---------------------------------" << endl; cerr << "Output aa[0]: " << aa[0] << endl; cerr << "Output aa[1]: " << aa[1] << endl; cerr << "Output aa[2]: " << aa[2] << endl; cerr << "====================================" << endl; cerr << endl; } ///////////////////////////////////////////////////////// // test polynomial fit order 3 static void testPolynomialOrder3(double a0, double a1, double a2, double a3, double noiseMean, double noiseSdev) { double aaa[4]; aaa[0] = a0; aaa[1] = a1; aaa[2] = a2; aaa[3] = a3; double xx[100]; double yy[100]; for (int ii = 0; ii < 100; ii++) { double noise = STATS_normal_gen(noiseMean, noiseSdev); double val = ii; xx[ii] = val; yy[ii] = aaa[0] + aaa[1] * val + aaa[2] * val * val + aaa[3] * val * val * val + noise; } cerr << "====================================" << endl; cerr << "Poly fit" << endl; cerr << "Input aaa[0]: " << aaa[0] << endl; cerr << "Input aaa[1]: " << aaa[1] << endl; cerr << "Input aaa[2]: " << aaa[2] << endl; cerr << "Input aaa[3]: " << aaa[3] << endl; cerr << "Noise mean: " << noiseMean << endl; cerr << "Noise sdev: " << noiseSdev << endl; double stdErr, rSquared; uPolyFit(100, xx, yy, aaa, 4, &stdErr, &rSquared); cerr << "---------------------------------" << endl; cerr << "Output aaa[0]: " << aaa[0] << endl; cerr << "Output aaa[1]: " << aaa[1] << endl; cerr << "Output aaa[2]: " << aaa[2] << endl; cerr << "Output aaa[3]: " << aaa[3] << endl; cerr << "====================================" << endl; cerr << endl; } ///////////////////////////////////////////////////////// // test linear fit static void testLinearFit(double mm, double cc, double noiseMean, double noiseSdev) { double xx[100]; double yy[100]; for (int ii = 0; ii < 100; ii++) { double noise = STATS_normal_gen(noiseMean, noiseSdev); double val = ii; xx[ii] = val; yy[ii] = mm * xx[ii] + cc + noise; } cerr << "====================================" << endl; cerr << "Linear fit" << endl; cerr << "Input mm: " << mm << endl; cerr << "Input cc: " << cc << endl; cerr << "Noise mean: " << noiseMean << endl; cerr << "Noise sdev: " << noiseSdev << endl; double aa[2]; double xmean, ymean, xsdev, ysdev; double corr; double stdErr, rSquared; uLinearFit(100, xx, yy, aa, &xmean, &ymean, &xsdev, &ysdev, &corr, &stdErr, &rSquared); cerr << "---------------------------------" << endl; cerr << "Fit mm: " << aa[1] << endl; cerr << "Fit cc: " << aa[0] << endl; cerr << " xmean: " << xmean << endl; cerr << " ymean: " << ymean << endl; cerr << " xsdev: " << xsdev << endl; cerr << " ysdev: " << ysdev << endl; cerr << " corr: " << corr << endl; cerr << " stdErr: " << stdErr << endl; cerr << " rSquared: " << rSquared << endl; uLinearFitPC(100, xx, yy, &mm, &cc); cerr << "---------------------------------" << endl; cerr << "PC Fit mm: " << mm << endl; cerr << "PC Fit cc: " << cc << endl; cerr << "====================================" << endl; cerr << endl; } ///////////////////////////////////////////////////////// // testing uNewtRaph // get roots of cubic static void testNewtRaph() { cerr << "Testing uNewtRaph" << endl; cerr << "Roots should be: -4, -1, 2" << endl; double root; uNewtRaph(cubic, -6, -3, 0.0001, &root); cerr << "========>> root: " << root << endl; uNewtRaph(cubic, -3, 0.0, 0.0001, &root); cerr << "========>> root: " << root << endl; uNewtRaph(cubic, 1.0, 300.0, 0.0001, &root); cerr << "========>> root: " << root << endl; } static void cubic(double x, double *val, double *deriv) { *val = (pow(x, 3.0) + 3 * pow(x, 2.0) - 6.0 * x - 8.0) / 4.0; *deriv = (3.0 * x * x + 6.0 * x - 6.0) / 4.0; }