amroc/amr/AMRStatistics.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) California Institute of Technology
// Carlos Pantano

#ifndef AMROC_AMRStatistics_H
#define AMROC_AMRStatistics_H

#include <iostream>
#include <fstream>
#include <string>
#include <stdlib.h>
#include <algorithm>
#include "DAGH.h"
#include "AMRBase.h"
#include "Interpolation.h"

#include "StatParser/StatParser.h"

#define DEBUG_PRINT_STATISTICS

#define MAXGROUPS 8
#define MAXSTATOPS 4
#define MAXPERVAR 4
#define AVGOP 0
#define MODEOP 1
#define MEDIANOP 2
#define STDDEVIATIONOP 3
#define MINMAX 4
#define GRADOP 5
#define CURLOP 6
#define FFTOP 7

#define NUMBINS 10
#define TOLSTAT 1e-8
#define MINFIN 0.9e+37

#define MAXGRADS 4

template <class VectorType, class InterpolationType, class OutputType, int dim>
class AMRStatistics : public AMRBase<VectorType,dim>,
    private StatParser<typename InterpolationType::point_type,
    typename VectorType::InternalDataType,dim> {
  typedef AMRBase<VectorType,dim> base;
  typedef typename VectorType::InternalDataType DataType;
  typedef typename InterpolationType::point_type point_type;
  typedef StatParser<point_type,DataType,dim> stparser;
  typedef typename stparser::symrec symrec;

  // bring members from templated base class into scope
  using StatParser<point_type,DataType,dim>::group_table;
  using StatParser<point_type,DataType,dim>::putsym;
  using StatParser<point_type,DataType,dim>::unlinksym;
  using StatParser<point_type,DataType,dim>::vectorize;
  using StatParser<point_type,DataType,dim>::buildcoords;
  using StatParser<point_type,DataType,dim>::clearcoords;
  //using StatParser<point_type,DataType,dim>::clonestack;

public:
  typedef typename base::vec_grid_fct_type vec_grid_fct_type;
  typedef GridFunction<DataType,dim> grid_fct_type;

public:
  AMRStatistics(InterpolationType *int_ref, OutputType* out_ref) : base() {
    _Step = 1;
    _FileOutput = out_ref;
    _Interpolation = int_ref;
    _OutputFile = "stats";
    _DumpMeshes = 0;
    _NumStats = 0;
    //_NumPerVar = 0;
    _NumGrads = 0;
    //gcount = 0;
    _SepFiles = 0;


    for (register int i = 0; i < 3+MAXPERVAR; i++)
      for (register int j = 0; j < MAXSTATOPS; j++)
        statOp[i][j] = -1;

    std::vector<double> d;
    std::vector<std::vector<double> > axis;
    axis.push_back(d);
    ctmp.resize(2*dim+1);
    ftmp.resize(2*dim+1);
    /*for (register int i = 0; i < 1+2*dim; i++) {
      ftmp.push_back(d);

    }*/
    for (register int i = 0; i < MAXGRADS; i++) {
      calcount[i]=1;
      gcount[i]=0;
      curltog[i]=0;
      wvtog[i]=0;
      uwtog[i]=0;
      vutog[i]=0;
    }
    /*for (register int i = 0; i < dim; i++) {
      dircheck[i]=0;
      ctmp[i].resize(2*dim+1);
      for (register int j = 0; j < 2*dim+1; j++) {
        ctmp[i][j].resize(1);
      }
    }*/

    OpNames.clear();
    OpNames.push_back(" avgerage ");
    OpNames.push_back(" mode ");
    OpNames.push_back(" median ");
    OpNames.push_back(" stddeviation ");
    OpNames.push_back(" min-max ");

  }

  void Setup(GridHierarchy *gh, const int& ghosts, int* shape, double* geom) {
    base::SetupData(gh, ghosts);

    int me = MY_PROC;
    if (me == VizServer) {
      for (register int j = 0; j < _NumStats; j++) {
        for (register int i = 0; i < 3+MAXPERVAR; i++)
          std::cout << "["<<i<<"]["<<j<<"] " << statOp[i][j]<< " : ";
        std::cout << " _NumPerVar " << _NumPerVar[j] << std::endl;
      }
      for (register int i = 0; i < _NumGrads; i++) {
        std:: cout << "calcount["<<i<<"] " << calcount[i] << "\n";
        std::cout << "Gradgroup["<<i<<"] " << GradGroup[i] << "\n";
        for (register int d = 0; d < dim; d++) {
          std::cout << "Grads["<<d<<"]["<<i<<"] " << Grads[d][i] << "\n";
        }
      }
      for (register int i = 0; i < MAXGROUPS; i++) {
        std:: cout << "GroupVars["<<i<<"] " << GroupVars[i] << "\n";
      }


    }

    int MaxLev;
    int nx, ny, nz;
    char s[36];

    MaxLev = MaxLevel(base::GH());

    nx = shape[0];
    putsym(TOK_INTEGER, STAT_NX)->vdata.ivalue = nx;
    // add level dependent point number
    sprintf(s, "%s(%d)", STAT_NX, 1);
    putsym(TOK_INTEGER, s)->vdata.ivalue = nx;
    for ( int i = 0; i < MaxLev ; i++ ) {
      sprintf(s, "%s(%d)", STAT_NX, i+2);
      nx *= RefineFactor(base::GH(),i);
      putsym(TOK_INTEGER, s)->vdata.ivalue = nx;
    }
    putsym(TOK_INTEGER, STAT_NXF)->vdata.ivalue = nx;
    putsym(TOK_DOUBLE, STAT_XMIN)->vdata.dvalue = geom[0];
    putsym(TOK_DOUBLE, STAT_XMAX)->vdata.dvalue = geom[1];

    if ( dim > 1 ) {
      ny = shape[1];
      putsym(TOK_INTEGER, STAT_NY)->vdata.ivalue = ny;
      // add level dependent point number
      sprintf(s, "%s(%d)", STAT_NY, 1);
      putsym(TOK_INTEGER, s)->vdata.ivalue = ny;
      for ( int i = 0; i < MaxLev ; i++ ) {
        sprintf(s, "%s(%d)", STAT_NY, i+2);
        ny *= RefineFactor(base::GH(),i);;
        putsym(TOK_INTEGER, s)->vdata.ivalue = ny;
      }
      putsym(TOK_INTEGER, STAT_NYF)->vdata.ivalue = ny;
      putsym(TOK_DOUBLE, STAT_YMIN)->vdata.dvalue = geom[2];
      putsym(TOK_DOUBLE, STAT_YMAX)->vdata.dvalue = geom[3];
    }

    if ( dim > 2 ) {
      nz = shape[2];
      putsym(TOK_INTEGER, STAT_NZ)->vdata.ivalue = nz;
      // add level dependent point number
      sprintf(s, "%s(%d)", STAT_NZ, 1);
      putsym(TOK_INTEGER, s)->vdata.ivalue = nz;
      for ( int i = 0; i < MaxLev ; i++ ) {
        sprintf(s, "%s(%d)", STAT_NZ, i+2);
        nz *= RefineFactor(base::GH(),i);;
        putsym(TOK_INTEGER, s)->vdata.ivalue = nz;
      }
      putsym(TOK_INTEGER, STAT_NZF)->vdata.ivalue = nz;
      putsym(TOK_DOUBLE, STAT_ZMIN)->vdata.dvalue = geom[4];
      putsym(TOK_DOUBLE, STAT_ZMAX)->vdata.dvalue = geom[5];
    }

    /* simulation variables (dimension dependent) */
    // by convention I take positive iptr's as conservative
    // variables and negative iptr's as primitive variables.
    symrec *ptr;

    putsym(TOK_VARIABLE, "rho")->vdata.iptr = 1;
    putsym(TOK_VARIABLE, "r")->vdata.iptr = 1;
    (ptr = putsym(TOK_VARIABLE, "u"))->vdata.iptr = 2; ptr->src = STAT_DERIVED;
    if ( dim > 1 ) {
      (ptr = putsym(TOK_VARIABLE, "v"))->vdata.iptr = 3; ptr->src = STAT_DERIVED;
    }
    if ( dim > 2 ) {
      (ptr = putsym(TOK_VARIABLE, "w"))->vdata.iptr = 4; ptr->src = STAT_DERIVED;
    }
    putsym(TOK_VARIABLE, "E")->vdata.iptr = 5;
    for ( int i = 0; i < NVARS - 5 ; i++ ) {
      char name[12];
      sprintf(name, "Y%d", i+1);
      (ptr = putsym(TOK_VARIABLE, name))->vdata.iptr =
          ((dim == 1 ? 7 : (dim == 2 ? 8 : 9))+i); ptr->src = STAT_DERIVED;
    }
    if ( dim == 3 ) {
      (ptr = putsym(TOK_VARIABLE, "sgske"))->vdata.iptr = (NVARS+5); ptr->src = STAT_DERIVED;
    }
    putsym(TOK_VARIABLE, "T")->vdata.iptr = NVARS+1;
    (ptr = putsym(TOK_VARIABLE, "p"))->vdata.iptr =
        (dim == 1 ? 5 : (dim == 2 ? 6 : 7)); ptr->src = STAT_DERIVED;
    (ptr = putsym(TOK_VARIABLE, "gamma"))->vdata.iptr =
        (dim == 1 ? 6 : (dim == 2 ? 7 : 8)); ptr->src = STAT_DERIVED;
    putsym(TOK_VARINARRAY, "q")->vdata.iptr = 0;
    (ptr = putsym(TOK_VARINARRAY, "qout"))->vdata.iptr = 0; ptr->src = STAT_DERIVED;

    if ( _DefFile.length() ) {
      stparser::parse(_DefFile.c_str());
      for ( typename stparser::grpvec::const_iterator g = group_table.begin() ;
            g != group_table.end(); g++)
        if ( (*g)->step == 0 ) (*g)->step = _Step;

      int me = MY_PROC;
      if (me == VizServer) {
        if ( _DumpMeshes == 1) {
          int ig = 1;
          for ( typename stparser::grpvec::const_iterator g = group_table.begin() ;
                g != group_table.end(); g++, ig++) {
            int ip = 1;
            for ( typename stparser::provec::const_iterator p = (*g)->probe_table.begin() ;
                  p != (*g)->probe_table.end() ; p ++, ip++ )
            {
              if ( buildcoords(*p) != RETURN_OK ) return;

              int i, j;
              FILE * fp;
              char name[128];
#if defined(WJSON) || defined(SJSON)
              sprintf(name, "group%d-probe%d.json", ig, ip);
#else
              sprintf(name, "group%d-probe%d.dat", ig, ip);
#endif
              fp = fopen(name, "w");
              if ( (*p)->type == TOK_THREAD ) {
                for ( i = 0 ; i < (*p)->options.thread.na ; i++ ) {
                  fprintf(fp, "%12.8lf %12.8lf %12.8lf\n", (*p)->coords[i](0),
                          (*p)->coords[i](1),
                          (*p)->coords[i](2));
                }
              } else {
                int na, nb;
                if ( (*p)->type == TOK_SURFACE ) {
                  na = (*p)->options.surface.na;
                  nb = (*p)->options.surface.nb;
                } else if ((*p)->type == TOK_PLANES ) {
                  na = (*p)->options.plane.na;
                  nb = (*p)->options.plane.nb;
                }
#if defined(WJSON) || defined(SJSON)
                fprintf(fp, "\{\n\t\"COORDS\": [\"X\", \"Y\", \"Z\"],\n");
                fprintf(fp, "\t\"NA\": %d, \"NB\": %d,\n", na, nb);
                fprintf(fp, "\t\"VALUES\": [\n");
#else
                fprintf(fp, "VARIABLES = \"X\", \"Y\", \"Z\"\n");
                fprintf(fp, "ZONE I=%d, J=%d, F=POINT\n", na, nb);
#endif
                for ( j = 0 ; j < nb ; j++ )
                  for ( i = 0 ; i < na ; i++ ) {
                    int l = i+j*na;
#if defined(WJSON) || defined(SJSON)
                    if (l < na*nb-1) fprintf(fp, "\t\t [ %12.8lf, %12.8lf, %12.8lf ],\n", (*p)->coords[l](0),
                                             (*p)->coords[l](1), (*p)->coords[l](2));
                        else fprintf(fp, "\t\t [ %12.8lf, %12.8lf, %12.8lf ]\n", (*p)->coords[l](0),
                                     (*p)->coords[l](1), (*p)->coords[l](2));
    #else
                    fprintf(fp, "%12.8lf %12.8lf %12.8lf\n", (*p)->coords[l](0),
                            (*p)->coords[l](1), (*p)->coords[l](2));
    #endif
                  }
#if defined(WJSON) || defined(SJSON)
                fprintf(fp, "\t]\n}");
#endif
              }
              fclose(fp);

              clearcoords(*p);
            }
          }
        }
      }
    }
  }

  virtual void register_at(ControlDevice& Ctrl) { register_at(Ctrl, ""); }
  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    LocCtrl = Ctrl.getSubDevice(prefix+"Statistics");

    RegisterAt(LocCtrl,"DefinitionFile",_DefFile);
    RegisterAt(LocCtrl,"OutputFile",_OutputFile);
    RegisterAt(LocCtrl,"Step",_Step);
    RegisterAt(LocCtrl,"WriteMeshes", _DumpMeshes);
    RegisterAt(LocCtrl,"SeparateFiles", _SepFiles);

    RegisterAt(LocCtrl,"NumStats", _NumStats);
    char VariableName[32];
    for (register int i = 0; i < MAXSTATOPS; i++) {
      std::sprintf(VariableName,"G(%02d)",i+1);
      RegisterAt(LocCtrl,std::string(VariableName),statOp[0][i]);
      std::sprintf(VariableName,"P(%02d)",i+1);
      RegisterAt(LocCtrl,std::string(VariableName),statOp[1][i]);
      std::sprintf(VariableName,"V(%02d)",i+1);
      RegisterAt(LocCtrl,std::string(VariableName),statOp[2][i]);
      std::sprintf(VariableName,"N(%02d)",i+1);
      RegisterAt(LocCtrl,std::string(VariableName),_NumPerVar[i]);
      for (register int j = 0; j < MAXPERVAR; j++) {
        std::sprintf(VariableName,"Op%02d(%02d)",j+1,i+1);
        RegisterAt(LocCtrl,std::string(VariableName),statOp[3+j][i]);
      }
    }

    RegisterAt(LocCtrl,"NumGrads", _NumGrads);
    for (register int j = 0; j < MAXGRADS; j++) {
      std::sprintf(VariableName,"GradGroup(%02d)",j+1);
      RegisterAt(LocCtrl,std::string(VariableName),GradGroup[j]);
      std::sprintf(VariableName,"GradVars(%02d)",j+1);
      RegisterAt(LocCtrl,std::string(VariableName),GradVars[j]);
      std::sprintf(VariableName,"Curl(%02d)",j+1);
      RegisterAt(LocCtrl,std::string(VariableName),curltog[j]);
      std::sprintf(VariableName,"wv(%02d)",j+1);
      RegisterAt(LocCtrl,std::string(VariableName),wvtog[j]);
      std::sprintf(VariableName,"uw(%02d)",j+1);
      RegisterAt(LocCtrl,std::string(VariableName),uwtog[j]);
      std::sprintf(VariableName,"vu(%02d)",j+1);
      RegisterAt(LocCtrl,std::string(VariableName),vutog[j]);
    }
    for (register int i = 0; i < dim; i++) {
      for (register int j = 0; j < MAXGRADS; j++) {
        if (i == 0)
          std::sprintf(VariableName,"GradX(%02d)",j+1);
        else if (i == 1)
          std::sprintf(VariableName,"GradY(%02d)",j+1);
        else
          std::sprintf(VariableName,"GradZ(%02d)",j+1);
        RegisterAt(LocCtrl,std::string(VariableName),Grads[i][j]);
      }
    }

    for (register int j = 0; j < MAXGROUPS; j++) {
      std::sprintf(VariableName,"GroupVars(%02d)",j+1);
      RegisterAt(LocCtrl,std::string(VariableName),GroupVars[j]);
    }

  }

  ~AMRStatistics() {
    _FileOutput = NULL;
    _Interpolation = NULL;
  }

public:

  virtual void CalcStat(int OpType, int Lpoints, double * adata, double * result) {

    double avg = 0., bin = 0., stddev = 0., st1 = 0., st2 = 0., min = MINFIN, max = -MINFIN;
    int intervals[NUMBINS+1], mfid = 0, mfreq = 0, count = 0, mini = 0, maxi = 0;
    std::vector<double> tmp( adata, adata + Lpoints );//*sizeof(double) );
    switch ( OpType ) {
    case AVGOP:  // 0
      for (register int i=0; i<Lpoints; i++) {
        avg += tmp[i];
      }
      avg /= Lpoints;
      *result = avg;
      break;
    case MODEOP: // 1
      tmp.resize(Lpoints);
      sort(tmp.begin(),tmp.end());
      for (register int i = 0; i < NUMBINS+1; i++ ) {
        intervals[i] = 0;
      }
      bin = (tmp[Lpoints-1] - tmp[0])/(double)NUMBINS;
      if (bin > TOLSTAT) {
        for (register int i = 0; i < Lpoints; i++ ) {
          intervals[ (int) std::floor( (tmp[i]-tmp[0])/bin )  ] += 1;
        }
        for (register int b = 0; b < NUMBINS+1; b++ ) {
          count = intervals[b];
          if (count > mfreq) {
            mfreq = count;
            mfid = b;
          }
        }
        *result = tmp[mfid];
      }
      else {
        *result = tmp[Lpoints/2];
      }
      tmp.clear();
      break;
    case MEDIANOP: // 2
      tmp.resize(Lpoints);
      sort(tmp.begin(),tmp.end());
      *result = tmp[Lpoints/2];
      tmp.clear();
      break;
    case STDDEVIATIONOP: // 3
      avg = 0.;
      st1 = 0.;
      st2 = 0.;
      for (register int i=0; i<Lpoints; i++) {
        stddev += (adata[i]-avg)*(adata[i]-avg);
        st1 += (adata[i]);
        st2 += (adata[i])*(adata[i]);
      }
      st1 = std::sqrt( ((double)Lpoints*st2 - st1*st1)/(double)(Lpoints*(Lpoints-1)) );
      *result = st1;
      break;
    case MINMAX:  // 4
      for (register int i=0; i<Lpoints; i++) {
        if (min > tmp[i]) { min = tmp[i]; mini = i; }
        if (max < tmp[i]) { max = tmp[i]; maxi = i; }
      }
      result[0] = min;  result[1] = (double)mini;
      result[2] = max;  result[3] = (double)maxi;
      break;

    }

  }

