Skip to content
Snippets Groups Projects
Commit 3c6da78e authored by Wuttke, Joachim's avatar Wuttke, Joachim
Browse files

scr+intp basically works

parent dc9972b4
No related branches found
No related tags found
No related merge requests found
......@@ -547,6 +547,13 @@ void globalPrintout( int n_par, const double* par, int m_dat, const void *data,
// iter : outer loop counter
// nfev : number of calls to *evaluate
{
printf( "fitmon %4i", iter );
for( uint ip=0; ip<n_par; ++ip )
printf( " %19.11g", par[ip] );
double msd=0;
for( uint i=0; i<m_dat; ++i )
msd += fvec[i]*fvec[i];
printf( " -> %19.9g\n", msd );
}
......@@ -615,7 +622,7 @@ void NCurveFile::Fit( bool _allow_slow_conv )
POlc fc;
while ( (fc=fiter()) ) {
data.k = fiter.k();
if( !fc->T->has_dummy() )
if( fc->evaMode==COlc::_EXPR && !fc->T->has_dummy() )
throw string( "curve has no formal argument t" );
if( fc->kd>=NOlm::MOM.size() )
throw string( "number of data file out of range" );
......@@ -660,7 +667,7 @@ void NCurveFile::Fit( bool _allow_slow_conv )
data.fc = fc;
data.fd = fd;
uint np = fc->T->npar();
uint np = fc->nPar();
control = lm_control_double;
control.epsilon = NFitTune::Epsilon;
......
......@@ -250,7 +250,6 @@ void COlc::curve_val_vec( vector<double>* ret, const vector<double>& vt,
ret->resize( vt.size() );
CResult val;
if ( evaMode==_EXPR ) {
cout << "DEB expr\n";
CContext ctx( k, j );
if ( !T->has_dummy() ) { // curve does not depend on t
T->tree_val( val, ctx );
......@@ -271,28 +270,26 @@ void COlc::curve_val_vec( vector<double>* ret, const vector<double>& vt,
// string args;
for( uint ip=0; ip<nPar(); ++ip )
fprintf( F, "%g\n", VC(j)->P[ip] );
//for( uint ip=0; ip<nPar(); ++ip )
// printf( "DEB p%i %g\n", ip, VC(j)->P[ip] );
fclose( F );
// args += " " + strg( VC(j)->P[ip] );
//cout << "DEB mkfifo\n";
//mkfifo( "/tmp/f2curve", 0666 );
// string cmd = "'" + expr + args + "' > /ram/f2curve";
string cmd = expr + "< /ram/pars > /ram/f2curve";
cout << "DEB system [" + cmd + "]\n";
system( cmd.c_str() );
cout << "DEB read\n";
//FILE *F;
if( !(F = fopen("/ram/f2curve", "r")) )
throw "BUG: cannot read results of system call";
CSpec S;
double x, y;
cout << "DEB reading ...\n";
while( fscanf( F, "%lg %lg\n", &x, &y )==2 ){
while( fscanf( F, "%lg %lg\n", &x, &y )==2 )
S.push_xy( x, y );
}
fclose( F );
system( "rm /ram/f2curve" );
cout << "DEB no lines: " << S.size() << "\n";
//system( "rm /ram/pars /ram/f2curve" );
// interpolation:
if( S.size()<2 )
throw "script returns only " + strg(S.size()) + " lines";
S.intpol( vt, *ret );
} else
throw "BUG: unexpected evaluation mode";
......
......@@ -137,15 +137,18 @@ uint CSpec::size() const
void CSpec::intpol( const vector<double>& xx, vector<double>& yy ) const
{
uint n = size();
if( n<2 )
throw "BUG: intpol n<2";
uint nn = xx.size();
yy.resize( nn );
if( nn<2 )
throw "BUG: intpol nn<2";
yy.resize( nn, 0. );
uint ii = 0;
while( xx[ii]<x[0] ){
if( ii>=nn )
if( ++ii>=nn )
return;
yy[ii] = 0;
++ii;
// yy[ii] = 0;
}
uint i = 0;
while ( i<n-1 && xx[ii]<x[n-1] ) {
......@@ -158,8 +161,8 @@ void CSpec::intpol( const vector<double>& xx, vector<double>& yy ) const
yy[ii] = y[i] + (xx[ii]-x[i])/(x[i+1]-x[i])*(y[i+1]-y[i]);
++ii;
}
for( ; ii<nn; ++ii )
yy[ii] = 0;
//for( ; ii<nn; ++ii )
// yy[ii] = 0;
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment