Add new analysis to calculate finite difference

.travis.yml

@@ -2,16 +2,14 @@ language: cpp
 dist: trusty
 sudo: false
 addons:
+  sources:
+    - ubuntu-toolchain-r-test
   apt:
-    sources:
-      - ubuntu-toolchain-r-test
-      - george-edison55-precise-backports
     packages:
       - gfortran
       - libbz2-dev
       - libblas-dev
       - liblapack-dev
-      - libarpack2-dev
       - libnetcdf-dev
       - libfftw3-dev
       - netcdf-bin

doc/cpptraj.lyx

@@ -38804,11 +38804,79 @@ integrate
 \end_layout
 
 \begin_layout LyX-Code
-integrate <dset0> [<dset1> ...] [out <outfile>] [name <outsetname>]
+integrate <dset0> [<dset1> ...] [out <outfile>] [intout <intfile>]
+\end_layout
+
+\begin_layout LyX-Code
+          [name <name>]
+\end_layout
+
+\begin_deeper
+\begin_layout Description
+<dset0>
+\begin_inset space ~
+\end_inset
+
+[<dset1>
+\begin_inset space ~
+\end_inset
+
+...] Data set(s) to integrate.
+\end_layout
+
+\begin_layout Description
+[out
+\begin_inset space ~
+\end_inset
+
+<outfile>] If specified, write cumulative sum curves to <outfile>.
+\end_layout
+
+\begin_layout Description
+[intout
+\begin_inset space ~
+\end_inset
+
+<intfile>] If specified, write final integral values to <intfile>.
+\end_layout
+
+\begin_layout Description
+[name
+\begin_inset space ~
+\end_inset
+
+<name>] Output data set(s) name.
 \end_layout
 
+\begin_layout Standard
+DataSets Created:
+\end_layout
+
+\begin_layout Description
+<name> Final integral values, 1 for each input data set (indexed from 0).
+\end_layout
+
+\begin_layout Description
+<name>[Sum]:<idx> Cumulative sum curves if 
+\series bold
+out
+\series default
+ was specified, 1 for each input data set (indexed from 0).
+\end_layout
+
+\end_deeper
 \begin_layout Standard
 Integrate specified data set(s) using trapezoid integration.
+ If 
+\series bold
+'out'
+\series default
+ is specified write cumulative sum curves to <outfile>.
+ If 
+\series bold
+'intout'
+\series default
+ is specified write final integral values for each set to <intfile>.
 \end_layout
 
 \begin_layout Subsection
@@ -42251,6 +42319,70 @@ Calculate running average over windows of given size for data in selected
  data set(s).
 \end_layout
 
+\begin_layout Subsection
+slope
+\end_layout
+
+\begin_layout LyX-Code
+slope <dset0> [<dset1> ...] [out <outfile>] [name <name>]
+\end_layout
+
+\begin_layout LyX-Code
+  	[type {forward|backward|central}]
+\end_layout
+
+\begin_deeper
+\begin_layout Description
+<dset0>
+\begin_inset space ~
+\end_inset
+
+[<dset1>
+\begin_inset space ~
+\end_inset
+
+...] Data set(s) to calculate finite difference for.
+\end_layout
+
+\begin_layout Description
+[out
+\begin_inset space ~
+\end_inset
+
+<outfile>] File to write finite difference curves to.
+\end_layout
+
+\begin_layout Description
+[name
+\begin_inset space ~
+\end_inset
+
+<name>] Output data set(s) name.
+\end_layout
+
+\begin_layout Description
+[type
+\begin_inset space ~
+\end_inset
+
+{forward|backward|central}] Specify type of finite difference to calculate
+ (default forward).
+\end_layout
+
+\begin_layout Standard
+DataSets generated:
+\end_layout
+
+\begin_layout Description
+<name>:<idx> Output finite difference curves for each input data set (indexed
+ from 0).
+\end_layout
+
+\end_deeper
+\begin_layout Standard
+Calculate finite differences for each input data set.
+\end_layout
+
 \begin_layout Subsection
 spline
 \end_layout
@@ -42310,7 +42442,7 @@ spline <dset0> [<dset1> ...] [out <outfile>] [meshsize <n> | meshfactor <x>]
 
 \end_deeper
 \begin_layout Standard
