Skip to content

Commit 8ad769d

Browse files
drroedrroe
authored
Add new analysis to calculate finite difference (#753)
* DRR - Cpptraj: Move integration routines from DataSet_Mesh to DataSet_1D * DRR - Cpptraj: Use forward declare. * DRR - Cpptraj: Add some finite difference calcs * DRR - Cpptraj: Add functions to header * DRR - Cpptraj: Warn when incoming X mesh doesnt have the right size * DRR - Cpptraj: Forward declare mesh * DRR - Cpptraj: Finite difference analysis class * DRR - Cpptraj: Add slope command. Have finite difference routines set X values directly. * DRR - Cpptraj: Add finite difference test * DRR - Cpptraj: Have integrate routine with cumulative sum set X values. * DRR - Cpptraj: Slight reorganization of loop for clarity * DRR - Cpptraj: Save final integral values to a data set; add 'intout' keyword to write that to a file. * DRR - Cpptraj: Compare final integral value. * DRR - Cpptraj: Forgot to actually put the test in. * DRR - Cpptraj: Better output. * DRR - Cpptraj: Revision bump for new slope (finite difference) analysis command and new integrate option 'intout'. * DRR - Cpptraj: Enable slope test * DRR - Cpptraj: Simplify help * DRR - Cpptraj: Update manual for new slope command and new integrate command option intout (also fix up the integrate command entry in general). * DRR - Cpptraj: Remove sources to try to fix travis builds. Try to add parallel netcdf to MPI build. * DRR - Cpptraj: Need to have the ubuntu-toolchain ppa in there. * DRR - Cpptraj: Turns out pnetcdf not part of the toolchain
1 parent c38a771 commit 8ad769d

24 files changed

+778
-96
lines changed

‎.travis.yml‎

Lines changed: 2 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -2,16 +2,14 @@ language: cpp
22
dist: trusty
33
sudo: false
44
addons:
5+
sources:
6+
- ubuntu-toolchain-r-test
57
apt:
6-
sources:
7-
- ubuntu-toolchain-r-test
8-
- george-edison55-precise-backports
98
packages:
109
- gfortran
1110
- libbz2-dev
1211
- libblas-dev
1312
- liblapack-dev
14-
- libarpack2-dev
1513
- libnetcdf-dev
1614
- libfftw3-dev
1715
- netcdf-bin

‎doc/CpptrajManual.pdf‎

1.55 KB
Binary file not shown.

‎doc/cpptraj.lyx‎

Lines changed: 134 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -38804,11 +38804,79 @@ integrate
3880438804
\end_layout
3880538805

3880638806
\begin_layout LyX-Code
38807-
integrate <dset0> [<dset1> ...] [out <outfile>] [name <outsetname>]
38807+
integrate <dset0> [<dset1> ...] [out <outfile>] [intout <intfile>]
38808+
\end_layout
38809+
38810+
\begin_layout LyX-Code
38811+
[name <name>]
38812+
\end_layout
38813+
38814+
\begin_deeper
38815+
\begin_layout Description
38816+
<dset0>
38817+
\begin_inset space ~
38818+
\end_inset
38819+
38820+
[<dset1>
38821+
\begin_inset space ~
38822+
\end_inset
38823+
38824+
...] Data set(s) to integrate.
38825+
\end_layout
38826+
38827+
\begin_layout Description
38828+
[out
38829+
\begin_inset space ~
38830+
\end_inset
38831+
38832+
<outfile>] If specified, write cumulative sum curves to <outfile>.
38833+
\end_layout
38834+
38835+
\begin_layout Description
38836+
[intout
38837+
\begin_inset space ~
38838+
\end_inset
38839+
38840+
<intfile>] If specified, write final integral values to <intfile>.
38841+
\end_layout
38842+
38843+
\begin_layout Description
38844+
[name
38845+
\begin_inset space ~
38846+
\end_inset
38847+
38848+
<name>] Output data set(s) name.
3880838849
\end_layout
3880938850

38851+
\begin_layout Standard
38852+
DataSets Created:
38853+
\end_layout
38854+
38855+
\begin_layout Description
38856+
<name> Final integral values, 1 for each input data set (indexed from 0).
38857+
\end_layout
38858+
38859+
\begin_layout Description
38860+
<name>[Sum]:<idx> Cumulative sum curves if
38861+
\series bold
38862+
out
38863+
\series default
38864+
was specified, 1 for each input data set (indexed from 0).
38865+
\end_layout
38866+
38867+
\end_deeper
3881038868
\begin_layout Standard
3881138869
Integrate specified data set(s) using trapezoid integration.
38870+
If
38871+
\series bold
38872+
'out'
38873+
\series default
38874+
is specified write cumulative sum curves to <outfile>.
38875+
If
38876+
\series bold
38877+
'intout'
38878+
\series default
38879+
is specified write final integral values for each set to <intfile>.
3881238880
\end_layout
3881338881

3881438882
\begin_layout Subsection
@@ -42251,6 +42319,70 @@ Calculate running average over windows of given size for data in selected
4225142319
data set(s).
4225242320
\end_layout
4225342321

42322+
\begin_layout Subsection
42323+
slope
42324+
\end_layout
42325+
42326+
\begin_layout LyX-Code
42327+
slope <dset0> [<dset1> ...] [out <outfile>] [name <name>]
42328+
\end_layout
42329+
42330+
\begin_layout LyX-Code
42331+
[type {forward|backward|central}]
42332+
\end_layout
42333+
42334+
\begin_deeper
42335+
\begin_layout Description
42336+
<dset0>
42337+
\begin_inset space ~
42338+
\end_inset
42339+
42340+
[<dset1>
42341+
\begin_inset space ~
42342+
\end_inset
42343+
42344+
...] Data set(s) to calculate finite difference for.
42345+
\end_layout
42346+
42347+
\begin_layout Description
42348+
[out
42349+
\begin_inset space ~
42350+
\end_inset
42351+
42352+
<outfile>] File to write finite difference curves to.
42353+
\end_layout
42354+
42355+
\begin_layout Description
42356+
[name
42357+
\begin_inset space ~
42358+
\end_inset
42359+
42360+
<name>] Output data set(s) name.
42361+
\end_layout
42362+
42363+
\begin_layout Description
42364+
[type
42365+
\begin_inset space ~
42366+
\end_inset
42367+
42368+
{forward|backward|central}] Specify type of finite difference to calculate
42369+
(default forward).
42370+
\end_layout
42371+
42372+
\begin_layout Standard
42373+
DataSets generated:
42374+
\end_layout
42375+
42376+
\begin_layout Description
42377+
<name>:<idx> Output finite difference curves for each input data set (indexed
42378+
from 0).
42379+
\end_layout
42380+
42381+
\end_deeper
42382+
\begin_layout Standard
42383+
Calculate finite differences for each input data set.
42384+
\end_layout
42385+
4225442386
\begin_layout Subsection
4225542387
spline
4225642388
\end_layout
@@ -42310,7 +42442,7 @@ spline <dset0> [<dset1> ...] [out <outfile>] [meshsize <n> | meshfactor <x>]
4231042442

4231142443
\end_deeper
4231242444
\begin_layout Standard
42313-
Cubic spline the given data sets.
42445+
Apply cubic splines to the given input data sets to create new data sets.
4231442446
\end_layout
4231542447

4231642448
\begin_layout Subsection

‎src/Action_VelocityAutoCorr.cpp‎

Lines changed: 1 addition & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,6 @@
11
#include "Action_VelocityAutoCorr.h"
22
#include "CpptrajStdio.h"
33
#include "ProgressBar.h"
4-
#include "DataSet_Mesh.h"
54
#include "DataSet_double.h"
65
#include "Constants.h"
76
#include "Corr.h"
@@ -281,9 +280,7 @@ void Action_VelocityAutoCorr::Print() {
281280
// Integration to get diffusion coefficient.
282281
VAC_->SetDim(Dimension::X, Dimension(0.0, tstep_, "Time (ps)"));
283282
mprintf("\tIntegrating data set %s, step is %f\n", VAC_->legend(), VAC_->Dim(0).Step());
284-
DataSet_Mesh mesh;
285-
mesh.SetMeshXY( static_cast<DataSet_1D const&>(*VAC_) );
286-
double total = mesh.Integrate_Trapezoid();
283+
double total = Ct.Integrate( DataSet_1D::TRAPEZOID );
287284
const double ANG2_PS_TO_CM2_S = 10.0; // Convert Ang^2/ps to 1E-5 cm^2/s
288285
const char* tab = "\t";
289286
if (!diffout_->IsStream()) {

‎src/Analysis_Integrate.cpp‎

Lines changed: 40 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -1,17 +1,22 @@
11
#include "Analysis_Integrate.h"
22
#include "CpptrajStdio.h"
3+
#include "DataSet_Mesh.h"
34

4-
Analysis_Integrate::Analysis_Integrate() : outfile_(0) {}
5+
Analysis_Integrate::Analysis_Integrate() : sumSet_(0) {}
56

67
void Analysis_Integrate::Help() const {
7-
mprintf("\t<dset0> [<dset1> ...] [out <outfile>] [name <outsetname>]\n"
8-
" Integrate given data sets.\n");
8+
mprintf("\t<dset0> [<dset1> ...] [out <outfile>] [intout <intfile>]\n"
9+
"\t[name <name>]\n"
10+
" Integrate given data sets. If 'out' is specified write cumulative\n"
11+
" sum curves to <outfile>. If 'intout' is specified write final\n"
12+
" integral values for each set to <intfile>.\n");
913
}
1014

1115
Analysis::RetType Analysis_Integrate::Setup(ArgList& analyzeArgs, AnalysisSetup& setup, int debugIn)
1216
{
1317
std::string setname = analyzeArgs.GetStringKey("name");
14-
outfile_ = setup.DFL().AddDataFile(analyzeArgs.GetStringKey("out"), analyzeArgs);
18+
DataFile* outfile = setup.DFL().AddDataFile(analyzeArgs.GetStringKey("out"), analyzeArgs);
19+
DataFile* intfile = setup.DFL().AddDataFile(analyzeArgs.GetStringKey("intout"), analyzeArgs);
1520
// Select datasets from remaining args
1621
if (input_dsets_.AddSetsFromArgs( analyzeArgs.RemainingArgs(), setup.DSL() )) {
1722
mprinterr("Error: Could not add data sets.\n");
@@ -23,48 +28,55 @@ Analysis::RetType Analysis_Integrate::Setup(ArgList& analyzeArgs, AnalysisSetup&
2328
}
2429

2530
// Set up output datasets
26-
if (outfile_ != 0) {
31+
if (setname.empty())
32+
setname = setup.DSL().GenerateDefaultName("Int");
33+
sumSet_ = setup.DSL().AddSet(DataSet::DOUBLE, setname);
34+
if (sumSet_ == 0) return Analysis::ERR;
35+
if (intfile != 0) intfile->AddDataSet( sumSet_ );
36+
if (outfile != 0) {
37+
int idx = 0;
38+
MetaData md(setname, "Sum");
2739
for (Array1D::const_iterator dsIn = input_dsets_.begin();
2840
dsIn != input_dsets_.end(); ++dsIn)
2941
{
30-
DataSet* ds = setup.DSL().AddSet(DataSet::XYMESH, setname, "Int");
42+
md.SetIdx(idx++);
43+
DataSet* ds = setup.DSL().AddSet(DataSet::XYMESH, md);
3144
if (ds == 0) return Analysis::ERR;
3245
ds->SetLegend( "Int(" + (*dsIn)->Meta().Legend() + ")" );
33-
outfile_->AddDataSet( ds );
46+
outfile->AddDataSet( ds );
3447
output_dsets_.push_back( (DataSet_Mesh*)ds );
3548
}
3649
}
3750

38-
mprintf(" INTEGRATE: Calculating integral of %zu data sets.\n",
51+
mprintf(" INTEGRATE: Calculating integral for %zu data sets.\n",
3952
input_dsets_.size());
40-
if (outfile_ != 0) {
53+
if (outfile != 0) {
4154
if (!setname.empty())
4255
mprintf("\tOutput set name: %s\n", setname.c_str());
43-
mprintf("\tOutfile name: %s\n", outfile_->DataFilename().base());
56+
mprintf("\tOutfile name: %s\n", outfile->DataFilename().base());
57+
}
58+
if (debugIn > 0) {
59+
for (Array1D::const_iterator set = input_dsets_.begin(); set != input_dsets_.end(); ++set)
60+
mprintf("\t%s\n", (*set)->legend());
4461
}
45-
//for (Array1D::const_iterator set = input_dsets_.begin(); set != input_dsets_.end(); ++set)
46-
// mprintf("\t%s\n", (*set)->legend());
4762
return Analysis::OK;
4863
}
4964

5065
Analysis::RetType Analysis_Integrate::Analyze() {
5166
double sum;
52-
int idx = 0;
53-
for (Array1D::const_iterator DS = input_dsets_.begin();
54-
DS != input_dsets_.end(); ++DS, ++idx)
55-
{
56-
if ( (*DS)->Size() < 1)
57-
mprintf("Warning: Set [%i] \"%s\" has no data.\n", idx, (*DS)->legend());
58-
else {
59-
DataSet_Mesh mesh;
60-
// Set XY mesh
61-
mesh.SetMeshXY( *(*DS) );
62-
if (outfile_ != 0)
63-
sum = mesh.Integrate_Trapezoid( *(output_dsets_[idx]) );
64-
else
65-
sum = mesh.Integrate_Trapezoid();
66-
mprintf("\tIntegral of %s is %g\n", (*DS)->legend(), sum);
67-
output_dsets_[idx]->SetDim(Dimension::X, (*DS)->Dim(0));
67+
for (unsigned int idx = 0; idx != input_dsets_.size(); idx++) {
68+
DataSet_1D const& inSet = *(input_dsets_[idx]);
69+
if (inSet.Size() < 1) {
70+
mprintf("Warning: Set '%s' has no data.\n", inSet.legend());
71+
} else {
72+
if (!output_dsets_.empty()) {
73+
DataSet_Mesh& outSet = static_cast<DataSet_Mesh&>( *(output_dsets_[idx]) );
74+
sum = inSet.Integrate( DataSet_1D::TRAPEZOID, outSet.SetMeshX(), outSet.SetMeshY() );
75+
outSet.SetDim(Dimension::X, inSet.Dim(0));
76+
} else
77+
sum = inSet.Integrate( DataSet_1D::TRAPEZOID );
78+
mprintf("\tIntegral of %s is %g\n", inSet.legend(), sum);
79+
sumSet_->Add(idx, &sum);
6880
}
6981
}
7082
return Analysis::OK;

‎src/Analysis_Integrate.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22
#define INC_ANALYSIS_INTEGRATE_H
33
#include "Analysis.h"
44
#include "Array1D.h"
5-
#include "DataSet_Mesh.h"
5+
class DataSet_Mesh;
66
class Analysis_Integrate : public Analysis {
77
public:
88
Analysis_Integrate();
@@ -12,8 +12,8 @@ class Analysis_Integrate : public Analysis {
1212
Analysis::RetType Setup(ArgList&, AnalysisSetup&, int);
1313
Analysis::RetType Analyze();
1414
private:
15-
DataFile* outfile_; // FIXME: May not need to be class var
1615
Array1D input_dsets_;
1716
std::vector<DataSet_Mesh*> output_dsets_;
17+
DataSet* sumSet_; ///< Hold final sum for each set
1818
};
1919
#endif

‎src/Analysis_Rotdif.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1590,7 +1590,7 @@ int Analysis_Rotdif::DetermineDeffs() {
15901590
// Calculate mesh Y values
15911591
spline.SetSplinedMeshY(pX, pY);
15921592
// Integrate
1593-
double integral = spline.Integrate_Trapezoid();
1593+
double integral = spline.Integrate( DataSet_1D::TRAPEZOID );
15941594
//mprintf("DEBUG: Vec %i integral= %g\n", nvec, integral);
15951595
// Solve for deff
15961596
D_eff_.push_back( calcEffectiveDiffusionConst(integral) );

0 commit comments

Comments
 (0)