// Author : Anthony Geay (CEA/DEN)
#include "InterpKernelMatrixTools.hxx"
+#include "InterpKernelException.hxx"
#include "InterpKernelAutoPtr.hxx"
+#include <sstream>
#include <algorithm>
namespace INTERP_KERNEL
void matrixProduct(const double *A, int n1, int p1, const double *B, int n2, int p2, double *C)
{
+ if(p1!=n2)
+ {
+ std::ostringstream oss; oss << "matrixProduct : the size of input matrix are not coherent the nb of cols of input matrix #0 is " << p1 << " whereas the number of rows of input matrix #1 is " << n2 << " !";
+ throw INTERP_KERNEL::Exception(oss.str().c_str());
+ }
for(int i=0;i<n1;i++)
{
for(int j=0;j<p2;j++)
DataArrayDouble *MEDCouplingFieldDiscretizationKriging::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfTargetPoints) const
{
- if(!mesh)
- throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::getValueOnMulti : NULL input mesh !");
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords=getLocalizationOfDiscValues(mesh);
- int nbOfPts=coords->getNumberOfTuples();
- int dimension=coords->getNumberOfComponents();
- //
- int delta=0;
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> KnewiK=computeVectorOfCoefficients(mesh,arr,delta);
+ if(!arr || !arr->isAllocated())
+ throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::getValueOnMulti : input array is null or not allocated !");
+ int nbOfRows(getNumberOfMeshPlaces(mesh));
+ if(arr->getNumberOfTuples()!=nbOfRows)
+ {
+ std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationKriging::getValueOnMulti : input array does not have correct number of tuples ! Excepted " << nbOfRows << " having " << arr->getNumberOfTuples() << " !";
+ throw INTERP_KERNEL::Exception(oss.str().c_str());
+ }
+ int nbCols(-1),nbCompo(arr->getNumberOfComponents());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> m(computeEvaluationMatrixOnGivenPts(mesh,loc,nbOfTargetPoints,nbCols));
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New());
+ ret->alloc(nbOfTargetPoints,nbCompo);
+ INTERP_KERNEL::matrixProduct(m->begin(),nbOfTargetPoints,nbCols,arr->begin(),nbOfRows,nbCompo,ret->getPointer());
+ return ret.retn();
+}
+
+void MEDCouplingFieldDiscretizationKriging::reprQuickOverview(std::ostream& stream) const throw(INTERP_KERNEL::Exception)
+{
+ stream << "Kriging spatial discretization.";
+}
+
+/*!
+ * Returns the matrix of size nbRows = \a nbOfTargetPoints and \a nbCols = \a nbCols. This matrix is useful if
+ *
+ * \return the new result matrix to be deallocated by the caller.
+ */
+DataArrayDouble *MEDCouplingFieldDiscretizationKriging::computeEvaluationMatrixOnGivenPts(const MEDCouplingMesh *mesh, const double *loc, int nbOfTargetPoints, int& nbCols) const
+{
+ int isDrift(-1),nbRows(-1);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixInv(computeInverseMatrix(mesh,isDrift,nbRows));
//
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords=getLocalizationOfDiscValues(mesh);
+ int nbOfPts(coords->getNumberOfTuples()),dimension(coords->getNumberOfComponents());
MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> locArr=DataArrayDouble::New();
locArr->useArray(loc,false,CPP_DEALLOC,nbOfTargetPoints,dimension);
+ nbCols=nbOfPts;
+ //
MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrix2=coords->buildEuclidianDistanceDenseMatrixWith(locArr);
- operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfPts*nbOfTargetPoints,matrix2->getPointer());
+ operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfTargetPoints*nbOfPts,matrix2->getPointer());
+ //
MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrix3=DataArrayDouble::New();
- matrix3->alloc((nbOfPts+delta)*nbOfTargetPoints,1);
+ matrix3->alloc(nbOfTargetPoints*nbRows,1);
double *work=matrix3->getPointer();
- const double *workCst=matrix2->getConstPointer();
- const double *workCst2=loc;
- for(int i=0;i<nbOfTargetPoints;i++,workCst+=nbOfPts,workCst2+=delta-1)
+ const double *workCst(matrix2->begin()),*workCst2(loc);
+ for(int i=0;i<nbOfTargetPoints;i++,workCst+=nbOfPts,workCst2+=isDrift-1)
{
for(int j=0;j<nbOfPts;j++)
- work[j*nbOfTargetPoints+i]=workCst[j];
- work[nbOfPts*nbOfTargetPoints+i]=1.0;
- for(int j=0;j<delta-1;j++)
- work[(nbOfPts+1+j)*nbOfTargetPoints+i]=workCst2[j];
+ work[i*nbRows+j]=workCst[j];
+ work[i*nbRows+nbOfPts]=1.0;
+ for(int j=0;j<isDrift-1;j++)
+ work[i*nbRows+(nbOfPts+1+j)]=workCst2[j];
}
- //
- int nbOfCompo=arr->getNumberOfComponents();
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
- ret->alloc(nbOfTargetPoints,nbOfCompo);
- INTERP_KERNEL::matrixProduct(KnewiK->getConstPointer(),1,nbOfPts+delta,matrix3->getConstPointer(),nbOfPts+delta,nbOfTargetPoints*nbOfCompo,ret->getPointer());
- return ret.retn();
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New());
+ ret->alloc(nbOfTargetPoints,nbRows);
+ INTERP_KERNEL::matrixProduct(matrix3->begin(),nbOfTargetPoints,nbRows,matrixInv->begin(),nbRows,nbRows,ret->getPointer());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret2(DataArrayDouble::New());
+ ret2->alloc(nbOfTargetPoints*nbOfPts,1);
+ workCst=ret->begin(); work=ret2->getPointer();
+ for(int i=0;i<nbOfTargetPoints;i++,workCst+=nbRows)
+ work=std::copy(workCst,workCst+nbOfPts,work);
+ return ret2.retn();
}
-void MEDCouplingFieldDiscretizationKriging::reprQuickOverview(std::ostream& stream) const throw(INTERP_KERNEL::Exception)
+/*!
+ * This method returns the square matrix of size \a matSz that is the inverse of the kriging matrix. The returned matrix can returned all the coeffs of kriging
+ * when multiplied by the vector of values attached to each point.
+ *
+ * \param [out] isDrift return if drift coefficients are present in the returned vector of coefficients. If different from 0 there is presence of drift coefficients.
+ * \param [out] matSz the size of returned square matrix
+ * \return the new result matrix to be deallocated by the caller.
+ */
+DataArrayDouble *MEDCouplingFieldDiscretizationKriging::computeInverseMatrix(const MEDCouplingMesh *mesh, int& isDrift, int& matSz) const
{
- stream << "Kriging spatial discretization.";
+ if(!mesh)
+ throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::computeVectorOfCoefficients : NULL input mesh !");
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords=getLocalizationOfDiscValues(mesh);
+ int nbOfPts=coords->getNumberOfTuples();
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrix=coords->buildEuclidianDistanceDenseMatrix();
+ operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfPts*nbOfPts,matrix->getPointer());
+ // Drift
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixWithDrift=performDrift(matrix,coords,isDrift);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixInv=DataArrayDouble::New();
+ matSz=nbOfPts+isDrift;
+ matrixInv->alloc(matSz*matSz,1);
+ INTERP_KERNEL::inverseMatrix(matrixWithDrift->getConstPointer(),matSz,matrixInv->getPointer());
+ return matrixInv.retn();
}
/*!
*
* \param [in] mesh is the sources of nodes on which kriging will be done regarding the parameters and the value of \c this->getSpaceDimension()
* \param [in] arr input field DataArrayDouble whose number of tuples must be equal to the number of nodes in \a mesh
- * \param [out] isDrift return if drift coefficients are present in the returned vector of coefficients, and if. If different from 0 there is presence of drift coefficients.
+ * \param [out] isDrift return if drift coefficients are present in the returned vector of coefficients. If different from 0 there is presence of drift coefficients.
* Whatever the value of \a isDrift the number of tuples of returned DataArrayDouble will be equal to \c arr->getNumberOfTuples() + \a isDrift.
* \return a newly allocated array containing coefficients including or not drift coefficient at the end depending the value of \a isDrift parameter.
*/
DataArrayDouble *MEDCouplingFieldDiscretizationKriging::computeVectorOfCoefficients(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, int& isDrift) const
{
- if(!mesh)
- throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::computeVectorOfCoefficients : NULL input mesh !");
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords=getLocalizationOfDiscValues(mesh);
- int nbOfPts=coords->getNumberOfTuples();
- //int dimension=coords->getNumberOfComponents();
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrix=coords->buildEuclidianDistanceDenseMatrix();
- operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfPts*nbOfPts,matrix->getPointer());
- // Drift
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixWithDrift=performDrift(matrix,coords,isDrift);
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixInv=DataArrayDouble::New();
- matrixInv->alloc((nbOfPts+isDrift)*(nbOfPts+isDrift),1);
- INTERP_KERNEL::inverseMatrix(matrixWithDrift->getConstPointer(),nbOfPts+isDrift,matrixInv->getPointer());
- //
+ int nbRows(-1);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixInv(computeInverseMatrix(mesh,isDrift,nbRows));
MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> KnewiK=DataArrayDouble::New();
- KnewiK->alloc((nbOfPts+isDrift)*1,1);
+ KnewiK->alloc(nbRows*1,1);
MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr2=DataArrayDouble::New();
- arr2->alloc((nbOfPts+isDrift)*1,1);
+ arr2->alloc(nbRows*1,1);
double *work=std::copy(arr->begin(),arr->end(),arr2->getPointer());
std::fill(work,work+isDrift,0.);
- INTERP_KERNEL::matrixProduct(matrixInv->getConstPointer(),nbOfPts+isDrift,nbOfPts+isDrift,arr2->getConstPointer(),nbOfPts+isDrift,1,KnewiK->getPointer());
+ INTERP_KERNEL::matrixProduct(matrixInv->getConstPointer(),nbRows,nbRows,arr2->getConstPointer(),nbRows,1,KnewiK->getPointer());
return KnewiK.retn();
}
DataArrayDouble *getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfPoints) const;
void reprQuickOverview(std::ostream& stream) const throw(INTERP_KERNEL::Exception);
public://specific part
+ DataArrayDouble *computeEvaluationMatrixOnGivenPts(const MEDCouplingMesh *mesh, const double *loc, int nbOfTargetPoints, int& nbCols) const;
+ DataArrayDouble *computeInverseMatrix(const MEDCouplingMesh *mesh, int& isDrift, int& matSz) const;
DataArrayDouble *computeVectorOfCoefficients(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, int& isDrift) const;
protected:
void operateOnDenseMatrix(int spaceDimension, int nbOfElems, double *matrixPtr) const;
self.assertTrue(bbt.computeNbOfInteractionsWith(bbs,1e-12).isEqual(wExp))
pass
-
def testKrSpatialDiscretization2(self):
srcPointCoordsXY=DataArrayDouble([0.8401877171547095,0.39438292681909304,0.7830992237586059,0.7984400334760733,0.9116473579367843,0.19755136929338396,0.335222755714889,0.768229594811904,0.2777747108031878,0.5539699557954305,0.47739705186216025,0.6288709247619244,0.36478447279184334,0.5134009101956155,0.9522297251747128,0.9161950680037007,0.6357117279599009,0.7172969294326831,0.14160255535580338,0.6069688762570586,0.01630057162432958,0.24288677062973696,0.13723157678601872,0.8041767542269904,0.15667908925408455,0.4009443942461835,0.12979044678145574,0.10880880202576929,0.998924518003559,0.21825690531090688,0.5129323944043984,0.8391122346926072,0.6126398325956612,0.29603161769734304,0.6375522677030192,0.5242871900667843,0.493582986990727,0.9727750238835695,0.29251678441302703,0.7713576977939148,0.5267449792133388,0.7699138362751873,0.4002286220901779,0.8915294520051822,0.2833147460051415,0.3524583472648907,0.8077245200088827,0.9190264739650424,0.06975527623191256,0.9493270753646861,0.5259953502221011,0.08605584785624214,0.19221384599442307,0.6632269270081198,0.8902326025488938,0.3488929352485076,0.06417132078864207,0.02002304886468828,0.4577017372742769,0.06309583832653977,0.23827995417559517,0.9706341316786754,0.9022080734848082,0.8509197867712563,0.2666657493760184,0.5397603407221662,0.3752069763723793,0.7602487363667454,0.5125353641400744,0.6677237607854063,0.5316064341606602,0.039280343353413204,0.4376375965949323,0.9318350562508382,0.9308097953585953,0.7209523430657351,0.28429340305006756,0.7385343149018168,0.6399788165651163,0.3540486797476414,0.687861390266503,0.16597416632155615,0.4401045276038835,0.880075236260926,0.829201093329676,0.3303371296871161,0.22896817104377232,0.8933724145839793,0.35036017855180435,0.6866699083180492,0.9564682529105192,0.5886401331930609,0.6573040395310633,0.8586763259296661,0.4395599194986559,0.9239697889070817,0.39843666665183225,0.8147668963366965,0.6842185252738271,0.9109720307919067,0.4824906566564416,0.21582495896882609,0.9502523741453198,0.9201282537170352,0.14766001475400292,0.8810621695039152,0.641080596317109,0.43195341826973177,0.6195964839400707,0.281059412416564,0.7860020980173732,0.3074578737409124,0.44703357920378145,0.22610662515559543,0.18753310953617705,0.27623467206779617,0.5564437553083728,0.4165012805799494,0.16960708618611428,0.9068039338601771,0.10317118843233734,0.1260753390966334,0.49544406658757667,0.7604752284290619,0.9847516650262995,0.9350039865518939,0.6844450168704823,0.3831883312124705,0.7497708824229291,0.36866354167864823,0.2941603620043771,0.2322615386137094,0.5844885006474743,0.24441273568403568,0.15238979186508328,0.7321485158671385,0.12547490472228962,0.7934703881821923,0.164101933671209,0.7450713891280216,0.07452980059875632,0.9501040316885822,0.05252926240327268,0.5215633798025378,0.1762106563785163,0.24006237240511102,0.797798051870334,0.732654411686889,0.6565636529850605,0.9674051385221095,0.6394583455470663,0.7597348418830591,0.09348047715308166,0.13490241166898162,0.5202100698464597,0.07823214171371988,0.06990639775521419,0.2046550862512808,0.4614204733918516,0.8196772801781433,0.5733186283955903,0.7555808353962288,0.05193881879185271,0.1578071285774033,0.9999935710802644,0.204328610656936,0.8899556444445419,0.12546847580255405,0.9977989993047895,0.054057577650089554,0.8705398649305757,0.07232879943788462,0.004161608873010431,0.9230691273338484,0.5938921792404224,0.180372265717188,0.16313149927329806,0.3916902306450951,0.9130266774040771,0.8196951527240198,0.35909536870154335,0.552485022485482,0.5794299941414176,0.452575845854625,0.687387434620125,0.09964006352221597,0.5308079880340062,0.7572938323753392,0.30429514977349675,0.9922284614258579,0.5769711125534824,0.877613778169087,0.7478092963564253,0.6289099313453351,0.03542090674649035,0.7478028669710285,0.8332385420022712,0.9253765511910322,0.8732713427735824,0.8310375408413995],100,2)
srcFieldValsOnPoints=DataArrayDouble([0.7643742528498438,-0.023507696856211995,1.1082895131907775,0.6299357452572031,0.8892623544912389,0.72212114810697,0.9196401044320336,-0.759961711221917,0.40801932617748826,0.8441134300809151,0.982483804252809,0.6752368914020778,0.9924403977479798,1.1063334970204484,0.9403055261137516,0.3624481886322733,1.1344772505996308,0.7522965618948239,0.17077741651388564,0.6504551671311436,0.45843479588425423,0.41098905950326753,1.0681420394050904,-0.3483587903820091,0.5620151050607809,1.384969776596035,0.7948875141132845,0.7931192000237167,1.062498042490183,1.3709072529577366,0.44929346605311893,-0.4469683401788374,0.9035857424514101,0.6137249300593463,0.6355610879026966,1.4318174829507697,0.3097567072129551,-0.20515052260807165,0.6922559820922779,1.0341638749443423,1.3072652153341024,0.38511367353000436,0.9160514929274943,0.54513408530581,0.722252267913328,0.06684522818576251,0.10571899758067793,0.3193844999960903,0.5213532270828706,-0.04834998649603944,1.2408805068350615,-0.7632951295676795,0.5980054665011202,0.9064738717547436,1.1541070755096696,1.008234260272265,1.2225806960553827,1.0788560195121106,0.9818990282104452,0.5621951325841853,1.0796757508374188,0.5082872315589883,-0.9153702001062469,0.9560418838920791,0.9251098559152824,1.1603063610984021,1.2122303611181837,0.7379539363312343,0.6877611899207183,0.723966552446608,0.5596025827162566,0.8849725005989729,1.0908363665075547,0.08956512916455672,-0.10247645571248344,0.3236718069555875,1.069478546398975,1.3900071080692746,1.0322398863403262,0.45315515354558034,0.4249870238786733,1.030226761858634,0.974024629584669,1.2838885424020365,1.3451943506525155,1.4029933267831995,0.6025539675442462,1.2947650597767038,1.0006061239483002,-0.4017336259949164,0.8771165113201297,0.9158909024218246,1.403798605551443,0.4742904006425974,0.3671787905896653,0.20646491720419674,0.40739337434288925,0.7341932402033597,-0.4295893651836911,-0.3187777570661546],100,1)
coeffs,isDrift=fd.computeVectorOfCoefficients(mesh,srcFieldValsOnPoints)
self.assertEqual(3,isDrift)
self.assertTrue(coeffsExpected.isEqual(coeffs,1e-8))
+ # testing matrix
+ pts3=[-0.5,-0.5,-0.5,-0.35,-0.35,-0.2]
+ mesh.setCoords(srcPointCoordsXY[:4])
+ m,nbCols=fd.computeEvaluationMatrixOnGivenPts(mesh,pts3)
+ self.assertTrue(m.isEqual(DataArrayDouble([0.05768877688524917,-4.438982030395039,1.9495386255911573,3.431754627918642,0.11803848510231275,-4.138339658420563,1.6630742187104417,3.357226954607818,0.14630203028580618,-3.5156045565871734,1.414680070737206,2.954622455564169]),1e-12))
+ if MEDCouplingHasNumPyBindings():
+ import numpy as np
+ m0=m.toNumPyArray() ; m0=m0.reshape(3,nbCols) ; m0=np.matrix(m0)
+ srcFieldValsOnPoints2=DataArrayDouble(4,2) ; srcFieldValsOnPoints2[:,0]=srcFieldValsOnPoints[:4] ; srcFieldValsOnPoints2[:,1]=2*srcFieldValsOnPoints[:4]
+ n0=srcFieldValsOnPoints2.toNumPyArray() ; n0=n0.reshape(4,2) ; n0=np.matrix(n0)
+ #
+ f=MEDCouplingFieldDouble.New(ON_NODES_KR,ONE_TIME) ; f.setMesh(mesh) ; f.setArray(srcFieldValsOnPoints2) ; f.checkCoherency()
+ self.assertTrue(DataArrayDouble(np.array((m0*n0))).isEqual(f.getValueOnMulti(pts3),1e-14))
+ pass
#
pass
- pass
def setUp(self):
pass
%extend ParaMEDMEM::MEDCouplingFieldDiscretizationKriging
{
- PyObject *computeVectorOfCoefficients(const MEDCouplingMesh *mesh, const DataArrayDouble *arr) const
+ PyObject *computeVectorOfCoefficients(const MEDCouplingMesh *mesh, const DataArrayDouble *arr) const throw(INTERP_KERNEL::Exception)
{
int ret1;
DataArrayDouble *ret0=self->computeVectorOfCoefficients(mesh,arr,ret1);
PyTuple_SetItem(ret,1,PyInt_FromLong(ret1));
return ret;
}
+
+ PyObject *computeInverseMatrix(const MEDCouplingMesh *mesh) const throw(INTERP_KERNEL::Exception)
+ {
+ int ret1(-1),ret2(-1);
+ DataArrayDouble *ret0=self->computeInverseMatrix(mesh,ret1,ret2);
+ PyObject *ret=PyTuple_New(3);
+ PyTuple_SetItem(ret,0,SWIG_NewPointerObj(SWIG_as_voidptr(ret0),SWIGTYPE_p_ParaMEDMEM__DataArrayDouble, SWIG_POINTER_OWN | 0 ));
+ PyTuple_SetItem(ret,1,PyInt_FromLong(ret1));
+ PyTuple_SetItem(ret,2,PyInt_FromLong(ret2));
+ return ret;
+ }
+
+ PyObject *computeEvaluationMatrixOnGivenPts(const MEDCouplingMesh *mesh, PyObject *locs) const throw(INTERP_KERNEL::Exception)
+ {
+ if(!mesh)
+ throw INTERP_KERNEL::Exception("wrap of MEDCouplingFieldDiscretizationKriging::computeEvaluationMatrixOnGivenPts : input mesh is empty !");
+ int sw,nbPts;
+ double v0; ParaMEDMEM::DataArrayDouble *v1(0); ParaMEDMEM::DataArrayDoubleTuple *v2(0); std::vector<double> v3;
+ const double *inp=convertObjToPossibleCpp5_Safe2(locs,sw,v0,v1,v2,v3,"wrap of MEDCouplingFieldDiscretizationKriging::computeEvaluationMatrixOnGivenPts",
+ mesh->getSpaceDimension(),true,nbPts);
+ //
+ int ret1(-1);
+ DataArrayDouble *ret0=self->computeEvaluationMatrixOnGivenPts(mesh,inp,nbPts,ret1);
+ PyObject *ret=PyTuple_New(2);
+ PyTuple_SetItem(ret,0,SWIG_NewPointerObj(SWIG_as_voidptr(ret0),SWIGTYPE_p_ParaMEDMEM__DataArrayDouble, SWIG_POINTER_OWN | 0 ));
+ PyTuple_SetItem(ret,1,PyInt_FromLong(ret1));
+ return ret;
+ }
}
namespace ParaMEDMEM
