Changes made to make mogi and cdm run properly

altar/bayesian/Metropolis.py

@@ -50,7 +50,7 @@ def initialize(self, application):
         # get the capsule of the random number generator
         rng = application.rng.rng
         # set up the distribution for building the sample multiplicities
-        self.uniform = altar.pdf.uniform(support=(0,1), rng=rng)
+        self.uniform = altar.pdf.uniform(support=(1.0e-15,1), rng=rng)
         # set up the distribution for the random walk displacement vectors
         self.uninormal = altar.pdf.ugaussian(rng=rng)
 

models/cdm/cdm/CDM.py

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

models/cdm/cdm/Fast.py

@@ -49,10 +49,8 @@ def initialize(self, model, **kwds):
         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.xIdx, model.dIdx, model.openingIdx, model.aXIdx, model.omegaXIdx,
                       model.offsetIdx)
-
         # nothing to do
         return self
 
@@ -90,7 +88,8 @@ def dataLikelihood(self, model, step):
             # get the residuals
             residuals = predicted.getRow(sample)
             # compute the norm, and normalize it
-            llk = normalization - norm.eval(v=residuals, sigma_inv=cd_inv) / 2
+            normeval = norm.eval(v=residuals, sigma_inv=cd_inv)
+            llk = normalization - normeval**2.0 / 2.0
             # store it
             dataLLK[sample] = llk
 

models/cdm/cdm/Native.py

@@ -50,20 +50,20 @@ 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
         dIdx = model.dIdx
+        openingIdx = model.openingIdx
         aXIdx = model.aXIdx
         aYIdx = model.aYIdx
         aZIdx = model.aZIdx
         omegaXIdx = model.omegaXIdx
         omegaYIdx = model.omegaYIdx
         omegaZIdx = model.omegaZIdx
-        openingIdx = model.openingIdx
         offsetIdx = model.offsetIdx
-
+        
         # get the observations
         los = model.los
         oid = model.oid
@@ -91,9 +91,10 @@ def dataLikelihood(self, model, step):
             omegaZ = parameters[omegaZIdx]
 
             # make a source using the sample parameters
-            cdm = source(x=x, y=y, d=d,
+            cdm = source(x=x, y=y, d=d, opening=opening,
                          ax=aX, ay=aY, az=aZ, omegaX=omegaX, omegaY=omegaY, omegaZ=omegaZ,
-                         opening=opening, v=model.nu)
+                         v=model.nu)
+
             # compute the expected displacement
             u = cdm.displacements(locations=locations, los=los)
 
@@ -103,11 +104,11 @@ 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
-            residual = 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 - residual/2
+            dataLLK[sample] = normalization - normeval**2.0 /2
 
         # all done
         return self

models/cdm/cdm/Source.py

@@ -78,9 +78,9 @@ def displacements(self, locations, los):
         # allocate space for the result
         u = altar.vector(shape=len(locations))
         # compute the displacements
