From 18e449fe2b8d5b3a2507664917675e07394b73e8 Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (o)" <j.wuttke@fz-juelich.de>
Date: Tue, 16 Dec 2014 10:42:02 +0100
Subject: [PATCH] Replace fit "Quality" by explicit outcome, chi2, R2.
---
TODO | 2 ++
pub/lib/fit.cpp | 26 +++++++-------------------
pub/lib/integrate.cpp | 14 +++++++-------
pub/lib/node.cpp | 26 +++++++++++++++-----------
pub/lib/var.cpp | 15 ++++++++++-----
pub/lib/var.hpp | 5 +++--
pub/lib/xax_lex.lpp | 2 +-
pub/lib/zentry.cpp | 12 +++++++-----
pub/lib/zentry.hpp | 5 +++--
9 files changed, 55 insertions(+), 52 deletions(-)
diff --git a/TODO b/TODO
index 9acf54b3..7dfa3ad1 100644
--- a/TODO
+++ b/TODO
@@ -4,6 +4,8 @@ cq0/q0 etc -> more explicit names throughout (cp,oi,..);
get rid of ci (geni 'integrate')
automatic mfj after oi
+update infoLine header
+
shibang mechanism to invoke Frida scripts
cp should show quality of last fit
latest action on pX should be part of history
diff --git a/pub/lib/fit.cpp b/pub/lib/fit.cpp
index fdffb062..f1d4e498 100644
--- a/pub/lib/fit.cpp
+++ b/pub/lib/fit.cpp
@@ -32,7 +32,6 @@ namespace NCurveFit {
double Factor=100.;
int nCall=200;
int verbosity=0;
- int mFitMetric=2;
double maxtime = 3;
}
@@ -42,8 +41,8 @@ void NCurveFit::setFitTuningPars( string which )
{
if ( which=="?" )
printf( "Fit setup:\n nCall=%d\n Verbosity=%d\n Factor=%g\n"
- " Eps=%g\n Metric=%i\n",
- nCall, verbosity, Factor, Epsilon, mFitMetric );
+ " Eps=%g\n",
+ nCall, verbosity, Factor, Epsilon );
else if ( which=="v" )
verbosity = iask( "Verbosity", verbosity );
else if ( which=="c" )
@@ -52,17 +51,7 @@ void NCurveFit::setFitTuningPars( string which )
Factor = dask( "NTuneFit::Factor", Factor );
else if ( which=="e" )
Epsilon = dask( "NTuneFit::Epsilon", Epsilon );
- else if ( which=="m" ){
- // metric (0) is the integer indicator returned by lmmin;
- // it cannot be chosen here
- cout << "Fit metric:\n"
- "(1) R^2 (coefficient of determination, for linear regression)\n"
- "(2) chi^2 (mean deviation, with weights as in fit)\n";
- int ret = iask( "Option", mFitMetric );
- if( ret<1 || ret>CCurve::mQuality )
- throw "invalid option";
- mFitMetric = ret;
- } else if ( which=="t" )
+ else if ( which=="t" )
maxtime = dask( "Max duration of fits (in s)", maxtime );
}
@@ -267,7 +256,7 @@ void NCurveFit::fit( bool _allow_slow_conv )
// Fit quality metrics Q:
// Q: status code
- C->Quality[0] = (double) status.outcome;
+ C->fitOutcome = (double) status.outcome;
// Q: coefficient of determination R^2 (for linear regression)
double ysum1 = 0, ysum2 = 0;
@@ -276,12 +265,12 @@ void NCurveFit::fit( bool _allow_slow_conv )
ysum2 += SQR( D->y[i] );
}
double ydev2 = ysum2 - SQR( ysum1 )/nd; // = nd * <(y-<y>)^2>
- C->Quality[1] = ydev2!=0 ? 1 - SQR(status.fnorm) / ydev2 : -1;
+ C->fitR2 = ydev2!=0 ? 1 - SQR(status.fnorm) / ydev2 : -1;
// Q: mean deviation (chi^2 if fvec is weighed with 1/dy)
int nfreedom = nd - npfree - 1;
if( nfreedom<1 ) nfreedom = 1;
- C->Quality[2] = SQR(status.fnorm) / nfreedom;
+ C->fitChi2 = SQR(status.fnorm) / nfreedom;
{
// print results:
@@ -291,7 +280,7 @@ void NCurveFit::fit( bool _allow_slow_conv )
for ( int i=0; i<npfree; i++ )
printf(" %12g", Par[i]);
printf( " (%3d->%9.7g)",
- status.nfev, C->Quality[mFitMetric] );
+ status.nfev, C->fitChi2 );
printf( "%s\n", lm_shortmsg[ status.outcome ] );
}
} catch (...) {
@@ -301,4 +290,3 @@ void NCurveFit::fit( bool _allow_slow_conv )
} // j
} // k
}
-
diff --git a/pub/lib/integrate.cpp b/pub/lib/integrate.cpp
index 2c874985..74945c86 100644
--- a/pub/lib/integrate.cpp
+++ b/pub/lib/integrate.cpp
@@ -92,7 +92,7 @@ void NCurveIntegral::integrate()
if( NRead::stack_empty() )
cout << "Generate integral file:\n"
" 0 | 1 | ... extract curve parameters\n"
- " q0 | q1 | q2 extract fit metrics\n"
+ " res | chi2 | R2 extract fit metrics\n"
" i numeric integral\n"
" s integral as sum over data file grid\n"
" <list> comma-separated list thereof\n"
@@ -191,12 +191,12 @@ void NCurveIntegral::integrate()
double val;
if ( mode=="p" )
val = fin->VC(j)->P[ip];
- else if ( mode=="q0" )
- val = fin->VC(j)->Quality[0];
- else if ( mode=="q1" )
- val = fin->VC(j)->Quality[1];
- else if ( mode=="q2" )
- val = fin->VC(j)->Quality[2];
+ else if ( mode=="res" )
+ val = fin->VC(j)->fitOutcome;
+ else if ( mode=="chi2" )
+ val = fin->VC(j)->fitChi2;
+ else if ( mode=="R2" )
+ val = fin->VC(j)->fitR2;
else if ( mode=="i" )
val = NumericIntegral( fin, fiter.k(), j,
numint_bound_low, numint_bound_hig );
diff --git a/pub/lib/node.cpp b/pub/lib/node.cpp
index 57315f89..fb461ec7 100644
--- a/pub/lib/node.cpp
+++ b/pub/lib/node.cpp
@@ -520,23 +520,22 @@ void CNodeIva::tree_val( CResult& ret, const CContext& ctx ) const
// References for curve file:
POlc fc = P2C( f );
if ( fc ) {
+ if ( ref->ti )
+ throw string("unexpected i dependence of curve property");
int j = ref->get_j( ctx, f->nJ() );
PCurve cj = fc->VC(j);
if ( var->typ == CVar::_CP ) {
- if ( ref->ti )
- throw string("P does not require 3 indices");
if ( var->num>= fc->nP )
throw string("invalid p ref(") + var->var_info() +
") in curve file";
ret.set_r( cj->P[var->num] );
- } else if ( var->typ == CVar::_CQ ) {
- if ( ref->ti )
- throw string("Q does not require 3 indices");
- if ( var->num >= CCurve::mQuality )
- throw string("invalid fm ref(") + var->var_info() +
- ") in curve file";
- ret.set_r( cj->Quality[var->num] );
+ } else if ( var->typ == CVar::_COUTCOME ) {
+ ret.set_r( cj->fitOutcome );
+ } else if ( var->typ == CVar::_CCHI2 ) {
+ ret.set_r( cj->fitChi2 );
+ } else if ( var->typ == CVar::_CR2 ) {
+ ret.set_r( cj->fitR2 );
} else
throw string("invalid ref(") + var->var_info() +
@@ -604,8 +603,13 @@ void CNodeIva::set_coord( CCoord& ret, int k_in ) const
if( var->num>=fc->nP )
throw "invalid reference p" + S(var->num);
ret = fc->PCo[var->num];
- } else if ( var->typ == CVar::_CQ )
- ret = CCoord("cq"+S(var->num), ""); // fit quality indicator
+ }
+ else if ( var->typ == CVar::_COUTCOME )
+ ret = CCoord("fit_outcome", "");
+ else if ( var->typ == CVar::_CCHI2 )
+ ret = CCoord("chi2", "");
+ else if ( var->typ == CVar::_CR2 )
+ ret = CCoord("R2", "");
else if ( var->typ == CVar::_C )
ret = CCoord(fc->expr,"");
else
diff --git a/pub/lib/var.cpp b/pub/lib/var.cpp
index e538f2cf..8e317d1c 100644
--- a/pub/lib/var.cpp
+++ b/pub/lib/var.cpp
@@ -39,10 +39,7 @@ CVar::CVar( string s )
if( s.size()<1 )
throw "Variable name missing";
- if ( s.substr(0,2)=="cq" ) { // curve quality
- num = requestNum( s, 2 );
- typ = _CQ;
- } else if ( s[0]=='p' ) { // curve parameter
+ if ( s[0]=='p' ) { // curve parameter
num = requestNum( s, 1 );
typ = _CP;
} else if ( s=="z+" ) { // new z value
@@ -67,6 +64,12 @@ CVar::CVar( string s )
typ = _DY;
} else if ( s=="c" ) {
typ = _C;
+ } else if ( s=="outcome" ) {
+ typ = _COUTCOME;
+ } else if ( s=="chi2" ) {
+ typ = _CCHI2;
+ } else if ( s=="R2" ) {
+ typ = _CR2;
} else if ( s=="k" ) {
typ = _K;
} else if ( s=="j" ) {
@@ -106,7 +109,9 @@ string CVar::var_info() const
case _Z: ret = "z"; break;
case _R: ret = "r"; break;
case _CP: ret = "p"; break;
- case _CQ: ret = "cq"; break;
+ case _COUTCOME: ret = "fit_outcome"; break;
+ case _CCHI2: ret = "chi2"; break;
+ case _CR2: ret = "R2"; break;
case _K: ret = "k"; break;
case _J: ret = "j"; break;
case _I: ret = "i"; break;
diff --git a/pub/lib/var.hpp b/pub/lib/var.hpp
index 44a54df1..b61e7d8b 100644
--- a/pub/lib/var.hpp
+++ b/pub/lib/var.hpp
@@ -16,7 +16,8 @@ class CVar {
_NOREF=0,
_X, _Y, _Z, _R, // data
_DY, // associated error
- _C, _CP, _CQ, // curve, curve parameters, curve quality
+ _C, _CP, // curve, curve parameters
+ _COUTCOME, _CCHI2, _CR2, // fit outcome
_K, _J, _I, // indices for file, spectrum, point
_NJ, _NI // number of spectra in file, points in spectrum
};
@@ -31,7 +32,7 @@ class CVar {
bool pointwise() const {
return typ==_X || typ==_Y || typ==_DY || typ==_I; };
bool numbered() const {
- return typ==_CP || typ==_CQ || typ==_Z; };
+ return typ==_CP || typ==_Z; };
CVar errvar() const;
diff --git a/pub/lib/xax_lex.lpp b/pub/lib/xax_lex.lpp
index dbdd81b1..2085485e 100644
--- a/pub/lib/xax_lex.lpp
+++ b/pub/lib/xax_lex.lpp
@@ -58,7 +58,7 @@ EXP [eE]"-"?[0-9]+
xaxlval->v = PVar( new CVar(xaxtext) );
return IDX; }
-[xy]|dy|n[ij]|([zpr]|cq){DIG}+ {
+[xy]|dy|n[ij]|([zpr]{DIG}+)|outcome|chi2|R2 {
xaxlval->v = PVar( new CVar(xaxtext) );
return REF; }
diff --git a/pub/lib/zentry.cpp b/pub/lib/zentry.cpp
index b3b437da..04c05885 100644
--- a/pub/lib/zentry.cpp
+++ b/pub/lib/zentry.cpp
@@ -275,9 +275,9 @@ void CCurve::clear(void)
{
P.clear();
z.clear();
- Quality[0] = -1;
- Quality[1] = -1;
- Quality[2] = -1;
+ fitOutcome = -1;
+ fitChi2 = 0;
+ fitR2 = 0;
}
//! Return line containing z and p values.
@@ -299,9 +299,11 @@ string CCurve::infoLine() const
if( ip<P.size()-1 )
out += " ";
}
- snprintf( wrd, LEN, "%-2li", lround(Quality[0]) );
+ snprintf( wrd, LEN, "%-2i", fitOutcome );
out += wrd;
- snprintf( wrd, LEN, "%-12.7g", Quality[2] );
+ snprintf( wrd, LEN, "%-12.7g", fitChi2 );
+ out += wrd;
+ snprintf( wrd, LEN, "%-12.7g", fitR2 );
out += wrd;
return out;
}
diff --git a/pub/lib/zentry.hpp b/pub/lib/zentry.hpp
index 5a343aaf..59260d53 100644
--- a/pub/lib/zentry.hpp
+++ b/pub/lib/zentry.hpp
@@ -26,10 +26,11 @@ class CZentry {
class CCurve: public CZentry {
public:
- static const int mQuality = 3;
vector<double> P;
vector<bool> fixed;
- double Quality[mQuality];
+ int fitOutcome;
+ double fitChi2;
+ double fitR2;
void clear( void );
string infoLine() const;
--
GitLab