style: fmt updates
All checks were successful
gitea-physics/pdme/pipeline/head This commit looks good

This commit is contained in:
2022-04-03 17:59:58 -05:00
parent b3adf33f59
commit d89d3585da
33 changed files with 956 additions and 358 deletions

View File

@@ -1,3 +1,3 @@
[flake8] [flake8]
ignore = W191, E501, W503 ignore = W191, E501, W503, E203
max-line-length = 120 max-line-length = 120

View File

@@ -1,4 +1,7 @@
from pdme.inputs.dot_inputs import inputs_with_frequency_range, input_pairs_with_frequency_range from pdme.inputs.dot_inputs import (
inputs_with_frequency_range,
input_pairs_with_frequency_range,
)
__all__ = ['inputs_with_frequency_range', 'input_pairs_with_frequency_range'] __all__ = ["inputs_with_frequency_range", "input_pairs_with_frequency_range"]

View File

@@ -4,10 +4,18 @@ import itertools
from typing import Sequence, Tuple from typing import Sequence, Tuple
def inputs_with_frequency_range(dots: Sequence[numpy.typing.ArrayLike], frequency_range: Sequence[float]) -> Sequence[Tuple[numpy.typing.ArrayLike, float]]: def inputs_with_frequency_range(
dots: Sequence[numpy.typing.ArrayLike], frequency_range: Sequence[float]
) -> Sequence[Tuple[numpy.typing.ArrayLike, float]]:
return list(itertools.chain(*[[(dot, f) for f in frequency_range] for dot in dots])) return list(itertools.chain(*[[(dot, f) for f in frequency_range] for dot in dots]))
def input_pairs_with_frequency_range(dots: Sequence[numpy.typing.ArrayLike], frequency_range: Sequence[float]) -> Sequence[Tuple[numpy.typing.ArrayLike, numpy.typing.ArrayLike, float]]: def input_pairs_with_frequency_range(
dots: Sequence[numpy.typing.ArrayLike], frequency_range: Sequence[float]
) -> Sequence[Tuple[numpy.typing.ArrayLike, numpy.typing.ArrayLike, float]]:
all_pairs = itertools.combinations(dots, 2) all_pairs = itertools.combinations(dots, 2)
return list(itertools.chain(*[[(dot1, dot2, f) for f in frequency_range] for (dot1, dot2) in all_pairs])) return list(
itertools.chain(
*[[(dot1, dot2, f) for f in frequency_range] for (dot1, dot2) in all_pairs]
)
)

View File

@@ -1,7 +1,22 @@
from pdme.measurement.dot_measure import DotMeasurement, DotRangeMeasurement from pdme.measurement.dot_measure import DotMeasurement, DotRangeMeasurement
from pdme.measurement.dot_pair_measure import DotPairMeasurement, DotPairRangeMeasurement from pdme.measurement.dot_pair_measure import (
from pdme.measurement.oscillating_dipole import OscillatingDipole, OscillatingDipoleArrangement DotPairMeasurement,
DotPairRangeMeasurement,
)
from pdme.measurement.oscillating_dipole import (
OscillatingDipole,
OscillatingDipoleArrangement,
)
from pdme.measurement.input_types import DotInput, DotPairInput from pdme.measurement.input_types import DotInput, DotPairInput
__all__ = ['DotMeasurement', 'DotRangeMeasurement', 'DotPairMeasurement', 'DotPairRangeMeasurement', 'OscillatingDipole', 'OscillatingDipoleArrangement', 'DotInput', 'DotPairInput'] __all__ = [
"DotMeasurement",
"DotRangeMeasurement",
"DotPairMeasurement",
"DotPairRangeMeasurement",
"OscillatingDipole",
"OscillatingDipoleArrangement",
"DotInput",
"DotPairInput",
]

View File

@@ -4,19 +4,20 @@ import numpy.typing
@dataclass @dataclass
class DotMeasurement(): class DotMeasurement:
''' """
Representation of a dot measuring oscillating dipoles. Representation of a dot measuring oscillating dipoles.
Parameters Parameters
---------- ----------
v : float v : float
The voltage measured at the dot. The voltage measured at the dot.
r : numpy.ndarray r : numpy.ndarray
The position of the dot. The position of the dot.
f : float f : float
The measurement frequency. The measurement frequency.
''' """
v: float v: float
r: numpy.ndarray r: numpy.ndarray
f: float f: float
@@ -26,21 +27,22 @@ class DotMeasurement():
@dataclass @dataclass
class DotRangeMeasurement(): class DotRangeMeasurement:
''' """
Representation of a dot measuring oscillating dipoles. Representation of a dot measuring oscillating dipoles.
Parameters Parameters
---------- ----------
v_low : float v_low : float
The lower range of voltage measured at the dot. The lower range of voltage measured at the dot.
v_high : float v_high : float
The upper range of voltage measured at the dot. The upper range of voltage measured at the dot.
r : numpy.ndarray r : numpy.ndarray
The position of the dot. The position of the dot.
f : float f : float
The measurement frequency. The measurement frequency.
''' """
v_low: float v_low: float
v_high: float v_high: float
r: numpy.ndarray r: numpy.ndarray

View File

@@ -4,21 +4,22 @@ import numpy.typing
@dataclass @dataclass
class DotPairMeasurement(): class DotPairMeasurement:
''' """
Representation of a dot measuring oscillating dipoles. Representation of a dot measuring oscillating dipoles.
Parameters Parameters
---------- ----------
v : float v : float
The voltage measured at the dot. The voltage measured at the dot.
r1 : numpy.ndarray r1 : numpy.ndarray
The position of the first dot. The position of the first dot.
r2 : numpy.ndarray r2 : numpy.ndarray
The position of the second dot. The position of the second dot.
f : float f : float
The measurement frequency. The measurement frequency.
''' """
v: float v: float
r1: numpy.ndarray r1: numpy.ndarray
r2: numpy.ndarray r2: numpy.ndarray
@@ -30,23 +31,24 @@ class DotPairMeasurement():
@dataclass @dataclass
class DotPairRangeMeasurement(): class DotPairRangeMeasurement:
''' """
Representation of a dot measuring oscillating dipoles. Representation of a dot measuring oscillating dipoles.
Parameters Parameters
---------- ----------
v_low : float v_low : float
The lower range of voltage measured at the dot. The lower range of voltage measured at the dot.
v_high : float v_high : float
The upper range of voltage measured at the dot. The upper range of voltage measured at the dot.
r1 : numpy.ndarray r1 : numpy.ndarray
The position of the first dot. The position of the first dot.
r2 : numpy.ndarray r2 : numpy.ndarray
The position of the second dot. The position of the second dot.
f : float f : float
The measurement frequency. The measurement frequency.
''' """
v_low: float v_low: float
v_high: float v_high: float
r1: numpy.ndarray r1: numpy.ndarray

View File

@@ -9,14 +9,28 @@ DotPairInput = Tuple[numpy.typing.ArrayLike, numpy.typing.ArrayLike, float]
def dot_inputs_to_array(dot_inputs: Sequence[DotInput]) -> numpy.ndarray: def dot_inputs_to_array(dot_inputs: Sequence[DotInput]) -> numpy.ndarray:
return numpy.array([numpy.append(numpy.array(input[0]), input[1]) for input in dot_inputs]) return numpy.array(
[numpy.append(numpy.array(input[0]), input[1]) for input in dot_inputs]
)
def dot_pair_inputs_to_array(pair_inputs: Sequence[DotPairInput]) -> numpy.ndarray: def dot_pair_inputs_to_array(pair_inputs: Sequence[DotPairInput]) -> numpy.ndarray:
return numpy.array([[numpy.append(numpy.array(input[0]), input[2]), numpy.append(numpy.array(input[1]), input[2])] for input in pair_inputs]) return numpy.array(
[
[
numpy.append(numpy.array(input[0]), input[2]),
numpy.append(numpy.array(input[1]), input[2]),
]
for input in pair_inputs
]
)
def dot_range_measurements_low_high_arrays(dot_range_measurements: Union[Sequence[DotRangeMeasurement], Sequence[DotPairRangeMeasurement]]) -> Tuple[numpy.ndarray, numpy.ndarray]: def dot_range_measurements_low_high_arrays(
dot_range_measurements: Union[
Sequence[DotRangeMeasurement], Sequence[DotPairRangeMeasurement]
]
) -> Tuple[numpy.ndarray, numpy.ndarray]:
lows = [measurement.v_low for measurement in dot_range_measurements] lows = [measurement.v_low for measurement in dot_range_measurements]
highs = [measurement.v_high for measurement in dot_range_measurements] highs = [measurement.v_high for measurement in dot_range_measurements]
return (numpy.array(lows), numpy.array(highs)) return (numpy.array(lows), numpy.array(highs))

View File

@@ -3,51 +3,55 @@ import numpy
import numpy.typing import numpy.typing
from typing import Sequence, List from typing import Sequence, List
from pdme.measurement.dot_measure import DotMeasurement, DotRangeMeasurement from pdme.measurement.dot_measure import DotMeasurement, DotRangeMeasurement
from pdme.measurement.dot_pair_measure import DotPairMeasurement, DotPairRangeMeasurement from pdme.measurement.dot_pair_measure import (
DotPairMeasurement,
DotPairRangeMeasurement,
)
from pdme.measurement.input_types import DotInput, DotPairInput from pdme.measurement.input_types import DotInput, DotPairInput
@dataclass @dataclass
class OscillatingDipole(): class OscillatingDipole:
''' """
Representation of an oscillating dipole, either known or guessed. Representation of an oscillating dipole, either known or guessed.
Parameters Parameters
---------- ----------
p : numpy.ndarray p : numpy.ndarray
The oscillating dipole moment, with overall sign arbitrary. The oscillating dipole moment, with overall sign arbitrary.
s : numpy.ndarray s : numpy.ndarray
The position of the dipole. The position of the dipole.
w : float w : float
The oscillation frequency. The oscillation frequency.
''' """
p: numpy.ndarray p: numpy.ndarray
s: numpy.ndarray s: numpy.ndarray
w: float w: float
def __post_init__(self) -> None: def __post_init__(self) -> None:
''' """
Coerce the inputs into numpy arrays. Coerce the inputs into numpy arrays.
''' """
self.p = numpy.array(self.p) self.p = numpy.array(self.p)
self.s = numpy.array(self.s) self.s = numpy.array(self.s)
def s_at_position(self, r: numpy.ndarray, f: float) -> float: def s_at_position(self, r: numpy.ndarray, f: float) -> float:
''' """
Returns the noise potential at a point r, at some frequency f. Returns the noise potential at a point r, at some frequency f.
Parameters Parameters
---------- ----------
r : numpy.ndarray r : numpy.ndarray
The position of the dot. The position of the dot.
f : float f : float
The dot frequency to sample. The dot frequency to sample.
''' """
return (self._alpha(r))**2 * self._b(f) return (self._alpha(r)) ** 2 * self._b(f)
def _alpha(self, r: numpy.ndarray) -> float: def _alpha(self, r: numpy.ndarray) -> float:
diff = r - self.s diff = r - self.s
return self.p.dot(diff) / (numpy.linalg.norm(diff)**3) return self.p.dot(diff) / (numpy.linalg.norm(diff) ** 3)
def _b(self, f: float) -> float: def _b(self, f: float) -> float:
return (1 / numpy.pi) * (self.w / (f**2 + self.w**2)) return (1 / numpy.pi) * (self.w / (f**2 + self.w**2))
@@ -59,65 +63,114 @@ class OscillatingDipole():
return numpy.concatenate([self.p, self.s, numpy.array([self.w])]) return numpy.concatenate([self.p, self.s, numpy.array([self.w])])
class OscillatingDipoleArrangement(): class OscillatingDipoleArrangement:
''' """
A collection of oscillating dipoles, which we are interested in being able to characterise. A collection of oscillating dipoles, which we are interested in being able to characterise.
Parameters Parameters
-------- --------
dipoles : Sequence[OscillatingDipole] dipoles : Sequence[OscillatingDipole]
''' """
def __init__(self, dipoles: Sequence[OscillatingDipole]): def __init__(self, dipoles: Sequence[OscillatingDipole]):
self.dipoles = dipoles self.dipoles = dipoles
def get_dot_measurement(self, dot_input: DotInput) -> DotMeasurement: def get_dot_measurement(self, dot_input: DotInput) -> DotMeasurement:
r = numpy.array(dot_input[0]) r = numpy.array(dot_input[0])
f = dot_input[1] f = dot_input[1]
return DotMeasurement(sum([dipole.s_at_position(r, f) for dipole in self.dipoles]), r, f) return DotMeasurement(
sum([dipole.s_at_position(r, f) for dipole in self.dipoles]), r, f
)
def get_dot_pair_measurement(self, dot_pair_input: DotPairInput) -> DotPairMeasurement: def get_dot_pair_measurement(
self, dot_pair_input: DotPairInput
) -> DotPairMeasurement:
r1 = numpy.array(dot_pair_input[0]) r1 = numpy.array(dot_pair_input[0])
r2 = numpy.array(dot_pair_input[1]) r2 = numpy.array(dot_pair_input[1])
f = dot_pair_input[2] f = dot_pair_input[2]
return DotPairMeasurement(sum([dipole.s_for_dot_pair(r1, r2, f) for dipole in self.dipoles]), r1, r2, f) return DotPairMeasurement(
sum([dipole.s_for_dot_pair(r1, r2, f) for dipole in self.dipoles]),
r1,
r2,
f,
)
def get_dot_measurements(self, dot_inputs: Sequence[DotInput]) -> List[DotMeasurement]: def get_dot_measurements(
''' self, dot_inputs: Sequence[DotInput]
) -> List[DotMeasurement]:
"""
For a series of points, each with three coordinates and a frequency, return a list of the corresponding DotMeasurements. For a series of points, each with three coordinates and a frequency, return a list of the corresponding DotMeasurements.
''' """
return [self.get_dot_measurement(dot_input) for dot_input in dot_inputs] return [self.get_dot_measurement(dot_input) for dot_input in dot_inputs]
def get_dot_pair_measurements(self, dot_pair_inputs: Sequence[DotPairInput]) -> List[DotPairMeasurement]: def get_dot_pair_measurements(
''' self, dot_pair_inputs: Sequence[DotPairInput]
) -> List[DotPairMeasurement]:
"""
For a series of pairs of points, each with three coordinates and a frequency, return a list of the corresponding DotPairMeasurements. For a series of pairs of points, each with three coordinates and a frequency, return a list of the corresponding DotPairMeasurements.
''' """
return [self.get_dot_pair_measurement(dot_pair_input) for dot_pair_input in dot_pair_inputs] return [
self.get_dot_pair_measurement(dot_pair_input)
for dot_pair_input in dot_pair_inputs
]
def get_percent_range_dot_measurement(self, dot_input: DotInput, low_percent: float, high_percent: float) -> DotRangeMeasurement: def get_percent_range_dot_measurement(
self, dot_input: DotInput, low_percent: float, high_percent: float
) -> DotRangeMeasurement:
r = numpy.array(dot_input[0]) r = numpy.array(dot_input[0])
f = dot_input[1] f = dot_input[1]
return DotRangeMeasurement(low_percent * sum([dipole.s_at_position(r, f) for dipole in self.dipoles]), high_percent * sum([dipole.s_at_position(r, f) for dipole in self.dipoles]), r, f) return DotRangeMeasurement(
low_percent * sum([dipole.s_at_position(r, f) for dipole in self.dipoles]),
high_percent * sum([dipole.s_at_position(r, f) for dipole in self.dipoles]),
r,
f,
)
def get_percent_range_dot_measurements(self, dot_inputs: Sequence[DotInput], low_percent: float, high_percent: float) -> List[DotRangeMeasurement]: def get_percent_range_dot_measurements(
''' self, dot_inputs: Sequence[DotInput], low_percent: float, high_percent: float
) -> List[DotRangeMeasurement]:
"""
For a series of pairs of points, each with three coordinates and a frequency, and also a lower error range and upper error range, return a list of the corresponding DotPairRangeMeasurements. For a series of pairs of points, each with three coordinates and a frequency, and also a lower error range and upper error range, return a list of the corresponding DotPairRangeMeasurements.
''' """
return [self.get_percent_range_dot_measurement(dot_input, low_percent, high_percent) for dot_input in dot_inputs] return [
self.get_percent_range_dot_measurement(dot_input, low_percent, high_percent)
for dot_input in dot_inputs
]
def get_percent_range_dot_pair_measurement(self, pair_input: DotPairInput, low_percent: float, high_percent: float) -> DotPairRangeMeasurement: def get_percent_range_dot_pair_measurement(
self, pair_input: DotPairInput, low_percent: float, high_percent: float
) -> DotPairRangeMeasurement:
r1 = numpy.array(pair_input[0]) r1 = numpy.array(pair_input[0])
r2 = numpy.array(pair_input[1]) r2 = numpy.array(pair_input[1])
f = pair_input[2] f = pair_input[2]
return DotPairRangeMeasurement(low_percent * sum([dipole.s_for_dot_pair(r1, r2, f) for dipole in self.dipoles]), high_percent * sum([dipole.s_for_dot_pair(r1, r2, f) for dipole in self.dipoles]), r1, r2, f) return DotPairRangeMeasurement(
low_percent
* sum([dipole.s_for_dot_pair(r1, r2, f) for dipole in self.dipoles]),
high_percent
* sum([dipole.s_for_dot_pair(r1, r2, f) for dipole in self.dipoles]),
r1,
r2,
f,
)
def get_percent_range_dot_pair_measurements(self, pair_inputs: Sequence[DotPairInput], low_percent: float, high_percent: float) -> List[DotPairRangeMeasurement]: def get_percent_range_dot_pair_measurements(
''' self,
pair_inputs: Sequence[DotPairInput],
low_percent: float,
high_percent: float,
) -> List[DotPairRangeMeasurement]:
"""
For a series of pairs of points, each with three coordinates and a frequency, and also a lower error range and upper error range, return a list of the corresponding DotPairRangeMeasurements. For a series of pairs of points, each with three coordinates and a frequency, and also a lower error range and upper error range, return a list of the corresponding DotPairRangeMeasurements.
''' """
return [self.get_percent_range_dot_pair_measurement(pair_input, low_percent, high_percent) for pair_input in pair_inputs] return [
self.get_percent_range_dot_pair_measurement(
pair_input, low_percent, high_percent
)
for pair_input in pair_inputs
]
def to_numpy_array(self) -> numpy.ndarray: def to_numpy_array(self) -> numpy.ndarray:
''' """
Returns a numpy array with the canonical representation of each dipole in a nx7 numpy array. Returns a numpy array with the canonical representation of each dipole in a nx7 numpy array.
''' """
return numpy.array([dipole.to_flat_array() for dipole in self.dipoles]) return numpy.array([dipole.to_flat_array() for dipole in self.dipoles])

View File

@@ -1,3 +1,3 @@
from importlib.metadata import version from importlib.metadata import version
__version__ = version('pdme') __version__ = version("pdme")

View File

@@ -4,4 +4,11 @@ from pdme.model.unrestricted_model import UnrestrictedModel
from pdme.model.fixed_dipole_model import FixedDipoleModel from pdme.model.fixed_dipole_model import FixedDipoleModel
from pdme.model.fixed_magnitude_model import FixedMagnitudeModel from pdme.model.fixed_magnitude_model import FixedMagnitudeModel
__all__ = ["Model", "Discretisation", "FixedZPlaneModel", "UnrestrictedModel", "FixedDipoleModel", "FixedMagnitudeModel"] __all__ = [
"Model",
"Discretisation",
"FixedZPlaneModel",
"UnrestrictedModel",
"FixedDipoleModel",
"FixedMagnitudeModel",
]

View File

@@ -4,21 +4,36 @@ from dataclasses import dataclass
from typing import Sequence, Tuple from typing import Sequence, Tuple
import scipy.optimize import scipy.optimize
from pdme.model.model import Model, Discretisation from pdme.model.model import Model, Discretisation
from pdme.measurement import DotMeasurement, OscillatingDipoleArrangement, OscillatingDipole from pdme.measurement import (
DotMeasurement,
OscillatingDipoleArrangement,
OscillatingDipole,
)
class FixedDipoleModel(Model): class FixedDipoleModel(Model):
''' """
Model of oscillating dipole with a fixed dipole moment. Model of oscillating dipole with a fixed dipole moment.
Parameters Parameters
---------- ----------
p : numpy.ndarray p : numpy.ndarray
The fixed dipole moment. The fixed dipole moment.
n : int n : int
The number of dipoles to assume. The number of dipoles to assume.
''' """
def __init__(self, xmin: float, xmax: float, ymin: float, ymax: float, zmin: float, zmax: float, p: numpy.ndarray, n: int) -> None:
def __init__(
self,
xmin: float,
xmax: float,
ymin: float,
ymax: float,
zmin: float,
zmax: float,
p: numpy.ndarray,
n: int,
) -> None:
self.xmin = xmin self.xmin = xmin
self.xmax = xmax self.xmax = xmax
self.ymin = ymin self.ymin = ymin
@@ -30,12 +45,20 @@ class FixedDipoleModel(Model):
self.rng = numpy.random.default_rng() self.rng = numpy.random.default_rng()
def __repr__(self) -> str: def __repr__(self) -> str:
return f'FixedDipoleModel({self.xmin}, {self.xmax}, {self.ymin}, {self.ymax}, {self.zmin}, {self.zmax}, {self.p}, {self.n()})' return f"FixedDipoleModel({self.xmin}, {self.xmax}, {self.ymin}, {self.ymax}, {self.zmin}, {self.zmax}, {self.p}, {self.n()})"
# TODO: this signature doesn't make sense. # TODO: this signature doesn't make sense.
def get_dipoles(self, frequency: float) -> OscillatingDipoleArrangement: def get_dipoles(self, frequency: float) -> OscillatingDipoleArrangement:
s_pts = numpy.array((self.rng.uniform(self.xmin, self.xmax), self.rng.uniform(self.ymin, self.ymax), self.rng.uniform(self.zmin, self.zmax))) s_pts = numpy.array(
return OscillatingDipoleArrangement([OscillatingDipole(self.p, s_pts, frequency)]) (
self.rng.uniform(self.xmin, self.xmax),
self.rng.uniform(self.ymin, self.ymax),
self.rng.uniform(self.zmin, self.zmax),
)
)
return OscillatingDipoleArrangement(
[OscillatingDipole(self.p, s_pts, frequency)]
)
def solution_single_dipole(self, pt: numpy.ndarray) -> OscillatingDipole: def solution_single_dipole(self, pt: numpy.ndarray) -> OscillatingDipole:
# assume length is 4. # assume length is 4.
@@ -44,10 +67,10 @@ class FixedDipoleModel(Model):
return OscillatingDipole(self.p, s, w) return OscillatingDipole(self.p, s, w)
def point_length(self) -> int: def point_length(self) -> int:
''' """
Dipole is constrained magnitude, but free orientation. Dipole is constrained magnitude, but free orientation.
Six degrees of freedom: (sx, sy, sz, w). Six degrees of freedom: (sx, sy, sz, w).
''' """
return 4 return 4
def n(self) -> int: def n(self) -> int:
@@ -58,45 +81,56 @@ class FixedDipoleModel(Model):
w = pt[3] w = pt[3]
diff = dot.r - s diff = dot.r - s
alpha = self.p.dot(diff) / (numpy.linalg.norm(diff)**3) alpha = self.p.dot(diff) / (numpy.linalg.norm(diff) ** 3)
b = (1 / numpy.pi) * (w / (w**2 + dot.f**2)) b = (1 / numpy.pi) * (w / (w**2 + dot.f**2))
return alpha**2 * b return alpha**2 * b
def jac_for_point_at_dot(self, dot: DotMeasurement, pt: numpy.ndarray) -> numpy.ndarray: def jac_for_point_at_dot(
self, dot: DotMeasurement, pt: numpy.ndarray
) -> numpy.ndarray:
s = pt[0:3] s = pt[0:3]
w = pt[3] w = pt[3]
diff = dot.r - s diff = dot.r - s
alpha = self.p.dot(diff) / (numpy.linalg.norm(diff)**3) alpha = self.p.dot(diff) / (numpy.linalg.norm(diff) ** 3)
b = (1 / numpy.pi) * (w / (w**2 + dot.f**2)) b = (1 / numpy.pi) * (w / (w**2 + dot.f**2))
r_divs = (-self.p / (numpy.linalg.norm(diff)**3) + 3 * self.p.dot(diff) * diff / (numpy.linalg.norm(diff)**5)) * 2 * alpha * b r_divs = (
(
-self.p / (numpy.linalg.norm(diff) ** 3)
+ 3 * self.p.dot(diff) * diff / (numpy.linalg.norm(diff) ** 5)
)
* 2
* alpha
* b
)
f2 = dot.f**2 f2 = dot.f**2
w2 = w**2 w2 = w**2
w_div = alpha**2 * (1 / numpy.pi) * ((f2 - w2) / ((f2 + w2)**2)) w_div = alpha**2 * (1 / numpy.pi) * ((f2 - w2) / ((f2 + w2) ** 2))
return numpy.concatenate((r_divs, w_div), axis=None) return numpy.concatenate((r_divs, w_div), axis=None)
@dataclass @dataclass
class FixedDipoleDiscretisation(Discretisation): class FixedDipoleDiscretisation(Discretisation):
''' """
Representation of a discretisation of a FixedDipoleDiscretisation. Representation of a discretisation of a FixedDipoleDiscretisation.
Also captures a rough maximum value of dipole. Also captures a rough maximum value of dipole.
Parameters Parameters
---------- ----------
model : FixedDipoleModel model : FixedDipoleModel
The parent model of the discretisation. The parent model of the discretisation.
num_x : int num_x : int
The number of partitions of the x axis. The number of partitions of the x axis.
num_y : int num_y : int
The number of partitions of the y axis. The number of partitions of the y axis.
num_z : int num_z : int
The number of partitions of the z axis. The number of partitions of the z axis.
''' """
model: FixedDipoleModel model: FixedDipoleModel
num_x: int num_x: int
num_y: int num_y: int
@@ -115,22 +149,32 @@ class FixedDipoleDiscretisation(Discretisation):
# We want to keep w unbounded, restrict sx, sy, sz, px and py based on step. # We want to keep w unbounded, restrict sx, sy, sz, px and py based on step.
return ( return (
[ [
xi * self.x_step + self.model.xmin, yi * self.y_step + self.model.ymin, zi * self.z_step + self.model.zmin, xi * self.x_step + self.model.xmin,
-numpy.inf yi * self.y_step + self.model.ymin,
zi * self.z_step + self.model.zmin,
-numpy.inf,
], ],
[ [
(xi + 1) * self.x_step + self.model.xmin, (yi + 1) * self.y_step + self.model.ymin, (zi + 1) * self.z_step + self.model.zmin, (xi + 1) * self.x_step + self.model.xmin,
numpy.inf (yi + 1) * self.y_step + self.model.ymin,
] (zi + 1) * self.z_step + self.model.zmin,
numpy.inf,
],
) )
def all_indices(self) -> numpy.ndindex: def all_indices(self) -> numpy.ndindex:
# see https://github.com/numpy/numpy/issues/20706 for why this is a mypy problem. # see https://github.com/numpy/numpy/issues/20706 for why this is a mypy problem.
return numpy.ndindex((self.num_x, self.num_y, self.num_z)) # type:ignore return numpy.ndindex((self.num_x, self.num_y, self.num_z)) # type:ignore
def solve_for_index(self, dots: Sequence[DotMeasurement], index: Tuple[float, ...]) -> scipy.optimize.OptimizeResult: def solve_for_index(
self, dots: Sequence[DotMeasurement], index: Tuple[float, ...]
) -> scipy.optimize.OptimizeResult:
bounds = self.bounds(index) bounds = self.bounds(index)
sx_mean = (bounds[0][0] + bounds[1][0]) / 2 sx_mean = (bounds[0][0] + bounds[1][0]) / 2
sy_mean = (bounds[0][1] + bounds[1][1]) / 2 sy_mean = (bounds[0][1] + bounds[1][1]) / 2
sz_mean = (bounds[0][2] + bounds[1][2]) / 2 sz_mean = (bounds[0][2] + bounds[1][2]) / 2
return self.model.solve(dots, initial_pt=numpy.array([sx_mean, sy_mean, sz_mean, .1]), bounds=bounds) return self.model.solve(
dots,
initial_pt=numpy.array([sx_mean, sy_mean, sz_mean, 0.1]),
bounds=bounds,
)

View File

@@ -4,21 +4,36 @@ from dataclasses import dataclass
from typing import Sequence, Tuple from typing import Sequence, Tuple
import scipy.optimize import scipy.optimize
from pdme.model.model import Model, Discretisation from pdme.model.model import Model, Discretisation
from pdme.measurement import DotMeasurement, OscillatingDipole, OscillatingDipoleArrangement from pdme.measurement import (
DotMeasurement,
OscillatingDipole,
OscillatingDipoleArrangement,
)
class FixedMagnitudeModel(Model): class FixedMagnitudeModel(Model):
''' """
Model of oscillating dipole with a fixed magnitude, but free rotation. Model of oscillating dipole with a fixed magnitude, but free rotation.
Parameters Parameters
---------- ----------
pfixed : float pfixed : float
The fixed dipole magnitude. The fixed dipole magnitude.
n : int n : int
The number of dipoles to assume. The number of dipoles to assume.
''' """
def __init__(self, xmin: float, xmax: float, ymin: float, ymax: float, zmin: float, zmax: float, pfixed: float, n: int) -> None:
def __init__(
self,
xmin: float,
xmax: float,
ymin: float,
ymax: float,
zmin: float,
zmax: float,
pfixed: float,
n: int,
) -> None:
self.xmin = xmin self.xmin = xmin
self.xmax = xmax self.xmax = xmax
self.ymin = ymin self.ymin = ymin
@@ -30,7 +45,7 @@ class FixedMagnitudeModel(Model):
self.rng = numpy.random.default_rng() self.rng = numpy.random.default_rng()
def __repr__(self) -> str: def __repr__(self) -> str:
return f'FixedMagnitudeModel({self.xmin}, {self.xmax}, {self.ymin}, {self.ymax}, {self.zmin}, {self.zmax}, {self.pfixed}, {self.n()})' 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: def solution_single_dipole(self, pt: numpy.ndarray) -> OscillatingDipole:
# assume length is 6, who needs error checking. # assume length is 6, who needs error checking.
@@ -39,18 +54,20 @@ class FixedMagnitudeModel(Model):
s = pt[2:5] s = pt[2:5]
w = pt[5] w = pt[5]
p = numpy.array([ 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.sin(p_theta) * numpy.cos(p_phi),
self.pfixed * numpy.cos(p_theta) self.pfixed * numpy.sin(p_theta) * numpy.sin(p_phi),
]) self.pfixed * numpy.cos(p_theta),
]
)
return OscillatingDipole(p, s, w) return OscillatingDipole(p, s, w)
def point_length(self) -> int: def point_length(self) -> int:
''' """
Dipole is constrained magnitude, but free orientation. Dipole is constrained magnitude, but free orientation.
Six degrees of freedom: (p_theta, p_phi, sx, sy, sz, w). Six degrees of freedom: (p_theta, p_phi, sx, sy, sz, w).
''' """
return 6 return 6
def get_dipoles(self, frequency: float) -> OscillatingDipoleArrangement: def get_dipoles(self, frequency: float) -> OscillatingDipoleArrangement:
@@ -59,8 +76,16 @@ class FixedMagnitudeModel(Model):
px = self.pfixed * numpy.sin(theta) * numpy.cos(phi) px = self.pfixed * numpy.sin(theta) * numpy.cos(phi)
py = self.pfixed * numpy.sin(theta) * numpy.sin(phi) py = self.pfixed * numpy.sin(theta) * numpy.sin(phi)
pz = self.pfixed * numpy.cos(theta) pz = self.pfixed * numpy.cos(theta)
s_pts = numpy.array((self.rng.uniform(self.xmin, self.xmax), self.rng.uniform(self.ymin, self.ymax), self.rng.uniform(self.zmin, self.zmax))) s_pts = numpy.array(
return OscillatingDipoleArrangement([OscillatingDipole(numpy.array([px, py, pz]), s_pts, frequency)]) (
self.rng.uniform(self.xmin, self.xmax),
self.rng.uniform(self.ymin, self.ymax),
self.rng.uniform(self.zmin, self.zmax),
)
)
return OscillatingDipoleArrangement(
[OscillatingDipole(numpy.array([px, py, pz]), s_pts, frequency)]
)
def get_n_single_dipoles(self, n: int, max_frequency: float) -> numpy.ndarray: def get_n_single_dipoles(self, n: int, max_frequency: float) -> numpy.ndarray:
# psw # psw
@@ -88,29 +113,35 @@ class FixedMagnitudeModel(Model):
s = pt[2:5] s = pt[2:5]
w = pt[5] w = pt[5]
p = numpy.array([ 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.sin(p_theta) * numpy.cos(p_phi),
self.pfixed * numpy.cos(p_theta) self.pfixed * numpy.sin(p_theta) * numpy.sin(p_phi),
]) self.pfixed * numpy.cos(p_theta),
]
)
diff = dot.r - s diff = dot.r - s
alpha = p.dot(diff) / (numpy.linalg.norm(diff)**3) alpha = p.dot(diff) / (numpy.linalg.norm(diff) ** 3)
b = (1 / numpy.pi) * (w / (w**2 + dot.f**2)) b = (1 / numpy.pi) * (w / (w**2 + dot.f**2))
return alpha**2 * b return alpha**2 * b
def jac_for_point_at_dot(self, dot: DotMeasurement, pt: numpy.ndarray) -> numpy.ndarray: def jac_for_point_at_dot(
self, dot: DotMeasurement, pt: numpy.ndarray
) -> numpy.ndarray:
p_theta = pt[0] p_theta = pt[0]
p_phi = pt[1] p_phi = pt[1]
s = pt[2:5] s = pt[2:5]
w = pt[5] w = pt[5]
p = numpy.array([ 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.sin(p_theta) * numpy.cos(p_phi),
self.pfixed * numpy.cos(p_theta) self.pfixed * numpy.sin(p_theta) * numpy.sin(p_phi),
]) self.pfixed * numpy.cos(p_theta),
]
)
diff = dot.r - s diff = dot.r - s
alpha = p.dot(diff) / (numpy.linalg.norm(diff)**3) alpha = p.dot(diff) / (numpy.linalg.norm(diff) ** 3)
b = (1 / numpy.pi) * (w / (w**2 + dot.f**2)) b = (1 / numpy.pi) * (w / (w**2 + dot.f**2))
theta_div_middle = self.pfixed * ( theta_div_middle = self.pfixed * (
@@ -118,45 +149,54 @@ class FixedMagnitudeModel(Model):
+ diff[1] * numpy.sin(p_phi) * numpy.cos(p_theta) + diff[1] * numpy.sin(p_phi) * numpy.cos(p_theta)
- diff[2] * numpy.sin(p_theta) - diff[2] * numpy.sin(p_theta)
) )
theta_div = 2 * alpha * (theta_div_middle) / (numpy.linalg.norm(diff)**3) * b theta_div = 2 * alpha * (theta_div_middle) / (numpy.linalg.norm(diff) ** 3) * b
phi_div_middle = self.pfixed * ( phi_div_middle = self.pfixed * (
diff[1] * numpy.sin(p_theta) * numpy.cos(p_phi) diff[1] * numpy.sin(p_theta) * numpy.cos(p_phi)
- diff[0] * numpy.sin(p_theta) * numpy.sin(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 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 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 f2 = dot.f**2
w2 = w**2 w2 = w**2
w_div = alpha**2 * (1 / numpy.pi) * ((f2 - w2) / ((f2 + w2)**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) return numpy.concatenate((theta_div, phi_div, r_divs, w_div), axis=None)
@dataclass @dataclass
class FixedMagnitudeDiscretisation(Discretisation): class FixedMagnitudeDiscretisation(Discretisation):
''' """
Representation of a discretisation of a FixedMagnitudeDiscretisation. Representation of a discretisation of a FixedMagnitudeDiscretisation.
Also captures a rough maximum value of dipole. Also captures a rough maximum value of dipole.
Parameters Parameters
---------- ----------
model : FixedMagnitudeModel model : FixedMagnitudeModel
The parent model of the discretisation. The parent model of the discretisation.
num_ptheta: int num_ptheta: int
The number of partitions of ptheta. The number of partitions of ptheta.
num_pphi: int num_pphi: int
The number of partitions of pphi. The number of partitions of pphi.
num_x : int num_x : int
The number of partitions of the x axis. The number of partitions of the x axis.
num_y : int num_y : int
The number of partitions of the y axis. The number of partitions of the y axis.
num_z : int num_z : int
The number of partitions of the z axis. The number of partitions of the z axis.
''' """
model: FixedMagnitudeModel model: FixedMagnitudeModel
num_ptheta: int num_ptheta: int
num_pphi: int num_pphi: int
@@ -179,15 +219,21 @@ class FixedMagnitudeDiscretisation(Discretisation):
# We want to keep w unbounded, restrict sx, sy, sz, px and py based on step. # We want to keep w unbounded, restrict sx, sy, sz, px and py based on step.
return ( return (
[ [
numpy.arccos(1 - pthetai * self.h_step), pphii * self.phi_step, numpy.arccos(1 - pthetai * self.h_step),
xi * self.x_step + self.model.xmin, yi * self.y_step + self.model.ymin, zi * self.z_step + self.model.zmin, pphii * self.phi_step,
-numpy.inf xi * self.x_step + self.model.xmin,
yi * self.y_step + self.model.ymin,
zi * self.z_step + self.model.zmin,
-numpy.inf,
], ],
[ [
numpy.arccos(1 - (pthetai + 1) * self.h_step), (pphii + 1) * self.phi_step, numpy.arccos(1 - (pthetai + 1) * self.h_step),
(xi + 1) * self.x_step + self.model.xmin, (yi + 1) * self.y_step + self.model.ymin, (zi + 1) * self.z_step + self.model.zmin, (pphii + 1) * self.phi_step,
numpy.inf (xi + 1) * self.x_step + self.model.xmin,
] (yi + 1) * self.y_step + self.model.ymin,
(zi + 1) * self.z_step + self.model.zmin,
numpy.inf,
],
) )
def get_model(self) -> Model: def get_model(self) -> Model:
@@ -195,13 +241,23 @@ class FixedMagnitudeDiscretisation(Discretisation):
def all_indices(self) -> numpy.ndindex: def all_indices(self) -> numpy.ndindex:
# see https://github.com/numpy/numpy/issues/20706 for why this is a mypy problem. # see https://github.com/numpy/numpy/issues/20706 for why this is a mypy problem.
return numpy.ndindex((self.num_ptheta, self.num_pphi, self.num_x, self.num_y, self.num_z)) # type:ignore return numpy.ndindex(
(self.num_ptheta, self.num_pphi, self.num_x, self.num_y, self.num_z)
) # type:ignore
def solve_for_index(self, dots: Sequence[DotMeasurement], index: Tuple[float, ...]) -> scipy.optimize.OptimizeResult: def solve_for_index(
self, dots: Sequence[DotMeasurement], index: Tuple[float, ...]
) -> scipy.optimize.OptimizeResult:
bounds = self.bounds(index) bounds = self.bounds(index)
ptheta_mean = (bounds[0][0] + bounds[1][0]) / 2 ptheta_mean = (bounds[0][0] + bounds[1][0]) / 2
pphi_mean = (bounds[0][1] + bounds[1][1]) / 2 pphi_mean = (bounds[0][1] + bounds[1][1]) / 2
sx_mean = (bounds[0][2] + bounds[1][2]) / 2 sx_mean = (bounds[0][2] + bounds[1][2]) / 2
sy_mean = (bounds[0][3] + bounds[1][3]) / 2 sy_mean = (bounds[0][3] + bounds[1][3]) / 2
sz_mean = (bounds[0][4] + bounds[1][4]) / 2 sz_mean = (bounds[0][4] + bounds[1][4]) / 2
return self.model.solve(dots, initial_pt=numpy.array([ptheta_mean, pphi_mean, sx_mean, sy_mean, sz_mean, .1]), bounds=bounds) return self.model.solve(
dots,
initial_pt=numpy.array(
[ptheta_mean, pphi_mean, sx_mean, sy_mean, sz_mean, 0.1]
),
bounds=bounds,
)

View File

@@ -8,26 +8,29 @@ from pdme.measurement import DotMeasurement
class FixedZPlaneModel(Model): class FixedZPlaneModel(Model):
''' """
Model of oscillating dipoles constrained to lie within a plane. Model of oscillating dipoles constrained to lie within a plane.
Additionally, each dipole is assumed to be orientated in the plus or minus z direction. Additionally, each dipole is assumed to be orientated in the plus or minus z direction.
Parameters Parameters
---------- ----------
z : float z : float
The z position of the plane where dipoles are constrained to lie. The z position of the plane where dipoles are constrained to lie.
xmin : float xmin : float
The minimum x value for dipoles. The minimum x value for dipoles.
xmax : float xmax : float
The maximum x value for dipoles. The maximum x value for dipoles.
ymin : float ymin : float
The minimum y value for dipoles. The minimum y value for dipoles.
ymax : float ymax : float
The maximum y value for dipoles. The maximum y value for dipoles.
n : int n : int
The number of dipoles to assume. The number of dipoles to assume.
''' """
def __init__(self, z: float, xmin: float, xmax: float, ymin: float, ymax: float, n: int) -> None:
def __init__(
self, z: float, xmin: float, xmax: float, ymin: float, ymax: float, n: int
) -> None:
self.z = z self.z = z
self.xmin = xmin self.xmin = xmin
self.xmax = xmax self.xmax = xmax
@@ -36,13 +39,13 @@ class FixedZPlaneModel(Model):
self._n = n self._n = n
def __repr__(self) -> str: def __repr__(self) -> str:
return f'FixedZPlaneModel({self.z}, {self.xmin}, {self.xmax}, {self.ymin}, {self.ymax}, {self.n()})' return f"FixedZPlaneModel({self.z}, {self.xmin}, {self.xmax}, {self.ymin}, {self.ymax}, {self.n()})"
def point_length(self) -> int: def point_length(self) -> int:
''' """
Dipole is constrained in this model to have (px, py, pz) = (0, 0, pz) and (sx, sy, sz) = (sx, sy, self.z). Dipole is constrained in this model to have (px, py, pz) = (0, 0, pz) and (sx, sy, sz) = (sx, sy, self.z).
With some frequency w, there are four degrees of freedom: (pz, sx, sy, w). With some frequency w, there are four degrees of freedom: (pz, sx, sy, w).
''' """
return 4 return 4
def n(self) -> int: def n(self) -> int:
@@ -54,48 +57,61 @@ class FixedZPlaneModel(Model):
w = pt[3] w = pt[3]
diff = dot.r - s diff = dot.r - s
alpha = p.dot(diff) / (numpy.linalg.norm(diff)**3) alpha = p.dot(diff) / (numpy.linalg.norm(diff) ** 3)
b = (1 / numpy.pi) * (w / (w**2 + dot.f**2)) b = (1 / numpy.pi) * (w / (w**2 + dot.f**2))
return alpha**2 * b return alpha**2 * b
def jac_for_point_at_dot(self, dot: DotMeasurement, pt: numpy.ndarray) -> numpy.ndarray: def jac_for_point_at_dot(
self, dot: DotMeasurement, pt: numpy.ndarray
) -> numpy.ndarray:
p = numpy.array([0, 0, pt[0]]) p = numpy.array([0, 0, pt[0]])
s = numpy.array([pt[1], pt[2], self.z]) s = numpy.array([pt[1], pt[2], self.z])
w = pt[3] w = pt[3]
diff = dot.r - s diff = dot.r - s
alpha = p.dot(diff) / (numpy.linalg.norm(diff)**3) alpha = p.dot(diff) / (numpy.linalg.norm(diff) ** 3)
b = (1 / numpy.pi) * (w / (w**2 + dot.f**2)) b = (1 / numpy.pi) * (w / (w**2 + dot.f**2))
p_divs = 2 * alpha * diff[2] / (numpy.linalg.norm(diff)**3) * b # only need the z component. p_divs = (
2 * alpha * diff[2] / (numpy.linalg.norm(diff) ** 3) * b
) # only need the z component.
r_divs = (-p[0:2] / (numpy.linalg.norm(diff)**3) + 3 * p.dot(diff) * diff[0:2] / (numpy.linalg.norm(diff)**5)) * 2 * alpha * b r_divs = (
(
-p[0:2] / (numpy.linalg.norm(diff) ** 3)
+ 3 * p.dot(diff) * diff[0:2] / (numpy.linalg.norm(diff) ** 5)
)
* 2
* alpha
* b
)
f2 = dot.f**2 f2 = dot.f**2
w2 = w**2 w2 = w**2
w_div = alpha**2 * (1 / numpy.pi) * ((f2 - w2) / ((f2 + w2)**2)) w_div = alpha**2 * (1 / numpy.pi) * ((f2 - w2) / ((f2 + w2) ** 2))
return numpy.concatenate((p_divs, r_divs, w_div), axis=None) return numpy.concatenate((p_divs, r_divs, w_div), axis=None)
@dataclass @dataclass
class FixedZPlaneDiscretisation(): class FixedZPlaneDiscretisation:
''' """
Representation of a discretisation of a FixedZPlaneModel. Representation of a discretisation of a FixedZPlaneModel.
Also captures a rough maximum value of dipole. Also captures a rough maximum value of dipole.
Parameters Parameters
---------- ----------
model : FixedZPlaneModel model : FixedZPlaneModel
The parent model of the discretisation. The parent model of the discretisation.
num_pz: int num_pz: int
The number of partitions of pz. The number of partitions of pz.
num_x : int num_x : int
The number of partitions of the x axis. The number of partitions of the x axis.
num_y : int num_y : int
The number of partitions of the y axis. The number of partitions of the y axis.
''' """
model: FixedZPlaneModel model: FixedZPlaneModel
num_pz: int num_pz: int
num_x: int num_x: int
@@ -108,24 +124,42 @@ class FixedZPlaneDiscretisation():
self.x_step = (self.model.xmax - self.model.xmin) / self.num_x self.x_step = (self.model.xmax - self.model.xmin) / self.num_x
self.y_step = (self.model.ymax - self.model.ymin) / self.num_y self.y_step = (self.model.ymax - self.model.ymin) / self.num_y
def bounds(self, index: Tuple[float, float, float]) -> Tuple[numpy.ndarray, numpy.ndarray]: def bounds(
self, index: Tuple[float, float, float]
) -> Tuple[numpy.ndarray, numpy.ndarray]:
pzi, xi, yi = index pzi, xi, yi = index
# For this model, a point is (pz, sx, sy, w). # For this model, a point is (pz, sx, sy, w).
# We want to keep w bounded, and restrict pz, sx and sy based on step. # We want to keep w bounded, and restrict pz, sx and sy based on step.
return ( return (
numpy.array((pzi * self.pz_step - self.max_pz, xi * self.x_step + self.model.xmin, yi * self.y_step + self.model.ymin, -numpy.inf)), numpy.array(
numpy.array(((pzi + 1) * self.pz_step - self.max_pz, (xi + 1) * self.x_step + self.model.xmin, (yi + 1) * self.y_step + self.model.ymin, numpy.inf)) (
pzi * self.pz_step - self.max_pz,
xi * self.x_step + self.model.xmin,
yi * self.y_step + self.model.ymin,
-numpy.inf,
)
),
numpy.array(
(
(pzi + 1) * self.pz_step - self.max_pz,
(xi + 1) * self.x_step + self.model.xmin,
(yi + 1) * self.y_step + self.model.ymin,
numpy.inf,
)
),
) )
def all_indices(self) -> numpy.ndindex: def all_indices(self) -> numpy.ndindex:
# see https://github.com/numpy/numpy/issues/20706 for why this is a mypy problem. # see https://github.com/numpy/numpy/issues/20706 for why this is a mypy problem.
return numpy.ndindex((self.num_pz, self.num_x, self.num_y)) # type:ignore return numpy.ndindex((self.num_pz, self.num_x, self.num_y)) # type:ignore
def solve_for_index(self, dots: Sequence[DotMeasurement], index: Tuple[float, float, float]) -> scipy.optimize.OptimizeResult: def solve_for_index(
self, dots: Sequence[DotMeasurement], index: Tuple[float, float, float]
) -> scipy.optimize.OptimizeResult:
bounds = self.bounds(index) bounds = self.bounds(index)
pz_mean = (bounds[0][0] + bounds[1][0]) / 2 pz_mean = (bounds[0][0] + bounds[1][0]) / 2
sx_mean = (bounds[0][1] + bounds[1][1]) / 2 sx_mean = (bounds[0][1] + bounds[1][1]) / 2
sy_mean = (bounds[0][2] + bounds[1][2]) / 2 sy_mean = (bounds[0][2] + bounds[1][2]) / 2
# I don't care about the typing here at the moment. # I don't care about the typing here at the moment.
return self.model.solve(dots, initial_pt=numpy.array((pz_mean, sx_mean, sy_mean, .1)), bounds=bounds) # type: ignore return self.model.solve(dots, initial_pt=numpy.array((pz_mean, sx_mean, sy_mean, 0.1)), bounds=bounds) # type: ignore

View File

@@ -1,7 +1,11 @@
import numpy import numpy
import scipy.optimize import scipy.optimize
from typing import Callable, Sequence, Tuple, List from typing import Callable, Sequence, Tuple, List
from pdme.measurement import DotMeasurement, OscillatingDipoleArrangement, OscillatingDipole from pdme.measurement import (
DotMeasurement,
OscillatingDipoleArrangement,
OscillatingDipole,
)
import pdme.util import pdme.util
import logging import logging
@@ -9,7 +13,7 @@ import logging
_logger = logging.getLogger(__name__) _logger = logging.getLogger(__name__)
class Model(): class Model:
""" """
Interface for models. Interface for models.
""" """
@@ -34,18 +38,20 @@ class Model():
def solution_as_dipoles(self, pts: numpy.ndarray) -> List[OscillatingDipole]: def solution_as_dipoles(self, pts: numpy.ndarray) -> List[OscillatingDipole]:
pt_length = self.point_length() pt_length = self.point_length()
chunked_pts = [pts[i: i + pt_length] for i in range(0, len(pts), pt_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] return [self.solution_single_dipole(pt) for pt in chunked_pts]
def cost_for_dot(self, dot: DotMeasurement, pts: numpy.ndarray) -> float: def cost_for_dot(self, dot: DotMeasurement, pts: numpy.ndarray) -> float:
# creates numpy.ndarrays in groups of self.point_length(). # creates numpy.ndarrays in groups of self.point_length().
# Will throw problems for irregular points, but that's okay for now. # Will throw problems for irregular points, but that's okay for now.
pt_length = self.point_length() pt_length = self.point_length()
chunked_pts = [pts[i: i + pt_length] for i in range(0, len(pts), pt_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 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]: 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. 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. Default implementation assumes a single dot cost from which to build the list.
@@ -56,7 +62,7 @@ class Model():
Returns Returns
---------- ----------
Returns the model's cost function. Returns the model's cost function.
''' """
_logger.debug(f"Constructing costs for dots: {dots}") _logger.debug(f"Constructing costs for dots: {dots}")
def costs_to_return(pts: numpy.ndarray) -> numpy.ndarray: def costs_to_return(pts: numpy.ndarray) -> numpy.ndarray:
@@ -64,18 +70,24 @@ class Model():
return costs_to_return return costs_to_return
def jac_for_point_at_dot(self, dot: DotMeasurement, pt: numpy.ndarray) -> numpy.ndarray: def jac_for_point_at_dot(
self, dot: DotMeasurement, pt: numpy.ndarray
) -> numpy.ndarray:
raise NotImplementedError raise NotImplementedError
def jac_for_dot(self, dot: DotMeasurement, pts: numpy.ndarray) -> numpy.ndarray: def jac_for_dot(self, dot: DotMeasurement, pts: numpy.ndarray) -> numpy.ndarray:
# creates numpy.ndarrays in groups of self.point_length(). # creates numpy.ndarrays in groups of self.point_length().
# Will throw problems for irregular points, but that's okay for now. # Will throw problems for irregular points, but that's okay for now.
pt_length = self.point_length() pt_length = self.point_length()
chunked_pts = [pts[i: i + pt_length] for i in range(0, len(pts), pt_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]) 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]: 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. 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. Default implementation assumes a single dot jacobian from which to build the list.
@@ -86,33 +98,53 @@ class Model():
Returns Returns
---------- ----------
Returns the model's cost function's Jacobian. Returns the model's cost function's Jacobian.
''' """
def jac_to_return(pts: numpy.ndarray) -> numpy.ndarray: def jac_to_return(pts: numpy.ndarray) -> numpy.ndarray:
return numpy.array([self.jac_for_dot(dot, pts) for dot in dots]) return numpy.array([self.jac_for_dot(dot, pts) for dot in dots])
return jac_to_return return jac_to_return
def solve(self, dots: Sequence[DotMeasurement], initial_pt: numpy.ndarray = None, bounds=(-numpy.inf, numpy.inf)) -> scipy.optimize.OptimizeResult: def solve(
self,
dots: Sequence[DotMeasurement],
initial_pt: numpy.ndarray = None,
bounds=(-numpy.inf, numpy.inf),
) -> scipy.optimize.OptimizeResult:
if initial_pt is None: if initial_pt is None:
initial = numpy.tile(.1, self.n() * self.point_length()) initial = numpy.tile(0.1, self.n() * self.point_length())
else: else:
if len(initial_pt) != self.point_length(): 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()}") 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()) 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 = scipy.optimize.least_squares(
result.normalised_x = pdme.util.normalise_point_list(result.x, self.point_length()) 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 return result
class Discretisation(): class Discretisation:
def bounds(self, index: Tuple[float, ...]) -> Tuple: def bounds(self, index: Tuple[float, ...]) -> Tuple:
raise NotImplementedError raise NotImplementedError
def all_indices(self) -> numpy.ndindex: def all_indices(self) -> numpy.ndindex:
raise NotImplementedError raise NotImplementedError
def solve_for_index(self, dots: Sequence[DotMeasurement], index: Tuple) -> scipy.optimize.OptimizeResult: def solve_for_index(
self, dots: Sequence[DotMeasurement], index: Tuple
) -> scipy.optimize.OptimizeResult:
raise NotImplementedError raise NotImplementedError
def get_model(self) -> Model: def get_model(self) -> Model:

View File

@@ -3,20 +3,35 @@ from dataclasses import dataclass
from typing import Sequence, Tuple from typing import Sequence, Tuple
import scipy.optimize import scipy.optimize
from pdme.model.model import Model, Discretisation from pdme.model.model import Model, Discretisation
from pdme.measurement import DotMeasurement, OscillatingDipoleArrangement, OscillatingDipole from pdme.measurement import (
DotMeasurement,
OscillatingDipoleArrangement,
OscillatingDipole,
)
class UnrestrictedModel(Model): class UnrestrictedModel(Model):
''' """
Model of oscillating dipoles with no restrictions. Model of oscillating dipoles with no restrictions.
Additionally, each dipole is assumed to be orientated in the plus or minus z direction. Additionally, each dipole is assumed to be orientated in the plus or minus z direction.
Parameters Parameters
---------- ----------
n : int n : int
The number of dipoles to assume. The number of dipoles to assume.
''' """
def __init__(self, xmin: float, xmax: float, ymin: float, ymax: float, zmin: float, zmax: float, max_p: float, n: int) -> None:
def __init__(
self,
xmin: float,
xmax: float,
ymin: float,
ymax: float,
zmin: float,
zmax: float,
max_p: float,
n: int,
) -> None:
self.xmin = xmin self.xmin = xmin
self.xmax = xmax self.xmax = xmax
self.ymin = ymin self.ymin = ymin
@@ -28,13 +43,13 @@ class UnrestrictedModel(Model):
self.rng = numpy.random.default_rng() self.rng = numpy.random.default_rng()
def __repr__(self) -> str: def __repr__(self) -> str:
return f'UnrestrictedModel({self.xmin}, {self.xmax}, {self.ymin}, {self.ymax}, {self.zmin}, {self.zmax}, {self.max_p}, {self.n()})' return f"UnrestrictedModel({self.xmin}, {self.xmax}, {self.ymin}, {self.ymax}, {self.zmin}, {self.zmax}, {self.max_p}, {self.n()})"
def point_length(self) -> int: def point_length(self) -> int:
''' """
Dipole is unconstrained in this model. Dipole is unconstrained in this model.
All seven degrees of freedom: (px, py, pz, sx, sy, sz, w). All seven degrees of freedom: (px, py, pz, sx, sy, sz, w).
''' """
return 7 return 7
def n(self) -> int: def n(self) -> int:
@@ -46,27 +61,37 @@ class UnrestrictedModel(Model):
w = pt[6] w = pt[6]
diff = dot.r - s diff = dot.r - s
alpha = p.dot(diff) / (numpy.linalg.norm(diff)**3) alpha = p.dot(diff) / (numpy.linalg.norm(diff) ** 3)
b = (1 / numpy.pi) * (w / (w**2 + dot.f**2)) b = (1 / numpy.pi) * (w / (w**2 + dot.f**2))
return alpha**2 * b return alpha**2 * b
def jac_for_point_at_dot(self, dot: DotMeasurement, pt: numpy.ndarray) -> numpy.ndarray: def jac_for_point_at_dot(
self, dot: DotMeasurement, pt: numpy.ndarray
) -> numpy.ndarray:
p = pt[0:3] p = pt[0:3]
s = pt[3:6] s = pt[3:6]
w = pt[6] w = pt[6]
diff = dot.r - s diff = dot.r - s
alpha = p.dot(diff) / (numpy.linalg.norm(diff)**3) alpha = p.dot(diff) / (numpy.linalg.norm(diff) ** 3)
b = (1 / numpy.pi) * (w / (w**2 + dot.f**2)) b = (1 / numpy.pi) * (w / (w**2 + dot.f**2))
p_divs = 2 * alpha * diff / (numpy.linalg.norm(diff)**3) * b p_divs = 2 * alpha * diff / (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 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 f2 = dot.f**2
w2 = w**2 w2 = w**2
w_div = alpha**2 * (1 / numpy.pi) * ((f2 - w2) / ((f2 + w2)**2)) w_div = alpha**2 * (1 / numpy.pi) * ((f2 - w2) / ((f2 + w2) ** 2))
return numpy.concatenate((p_divs, r_divs, w_div), axis=None) return numpy.concatenate((p_divs, r_divs, w_div), axis=None)
@@ -77,35 +102,44 @@ class UnrestrictedModel(Model):
px = p * numpy.sin(theta) * numpy.cos(phi) px = p * numpy.sin(theta) * numpy.cos(phi)
py = p * numpy.sin(theta) * numpy.sin(phi) py = p * numpy.sin(theta) * numpy.sin(phi)
pz = p * numpy.cos(theta) pz = p * numpy.cos(theta)
s_pts = numpy.array((self.rng.uniform(self.xmin, self.xmax), self.rng.uniform(self.ymin, self.ymax), self.rng.uniform(self.zmin, self.zmax))) s_pts = numpy.array(
return OscillatingDipoleArrangement([OscillatingDipole(numpy.array([px, py, pz]), s_pts, frequency)]) (
self.rng.uniform(self.xmin, self.xmax),
self.rng.uniform(self.ymin, self.ymax),
self.rng.uniform(self.zmin, self.zmax),
)
)
return OscillatingDipoleArrangement(
[OscillatingDipole(numpy.array([px, py, pz]), s_pts, frequency)]
)
@dataclass @dataclass
class UnrestrictedDiscretisation(Discretisation): class UnrestrictedDiscretisation(Discretisation):
''' """
Representation of a discretisation of a UnrestrictedModel. Representation of a discretisation of a UnrestrictedModel.
Also captures a rough maximum value of dipole. Also captures a rough maximum value of dipole.
Parameters Parameters
---------- ----------
model : UnrestrictedModel model : UnrestrictedModel
The parent model of the discretisation. The parent model of the discretisation.
num_px: int num_px: int
The number of partitions of the px. The number of partitions of the px.
num_py: int num_py: int
The number of partitions of the py. The number of partitions of the py.
num_pz: int num_pz: int
The number of partitions of pz. The number of partitions of pz.
num_x : int num_x : int
The number of partitions of the x axis. The number of partitions of the x axis.
num_y : int num_y : int
The number of partitions of the y axis. The number of partitions of the y axis.
num_z : int num_z : int
The number of partitions of the z axis. The number of partitions of the z axis.
max_p : int max_p : int
The maximum p coordinate in any direction. The maximum p coordinate in any direction.
''' """
model: UnrestrictedModel model: UnrestrictedModel
num_px: int num_px: int
num_py: int num_py: int
@@ -131,22 +165,34 @@ class UnrestrictedDiscretisation(Discretisation):
# We want to keep w unbounded, restrict sx, sy, sz, px and py based on step. # We want to keep w unbounded, restrict sx, sy, sz, px and py based on step.
return ( return (
[ [
pxi * self.px_step - self.max_p, pyi * self.py_step - self.max_p, pzi * self.pz_step - self.max_p, pxi * self.px_step - self.max_p,
xi * self.x_step + self.model.xmin, yi * self.y_step + self.model.ymin, zi * self.z_step + self.model.zmin, pyi * self.py_step - self.max_p,
-numpy.inf pzi * self.pz_step - self.max_p,
xi * self.x_step + self.model.xmin,
yi * self.y_step + self.model.ymin,
zi * self.z_step + self.model.zmin,
-numpy.inf,
], ],
[ [
(pxi + 1) * self.px_step - self.max_p, (pyi + 1) * self.py_step - self.max_p, (pzi + 1) * self.pz_step - self.max_p, (pxi + 1) * self.px_step - self.max_p,
(xi + 1) * self.x_step + self.model.xmin, (yi + 1) * self.y_step + self.model.ymin, (zi + 1) * self.z_step + self.model.zmin, (pyi + 1) * self.py_step - self.max_p,
numpy.inf (pzi + 1) * self.pz_step - self.max_p,
] (xi + 1) * self.x_step + self.model.xmin,
(yi + 1) * self.y_step + self.model.ymin,
(zi + 1) * self.z_step + self.model.zmin,
numpy.inf,
],
) )
def all_indices(self) -> numpy.ndindex: def all_indices(self) -> numpy.ndindex:
# see https://github.com/numpy/numpy/issues/20706 for why this is a mypy problem. # see https://github.com/numpy/numpy/issues/20706 for why this is a mypy problem.
return numpy.ndindex((self.num_px, self.num_py, self.num_pz, self.num_x, self.num_y, self.num_z)) # type:ignore return numpy.ndindex(
(self.num_px, self.num_py, self.num_pz, self.num_x, self.num_y, self.num_z)
) # type:ignore
def solve_for_index(self, dots: Sequence[DotMeasurement], index: Tuple[float, ...]) -> scipy.optimize.OptimizeResult: def solve_for_index(
self, dots: Sequence[DotMeasurement], index: Tuple[float, ...]
) -> scipy.optimize.OptimizeResult:
bounds = self.bounds(index) bounds = self.bounds(index)
px_mean = (bounds[0][0] + bounds[1][0]) / 2 px_mean = (bounds[0][0] + bounds[1][0]) / 2
py_mean = (bounds[0][1] + bounds[1][1]) / 2 py_mean = (bounds[0][1] + bounds[1][1]) / 2
@@ -154,4 +200,10 @@ class UnrestrictedDiscretisation(Discretisation):
sx_mean = (bounds[0][3] + bounds[1][3]) / 2 sx_mean = (bounds[0][3] + bounds[1][3]) / 2
sy_mean = (bounds[0][4] + bounds[1][4]) / 2 sy_mean = (bounds[0][4] + bounds[1][4]) / 2
sz_mean = (bounds[0][5] + bounds[1][5]) / 2 sz_mean = (bounds[0][5] + bounds[1][5]) / 2
return self.model.solve(dots, initial_pt=numpy.array([px_mean, py_mean, pz_mean, sx_mean, sy_mean, sz_mean, .1]), bounds=bounds) return self.model.solve(
dots,
initial_pt=numpy.array(
[px_mean, py_mean, pz_mean, sx_mean, sy_mean, sz_mean, 0.1]
),
bounds=bounds,
)

View File

@@ -5,10 +5,12 @@ import logging
_logger = logging.getLogger(__name__) _logger = logging.getLogger(__name__)
def fast_s_nonlocal(dot_pair_inputs: numpy.ndarray, dipoles: numpy.ndarray) -> numpy.ndarray: def fast_s_nonlocal(
''' dot_pair_inputs: numpy.ndarray, dipoles: numpy.ndarray
No error correction here baby. ) -> numpy.ndarray:
''' """
No error correction here baby.
"""
ps = dipoles[:, 0:3] ps = dipoles[:, 0:3]
ss = dipoles[:, 3:6] ss = dipoles[:, 3:6]
ws = dipoles[:, 6] ws = dipoles[:, 6]
@@ -31,19 +33,19 @@ def fast_s_nonlocal(dot_pair_inputs: numpy.ndarray, dipoles: numpy.ndarray) -> n
_logger.debug(f"diffses1: {diffses1}") _logger.debug(f"diffses1: {diffses1}")
_logger.debug(f"diffses2: {diffses2}") _logger.debug(f"diffses2: {diffses2}")
norms1 = numpy.linalg.norm(diffses1, axis=2)**3 norms1 = numpy.linalg.norm(diffses1, axis=2) ** 3
norms2 = numpy.linalg.norm(diffses2, axis=2)**3 norms2 = numpy.linalg.norm(diffses2, axis=2) ** 3
if _logger.isEnabledFor(logging.DEBUG): if _logger.isEnabledFor(logging.DEBUG):
_logger.debug(f"norms1: {norms1}") _logger.debug(f"norms1: {norms1}")
_logger.debug(f"norms2: {norms2}") _logger.debug(f"norms2: {norms2}")
alphses1 = numpy.einsum('...ji, ...i', diffses1, ps) / norms1 alphses1 = numpy.einsum("...ji, ...i", diffses1, ps) / norms1
alphses2 = numpy.einsum('...ji, ...i', diffses2, ps) / norms2 alphses2 = numpy.einsum("...ji, ...i", diffses2, ps) / norms2
if _logger.isEnabledFor(logging.DEBUG): if _logger.isEnabledFor(logging.DEBUG):
_logger.debug(f"alphses1: {alphses1}") _logger.debug(f"alphses1: {alphses1}")
_logger.debug(f"alphses2: {alphses2}") _logger.debug(f"alphses2: {alphses2}")
bses = (1 / numpy.pi) * (ws[:, None] / (f1s**2 + ws[:, None]**2)) bses = (1 / numpy.pi) * (ws[:, None] / (f1s**2 + ws[:, None] ** 2))
if _logger.isEnabledFor(logging.DEBUG): if _logger.isEnabledFor(logging.DEBUG):
_logger.debug(f"bses: {bses}") _logger.debug(f"bses: {bses}")

View File

@@ -5,10 +5,12 @@ import logging
_logger = logging.getLogger(__name__) _logger = logging.getLogger(__name__)
def fast_vs_for_dipoles(dot_inputs: numpy.ndarray, dipoles: numpy.ndarray) -> numpy.ndarray: def fast_vs_for_dipoles(
''' dot_inputs: numpy.ndarray, dipoles: numpy.ndarray
No error correction here baby. ) -> numpy.ndarray:
''' """
No error correction here baby.
"""
ps = dipoles[:, 0:3] ps = dipoles[:, 0:3]
ss = dipoles[:, 3:6] ss = dipoles[:, 3:6]
ws = dipoles[:, 6] ws = dipoles[:, 6]
@@ -23,18 +25,18 @@ def fast_vs_for_dipoles(dot_inputs: numpy.ndarray, dipoles: numpy.ndarray) -> nu
diffses = rs - ss[:, None] diffses = rs - ss[:, None]
_logger.debug(f"diffses: {diffses}") _logger.debug(f"diffses: {diffses}")
norms = numpy.linalg.norm(diffses, axis=2)**3 norms = numpy.linalg.norm(diffses, axis=2) ** 3
_logger.debug(f"norms: {norms}") _logger.debug(f"norms: {norms}")
ases = (numpy.einsum('...ji, ...i', diffses, ps) / norms)**2 ases = (numpy.einsum("...ji, ...i", diffses, ps) / norms) ** 2
_logger.debug(f"ases: {ases}") _logger.debug(f"ases: {ases}")
bses = (1 / numpy.pi) * (ws[:, None] / (fs**2 + ws[:, None]**2)) bses = (1 / numpy.pi) * (ws[:, None] / (fs**2 + ws[:, None] ** 2))
_logger.debug(f"bses: {bses}") _logger.debug(f"bses: {bses}")
return ases * bses return ases * bses
def between(a: numpy.ndarray, low: numpy.ndarray, high: numpy.ndarray) -> numpy.ndarray: def between(a: numpy.ndarray, low: numpy.ndarray, high: numpy.ndarray) -> numpy.ndarray:
''' """
Intended specifically for the case where a is a list of arrays, and each array must be between the single array low and high, but without error checking. Intended specifically for the case where a is a list of arrays, and each array must be between the single array low and high, but without error checking.
''' """
return numpy.all(numpy.logical_and(low < a, high > a), axis=1) return numpy.all(numpy.logical_and(low < a, high > a), axis=1)

View File

@@ -26,7 +26,20 @@ def flip_chunk_to_positive_px(pt: numpy.ndarray) -> numpy.ndarray:
def normalise_point_list(pts: numpy.ndarray, pt_length) -> numpy.ndarray: 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)] 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)) range_to_length = list(range(pt_length))
rotated_range = range_to_length[pt_length - 1:] + range_to_length[0:pt_length - 1] rotated_range = (
return numpy.concatenate(sorted(chunked_pts, key=lambda x: tuple(round(val, 3) for val in operator.itemgetter(*rotated_range)(x))), axis=None) 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,
)

View File

@@ -8,7 +8,12 @@ def test_inputs_with_frequency_range():
frequencies = [5, 7, 9] frequencies = [5, 7, 9]
expected = [ expected = [
([1, 2, 3], 5), ([1, 2, 3], 7), ([1, 2, 3], 9), ([2, 4, 6], 7), ([2, 4, 6], 9), ([2, 4, 6], 5) ([1, 2, 3], 5),
([1, 2, 3], 7),
([1, 2, 3], 9),
([2, 4, 6], 7),
([2, 4, 6], 9),
([2, 4, 6], 5),
] ]
actual = pdme.inputs.inputs_with_frequency_range([dot1, dot2], frequencies) actual = pdme.inputs.inputs_with_frequency_range([dot1, dot2], frequencies)
@@ -23,7 +28,9 @@ def test_input_pairs_with_frequency_range():
frequencies = [5, 7, 9] frequencies = [5, 7, 9]
expected = [ expected = [
([1, 2, 3], [2, 4, 6], 5), ([1, 2, 3], [2, 4, 6], 7), ([1, 2, 3], [2, 4, 6], 9) ([1, 2, 3], [2, 4, 6], 5),
([1, 2, 3], [2, 4, 6], 7),
([1, 2, 3], [2, 4, 6], 9),
] ]
actual = pdme.inputs.input_pairs_with_frequency_range([dot1, dot2], frequencies) actual = pdme.inputs.input_pairs_with_frequency_range([dot1, dot2], frequencies)

View File

@@ -12,14 +12,30 @@ def test_static_dipole():
expected_v1 = 0.00001421963287022476 expected_v1 = 0.00001421963287022476
expected_v2 = 0.00001107180225755457 expected_v2 = 0.00001107180225755457
numpy.testing.assert_allclose(d1.s_at_position(dot_position1, dot_frequency1), expected_v1, err_msg="Voltage at dot isn't as expected.") numpy.testing.assert_allclose(
d1.s_at_position(dot_position1, dot_frequency1),
expected_v1,
err_msg="Voltage at dot isn't as expected.",
)
dot_measurements = dipoles.get_dot_measurements([(dot_position1, dot_frequency1)]) dot_measurements = dipoles.get_dot_measurements([(dot_position1, dot_frequency1)])
assert len(dot_measurements) == 1, "Should have only had one dot measurement." assert len(dot_measurements) == 1, "Should have only had one dot measurement."
measurement = dot_measurements[0] measurement = dot_measurements[0]
numpy.testing.assert_allclose(measurement.r, dot_position1, err_msg="Dot position should have been passed over") numpy.testing.assert_allclose(
numpy.testing.assert_allclose(measurement.f, dot_frequency1, err_msg="Dot frequency should have been passed over") measurement.r,
numpy.testing.assert_allclose(measurement.v, expected_v1 + expected_v2, err_msg="Voltage at dot isn't as expected.") dot_position1,
err_msg="Dot position should have been passed over",
)
numpy.testing.assert_allclose(
measurement.f,
dot_frequency1,
err_msg="Dot frequency should have been passed over",
)
numpy.testing.assert_allclose(
measurement.v,
expected_v1 + expected_v2,
err_msg="Voltage at dot isn't as expected.",
)
def test_dipole_dot_pair(): def test_dipole_dot_pair():
@@ -31,8 +47,21 @@ def test_dipole_dot_pair():
dot_frequency = 8 dot_frequency = 8
expected_sij = 0.000083328037100902801698 expected_sij = 0.000083328037100902801698
numpy.testing.assert_allclose(d1.s_for_dot_pair(dot_position1, dot_position2, dot_frequency), expected_sij, err_msg="Sij for dot pair isn't as expected.") numpy.testing.assert_allclose(
numpy.testing.assert_allclose([m.v for m in dipoles.get_dot_pair_measurements([(dot_position1, dot_position2, dot_frequency)])], [expected_sij], err_msg="Sij for dot pair isn't as expected via dipole.") d1.s_for_dot_pair(dot_position1, dot_position2, dot_frequency),
expected_sij,
err_msg="Sij for dot pair isn't as expected.",
)
numpy.testing.assert_allclose(
[
m.v
for m in dipoles.get_dot_pair_measurements(
[(dot_position1, dot_position2, dot_frequency)]
)
],
[expected_sij],
err_msg="Sij for dot pair isn't as expected via dipole.",
)
def test_range_pairs(): def test_range_pairs():
@@ -44,10 +73,16 @@ def test_range_pairs():
dot_frequency = 8 dot_frequency = 8
expected_sij = 0.000083328037100902801698 expected_sij = 0.000083328037100902801698
actuals = dipoles.get_percent_range_dot_pair_measurements([(dot_position1, dot_position2, dot_frequency)], 0.5, 1.5) actuals = dipoles.get_percent_range_dot_pair_measurements(
[(dot_position1, dot_position2, dot_frequency)], 0.5, 1.5
)
assert len(actuals) == 1, "should have only been one pair" assert len(actuals) == 1, "should have only been one pair"
actual = actuals[0] actual = actuals[0]
numpy.testing.assert_allclose([actual.v_low, actual.v_high], expected_sij * numpy.array([0.5, 1.5]), err_msg="Sij for dot pair isn't as expected via dipole with range.") numpy.testing.assert_allclose(
[actual.v_low, actual.v_high],
expected_sij * numpy.array([0.5, 1.5]),
err_msg="Sij for dot pair isn't as expected via dipole with range.",
)
def test_range_dipole_measurements(): def test_range_dipole_measurements():
@@ -60,12 +95,34 @@ def test_range_dipole_measurements():
expected_v1 = 0.00001421963287022476 expected_v1 = 0.00001421963287022476
expected_v2 = 0.00001107180225755457 expected_v2 = 0.00001107180225755457
numpy.testing.assert_allclose(d1.s_at_position(dot_position1, dot_frequency1), expected_v1, err_msg="Voltage at dot isn't as expected.") numpy.testing.assert_allclose(
d1.s_at_position(dot_position1, dot_frequency1),
expected_v1,
err_msg="Voltage at dot isn't as expected.",
)
range_dot_measurements = dipoles.get_percent_range_dot_measurements([(dot_position1, dot_frequency1)], 0.5, 1.5) range_dot_measurements = dipoles.get_percent_range_dot_measurements(
[(dot_position1, dot_frequency1)], 0.5, 1.5
)
assert len(range_dot_measurements) == 1, "Should have only had one dot measurement." assert len(range_dot_measurements) == 1, "Should have only had one dot measurement."
range_measurement = range_dot_measurements[0] range_measurement = range_dot_measurements[0]
numpy.testing.assert_allclose(range_measurement.r, dot_position1, err_msg="Dot position should have been passed over") numpy.testing.assert_allclose(
numpy.testing.assert_allclose(range_measurement.f, dot_frequency1, err_msg="Dot frequency should have been passed over") range_measurement.r,
numpy.testing.assert_allclose(range_measurement.v_low, (expected_v1 + expected_v2) / 2, err_msg="Lower voltage at dot isn't as expected.") dot_position1,
numpy.testing.assert_allclose(range_measurement.v_high, (expected_v1 + expected_v2) * 3 / 2, err_msg="Lower oltage at dot isn't as expected.") err_msg="Dot position should have been passed over",
)
numpy.testing.assert_allclose(
range_measurement.f,
dot_frequency1,
err_msg="Dot frequency should have been passed over",
)
numpy.testing.assert_allclose(
range_measurement.v_low,
(expected_v1 + expected_v2) / 2,
err_msg="Lower voltage at dot isn't as expected.",
)
numpy.testing.assert_allclose(
range_measurement.v_high,
(expected_v1 + expected_v2) * 3 / 2,
err_msg="Lower oltage at dot isn't as expected.",
)

View File

@@ -18,5 +18,7 @@ def test_dipole_to_array():
arrangement = pdme.measurement.OscillatingDipoleArrangement([d1, d2]) arrangement = pdme.measurement.OscillatingDipoleArrangement([d1, d2])
numpy.testing.assert_array_equal( numpy.testing.assert_array_equal(
[expected1, expected2], arrangement.to_numpy_array(), err_msg="Didn't convert multiple dipoles right" [expected1, expected2],
arrangement.to_numpy_array(),
err_msg="Didn't convert multiple dipoles right",
) )

View File

@@ -9,8 +9,16 @@ def test_swap_high_low():
m1 = DotRangeMeasurement(actual_high, actual_low, 100, 1000) m1 = DotRangeMeasurement(actual_high, actual_low, 100, 1000)
m2 = DotRangeMeasurement(actual_low, actual_high, 100, 1000) m2 = DotRangeMeasurement(actual_low, actual_high, 100, 1000)
numpy.testing.assert_array_equal([m1.v_low, m1.v_high], [actual_low, actual_high], err_msg="Highs were wrong with swap") numpy.testing.assert_array_equal(
numpy.testing.assert_array_equal([m2.v_low, m2.v_high], [actual_low, actual_high], err_msg="Highs were wrong without swap") [m1.v_low, m1.v_high],
[actual_low, actual_high],
err_msg="Highs were wrong with swap",
)
numpy.testing.assert_array_equal(
[m2.v_low, m2.v_high],
[actual_low, actual_high],
err_msg="Highs were wrong without swap",
)
def test_swap_high_low_negative(): def test_swap_high_low_negative():
@@ -19,8 +27,16 @@ def test_swap_high_low_negative():
m1 = DotRangeMeasurement(actual_high, actual_low, 100, 1000) m1 = DotRangeMeasurement(actual_high, actual_low, 100, 1000)
m2 = DotRangeMeasurement(actual_low, actual_high, 100, 1000) m2 = DotRangeMeasurement(actual_low, actual_high, 100, 1000)
numpy.testing.assert_array_equal([m1.v_low, m1.v_high], [actual_low, actual_high], err_msg="Highs were wrong with swap, negative") numpy.testing.assert_array_equal(
numpy.testing.assert_array_equal([m2.v_low, m2.v_high], [actual_low, actual_high], err_msg="Highs were wrong without swap, negative") [m1.v_low, m1.v_high],
[actual_low, actual_high],
err_msg="Highs were wrong with swap, negative",
)
numpy.testing.assert_array_equal(
[m2.v_low, m2.v_high],
[actual_low, actual_high],
err_msg="Highs were wrong without swap, negative",
)
def test_swap_high_low_pair(): def test_swap_high_low_pair():
@@ -29,8 +45,16 @@ def test_swap_high_low_pair():
m1 = DotPairRangeMeasurement(actual_high, actual_low, 100, 1000, 10000) m1 = DotPairRangeMeasurement(actual_high, actual_low, 100, 1000, 10000)
m2 = DotPairRangeMeasurement(actual_low, actual_high, 100, 1000, 10000) m2 = DotPairRangeMeasurement(actual_low, actual_high, 100, 1000, 10000)
numpy.testing.assert_array_equal([m1.v_low, m1.v_high], [actual_low, actual_high], err_msg="Highs were wrong with swap") numpy.testing.assert_array_equal(
numpy.testing.assert_array_equal([m2.v_low, m2.v_high], [actual_low, actual_high], err_msg="Highs were wrong without swap") [m1.v_low, m1.v_high],
[actual_low, actual_high],
err_msg="Highs were wrong with swap",
)
numpy.testing.assert_array_equal(
[m2.v_low, m2.v_high],
[actual_low, actual_high],
err_msg="Highs were wrong without swap",
)
def test_swap_high_low_pair_negative(): def test_swap_high_low_pair_negative():
@@ -39,5 +63,13 @@ def test_swap_high_low_pair_negative():
m1 = DotPairRangeMeasurement(actual_high, actual_low, 100, 1000, 10000) m1 = DotPairRangeMeasurement(actual_high, actual_low, 100, 1000, 10000)
m2 = DotPairRangeMeasurement(actual_low, actual_high, 100, 1000, 10000) m2 = DotPairRangeMeasurement(actual_low, actual_high, 100, 1000, 10000)
numpy.testing.assert_array_equal([m1.v_low, m1.v_high], [actual_low, actual_high], err_msg="Highs were wrong with swap, negative") numpy.testing.assert_array_equal(
numpy.testing.assert_array_equal([m2.v_low, m2.v_high], [actual_low, actual_high], err_msg="Highs were wrong without swap, negative") [m1.v_low, m1.v_high],
[actual_low, actual_high],
err_msg="Highs were wrong with swap, negative",
)
numpy.testing.assert_array_equal(
[m2.v_low, m2.v_high],
[actual_low, actual_high],
err_msg="Highs were wrong without swap, negative",
)

View File

@@ -11,7 +11,9 @@ def test_inputs_to_array():
actual = pdme.measurement.input_types.dot_inputs_to_array([i1, i2]) actual = pdme.measurement.input_types.dot_inputs_to_array([i1, i2])
expected = numpy.array([[1, 2, 3, 5], [-1, 4, -2, 9]]) expected = numpy.array([[1, 2, 3, 5], [-1, 4, -2, 9]])
numpy.testing.assert_allclose(actual, expected, err_msg="Didn't convert to array properly") numpy.testing.assert_allclose(
actual, expected, err_msg="Didn't convert to array properly"
)
def test_pair_inputs_to_array(): def test_pair_inputs_to_array():
@@ -19,19 +21,26 @@ def test_pair_inputs_to_array():
i2 = ([-1, 4, -2], [6, 7, 8], 9) i2 = ([-1, 4, -2], [6, 7, 8], 9)
actual = pdme.measurement.input_types.dot_pair_inputs_to_array([i1, i2]) actual = pdme.measurement.input_types.dot_pair_inputs_to_array([i1, i2])
expected = numpy.array([ expected = numpy.array(
[[1, 2, 3, 5], [-1, 4, -2, 5]], [
[[-1, 4, -2, 9], [6, 7, 8, 9]], [[1, 2, 3, 5], [-1, 4, -2, 5]],
]) [[-1, 4, -2, 9], [6, 7, 8, 9]],
]
)
numpy.testing.assert_allclose(actual, expected, err_msg="Didn't convert to array properly") numpy.testing.assert_allclose(
actual, expected, err_msg="Didn't convert to array properly"
)
def test_ranges_to_array(): def test_ranges_to_array():
m1 = DotRangeMeasurement(1, 2, 100, 1000) m1 = DotRangeMeasurement(1, 2, 100, 1000)
m2 = DotRangeMeasurement(0.5, 3, 100, 1000) m2 = DotRangeMeasurement(0.5, 3, 100, 1000)
actual_lows, actual_highs = pdme.measurement.input_types.dot_range_measurements_low_high_arrays([m1, m2]) (
actual_lows,
actual_highs,
) = pdme.measurement.input_types.dot_range_measurements_low_high_arrays([m1, m2])
numpy.testing.assert_allclose(actual_lows, [1, 0.5], err_msg="Lows were wrong") numpy.testing.assert_allclose(actual_lows, [1, 0.5], err_msg="Lows were wrong")
numpy.testing.assert_allclose(actual_highs, [2, 3], err_msg="Highs were wrong") numpy.testing.assert_allclose(actual_highs, [2, 3], err_msg="Highs were wrong")
@@ -40,6 +49,9 @@ def test_pair_ranges_to_array():
m1 = DotPairRangeMeasurement(1, 2, 100, 1000, 10000) m1 = DotPairRangeMeasurement(1, 2, 100, 1000, 10000)
m2 = DotPairRangeMeasurement(0.5, 3, 100, 1000, 10000) m2 = DotPairRangeMeasurement(0.5, 3, 100, 1000, 10000)
actual_lows, actual_highs = pdme.measurement.input_types.dot_range_measurements_low_high_arrays([m1, m2]) (
actual_lows,
actual_highs,
) = pdme.measurement.input_types.dot_range_measurements_low_high_arrays([m1, m2])
numpy.testing.assert_allclose(actual_lows, [1, 0.5], err_msg="Lows were wrong") numpy.testing.assert_allclose(actual_lows, [1, 0.5], err_msg="Lows were wrong")
numpy.testing.assert_allclose(actual_highs, [2, 3], err_msg="Highs were wrong") numpy.testing.assert_allclose(actual_highs, [2, 3], err_msg="Highs were wrong")

View File

@@ -10,10 +10,26 @@ import pytest
def test_fixed_dipole_model_solve_basic(): def test_fixed_dipole_model_solve_basic():
# Initialise our dipole arrangement and create dot measurements along a square. # Initialise our dipole arrangement and create dot measurements along a square.
fixed_dipole_moment = numpy.array((2, 1, 2)) fixed_dipole_moment = numpy.array((2, 1, 2))
dipoles = OscillatingDipoleArrangement([OscillatingDipole(fixed_dipole_moment, (1, 2, 4), 6)]) dipoles = OscillatingDipoleArrangement(
dot_inputs = list(itertools.chain.from_iterable( [OscillatingDipole(fixed_dipole_moment, (1, 2, 4), 6)]
(([1, 2, 1], f), ([1, 1, -2], f), ([1.5, 2, 0.1], f), ([1.5, 1, -0.2], f), ([2, 1, 0], f), ([2, 2, 0], f), ([0, 2, -.1], f), ([0, 1, 0.04], f), ([2, 0, 2], f), ([1, 0, 0], f)) for f in numpy.arange(1, 10, 1) )
)) dot_inputs = list(
itertools.chain.from_iterable(
(
([1, 2, 1], f),
([1, 1, -2], f),
([1.5, 2, 0.1], f),
([1.5, 1, -0.2], f),
([2, 1, 0], f),
([2, 2, 0], f),
([0, 2, -0.1], f),
([0, 1, 0.04], f),
([2, 0, 2], f),
([1, 0, 0], f),
)
for f in numpy.arange(1, 10, 1)
)
)
dots = dipoles.get_dot_measurements(dot_inputs) dots = dipoles.get_dot_measurements(dot_inputs)
model = FixedDipoleModel(-3, 3, -3, 3, 3, 5, fixed_dipole_moment, 1) model = FixedDipoleModel(-3, 3, -3, 3, 3, 5, fixed_dipole_moment, 1)
@@ -24,4 +40,10 @@ def test_fixed_dipole_model_solve_basic():
result = model.solve(dots) result = model.solve(dots)
logging.info(result) logging.info(result)
assert result.success assert result.success
numpy.testing.assert_allclose(result.normalised_x, expected_solution, err_msg="Even well specified problem solution was wrong.", rtol=1e-6, atol=1e-11) numpy.testing.assert_allclose(
result.normalised_x,
expected_solution,
err_msg="Even well specified problem solution was wrong.",
rtol=1e-6,
atol=1e-11,
)

View File

@@ -8,9 +8,23 @@ import itertools
def test_fixed_magnitude_model_solve_basic(): def test_fixed_magnitude_model_solve_basic():
# Initialise our dipole arrangement and create dot measurements along a square. # Initialise our dipole arrangement and create dot measurements along a square.
dipoles = OscillatingDipoleArrangement([OscillatingDipole((2, 0, 0), (1, 2, 4), 1)]) dipoles = OscillatingDipoleArrangement([OscillatingDipole((2, 0, 0), (1, 2, 4), 1)])
dot_inputs = list(itertools.chain.from_iterable( dot_inputs = list(
(([1, 2, 0.01], f), ([1, 1, -0.2], f), ([1.5, 2, 0.01], f), ([1.5, 1, -0.2], f), ([2, 1, 0], f), ([2, 2, 0], f), ([0, 2, -.1], f), ([0, 1, 0.04], f), ([2, 0, 0], f), ([1, 0, 0], f)) for f in numpy.arange(1, 10, 2) itertools.chain.from_iterable(
)) (
([1, 2, 0.01], f),
([1, 1, -0.2], f),
([1.5, 2, 0.01], f),
([1.5, 1, -0.2], f),
([2, 1, 0], f),
([2, 2, 0], f),
([0, 2, -0.1], f),
([0, 1, 0.04], f),
([2, 0, 0], f),
([1, 0, 0], f),
)
for f in numpy.arange(1, 10, 2)
)
)
dots = dipoles.get_dot_measurements(dot_inputs) dots = dipoles.get_dot_measurements(dot_inputs)
model = FixedMagnitudeModel(-5, 5, -5, 5, -5, 5, 2, 1) model = FixedMagnitudeModel(-5, 5, -5, 5, -5, 5, 2, 1)
@@ -21,14 +35,38 @@ def test_fixed_magnitude_model_solve_basic():
result = model.solve(dots) result = model.solve(dots)
logging.info(result) logging.info(result)
assert result.success assert result.success
numpy.testing.assert_allclose(result.normalised_x, expected_solution, err_msg="Even well specified problem solution was wrong.", rtol=1e-6, atol=1e-11) numpy.testing.assert_allclose(
result.normalised_x,
expected_solution,
err_msg="Even well specified problem solution was wrong.",
rtol=1e-6,
atol=1e-11,
)
solved_dipoles = model.solution_as_dipoles(result.normalised_x) solved_dipoles = model.solution_as_dipoles(result.normalised_x)
assert len(solved_dipoles) == 1 assert len(solved_dipoles) == 1
numpy.testing.assert_allclose(solved_dipoles[0].p, (2, 0, 0), err_msg="Shove it in a dipole correctly.", rtol=1e-6, atol=1e-11) numpy.testing.assert_allclose(
numpy.testing.assert_allclose(solved_dipoles[0].s, (1, 2, 4), err_msg="Shove it in a dipole correctly.", rtol=1e-6, atol=1e-11) solved_dipoles[0].p,
numpy.testing.assert_allclose(solved_dipoles[0].w, 1, err_msg="Shove it in a dipole correctly.", rtol=1e-6, atol=1e-11) (2, 0, 0),
err_msg="Shove it in a dipole correctly.",
rtol=1e-6,
atol=1e-11,
)
numpy.testing.assert_allclose(
solved_dipoles[0].s,
(1, 2, 4),
err_msg="Shove it in a dipole correctly.",
rtol=1e-6,
atol=1e-11,
)
numpy.testing.assert_allclose(
solved_dipoles[0].w,
1,
err_msg="Shove it in a dipole correctly.",
rtol=1e-6,
atol=1e-11,
)
def test_fixed_magnitude_model_repr(): def test_fixed_magnitude_model_repr():

View File

