From 3578394e40ae955a2f463e081e1368b08aa786f1 Mon Sep 17 00:00:00 2001 From: gracebato Date: Wed, 18 Sep 2019 21:25:45 -0700 Subject: [PATCH] CDM-only related PR --- models/cdm/cdm/CDM.py | 24 +++-- models/cdm/cdm/Fast.py | 31 +++--- models/cdm/cdm/Native.py | 13 ++- models/cdm/cdm/Source.py | 11 +- models/cdm/cdm/libcdm.py | 8 +- models/cdm/ext/cdm/source.cc | 11 +- models/cdm/lib/libcdm/Source.cc | 52 ++++++++-- models/cdm/lib/libcdm/Source.h | 8 +- models/cdm/lib/libcdm/Source.icc | 16 +-- models/cdm/lib/libcdm/cdm.cc | 169 +++++++++++++++++++++++-------- models/cdm/lib/libcdm/cdm.h | 15 ++- models/cdm/tests/libcdm.py | 3 +- 12 files changed, 259 insertions(+), 102 deletions(-) diff --git a/models/cdm/cdm/CDM.py b/models/cdm/cdm/CDM.py index ccc88747..d8fad827 100644 --- a/models/cdm/cdm/CDM.py +++ b/models/cdm/cdm/CDM.py @@ -16,6 +16,7 @@ # the encapsulation of the layout of the data in a file from .Data import Data as datasheet +import gc # declaration class CDM(altar.models.bayesian, family="altar.models.cdm"): @@ -62,8 +63,7 @@ class CDM(altar.models.bayesian, family="altar.models.cdm"): # public data parameters = 0 # adjusted during model initialization strategy = None # the strategy for computing the data log likelihood - - + # protocol obligations @altar.export def initialize(self, application): @@ -106,7 +106,7 @@ def initialize(self, application): theta = record.theta phi = record.phi # form the projection vectors and store them - self.los[obs, 0] = sin(theta) * cos(phi) + self.los[obs, 0] = -sin(theta) * cos(phi) self.los[obs, 1] = sin(theta) * sin(phi) self.los[obs, 2] = cos(theta) @@ -237,12 +237,18 @@ def initializeParameterSets(self): self.yIdx = self.xIdx + 1 self.dIdx = psets["depth"].offset self.openingIdx = psets["opening"].offset - self.aXIdx = psets["a"].offset - self.aYIdx = self.aXIdx + 1 - self.aZIdx = self.aXIdx + 2 - self.omegaXIdx = psets["omega"].offset - self.omegaYIdx = self.omegaXIdx + 1 - self.omegaZIdx = self.omegaXIdx + 2 +# self.aXIdx = psets["a"].offset +# self.aYIdx = self.aXIdx + 1 +# self.aZIdx = self.aXIdx + 2 + self.aXIdx = psets["aX"].offset + self.aYIdx = psets["aY"].offset + self.aZIdx = psets["aZ"].offset +# self.omegaXIdx = psets["omega"].offset +# self.omegaYIdx = self.omegaXIdx + 1 +# self.omegaZIdx = self.omegaXIdx + 2 + self.omegaXIdx = psets["omegaX"].offset + self.omegaYIdx = psets["omegaY"].offset + self.omegaZIdx = psets["omegaZ"].offset self.offsetIdx = psets["offsets"].offset # all done diff --git a/models/cdm/cdm/Fast.py b/models/cdm/cdm/Fast.py index 83bc63c4..1ccfad57 100644 --- a/models/cdm/cdm/Fast.py +++ b/models/cdm/cdm/Fast.py @@ -13,7 +13,7 @@ import altar # the pure python implementation of the CDM source from altar.models.cdm.ext import libcdm - +import gc # declaration class Fast: @@ -48,11 +48,13 @@ def initialize(self, model, **kwds): # attach the map of observations to their set libcdm.oid(source, oid) # inform the source about the parameter layout; assumes contiguous parameter sets +# libcdm.layout(source, +# model.xIdx, model.dIdx, model.openingIdx, model.aXIdx, model.omegaXIdx, +# model.offsetIdx) libcdm.layout(source, - model.xIdx, model.dIdx, - model.openingIdx, model.aXIdx, model.omegaXIdx, - model.offsetIdx) - + model.xIdx, model.dIdx, model.openingIdx, model.aXIdx, model.aYIdx, model.aZIdx, + model.omegaXIdx, model.omegaYIdx, model.omegaZIdx, model.offsetIdx) + # nothing to do return self @@ -89,19 +91,24 @@ def dataLikelihood(self, model, step): for sample in range(samples): # get the residuals residuals = predicted.getRow(sample) - # compute the norm - nrm = norm.eval(v=residuals, sigma_inv=cd_inv) - # and normalize it - llk = normalization - nrm**2 / 2 + + # compute the norm, and normalize it + normeval = norm.eval(v=residuals, sigma_inv=cd_inv) + llk = normalization - normeval**2.0 / 2.0 + # store it dataLLK[sample] = llk - + llk = None # + residuals = None # + normeval = None # + + cd_inv = None # + normalization = None # # all done return self - # private data source = None - + print('Garbage collected (Fast.py--bottom): ', gc.collect()) # end of file diff --git a/models/cdm/cdm/Native.py b/models/cdm/cdm/Native.py index 6302fd7a..5f9d6ba5 100644 --- a/models/cdm/cdm/Native.py +++ b/models/cdm/cdm/Native.py @@ -50,7 +50,7 @@ def dataLikelihood(self, model, step): samples = θ.rows # get the parameter sets psets = model.psets - + # get the offsets of the various parameter sets xIdx = model.xIdx yIdx = model.yIdx @@ -63,7 +63,7 @@ def dataLikelihood(self, model, step): omegaYIdx = model.omegaYIdx omegaZIdx = model.omegaZIdx offsetIdx = model.offsetIdx - + # get the observations los = model.los oid = model.oid @@ -90,10 +90,12 @@ def dataLikelihood(self, model, step): omegaY = parameters[omegaYIdx] omegaZ = parameters[omegaZIdx] + #print(x,' ',y,' ',d,' ',opening,' ',aX,' ',aY,' ',aZ,' ',omegaX,' ',omegaY,' ',omegaZ) # make a source using the sample parameters cdm = source(x=x, y=y, d=d, opening=opening, ax=aX, ay=aY, az=aZ, omegaX=omegaX, omegaY=omegaY, omegaZ=omegaZ, v=model.nu) + # compute the expected displacement u = cdm.displacements(locations=locations, los=los) @@ -103,11 +105,12 @@ def dataLikelihood(self, model, step): for obs in range(observations): # appropriate for the corresponding dataset u[obs] -= parameters[offsetIdx + oid[obs]] - + # compute the norm of the displacements - nrm = norm.eval(v=u, sigma_inv=cd_inv) + normeval = norm.eval(v=u, sigma_inv=cd_inv) # normalize and store it as the data log likelihood - dataLLK[sample] = normalization - nrm**2 /2 + dataLLK[sample] = normalization - normeval**2.0 /2 + #print(dataLLK[sample]) # all done return self diff --git a/models/cdm/cdm/Source.py b/models/cdm/cdm/Source.py index 448ac774..5906113c 100644 --- a/models/cdm/cdm/Source.py +++ b/models/cdm/cdm/Source.py @@ -78,10 +78,11 @@ def displacements(self, locations, los): # allocate space for the result u = altar.vector(shape=len(locations)) # compute the displacements - ue, un, uv = CDM(X=Xf, Y=Yf, X0=x_src, Y0=y_src, depth=d_src, + ue, un, uv = CDM(X=Xf, Y=Yf, X0=x_src, Y0=y_src, depth=d_src, opening=opening, ax=ax_src, ay=ay_src, az=az_src, omegaX=omegaX_src, omegaY=omegaY_src, omegaZ=omegaZ_src, - opening=opening, nu=v) + nu=v) + # go through each observation location for idx, (ux,uy,uz) in enumerate(zip(ue, un, uv)): # project the expected displacement along LOS and store @@ -92,9 +93,9 @@ def displacements(self, locations, los): # meta-methods - def __init__(self, x=x, y=y, d=d, - ax=ax, ay=ay, az=az, omegaX=omegaX, omegaY=omegaY, omegaZ=omegaZ, - opening=opening, v=v, **kwds): + def __init__(self, x=x, y=y, d=d, opening=opening, ax=ax, ay=ay, az=az, + omegaX=omegaX, omegaY=omegaY, omegaZ=omegaZ, v=v, **kwds): + # chain up super().__init__(**kwds) # store the location diff --git a/models/cdm/cdm/libcdm.py b/models/cdm/cdm/libcdm.py index 0ec13bbd..66aad7f4 100644 --- a/models/cdm/cdm/libcdm.py +++ b/models/cdm/cdm/libcdm.py @@ -35,7 +35,7 @@ def norm(v): return numpy.sqrt(v.dot(v)) -def CDM(X, Y, X0, Y0, depth, ax, ay, az, omegaX, omegaY, omegaZ, opening, nu): +def CDM(X, Y, X0, Y0, depth, opening, ax, ay, az, omegaX, omegaY, omegaZ, nu): """ CDM calculates the surface displacements and potency associated with a CDM @@ -219,8 +219,10 @@ def CDM(X, Y, X0, Y0, depth, ax, ay, az, omegaX, omegaY, omegaZ, opening, nu): R3[2], R4[2]]) if numpy.any(VertVec>0): - raise ValueError('Half-space solution: The CDM must be under the free surface!' + - ' VertVec={}'.format(VertVec)) +# raise ValueError('Half-space solution: The CDM must be under the free surface!' + +# ' VertVec={}'.format(VertVec)) + print('Half-space solution: The CDM must be under the free surface!' + + ' VertVec={}'.format(VertVec)) if ax == 0 and ay == 0 and az == 0: ue = numpy.zeros(numpy.shape(X)) diff --git a/models/cdm/ext/cdm/source.cc b/models/cdm/ext/cdm/source.cc index b1c3753f..3a5df0a5 100644 --- a/models/cdm/ext/cdm/source.cc +++ b/models/cdm/ext/cdm/source.cc @@ -366,14 +366,16 @@ layout(PyObject *, PyObject * args) { // storage PyObject * pySource; - std::size_t xIdx, dIdx, openingIdx, aXIdx, omegaXIdx, offsetIdx; + //std::size_t xIdx, dIdx, openingIdx, aXIdx, omegaXIdx, offsetIdx; + std::size_t xIdx, dIdx, openingIdx, aXIdx, aYIdx, aZIdx, omegaXIdx, omegaYIdx, omegaZIdx, offsetIdx; // unpack the arguments int status = PyArg_ParseTuple(args, - "O!kkkkkk:layout", + "O!kkkkkkkkkk:layout", &PyCapsule_Type, &pySource, &xIdx, &dIdx, - &openingIdx, &aXIdx, &omegaXIdx, + &openingIdx, &aXIdx, &aYIdx, &aZIdx, + &omegaXIdx, &omegaYIdx, &omegaZIdx, &offsetIdx); // if something went wrong if (!status) { @@ -401,7 +403,8 @@ layout(PyObject *, PyObject * args) << pyre::journal::endl; // attach the map - source->layout(xIdx, dIdx, openingIdx, aXIdx, omegaXIdx, offsetIdx); + //source->layout(xIdx, dIdx, openingIdx, aXIdx, omegaXIdx, offsetIdx); + source->layout(xIdx, dIdx, openingIdx, aXIdx, aYIdx, aZIdx, omegaXIdx, omegaYIdx, omegaZIdx, offsetIdx); // all done Py_INCREF(Py_None); diff --git a/models/cdm/lib/libcdm/Source.cc b/models/cdm/lib/libcdm/Source.cc index 855b48b2..ffd4bb22 100644 --- a/models/cdm/lib/libcdm/Source.cc +++ b/models/cdm/lib/libcdm/Source.cc @@ -72,6 +72,9 @@ displacements(gsl_matrix_view * samples, gsl_matrix * predicted) const { // clean up the resulting matrix gsl_matrix_set_zero(predicted); + // allocate storage for the displacement vectors; we reuse this for all samples + gsl_matrix * disp = gsl_matrix_alloc(_locations->size1, 3); + // go through all the samples for (auto sample=0; samplematrix, sample, _omegaYIdx); auto omegaZ = gsl_matrix_get(&samples->matrix, sample, _omegaZIdx); - // compute the displacements - cdm(sample, _locations, _los, + // // compute the displacements + // cdm(sample, _locations, _los, + // xSrc, ySrc, dSrc, + // aX, aY, aZ, + // omegaX, omegaY, omegaZ, + // openingSrc, + // _nu, + // predicted); + + // // apply the location specific projection to LOS vector and dataset shift + // for (auto loc=0; loc<_locations->size1; ++loc) { + // // get the current value + // auto u = gsl_matrix_get(predicted, sample, loc); + // // find the shift that corresponds to this observation + // auto shift = gsl_matrix_get(&samples->matrix, sample, _offsetIdx+_oids[loc]); + // // and apply it to the projected displacement + // u -= shift; + // // save + // gsl_matrix_set(predicted, sample, loc, u); + // } + // } + cdm(_locations, xSrc, ySrc, dSrc, + openingSrc, aX, aY, aZ, omegaX, omegaY, omegaZ, - openingSrc, _nu, - predicted); - + disp); + // apply the location specific projection to LOS vector and dataset shift for (auto loc=0; loc<_locations->size1; ++loc) { - // get the current value - auto u = gsl_matrix_get(predicted, sample, loc); + // compute the components of the unit LOS vector + auto nx = gsl_matrix_get(_los, loc, 0); + auto ny = gsl_matrix_get(_los, loc, 1); + auto nz = gsl_matrix_get(_los, loc, 2); + + // get the three components of the predicted displacement for this location + auto ux = gsl_matrix_get(disp, loc, 0); + auto uy = gsl_matrix_get(disp, loc, 1); + auto ud = gsl_matrix_get(disp, loc, 2); + + // project + auto u = ux*nx + uy*ny + ud*nz; // find the shift that corresponds to this observation auto shift = gsl_matrix_get(&samples->matrix, sample, _offsetIdx+_oids[loc]); // and apply it to the projected displacement u -= shift; + // save gsl_matrix_set(predicted, sample, loc, u); } } + // clean up + gsl_matrix_free(disp); + // all done return; } @@ -154,3 +191,4 @@ residuals(gsl_matrix * predicted) const { // end of file + diff --git a/models/cdm/lib/libcdm/Source.h b/models/cdm/lib/libcdm/Source.h index 07635749..f6e17a32 100644 --- a/models/cdm/lib/libcdm/Source.h +++ b/models/cdm/lib/libcdm/Source.h @@ -44,10 +44,14 @@ class altar::models::cdm::Source { inline void locations(gsl_matrix * locations); inline void los(gsl_matrix * los); inline void oids(const oids_type & oids); +// inline void layout(size_type xIdx, size_type dIdx, +// size_type openingIdx, size_type aXIdx, size_type omegaXIdx, +// size_type offsetIdx); inline void layout(size_type xIdx, size_type dIdx, - size_type openingIdx, size_type aXIdx, size_type omegaXIdx, - size_type offsetIdx); + size_type openingIdx, size_type aXIdx, size_type aYIdx, size_type aZIdx, + size_type omegaXIdx, size_type omegaYIdx, size_type omegaZIdx, size_type offsetIdx); + void displacements(gsl_matrix_view * samples, gsl_matrix * predicted) const; void residuals(gsl_matrix * predicted) const; diff --git a/models/cdm/lib/libcdm/Source.icc b/models/cdm/lib/libcdm/Source.icc index 4d299a8f..3342e108 100644 --- a/models/cdm/lib/libcdm/Source.icc +++ b/models/cdm/lib/libcdm/Source.icc @@ -110,19 +110,23 @@ oids(const oids_type & oids) { void altar::models::cdm::Source:: layout(size_type xIdx, size_type dIdx, - size_type openingIdx, size_type aXIdx, size_type omegaXIdx, - size_type offsetIdx) { + size_type openingIdx, size_type aXIdx, size_type aYIdx, size_type aZIdx, size_type omegaXIdx, + size_type omegaYIdx, size_type omegaZIdx, size_type offsetIdx) { // assign _xIdx = xIdx; _yIdx = xIdx + 1; _dIdx = dIdx; _openingIdx = openingIdx; _aXIdx = aXIdx; - _aYIdx = aXIdx + 1; - _aZIdx = aXIdx + 2; + //_aYIdx = aXIdx + 1; + //_aZIdx = aXIdx + 2; + _aYIdx = aYIdx; + _aZIdx = aZIdx; _omegaXIdx = omegaXIdx; - _omegaYIdx = omegaXIdx + 1; - _omegaZIdx = omegaXIdx + 2; + //_omegaYIdx = omegaXIdx + 1; + //_omegaZIdx = omegaXIdx + 2; + _omegaYIdx = omegaYIdx; + _omegaZIdx = omegaZIdx; _offsetIdx = offsetIdx; // make a channel diff --git a/models/cdm/lib/libcdm/cdm.cc b/models/cdm/lib/libcdm/cdm.cc index 6bd26abe..03c598b5 100644 --- a/models/cdm/lib/libcdm/cdm.cc +++ b/models/cdm/lib/libcdm/cdm.cc @@ -29,8 +29,12 @@ namespace altar { // local helpers static void - RDdispSurf(int sample, - const gsl_matrix * locations, const gsl_matrix * los, + // RDdispSurf(int sample, + // const gsl_matrix * locations, const gsl_matrix * los, + // const vec_t & P1, const vec_t & P2, const vec_t & P3, const vec_t & P4, + // double opening, double nu, + // gsl_matrix * results); + RDdispSurf(const gsl_matrix * locations, const vec_t & P1, const vec_t & P2, const vec_t & P3, const vec_t & P4, double opening, double nu, gsl_matrix * results); @@ -76,14 +80,21 @@ namespace altar { // definitions void altar::models::cdm:: -cdm(int sample, - const gsl_matrix * locations, const gsl_matrix * los, +cdm(const gsl_matrix * locations, double x, double y, double depth, + double opening, double aX, double aY, double aZ, double omegaX, double omegaY, double omegaZ, - double opening, double nu, gsl_matrix * predicted) +// cdm(int sample, +// const gsl_matrix * locations, const gsl_matrix * los, +// double x, double y, double depth, +// double aX, double aY, double aZ, +// double omegaX, double omegaY, double omegaZ, +// double opening, +// double nu, +// gsl_matrix * predicted) { // convert semi-axes to axes aX *= 2; @@ -91,12 +102,12 @@ cdm(int sample, aZ *= 2; // short circuit the trivial case - if (std::abs(aX) < eps && std::abs(aY) < eps && std::abs(aZ) < eps) { - // no displacements - gsl_matrix_set_zero(predicted); - // all done - return; - } + // if (std::abs(aX) < eps && std::abs(aY) < eps && std::abs(aZ) < eps) { + // // no displacements + // gsl_matrix_set_zero(predicted); + // // all done + // return; + // } // the axis specific coordinate matrices mat_t Rx = {1., 0., 0., @@ -141,22 +152,70 @@ cdm(int sample, Q1[2] > 0 || Q2[2] > 0 || Q3[2] > 0 || Q4[2] > 0 || R1[2] > 0 || R2[2] > 0 || R3[2] > 0 || R4[2] > 0) { // complain... - throw std::domain_error("the CDM must be below the surface"); + // throw std::domain_error("the CDM must be below the surface"); + printf("WARNING: The CDM must be below the surface \n"); } + // // dispatch the various cases + // if (std::abs(aX) < eps && std::abs(aY) > eps && std::abs(aZ) > eps) { + // RDdispSurf(sample, locations, los, P1, P2, P3, P4, opening, nu, predicted); + // } else if (std::abs(aX) > eps && std::abs(aY) < eps && std::abs(aZ) > eps) { + // RDdispSurf(sample, locations, los, Q1, Q2, Q3, Q4, opening, nu, predicted); + // } else if (std::abs(aX) > eps && std::abs(aY) > eps && std::abs(aZ) < eps) { + // RDdispSurf(sample, locations, los, R1, R2, R3, R4, opening, nu, predicted); + // } else { + // RDdispSurf(sample, locations, los, P1, P2, P3, P4, opening, nu, predicted); + // RDdispSurf(sample, locations, los, Q1, Q2, Q3, Q4, opening, nu, predicted); + // RDdispSurf(sample, locations, los, R1, R2, R3, R4, opening, nu, predicted); + // } + // dispatch the various cases - if (std::abs(aX) < eps && std::abs(aY) > eps && std::abs(aZ) > eps) { - RDdispSurf(sample, locations, los, P1, P2, P3, P4, opening, nu, predicted); - } else if (std::abs(aX) > eps && std::abs(aY) < eps && std::abs(aZ) > eps) { - RDdispSurf(sample, locations, los, Q1, Q2, Q3, Q4, opening, nu, predicted); - } else if (std::abs(aX) > eps && std::abs(aY) > eps && std::abs(aZ) < eps) { - RDdispSurf(sample, locations, los, R1, R2, R3, R4, opening, nu, predicted); + // short circuit the trivial case + if (std::abs(aX) == 0 && std::abs(aY) == 0 && std::abs(aZ) == 0) { + // no displacements + printf("I am if\n"); + gsl_matrix_set_zero(predicted); + } + else if (std::abs(aX) == 0 && std::abs(aY) > 0 && std::abs(aZ) > 0) { + printf("I am elseif1 \n"); + RDdispSurf(locations, P1, P2, P3, P4, opening, nu, predicted); + } else if (std::abs(aX) > 0 && std::abs(aY) == 0 && std::abs(aZ) > 0) { + printf("I am elseif2 \n"); + RDdispSurf(locations, Q1, Q2, Q3, Q4, opening, nu, predicted); + } else if (std::abs(aX) > 0 && std::abs(aY) > 0 && std::abs(aZ) == 0) { + printf("I am elseif3 \n"); + RDdispSurf(locations, R1, R2, R3, R4, opening, nu, predicted); } else { - RDdispSurf(sample, locations, los, P1, P2, P3, P4, opening, nu, predicted); - RDdispSurf(sample, locations, los, Q1, Q2, Q3, Q4, opening, nu, predicted); - RDdispSurf(sample, locations, los, R1, R2, R3, R4, opening, nu, predicted); + + // allocate room for the partial results + gsl_matrix * P = gsl_matrix_alloc(locations->size1, 3); + gsl_matrix * Q = gsl_matrix_alloc(locations->size1, 3); + gsl_matrix * R = gsl_matrix_alloc(locations->size1, 3); + + // compute + RDdispSurf(locations, P1, P2, P3, P4, opening, nu, P); + RDdispSurf(locations, Q1, Q2, Q3, Q4, opening, nu, Q); + RDdispSurf(locations, R1, R2, R3, R4, opening, nu, R); + + // assemble + for (auto loc=0; locsize1; ++loc) { + for (auto axis=0; axis<3; ++axis) { + // combine + auto result = + gsl_matrix_get(P, loc, axis) + + gsl_matrix_get(Q, loc, axis) + + gsl_matrix_get(R, loc, axis); + + // assign + gsl_matrix_set(predicted, loc, axis, result); + } + } + gsl_matrix_free(P); + gsl_matrix_free(Q); + gsl_matrix_free(R); } + // all done return; } @@ -164,11 +223,16 @@ cdm(int sample, // implementations static void altar::models::cdm:: -RDdispSurf(int sample, - const gsl_matrix * locations, const gsl_matrix * los, +// RDdispSurf(int sample, +// const gsl_matrix * locations, const gsl_matrix * los, +// const vec_t & P1, const vec_t & P2, const vec_t & P3, const vec_t & P4, +// double opening, double nu, +// gsl_matrix * results) { +RDdispSurf(const gsl_matrix * locations, const vec_t & P1, const vec_t & P2, const vec_t & P3, const vec_t & P4, double opening, double nu, gsl_matrix * results) { + // cross auto V = cross(P2-P1, P4-P1); auto b = opening * V/norm(V); @@ -184,23 +248,30 @@ RDdispSurf(int sample, auto u3 = AngSetupFSC(x,y, b, P3, P4, nu); auto u4 = AngSetupFSC(x,y, b, P4, P1, nu); + // // assemble + // auto u = u1 + u2 + u3 + u4; + // // compute the unit LOS vector + // vec_t n = { gsl_matrix_get(los, loc, 0), + // gsl_matrix_get(los, loc, 1), + // gsl_matrix_get(los, loc, 2) }; + + // // project the displacement to the LOS + // auto uLOS = dot(u, n); + // // save by accumulating my contribution to the slot + // // N.B.: note the "+=": the general case call this function three times + // // get the current value + // auto current = gsl_matrix_get(results, sample, loc); + // // update + // current += uLOS; + // // save + // gsl_matrix_set(results, sample, loc, current); + // } + // assemble - auto u = u1 + u2 + u3 + u4; - // compute the unit LOS vector - vec_t n = { gsl_matrix_get(los, loc, 0), - gsl_matrix_get(los, loc, 1), - gsl_matrix_get(los, loc, 2) }; - - // project the displacement to the LOS - auto uLOS = dot(u, n); - // save by accumulating my contribution to the slot - // N.B.: note the "+=": the general case call this function three times - // get the current value - auto current = gsl_matrix_get(results, sample, loc); - // update - current += uLOS; - // save - gsl_matrix_set(results, sample, loc, current); + for (auto axis=0; axis<3; ++axis) { + auto u = u1[axis] + u2[axis] + u3[axis] + u4[axis]; + gsl_matrix_set(results, loc, axis, u); + } } // all done @@ -215,7 +286,8 @@ AngSetupFSC(double x, double y, double nu) { vec_t SideVec = PB - PA; vec_t eZ = {0, 0, 1}; - auto beta = std::acos(dot(SideVec, eZ) / norm(SideVec)); + // auto beta = std::acos(dot(SideVec, eZ) / norm(SideVec)); + auto beta = std::acos(dot(-SideVec, eZ) / norm(SideVec)); if (std::abs(beta) < eps || std::abs(pi - beta) < eps) { return { 0,0,0 }; @@ -239,14 +311,20 @@ AngSetupFSC(double x, double y, vec_t vA, vB; // distinguish the two configurations - if (beta*adcsA[0] > 0) { + // if (beta*adcsA[0] > 0) { + if (beta*adcsA[0] >= 0) { // configuration I - vA = AngDisDispSurf(adcsA, -pi+beta, b, nu, -PA[2]); - vB = AngDisDispSurf(adcsB, -pi+beta, b, nu, -PB[2]); + // vA = AngDisDispSurf(adcsA, -pi+beta, b, nu, -PA[2]); + // vB = AngDisDispSurf(adcsB, -pi+beta, b, nu, -PB[2]); + vA = AngDisDispSurf(adcsA, -pi+beta, bADCS, nu, -PA[2]); + vB = AngDisDispSurf(adcsB, -pi+beta, bADCS, nu, -PB[2]); + } else { // configuration II - vA = AngDisDispSurf(adcsA, beta, b, nu, -PB[2]); - vB = AngDisDispSurf(adcsB, beta, b, nu, -PB[2]); + // vA = AngDisDispSurf(adcsA, beta, b, nu, -PB[2]); + // vB = AngDisDispSurf(adcsB, beta, b, nu, -PB[2]); + vA = AngDisDispSurf(adcsA, beta, bADCS, nu, -PA[2]); + vB = AngDisDispSurf(adcsB, beta, bADCS, nu, -PB[2]); } vec_t v = xform(transpose(A), vB - vA); @@ -432,3 +510,4 @@ cos(double omega) { }; // end-of-file + diff --git a/models/cdm/lib/libcdm/cdm.h b/models/cdm/lib/libcdm/cdm.h index bfd6eb9a..c5526a54 100644 --- a/models/cdm/lib/libcdm/cdm.h +++ b/models/cdm/lib/libcdm/cdm.h @@ -12,12 +12,20 @@ namespace altar { namespace models { namespace cdm { - void cdm(int sample, - const gsl_matrix * locations, const gsl_matrix * los, + // void cdm(int sample, + // const gsl_matrix * locations, const gsl_matrix * los, + // double X0, double Y0, double depth, + // double ax, double ay, double az, + // double omegaX, double omegaY, double omegaZ, + // double opening, + // double nu, + // gsl_matrix * predicted + // ); + void cdm(const gsl_matrix * locations, + double opening, double X0, double Y0, double depth, double ax, double ay, double az, double omegaX, double omegaY, double omegaZ, - double opening, double nu, gsl_matrix * predicted ); @@ -28,3 +36,4 @@ namespace altar { #endif // end of file + diff --git a/models/cdm/tests/libcdm.py b/models/cdm/tests/libcdm.py index 05d6006b..a0eb8f38 100755 --- a/models/cdm/tests/libcdm.py +++ b/models/cdm/tests/libcdm.py @@ -17,7 +17,8 @@ import sys import numpy -from altar.models.cdm.libcdm import CDM +#from altar.models.cdm.libcdm import CDM +from altar.models.cdm.ext import libcdm # Matlab test coordinates and displacements in EFCS (East, North, Vertical). # Sampled from computed output from original Matlab code.