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

oio seems to work 

18 TODO's left
parent 1839717f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
pub/src/STYLE 0 → 100644
+ 21
− 0
View file @ 169ba463
== function arguments and "const" qualifier ==
- never use "const" with simple data types
double f( double x ); // and not ( const double x )
- prefer passing complex data types (including string) by reference:
double sum( const vector<int>& row );
- generally, when passing arguments by & reference: not putting "const"
means that mutation of the argument is intended
- omitting "const" has the side effect, that the function cannot be
called with a constructor in place of the argument:
out = sum( vector<int>(7) ); // only allowed when arg is declared const
- "const" with pointer arguments is poorly understood (even less so with
shared_ptr)
pub/src/opr.cpp
+ 55
− 62
View file @ 169ba463
...@@ -820,7 +820,7 @@ void NOperate::FunctionalAutocorrelate() ...@@ -820,7 +820,7 @@ void NOperate::FunctionalAutocorrelate()
namespace NIntOld { namespace NIntOld {
static int mod; static int mod;
double Eval(PSpec& S); double evaluate( const PSpec& S );
static double arg2; static double arg2;
static int iarg2; static int iarg2;
}; };
...@@ -870,17 +870,12 @@ void NIntOld::Opr() ...@@ -870,17 +870,12 @@ void NIntOld::Opr()
int nz = fin->ZCo.size(); int nz = fin->ZCo.size();
bool savable = nz > 0; bool savable = nz > 0;
POld fout; // TODO fout = *fin; POld fout( new COld( *fin ) );
cerr << "TODO probably broken after V rewrite\n";
fout->V.clear(); fout->V.clear();
CCoord zco; fout->ZCo.pop_back();
double xval, zval, result; fout->xco = fin->ZCo.back();
if (savable) {
zco = fin->ZCo.back();
fout->ZCo.pop_back();
fout->xco = zco;
}
if (mod== 2) { if (mod== 2) {
fout->yco.name = fin->yco.name + " (#" fout->yco.name = fin->yco.name + " (#"
+ strg(iarg2) + ")"; + strg(iarg2) + ")";
...@@ -931,164 +926,162 @@ void NIntOld::Opr() ...@@ -931,164 +926,162 @@ void NIntOld::Opr()
"|" + fin->yco.name + ")"; "|" + fin->yco.name + ")";
fout->yco.unit = "" ; fout->yco.unit = "" ;
} else { } else {
cout << "PROGRAM ERROR not implemented\n"; throw string( "invalid mode" );
return;
} }
fout->lDoc.push_back( "oio "+strg(mod)+" # y = "+fout->yco.name ); fout->lDoc.push_back( "oio "+strg(mod)+" # y = "+fout->yco.name );
PSpec S( new CSpec ); PSpec eout( new CSpec );
S->z = fin->V[0]->z; eout->z = fin->V[0]->z;
S->z.pop_back(); eout->z.pop_back();
for (uint j=0; j<fin->nJ(); j++) { for (uint j=0; j<fin->nJ(); j++) {
throw string( "TODO broken after V rewrite" ); double xval, zval, result;
// ?????? result = Eval( fin->V[j] ); result = NIntOld::evaluate( fin->VS(j) );
if (!savable) { if (!savable) {
cout << "results: " << result << "\n"; cout << "results: " << result << "\n";
break; continue;
} }
xval = fin->V[j]->z[nz-1]; xval = fin->V[j]->z[nz-1];
S->push_xy(xval, result); eout->push_xy(xval, result);
if (nz>=2) { // new spectrum if jump in other z values if (nz>=2) { // new spectrum if jump in other z values
zval = fin->V[j]->z[nz-2]; zval = fin->V[j]->z[nz-2];
if ((j+1)<fin->nJ() && if ((j+1)<fin->nJ() &&
fin->V[j+1]->z[nz-2]!=zval) { fin->V[j+1]->z[nz-2]!=zval) {
fout->V.push_back(S); fout->V.push_back(eout);
S = PSpec( new CSpec ); eout = PSpec( new CSpec );
S->z = fin->V[j+1]->z; eout->z = fin->V[j+1]->z;
S->z.pop_back(); eout->z.pop_back();
} }
} }
} }
if (savable) { if (savable) {
fout->V.push_back(S); fout->V.push_back(eout);
NOlm::OloAdd( fout ); NOlm::OloAdd( fout );
} }
} }
}; };
double NIntOld::Eval( PSpec& S ) double NIntOld::evaluate( const PSpec& ein )
{ {
int i, i0; int i, i0;
double r0, r1, r2, s0, s1, v, dr, area, triangle, p, center; double r0, r1, r2, s0, s1, v, dr, area, triangle, p, center;
int n = S->size(); int n = ein->size();
switch(mod) { switch(mod) {
case 2: case 2:
if (iarg2<0 || iarg2>=n) if (iarg2<0 || iarg2>=n)
throw string( "i out of range" ); throw string( "i out of range" );
return S->y[iarg2]; return ein->y[iarg2];
case 3: case 3:
r0 = 0; r0 = 0;
for (i=0; i<n; i++) { for (i=0; i<n; i++) {
r0 += S->y[i]; r0 += ein->y[i];
} }
return r0/n; return r0/n;
case 4: case 4:
r0 = 0; s0 = 0; r0 = 0; s0 = 0;
for (i=0; i<n; i++) { for (i=0; i<n; i++) {
r0 += S->y[i] / n; r0 += ein->y[i] / n;
} }
for (i=0; i<n; i++) { for (i=0; i<n; i++) {
s0 += SQR(S->y[i]-r0) / n; s0 += SQR(ein->y[i]-r0) / n;
} }
return sqrt(s0); return sqrt(s0);
case 5: case 5:
i0 = 0; i0 = 0;
for (i=1; i<n; i++) { for (i=1; i<n; i++) {
if( S->y[i]>S->y[i0] ) i0 = i; if( ein->y[i]>ein->y[i0] ) i0 = i;
} }
return S->x[i0]; return ein->x[i0];
case 6: case 6:
i0 = 0; i0 = 0;
for (i=1; i<n; i++) { for (i=1; i<n; i++) {
if( S->y[i]>S->y[i0] ) i0 = i; if( ein->y[i]>ein->y[i0] ) i0 = i;
} }
return S->y[i0]; return ein->y[i0];
case 7: case 7:
i0 = 0; i0 = 0;
for (i=1; i<n; i++) { for (i=1; i<n; i++) {
if( S->y[i]<S->y[i0] ) i0 = i; if( ein->y[i]<ein->y[i0] ) i0 = i;
} }
return S->x[i0]; return ein->x[i0];
case 8: case 8:
i0 = 0; i0 = 0;
for (i=1; i<n; i++) { for (i=1; i<n; i++) {
if( S->y[i]<S->y[i0] ) i0 = i; if( ein->y[i]<ein->y[i0] ) i0 = i;
} }
return S->y[i0]; return ein->y[i0];
case 9: // sum_i y case 9: // sum_i y
r0 = 0; r0 = 0;
for (i=0; i<n; i++) for (i=0; i<n; i++)
r0 += S->y[i]; r0 += ein->y[i];
return r0; return r0;
case 10: // integral dx y case 10: // integral dx y
r0 = 0; r0 = 0;
for (i=0; i<n-1; i++) { for (i=0; i<n-1; i++) {
if (S->x[i+1]<=S->x[i]) if (ein->x[i+1]<=ein->x[i])
throw string( "x not sorted" ); throw string( "x not sorted" );
r0 += (S->x[i+1]-S->x[i])*(S->y[i+1]+S->y[i])/2; r0 += (ein->x[i+1]-ein->x[i])*(ein->y[i+1]+ein->y[i])/2;
} }
return r0; return r0;
case 11: case 11:
r0 = 0; r1 = 0; r0 = 0; r1 = 0;
for (i=0; i<n-1; i++) { for (i=0; i<n-1; i++) {
if (S->x[i+1]<=S->x[i]) if (ein->x[i+1]<=ein->x[i])
throw string( "x not sorted" ); throw string( "x not sorted" );
r0 += (S->x[i+1]-S->x[i])*(S->y[i+1]+S->y[i])/2; r0 += (ein->x[i+1]-ein->x[i])*(ein->y[i+1]+ein->y[i])/2;
} }
if (r0<=0) if (r0<=0)
throw string( "int <= 0" ); throw string( "int <= 0" );
for (i=0; i<n-1; i++) { for (i=0; i<n-1; i++) {
dr = (S->x[i+1]-S->x[i])*(S->y[i+1]+S->y[i])/2; dr = (ein->x[i+1]-ein->x[i])*(ein->y[i+1]+ein->y[i])/2;
if ((r1+dr)/r0>=arg2) break; if ((r1+dr)/r0>=arg2) break;
r1 += dr; r1 += dr;
} }
area = arg2*r0 - r1; area = arg2*r0 - r1;
triangle = (S->y[i+1]-S->y[i])/(S->x[i+1]-S->x[i])/2; triangle = (ein->y[i+1]-ein->y[i])/(ein->x[i+1]-ein->x[i])/2;
if (triangle==0) { if (triangle==0) {
return S->x[i]+area/dr; return ein->x[i]+area/dr;
} else { } else {
p = S->y[i]/2/triangle; p = ein->y[i]/2/triangle;
return S->x[i] - p + return ein->x[i] - p +
(triangle>0?+1:-1) * sqrt(p*p + area/triangle); (triangle>0?+1:-1) * sqrt(p*p + area/triangle);
} }
case 13: case 13:
for (i=n-1; i>=0; i--) for (i=n-1; i>=0; i--)
if (S->y[i]>0) break; if (ein->y[i]>0) break;
return S->x[i]; return ein->x[i];
case 15: // center of gravity case 15: // center of gravity
r0 = 0; r1 = 0; r0 = 0; r1 = 0;
for (i=0; i<n; i++) { for (i=0; i<n; i++) {
r0 += S->y[i]; r0 += ein->y[i];
r1 += S->y[i] * S->x[i]; r1 += ein->y[i] * ein->x[i];
} }
return r0>0 ? r1/r0 : 0; return r0>0 ? r1/r0 : 0;
case 16: // stddev in x case 16: // stddev in x
r0 = 0; r1 = 0; r2 = 0; r0 = 0; r1 = 0; r2 = 0;
for (i=0; i<n; i++) { for (i=0; i<n; i++) {
r0 += S->y[i]; r0 += ein->y[i];
r1 += S->y[i] * S->x[i]; r1 += ein->y[i] * ein->x[i];
} }
if( r0<=0 ) if( r0<=0 )
throw string( "r0 <= 0" ); throw string( "r0 <= 0" );
center = r1/r0; center = r1/r0;
for (i=0; i<n; i++) { for (i=0; i<n; i++) {
r2 += S->y[i] * SQR(S->x[i]-center); r2 += ein->y[i] * SQR(ein->x[i]-center);
} }
return sqrt( r2/r0 ); return sqrt( r2/r0 );
case 17: // correlation coefficient case 17: // correlation coefficient
r0 = 0; r1 = 0; s0 = 0; s1 = 0; v = 0; r0 = 0; r1 = 0; s0 = 0; s1 = 0; v = 0;
for (i=0; i<n; i++) { for (i=0; i<n; i++) {
r0 += S->x[i] / n; r0 += ein->x[i] / n;
r1 += S->y[i] / n; r1 += ein->y[i] / n;
} }
for (i=0; i<n; i++) { for (i=0; i<n; i++) {
s0 += SQR(S->x[i]-r0) / n; s0 += SQR(ein->x[i]-r0) / n;
s1 += SQR(S->y[i]-r1) / n; s1 += SQR(ein->y[i]-r1) / n;
v += (S->x[i]-r0) * (S->y[i]-r1) / n; v += (ein->x[i]-r0) * (ein->y[i]-r1) / n;
} }
return v / sqrt(s0*s1); return v / sqrt(s0*s1);
default: default:
throw string( "PROGRAM ERROR Integral not implemented" ); throw string( "BUG: this oio mode not implemented" );
} }
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment