feat!: Guts the model interface for only the things useful for monte carlo bayes stuff
This commit is contained in:
@@ -130,52 +130,3 @@ class FixedMagnitudeModel(Model):
|
||||
alpha = p.dot(diff) / (numpy.linalg.norm(diff) ** 3)
|
||||
b = (1 / numpy.pi) * (w / (w**2 + dot.f**2))
|
||||
return alpha**2 * b
|
||||
|
||||
def jac_for_point_at_dot(
|
||||
self, dot: DotMeasurement, pt: numpy.ndarray
|
||||
) -> numpy.ndarray:
|
||||
p_theta = pt[0]
|
||||
p_phi = pt[1]
|
||||
s = pt[2:5]
|
||||
w = pt[5]
|
||||
|
||||
p = numpy.array(
|
||||
[
|
||||
self.pfixed * numpy.sin(p_theta) * numpy.cos(p_phi),
|
||||
self.pfixed * numpy.sin(p_theta) * numpy.sin(p_phi),
|
||||
self.pfixed * numpy.cos(p_theta),
|
||||
]
|
||||
)
|
||||
diff = dot.r - s
|
||||
alpha = p.dot(diff) / (numpy.linalg.norm(diff) ** 3)
|
||||
b = (1 / numpy.pi) * (w / (w**2 + dot.f**2))
|
||||
|
||||
theta_div_middle = self.pfixed * (
|
||||
diff[0] * numpy.cos(p_phi) * numpy.cos(p_theta)
|
||||
+ diff[1] * numpy.sin(p_phi) * numpy.cos(p_theta)
|
||||
- diff[2] * numpy.sin(p_theta)
|
||||
)
|
||||
theta_div = 2 * alpha * (theta_div_middle) / (numpy.linalg.norm(diff) ** 3) * b
|
||||
|
||||
phi_div_middle = self.pfixed * (
|
||||
diff[1] * numpy.sin(p_theta) * numpy.cos(p_phi)
|
||||
- diff[0] * numpy.sin(p_theta) * numpy.sin(p_phi)
|
||||
)
|
||||
phi_div = 2 * alpha * (phi_div_middle) / (numpy.linalg.norm(diff) ** 3) * b
|
||||
|
||||
r_divs = (
|
||||
(
|
||||
-p / (numpy.linalg.norm(diff) ** 3)
|
||||
+ 3 * p.dot(diff) * diff / (numpy.linalg.norm(diff) ** 5)
|
||||
)
|
||||
* 2
|
||||
* alpha
|
||||
* b
|
||||
)
|
||||
|
||||
f2 = dot.f**2
|
||||
w2 = w**2
|
||||
|
||||
w_div = alpha**2 * (1 / numpy.pi) * ((f2 - w2) / ((f2 + w2) ** 2))
|
||||
|
||||
return numpy.concatenate((theta_div, phi_div, r_divs, w_div), axis=None)
|
||||
|
||||
@@ -1,13 +1,10 @@
|
||||
import numpy
|
||||
import numpy.random
|
||||
import scipy.optimize
|
||||
from typing import Callable, Sequence, List
|
||||
from typing import Callable, Sequence
|
||||
from pdme.measurement import (
|
||||
DotMeasurement,
|
||||
OscillatingDipoleArrangement,
|
||||
OscillatingDipole,
|
||||
)
|
||||
import pdme.util
|
||||
import logging
|
||||
|
||||
|
||||
@@ -36,14 +33,6 @@ class Model:
|
||||
) -> numpy.ndarray:
|
||||
raise NotImplementedError
|
||||
|
||||
def solution_single_dipole(self, pt: numpy.ndarray) -> OscillatingDipole:
|
||||
raise NotImplementedError
|
||||
|
||||
def solution_as_dipoles(self, pts: numpy.ndarray) -> List[OscillatingDipole]:
|
||||
pt_length = self.point_length()
|
||||
chunked_pts = [pts[i : i + pt_length] for i in range(0, len(pts), pt_length)]
|
||||
return [self.solution_single_dipole(pt) for pt in chunked_pts]
|
||||
|
||||
def cost_for_dot(self, dot: DotMeasurement, pts: numpy.ndarray) -> float:
|
||||
# creates numpy.ndarrays in groups of self.point_length().
|
||||
# Will throw problems for irregular points, but that's okay for now.
|
||||
@@ -72,67 +61,3 @@ class Model:
|
||||
return numpy.array([self.cost_for_dot(dot, pts) for dot in dots])
|
||||
|
||||
return costs_to_return
|
||||
|
||||
def jac_for_point_at_dot(
|
||||
self, dot: DotMeasurement, pt: numpy.ndarray
|
||||
) -> numpy.ndarray:
|
||||
raise NotImplementedError
|
||||
|
||||
def jac_for_dot(self, dot: DotMeasurement, pts: numpy.ndarray) -> numpy.ndarray:
|
||||
# creates numpy.ndarrays in groups of self.point_length().
|
||||
# Will throw problems for irregular points, but that's okay for now.
|
||||
pt_length = self.point_length()
|
||||
chunked_pts = [pts[i : i + pt_length] for i in range(0, len(pts), pt_length)]
|
||||
return numpy.append(
|
||||
[], [self.jac_for_point_at_dot(dot, pt) for pt in chunked_pts]
|
||||
)
|
||||
|
||||
def jac(
|
||||
self, dots: Sequence[DotMeasurement]
|
||||
) -> Callable[[numpy.ndarray], numpy.ndarray]:
|
||||
"""
|
||||
Returns a function that returns the cost function's Jacobian for the given list of DotMeasurements for a particular model-dependent phase space point.
|
||||
Default implementation assumes a single dot jacobian from which to build the list.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
dots: A list of dot measurements to use to find the cost functions and their Jacobian.
|
||||
|
||||
Returns
|
||||
----------
|
||||
Returns the model's cost function's Jacobian.
|
||||
"""
|
||||
|
||||
def jac_to_return(pts: numpy.ndarray) -> numpy.ndarray:
|
||||
return numpy.array([self.jac_for_dot(dot, pts) for dot in dots])
|
||||
|
||||
return jac_to_return
|
||||
|
||||
def solve(
|
||||
self,
|
||||
dots: Sequence[DotMeasurement],
|
||||
initial_pt: numpy.ndarray = None,
|
||||
bounds=(-numpy.inf, numpy.inf),
|
||||
) -> scipy.optimize.OptimizeResult:
|
||||
if initial_pt is None:
|
||||
initial = numpy.tile(0.1, self.n() * self.point_length())
|
||||
else:
|
||||
if len(initial_pt) != self.point_length():
|
||||
raise ValueError(
|
||||
f"The initial point {initial_pt} does not have the model's expected length: {self.point_length()}"
|
||||
)
|
||||
initial = numpy.tile(initial_pt, self.n())
|
||||
|
||||
result = scipy.optimize.least_squares(
|
||||
self.costs(dots),
|
||||
initial,
|
||||
jac=self.jac(dots),
|
||||
ftol=1e-15,
|
||||
gtol=3e-16,
|
||||
xtol=None,
|
||||
bounds=bounds,
|
||||
)
|
||||
result.normalised_x = pdme.util.normalise_point_list(
|
||||
result.x, self.point_length()
|
||||
)
|
||||
return result
|
||||
|
||||
Reference in New Issue
Block a user