Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

Refs #5133 Removed tied variables to get sensible errors

1. Errors reported when variables are tied are way too big
2. Added constraints to restrict range for the variables
3. Default now is to NOT calculate variances( Used variables tied with a complex formula)
4. Improved initial parameters a bit.
5. Updated test programs to reflect new values.
6. Improved derivatives in the case (co)variances were calculated.

7. Added a new column to the output table for IntegratePeakTimeSlices to list the
   specID's used to calcualte the integrated intensities.
  • Loading branch information...
commit 98b8459ce1b2dcc18f12764a5cb8bddec15fde35 1 parent 2563eb1
@RuthFromDuluth RuthFromDuluth authored
View
18 Code/Mantid/Framework/Crystal/inc/MantidCrystal/IntegratePeakTimeSlices.h
@@ -126,7 +126,8 @@ class DLLExport IntegratePeakTimeSlices: public Mantid::API::Algorithm
Kernel::V3D &CentNghbr,
double &neighborRadius,//from CentDet
- double Radius);
+ double Radius,
+ std::string &spec_idList);
@@ -148,12 +149,12 @@ class DLLExport IntegratePeakTimeSlices: public Mantid::API::Algorithm
void InitializeColumnNamesInTableWorkspace( DataObjects::TableWorkspace_sptr &TabWS) ;
- void SetUpData1( API::MatrixWorkspace_sptr &Data,
- API::MatrixWorkspace_sptr const &inpWkSpace,
- const int chan,
-
- double Radius,
- Kernel::V3D CentPos ///< Center on Plane
+ void SetUpData1( API::MatrixWorkspace_sptr &Data,
+ API::MatrixWorkspace_sptr const &inpWkSpace,
+ const int chan,
+ double Radius,
+ Kernel::V3D CentPos , ///< Center on Plane,
+ std::string &spec_idList
) ;
@@ -173,7 +174,8 @@ class DLLExport IntegratePeakTimeSlices: public Mantid::API::Algorithm
std::vector<double > const &errs,
std::vector<std::string> const &names,
const double chisq,
- const double time) ;
+ const double time,
+ std::string spec_idList) ;
void updatePeakInformation( std::vector<double > const &params,
View
57 Code/Mantid/Framework/Crystal/src/IntegratePeakTimeSlices.cpp
@@ -56,14 +56,14 @@ This Algorithm is also used by the [[PeakIntegration]] algorithm when the Fit ta
#include "MantidKernel/Unit.h"
#include "MantidKernel/V3D.h"
#include "MantidGeometry/Surfaces/Surface.h"
-//
+#include <boost/lexical_cast.hpp>
#include <vector>
#include "MantidAPI/Algorithm.h"
#include <algorithm>
#include <math.h>
#include <cstdio>
-#include <boost/random/poisson_distribution.hpp>
-
+//#include <boost/random/poisson_distribution.hpp>
+#include "MantidAPI/ISpectrum.h"
using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::DataObjects;
@@ -161,6 +161,8 @@ namespace Mantid
for (int i = 0; i < NParameters; i++)
ParameterValues[i] = 0;
+
+
}
/// Destructor
@@ -418,6 +420,7 @@ namespace Mantid
double Row0 = lastRow;
double lastCol = COL;
double Col0 = lastCol;
+ string spec_idList="";
// For quickly looking up workspace index from det id
wi_to_detid_map = inpWkSpace->getDetectorIDToWorkspaceIndexMap( false );
@@ -505,7 +508,7 @@ namespace Mantid
std::string("Workspace2D"), 3,NN,NN);
Kernel::V3D CentPos = center+yvec*(Centy-ROW)*CellHeight + xvec*(Centx-COL)*CellWidth;
- SetUpData1(Data, inpWkSpace, Chan+dir*t, R ,CenterDet->getPos());
+ SetUpData1(Data, inpWkSpace, Chan+dir*t, R ,CenterDet->getPos(), spec_idList);
if (AttributeValues[ISSIxx] > 0)
{
if (AttributeValues[IIntensities] > MaxCounts)
@@ -598,7 +601,7 @@ namespace Mantid
SetUpData(Data, inpWkSpace, panel, xchan, lastCol,lastRow,Cent,//CentDetspec,//NeighborIDs,
- neighborRadius, Radius);
+ neighborRadius, Radius, spec_idList);
Data->setName("index0");
@@ -647,7 +650,7 @@ namespace Mantid
double chisq = fit_alg->getProperty("OutputChi2overDoF");
ITableWorkspace_sptr RRes = fit_alg->getProperty( "OutputParameters");
- for( int prm = 0; prm < (int)RRes->rowCount(); prm++ )
+ for( int prm = 0; prm < (int)RRes->rowCount()-1; prm++ )
{
names.push_back( RRes->getRef< string >( "Name", prm ) );
params.push_back( RRes->getRef< double >( "Value", prm ));
@@ -656,11 +659,14 @@ namespace Mantid
ostringstream res;
res << " Thru Algorithm: chiSq=" << setw(7) << chisq << endl;
- res<<" Row,Col Radius="<<lastRow<<","<<lastCol<<","<<neighborRadius<<std::endl;
+ res<<" Row,Col Radius=" << lastRow << "," << lastCol <<"," <<neighborRadius << std::endl;
double sqrtChisq = -1;
if (chisq >= 0)
- sqrtChisq = sqrt(chisq);
+ sqrtChisq =(chisq);
+
+ sqrtChisq = max< double >( sqrtChisq,AttributeValues[IIntensities]/AttributeValues[ISS1] );
+ sqrtChisq = sqrt( sqrtChisq);
for (size_t kk = 0; kk < params.size(); kk++)
{
@@ -674,7 +680,7 @@ namespace Mantid
if (isGoodFit(params, errs, names, chisq))
{
- UpdateOutputWS(TabWS, dir, xchan, params, errs, names, chisq, time);
+ UpdateOutputWS(TabWS, dir, xchan, params, errs, names, chisq, time, spec_idList);
double TotSliceIntensity = AttributeValues[IIntensities];
double TotSliceVariance = AttributeValues[IVariance];
@@ -851,6 +857,7 @@ namespace Mantid
TabWS->addColumn("double", "Start Col");
TabWS->addColumn("double", "End Col");
TabWS->addColumn("double", "TotIntensityError");
+ TabWS->addColumn("str", "SpecIDs");
}
@@ -1024,7 +1031,8 @@ namespace Mantid
//int* nghbrs,//Not used
// std::map< specid_t, V3D > &neighbors,
double &neighborRadius,//from CentDetspec
- double Radius)
+ double Radius,
+ string &spec_idList)
{
// size_t wsIndx = inpWkSpace->getIndexFromSpectrumNumber(CentDetspec);
@@ -1034,13 +1042,13 @@ namespace Mantid
+ yvec*(CentY-ROW)*CellHeight;
SetUpData1(Data , inpWkSpace, chan,
- Radius , CentPos1);
+ Radius , CentPos1, spec_idList);
if( AttributeValues[ISSIxx]< 0)// Not enough data
return;
double DD = max<double>( sqrt(ParameterValues[IVYY])*CellHeight, sqrt(ParameterValues[IVXX])*CellWidth);
- double NewRadius= 1.1* max<double>( 5*max<double>(CellWidth,CellHeight), 4*DD);
+ double NewRadius= 1.4* max<double>( 5*max<double>(CellWidth,CellHeight), 4*DD);
//1.4 is needed to get more background cells. In rectangle the corners were background
NewRadius = min<double>(30*max<double>(CellWidth,CellHeight),NewRadius);
@@ -1081,7 +1089,7 @@ namespace Mantid
neighborRadius -= DD;
SetUpData1(Data, inpWkSpace, chan,
- NewRadius, CentPos );
+ NewRadius, CentPos, spec_idList );
}
@@ -1091,7 +1099,8 @@ namespace Mantid
API::MatrixWorkspace_sptr const &inpWkSpace,
const int chan,
double Radius,
- Kernel::V3D CentPos
+ Kernel::V3D CentPos,
+ string &spec_idList
)
{
UNUSED_ARG(g_log);
@@ -1104,7 +1113,7 @@ namespace Mantid
boost::shared_ptr<Workspace2D> ws = boost::shared_dynamic_cast<Workspace2D>(Data);
std::vector<double> StatBase;
-
+ spec_idList.clear();
for (int i = 0; i < NAttributes + 2; i++)
StatBase.push_back(0);
@@ -1137,6 +1146,10 @@ namespace Mantid
// }else
// std::cout<<"-------- data-------"<<std::endl;
+ int jj=0;
+ Mantid::MantidVecPtr pX;
+
+ Mantid::MantidVec& xRef = pX.access();
for (int i = 2; i < NeighborIDs[1]; i++)
{
int DetID = NeighborIDs[i];
@@ -1213,15 +1226,6 @@ namespace Mantid
}
}
-<<<<<<< Updated upstream
- Mantid::MantidVecPtr pX;
-
- Mantid::MantidVec& xRef = pX.access();
- for (int j = 0; j < N+1; j++)
- {
- xRef.push_back((double)j);
- }
-
Data->setX(0, pX);
Data->setX(1, pX);
@@ -1419,7 +1423,8 @@ namespace Mantid
std::vector<double > const &errs,
std::vector<std::string> const &names,
const double Chisq,
- const double time)
+ const double time,
+ string spec_idList)
{
int Ibk = find("Background", names);
int IIntensity = find("Intensity", names);
@@ -1477,7 +1482,7 @@ namespace Mantid
TabWS->getRef<double> (std::string("End Col"), TableRow) = AttributeValues[IStartCol]
+ AttributeValues[INCol] - 1;
TabWS->getRef<double> (std::string("TotIntensityError"), TableRow) = sqrt(AttributeValues[IVariance]);
-
+ TabWS->getRef<string>( std::string("SpecIDs"), TableRow ) = spec_idList;
}
View
59 Code/Mantid/Framework/Crystal/test/IntegratePeakTimeSlicesTest.h
@@ -194,9 +194,9 @@ class IntegratePeakTimeSlicesTest: public CxxTest::TestSuite
TableWorkspace_sptr Twk = algP.getProperty("OutputWorkspace");
- TS_ASSERT_LESS_THAN(fabs(intensity -59869.3), 100.0);
+ TS_ASSERT_LESS_THAN(fabs(intensity -60340.5), 100.0);
//Not sure why this reduced the error so much in the test
- TS_ASSERT_LESS_THAN(fabs(sigma -538), 1.0);
+ TS_ASSERT_LESS_THAN(fabs(sigma -466), 1.0);
TS_ASSERT_EQUALS( Twk->rowCount(), 7);
@@ -206,23 +206,23 @@ class IntegratePeakTimeSlicesTest: public CxxTest::TestSuite
TS_ASSERT_LESS_THAN(fabs(Twk->getRef<double> (std::string("Time"), 0) - 19250), 20);
- TS_ASSERT_LESS_THAN(fabs(Twk->getRef<double> (std::string("Background"), 1) - 1.40532 ), .2);
+ TS_ASSERT_LESS_THAN(fabs(Twk->getRef<double> (std::string("Background"), 1) - 1.23611 ), .2);
- TS_ASSERT_LESS_THAN(fabs(Twk->getRef<double> ("Intensity", 2) - 11238.6 ), 20);
+ TS_ASSERT_LESS_THAN(fabs(Twk->getRef<double> ("Intensity", 2) - 11331), 20);
- TS_ASSERT_LESS_THAN(fabs(Twk->getRef<double> ("NCells", 3) - 457), 5);
+ TS_ASSERT_LESS_THAN(fabs(Twk->getRef<double> ("NCells", 3) - 437), 5);
- TS_ASSERT_LESS_THAN(fabs(Twk->getRef<double> ("ChiSqrOverDOF", 4) - 73.9528 ), 1.5);
+ TS_ASSERT_LESS_THAN(fabs(Twk->getRef<double> ("ChiSqrOverDOF", 4) - 76.7824), 1.5);
- TS_ASSERT_LESS_THAN(fabs(Twk->getRef<double> ("TotIntensity", 0) - 4389.8), 10);
+ TS_ASSERT_LESS_THAN(fabs(Twk->getRef<double> ("TotIntensity", 0) - 4361.8 ), 10);
- /* std::vector<std::string> names = Twk->getColumnNames();
+ /* std::vector<std::string> names = Twk->getColumnNames();
std::cout<<"Intensitty="<<intensity<<" sigma="<<sigma<<
" Theoret intensity="<<TotIntensity<<std::endl;
@@ -248,6 +248,7 @@ class IntegratePeakTimeSlicesTest: public CxxTest::TestSuite
}
+
Intensitty=59870.7 sigma=665.098
Time 19250 19350 19450 19550 19650 19750 19850
Channel 12 13 14 15 16 17 18
@@ -324,31 +325,31 @@ ISAWIntensityError 108.136 189.481 271.573 353.892 271.573
End Col 38 38 38 38 38 38 38
==================================================================================
Circular Regions 4sigs, CalcVar=0(redo of optimization). Constraints on variables added
-Intensitty=59869.3 sigma=537.644 Theoret intensity=60000
+Intensitty=60340.5 sigma=466.885 Theoret intensity=60000
Act Int 3750 7500 11250 15000 11250 7500 3750
Time 19250 19350 19450 19550 19650 19750 19850
Channel 12 13 14 15 16 17 18
- Background 1.41788 1.43576 1.45363 1.47146 1.45363 1.43576 1.41788
- Intensity 3746.19 7492.38 11238.6 14984.8 11238.6 7492.38 3746.19
+ Background 1.41953 1.23611 1.23688 1.23589 1.23688 1.23611 1.41953
+ Intensity 3745.82 7568.01 11331 15097.3 11331 7568.01 3745.82
Mcol 27 27 27 27 27 27 27
- Mrow 21.9999 21.9999 21.9999 21.9999 21.9999 21.9999 21.9999
- SScol 4.44989 4.44989 4.4499 4.44995 4.4499 4.44989 4.44989
- SSrow 4.44998 4.44998 4.45 4.45004 4.45 4.44998 4.44998
- SSrc 0.000104969 0.000209915 0.000314833 0.000419728 0.000314833 0.000209915 0.000104969
- NCells 457 457 457 457 457 457 457
- ChiSqrOverDOF 9.60569 32.8682 73.9528 131.467 73.9528 32.8682 9.60569
- TotIntensity 4389.8 8139.8 11889.8 15639.8 11889.8 8139.8 4389.8
-BackgroundError 0.166822 0.308586 0.462878 0.615341 0.462878 0.308586 0.166822
-FitIntensityError 37.7145 69.764 104.646 137.837 104.646 69.764 37.7145
- ISAWIntensity 3741.83 7483.66 11225.5 14967.3 11225.5 7483.66 3741.83
-ISAWIntensityError 104.163 169.363 239.377 308.855 239.377 169.363 104.163
- TotalBoundary 151.2 151.2 151.2 151.2 151.2 151.2 151.2
- NBoundaryCells 108 108 108 108 108 108 108
- Start Row 10 10 10 10 10 10 10
- End Row 34 34 34 34 34 34 34
- Start Col 15 15 15 15 15 15 15
- End Col 39 39 39 39 39 39 39
-TotIntensityError 66.2556 90.2208 109.04 125.059 109.04 90.2208 66.2556
+ Mrow 22 22 22 22 22 22 22
+ SScol 4.44865 4.45006 4.45015 4.45 4.45015 4.45006 4.44865
+ SSrow 4.44851 4.45003 4.45003 4.45013 4.45003 4.45003 4.44851
+ SSrc 0.000101854 0.000169181 0.000257058 0.000347122 0.000257058 0.000169181 0.000101854
+ NCells 437 437 437 437 437 437 437
+ ChiSqrOverDOF 9.98124 34.0985 76.7824 136.568 76.7824 34.0985 9.98124
+ TotIntensity 4361.8 8111.8 11861.8 15611.8 11861.8 8111.8 4361.8
+BackgroundError 0.175194 0.317447 0.480952 0.242143 0.480952 0.317447 0.175194
+FitIntensityError 38.725 65.9045 103.054 124.116 103.054 65.9045 38.725
+ ISAWIntensity 3741.47 7571.62 11321.3 15071.7 11321.3 7571.62 3741.47
+ISAWIntensityError 104.132 167.022 237.858 165.375 237.858 167.022 104.132
+ TotalBoundary 123.2 123.2 123.2 123.2 123.2 123.2 123.2
+ NBoundaryCells 88 88 88 88 88 88 88
+ Start Row 11 11 11 11 11 11 11
+ End Row 33 33 33 33 33 33 33
+ Start Col 16 16 16 16 16 16 16
+ End Col 38 38 38 38 38 38 38
+TotIntensityError 66.0439 90.0655 108.912 124.947 108.912 90.0655 66.0439
*/
View
2  Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/BivariateNormal.h
@@ -144,7 +144,7 @@ namespace Mantid
return false;
}
-
+ bool CalcVxx, CalcVyy, CalcVxy;
protected:
void init();
View
49 Code/Mantid/Framework/CurveFitting/src/BivariateNormal.cpp
@@ -52,10 +52,11 @@ BivariateNormal::BivariateNormal() : API::IFunction1D(), expVals(NULL)
{
LastParams[IVXX] = -1;
- // g_log.setLevel(7);
-
+// g_log.setLevel(7);
+ CalcVxx= CalcVyy= CalcVxy= false;
//BackConstraint = IntensityConstraint= MeanxConstraint = MeanyConstraint = NULL;
expVals = NULL;
+ CalcVxx = CalcVxy = CalcVyy = false;
}
BivariateNormal::~BivariateNormal()
@@ -194,7 +195,6 @@ void BivariateNormal::functionDeriv1D(API::Jacobian *out, const double *xValues,
double C = -LastParams[IVYY] / 2 / uu;
-
double SIVXX = coefExp * expVals[x] * (C +
-LastParams[IVYY]/uu* (coefx2 * (c - LastParams[IXMEAN])* (c- LastParams[IXMEAN]) +
coefxy * (r - LastParams[IYMEAN]) * (c - LastParams[IXMEAN]) + coefy2
@@ -226,6 +226,7 @@ void BivariateNormal::functionDeriv1D(API::Jacobian *out, const double *xValues,
* (r - LastParams[IYMEAN]) * (r - LastParams[IYMEAN]))
+(r - LastParams[IYMEAN]) * (c - LastParams[IXMEAN])/uu) ;
+
if( !CalcVxx)
out->set(x, IVXX, SIVXX);
else
@@ -379,9 +380,6 @@ void BivariateNormal::initCommon()
}
- double mIx, mx, mIy, my;
-
- double SIxx, SIyy, SIxy, Sxx, Syy, Sxy;
mIx = Attrib[S_xint] / Attrib[S_int];
@@ -422,9 +420,7 @@ void BivariateNormal::initCommon()
if( getConstraint(1) == NULL)
- { //boost::shared_ptr<BoundaryConstraint> ic(new BoundaryConstraint(this, "Intensity", 0, maxIntensity));
- //IntensityConstraint =ic;// new BoundaryConstraint(this, "Intensity", 0, maxIntensity);
-
+ {
addConstraint(new BoundaryConstraint(this, "Intensity", 0, maxIntensity));
}
@@ -433,8 +429,7 @@ void BivariateNormal::initCommon()
double maxMeany = MinY * .1 + .9 * MaxY;
if( getConstraint(3) == NULL)
- { //boost::shared_ptr<BoundaryConstraint>my(new BoundaryConstraint(this, "Mrow", minMeany, maxMeany));
- //MeanyConstraint = my;// new BoundaryConstraint(this, "Mrow", minMeany, maxMeany);
+ {
addConstraint(new BoundaryConstraint(this, "Mrow", minMeany, maxMeany));
}
@@ -442,8 +437,7 @@ void BivariateNormal::initCommon()
double minMeanx = MinX * .9 + .1 * MaxX;
double maxMeanx = MinX * .1 + .9 * MaxX;
if( getConstraint(2)==NULL)
- {//boost::shared_ptr<BoundaryConstraint>mx(new BoundaryConstraint(this, "Mcol", minMeanx, maxMeanx));
- // MeanxConstraint =mx;// new BoundaryConstraint(this, "Mcol", minMeanx, maxMeanx);
+ {
addConstraint(new BoundaryConstraint(this, "Mcol", minMeanx, maxMeanx));
}
@@ -475,6 +469,7 @@ void BivariateNormal::initCommon()
if (getTie(IVYY) == NULL)
{
tie("SSrow", ssyy.str());
+ CalcVxx=true;
}
@@ -486,7 +481,7 @@ void BivariateNormal::initCommon()
if (getTie(IVXX) == NULL)
{
tie("SScol", ssxx.str());
-
+ CalcVyy=true;
}
ssxy << std::string("(") << (SIxy) << "+(Mcol-" << (mIx) << ")*(Mrow-" << (mIy) << ")*"
@@ -497,7 +492,7 @@ void BivariateNormal::initCommon()
if (getTie(IVXY) == NULL)
{
tie("SSrc", ssxy.str());
-
+ CalcVxy = true;
}
}
CommonsOK = true;
@@ -518,25 +513,22 @@ void BivariateNormal::initCommon()
if (!ParamsOK)
{
-
- for (size_t i = 0; i < nParams(); i++)
- LastParams[i] = getParameter(i);
-
-
+ for (size_t i = 0; i < nParams(); i++)
+ LastParams[i] = getParameter(i);
double Varxx = LastParams[IVXX];
if (CalcVxx)
{
- Varxx = (SIxx + (getParameter("Mcol") - mIx) * (getParameter("Mcol") - mIx) *
- TotI - getParameter("Background") * Sxx - getParameter("Background") * (getParameter("Mcol")
+ Varxx = (SIxx + (getParameter("Mcol") - mIx) * (getParameter("Mcol") - mIx) * TotI
+ - getParameter("Background") * Sxx - getParameter("Background") * (getParameter("Mcol")
- (mx)) * (getParameter("Mcol") - (mx)) * TotN) / (TotI - getParameter("Background")
* TotN);
- LastParams[IVXX] =Varxx;
- }
+ LastParams[IVXX] = Varxx;
- double Varyy = LastParams[IVYY];
+ }
+ double Varyy = LastParams[IVYY];
if (CalcVyy)
{
@@ -544,11 +536,10 @@ void BivariateNormal::initCommon()
- getParameter("Background") * (Syy) - getParameter("Background") * (getParameter("Mrow")
- (my)) * (getParameter("Mrow") - (my)) * TotN) / (TotI - getParameter("Background")
* TotN);
- LastParams[IVYY] =Varyy;
+ LastParams[IVYY] = Varyy;
}
double Varxy = LastParams[IVXY];
-
if (CalcVxy)
{
Varxy = ((SIxy) + (getParameter("Mcol") - (mIx)) * (getParameter("Mrow") - (mIy)) * TotI
@@ -556,10 +547,10 @@ void BivariateNormal::initCommon()
- (mx)) * (getParameter("Mrow") - (my)) * TotN) / (TotI - getParameter("Background")
* TotN);
- LastParams[IVXY] =Varxy;
+ LastParams[IVXY] = Varxy;
}
- }
+ }
if (!CommonsOK || !ParamsOK)
View
14 Code/Mantid/Framework/CurveFitting/test/BivariateNormalTest.h
@@ -115,8 +115,10 @@ class BivariateNormalTest: public CxxTest::TestSuite
else if(sgn2>0)sgn2=-sgn2;else {sgn1=-sgn1+1; sgn2=sgn1;}
xvals.access().push_back( x );
yvals.access().push_back( y );
- data.access().push_back( NormVal(background,intensity,Mcol,
- Mrow,Vx,Vy, Vxy, y, x));
+ double val = NormVal(background,intensity,Mcol,
+ Mrow,Vx,Vy, Vxy, y, x);
+ data.access().push_back( val );
+
}
@@ -147,7 +149,7 @@ class BivariateNormalTest: public CxxTest::TestSuite
NormalFit.setParameter("SSrc", 2.243584414, true);
- NormalFit.setAttributeValue("CalcVariances", 1);
+ NormalFit.setAttributeValue("CalcVariances",1);
std::vector<double> out(nCells);
NormalFit.function1D(out.data(), xx, nCells);
@@ -177,7 +179,7 @@ class BivariateNormalTest: public CxxTest::TestSuite
NormalFit.functionDeriv1D(Jac.get(), xx, nCells);
-/* for( int i=0; i< nCells; i+=6)//points
+ /* for( int i=0; i< nCells; i+=6)//points
{
std::cout<<i;
for( int j=0;j< 7;j++)//params
@@ -200,8 +202,8 @@ class BivariateNormalTest: public CxxTest::TestSuite
}
}
-
- /* std::vector<double>out1(nCells);
+/*
+ std::vector<double>out1(nCells);
std::vector<double>out2(nCells);
for( size_t param = 0; param < 4; param++ )
{
Please sign in to comment.
Something went wrong with that request. Please try again.