double compute_deviation()

in Documents/NamdSample/namd210/lib/abf_integrate/abf_integrate.cpp [201:293]

• 78 lines of code
• 19 McCabe index (conditional complexity)

double compute_deviation(ABFdata * data, bool meta, double kT)
{
    // Computing deviation between gradients differentiated from pmf
    // and input data
    // NOTE: this is mostly for output, hence NOT performance-critical
    double        *dev = data->deviation;
    double        *est = data->estimate;
    const double  *grad = data->gradients;
    int           *pos, *dpos, *newpos;
    double        rmsd = 0.0;
    unsigned int  offset, newoffset;
    double        sum;
    int           c;
    bool          moved;
    unsigned int  norm = 0; // number of data points summmed

    pos = new int[data->Nvars];
    dpos = new int[data->Nvars];
    newpos = new int[data->Nvars];

    for (int i = 0; i < data->Nvars; i++)
        pos[i] = 0;

    for (offset = 0; offset < data->scalar_dim; offset++) {
        for (int i = data->Nvars - 1; i > 0; i--) {
            if (pos[i] == data->sizes[i]) {
                pos[i] = 0;
                pos[i - 1]++;
            }
        }

        if (data->allowed (offset)) {
          for (int i = 0; i < data->Nvars; i++)
              newpos[i] = pos[i];

          for (int i = 0; i < data->Nvars; i++) {
              est[i] = 0.0;
              sum = 0.0;          // sum of finite differences on two sides (if not on edge of the grid)
              c = 0;              // count of summed values

              newpos[i] = pos[i] - 1;
              moved = data->wrap(newpos[i], i);
              newoffset = data->offset(newpos);
              if ( moved && data->allowed (newoffset) ) {
                  if (meta) {
                      sum = (data->bias[newoffset] - data->bias[offset]) / data->widths[i];
                      c++;
                  } else {
                      if (data->histogram[offset] && data->histogram[newoffset]) {
                          sum = kT * log(double (data->histogram[newoffset]) /
                                            double (data->histogram[offset])) / data->widths[i];
                          c++;
                      }
                  }
              }

              newpos[i] = pos[i] + 1;
              moved = data->wrap(newpos[i], i);
              newoffset = data->offset(newpos);
              if ( moved && data->allowed (newoffset) ) {
                  if (meta) {
                      sum += (data->bias[offset] - data->bias[newoffset]) / data->widths[i];
                      c++;
                  } else {
                      if (data->histogram[offset] && data->histogram[newoffset]) {
                          sum += kT * log(double (data->histogram[offset]) /
                                             double (data->histogram[newoffset])) / data->widths[i];
                          c++;
                      }
                  }
              }

              newpos[i] = pos[i]; // Go back to initial position for next dimension

              est[i] = (c ? sum/double(c) : 0.0);
              dev[i] = grad[i] - est[i];
              rmsd += dev[i] * dev[i];
              norm++;
          }
        }

        pos[data->Nvars - 1]++; // move on to next point
        est += data->Nvars;
        dev += data->Nvars;
        grad += data->Nvars;
    }

    delete [] pos;
    delete [] newpos;
    delete [] dpos;

    return sqrt(rmsd / norm);
}