-Cubic spline the given data sets.
+Apply cubic splines to the given input data sets to create new data sets.
 \end_layout
 
 \begin_layout Subsection

src/Action_VelocityAutoCorr.cpp

@@ -1,7 +1,6 @@
 #include "Action_VelocityAutoCorr.h"
 #include "CpptrajStdio.h"
 #include "ProgressBar.h"
-#include "DataSet_Mesh.h"
 #include "DataSet_double.h"
 #include "Constants.h"
 #include "Corr.h"
@@ -281,9 +280,7 @@ void Action_VelocityAutoCorr::Print() {
   // Integration to get diffusion coefficient.
   VAC_->SetDim(Dimension::X, Dimension(0.0, tstep_, "Time (ps)"));
   mprintf("\tIntegrating data set %s, step is %f\n", VAC_->legend(), VAC_->Dim(0).Step());
-  DataSet_Mesh mesh;
-  mesh.SetMeshXY( static_cast<DataSet_1D const&>(*VAC_) );
-  double total = mesh.Integrate_Trapezoid();
+  double total = Ct.Integrate( DataSet_1D::TRAPEZOID );
   const double ANG2_PS_TO_CM2_S = 10.0; // Convert Ang^2/ps to 1E-5 cm^2/s
   const char* tab = "\t";
   if (!diffout_->IsStream()) {

src/Analysis_Integrate.cpp

@@ -1,17 +1,22 @@
 #include "Analysis_Integrate.h"
 #include "CpptrajStdio.h"
+#include "DataSet_Mesh.h"
 
-Analysis_Integrate::Analysis_Integrate() : outfile_(0) {}
+Analysis_Integrate::Analysis_Integrate() : sumSet_(0) {}
 
 void Analysis_Integrate::Help() const {
-  mprintf("\t<dset0> [<dset1> ...] [out <outfile>] [name <outsetname>]\n"
-          "  Integrate given data sets.\n");
+  mprintf("\t<dset0> [<dset1> ...] [out <outfile>] [intout <intfile>]\n"
+          "\t[name <name>]\n"
+          "  Integrate given data sets. If 'out' is specified write cumulative\n"
+          "  sum curves to <outfile>. If 'intout' is specified write final\n"
+          "  integral values for each set to <intfile>.\n");
 }
 
 Analysis::RetType Analysis_Integrate::Setup(ArgList& analyzeArgs, AnalysisSetup& setup, int debugIn)
 {
   std::string setname = analyzeArgs.GetStringKey("name");
-  outfile_ = setup.DFL().AddDataFile(analyzeArgs.GetStringKey("out"), analyzeArgs);
+  DataFile* outfile = setup.DFL().AddDataFile(analyzeArgs.GetStringKey("out"), analyzeArgs);
+  DataFile* intfile = setup.DFL().AddDataFile(analyzeArgs.GetStringKey("intout"), analyzeArgs);
   // Select datasets from remaining args
   if (input_dsets_.AddSetsFromArgs( analyzeArgs.RemainingArgs(), setup.DSL() )) {
     mprinterr("Error: Could not add data sets.\n");
@@ -23,48 +28,55 @@ Analysis::RetType Analysis_Integrate::Setup(ArgList& analyzeArgs, AnalysisSetup&
   }
 
   // Set up output datasets
-  if (outfile_ != 0) {
+  if (setname.empty())
+    setname = setup.DSL().GenerateDefaultName("Int");
+  sumSet_ = setup.DSL().AddSet(DataSet::DOUBLE, setname);
+  if (sumSet_ == 0) return Analysis::ERR;
+  if (intfile != 0) intfile->AddDataSet( sumSet_ );
+  if (outfile != 0) {
+    int idx = 0;
+    MetaData md(setname, "Sum");
     for (Array1D::const_iterator dsIn = input_dsets_.begin();
                                  dsIn != input_dsets_.end(); ++dsIn)
     {
-      DataSet* ds = setup.DSL().AddSet(DataSet::XYMESH, setname, "Int");
+      md.SetIdx(idx++);
+      DataSet* ds = setup.DSL().AddSet(DataSet::XYMESH, md);
       if (ds == 0) return Analysis::ERR;
       ds->SetLegend( "Int(" + (*dsIn)->Meta().Legend() + ")" );
-      outfile_->AddDataSet( ds );
+      outfile->AddDataSet( ds );
       output_dsets_.push_back( (DataSet_Mesh*)ds );
     }
   }
 
-  mprintf("    INTEGRATE: Calculating integral of %zu data sets.\n",
+  mprintf("    INTEGRATE: Calculating integral for %zu data sets.\n",
           input_dsets_.size());
-  if (outfile_ != 0) {
+  if (outfile != 0) {
     if (!setname.empty())
       mprintf("\tOutput set name: %s\n", setname.c_str());
-    mprintf("\tOutfile name: %s\n", outfile_->DataFilename().base());
+    mprintf("\tOutfile name: %s\n", outfile->DataFilename().base());
+  }
+  if (debugIn > 0) {
+    for (Array1D::const_iterator set = input_dsets_.begin(); set != input_dsets_.end(); ++set)
+      mprintf("\t%s\n", (*set)->legend());
   }
-  //for (Array1D::const_iterator set = input_dsets_.begin(); set != input_dsets_.end(); ++set)
-  //  mprintf("\t%s\n", (*set)->legend());
   return Analysis::OK;
 }
 
 Analysis::RetType Analysis_Integrate::Analyze() {
   double sum;
-  int idx = 0;
-  for (Array1D::const_iterator DS = input_dsets_.begin();
-                               DS != input_dsets_.end(); ++DS, ++idx)
-  {
-    if ( (*DS)->Size() < 1)
-      mprintf("Warning: Set [%i] \"%s\" has no data.\n", idx, (*DS)->legend());
-    else {
-      DataSet_Mesh mesh;
-      // Set XY mesh
-      mesh.SetMeshXY( *(*DS) );
-      if (outfile_ != 0)
-        sum = mesh.Integrate_Trapezoid( *(output_dsets_[idx]) );
-      else
-        sum = mesh.Integrate_Trapezoid();
-      mprintf("\tIntegral of %s is %g\n", (*DS)->legend(), sum);
-      output_dsets_[idx]->SetDim(Dimension::X, (*DS)->Dim(0));
+  for (unsigned int idx = 0; idx != input_dsets_.size(); idx++) {
+    DataSet_1D const& inSet = *(input_dsets_[idx]);
+    if (inSet.Size() < 1) {
+      mprintf("Warning: Set '%s' has no data.\n", inSet.legend());
+    } else {
+      if (!output_dsets_.empty()) {
+        DataSet_Mesh& outSet = static_cast<DataSet_Mesh&>( *(output_dsets_[idx]) );
+        sum = inSet.Integrate( DataSet_1D::TRAPEZOID, outSet.SetMeshX(), outSet.SetMeshY() );
+        outSet.SetDim(Dimension::X, inSet.Dim(0));
+      } else
+        sum = inSet.Integrate( DataSet_1D::TRAPEZOID );
+      mprintf("\tIntegral of %s is %g\n", inSet.legend(), sum);
+      sumSet_->Add(idx, &sum);
     }
   }
   return Analysis::OK;

src/Analysis_Integrate.h

@@ -2,7 +2,7 @@
 #define INC_ANALYSIS_INTEGRATE_H
 #include "Analysis.h"
 #include "Array1D.h"
-#include "DataSet_Mesh.h"
+class DataSet_Mesh;
 class Analysis_Integrate : public Analysis {
   public:
     Analysis_Integrate();
@@ -12,8 +12,8 @@ class Analysis_Integrate : public Analysis {
     Analysis::RetType Setup(ArgList&, AnalysisSetup&, int);
     Analysis::RetType Analyze();
   private:
-    DataFile* outfile_; // FIXME: May not need to be class var
     Array1D input_dsets_;
     std::vector<DataSet_Mesh*> output_dsets_;
+    DataSet* sumSet_; ///< Hold final sum for each set
 };
 #endif

src/Analysis_Rotdif.cpp

@@ -1590,7 +1590,7 @@ int Analysis_Rotdif::DetermineDeffs() {
     // Calculate mesh Y values
     spline.SetSplinedMeshY(pX, pY);
     // Integrate
-    double integral = spline.Integrate_Trapezoid();
+    double integral = spline.Integrate( DataSet_1D::TRAPEZOID );
     //mprintf("DEBUG: Vec %i integral= %g\n", nvec, integral);
     // Solve for deff
     D_eff_.push_back( calcEffectiveDiffusionConst(integral) );