@@ -17,9 +17,12 @@ def test_fixed_z_plane_model_solve_error_initial():
def test_fixed_z_plane_model_solve_basic(): def test_fixed_z_plane_model_solve_basic():
# Initialise our dipole arrangement and create dot measurements along a square. # Initialise our dipole arrangement and create dot measurements along a square.
dipoles = OscillatingDipoleArrangement([OscillatingDipole((0, 0, 2), (1, 2, 4), 1)]) dipoles = OscillatingDipoleArrangement([OscillatingDipole((0, 0, 2), (1, 2, 4), 1)])
dot_inputs = list(itertools.chain.from_iterable( dot_inputs = list(
(([1, 2, 0], f), ([1, 1, 0], f), ([2, 1, 0], f), ([2, 2, 0], f)) for f in numpy.arange(1, 10, 2) itertools.chain.from_iterable(
)) (([1, 2, 0], f), ([1, 1, 0], f), ([2, 1, 0], f), ([2, 2, 0], f))
for f in numpy.arange(1, 10, 2)
)
)
dots = dipoles.get_dot_measurements(dot_inputs) dots = dipoles.get_dot_measurements(dot_inputs)
model = FixedZPlaneModel(4, -10, 10, -10, 10, 1) model = FixedZPlaneModel(4, -10, 10, -10, 10, 1)
@@ -30,10 +33,22 @@ def test_fixed_z_plane_model_solve_basic():
result = model.solve(dots) result = model.solve(dots)
logging.info(result) logging.info(result)
assert result.success assert result.success
numpy.testing.assert_allclose(result.x, expected_solution, err_msg="Even well specified problem solution was wrong.", rtol=1e-6, atol=1e-11) numpy.testing.assert_allclose(
result.x,
expected_solution,
err_msg="Even well specified problem solution was wrong.",
rtol=1e-6,
atol=1e-11,
)
# Do it again with an initial point # Do it again with an initial point
result = model.solve(dots, initial_pt=[2, 2, 2, 2]) result = model.solve(dots, initial_pt=[2, 2, 2, 2])
logging.info(result) logging.info(result)
assert result.success assert result.success
numpy.testing.assert_allclose(result.x, expected_solution, err_msg="Even well specified problem solution was wrong.", rtol=1e-6, atol=1e-11) numpy.testing.assert_allclose(
result.x,
expected_solution,
err_msg="Even well specified problem solution was wrong.",
rtol=1e-6,
atol=1e-11,
)

View File

@@ -11,5 +11,10 @@ def test_fixed_z_plane_model_discretization():
assert discretisation.pz_step == 30 assert discretisation.pz_step == 30
assert discretisation.x_step == 10 assert discretisation.x_step == 10
assert discretisation.y_step == 4 assert discretisation.y_step == 4
numpy.testing.assert_allclose(discretisation.bounds((0, 0, 0)), ((-15, -10, -10, -numpy.inf), (15, 0, -6, numpy.inf))) numpy.testing.assert_allclose(
numpy.testing.assert_allclose(list(discretisation.all_indices()), list(numpy.ndindex((1, 2, 5)))) discretisation.bounds((0, 0, 0)),
((-15, -10, -10, -numpy.inf), (15, 0, -6, numpy.inf)),
)
numpy.testing.assert_allclose(
list(discretisation.all_indices()), list(numpy.ndindex((1, 2, 5)))
)

View File

@@ -20,13 +20,27 @@ def test_fixed_z_plane_model_cost_and_jac_single():
expected_cost = [0.0000946746] expected_cost = [0.0000946746]
actual_cost = cost_function(pt) actual_cost = cost_function(pt)
numpy.testing.assert_allclose(actual_cost, expected_cost, err_msg="Cost wasn't as expected.", rtol=1e-6, atol=1e-11) numpy.testing.assert_allclose(
actual_cost,
expected_cost,
err_msg="Cost wasn't as expected.",
rtol=1e-6,
atol=1e-11,
)
jac_function = model.jac([dot]) jac_function = model.jac([dot])
expected_jac = [[0.0002859666151836802, -0.0001009293935942401, 0, 0.0001035396365320221]] expected_jac = [
[0.0002859666151836802, -0.0001009293935942401, 0, 0.0001035396365320221]
]
actual_jac = jac_function(pt) actual_jac = jac_function(pt)
logging.warning(actual_jac) logging.warning(actual_jac)
numpy.testing.assert_allclose(actual_jac, expected_jac, err_msg="Jac wasn't as expected.", rtol=1e-6, atol=1e-11) numpy.testing.assert_allclose(
actual_jac,
expected_jac,
err_msg="Jac wasn't as expected.",
rtol=1e-6,
atol=1e-11,
)

View File

@@ -7,10 +7,26 @@ import itertools
def test_unrestricted_model_solve_basic(): def test_unrestricted_model_solve_basic():
# Initialise our dipole arrangement and create dot measurements along a square. # Initialise our dipole arrangement and create dot measurements along a square.
dipoles = OscillatingDipoleArrangement([OscillatingDipole((.2, 0, 2), (1, 2, 4), 1)]) dipoles = OscillatingDipoleArrangement(
dot_inputs = list(itertools.chain.from_iterable( [OscillatingDipole((0.2, 0, 2), (1, 2, 4), 1)]
(([1, 2, 0.01], f), ([1, 1, -0.2], f), ([1.5, 2, 0.01], f), ([1.5, 1, -0.2], f), ([2, 1, 0], f), ([2, 2, 0], f), ([0, 2, -.1], f), ([0, 1, 0.04], f), ([2, 0, 0], f), ([1, 0, 0], f)) for f in numpy.arange(1, 10, 2) )
)) dot_inputs = list(
itertools.chain.from_iterable(
(
([1, 2, 0.01], f),
([1, 1, -0.2], f),
([1.5, 2, 0.01], f),
([1.5, 1, -0.2], f),
([2, 1, 0], f),
([2, 2, 0], f),
([0, 2, -0.1], f),
([0, 1, 0.04], f),
([2, 0, 0], f),
([1, 0, 0], f),
)
for f in numpy.arange(1, 10, 2)
)
)
dots = dipoles.get_dot_measurements(dot_inputs) dots = dipoles.get_dot_measurements(dot_inputs)
model = UnrestrictedModel(1, -1, 1, -1, 1, -1, 5, 1) model = UnrestrictedModel(1, -1, 1, -1, 1, -1, 5, 1)
@@ -21,4 +37,10 @@ def test_unrestricted_model_solve_basic():
result = model.solve(dots) result = model.solve(dots)
logging.info(result) logging.info(result)
assert result.success assert result.success
numpy.testing.assert_allclose(result.normalised_x, expected_solution, err_msg="Even well specified problem solution was wrong.", rtol=1e-6, atol=1e-11) numpy.testing.assert_allclose(
result.normalised_x,
expected_solution,
err_msg="Even well specified problem solution was wrong.",
rtol=1e-6,
atol=1e-11,
)

View File

@@ -14,5 +14,10 @@ def test_unrestricted_model_discretization():
assert discretisation.x_step == 10 assert discretisation.x_step == 10
assert discretisation.y_step == 4 assert discretisation.y_step == 4
assert discretisation.z_step == 20 assert discretisation.z_step == 20
numpy.testing.assert_allclose(discretisation.bounds((0, 0, 0, 0, 0, 0)), ((-15, -15, -15, -10, -10, -10, -numpy.inf), (15, 15, 0, 0, -6, 10, numpy.inf))) numpy.testing.assert_allclose(
numpy.testing.assert_allclose(list(discretisation.all_indices()), list(numpy.ndindex((1, 1, 2, 2, 5, 1)))) discretisation.bounds((0, 0, 0, 0, 0, 0)),
((-15, -15, -15, -10, -10, -10, -numpy.inf), (15, 15, 0, 0, -6, 10, numpy.inf)),
)
numpy.testing.assert_allclose(
list(discretisation.all_indices()), list(numpy.ndindex((1, 1, 2, 2, 5, 1)))
)

View File

@@ -20,19 +20,35 @@ def test_unrestricted_model_cost_and_jac_single():
expected_cost = [0.0000946746] expected_cost = [0.0000946746]
actual_cost = cost_function(pt) actual_cost = cost_function(pt)
numpy.testing.assert_allclose(actual_cost, expected_cost, err_msg="Cost wasn't as expected.", rtol=1e-6, atol=1e-11) numpy.testing.assert_allclose(
actual_cost,
expected_cost,
err_msg="Cost wasn't as expected.",
rtol=1e-6,
atol=1e-11,
)
jac_function = model.jac([dot]) jac_function = model.jac([dot])
expected_jac = [ expected_jac = [
[ [
0.00007149165379592005, 0, 0.0002859666151836802, 0.00007149165379592005,
-0.0001009293935942401, 0, -0.0002607342667851202, 0,
0.0001035396365320221 0.0002859666151836802,
-0.0001009293935942401,
0,
-0.0002607342667851202,
0.0001035396365320221,
] ]
] ]
actual_jac = jac_function(pt) actual_jac = jac_function(pt)
logging.warning(actual_jac) logging.warning(actual_jac)
numpy.testing.assert_allclose(actual_jac, expected_jac, err_msg="Jac wasn't as expected.", rtol=1e-6, atol=1e-11) numpy.testing.assert_allclose(
actual_jac,
expected_jac,
err_msg="Jac wasn't as expected.",
rtol=1e-6,
atol=1e-11,
)

View File

@@ -12,10 +12,9 @@ def test_fast_nonlocal_calc():
dipoles = numpy.array([d1, d2, d3, d4]) dipoles = numpy.array([d1, d2, d3, d4])
dot_pairs = numpy.array([ dot_pairs = numpy.array(
[[-1, -2, -3, 11], [-1, 2, 5, 11]], [[[-1, -2, -3, 11], [-1, 2, 5, 11]], [[-1, -2, -3, 6], [2, 4, 6, 6]]]
[[-1, -2, -3, 6], [2, 4, 6, 6]] )
])
# expected_ij is for pair i, dipole j # expected_ij is for pair i, dipole j
expected_11 = 0.000021124454334947546213 expected_11 = 0.000021124454334947546213
expected_12 = 0.000022184755131682365135 expected_12 = 0.000022184755131682365135
@@ -26,12 +25,23 @@ def test_fast_nonlocal_calc():
expected_23 = 0.000017558321044891869169 expected_23 = 0.000017558321044891869169
expected_24 = -0.000034714318479634499683 expected_24 = -0.000034714318479634499683
expected = numpy.array([[expected_11, expected_21], [expected_12, expected_22], [expected_13, expected_23], [expected_14, expected_24]]) expected = numpy.array(
[
[expected_11, expected_21],
[expected_12, expected_22],
[expected_13, expected_23],
[expected_14, expected_24],
]
)
# this is a bit silly but just set the logger to debug so that the coverage stats don't get affected by the debug statements. # this is a bit silly but just set the logger to debug so that the coverage stats don't get affected by the debug statements.
pdme.util.fast_nonlocal_spectrum._logger.setLevel(logging.DEBUG) pdme.util.fast_nonlocal_spectrum._logger.setLevel(logging.DEBUG)
numpy.testing.assert_allclose(pdme.util.fast_nonlocal_spectrum.fast_s_nonlocal(dot_pairs, dipoles), expected, err_msg="nonlocal voltages at dot aren't as expected.") numpy.testing.assert_allclose(
pdme.util.fast_nonlocal_spectrum.fast_s_nonlocal(dot_pairs, dipoles),
expected,
err_msg="nonlocal voltages at dot aren't as expected.",
)
def test_fast_nonlocal_frequency_check(): def test_fast_nonlocal_frequency_check():
@@ -39,9 +49,7 @@ def test_fast_nonlocal_frequency_check():
dipoles = numpy.array([d1]) dipoles = numpy.array([d1])
dot_pairs = numpy.array([ dot_pairs = numpy.array([[[-1, -2, -3, 11], [-1, 2, 5, 10]]])
[[-1, -2, -3, 11], [-1, 2, 5, 10]]
])
with pytest.raises(ValueError): with pytest.raises(ValueError):
pdme.util.fast_nonlocal_spectrum.fast_s_nonlocal(dot_pairs, dipoles) pdme.util.fast_nonlocal_spectrum.fast_s_nonlocal(dot_pairs, dipoles)

View File

@@ -17,14 +17,18 @@ def test_fast_v_calc():
expected = numpy.array([[expected_11, expected_21], [expected_12, expected_22]]) expected = numpy.array([[expected_11, expected_21], [expected_12, expected_22]])
numpy.testing.assert_allclose(pdme.util.fast_v_calc.fast_vs_for_dipoles(dot_inputs, dipoles), expected, err_msg="Voltages at dot aren't as expected.") numpy.testing.assert_allclose(
pdme.util.fast_v_calc.fast_vs_for_dipoles(dot_inputs, dipoles),
expected,
err_msg="Voltages at dot aren't as expected.",
)
def test_between(): def test_between():
low = numpy.array([1, 2, 3]) low = numpy.array([1, 2, 3])
high = numpy.array([6, 7, 8]) high = numpy.array([6, 7, 8])
# FALSE FALSE TRUE # FALSE FALSE TRUE
a = [[0, 1, 2], [0, 9, 5], [4, 5, 6]] a = [[0, 1, 2], [0, 9, 5], [4, 5, 6]]
actual = pdme.util.fast_v_calc.between(a, low, high) actual = pdme.util.fast_v_calc.between(a, low, high)