CDM-only related PR

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

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

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

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

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))

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);

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; sample<nSamples; ++sample) {
         // unpack the parameters
@@ -90,28 +93,62 @@ displacements(gsl_matrix_view * samples, gsl_matrix * predicted) const {
         auto omegaY = gsl_matrix_get(&samples->matrix, 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
+

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;
 

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