#ifdef SJSON
  virtual void CalcGradient(int groupNum, int * dir, int nvars, int Lpoints, int na, int nb, std::vector<std::vector<point_type> > cdata,
                            std::vector<std::vector<double> > vdata, int nstep, double * t, int curl, int wv, int uw, int vu) {
#else
  virtual void CalcGradient(int groupNum, int * dir, int nvars, int Lpoints, int na, int nb, std::vector<std::vector<point_type> > cdata,
                            std::vector<std::vector<double> > vdata, int nstep, double * t, std::ostream* out, int curl, int wv, int uw, int vu) {
#endif

    std::vector<std::vector<double> > gradi;
    std::vector<double> res;
    res.resize(Lpoints);
    int pos,posi;

#ifdef STATS_DEBUG
    std::cout << "******* *CalgGrad()* ********\n";
    std::cout << "vdata.size() " << vdata.size() << std::endl;
    std::cout << "vdata[0].size() " << vdata[0].size() << std::endl;
#endif

#if defined(WJSON)
    *out << "\t\t\t\"Gradient\": {\n";
#elif !defined(WJSON) && !defined(SJSON)
    *out << "# Gradient\n";
#endif
    for (register int nv=0; nv<nvars; nv++) {
#if defined(WJSON)
      *out << "\t\t\t\t\"VAR" << nv+1 << "\": {\n";
#elif !defined(WJSON) && !defined(SJSON)
      *out << "# NVAR " << nv+1 << "\n";
#endif
      pos = nv*Lpoints;
      for ( register int d=0; d<dim; d++) {
        if (dir[d] == 1) {
#if defined(WJSON)
          *out << "\t\t\t\t\t\"dim" << d << "\": [";
#elif !defined(WJSON) && !defined(SJSON)
          *out << "# dim " << d << "\n";
#endif
          for (register int i=0; i<Lpoints; i++) {
            posi = pos+i;
            if (vdata[0][posi] < -MINFIN ) {
              res[i] = MINFIN;
            } else if (vdata[2*d+1][posi] < -MINFIN ) {
              res[i] = (vdata[0][posi] -vdata[2*d+2][posi])/(cdata[0][i](d)-cdata[2*d+2][i](d));
            } else if (vdata[2*d+2][posi] < -MINFIN ) {
              res[i] = (vdata[2*d+1][posi] -vdata[0][posi])/(cdata[2*d+1][i](d)-cdata[0][i](d));
            } else if (vdata[2*d][posi] > -MINFIN && vdata[2*d+1][posi] > -MINFIN &&
                       vdata[2*d+2][posi] > -MINFIN) {
              res[i] = (vdata[2*d+1][posi] -vdata[2*d+2][posi])/(cdata[2*d+1][i](d)-cdata[2*d+2][i](d));
            }
#if defined(WJSON)
            if (i<Lpoints-1) *out << res[i] << ", ";
            else {
              if (d<dim-1) *out << res[i] << " ],";
              else *out << res[i] << " ]";
            }
#elif !defined(WJSON) && !defined(SJSON)
            *out << res[i] << " ";
#endif
          }
#ifndef SJSON
          *out << std::endl;
#endif

          gradi.push_back(res);

        }
      }
#if defined(WJSON)
      if (nv<nvars-1) *out << "\t\t\t\t},\n";
      else *out << "\t\t\t\t}\n";
#endif
    }
#if defined(WJSON)
    if (curl) {
      *out << "\t\t\t},\n";
    } else *out << "\t\t\t}\n";
#endif

    if (curl > 0) {
#if defined(WJSON)
      *out << "\t\t\t\"CURL\": { \n";
#elif defined SJSON
      if (_SepFiles == 1) {
        char sfname[256];
        std::ofstream sepOut;
        if (wv) {
          std::sprintf(sfname,"%s%d_wv_%d.json",_OutputFile.c_str(),groupNum,CurrentTime(base::GH(), 0) );
          std::cout << " *** Writing " << sfname << std::endl;
          sepOut.open(sfname, std::ios::out );
          sepOut << "{\n\t\"iter\": " << CurrentTime(base::GH(), 0) << ", \"time\": " << t[0] << "," << std::endl;
          sepOut << "\t\"wv\": [ ";
          for (register int i=0; i<Lpoints-1; i++) {
            if (i<Lpoints-1) sepOut << gradi[3][i] - gradi[2][i] << ", ";
          }
          sepOut << gradi[3][Lpoints-1] - gradi[2][Lpoints-1] << " ]\n}";
          sepOut.close();
        }
        if (uw) {
          std::sprintf(sfname,"%s%d_uw_%d.json",_OutputFile.c_str(),groupNum,CurrentTime(base::GH(), 0) );
          std::cout << " *** Writing "  << sfname << std::endl;
          sepOut.open(sfname, std::ios::out );
          sepOut << "{\n\t\"iter\": " << CurrentTime(base::GH(), 0) << ", \"time\": " << t[0] << "," << std::endl;
          sepOut << "\t\"uw\": [ ";
          for (register int i=0; i<Lpoints-1; i++) {
            if (i<Lpoints-1) sepOut << gradi[1][i] - gradi[3][i] << ", ";
          }
          sepOut << gradi[1][Lpoints-1] - gradi[3][Lpoints-1] << " ]\n}";
          sepOut.close();
        }
        if (vu) {
          std::sprintf(sfname,"%s%d_vu_%d.json",_OutputFile.c_str(),groupNum,CurrentTime(base::GH(), 0) );
          std::cout << " *** Writing "  << sfname << std::endl;
          sepOut.open(sfname, std::ios::out );
          sepOut << "{\n\t\"iter\": " << CurrentTime(base::GH(), 0) << ", \"time\": " << t[0] << "," << std::endl;
          sepOut << "\t\"vu\": [ ";
          for (register int i=0; i<Lpoints-1; i++) {
            if (i<Lpoints-1) sepOut << gradi[2][i] - gradi[1][i] << ", ";
          }
          sepOut << gradi[2][Lpoints-1] - gradi[1][Lpoints-1] << " ]\n}";
          sepOut.close();
        }
      }
#else
      *out << "# CURL\n";
#endif

#if defined(WJSON)
      *out << "\t\t\t\t\"wv\": [ ";
#endif
      for (register int i=0; i<Lpoints; i++) {
#if defined(WJSON)
        if (i<Lpoints-1) *out << gradi[3][i] - gradi[2][i] << ", ";
        else *out << gradi[3][i] - gradi[2][i] << " ],";
#elif !defined(SJSON)
        *out << gradi[3][i] - gradi[2][i] << " ";
#endif
      }
#if defined(WJSON)
      *out << std::endl;
      *out << "\t\t\t\t\"uw\": [ ";
#endif
      for (register int i=0; i<Lpoints; i++) {
#if defined(WJSON)
        if (i<Lpoints-1) *out << gradi[1][i] - gradi[3][i] << ", ";
        else *out << gradi[1][i] - gradi[3][i] << "], ";
#elif !defined(SJSON)
        *out << gradi[1][i] - gradi[3][i] << " ";
#endif
      }
#if defined(WJSON)
      *out << std::endl;
      *out << "\t\t\t\t\"vu\": [ ";
#endif
      for (register int i=0; i<Lpoints; i++) {
#if defined(WJSON)
        if (i<Lpoints-1) *out  << gradi[2][i] - gradi[1][i] << ", ";
        else *out  << gradi[2][i] - gradi[1][i] << "]";
#elif !defined(SJSON)
        *out  << gradi[2][i] - gradi[1][i] << " ";
#endif
      }
#if defined(WJSON)
      *out << std::endl;
#endif
    }
#if defined(WJSON)
    *out << "\t\t\t},\n";
#endif
  }

  void Evaluate(double * t, vec_grid_fct_type& U, grid_fct_type& Work)
  {
    int me;
    int Npoints;

    char fname[256];
    std::ofstream outfile;
    std::ostream* out;

    if ( group_table.size() == 0 ) return;

    int TimeC = StepsTaken(base::GH(),0);

    // check step sampling rate
    int gcheck = 0;
    for ( typename stparser::grpvec::const_iterator g = group_table.begin() ;
          g != group_table.end(); g++) {
      if ( TimeC%((*g)->step) == 0 ) {
        gcheck = 1;
        break;
      }
    }

    if ( gcheck == 0 ) return;

    me = MY_PROC;

    if (me == VizServer) {
#if defined(WJSON)
      std::sprintf(fname,"%s%d.json",_OutputFile.c_str(),CurrentTime(base::GH(), 0));
      outfile.open(fname, std::ios::out );
      out = new std::ostream(outfile.rdbuf());
#elif !defined(WJSON) && !defined(SJSON)
      std::sprintf(fname,"%s%d.dat",_OutputFile.c_str(),CurrentTime(base::GH(), 0));
      outfile.open(fname, std::ios::out );
      out = new std::ostream(outfile.rdbuf());
#endif

#if defined(WJSON)
      *out << "{\n\t\"iter\": " << CurrentTime(base::GH(), 0) << ", \"time\": " << t[0] << "," << std::endl; //newformat
#elif !defined(WJSON) && !defined(SJSON)
      *out << "# time " << t[0] << std::endl;
#endif
    }

    int ig = 1, iptr;
    for ( typename stparser::grpvec::const_iterator g = group_table.begin() ;
          g != group_table.end(); g++, ig++ ) {

      int ipv = 0;
      if ( TimeC%((*g)->step) == 0 ) {
        if (me == VizServer) {
#if defined(WJSON)
          *out << "\t\t\"group" << ig << "\": {"<< std::endl;
#elif !defined(WJSON) && !defined(SJSON)
          *out << "# group " << ig << std::endl;
#endif

        }

        Npoints = 0;
        for ( typename stparser::provec::const_iterator p = (*g)->probe_table.begin() ;
              p != (*g)->probe_table.end() ; p ++) {
          ipv = 0;
          switch ( (*p)->type ) {
          case TOK_THREAD:
            Npoints += (*p)->options.thread.na;
            break;
          case TOK_SURFACE:
            Npoints += (*p)->options.surface.na*(*p)->options.surface.nb;
            break;
          case TOK_PLANES:
            Npoints += (*p)->options.plane.na*(*p)->options.plane.nb;
            break;
          case TOK_VOLUME:
            /* not implemented yet */
            break;
          }
        }

        point_type* CList = new point_type[Npoints];
        int cli=0;

        for ( typename stparser::provec::const_iterator p = (*g)->probe_table.begin() ;
              p != (*g)->probe_table.end() ; p ++) {
          if ( buildcoords(*p) != RETURN_OK ) return;
          int Lpoints;
          switch ( (*p)->type ) {
          case TOK_THREAD:
            Lpoints = (*p)->options.thread.na;
            break;
          case TOK_SURFACE:
            Lpoints = (*p)->options.surface.na*(*p)->options.surface.nb;
            break;
          case TOK_PLANES:
            Lpoints = (*p)->options.plane.na*(*p)->options.plane.nb;
            break;
          case TOK_VOLUME:
            /* not implemented yet */
            break;
          }
          for (register int li=0; li<Lpoints; li++, cli++)
            CList[cli] = (*p)->coords[li];
          if (me == VizServer)
            (*p)->out = new std::stringstream(std::stringstream::in | std::stringstream::out);
        }

        bool *lused = new bool[Npoints];
        int Level = 0;
        int Time = CurrentTime(base::GH(),Level);
        int NpUsed =_Interpolation->LocalPoints(Work,Time,Level,Npoints,CList,lused);
        point_type* CUList = new point_type[NpUsed];
        cli = 0;
        for (register int i=0; i<Npoints; i++)
          if (lused[i]) CUList[cli++] = CList[i];
        delete [] CList;

        symrec *ptr, *ptr_ref;

        /* first, make a copy of the stack */
        ptr = stparser::clonestack((*g)->sym_table);

        /* traverse the stack and evaluate */
        while ( ptr ) {
#ifdef DEBUG_PRINT_STATISTICS
          (comm_service::log() << stparser::token(ptr->type) << "\n").flush();
#endif
          // vectorize the data array that will be used for the result
          vectorize(ptr, NpUsed);
          // the arguments of the operator do not need to be vectorized
          // for all processors since they are not used. Only VizServer
          // performs operations with them.

          switch( ptr->type ) {
          case TOK_OP_PLUS:
            vectorize(ptr->prev, NpUsed);
            vectorize(ptr->prev->prev, NpUsed);
            for ( register int l = 0 ; l < NpUsed ; l++ )
              ptr->pdata[l] = ptr->prev->prev->pdata[l] + ptr->prev->pdata[l];
            unlinksym(ptr->prev->prev);
            unlinksym(ptr->prev);
            break;
          case TOK_OP_MINUS:
            vectorize(ptr->prev, NpUsed);
            vectorize(ptr->prev->prev, NpUsed);
            for ( register int l = 0 ; l < NpUsed ; l++ )
              ptr->pdata[l] = ptr->prev->prev->pdata[l] - ptr->prev->pdata[l];
            unlinksym(ptr->prev->prev);
            unlinksym(ptr->prev);
            break;
          case TOK_OP_PROD:
            vectorize(ptr->prev, NpUsed);
            vectorize(ptr->prev->prev, NpUsed);
            for ( register int l = 0 ; l < NpUsed ; l++ )
              ptr->pdata[l] = ptr->prev->prev->pdata[l] * ptr->prev->pdata[l];
            unlinksym(ptr->prev->prev);
            unlinksym(ptr->prev);
            break;
          case TOK_OP_DIV:
            vectorize(ptr->prev, NpUsed);
            vectorize(ptr->prev->prev, NpUsed);
            for ( register int l = 0 ; l < NpUsed ; l++ )
              ptr->pdata[l] = ptr->prev->prev->pdata[l] / ptr->prev->pdata[l];
            unlinksym(ptr->prev->prev);
            unlinksym(ptr->prev);
            break;
          case TOK_OP_NEG:
            vectorize(ptr->prev, NpUsed);
            for ( register int l = 0 ; l < NpUsed ; l++ )
              ptr->pdata[l] = - ptr->prev->pdata[l];
            unlinksym(ptr->prev);
            break;
          case TOK_OP_POWER:
            vectorize(ptr->prev, NpUsed);
            vectorize(ptr->prev->prev, NpUsed);
            for ( register int l = 0 ; l < NpUsed ; l++ )
              ptr->pdata[l] = pow(ptr->prev->prev->pdata[l], ptr->prev->pdata[l]);
            unlinksym(ptr->prev->prev);
            unlinksym(ptr->prev);
            break;
          case TOK_POINTER:
            ptr_ref = ptr->ref;
            switch(ptr_ref->type) {
            case TOK_FUNCTION1:
              vectorize(ptr->prev, NpUsed);
              for ( register int l = 0 ; l < NpUsed ; l++ )
                ptr->pdata[l] = ((func_t1)ptr_ref->vdata.ptr)(ptr->prev->pdata[l]);
              unlinksym(ptr->prev);
              break;
            case TOK_FUNCTION2:
              vectorize(ptr->prev, NpUsed);
              vectorize(ptr->prev->prev, NpUsed);
              for ( register int l = 0 ; l < NpUsed ; l++ )
                ptr->pdata[l] = ((func_t2)ptr_ref->vdata.ptr)(ptr->prev->prev->pdata[l],
                                                              ptr->prev->pdata[l]);
              unlinksym(ptr->prev->prev);
              unlinksym(ptr->prev);
              break;
            case TOK_VARINARRAY:
            case TOK_VARIABLE:

              // for variables in one of the AMR arrays, the pointer is located in the
              // current token, not the defered one (ref). For variables in regular
              // variables, the pointer is stored in the pointed variable.
              if ( ptr_ref->type == TOK_VARINARRAY )
                iptr = ptr->vdata.iptr;
              else
                iptr = ptr_ref->vdata.iptr;

              // process the variables here
              if ( ptr_ref->src == STAT_DERIVED ) {
                // output variable type case
                for (register int l=0; l<=FineLevel(base::GH()); l++) {
                  int Time = CurrentTime(base::GH(),l);
                  // if derived variable (call Transform)
                  _FileOutput->Transform(U, Work, Time, l, iptr, t[l]);
                }
#ifdef DEBUG_PRINT_STATISTICS
                (comm_service::log() << " convert output\n" ).flush();
#endif
              } else {
                int _dim = dim;
                // conservative vector of state variable case
                if ( _dim == 2 ) {
                  for ( int Level = FineLevel(base::GH()) ; Level >= 0 ; Level --) {
                    int Time = CurrentTime(base::GH(),Level);
                    forall (U,Time,Level,c)
                        Coords ss(U(Time,Level,c).bbox().stepsize());
                    BBox bbi(U.boundingbbox(Time,Level,c));
                    BeginFastIndex2(U, U(Time,Level,c).bbox(),
                                    U(Time,Level,c).data(), VectorType);
                    BeginFastIndex2(Work, Work(Time,Level,c).bbox(),
                                    Work(Time,Level,c).data(), DataType);
                    for_2 (i, j, bbi, ss)
                        FastIndex2(Work,i,j) = FastIndex2(U,i,j)(iptr-1);
                    end_for
                        EndFastIndex2(U);
                    EndFastIndex2(Work);
                    end_forall
                  }
                } else if ( _dim == 3 ) {
                  for ( int Level = FineLevel(base::GH()) ; Level >= 0 ; Level --) {
                    int Time = CurrentTime(base::GH(),Level);
                    forall (U,Time,Level,c)
                        Coords ss(U(Time,Level,c).bbox().stepsize());
                    BBox bbi(U.boundingbbox(Time,Level,c));
                    BeginFastIndex3(U, U(Time,Level,c).bbox(),
                                    U(Time,Level,c).data(), VectorType);
                    BeginFastIndex3(Work, Work(Time,Level,c).bbox(),
                                    Work(Time,Level,c).data(), DataType);
                    for_3 (i, j, k, bbi, ss)
                        FastIndex3(Work,i,j,k) = FastIndex3(U,i,j,k)(iptr-1);
                    end_for
                        EndFastIndex3(U);
                    EndFastIndex3(Work);
                    end_forall
                  }
                }
#ifdef DEBUG_PRINT_STATISTICS
                (comm_service::log() << " convert array\n" ).flush();
#endif
              }
              // note that we use data so the pointer token will be converted.
              _Interpolation->PointsValues(Work,t[0],NpUsed,CUList,ptr->pdata);
#ifdef DEBUG_PRINT_STATISTICS
              (comm_service::log() << "after\n" ).flush();
#endif
              break;
            case TOK_COORDINATE:
            {
              int iptr = (ptr_ref->vdata.iptr)-1;
              if ( iptr >= dim ) {
                staterror("invalid coordinate in this stack");
                return;
              }
              for ( register int l = 0 ; l < NpUsed ; l++ )
                ptr->pdata[l] = CUList[l](iptr);
            }
              break;
            default:
              staterror("there should not be any other pointer type left in this stack");
              return;
            }

            break;

          case TOK_SEPARATOR:
            // output
          {
            int na, nb, avea, aveb;
            vectorize(ptr->prev, NpUsed);

            cli=0;
            int ip = 1;
            int ccli=0;
            for ( typename stparser::provec::const_iterator p = (*g)->probe_table.begin() ;
                  p != (*g)->probe_table.end() ; p ++, ip++ ) {
              // write probe data
              int Lpoints;
              //ipv=0;
              switch((*p)->type) {
              case TOK_THREAD:
                Lpoints = (*p)->options.thread.na;
                if ( (*p)->options.thread.ave == STAT_AVERAGED ) {
                  double sum = 0.0, rsum;
                  for ( register int ipl = 0 ; ipl < Lpoints ; ipl ++ )
                    if (lused[ipl+cli])
                      sum += ptr->prev->pdata[ccli++];
                  MPI_Reduce(&sum, &rsum, 1, MPI_DOUBLE, MPI_SUM, VizServer, comm_service::comm());
                  if (me == VizServer)
                    *(*p)->out << (rsum/Lpoints) << std::endl;
                } else {
                  double *adata = new double[Lpoints];
                  for ( register int ipl = 0 ; ipl < Lpoints ; ipl ++ )
                    if (lused[ipl+cli])
                      adata[ipl] = ptr->prev->pdata[ccli++];
                    else
                      adata[ipl] = _Interpolation->ErrorValue();
                  _Interpolation->ArrayCombine(VizServer, Lpoints, adata);
                  if (me == VizServer) {
                    double *result = 0;// = new double[3];
                    int ti = 0;
                    for (register int i = 0; i < MAXSTATOPS; i++)
                      if ( statOp[0][i] == ig ) {
                        //std::cout << " DO a Stat on group " << ig << std::endl;
                        if ( statOp[1][i] == ip) {
                          //std::cout << " DO a Stat on probe " << ip << std::endl;
                          if ( statOp[2][i] == ipv) {
                            //std::cout << " DO a Stat on group " << ig; // << std::endl;
                            //std::cout << " on probe " << ip << std::endl;
                            //std::cout << "  __ operation " << statOp[3][i] << " Stat on variable " << ipv << std::endl;
                            ti = i;
                            if (statOp[3][i] < 4) {
                              result = new double;//[3];
                              CalcStat(statOp[3][i], Lpoints, adata, result);
                            }
                            else if (statOp[3][i] == 4) {
                              //double *result = new double[3];
                              //result = new double;//[3];
                              //CalcStat(statOp[3][i], Lpoints, adata, result);
                            }

                          }
                        }
                      }

                    for ( register int ipl = 0 ; ipl < Lpoints ; ipl ++ ) {

                      *(*p)->out << adata[ipl] << " ";
                    }
#if defined SJSON
                    //std::cout << " _SepFiles = " << _SepFiles << " _NumGrads " << _NumGrads << std::endl;
                    if (_SepFiles == 1) {
                      int write = 1;
                      int gtest = 1;
                      if (_NumGrads > 1 ) gtest = _NumGrads-1;
                      for (register int n = 0; n < gtest; n++) {
                        if (ig == GradGroup[n]) {
                          write = 0;
                        }
                        if (write == 1) {
                          char sfname[256];
                          std::ofstream sepOut;
                          std::sprintf(sfname,"%s%d_%d_%d_%d.json",_OutputFile.c_str(), ig, ip, ipv, CurrentTime(base::GH(), 0) );
                          std::cout << " *** Writing " << sfname << std::endl;
                          sepOut.open(sfname, std::ios::out );
                          sepOut << "{\n\t\"iter\": " << CurrentTime(base::GH(), 0) << ", \"time\": " << t[0] << "," << std::endl;
                          sepOut << "\t\"VAR"<<ipv<<"\": [ ";
                          for ( register int ipl = 0 ; ipl < Lpoints-1 ; ipl ++ ) {
                            sepOut << adata[ipl] << ", ";
                            //std::cout << "adata["<<ipl<<"] "<< adata[ipl] << "\n";
                          }
                          sepOut << adata[Lpoints-1] << " ]\n}";
                          sepOut.close();
                        }
                      }
                    }
#endif

                    if (result) {
                      if (statOp[3][ti] < 4) {
                        *(*p)->out << *result << " ";
                        //std::cout << "result"<< *result << "\n";
                      }
                      else if (statOp[3][ti] == 4) {
                        *(*p)->out << *result << " ";
                        //std::cout << "result"<< *result << "\n";
                        //*(*p)->out << result[0] << " " << result[1] << " " << result[2] << " ";
                        //std::cout << "result "<< result[0] << "\t"<< result[1] << "\t"<< result[2] << "\n";
                      }
                    }
#ifndef SJSON
                    *(*p)->out << std::endl;
#endif
                    ipv++;
                    delete [] result;
                  }
                  delete [] adata;
                }
                break;
              case TOK_SURFACE:
                Lpoints = (*p)->options.surface.na*(*p)->options.surface.nb;
                avea = (*p)->options.surface.avea;
                aveb = (*p)->options.surface.aveb;
                na = (*p)->options.surface.na;
                nb = (*p)->options.surface.nb;
              case TOK_PLANES:
                if ( (*p)->type == TOK_PLANES ) {
                  Lpoints = (*p)->options.plane.na*(*p)->options.plane.nb;
                  avea = STAT_AVERAGED;
                  aveb = STAT_AVERAGED;
                  na = (*p)->options.plane.na;
                  nb = (*p)->options.plane.nb;
                }

                if ( avea == STAT_AVERAGED && aveb == STAT_AVERAGED ) {
                  double sum = 0.0, rsum;
                  for ( register int ip = 0 ; ip < Lpoints ; ip ++ )
                    if (lused[ip+cli])
                      sum += ptr->prev->pdata[ccli++];
                  MPI_Reduce(&sum, &rsum, 1, MPI_DOUBLE, MPI_SUM, VizServer, comm_service::comm());
                  if (me == VizServer)
                    *(*p)->out << (rsum/Lpoints) << std::endl;
                } else {
                  double *adata = new double[Lpoints];
                  for ( register int j = 0 ; j < nb ; j++ )
                    for ( register int i = 0 ; i < na ; i++ )
                      if (lused[i+j*na+cli])
                        adata[i+j*na] = ptr->prev->pdata[ccli++];
                      else
                        adata[i+j*na] = _Interpolation->ErrorValue();

                  if ( avea == STAT_RAW && aveb == STAT_AVERAGED ) {
                    for ( register int i = 0 ; i < na ; i++ ) {
                      double sum=0.0, rsum;
                      for ( register int j = 0 ; j < nb ; j++ )
                        if (adata[i+j*na]!=_Interpolation->ErrorValue())
                          sum += adata[i+j*na];
                      MPI_Reduce(&sum, &rsum, 1, MPI_DOUBLE, MPI_SUM, VizServer, comm_service::comm());
                      if (me == VizServer)
                        *(*p)->out << (rsum/nb) << " ";
                    }
#if !defined(WJSON) && !defined(SJSON)
                    if (me == VizServer)
                      *(*p)->out << std::endl;
#endif
                  } else if ( avea == STAT_AVERAGED && aveb == STAT_RAW ) {
                    for ( register int j = 0 ; j < nb ; j++ ) {
                      double sum=0.0, rsum;
                      for ( register int i = 0 ; i < na ; i++ )
                        if (adata[i+j*na]!=_Interpolation->ErrorValue())
                          sum += adata[i+j*na];
                      MPI_Reduce(&sum, &rsum, 1, MPI_DOUBLE, MPI_SUM, VizServer, comm_service::comm());
                      if (me == VizServer)
                        *(*p)->out << (rsum/na) << " ";
                    }
#if !defined(WJSON) && !defined(SJSON)
                    if (me == VizServer)
                      *(*p)->out << std::endl;
#endif
                  } else {
                    _Interpolation->ArrayCombine(VizServer, Lpoints, adata);
                    if (me == VizServer) {
                      //CHANGE VARipv wo actual variable name from key
#ifdef WJSON
                      //if (ipv>0) *(*p)->out << "\t\t\t\t,\"VAR"<<ipv<<"\": [ ";
                      //else *(*p)->out << "\t\t\t\t\"VAR"<<ipv<<"\": [ ";
                      *(*p)->out << "\t\t\t\t\"VAR"<<ipv<<"\": [ ";
                      for ( register int ipl = 0 ; ipl < Lpoints ; ipl ++ ) {
                        *(*p)->out << adata[ipl];
                        if (ipl<Lpoints-1) *(*p)->out << ", ";
                        //std::cout << "adata["<<ipl<<"] "<< adata[ipl] << "\n";
                      }
                      if (ipv<GroupVars[ig-1]-1 ) *(*p)->out << " ],";//<<ipv<<" "<<GroupVars[ig-1]-1;
                      else *(*p)->out << " ]";//<<ipv<<" "<<GroupVars[ig-1]-1;
#elif defined(SJSON)
                      //std::cout << " _SepFiles = " << _SepFiles << " _NumGrads " << _NumGrads << std::endl;
                      if (_SepFiles == 1) {
                        int write = 1;
                        int gtest = 1;
                        if (_NumGrads > 1 ) gtest = _NumGrads-1;
                        for (register int n = 0; n < gtest; n++) {
                          //for (register int n = 0; n <_NumGrads+1; n++) {
                          if (ig == GradGroup[n]) {
                            write = 0;
                          }
                          if (write == 1) {
                            char sfname[256];
                            std::ofstream sepOut;
                            std::sprintf(sfname,"%s%d_%d_%d_%d.json",_OutputFile.c_str(), ig, ip, ipv, CurrentTime(base::GH(), 0) );
                            std::cout << " *** Writing " << sfname << std::endl;
                            sepOut.open(sfname, std::ios::out );
                            sepOut << "{\n\t\"iter\": " << CurrentTime(base::GH(), 0) << ", \"time\": " << t[0] << "," << std::endl;
                            sepOut << "\t\"VAR"<<ipv<<"\": [ ";
                            for ( register int ipl = 0 ; ipl < Lpoints-1 ; ipl ++ ) {
                              sepOut << adata[ipl] << ", ";
                              //std::cout << "adata["<<ipl<<"] "<< adata[ipl] << "\n";
                            }
                            sepOut << adata[Lpoints-1] << " ]\n}";
                            sepOut.close();
                          }
                        }
                      }

#else
                      for ( register int ipl = 0 ; ipl < Lpoints ; ipl ++ ) {
                        *(*p)->out << adata[ipl] << " ";
                        //std::cout << "adata["<<ipl<<"] "<< adata[ipl] << "\n";
                      }
#endif
                      double * result = 0;// = new double[3];
                      std::vector<double>stats;

                      for (register int i = 0; i < _NumStats; i++) {
                        if ( statOp[0][i] == ig ) {
                          //std::cout << " DO a Stat on group " << ig << std::endl;
                          if ( statOp[1][i] == ip) {
                            //std::cout << " DO a Stat on probe " << ip << std::endl;
                            if ( statOp[2][i] == ipv) {
#ifdef STATS_DEBUG
                              std::cout << " DO a Stat on group " << ig; // << std::endl;
                              std::cout << " on probe " << ip << std::endl;
#endif
                              //std::cout << "  __ operation " << statOp[3][i] << " Stat on variable " << ipv << std::endl;
                              for (register int j = 0; j < _NumPerVar[i]; j++) {
                                if (statOp[3+j][i] >= 0) {
                                  if (statOp[3+j][i] < 4) {
                                    //std::cout << "  __ operation " << statOp[3+j][i] << " Stat on variable " << ipv << std::endl;
                                    result = new double;//[3];
                                    CalcStat(statOp[3+j][i], Lpoints, adata, result);
#if defined(WJSON) || defined(SJSON)
                                    if (j==0 && ipv==GroupVars[ig-1]-1 ) *(*p)->out << ",\n\t\t\t\t\"VAR"<< ipv <<"STAT" << j << "\": { \"type\": \"" << OpNames[statOp[3+j][i]] << "\", ";
                                    else *(*p)->out << "\n\t\t\t\t\"VAR"<< ipv <<"STAT" << j << "\": { \"type\": \"" << OpNames[statOp[3+j][i]] << "\", ";

                                    if (j==_NumPerVar[i]-1 && ipv==GroupVars[ig-1]-1 ) *(*p)->out << "\"value\": " << result[0] << " }";
                                    else *(*p)->out << "\"value\": " << result[0] << " },";
#else
                                    *(*p)->out << "\n# STAT " << statOp[3+j][i]  << " :" << OpNames[statOp[3+j][i]] << "\n";
                                    *(*p)->out << result[0] << " ";
#endif

                                  }
                                  else if (statOp[3+j][i] == 4) {
                                    //std::cout << "  __ operation " << statOp[3+j][i] << " Stat on variable " << ipv << std::endl;
                                    result = new double[4];
                                    CalcStat(statOp[3+j][i], Lpoints, adata, result);
#if defined(WJSON) || defined(SJSON)
                                    if (j==0 && ipv==GroupVars[ig-1]-1 ) *(*p)->out << ",\n\t\t\t\t\"VAR"<< ipv <<"STAT" << j << "\": { \"type\": \""<< OpNames[statOp[3+j][i]] << "\",";
                                    else *(*p)->out << "\n\t\t\t\t\"VAR"<< ipv <<"STAT" << j << "\": { \"type\": \""<< OpNames[statOp[3+j][i]] << "\",";

                                    if (j==_NumPerVar[i]-1 && ipv==GroupVars[ig-1]-1 ) *(*p)->out  << "\"value\": [ " << result[0] << ", " << result[1] << ", " << result[2] << ", " << result[3] << " ] }";
                                    else *(*p)->out  << "\"value\": [ " << result[0] << ", " << result[1] << ", " << result[2] << ", " << result[3] << " ] },";
#else
                                    *(*p)->out << "\n# STAT " << statOp[3+j][i] << " :" << OpNames[statOp[3+j][i]] << "\n";
                                    *(*p)->out << result[0] << " " << result[1] << " " << result[2] << " " << result[3];
#endif
                                  }
                                }
                              }
                            }
                          }
                        }

                      }

                      //--- Gradients
                      for (register int n = 0; n <_NumGrads; n++) {
                        int na, nb;
                        /*std::cout << "CHECKING GRADIENT " << n << " ig " << ig
                                  << " GradGroup[n] " << GradGroup[n]
                                  << " ipv " << ipv
                                  << " ip " << ip
                                  << " (*g)->probe_table.size() " << (*g)->probe_table.size() << std::endl;
                                  */
                        if (GradGroup[n] == ig  && ip == 1 && ipv == 0) {  // once per group
                          for (register int i = 0; i < dim; i++) {
                            if (Grads[i][n]) {
                              calcount[n] += 2;
                            }
                          }

                        }
                        //std::cout << " ========== 0\n\n";

                        if (GradGroup[n] == ig  && ipv == 0) { // once per probe in group
                          if ( (*p)->type == TOK_SURFACE ) {
                            na = (*p)->options.surface.na;
                            nb = (*p)->options.surface.nb;
                          } else if ((*p)->type == TOK_PLANES ) {
                            na = (*p)->options.plane.na;
                            nb = (*p)->options.plane.nb;
                          }
                          /*printf( "VARIABLES = \"X\", \"Y\", \"Z\"\n");
                          printf( "ZONE I=%d, J=%d, F=POINT\n", na, nb);
                          for ( int j = 0 ; j < nb ; j++ )
                            for ( int i = 0 ; i < na ; i++ ) {
                              int l = i+j*na;
                              printf( "%12.8lf %12.8lf %12.8lf\n", (*p)->coords[l](0),
                                      (*p)->coords[l](1), (*p)->coords[l](2));
                              *out  << (*p)->coords[l](0) << " "
                                    << (*p)->coords[l](1) << " "
                                    << (*p)->coords[l](2) << "\n";
                            }*/

                          ctmp[ip-1].clear();
                          std::copy( &(*p)->coords[0], &(*p)->coords[na*nb-1], back_inserter(ctmp[ip-1]) );
                        }
                        //std::cout << " ========== 1\n\n";


                        if (Grads[(ip)/2][n] && GradGroup[n] == ig ) { //&& ipv > dim) { // once per variable per probe
                          if (ip==1) {
                            //std::cout << "+ ip "<< ip << " Grads["<<0<<"]["<<n<<"] " << Grads[0][n]
                            //        << " ipv " << ipv << "\n";
                          }
                          else {
                            dircheck[(ip-2)/2] = 1;
                            //std::cout << "+ ip "<< ip << " Grads["<<(ip-2)/2<<"]["<<n<<"] " << Grads[(ip-2)/2][n]
                            //        << " ipv " << ipv << "\n";
                          }

                          std::copy( adata, adata + Lpoints, back_inserter(ftmp[ip-1]) );
                          gcount[n]++;
                          //std::cout << "\n calcount["<<n<<"] " << calcount[n]
                          //        << " gcount["<<n<<"] " << gcount[n]
                          //      << " ip " << ip << " ipv " << ipv
                          //    << " GradVars[n] " << GradVars[n] << std::endl;

                          if ( gcount[n] == GradVars[n]*calcount[n] && ip == (*g)->probe_table.size() && ipv+1 == GradVars[n] ) {
#ifdef SJSON
                            CalcGradient(ig,dircheck,GradVars[n],Lpoints,na,nb,ctmp,ftmp,TimeC,t,curltog[n],wvtog[n],uwtog[n],vutog[n]);
#else
                            CalcGradient(ig,dircheck,GradVars[n],Lpoints,na,nb,ctmp,ftmp,TimeC,t,out,curltog[n],wvtog[n],uwtog[n],vutog[n]);
#endif
                            calcount[n] = 1;
                            gcount[n] = 0;
                            ipv=0;
                            for (register int cp=0; cp<2*dim+1; cp++)
                              ftmp[cp].clear();

                          }

                        }
                      }

#ifndef SJSON
                      *(*p)->out << std::endl;
#endif

                      delete [] result;
                    }
                  }
                  delete [] adata;
                }
                break;
              case TOK_VOLUME:
                break;
              }
              cli += Lpoints;
              clearcoords(*p);
              //ipv++;
            }
            ipv++;
            unlinksym(ptr->prev);
          }
#ifndef DAGH_NO_MPI
            MPI_Barrier(comm_service::comm());
#endif
            break;
          case TOK_FUNCTION1:
          case TOK_FUNCTION2:
          case TOK_VARINARRAY:
          case TOK_VARIABLE:
          case TOK_COORDINATE:
            staterror("there should not be any pointer type tokens"
                      " left in the stack");
            return;
          default:
            if ( ptr->type != TOK_DOUBLE ) {
              staterror("unhandled stack token");
              return;
            }
            break;
          }
          // need special code to handle the separator token, to flush the queue
          if ( ptr->type == TOK_SEPARATOR ) {
            ptr_ref = ptr->next;
            unlinksym(ptr);
            ptr = ptr_ref;
          } else {
            ptr->type = TOK_DOUBLE;
            ptr = ptr->next;
          }

#ifdef DEBUG_PRINT_STATISTICS
          (comm_service::log() << " end of token switch\n" ).flush();
#endif
        }
        delete [] lused;
        delete [] CUList;

        if (me == VizServer) {
          int ip = 1;
          for ( typename stparser::provec::const_iterator p = (*g)->probe_table.begin() ;
                p != (*g)->probe_table.end() ; p ++, ip++) {
#ifdef WJSON
            *out << "\t\t\t\"probe" << stparser::token((*p)->type) << ip << "\": {" << std::endl;
#elif !defined(WJSON) && !defined(SJSON)
            *out << "# probe " << stparser::token((*p)->type) << " " << ip << std::endl;
            ( *out << (*p)->out->str() ).flush();;
#endif
#ifdef WJSON
            if (ip<(*g)->probe_table.size()) *out << "\t\t\t},\n";
            else {
              *out << "\t\t\t}\n";
              if (ig<group_table.size()) *out << "\t\t},\n";
              //else *out << "\t\t}\n\t}\n}";
              else *out << "\t\t}\n}"; //newformat
            }
#endif
            delete (*p)->out;
          }
        }
      }
    }

    // close output file
#ifndef SJSON
    if (me == VizServer) {
      outfile.close();
      delete out;
    }
#endif
  }

private:
  int _DumpMeshes;
  int _Step;
  ControlDevice LocCtrl;
  std::string _OutputFile;
  std::string _DefFile;

  OutputType* _FileOutput;
  InterpolationType* _Interpolation;

  int _NumStats, _NumPerVar[MAXSTATOPS];
  int statOp[3+MAXPERVAR][MAXSTATOPS];

  int _NumGrads, GradGroup[MAXGRADS], GradVars[MAXGRADS], Grads[dim][MAXGRADS];
  std::vector<std::vector <double> > ftmp;
  std::vector<std::vector <point_type> > ctmp;
  int gcount[MAXGRADS], calcount[MAXGRADS], curltog[MAXGRADS], wvtog[MAXGRADS], uwtog[MAXGRADS], vutog[MAXGRADS];
  int dircheck[dim];

  std::vector<std::string> OpNames;
  int GroupVars[MAXGROUPS];

  int _SepFiles;
};

#endif

---