From 67ffac194801dea96754b6d96332a57fa85259fa Mon Sep 17 00:00:00 2001 From: ageay Date: Tue, 20 Aug 2013 08:42:37 +0000 Subject: [PATCH] Externalize Kriging matrix for advanced users. --- src/INTERP_KERNEL/InterpKernelMatrixTools.cxx | 7 + .../MEDCouplingFieldDiscretization.cxx | 122 ++++++++++++------ .../MEDCouplingFieldDiscretization.hxx | 2 + src/MEDCoupling_Swig/MEDCouplingBasicsTest.py | 16 ++- src/MEDCoupling_Swig/MEDCouplingCommon.i | 30 ++++- 5 files changed, 132 insertions(+), 45 deletions(-) diff --git a/src/INTERP_KERNEL/InterpKernelMatrixTools.cxx b/src/INTERP_KERNEL/InterpKernelMatrixTools.cxx index ada00883d..7b228af4f 100644 --- a/src/INTERP_KERNEL/InterpKernelMatrixTools.cxx +++ b/src/INTERP_KERNEL/InterpKernelMatrixTools.cxx @@ -19,8 +19,10 @@ // Author : Anthony Geay (CEA/DEN) #include "InterpKernelMatrixTools.hxx" +#include "InterpKernelException.hxx" #include "InterpKernelAutoPtr.hxx" +#include #include namespace INTERP_KERNEL @@ -396,6 +398,11 @@ 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 coords=getLocalizationOfDiscValues(mesh); - int nbOfPts=coords->getNumberOfTuples(); - int dimension=coords->getNumberOfComponents(); - // - int delta=0; - MEDCouplingAutoRefCountObjectPtr 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 m(computeEvaluationMatrixOnGivenPts(mesh,loc,nbOfTargetPoints,nbCols)); + MEDCouplingAutoRefCountObjectPtr 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 matrixInv(computeInverseMatrix(mesh,isDrift,nbRows)); // + MEDCouplingAutoRefCountObjectPtr coords=getLocalizationOfDiscValues(mesh); + int nbOfPts(coords->getNumberOfTuples()),dimension(coords->getNumberOfComponents()); MEDCouplingAutoRefCountObjectPtr locArr=DataArrayDouble::New(); locArr->useArray(loc,false,CPP_DEALLOC,nbOfTargetPoints,dimension); + nbCols=nbOfPts; + // MEDCouplingAutoRefCountObjectPtr matrix2=coords->buildEuclidianDistanceDenseMatrixWith(locArr); - operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfPts*nbOfTargetPoints,matrix2->getPointer()); + operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfTargetPoints*nbOfPts,matrix2->getPointer()); + // MEDCouplingAutoRefCountObjectPtr 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;ibegin()),*workCst2(loc); + for(int i=0;igetNumberOfComponents(); - MEDCouplingAutoRefCountObjectPtr 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 ret(DataArrayDouble::New()); + ret->alloc(nbOfTargetPoints,nbRows); + INTERP_KERNEL::matrixProduct(matrix3->begin(),nbOfTargetPoints,nbRows,matrixInv->begin(),nbRows,nbRows,ret->getPointer()); + MEDCouplingAutoRefCountObjectPtr ret2(DataArrayDouble::New()); + ret2->alloc(nbOfTargetPoints*nbOfPts,1); + workCst=ret->begin(); work=ret2->getPointer(); + for(int i=0;i coords=getLocalizationOfDiscValues(mesh); + int nbOfPts=coords->getNumberOfTuples(); + MEDCouplingAutoRefCountObjectPtr matrix=coords->buildEuclidianDistanceDenseMatrix(); + operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfPts*nbOfPts,matrix->getPointer()); + // Drift + MEDCouplingAutoRefCountObjectPtr matrixWithDrift=performDrift(matrix,coords,isDrift); + MEDCouplingAutoRefCountObjectPtr matrixInv=DataArrayDouble::New(); + matSz=nbOfPts+isDrift; + matrixInv->alloc(matSz*matSz,1); + INTERP_KERNEL::inverseMatrix(matrixWithDrift->getConstPointer(),matSz,matrixInv->getPointer()); + return matrixInv.retn(); } /*! @@ -2624,32 +2673,21 @@ void MEDCouplingFieldDiscretizationKriging::reprQuickOverview(std::ostream& stre * * \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 coords=getLocalizationOfDiscValues(mesh); - int nbOfPts=coords->getNumberOfTuples(); - //int dimension=coords->getNumberOfComponents(); - MEDCouplingAutoRefCountObjectPtr matrix=coords->buildEuclidianDistanceDenseMatrix(); - operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfPts*nbOfPts,matrix->getPointer()); - // Drift - MEDCouplingAutoRefCountObjectPtr matrixWithDrift=performDrift(matrix,coords,isDrift); - MEDCouplingAutoRefCountObjectPtr matrixInv=DataArrayDouble::New(); - matrixInv->alloc((nbOfPts+isDrift)*(nbOfPts+isDrift),1); - INTERP_KERNEL::inverseMatrix(matrixWithDrift->getConstPointer(),nbOfPts+isDrift,matrixInv->getPointer()); - // + int nbRows(-1); + MEDCouplingAutoRefCountObjectPtr matrixInv(computeInverseMatrix(mesh,isDrift,nbRows)); MEDCouplingAutoRefCountObjectPtr KnewiK=DataArrayDouble::New(); - KnewiK->alloc((nbOfPts+isDrift)*1,1); + KnewiK->alloc(nbRows*1,1); MEDCouplingAutoRefCountObjectPtr 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(); } diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx index 84b39ca47..2cbd5511f 100644 --- a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx @@ -383,6 +383,8 @@ namespace ParaMEDMEM 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; diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py index 739813d1c..61ebab008 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py @@ -13752,7 +13752,6 @@ class MEDCouplingBasicsTest(unittest.TestCase): 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) @@ -13783,9 +13782,22 @@ class MEDCouplingBasicsTest(unittest.TestCase): 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 diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index 3c3e6a0be..3da4c59fb 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -3243,7 +3243,7 @@ namespace ParaMEDMEM %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); @@ -3252,6 +3252,34 @@ namespace ParaMEDMEM 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 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 -- 2.39.2