-        ue, un, uv =  CDM(Xf, Yf, x_src, y_src, d_src,
-                          omegaX_src, omegaY_src, omegaZ_src, ax_src, ay_src, az_src,
-                          opening, v)
+        ue, un, uv =  CDM(Xf, Yf, x_src, y_src, d_src, opening,
+                          ax_src, ay_src, az_src, omegaX_src, omegaY_src, omegaZ_src,
+                          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
@@ -91,8 +91,8 @@ def displacements(self, locations, los):
 
 
     # meta-methods
-    def __init__(self, x=x, y=y, d=d, omegaX=omegaX, omegaY=omegaY, omegaZ=omegaZ,
-                 ax=ax, ay=ay, az=az, 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, omegaX, omegaY, omegaZ, ax, ay, az, 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

models/cdm/examples/synthetic/cdm.py

@@ -17,7 +17,6 @@
 # my model
 import altar.models.cdm
 
-
 # app
 class CDM(altar.application, family="altar.applications.cdm"):
     """
@@ -46,13 +45,13 @@ class CDM(altar.application, family="altar.applications.cdm"):
     omegaX = altar.properties.float(default=0)
     omegaX.doc = "the CDM rotation about the x axis"
 
-    omegaY = altar.properties.float(default=0)
+    omegaY = altar.properties.float(default=-45)
     omegaY.doc = "the CDM rotation about the y axis"
 
     omegaZ = altar.properties.float(default=0)
     omegaZ.doc = "the CDM rotation about the z axis"
 
-    opening = altar.properties.float(default=1e-3)
+    opening = altar.properties.float(default=100.0)
     opening.doc = "the tensile component of the Burgers vector of the dislocation"
 
     nu = altar.properties.float(default=.25)
@@ -97,16 +96,16 @@ def cdm(self):
         observations = len(stations)
         # make a source
         source = altar.models.cdm.source(
-            x=self.x, y=self.y, d=self.d,
+            x=self.x, y=self.y, d=self.d, opening=self.opening,
             ax=self.aX, ay=self.aY, az=self.aZ,
             omegaX=self.omegaX, omegaY=self.omegaY, omegaZ=self.omegaZ,
-            opening=self.opening, v=self.nu)
+            v=self.nu)
 
         # observe all displacements from the same angle for now
         theta = π/4 # the azimuthal angle
         phi = π     # the polar angle
         # build the common projection vector
-        s = sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta)
+        s = -sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta)
 
         # allocate a matrix to hold the projections
         los = altar.matrix(shape=(observations,3))
@@ -131,14 +130,7 @@ def cdm(self):
         for idx, (x,y) in enumerate(stations):
             # make a new entry in the data sheet
             observation = data.pyre_new()
-            # western stations
-            if x < 0:
-                # come from a different data set
-                observation.oid = 1
-            # than
-            else:
-                # eastern stations
-                observation.oid = 0
+            observation.oid = 0
 
             # record the location of this observation
             observation.x = x
@@ -158,7 +150,7 @@ def cdm(self):
         # go through the observations
         for idx, observation in enumerate(data):
             # set the covariance to a fraction of the "observed" displacement
-            correlation[idx,idx] = 1.0 #.01 * observation.d
+            correlation[idx,idx] = 0.0001 #1.0 #.01 * observation.d
 
         # all done
         return data, correlation
@@ -169,11 +161,16 @@ def makeStations(self):
         Create a set of station coordinate
         """
         # get some help
-        import itertools
-        # build a set of points on a grid
-        stations = itertools.product(range(-5,6), range(-5,6))
+        import numpy as np
+        xx = np.linspace(-6000,6000, 15)
+        yy = np.linspace(-6000,6000, 15)
+        x,y = np.meshgrid(xx, yy)
+        x=x.flatten()
+        y = y.flatten()
+        stations = [tuple((y[i], x[i])) for i in range(0,len(x))]
+
         # and return it
-        return tuple((x*1000, y*1000) for x,y in stations)
+        return stations
 
 
 # bootstrap

models/cdm/lib/libcdm/Source.cc

@@ -90,16 +90,14 @@ 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(_locations,
             xSrc, ySrc, dSrc,
+            openingSrc,
             aX, aY, aZ,
             omegaX, omegaY, omegaZ,
-            openingSrc,
             _nu,
             disp);
-
+        
         // apply the location specific projection to LOS vector and dataset shift
         for (auto loc=0; loc<_locations->size1; ++loc) {
             // compute the components of the unit LOS vector
@@ -109,16 +107,16 @@ displacements(gsl_matrix_view * samples, gsl_matrix * predicted) const {
 
             // 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, 0);
-            auto ud = gsl_matrix_get(disp, loc, 0);
+            auto uy = gsl_matrix_get(disp, loc, 1);
+            auto ud = gsl_matrix_get(disp, loc, 2);
 
-            // project; don't forget {ud} is positive into the ground
-            auto u = ux*nx + uy*ny - ud*nz;
+            // 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);
         }
@@ -147,7 +145,6 @@ residuals(gsl_matrix * predicted) const {
     // unpack the number of samples and number of observations
     auto nSamples = predicted->size1;
     auto nObservations = predicted->size2;
-
     // go through all observations
     for (auto obs=0; obs < nObservations; ++obs) {
         // get the corresponding measurement