diff --git a/altar/bayesian/Metropolis.py b/altar/bayesian/Metropolis.py index c47e7a57..2a3da331 100644 --- a/altar/bayesian/Metropolis.py +++ b/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) diff --git a/models/cdm/cdm/CDM.py b/models/cdm/cdm/CDM.py index ad5544e8..768371fa 100644 --- a/models/cdm/cdm/CDM.py +++ b/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) diff --git a/models/cdm/cdm/Fast.py b/models/cdm/cdm/Fast.py index bca6ab59..a57c8338 100644 --- a/models/cdm/cdm/Fast.py +++ b/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 diff --git a/models/cdm/cdm/Native.py b/models/cdm/cdm/Native.py index 2a86963c..b504e94b 100644 --- a/models/cdm/cdm/Native.py +++ b/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 diff --git a/models/cdm/cdm/Source.py b/models/cdm/cdm/Source.py index 5cb4c248..0630126a 100644 --- a/models/cdm/cdm/Source.py +++ b/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 diff --git a/models/cdm/cdm/libcdm.py b/models/cdm/cdm/libcdm.py index 4244884e..a247cafc 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, 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 diff --git a/models/cdm/examples/synthetic/cdm.py b/models/cdm/examples/synthetic/cdm.py index df7251b8..f670ed23 100755 --- a/models/cdm/examples/synthetic/cdm.py +++ b/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 diff --git a/models/cdm/lib/libcdm/Source.cc b/models/cdm/lib/libcdm/Source.cc index aa22c07f..adda5aed 100644 --- a/models/cdm/lib/libcdm/Source.cc +++ b/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 diff --git a/models/cdm/lib/libcdm/cdm.cc b/models/cdm/lib/libcdm/cdm.cc index 12c16789..6f9eade6 100644 --- a/models/cdm/lib/libcdm/cdm.cc +++ b/models/cdm/lib/libcdm/cdm.cc @@ -77,25 +77,18 @@ void altar::models::cdm:: 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) { + // convert semi-axes to axes aX *= 2; aY *= 2; 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; - } - // the axis specific coordinate matrices mat_t Rx = {1., 0., 0., 0., cos(omegaX), sin(omegaX), @@ -108,9 +101,10 @@ cdm(const gsl_matrix * locations, mat_t Rz = { cos(omegaZ), sin(omegaZ), 0., -sin(omegaZ), cos(omegaZ), 0., 0., 0., 1.}; - + // the coordinate rotation matrix mat_t R = Rz * (Ry * Rx); + // extract its three columns vec_t R_0 = { R[0], R[3], R[6] }; vec_t R_1 = { R[0+1], R[3+1], R[6+1] }; @@ -118,46 +112,58 @@ cdm(const gsl_matrix * locations, // the centroid vec_t P0 = { x, y, -depth }; - - vec_t P1 = P0 + (aY*R_1 + aZ*R_2)/2; + vec_t P1 = P0 + aY*R_1/2 + aZ*R_2/2; vec_t P2 = P1 - aY*R_1; vec_t P3 = P2 - aZ*R_2; vec_t P4 = P1 - aZ*R_2; - - vec_t Q1 = P0 + (aZ*R_2 - aX*R_0)/2; + + vec_t Q1 = P0 - aX*R_0/2 + aZ*R_2/2; vec_t Q2 = Q1 + aX*R_0; vec_t Q3 = Q2 - aZ*R_2; vec_t Q4 = Q1 - aZ*R_2; - vec_t R1 = P0 + (aX*R_0 + aY*R_1)/2; + vec_t R1 = P0 + aX*R_0/2 + aY*R_1/2; vec_t R2 = R1 - aX*R_0; vec_t R3 = R2 - aY*R_1; vec_t R4 = R1 - aY*R_1; - + // check that all z components are negative if (P1[2] > 0 || P2[2] > 0 || P3[2] > 0 || P4[2] > 0 || 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) { + // 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) > eps && std::abs(aY) < eps && std::abs(aZ) > eps) { + } 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) > eps && std::abs(aY) > eps && std::abs(aZ) < eps) { + } 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 { + // 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) { @@ -166,6 +172,7 @@ cdm(const gsl_matrix * locations, 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); } @@ -186,7 +193,7 @@ RDdispSurf(const gsl_matrix * locations, // cross auto V = cross(P2-P1, P4-P1); auto b = opening * V/norm(V); - + // go through each location for (auto loc=0; locsize1; ++loc) { // unpack the observation point coordinates @@ -218,43 +225,45 @@ AngSetupFSC(double x, double y, 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)); if (std::abs(beta) < eps || std::abs(pi - beta) < eps) { return { 0,0,0 }; } - vec_t ey1 = { SideVec[0], SideVec[1], 0 }; - ey1 = ey1 / norm(ey1); - vec_t ey3 = -eZ; - vec_t ey2 = cross(ey3, ey1); - - mat_t A = { ey1[0], ey1[1], ey1[2], - ey2[0], ey2[1], ey2[2], - ey3[0], ey3[1], ey3[2]}; - - vec_t adcsA = xform(A, {x-PA[0], y-PA[1], -PA[2]}); - vec_t adcsAB = xform(A, SideVec); - vec_t adcsB = adcsA - adcsAB; - - // transform the slip vector - vec_t bADCS = xform(A, b); + else { + vec_t ey1 = { SideVec[0], SideVec[1], 0 }; + ey1 = ey1 / norm(ey1); + vec_t ey3 = -eZ; + vec_t ey2 = cross(ey3, ey1); + + mat_t A = { ey1[0], ey1[1], ey1[2], + ey2[0], ey2[1], ey2[2], + ey3[0], ey3[1], ey3[2]}; + + vec_t adcsA = xform(A, {x-PA[0], y-PA[1], -PA[2]}); + vec_t adcsAB = xform(A, SideVec); + vec_t adcsB = adcsA - adcsAB; + + // transform the slip vector + vec_t bADCS = xform(A, b); + vec_t vA, vB; + + // distinguish the two configurations + if (beta*adcsA[0] >= 0) { + // configuration I + 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, bADCS, nu, -PA[2]); + vB = AngDisDispSurf(adcsB, beta, bADCS, nu, -PB[2]); + } - vec_t vA, vB; - // distinguish the two configurations - 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]); - } else { - // configuration II - vA = AngDisDispSurf(adcsA, beta, b, nu, -PB[2]); - vB = AngDisDispSurf(adcsB, beta, b, nu, -PB[2]); + vec_t v = xform(transpose(A), vB - vA); + + return v; } - - vec_t v = xform(transpose(A), vB - vA); - - return v; } static @@ -281,21 +290,21 @@ AngDisDispSurf(const vec_t & y, double beta, const vec_t & b, // the Burgers function auto Fi = 2*std::atan2(y2, (r+a)/std::tan(beta/2) - y1); - auto v1b1 = b1/2/pi*((1-(1-2*nu)*cotB*cotB)*Fi + + auto v1b1 = b1/2/pi*((1-(1-2*nu)*(cotB*cotB))*Fi + y2/(r+a)*((1-2*nu)*(cotB+y1/2/(r+a))-y1/r) - y2*(r*sinB-y1)*cosB/r/(r-z3)); - auto v2b1 = b1/2/pi*((1-2*nu)*((.5+cotB*cotB)*std::log(r+a)-cotB/sinB*std::log(r-z3)) - - 1./(r+a)*((1-2*nu)*(y1*cotB-a/2-y2*y2/2/(r+a))+y2*y2/r) + - y2*y2*cosB/r/(r-z3)); + auto v2b1 = b1/2/pi*((1-2*nu)*((.5+(cotB*cotB))*std::log(r+a)-cotB/sinB*std::log(r-z3)) - + 1./(r+a)*((1-2*nu)*(y1*cotB-a/2-(y2*y2)/2/(r+a))+(y2*y2)/r) + + (y2*y2)*cosB/r/(r-z3)); auto v3b1 = b1/2/pi*((1-2*nu)*Fi*cotB+y2/(r+a)*(2*nu+a/r) - y2*cosB/(r-z3)*(cosB+a/r)); - auto v1b2 = b2/2/pi*(-(1-2*nu)*((.5-cotB*cotB)*std::log(r+a) + cotB*cotB*cosB*std::log(r-z3) ) - - 1/(r+a)*((1-2*nu)*(y1*cotB+.5*a+y1*y1/2/(r+a)) - y1*y1/r) + + auto v1b2 = b2/2/pi*(-(1-2*nu)*((.5-(cotB*cotB))*std::log(r+a) + (cotB*cotB)*cosB*std::log(r-z3) ) - + 1/(r+a)*((1-2*nu)*(y1*cotB+.5*a+(y1*y1)/2/(r+a)) - (y1*y1)/r) + z1*(r*sinB-y1)/r/(r-z3)); - auto v2b2 = b2/2/pi*((1+(1-2*nu)*cotB*cotB)*Fi - + auto v2b2 = b2/2/pi*((1+(1-2*nu)*(cotB*cotB))*Fi - y2/(r+a)*((1-2*nu)*(cotB+y1/2/(r+a))-y1/r) - y2*z1/r/(r-z3)); @@ -303,13 +312,13 @@ AngDisDispSurf(const vec_t & y, double beta, const vec_t & b, y1/(r+a)*(2*nu+a/r) + z1/(r-z3)*(cosB+a/r)); auto v1b3 = b3/2/pi*(y2*(r*sinB-y1)*sinB/r/(r-z3)); - auto v2b3 = b3/2/pi*(-y2*y2*sinB/r/(r-z3)); + auto v2b3 = b3/2/pi*(-(y2*y2)*sinB/r/(r-z3)); auto v3b3 = b3/2/pi*(Fi + y2*(r*cosB+a)*sinB/r/(r-z3)); auto v1 = v1b1 + v1b2 + v1b3; auto v2 = v2b1 + v2b2 + v2b3; auto v3 = v3b1 + v3b2 + v3b3; - + return {v1, v2, v3}; } @@ -367,9 +376,9 @@ inline static altar::models::cdm::mat_t altar::models::cdm:: operator*(const mat_t & m1, const mat_t & m2) { - return { m1[0]*m2[0] + m1[2]*m2[3] + m1[2]*m2[6], - m1[0]*m2[1] + m1[2]*m2[4] + m1[2]*m2[7], - m1[0]*m2[2] + m1[2]*m2[5] + m1[2]*m2[8], + return { m1[0]*m2[0] + m1[1]*m2[3] + m1[2]*m2[6], + m1[0]*m2[1] + m1[1]*m2[4] + m1[2]*m2[7], + m1[0]*m2[2] + m1[1]*m2[5] + m1[2]*m2[8], m1[3]*m2[0] + m1[4]*m2[3] + m1[5]*m2[6], m1[3]*m2[1] + m1[4]*m2[4] + m1[5]*m2[7], @@ -414,9 +423,9 @@ inline static altar::models::cdm::vec_t altar::models::cdm:: xform(const mat_t & m, const vec_t & v) { - return {m[0]*v[0] + m[1]*v[2] + m[2]*v[3], - m[3]*v[0] + m[4]*v[2] + m[5]*v[3], - m[6]*v[0] + m[7]*v[2] + m[8]*v[3]}; + return {m[0]*v[0] + m[1]*v[1] + m[2]*v[2], + m[3]*v[0] + m[4]*v[1] + m[5]*v[2], + m[6]*v[0] + m[7]*v[1] + m[8]*v[2]}; }; // trig @@ -424,14 +433,14 @@ inline double altar::models::cdm:: sin(double omega) { - return std::sin(omega * 180/pi); + return std::sin(omega * pi/180); }; inline double altar::models::cdm:: cos(double omega) { - return std::cos(omega * 180/pi); + return std::cos(omega * pi/180); }; // end-of-file diff --git a/models/cdm/lib/libcdm/cdm.h b/models/cdm/lib/libcdm/cdm.h index d709d4b2..a674f726 100644 --- a/models/cdm/lib/libcdm/cdm.h +++ b/models/cdm/lib/libcdm/cdm.h @@ -14,9 +14,9 @@ namespace altar { namespace cdm { void cdm(const gsl_matrix * locations, double X0, double Y0, double depth, + double opening, double ax, double ay, double az, double omegaX, double omegaY, double omegaZ, - double opening, double nu, gsl_matrix * predicted ); diff --git a/models/mogi/lib/libmogi/Source.cc b/models/mogi/lib/libmogi/Source.cc index 2b41b7e3..6a5a81c4 100644 --- a/models/mogi/lib/libmogi/Source.cc +++ b/models/mogi/lib/libmogi/Source.cc @@ -75,7 +75,8 @@ displacements(gsl_matrix_view * samples, gsl_matrix * predicted) const { auto xSrc = gsl_matrix_get(&samples->matrix, sample, _xIdx); auto ySrc = gsl_matrix_get(&samples->matrix, sample, _yIdx); auto dSrc = gsl_matrix_get(&samples->matrix, sample, _dIdx); - auto sSrc = std::pow(10, gsl_matrix_get(&samples->matrix, sample, _sIdx)); + auto sSrc = 1000000*(gsl_matrix_get(&samples->matrix, sample, _sIdx)); + //auto sSrc = std::pow(10, gsl_matrix_get(&samples->matrix, sample, _sIdx)); // go through the locations for (auto loc=0; loc<_locations->size1; ++loc) { diff --git a/models/mogi/mogi/CUDA.py b/models/mogi/mogi/CUDA.py index 8f9c3d5e..696203aa 100644 --- a/models/mogi/mogi/CUDA.py +++ b/models/mogi/mogi/CUDA.py @@ -96,7 +96,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 diff --git a/models/mogi/mogi/Fast.py b/models/mogi/mogi/Fast.py index e7b53392..b410c4d8 100644 --- a/models/mogi/mogi/Fast.py +++ b/models/mogi/mogi/Fast.py @@ -87,7 +87,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 diff --git a/models/mogi/mogi/Mogi.py b/models/mogi/mogi/Mogi.py index 76e7d98d..db850bb2 100644 --- a/models/mogi/mogi/Mogi.py +++ b/models/mogi/mogi/Mogi.py @@ -112,7 +112,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) diff --git a/models/mogi/mogi/Native.py b/models/mogi/mogi/Native.py index 635dad08..e4f1d8b9 100644 --- a/models/mogi/mogi/Native.py +++ b/models/mogi/mogi/Native.py @@ -74,7 +74,8 @@ def dataLikelihood(self, model, step): # its depth d = parameters[dIdx] # and its strength; we model the logarithm of this one, so we have to exponentiate - dV = 10**parameters[sIdx] + #dV = 10**parameters[sIdx] + dV = 1.e6*parameters[sIdx] # make a source using the sample parameters mogi = source(x=x, y=y, d=d, dV=dV) @@ -89,9 +90,9 @@ def dataLikelihood(self, model, step): 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.0 # all done return self