// This file is part of fityk program. Copyright (C) Marcin Wojdyr
// Licence: GNU General Public License version 2
// $Id: datatrans.cpp 288 2007-04-22 02:14:02Z wojdyr $
// the idea of VM is based on one of boost::spirit samples - vmachine_calc
// this file can be compiled to stand-alone test program:
// $ g++ -I../3rdparty -DSTANDALONE_DATATRANS datatrans*.cpp numfuncs.cpp -o dt
// $ ./dt
// right hand variables:
// arrays (index in square brackets can be ommitted and n is assumed):
// x[] - x coordinate of point before transformation (old x)
// X[] - x coordinate after transformation (new x),
// y[], Y[] - y coordinate of point (before/after transformation)
// s[], S[] - std. dev. of y,
// a[], A[] - active point, either 0 or 1
// scalars:
// M - size of arrays,
// n - current index (index of currently transformed point)
// upper-case variables are assignable
// assignment syntax:
// or: assignable_var = expression --> executed for all points
// or: assignable_var[k] = expression --> executed for k-th point
// or: assignable_var[k...m] = expression
// or: assignable_var[k...] = expression
// or: assignable_var[...k] = expression
// where k and m are integers (negative are counted from the end:
// -1 is last, -2 one before last, and so on)
// FIXME: should x[-1] also mean x[M-1]?
// k...m means that expression is executed for k, k+1, ... m-1. Not for m.
//
// Assignments are executed for n = 0, ..., M-1 (if assignments are joined
// with '&', all are executed for n=0, then all for n=1, and so on).
// Before transformation the new point coordinates are a copy of the old ones.
//
// There is a special function sum(), which can be used as a real value.
// and returns a sum over all points of value of expression
// eg. sum(n > 0 ? (x[n] - x[n-1]) * (y[n-1] + y[n])/2 : 0) -> area.
// The value of sum() is computed before transformation.
//
//
// There are statements, that are executed only once:
// M=...,
// order=...
// and assignments T[k]=..., if these statements are joined with other
// assignments using ',', they are executed first, before iteration over all
// indices.
//
// M=length (eg. M=M+5 or M=100) changes number of points:
// either appends new data points with x=y=sigma=0, is_active=false,
// or deletes last points.
// order=[-]coordinate (eg. order=x, order=-x)
// sort data points; sorting is stable. After finishing transform
// command, points are sorted using x.
// delete[k] (or delete[k...m])
// delete(condition)
// deletion and change of M is postponed to the end of computation
//
// operators: binary +, -, *, /, % , ^
// ternary operator: condition ? expression1 : expression2
//
//functions: sqrt exp log10 ln sin cos tan atan asin acos abs min2 max2 round
// boolean: AND, OR, >, >=, <, <=, = (or ==), != (or <>), TRUE, FALSE, NOT
//
//parametrized functions: spline, interpolate //TODO->polyline
// The general syntax is: pfunc[param1, param2,...](expression),
// eg. spline[22.1, 37.9, 48.1, 17.2, 93.0, 20.7](x)
// spline - cubic spline interpolation, parameters are x1, y1, x2, y2, ...
// interpolate - polyline interpolation, parameters are x1, y1, x2, y2, ...
//
// All computations are performed using real numbers, but using round for
// comparisions should not be neccessary. Two numbers that differ less
// than epsilon (see: common.h) ie. abs(a-b)<epsilon, are considered equal.
// Indices are also computed in real number domain,
// and if they are not integers, interpolation of two values
// is taken (i.e of values with indices floor(idx) and ceil(idx)
//
//
// Syntax examples:
// set standard deviation as sqrt(max2(1,y))
// s = sqrt(max2(1,y))
// or (the same): s = sqrt(max2(1,y[n]))
// or: S = sqrt(max2(1,y[n]))
// or: S = sqrt(max2(1,Y[n]))
// or: S = sqrt(max2(1,Y))
//
// integration: Y[1...] = Y[n-1] + y[n]
//
// swaping x and y axes: y=x , x=y , s=sqrt(Y)
// or: y=x , x=y , s=sqrt(x)
//
// smoothing: Y[1...-1] = (y[n-1] + y[n+1])/2
//
// reducing: order = x, # not neccessary, points are always in this order
// x[...-1] = (x[n]+x[n+1])/2,
// y[...-1] = y[n]+y[n+1],
// delete(n%2==1)
//
// delete inactive points: delete(not a)
// normalize area:
// y = y / sum(n > 0 ? (x[n] - x[n-1]) * (y[n-1] + y[n])/2 : 0)
//
// substract spline baseline given by points (22.17, 37.92), (48.06, 17.23),
// (93.03, 20.68).
// y = y - spline[22.17 37.92 48.06 17.23 93.03 20.68](x)
//
//-----------------------------------------------------------------
// how it works:
// * parse string and prepare VM code (list of operators)
// and data (list of numbers). There is one VM code and
// data array for both do-once operations
// for for-many-points operations (perhaps it should be separated).
//
// * execute do-once operations (eg. order=x, x[3]=4)
// and computes value of sum()
//
// * execute all assignments for every point (from first to last),
// unless the point is outside of specified range.
//
//#define BOOST_SPIRIT_DEBUG
//#include "datatrans.h"
#include "datatrans2.h"
#include "sum.h"
#include "voigt.h"
//#include "common.h"
//#include "data.h"
//#include "var.h"
//#include "numfuncs.h"
//#include "logic.h"
//#include <boost/spirit/core.hpp>
//using namespace std;
//using namespace boost::spirit;
using namespace datatrans;
//#include "ui.h"
//#define DT_DEBUG(x) msg(x);
#ifdef STANDALONE_DATATRANS
#include <iostream>
bool dt_verbose = false;
#define DT_DEBUG(x) if (dt_verbose) std::cout << (x) << std::endl;
//------------------------- Main program ---------------------------------
#include <string.h>
int main(int argc, char **argv)
{
cout << "DataTransformVMachine test started...\n";
if (argc > 1 && !strcmp(argv[1], "-v")) {
cout << "[verbose mode]" << endl;
dt_verbose=true;
}
cout << "==> ";
string str;
while (getline(cin, str))
{
vector<Point> points;
points.push_back(Point(0.1, 4, 1));
points.push_back(Point(0.2, 2, 1));
points.push_back(Point(0.3, -2, 1));
points.push_back(Point(0.4, 3, 1));
points.push_back(Point(0.5, 3, 1));
try {
vector<Point> transformed_points = transform_data(str, points);
int len = (int) points.size();
for (int n = 0; n != (int) transformed_points.size(); n++)
cout << "point " << n << ": "
<< (n < len ? points[n] : Point()).str()
<< " -> " << transformed_points[n].str() << endl;
} catch (ExecuteError& e) {
cout << e.what() << endl;
}
cout << "==> ";
}
return 0;
}
#else
# ifndef DT_DEBUG
# define DT_DEBUG(x) ;
# endif
#endif //STANDALONE_DATATRANS
//---------------------------- grammar ----------------------------------
template <typename ScannerT>
DataTransformGrammar::definition<ScannerT>::definition(
DataTransformGrammar const& /*self*/)
{
range
= '[' >>
( eps_p(int_p >> ']')[push_op(OP_DO_ONCE)]
>> int_p [push_double()] [push_op(OP_INDEX)]
| (
ch_p('n') [push_the_double(0.)] [push_the_double(0.)]
| (
( int_p [push_double()]
| eps_p [push_the_double(0.)]
)
>> (str_p("...") | "..")
>> ( int_p [push_double()]
| eps_p [push_the_double(0)]
)
)
) [push_op(OP_RANGE)]
)
>> ']'
;
order
= ('-' >> as_lower_d["x"]) [push_the_double(-1.)]
| (!ch_p('+') >> as_lower_d["x"]) [push_the_double(+1.)]
| ('-' >> as_lower_d["y"]) [push_the_double(-2.)]
| (!ch_p('+') >> as_lower_d["y"]) [push_the_double(+2.)]
| ('-' >> as_lower_d["s"]) [push_the_double(-3.)]
| (!ch_p('+') >> as_lower_d["s"]) [push_the_double(+3.)]
;
assignment //not only assignments
= (as_lower_d["x"] >> !range >> '=' >> DataExpressionG)
[push_op(OP_ASSIGN_X)]
| (as_lower_d["y"] >> !range >> '=' >> DataExpressionG)
[push_op(OP_ASSIGN_Y)]
| (as_lower_d["s"] >> !range >> '=' >> DataExpressionG)
[push_op(OP_ASSIGN_S)]
| (as_lower_d["a"] >> !range >> '=' >> DataExpressionG)
[push_op(OP_ASSIGN_A)]
| ((ch_p('M') >> '=') [push_op(OP_DO_ONCE)]
>> DataExpressionG) [push_op(OP_RESIZE)]
| ((as_lower_d["order"] >> '=') [push_op(OP_DO_ONCE)]
>> order) [push_op(OP_ORDER)]
| (as_lower_d["delete"] >> eps_p('[') [push_op(OP_DO_ONCE)]
>> range) [push_op(OP_DELETE)]
| (as_lower_d["delete"] >> '(' >> DataExpressionG >> ')')
[push_op(OP_DELETE_COND)]
;
statement
= (eps_p[push_op(OP_BEGIN)] >> assignment >> eps_p[push_op(OP_END)])
% ch_p(',')
;
}
// explicit template instantiations
template DataTransformGrammar::definition<scanner<char const*, scanner_policies<skipper_iteration_policy<iteration_policy>, match_policy, no_actions_action_policy<action_policy> > > >::definition(DataTransformGrammar const&);
template DataTransformGrammar::definition<scanner<char const*, scanner_policies<skipper_iteration_policy<iteration_policy>, match_policy, no_actions_action_policy<no_actions_action_policy<action_policy> > > > >::definition(DataTransformGrammar const&);
DataTransformGrammar DataTransformG;
namespace datatrans {
/// debuging utility
#define OP_(x) \
if (op == OP_##x) return #x;
string dt_op(int op)
{
OP_(NEG) OP_(EXP) OP_(SIN) OP_(COS) OP_(ATAN) OP_(ABS) OP_(ROUND)
OP_(TAN) OP_(ASIN) OP_(ACOS)
OP_(LOG10) OP_(LN) OP_(SQRT) OP_(POW)
OP_(GAMMA) OP_(LGAMMA) OP_(VOIGT)
OP_(ADD) OP_(SUB) OP_(MUL) OP_(DIV) OP_(MOD)
OP_(MIN2) OP_(MAX2) OP_(RANDNORM) OP_(RANDU)
OP_(VAR_X) OP_(VAR_Y) OP_(VAR_S) OP_(VAR_A)
OP_(VAR_x) OP_(VAR_y) OP_(VAR_s) OP_(VAR_a)
OP_(VAR_n) OP_(VAR_M) OP_(NUMBER)
OP_(OR) OP_(AFTER_OR) OP_(AND) OP_(AFTER_AND) OP_(NOT)
OP_(TERNARY) OP_(TERNARY_MID) OP_(AFTER_TERNARY) OP_(DELETE_COND)
OP_(GT) OP_(GE) OP_(LT) OP_(LE) OP_(EQ) OP_(NEQ) OP_(NCMP_HACK)
OP_(RANGE) OP_(INDEX) OP_(x_IDX)
OP_(ASSIGN_X) OP_(ASSIGN_Y) OP_(ASSIGN_S) OP_(ASSIGN_A)
OP_(DO_ONCE) OP_(RESIZE) OP_(ORDER) OP_(DELETE) OP_(BEGIN) OP_(END)
OP_(END_AGGREGATE) OP_(AGCONDITION)
OP_(AGSUM) OP_(AGMIN) OP_(AGMAX) OP_(AGAREA) OP_(AGAVG) OP_(AGSTDDEV)
OP_(PARAMETERIZED) OP_(PLIST_BEGIN) OP_(PLIST_SEP) OP_(PLIST_END)
OP_(FUNC) OP_(SUM_F) OP_(SUM_Z) OP_(NUMAREA) OP_(FINDX) OP_(FIND_EXTR)
return S(op);
};
/// debuging utility
string dt_ops(vector<int> const& code)
{
string r;
for (std::vector<int>::const_iterator i=code.begin(); i != code.end(); ++i)
if (*i < 0)
r += dt_op(*i) + " ";
else
r += "[" + S(*i) + "] ";
return r;
}
//------------------------ Virtual Machine --------------------------------
vector<int> code; // VM code
vector<fp> numbers; // VM data (numeric values)
vector<ParameterizedFunction*> parameterized; // also used by VM
const int stack_size = 128; //should be enough,
//there are no checks for stack overflow
void clear_parse_vecs()
{
code.clear();
numbers.clear();
//TODO shared_ptr
purge_all_elements(parameterized);
}
vector<int>::const_iterator
skip_code(vector<int>::const_iterator i, int start_op, int finish_op)
{
int counter = 1;
while (counter) {
++i;
if (*i == finish_op) counter--;
else if (*i == start_op) counter++;
}
return i;
}
void skip_to_end(vector<int>::const_iterator &i)
{
DT_DEBUG("SKIPing")
while (*i != OP_END)
++i;
}
vector<int>::const_iterator find_end(vector<int>::const_iterator i)
{
while (*i != OP_END)
++i;
return i;
}
template<typename T>
fp get_var_with_idx(fp idx, vector<Point> const& points, T Point::*t)
{
if (idx < 0 || idx+1 > points.size())
return 0.;
else if (is_eq(idx, iround(idx)))
return points[iround(idx)].*t;
else {
int flo = int(floor(idx));
fp fra = idx - flo;
return (1-fra) * fp(points[flo].*t) + fra * fp(points[flo+1].*t);
}
}
/// returns floating-point "index" of x in sorted vec of points
fp find_idx_in_sorted(vector<Point> const& pp, fp x)
{
if (pp.empty())
return 0;
else if (x <= pp.front().x)
return 0;
else if (x >= pp.back().x)
return pp.size() - 1;
vector<Point>::const_iterator i = lower_bound(pp.begin(), pp.end(),
Point(x, 0));
assert (i > pp.begin() && i < pp.end());
if (is_eq(x, i->x))
return i - pp.begin();
else
return i - pp.begin() - (i->x - x) / (i->x - (i-1)->x);
}
// Return value:
// if once==true (one-time pass).
// true means: it is neccessary to call this function for every point.
// if iterative pass:
// true means: delete this point.
// n: index of point; // special value: n==M: one-time operations
// n(==M) can be changed in OP_INDEX
// M: number of all points (==new_points.size())
bool execute_code(int n, int &M, vector<fp>& stack,
vector<Point> const& old_points, vector<Point>& new_points,
vector<int> const& code)
{
DT_DEBUG("executing code; n=" + S(n) + " M=" + S(M))
assert(M == size(new_points));
bool once = (n == M);
bool return_value = false;
vector<fp>::iterator stackPtr = stack.begin() - 1;//will be ++'ed first
#define STACK_OP
#define STACK_OFFSET_CHANGE(ch) stackPtr+=(ch);
for (vector<int>::const_iterator i=code.begin(); i != code.end(); i++) {
//---8<--- START OF DT COMMON BLOCK --------------------------------------
DT_DEBUG("op " + dt_op(*i))
switch (*i) {
//unary-operators
case OP_NEG:
STACK_OP *stackPtr = - *stackPtr;
break;
case OP_SQRT:
STACK_OP *stackPtr = sqrt(*stackPtr);
break;
case OP_GAMMA:
STACK_OP *stackPtr = gammafn(*stackPtr);
break;
case OP_LGAMMA:
STACK_OP *stackPtr = lgammafn(*stackPtr);
break;
case OP_EXP:
STACK_OP *stackPtr = exp(*stackPtr);
break;
case OP_ERFC:
STACK_OP *stackPtr = erfc(*stackPtr);
break;
case OP_ERF:
STACK_OP *stackPtr = erf(*stackPtr);
break;
case OP_LOG10:
STACK_OP *stackPtr = log10(*stackPtr);
break;
case OP_LN:
STACK_OP *stackPtr = log(*stackPtr);
break;
case OP_SIN:
STACK_OP *stackPtr = sin(*stackPtr);
break;
case OP_COS:
STACK_OP *stackPtr = cos(*stackPtr);
break;
case OP_TAN:
STACK_OP *stackPtr = tan(*stackPtr);
break;
case OP_ATAN:
STACK_OP *stackPtr = atan(*stackPtr);
break;
case OP_ASIN:
STACK_OP *stackPtr = asin(*stackPtr);
break;
case OP_ACOS:
STACK_OP *stackPtr = acos(*stackPtr);
break;
case OP_ABS:
STACK_OP *stackPtr = fabs(*stackPtr);
break;
case OP_ROUND:
STACK_OP *stackPtr = floor(*stackPtr + 0.5);
break;
case OP_x_IDX:
STACK_OP *stackPtr = find_idx_in_sorted(old_points, *stackPtr);
break;
#ifndef STANDALONE_DATATRANS
case OP_FUNC:
i++;
STACK_OP
*stackPtr = AL->get_function(*i)->calculate_value(*stackPtr);
break;
case OP_SUM_F:
i++;
STACK_OP *stackPtr = AL->get_sum(*i)->value(*stackPtr);
break;
case OP_SUM_Z:
i++;
STACK_OP *stackPtr = AL->get_sum(*i)->zero_shift(*stackPtr);
break;
case OP_NUMAREA:
i += 2;
STACK_OFFSET_CHANGE(-2)
if (*(i-1) == OP_FUNC) {
STACK_OP
*stackPtr = AL->get_function(*i)->numarea(*stackPtr,
*(stackPtr+1), iround(*(stackPtr+2)));
}
else if (*(i-1) == OP_SUM_F) {
STACK_OP
*stackPtr = AL->get_sum(*i)->numarea(*stackPtr,
*(stackPtr+1), iround(*(stackPtr+2)));
}
else // OP_SUM_Z
throw ExecuteError("numarea(Z,...) is not implemented."
"Does anyone need it?");
break;
case OP_FINDX:
i += 2;
STACK_OFFSET_CHANGE(-2)
if (*(i-1) == OP_FUNC) {
STACK_OP
*stackPtr = AL->get_function(*i)->find_x_with_value(
*stackPtr, *(stackPtr+1), *(stackPtr+2));
}
else if (*(i-1) == OP_SUM_F) {
throw ExecuteError("findx(F,...) is not implemented. "
"Does anyone need it?");
}
else // OP_SUM_Z
throw ExecuteError("findx(Z,...) is not implemented. "
"Does anyone need it?");
break;
case OP_FIND_EXTR:
i += 2;
STACK_OFFSET_CHANGE(-1)
if (*(i-1) == OP_FUNC) {
STACK_OP
*stackPtr = AL->get_function(*i)->find_extremum(*stackPtr,
*(stackPtr+1));
}
else if (*(i-1) == OP_SUM_F) {
throw ExecuteError("extremum(F,...) is not implemented. "
"Does anyone need it?");
}
else // OP_SUM_Z
throw ExecuteError("extremum(Z,...) is not implemented. "
"Does anyone need it?");
break;
#endif //not STANDALONE_DATATRANS
case OP_PARAMETERIZED:
i++;
STACK_OP *stackPtr = parameterized[*i]->calculate(*stackPtr);
break;
//binary-operators
case OP_MIN2:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = min(*stackPtr, *(stackPtr+1));
break;
case OP_MAX2:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = max(*stackPtr, *(stackPtr+1));
break;
case OP_RANDU:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = rand_uniform(*stackPtr, *(stackPtr+1));
break;
case OP_RANDNORM:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr += rand_gauss() * *(stackPtr+1);
break;
case OP_ADD:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr += *(stackPtr+1);
break;
case OP_SUB:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr -= *(stackPtr+1);
break;
case OP_MUL:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr *= *(stackPtr+1);
break;
case OP_DIV:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr /= *(stackPtr+1);
break;
case OP_MOD:
STACK_OFFSET_CHANGE(-1)
STACK_OP
*stackPtr -= floor(*stackPtr / *(stackPtr+1)) * *(stackPtr+1);
break;
case OP_POW:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = pow(*stackPtr, *(stackPtr+1));
break;
case OP_VOIGT:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = humlik(*stackPtr, *(stackPtr+1))
/ sqrt(M_PI);
break;
// comparisions
case OP_LT:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = is_lt(*stackPtr, *(stackPtr+1));
break;
case OP_GT:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = is_gt(*stackPtr, *(stackPtr+1));
break;
case OP_LE:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = is_le(*stackPtr, *(stackPtr+1));
break;
case OP_GE:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = is_ge(*stackPtr, *(stackPtr+1));
break;
case OP_EQ:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = is_eq(*stackPtr, *(stackPtr+1));
break;
case OP_NEQ:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = is_neq(*stackPtr, *(stackPtr+1));
break;
// next comparision hack, see rbool rule for more...
case OP_NCMP_HACK:
STACK_OFFSET_CHANGE(+1)
// put number that is accidentally in unused part of the stack
STACK_OP *stackPtr = *(stackPtr+1);
break;
// putting-number-to-stack-operators
// stack overflow not checked
case OP_NUMBER:
STACK_OFFSET_CHANGE(+1)
i++;
STACK_OP *stackPtr = numbers[*i];
break;
case OP_VAR_n:
STACK_OFFSET_CHANGE(+1)
STACK_OP *stackPtr = static_cast<fp>(n);
break;
case OP_VAR_M:
STACK_OFFSET_CHANGE(+1)
STACK_OP *stackPtr = static_cast<fp>(new_points.size());
break;
case OP_VAR_x:
STACK_OP
*stackPtr = get_var_with_idx(*stackPtr, old_points, &Point::x);
break;
case OP_VAR_y:
STACK_OP
*stackPtr = get_var_with_idx(*stackPtr, old_points, &Point::y);
break;
case OP_VAR_s:
STACK_OP
*stackPtr = get_var_with_idx(*stackPtr, old_points,
&Point::sigma);
break;
case OP_VAR_a:
STACK_OP
*stackPtr = bool(iround(get_var_with_idx(*stackPtr, old_points,
&Point::is_active)));
break;
case OP_VAR_X:
STACK_OP
*stackPtr = get_var_with_idx(*stackPtr, new_points, &Point::x);
break;
case OP_VAR_Y:
STACK_OP
*stackPtr = get_var_with_idx(*stackPtr, new_points, &Point::y);
break;
case OP_VAR_S:
STACK_OP
*stackPtr = get_var_with_idx(*stackPtr, new_points,
&Point::sigma);
break;
case OP_VAR_A:
STACK_OP
*stackPtr = bool(iround(get_var_with_idx(*stackPtr, new_points,
&Point::is_active)));
break;
//assignment-operators
case OP_ASSIGN_X:
STACK_OP new_points[n].x = *stackPtr;
STACK_OFFSET_CHANGE(-1)
break;
case OP_ASSIGN_Y:
STACK_OP new_points[n].y = *stackPtr;
STACK_OFFSET_CHANGE(-1)
break;
case OP_ASSIGN_S:
STACK_OP new_points[n].sigma = *stackPtr;
STACK_OFFSET_CHANGE(-1)
break;
case OP_ASSIGN_A:
STACK_OP new_points[n].is_active = is_neq(*stackPtr, 0.);
STACK_OFFSET_CHANGE(-1)
break;
// logical; can skip part of VM code !
case OP_NOT:
STACK_OP *stackPtr = is_eq(*stackPtr, 0.);
break;
//---8<--- END OF DT COMMON BLOCK --------------------------------------
case OP_AND:
if (is_neq(*stackPtr, 0)) //return second
stackPtr--;
else // return first
i = skip_code(i, OP_AND, OP_AFTER_AND);
break;
case OP_OR:
if (is_neq(*stackPtr, 0)) //return first
i = skip_code(i, OP_OR, OP_AFTER_OR);
else // return second
stackPtr--;
break;
case OP_TERNARY:
if (! *stackPtr)
i = skip_code(i, OP_TERNARY, OP_TERNARY_MID);
stackPtr--;
break;
case OP_TERNARY_MID:
//if we are here, condition was true. Skip.
i = skip_code(i, OP_TERNARY_MID, OP_AFTER_TERNARY);
break;
case OP_AFTER_AND: //do nothing
case OP_AFTER_OR:
case OP_AFTER_TERNARY:
break;
case OP_DELETE_COND:
if (!is_zero(*stackPtr))
return_value = true;
stackPtr--;
break;
//transformation condition
case OP_INDEX:
assert(once); //x[n]= or delete[n]
n = iround(*stackPtr);//changing n(!!) for use in OP_ASSIGN_
stackPtr--;
if (n < 0)
n += M;
if (n < 0 || n >= M)
skip_to_end(i);
if (*(i+1) == OP_DELETE) {
new_points.erase(new_points.begin() + n);
M--;
skip_to_end(i);
}
break;
case OP_RANGE:
{
//x[i...j]= or delete[i...j]
int right = iround(*stackPtr); //Last In First Out
stackPtr--;
if (right <= 0)
right += M;
int left = iround(*stackPtr);
stackPtr--;
if (left < 0)
left += M;
if (*(i+1) == OP_DELETE) {
if (0 < left && left < right && right <= M) {
new_points.erase(new_points.begin()+left,
new_points.begin()+right);
M = size(new_points);
}
skip_to_end(i);
}
else {
//if n not in [i...j] then skip to prevent OP_ASSIGN_.
bool n_between = (left <= n && n < right);
if (!n_between)
skip_to_end(i);
}
break;
}
case OP_BEGIN:
{
bool next_op_once = (*(i+1) == OP_DO_ONCE);
if (next_op_once != once) {
skip_to_end(i);
if (once)
return_value=true;
}
break;
}
case OP_END: //nothing -- it is used only for skipping assignment
break;
// once - operators
case OP_DO_ONCE: //do nothing, it is only a flag
assert(once);
break;
case OP_RESIZE:
assert(once);
M = iround(*stackPtr);
stackPtr--;
new_points.resize(M);
break;
case OP_ORDER:
{
assert(once);
int ord = iround(*stackPtr);
stackPtr--;
DT_DEBUG("in OP_ORDER with " + S(ord))
if (ord == 1) {
DT_DEBUG("sort x_lt")
stable_sort(new_points.begin(), new_points.end(), x_lt);
}
else if(ord == -1) {
DT_DEBUG("sort x_gt")
stable_sort(new_points.begin(), new_points.end(), x_gt);
}
else if(ord == 2) {
DT_DEBUG("sort y_lt")
stable_sort(new_points.begin(), new_points.end(), y_lt);
}
else if(ord == -2) {
DT_DEBUG("sort y_gt")
stable_sort(new_points.begin(), new_points.end(), y_gt);
}
break;
}
case OP_DELETE:
assert(0); //OP_DELETE is processed in OP_INDEX or OP_RANGE
break;
default:
DT_DEBUG("Unknown operator in VM code: " + S(*i))
assert(0);
}
}
//assert(stackPtr == stack.begin() - 1 //DataTransformGrammar
// || (stackPtr == stack.begin() && once)); //DataExpressionGrammar
return return_value;
}
#if 0
// not tested !!!!!!
// this function could be used for optimization of data transformations
// but optimization is not neccessary (?)
void multipoint_execute_code(int M, vector<fp>& stack,
vector<Point> const& old_points, vector<Point>& new_points,
vector<int> const& code)
{
DT_DEBUG("executing code (multipoint), M=" + S(M))
assert(M == size(new_points));
int offset = -1; //stack offset, will be ++'ed first
vector<vector<int>::const_iterator> skips(M, code.begin());
#undef STACK_OP
#define STACK_OP \
for (vector<fp>::iterator stackPtr = stack.begin() + M * offset; \
stackPtr!=stack.end(); stackPtr=stack.end()) \
for (int n = 0; n < M; ++n, ++stackPtr) \
// if (i >= skips[n])
#undef STACK_OFFSET_CHANGE
#define STACK_OFFSET_CHANGE(ch) \
offset += (ch);
for (vector<int>::const_iterator i=code.begin(); i != code.end(); i++) {
if (size(stack) < (offset+2) * M)
stack.resize((offset+2) * M);
// the exact copy ...
//---8<--- START OF DT COMMON BLOCK --------------------------------------
DT_DEBUG("op " + dt_op(*i))
switch (*i) {
//unary-operators
case OP_NEG:
STACK_OP *stackPtr = - *stackPtr;
break;
case OP_SQRT:
STACK_OP *stackPtr = sqrt(*stackPtr);
break;
case OP_GAMMA:
STACK_OP *stackPtr = gammafn(*stackPtr);
break;
case OP_LGAMMA:
STACK_OP *stackPtr = lgammafn(*stackPtr);
break;
case OP_EXP:
STACK_OP *stackPtr = exp(*stackPtr);
break;
case OP_ERFC:
STACK_OP *stackPtr = erfc(*stackPtr);
break;
case OP_ERF:
STACK_OP *stackPtr = erf(*stackPtr);
break;
case OP_LOG10:
STACK_OP *stackPtr = log10(*stackPtr);
break;
case OP_LN:
STACK_OP *stackPtr = log(*stackPtr);
break;
case OP_SIN:
STACK_OP *stackPtr = sin(*stackPtr);
break;
case OP_COS:
STACK_OP *stackPtr = cos(*stackPtr);
break;
case OP_TAN:
STACK_OP *stackPtr = tan(*stackPtr);
break;
case OP_ATAN:
STACK_OP *stackPtr = atan(*stackPtr);
break;
case OP_ASIN:
STACK_OP *stackPtr = asin(*stackPtr);
break;
case OP_ACOS:
STACK_OP *stackPtr = acos(*stackPtr);
break;
case OP_ABS:
STACK_OP *stackPtr = fabs(*stackPtr);
break;
case OP_ROUND:
STACK_OP *stackPtr = floor(*stackPtr + 0.5);
break;
case OP_x_IDX:
STACK_OP *stackPtr = find_idx_in_sorted(old_points, *stackPtr);
break;
#ifndef STANDALONE_DATATRANS
case OP_FUNC:
i++;
STACK_OP
*stackPtr = AL->get_function(*i)->calculate_value(*stackPtr);
break;
case OP_SUM_F:
i++;
STACK_OP *stackPtr = AL->get_sum(*i)->value(*stackPtr);
break;
case OP_SUM_Z:
i++;
STACK_OP *stackPtr = AL->get_sum(*i)->zero_shift(*stackPtr);
break;
case OP_NUMAREA:
i += 2;
STACK_OFFSET_CHANGE(-2)
if (*(i-1) == OP_FUNC) {
STACK_OP
*stackPtr = AL->get_function(*i)->numarea(*stackPtr,
*(stackPtr+1), iround(*(stackPtr+2)));
}
else if (*(i-1) == OP_SUM_F) {
STACK_OP
*stackPtr = AL->get_sum(*i)->numarea(*stackPtr,
*(stackPtr+1), iround(*(stackPtr+2)));
}
else // OP_SUM_Z
throw ExecuteError("numarea(Z,...) is not implemented."
"Does anyone need it?");
break;
case OP_FINDX:
i += 2;
STACK_OFFSET_CHANGE(-2)
if (*(i-1) == OP_FUNC) {
STACK_OP
*stackPtr = AL->get_function(*i)->find_x_with_value(
*stackPtr, *(stackPtr+1), *(stackPtr+2));
}
else if (*(i-1) == OP_SUM_F) {
throw ExecuteError("findx(F,...) is not implemented. "
"Does anyone need it?");
}
else // OP_SUM_Z
throw ExecuteError("findx(Z,...) is not implemented. "
"Does anyone need it?");
break;
case OP_FIND_EXTR:
i += 2;
STACK_OFFSET_CHANGE(-1)
if (*(i-1) == OP_FUNC) {
STACK_OP
*stackPtr = AL->get_function(*i)->find_extremum(*stackPtr,
*(stackPtr+1));
}
else if (*(i-1) == OP_SUM_F) {
throw ExecuteError("extremum(F,...) is not implemented. "
"Does anyone need it?");
}
else // OP_SUM_Z
throw ExecuteError("extremum(Z,...) is not implemented. "
"Does anyone need it?");
break;
#endif //not STANDALONE_DATATRANS
case OP_PARAMETERIZED:
i++;
STACK_OP *stackPtr = parameterized[*i]->calculate(*stackPtr);
break;
//binary-operators
case OP_MIN2:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = min(*stackPtr, *(stackPtr+1));
break;
case OP_MAX2:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = max(*stackPtr, *(stackPtr+1));
break;
case OP_RANDU:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = rand_uniform(*stackPtr, *(stackPtr+1));
break;
case OP_RANDNORM:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr += rand_gauss() * *(stackPtr+1);
break;
case OP_ADD:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr += *(stackPtr+1);
break;
case OP_SUB:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr -= *(stackPtr+1);
break;
case OP_MUL:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr *= *(stackPtr+1);
break;
case OP_DIV:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr /= *(stackPtr+1);
break;
case OP_MOD:
STACK_OFFSET_CHANGE(-1)
STACK_OP
*stackPtr -= floor(*stackPtr / *(stackPtr+1)) * *(stackPtr+1);
break;
case OP_POW:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = pow(*stackPtr, *(stackPtr+1));
break;
// comparisions
case OP_LT:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = is_lt(*stackPtr, *(stackPtr+1));
break;
case OP_GT:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = is_gt(*stackPtr, *(stackPtr+1));
break;
case OP_LE:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = is_le(*stackPtr, *(stackPtr+1));
break;
case OP_GE:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = is_ge(*stackPtr, *(stackPtr+1));
break;
case OP_EQ:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = is_eq(*stackPtr, *(stackPtr+1));
break;
case OP_NEQ:
STACK_OFFSET_CHANGE(-1)
STACK_OP *stackPtr = is_neq(*stackPtr, *(stackPtr+1));
break;
// next comparision hack, see rbool rule for more...
case OP_NCMP_HACK:
STACK_OFFSET_CHANGE(+1)
// put number that is accidentally in unused part of the stack
STACK_OP *stackPtr = *(stackPtr+1);
break;
// putting-number-to-stack-operators
// stack overflow not checked
case OP_NUMBER:
STACK_OFFSET_CHANGE(+1)
i++;
STACK_OP *stackPtr = numbers[*i];
break;
case OP_VAR_n:
STACK_OFFSET_CHANGE(+1)
STACK_OP *stackPtr = static_cast<fp>(n);
break;
case OP_VAR_M:
STACK_OFFSET_CHANGE(+1)
STACK_OP *stackPtr = static_cast<fp>(new_points.size());
break;
case OP_VAR_x:
STACK_OP
*stackPtr = get_var_with_idx(*stackPtr, old_points, &Point::x);
break;
case OP_VAR_y:
STACK_OP
*stackPtr = get_var_with_idx(*stackPtr, old_points, &Point::y);
break;
case OP_VAR_s:
STACK_OP
*stackPtr = get_var_with_idx(*stackPtr, old_points,
&Point::sigma);
break;
case OP_VAR_a:
STACK_OP
*stackPtr = bool(iround(get_var_with_idx(*stackPtr, old_points,
&Point::is_active)));
break;
case OP_VAR_X:
STACK_OP
*stackPtr = get_var_with_idx(*stackPtr, new_points, &Point::x);
break;
case OP_VAR_Y:
STACK_OP
*stackPtr = get_var_with_idx(*stackPtr, new_points, &Point::y);
break;
case OP_VAR_S:
STACK_OP
*stackPtr = get_var_with_idx(*stackPtr, new_points,
&Point::sigma);
break;
case OP_VAR_A:
STACK_OP
*stackPtr = bool(iround(get_var_with_idx(*stackPtr, new_points,
&Point::is_active)));
break;
//assignment-operators
case OP_ASSIGN_X:
STACK_OP new_points[n].x = *stackPtr;
STACK_OFFSET_CHANGE(-1)
break;
case OP_ASSIGN_Y:
STACK_OP new_points[n].y = *stackPtr;
STACK_OFFSET_CHANGE(-1)
break;
case OP_ASSIGN_S:
STACK_OP new_points[n].sigma = *stackPtr;
STACK_OFFSET_CHANGE(-1)
break;
case OP_ASSIGN_A:
STACK_OP new_points[n].is_active = is_neq(*stackPtr, 0.);
STACK_OFFSET_CHANGE(-1)
break;
// logical; can skip part of VM code !
case OP_NOT:
STACK_OP *stackPtr = is_eq(*stackPtr, 0.);
break;
//---8<--- END OF DT COMMON BLOCK --------------------------------------
case OP_AND:
STACK_OP
if (is_zero(*stackPtr))
skips[n] = skip_code(i, OP_AND, OP_AFTER_AND) + 1;
STACK_OFFSET_CHANGE(-1)
break;
case OP_OR:
STACK_OP
if (!is_zero(*stackPtr))
skips[n] = skip_code(i, OP_OR, OP_AFTER_OR) + 1;
STACK_OFFSET_CHANGE(-1)
break;
case OP_TERNARY:
STACK_OP
if (is_zero(*stackPtr))
skips[n] = skip_code(i, OP_TERNARY, OP_TERNARY_MID) + 1;
STACK_OFFSET_CHANGE(-1)
break;
case OP_TERNARY_MID:
//if we are here, condition was true. Skip.
STACK_OP
skips[n] = skip_code(i, OP_TERNARY_MID, OP_AFTER_TERNARY) + 1;
break;
case OP_AFTER_AND: //do nothing
case OP_AFTER_OR:
case OP_AFTER_TERNARY:
break;
case OP_DELETE_COND:
for (int n = M-1; n >= 0; --n)
if (!is_zero(stack[n + M * offset]))
new_points.erase(new_points.begin() + n);
offset--;
break;
//transformation condition
case OP_RANGE:
{
//x[i...j]= or delete[i...j]
STACK_OP
{
int right = iround(*stackPtr); //Last In First Out
if (right <= 0)
right += M;
int left = iround(*(stackPtr-1));
stackPtr-=2;
if (left < 0)
left += M;
assert (*(i+1) != OP_DELETE);
//if n not in [i...j] then skip to prevent OP_ASSIGN_.
bool n_between = (left <= n && n < right);
if (!n_between)
skips[n] = find_end(i) + 1;
}
break;
}
case OP_BEGIN:
if (*(i+1) == OP_DO_ONCE)
skip_to_end(i);
break;
case OP_END: //nothing -- it is used only for skipping assignment
break;
default:
DT_DEBUG("Unknown operator in VM code: " + S(*i))
assert(0);
}
}
}
#endif
/// change AGSUM X ... X END_AGGREGATE to NUMBER INDEX
void replace_aggregates(int M, vector<Point> const& old_points,
vector<int>& code, vector<int>::iterator cb)
{
vector<fp> stack(stack_size);
for (vector<int>::iterator i = cb; i != code.end(); ++i) {
if (*i == OP_AGMIN || *i == OP_AGMAX || *i == OP_AGSUM
|| *i == OP_AGAREA || *i == OP_AGAVG || *i == OP_AGSTDDEV) {
int op = *i;
vector<int>::iterator const start = i;
DT_DEBUG("code before replace: " + dt_ops(code));
replace_aggregates(M, old_points, code, start+1);
fp result = 0.;
fp mean = 0.; //needed for OP_AGSTDDEV
int counter = 0;
vector<Point> fake_new_points(M);
++i;
while (*i != OP_AGCONDITION && *i != OP_END_AGGREGATE)
++i;
vector<int> acode(start+1, i); //code for aggragate function
vector<int> ccode; //code for condition, empty if no condition
if (*i == OP_AGCONDITION) {
vector<int>::iterator const start_cond = i;
while (*i != OP_END_AGGREGATE)
++i;
ccode = vector<int>(start_cond+1, i);
}
for (int n = 0; n != M; n++) {
if (!ccode.empty()) {
execute_code(n,M,stack, old_points, fake_new_points, ccode);
if (is_eq(stack.front(), 0))
continue;
}
++counter;
execute_code(n, M, stack, old_points, fake_new_points, acode);
if (op == OP_AGSUM)
result += stack.front();
else if (op == OP_AGMIN) {
if (counter == 1)
result = stack.front();
else if (result > stack.front())
result = stack.front();
}
else if (op == OP_AGMAX) {
if (counter == 1)
result = stack.front();
else if (result < stack.front())
result = stack.front();
}
else if (op == OP_AGAREA) {
fp dx = (old_points[min(n+1, M-1)].x
- old_points[max(n-1, 0)].x) / 2.;
result += stack.front() * dx;
}
else if (op == OP_AGAVG) {
result += (stack.front() - result) / counter;
}
else if (op == OP_AGSTDDEV) {
// see: http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
fp x = stack.front();
fp delta = x - mean;
mean += delta / counter;
result += delta * (x - mean);
}
else
assert(0);
DT_DEBUG("n=" + S(n) + " stack.front() = " + S(stack.front()));
}
if (op == OP_AGSTDDEV)
result = sqrt(result / (counter - 1));
*start = OP_NUMBER;
*(start+1) = size(numbers);
numbers.push_back(result);
code.erase(start+2, i+1);
i = start+1;
DT_DEBUG("code after replace: " + dt_ops(code));
}
}
}
//-------------------------------------------------------------------------
void execute_vm_code(const vector<Point> &old_points, vector<Point> &new_points)
{
vector<fp> stack(stack_size);
int M = (int) new_points.size();
for (vector<ParameterizedFunction*>::iterator i = parameterized.begin();
i != parameterized.end(); ++i)
(*i)->prepare_parameters(old_points);
replace_aggregates(M, old_points, code, code.begin());
// first execute one-time operations: sorting, x[15]=3, etc.
// n==M => one-time op.
bool t = execute_code(M, M, stack, old_points, new_points, code);
if (!t)
return;
#if 1
vector<int> to_be_deleted;
for (int n = 0; n != M; n++) {
bool r = execute_code(n, M, stack, old_points, new_points, code);
if (r)
to_be_deleted.push_back(n);
}
if (!to_be_deleted.empty())
for (vector<int>::const_iterator i = to_be_deleted.end() - 1;
i >= to_be_deleted.begin(); --i)
new_points.erase(new_points.begin() + *i);
#else
//multipoint_execute_code(M, stack, old_points, new_points, code);
#endif
}
} //namespace
vector<Point> transform_data(string const& str, vector<Point> const& old_points)
{
clear_parse_vecs();
// First compile string...
parse_info<> result = parse(str.c_str(), DataTransformG, space_p);
if (!result.full) {
throw ExecuteError("Syntax error in data transformation formula.");
}
// and then execute compiled code.
vector<Point> new_points = old_points; //initial values of new_points
execute_vm_code(old_points, new_points);
return new_points;
}
bool validate_transformation(string const& str)
{
clear_parse_vecs();
parse_info<> result = parse(str.c_str(), DataTransformG, space_p);
return (bool) result.full;
}
bool validate_data_expression(string const& str)
{
clear_parse_vecs();
parse_info<> result = parse(str.c_str(), DataExpressionG, space_p);
return (bool) result.full;
}
bool is_data_dependent_code(vector<int> const& code)
{
for (vector<int>::const_iterator i = code.begin(); i != code.end(); ++i)
if ((*i >= OP_VAR_FIRST_OP && *i <= OP_VAR_LAST_OP)
|| *i == OP_END_AGGREGATE)
return true;
return false;
}
bool is_data_dependent_expression(string const& s)
{
if (!validate_data_expression(s)) //it fills `code'
return false;
return is_data_dependent_code(code);
}
fp get_transform_expression_value(string const &s, Data const* data)
{
clear_parse_vecs();
// First compile string... puts result into code etc.
parse_info<> result = parse(s.c_str(), DataExpressionG, space_p);
if (!result.full)
throw ExecuteError("Syntax error in expression: " + s);
if (!data && is_data_dependent_code(code))
throw ExecuteError("Dataset is not specified and the expression "
"depends on it: " + s);
vector<Point> const no_points;
return get_transform_expr_value(code, data ? data->points() : no_points);
}
fp get_transform_expr_value(vector<int>& code_, vector<Point> const& points)
{
static vector<fp> stack(stack_size);
int M = (int) points.size();
vector<Point> new_points = points;
for (vector<int>::const_iterator i = code_.begin(); i != code_.end(); ++i)
if (*i == OP_PARAMETERIZED)
parameterized[*(i+1)]->prepare_parameters(points);
replace_aggregates(M, points, code_, code_.begin());
// n==M => one-time op.
bool t = execute_code(M, M, stack, points, new_points, code_);
if (t)
throw ExecuteError("Expression depends on undefined `n' index.");
return stack.front();
}
bool get_dt_code(string const& s, vector<int>& code_, vector<fp>& numbers_)
{
clear_parse_vecs();
// First compile string... puts result into code etc.
parse_info<> result = parse(s.c_str(), DataExpressionG, space_p);
if (!result.full)
return false;
for (vector<int>::iterator i = code.begin(); i != code.end(); ++i)
if (*i == OP_PARAMETERIZED
|| *i == OP_AGMIN || *i == OP_AGMAX || *i == OP_AGSUM
|| *i == OP_AGAREA || *i == OP_AGAVG || *i == OP_AGSTDDEV)
return false;
code_ = code;
numbers_ = numbers;
return true;
}
fp get_value_for_point(vector<int> const& code_, vector<fp> const& numbers_,
fp x, fp y)
{
static vector<fp> stack(stack_size);
static vector<Point> points(1);
static vector<Point> new_points(1);
numbers = numbers_;
points[0] = new_points[0] = Point(x, y);
int M = 1;
execute_code(0, M, stack, points, new_points, code_);
return stack.front();
}
vector<fp> get_all_point_expressions(string const &s, Data const* data,
bool only_active)
{
vector<fp> values;
vector<Point> const& points = data->points();
clear_parse_vecs();
// First compile string... puts result into code etc.
parse_info<> result = parse(s.c_str(), DataExpressionG, space_p);
if (!result.full)
throw ExecuteError("Syntax error in expression: " + s);
int M = (int) points.size();
vector<Point> new_points = points;
vector<fp> stack(stack_size);
for (vector<ParameterizedFunction*>::iterator i = parameterized.begin();
i != parameterized.end(); ++i)
(*i)->prepare_parameters(points);
replace_aggregates(M, points, code, code.begin());
for (int i = 0; i < M; ++i) {
if (only_active && !points[i].is_active)
continue;
execute_code(i, M, stack, points, new_points, code);
values.push_back(stack.front());
}
return values;
}
syntax highlighted by Code2HTML, v. 0.9.1