feat!: reduces model to minimal bayes needed stuff
All checks were successful
gitea-physics/pdme/pipeline/head This commit looks good
All checks were successful
gitea-physics/pdme/pipeline/head This commit looks good
This commit is contained in:
@@ -2,7 +2,6 @@ import numpy
|
||||
import numpy.random
|
||||
from pdme.model.model import Model
|
||||
from pdme.measurement import (
|
||||
DotMeasurement,
|
||||
OscillatingDipole,
|
||||
OscillatingDipoleArrangement,
|
||||
)
|
||||
@@ -30,7 +29,6 @@ class FixedMagnitudeModel(Model):
|
||||
zmin: float,
|
||||
zmax: float,
|
||||
pfixed: float,
|
||||
n: int,
|
||||
) -> None:
|
||||
self.xmin = xmin
|
||||
self.xmax = xmax
|
||||
@@ -39,34 +37,10 @@ class FixedMagnitudeModel(Model):
|
||||
self.zmin = zmin
|
||||
self.zmax = zmax
|
||||
self.pfixed = pfixed
|
||||
self._n = n
|
||||
self.rng = numpy.random.default_rng()
|
||||
|
||||
def __repr__(self) -> str:
|
||||
return f"FixedMagnitudeModel({self.xmin}, {self.xmax}, {self.ymin}, {self.ymax}, {self.zmin}, {self.zmax}, {self.pfixed}, {self.n()})"
|
||||
|
||||
def solution_single_dipole(self, pt: numpy.ndarray) -> OscillatingDipole:
|
||||
# assume length is 6, who needs error checking.
|
||||
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),
|
||||
]
|
||||
)
|
||||
return OscillatingDipole(p, s, w)
|
||||
|
||||
def point_length(self) -> int:
|
||||
"""
|
||||
Dipole is constrained magnitude, but free orientation.
|
||||
Six degrees of freedom: (p_theta, p_phi, sx, sy, sz, w).
|
||||
"""
|
||||
return 6
|
||||
return f"FixedMagnitudeModel({self.xmin}, {self.xmax}, {self.ymin}, {self.ymax}, {self.zmin}, {self.zmax}, {self.pfixed})"
|
||||
|
||||
def get_dipoles(self, frequency: float) -> OscillatingDipoleArrangement:
|
||||
theta = numpy.arccos(self.rng.uniform(-1, 1))
|
||||
@@ -109,24 +83,3 @@ class FixedMagnitudeModel(Model):
|
||||
w = rng.uniform(1, max_frequency, n)
|
||||
|
||||
return numpy.array([px, py, pz, sx, sy, sz, w]).T
|
||||
|
||||
def n(self) -> int:
|
||||
return self._n
|
||||
|
||||
def v_for_point_at_dot(self, dot: DotMeasurement, pt: numpy.ndarray) -> float:
|
||||
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))
|
||||
return alpha**2 * b
|
||||
|
||||
@@ -1,8 +1,6 @@
|
||||
import numpy
|
||||
import numpy.random
|
||||
from typing import Callable, Sequence
|
||||
from pdme.measurement import (
|
||||
DotMeasurement,
|
||||
OscillatingDipoleArrangement,
|
||||
)
|
||||
import logging
|
||||
@@ -16,15 +14,6 @@ class Model:
|
||||
Interface for models.
|
||||
"""
|
||||
|
||||
def point_length(self) -> int:
|
||||
raise NotImplementedError
|
||||
|
||||
def n(self) -> int:
|
||||
raise NotImplementedError
|
||||
|
||||
def v_for_point_at_dot(self, dot: DotMeasurement, pt: numpy.ndarray) -> float:
|
||||
raise NotImplementedError
|
||||
|
||||
def get_dipoles(self, frequency: float) -> OscillatingDipoleArrangement:
|
||||
raise NotImplementedError
|
||||
|
||||
@@ -32,32 +21,3 @@ class Model:
|
||||
self, n: int, max_frequency: float, rng: numpy.random.Generator = None
|
||||
) -> numpy.ndarray:
|
||||
raise NotImplementedError
|
||||
|
||||
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.
|
||||
pt_length = self.point_length()
|
||||
chunked_pts = [pts[i : i + pt_length] for i in range(0, len(pts), pt_length)]
|
||||
return sum(self.v_for_point_at_dot(dot, pt) for pt in chunked_pts) - dot.v
|
||||
|
||||
def costs(
|
||||
self, dots: Sequence[DotMeasurement]
|
||||
) -> Callable[[numpy.ndarray], numpy.ndarray]:
|
||||
"""
|
||||
Returns a function that returns the cost for the given list of DotMeasurements for a particular model-dependent phase space point.
|
||||
Default implementation assumes a single dot cost from which to build the list.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
dots: A list of dot measurements to use to find the cost functions.
|
||||
|
||||
Returns
|
||||
----------
|
||||
Returns the model's cost function.
|
||||
"""
|
||||
_logger.debug(f"Constructing costs for dots: {dots}")
|
||||
|
||||
def costs_to_return(pts: numpy.ndarray) -> numpy.ndarray:
|
||||
return numpy.array([self.cost_for_dot(dot, pts) for dot in dots])
|
||||
|
||||
return costs_to_return
|
||||
|
||||
@@ -1,3 +0,0 @@
|
||||
from pdme.util.normal_form import normalise_point_list
|
||||
|
||||
__all__ = ["normalise_point_list"]
|
||||
|
||||
@@ -1,45 +0,0 @@
|
||||
import numpy
|
||||
import operator
|
||||
import logging
|
||||
|
||||
|
||||
# flips px, py, pz
|
||||
SIGN_ARRAY_7 = numpy.array((-1, -1, -1, 1, 1, 1, 1))
|
||||
SIGN_ARRAY_4 = numpy.array((-1, 1, 1, 1))
|
||||
|
||||
|
||||
_logger = logging.getLogger(__name__)
|
||||
|
||||
|
||||
def flip_chunk_to_positive_px(pt: numpy.ndarray) -> numpy.ndarray:
|
||||
if pt[0] > 0:
|
||||
return pt
|
||||
else:
|
||||
# godawful hack.
|
||||
if len(pt) == 7:
|
||||
return SIGN_ARRAY_7 * pt
|
||||
elif len(pt) == 4:
|
||||
return SIGN_ARRAY_4 * pt
|
||||
else:
|
||||
_logger.warning(f"Could not normalise pt: {pt}. Returning as is...")
|
||||
return pt
|
||||
|
||||
|
||||
def normalise_point_list(pts: numpy.ndarray, pt_length) -> numpy.ndarray:
|
||||
chunked_pts = [
|
||||
flip_chunk_to_positive_px(pts[i : i + pt_length])
|
||||
for i in range(0, len(pts), pt_length)
|
||||
]
|
||||
range_to_length = list(range(pt_length))
|
||||
rotated_range = (
|
||||
range_to_length[pt_length - 1 :] + range_to_length[0 : pt_length - 1]
|
||||
)
|
||||
return numpy.concatenate(
|
||||
sorted(
|
||||
chunked_pts,
|
||||
key=lambda x: tuple(
|
||||
round(val, 3) for val in operator.itemgetter(*rotated_range)(x)
|
||||
),
|
||||
),
|
||||
axis=None,
|
||||
)
|
||||
Reference in New Issue
Block a user