CHANGELOG.md

@@ -2,6 +2,50 @@
 
 All notable changes to this project will be documented in this file. See [standard-version](https://github.com/conventional-changelog/standard-version) for commit guidelines.
 
+### [0.8.1](https://gitea.deepak.science:2222/physics/pdme/compare/0.8.0...0.8.1) (2022-04-30)
+
+
+### Features
+
+* adds multidipole nonlocal spectrum calculations ([8264ca3](https://gitea.deepak.science:2222/physics/pdme/commit/8264ca3d6edb703422229fde57bbb1d726ce5139))
+
+## [0.8.0](https://gitea.deepak.science:2222/physics/pdme/compare/0.7.0...0.8.0) (2022-04-30)
+
+
+### ⚠ BREAKING CHANGES
+
+* single dipole still outputs collection now to make interface consistent
+
+### Features
+
+* single dipole still outputs collection now to make interface consistent ([cb3c280](https://gitea.deepak.science:2222/physics/pdme/commit/cb3c280464aa2c4b0ec9320e6a319f4a454a0e9f))
+
+## [0.7.0](https://gitea.deepak.science:2222/physics/pdme/compare/0.6.2...0.7.0) (2022-04-30)
+
+
+### ⚠ BREAKING CHANGES
+
+* Changes names of classes to make more clear for dipole model and single dipole fixed magnitude model
+* reduces model to minimal bayes needed stuff
+* Guts the model interface for only the things useful for monte carlo bayes stuff
+* Removes unused models to make refactoring a bit easier
+
+### Features
+
+* adds fast method for calculating multiple dipole calculations ([1bad2f7](https://gitea.deepak.science:2222/physics/pdme/commit/1bad2f743af6a95189448d8929c70a036e6c21ab))
+* adds multidipole fixed mag model ([0c64ed0](https://gitea.deepak.science:2222/physics/pdme/commit/0c64ed02f0d42cd0957be0511a35dc49153a256f))
+* Changes names of classes to make more clear for dipole model and single dipole fixed magnitude model ([7eba331](https://gitea.deepak.science:2222/physics/pdme/commit/7eba3311c7ba90e4404b1a7b7c33df1585a800ba))
+* Guts the model interface for only the things useful for monte carlo bayes stuff ([946945d](https://gitea.deepak.science:2222/physics/pdme/commit/946945d791ed763e435afff7a6ae8c2b1c0e1711))
+* reduces model to minimal bayes needed stuff ([0d5508a](https://gitea.deepak.science:2222/physics/pdme/commit/0d5508a0b5641d98ee6f3484f58c923440c4c2c1))
+* Removes unused models to make refactoring a bit easier ([c04d863](https://gitea.deepak.science:2222/physics/pdme/commit/c04d863d7fa22c98d8542a5f72791b542358ff61))
+
+
+### Bug Fixes
+
+* makes name of method match interface ([029021e](https://gitea.deepak.science:2222/physics/pdme/commit/029021e39328c7ca7a73afcf23d86ad17516982f))
+* makes repr return actual name ([0f78c7c](https://gitea.deepak.science:2222/physics/pdme/commit/0f78c7c2db742d59f2a5ee9da767f35aafa49b78))
+* uses rng passed in correctly ([70d95c8](https://gitea.deepak.science:2222/physics/pdme/commit/70d95c8d6d5af10f926d8729a5a0eaed9ad8259d))
+
 ### [0.6.2](https://gitea.deepak.science:2222/physics/pdme/compare/0.6.1...0.6.2) (2022-04-18)
 
 

pdme/model/__init__.py

@@ -1,14 +1,11 @@
-from pdme.model.model import Model, Discretisation
+from pdme.model.model import DipoleModel
-from pdme.model.fixed_z_plane_model import FixedZPlaneModel
+from pdme.model.fixed_magnitude_model import SingleDipoleFixedMagnitudeModel
-from pdme.model.unrestricted_model import UnrestrictedModel
+from pdme.model.multidipole_fixed_magnitude_model import (
-from pdme.model.fixed_dipole_model import FixedDipoleModel
+	MultipleDipoleFixedMagnitudeModel,
-from pdme.model.fixed_magnitude_model import FixedMagnitudeModel
+)
 
 __all__ = [
-	"Model",
+	"DipoleModel",
-	"Discretisation",
+	"SingleDipoleFixedMagnitudeModel",
-	"FixedZPlaneModel",
+	"MultipleDipoleFixedMagnitudeModel",
-	"UnrestrictedModel",
-	"FixedDipoleModel",
-	"FixedMagnitudeModel",
 ]

pdme/model/fixed_dipole_model.py

@@ -1,184 +0,0 @@
-import numpy
-import numpy.random
-from dataclasses import dataclass
-from typing import Sequence, Tuple
-import scipy.optimize
-from pdme.model.model import Model, Discretisation
-from pdme.measurement import (
-	DotMeasurement,
-	OscillatingDipoleArrangement,
-	OscillatingDipole,
-)
-
-
-class FixedDipoleModel(Model):
-	"""
-	Model of oscillating dipole with a fixed dipole moment.
-
-	Parameters
-	----------
-	p : numpy.ndarray
-	The fixed dipole moment.
-
-	n : int
-	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:
-		self.xmin = xmin
-		self.xmax = xmax
-		self.ymin = ymin
-		self.ymax = ymax
-		self.zmin = zmin
-		self.zmax = zmax
-		self.p = p
-		self._n = n
-		self.rng = numpy.random.default_rng()
-
-	def __repr__(self) -> str:
-		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.
-	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),
-			)
-		)
-		return OscillatingDipoleArrangement(
-			[OscillatingDipole(self.p, s_pts, frequency)]
-		)
-
-	def solution_single_dipole(self, pt: numpy.ndarray) -> OscillatingDipole:
-		# assume length is 4.
-		s = pt[0:3]
-		w = pt[3]
-		return OscillatingDipole(self.p, s, w)
-
-	def point_length(self) -> int:
-		"""
-		Dipole is constrained magnitude, but free orientation.
-		Six degrees of freedom: (sx, sy, sz, w).
-		"""
-		return 4
-
-	def n(self) -> int:
-		return self._n
-
-	def v_for_point_at_dot(self, dot: DotMeasurement, pt: numpy.ndarray) -> float:
-		s = pt[0:3]
-		w = pt[3]
-
-		diff = dot.r - s
-		alpha = self.p.dot(diff) / (numpy.linalg.norm(diff) ** 3)
-		b = (1 / numpy.pi) * (w / (w**2 + dot.f**2))
-		return alpha**2 * b
-
-	def jac_for_point_at_dot(
-		self, dot: DotMeasurement, pt: numpy.ndarray
-	) -> numpy.ndarray:
-		s = pt[0:3]
-		w = pt[3]
-
-		diff = dot.r - s
-		alpha = self.p.dot(diff) / (numpy.linalg.norm(diff) ** 3)
-		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
-		)
-
-		f2 = dot.f**2
-		w2 = w**2
-
-		w_div = alpha**2 * (1 / numpy.pi) * ((f2 - w2) / ((f2 + w2) ** 2))
-
-		return numpy.concatenate((r_divs, w_div), axis=None)
-
-
-@dataclass
-class FixedDipoleDiscretisation(Discretisation):
-	"""
-	Representation of a discretisation of a FixedDipoleDiscretisation.
-	Also captures a rough maximum value of dipole.
-
-	Parameters
-	----------
-	model : FixedDipoleModel
-	The parent model of the discretisation.
-
-	num_x : int
-	The number of partitions of the x axis.
-
-	num_y : int
-	The number of partitions of the y axis.
-
-	num_z : int
-	The number of partitions of the z axis.
-	"""
-
-	model: FixedDipoleModel
-	num_x: int
-	num_y: int
-	num_z: int
-
-	def __post_init__(self):
-		self.cell_count = self.num_x * self.num_y * self.num_z
-		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.z_step = (self.model.zmax - self.model.zmin) / self.num_z
-
-	def bounds(self, index: Tuple[float, ...]) -> Tuple:
-		xi, yi, zi = index
-
-		# For this model, a point is (sx, sx, sy, w).
-		# We want to keep w unbounded, restrict sx, sy, sz, px and py based on step.
-		return (
-			[
-				xi * self.x_step + self.model.xmin,
-				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,
-				numpy.inf,
-			],
-		)
-
-	def all_indices(self) -> numpy.ndindex:
-		# 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
-
-	def solve_for_index(
-		self, dots: Sequence[DotMeasurement], index: Tuple[float, ...]
-	) -> scipy.optimize.OptimizeResult:
-		bounds = self.bounds(index)
-		sx_mean = (bounds[0][0] + bounds[1][0]) / 2
-		sy_mean = (bounds[0][1] + bounds[1][1]) / 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, 0.1]),
-			bounds=bounds,
-		)

pdme/model/fixed_magnitude_model.py

@@ -1,27 +1,20 @@
 import numpy
 import numpy.random
-from dataclasses import dataclass
+from pdme.model.model import DipoleModel
-from typing import Sequence, Tuple
-import scipy.optimize
-from pdme.model.model import Model, Discretisation
 from pdme.measurement import (
-	DotMeasurement,
 	OscillatingDipole,
 	OscillatingDipoleArrangement,
 )
 
 
-class FixedMagnitudeModel(Model):
+class SingleDipoleFixedMagnitudeModel(DipoleModel):
 	"""
-	Model of oscillating dipole with a fixed magnitude, but free rotation.
+	Model of single oscillating dipole with a fixed magnitude, but free rotation.
 
 	Parameters
 	----------
 	pfixed : float
 	The fixed dipole magnitude.
-
-	n : int
-	The number of dipoles to assume.
 	"""
 
 	def __init__(
@@ -33,7 +26,6 @@ class FixedMagnitudeModel(Model):
 		zmin: float,
 		zmax: float,
 		pfixed: float,
-		n: int,
 	) -> None:
 		self.xmin = xmin
 		self.xmax = xmax
@@ -42,56 +34,40 @@ class FixedMagnitudeModel(Model):
 		self.zmin = zmin
 		self.zmax = zmax
 		self.pfixed = pfixed
-		self._n = n
 		self.rng = numpy.random.default_rng()
 
 	def __repr__(self) -> str:
-		return f"FixedMagnitudeModel({self.xmin}, {self.xmax}, {self.ymin}, {self.ymax}, {self.zmin}, {self.zmax}, {self.pfixed}, {self.n()})"
+		return f"SingleDipoleFixedMagnitudeModel({self.xmin}, {self.xmax}, {self.ymin}, {self.ymax}, {self.zmin}, {self.zmax}, {self.pfixed})"
 
-	def solution_single_dipole(self, pt: numpy.ndarray) -> OscillatingDipole:
+	def get_dipoles(
-		# assume length is 6, who needs error checking.
+		self, max_frequency: float, rng_to_use: numpy.random.Generator = None
-		p_theta = pt[0]
+	) -> OscillatingDipoleArrangement:
-		p_phi = pt[1]
+		rng: numpy.random.Generator
-		s = pt[2:5]
+		if rng_to_use is None:
-		w = pt[5]
+			rng = self.rng
+		else:
+			rng = rng_to_use
 
-		p = numpy.array(
+		theta = numpy.arccos(rng.uniform(-1, 1))
-			[
+		phi = rng.uniform(0, 2 * numpy.pi)
-				self.pfixed * numpy.sin(p_theta) * numpy.cos(p_phi),
-				self.pfixed * numpy.sin(p_theta) * numpy.sin(p_phi),
-				self.pfixed * numpy.cos(p_theta),
-			]
-		)
-		return OscillatingDipole(p, s, w)
-
-	def point_length(self) -> int:
-		"""
-		Dipole is constrained magnitude, but free orientation.
-		Six degrees of freedom: (p_theta, p_phi, sx, sy, sz, w).
-		"""
-		return 6
-
-	def get_dipoles(self, frequency: float) -> OscillatingDipoleArrangement:
-		theta = numpy.arccos(self.rng.uniform(-1, 1))
-		phi = self.rng.uniform(0, 2 * numpy.pi)
 		px = self.pfixed * numpy.sin(theta) * numpy.cos(phi)
 		py = self.pfixed * numpy.sin(theta) * numpy.sin(phi)
 		pz = self.pfixed * numpy.cos(theta)
 		s_pts = numpy.array(
 			(
-				self.rng.uniform(self.xmin, self.xmax),
+				rng.uniform(self.xmin, self.xmax),
-				self.rng.uniform(self.ymin, self.ymax),
+				rng.uniform(self.ymin, self.ymax),
-				self.rng.uniform(self.zmin, self.zmax),
+				rng.uniform(self.zmin, self.zmax),
 			)
 		)
+		frequency = rng.uniform(0, max_frequency)
 		return OscillatingDipoleArrangement(
 			[OscillatingDipole(numpy.array([px, py, pz]), s_pts, frequency)]
 		)
 
-	def get_n_single_dipoles(
+	def get_monte_carlo_dipole_inputs(
 		self, n: int, max_frequency: float, rng_to_use: numpy.random.Generator = None
 	) -> numpy.ndarray:
-		# psw
 
 		rng: numpy.random.Generator
 		if rng_to_use is None:
@@ -99,179 +75,17 @@ class FixedMagnitudeModel(Model):
 		else:
 			rng = rng_to_use
 
-		theta = 2 * numpy.pi * rng.random(n)
+		shape = (n, 1)
-		phi = numpy.arccos(2 * rng.random(n) - 1)
+		theta = 2 * numpy.pi * rng.random(shape)
+		phi = numpy.arccos(2 * rng.random(shape) - 1)
 		px = self.pfixed * numpy.cos(theta) * numpy.sin(phi)
 		py = self.pfixed * numpy.sin(theta) * numpy.sin(phi)
 		pz = self.pfixed * numpy.cos(phi)
 
-		sx = rng.uniform(self.xmin, self.xmax, n)
+		sx = rng.uniform(self.xmin, self.xmax, shape)
-		sy = rng.uniform(self.ymin, self.ymax, n)
+		sy = rng.uniform(self.ymin, self.ymax, shape)
-		sz = rng.uniform(self.zmin, self.zmax, n)
+		sz = rng.uniform(self.zmin, self.zmax, shape)
 
-		w = rng.uniform(1, max_frequency, n)
+		w = rng.uniform(1, max_frequency, shape)
 
-		return numpy.array([px, py, pz, sx, sy, sz, w]).T
+		return numpy.stack([px, py, pz, sx, sy, sz, w], axis=-1)
-
-	def n(self) -> int:
-		return self._n
-
-	def v_for_point_at_dot(self, dot: DotMeasurement, pt: numpy.ndarray) -> float:
-		p_theta = pt[0]
-		p_phi = pt[1]
-		s = pt[2:5]
-		w = pt[5]
-
-		p = numpy.array(
-			[
-				self.pfixed * numpy.sin(p_theta) * numpy.cos(p_phi),
-				self.pfixed * numpy.sin(p_theta) * numpy.sin(p_phi),
-				self.pfixed * numpy.cos(p_theta),
-			]
-		)
-		diff = dot.r - s
-		alpha = p.dot(diff) / (numpy.linalg.norm(diff) ** 3)
-		b = (1 / numpy.pi) * (w / (w**2 + dot.f**2))
-		return alpha**2 * b
-
-	def jac_for_point_at_dot(
-		self, dot: DotMeasurement, pt: numpy.ndarray
-	) -> numpy.ndarray:
-		p_theta = pt[0]
-		p_phi = pt[1]
-		s = pt[2:5]
-		w = pt[5]
-
-		p = numpy.array(
-			[
-				self.pfixed * numpy.sin(p_theta) * numpy.cos(p_phi),
-				self.pfixed * numpy.sin(p_theta) * numpy.sin(p_phi),
-				self.pfixed * numpy.cos(p_theta),
-			]
-		)
-		diff = dot.r - s
-		alpha = p.dot(diff) / (numpy.linalg.norm(diff) ** 3)
-		b = (1 / numpy.pi) * (w / (w**2 + dot.f**2))
-
-		theta_div_middle = self.pfixed * (
-			diff[0] * numpy.cos(p_phi) * numpy.cos(p_theta)
-			+ diff[1] * numpy.sin(p_phi) * numpy.cos(p_theta)
-			- diff[2] * numpy.sin(p_theta)
-		)
-		theta_div = 2 * alpha * (theta_div_middle) / (numpy.linalg.norm(diff) ** 3) * b
-
-		phi_div_middle = self.pfixed * (
-			diff[1] * numpy.sin(p_theta) * numpy.cos(p_phi)
-			- diff[0] * numpy.sin(p_theta) * numpy.sin(p_phi)
-		)
-		phi_div = 2 * alpha * (phi_div_middle) / (numpy.linalg.norm(diff) ** 3) * b
-
-		r_divs = (
-			(
-				-p / (numpy.linalg.norm(diff) ** 3)
-				+ 3 * p.dot(diff) * diff / (numpy.linalg.norm(diff) ** 5)
-			)
-			* 2
-			* alpha
-			* b
-		)
-
-		f2 = dot.f**2
-		w2 = w**2
-
-		w_div = alpha**2 * (1 / numpy.pi) * ((f2 - w2) / ((f2 + w2) ** 2))
-
-		return numpy.concatenate((theta_div, phi_div, r_divs, w_div), axis=None)
-
-
-@dataclass
-class FixedMagnitudeDiscretisation(Discretisation):
-	"""
-	Representation of a discretisation of a FixedMagnitudeDiscretisation.
-	Also captures a rough maximum value of dipole.
-
-	Parameters
-	----------
-	model : FixedMagnitudeModel
-	The parent model of the discretisation.
-
-	num_ptheta: int
-	The number of partitions of ptheta.
-
-	num_pphi: int
-	The number of partitions of pphi.
-
-	num_x : int
-	The number of partitions of the x axis.
-
-	num_y : int
-	The number of partitions of the y axis.
-
-	num_z : int
-	The number of partitions of the z axis.
-	"""
-
-	model: FixedMagnitudeModel
-	num_ptheta: int
-	num_pphi: int
-	num_x: int
-	num_y: int
-	num_z: int
-
-	def __post_init__(self):
-		self.cell_count = self.num_x * self.num_y * self.num_z
-		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.z_step = (self.model.zmax - self.model.zmin) / self.num_z
-		self.h_step = 2 / self.num_ptheta
-		self.phi_step = 2 * numpy.pi / self.num_pphi
-
-	def bounds(self, index: Tuple[float, ...]) -> Tuple:
-		pthetai, pphii, xi, yi, zi = index
-
-		# For this model, a point is (p_theta, p_phi, sx, sx, sy, w).
-		# We want to keep w unbounded, restrict sx, sy, sz, px and py based on step.
-		return (
-			[
-				numpy.arccos(1 - pthetai * self.h_step),
-				pphii * self.phi_step,
-				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,
-				(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:
-		return self.model
-
-	def all_indices(self) -> numpy.ndindex:
-		# 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
-
-	def solve_for_index(
-		self, dots: Sequence[DotMeasurement], index: Tuple[float, ...]
-	) -> scipy.optimize.OptimizeResult:
-		bounds = self.bounds(index)
-		ptheta_mean = (bounds[0][0] + bounds[1][0]) / 2
-		pphi_mean = (bounds[0][1] + bounds[1][1]) / 2
-		sx_mean = (bounds[0][2] + bounds[1][2]) / 2
-		sy_mean = (bounds[0][3] + bounds[1][3]) / 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, 0.1]
-			),
-			bounds=bounds,
-		)

pdme/model/fixed_z_plane_model.py

@@ -1,173 +0,0 @@
-import numpy
-from dataclasses import dataclass
-from typing import Tuple, Sequence
-import scipy.optimize
-
-from pdme.model.model import Model
-from pdme.measurement import DotMeasurement
-
-
-class FixedZPlaneModel(Model):
-	"""
-	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.
-
-	Parameters
-	----------
-	z : float
-	The z position of the plane where dipoles are constrained to lie.
-
-	xmin : float
-	The minimum x value for dipoles.
-
-	xmax : float
-	The maximum x value for dipoles.
-
-	ymin : float
-	The minimum y value for dipoles.
-
-	ymax : float
-	The maximum y value for dipoles.
-
-	n : int
-	The number of dipoles to assume.
-	"""
-
-	def __init__(
-		self, z: float, xmin: float, xmax: float, ymin: float, ymax: float, n: int
-	) -> None:
-		self.z = z
-		self.xmin = xmin
-		self.xmax = xmax
-		self.ymin = ymin
-		self.ymax = ymax
-		self._n = n
-
-	def __repr__(self) -> str:
-		return f"FixedZPlaneModel({self.z}, {self.xmin}, {self.xmax}, {self.ymin}, {self.ymax}, {self.n()})"
-
-	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).
-		With some frequency w, there are four degrees of freedom: (pz, sx, sy, w).
-		"""
-		return 4
-
-	def n(self) -> int:
-		return self._n
-
-	def v_for_point_at_dot(self, dot: DotMeasurement, pt: numpy.ndarray) -> float:
-		p = numpy.array([0, 0, pt[0]])
-		s = numpy.array([pt[1], pt[2], self.z])
-		w = pt[3]
-
-		diff = dot.r - s
-		alpha = p.dot(diff) / (numpy.linalg.norm(diff) ** 3)
-		b = (1 / numpy.pi) * (w / (w**2 + dot.f**2))
-		return alpha**2 * b
-
-	def jac_for_point_at_dot(
-		self, dot: DotMeasurement, pt: numpy.ndarray
-	) -> numpy.ndarray:
-		p = numpy.array([0, 0, pt[0]])
-		s = numpy.array([pt[1], pt[2], self.z])
-		w = pt[3]
-
-		diff = dot.r - s
-		alpha = p.dot(diff) / (numpy.linalg.norm(diff) ** 3)
-		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.
-
-		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
-		w2 = w**2
-
-		w_div = alpha**2 * (1 / numpy.pi) * ((f2 - w2) / ((f2 + w2) ** 2))
-
-		return numpy.concatenate((p_divs, r_divs, w_div), axis=None)
-
-
-@dataclass
-class FixedZPlaneDiscretisation:
-	"""
-	Representation of a discretisation of a FixedZPlaneModel.
-	Also captures a rough maximum value of dipole.
-
-	Parameters
-	----------
-	model : FixedZPlaneModel
-	The parent model of the discretisation.
-
-	num_pz: int
-	The number of partitions of pz.
-
-	num_x : int
-	The number of partitions of the x axis.
-
-	num_y : int
-	The number of partitions of the y axis.
-	"""
-
-	model: FixedZPlaneModel
-	num_pz: int
-	num_x: int
-	num_y: int
-	max_pz: int
-
-	def __post_init__(self):
-		self.cell_count = self.num_x * self.num_y
-		self.pz_step = (2 * self.max_pz) / self.num_pz
-		self.x_step = (self.model.xmax - self.model.xmin) / self.num_x
-		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]:
-		pzi, xi, yi = index
-
-		# 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.
-		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(
-				(
-					(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:
-		# 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
-
-	def solve_for_index(
-		self, dots: Sequence[DotMeasurement], index: Tuple[float, float, float]
-	) -> scipy.optimize.OptimizeResult:
-		bounds = self.bounds(index)
-		pz_mean = (bounds[0][0] + bounds[1][0]) / 2
-		sx_mean = (bounds[0][1] + bounds[1][1]) / 2
-		sy_mean = (bounds[0][2] + bounds[1][2]) / 2
-		# 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, 0.1)), bounds=bounds)  # type: ignore

pdme/model/model.py

@@ -1,154 +1,34 @@
 import numpy
 import numpy.random
-import scipy.optimize
-from typing import Callable, Sequence, Tuple, List
 from pdme.measurement import (
-	DotMeasurement,
 	OscillatingDipoleArrangement,
-	OscillatingDipole,
 )
-import pdme.util
 import logging
 
 
 _logger = logging.getLogger(__name__)
 
 
-class Model:
+class DipoleModel:
 	"""
-	Interface for models.
+	Interface for models based on dipoles.
+	Some concepts are kept specific for dipole-based models, even though other types of models could be useful later on.
 	"""
 
-	def point_length(self) -> int:
+	def get_dipoles(
+		self, max_frequency: float, rng: numpy.random.Generator = None
+	) -> OscillatingDipoleArrangement:
+		"""
+		For a particular maximum frequency, gets a dipole arrangement based on the model that uniformly distributes its choices according to the model.
+		If no rng is passed in, uses some default, but you might not want that.
+		Frequencies should be chosen uniformly on range of (0, max_frequency).
+		"""
 		raise NotImplementedError
 
-	def n(self) -> int:
+	def get_monte_carlo_dipole_inputs(
-		raise NotImplementedError
-
-	def v_for_point_at_dot(self, dot: DotMeasurement, pt: numpy.ndarray) -> float:
-		raise NotImplementedError
-
-	def get_dipoles(self, frequency: float) -> OscillatingDipoleArrangement:
-		raise NotImplementedError
-
-	def get_n_single_dipoles(
 		self, n: int, max_frequency: float, rng: numpy.random.Generator = None
 	) -> numpy.ndarray:
-		raise NotImplementedError
-
-	def solution_single_dipole(self, pt: numpy.ndarray) -> OscillatingDipole:
-		raise NotImplementedError
-
-	def solution_as_dipoles(self, pts: numpy.ndarray) -> List[OscillatingDipole]:
-		pt_length = self.point_length()
-		chunked_pts = [pts[i : i + pt_length] for i in range(0, len(pts), pt_length)]
-		return [self.solution_single_dipole(pt) for pt in chunked_pts]
-
-	def cost_for_dot(self, dot: DotMeasurement, pts: numpy.ndarray) -> float:
-		# creates numpy.ndarrays in groups of self.point_length().
-		# Will throw problems for irregular points, but that's okay for now.
-		pt_length = self.point_length()
-		chunked_pts = [pts[i : i + pt_length] for i in range(0, len(pts), pt_length)]
-		return sum(self.v_for_point_at_dot(dot, pt) for pt in chunked_pts) - dot.v
-
-	def costs(
-		self, dots: Sequence[DotMeasurement]
-	) -> Callable[[numpy.ndarray], numpy.ndarray]:
 		"""
-		Returns a function that returns the cost for the given list of DotMeasurements for a particular model-dependent phase space point.
+		For a given DipoleModel, gets a set of dipole collections as a monte_carlo_n x dipole_count x 7 numpy array.
-		Default implementation assumes a single dot cost from which to build the list.
-
-		Parameters
-		----------
-		dots: A list of dot measurements to use to find the cost functions.
-
-		Returns
-		----------
-		Returns the model's cost function.
 		"""
-		_logger.debug(f"Constructing costs for dots: {dots}")
-
-		def costs_to_return(pts: numpy.ndarray) -> numpy.ndarray:
-			return numpy.array([self.cost_for_dot(dot, pts) for dot in dots])
-
-		return costs_to_return
-
-	def jac_for_point_at_dot(
-		self, dot: DotMeasurement, pt: numpy.ndarray
-	) -> numpy.ndarray:
-		raise NotImplementedError
-
-	def jac_for_dot(self, dot: DotMeasurement, pts: numpy.ndarray) -> numpy.ndarray:
-		# creates numpy.ndarrays in groups of self.point_length().
-		# Will throw problems for irregular points, but that's okay for now.
-		pt_length = self.point_length()
-		chunked_pts = [pts[i : i + pt_length] for i in range(0, len(pts), pt_length)]
-		return numpy.append(
-			[], [self.jac_for_point_at_dot(dot, pt) for pt in chunked_pts]
-		)
-
-	def jac(
-		self, dots: Sequence[DotMeasurement]
-	) -> Callable[[numpy.ndarray], numpy.ndarray]:
-		"""
-		Returns a function that returns the cost function's Jacobian for the given list of DotMeasurements for a particular model-dependent phase space point.
-		Default implementation assumes a single dot jacobian from which to build the list.
-
-		Parameters
-		----------
-		dots: A list of dot measurements to use to find the cost functions and their Jacobian.
-
-		Returns
-		----------
-		Returns the model's cost function's Jacobian.
-		"""
-
-		def jac_to_return(pts: numpy.ndarray) -> numpy.ndarray:
-			return numpy.array([self.jac_for_dot(dot, pts) for dot in dots])
-
-		return jac_to_return
-
-	def solve(
-		self,
-		dots: Sequence[DotMeasurement],
-		initial_pt: numpy.ndarray = None,
-		bounds=(-numpy.inf, numpy.inf),
-	) -> scipy.optimize.OptimizeResult:
-		if initial_pt is None:
-			initial = numpy.tile(0.1, self.n() * self.point_length())
-		else:
-			if len(initial_pt) != self.point_length():
-				raise ValueError(
-					f"The initial point {initial_pt} does not have the model's expected length: {self.point_length()}"
-				)
-			initial = numpy.tile(initial_pt, self.n())
-
-		result = scipy.optimize.least_squares(
-			self.costs(dots),
-			initial,
-			jac=self.jac(dots),
-			ftol=1e-15,
-			gtol=3e-16,
-			xtol=None,
-			bounds=bounds,
-		)
-		result.normalised_x = pdme.util.normalise_point_list(
-			result.x, self.point_length()
-		)
-		return result
-
-
-class Discretisation:
-	def bounds(self, index: Tuple[float, ...]) -> Tuple:
-		raise NotImplementedError
-
-	def all_indices(self) -> numpy.ndindex:
-		raise NotImplementedError
-
-	def solve_for_index(
-		self, dots: Sequence[DotMeasurement], index: Tuple
-	) -> scipy.optimize.OptimizeResult:
-		raise NotImplementedError
-
-	def get_model(self) -> Model:
 		raise NotImplementedError

pdme/model/multidipole_fixed_magnitude_model.py

@@ -0,0 +1,104 @@
+import numpy
+import numpy.random
+from pdme.model.model import DipoleModel
+from pdme.measurement import (
+	OscillatingDipole,
+	OscillatingDipoleArrangement,
+)
+
+
+class MultipleDipoleFixedMagnitudeModel(DipoleModel):
+	"""
+	Model of multiple oscillating dipoles with a fixed magnitude, but free rotation.
+
+	Parameters
+	----------
+	pfixed : float
+	The fixed dipole magnitude.
+
+	n : int
+	The number of dipoles.
+	"""
+
+	def __init__(
+		self,
+		xmin: float,
+		xmax: float,
+		ymin: float,
+		ymax: float,
+		zmin: float,
+		zmax: float,
+		pfixed: float,
+		n: int,
+	) -> None:
+		self.xmin = xmin
+		self.xmax = xmax
+		self.ymin = ymin
+		self.ymax = ymax
+		self.zmin = zmin
+		self.zmax = zmax
+		self.pfixed = pfixed
+		self.rng = numpy.random.default_rng()
+		self.n = n
+
+	def __repr__(self) -> str:
+		return f"MultipleDipoleFixedMagnitudeModel({self.xmin}, {self.xmax}, {self.ymin}, {self.ymax}, {self.zmin}, {self.zmax}, {self.pfixed}, {self.n})"
+
+	def get_dipoles(
+		self, max_frequency: float, rng_to_use: numpy.random.Generator = None
+	) -> OscillatingDipoleArrangement:
+		rng: numpy.random.Generator
+		if rng_to_use is None:
+			rng = self.rng
+		else:
+			rng = rng_to_use
+
+		dipoles = []
+
+		for i in range(self.n):
+			theta = numpy.arccos(rng.uniform(-1, 1))
+			phi = rng.uniform(0, 2 * numpy.pi)
+			px = self.pfixed * numpy.sin(theta) * numpy.cos(phi)
+			py = self.pfixed * numpy.sin(theta) * numpy.sin(phi)
+			pz = self.pfixed * numpy.cos(theta)
+			s_pts = numpy.array(
+				(
+					rng.uniform(self.xmin, self.xmax),
+					rng.uniform(self.ymin, self.ymax),
+					rng.uniform(self.zmin, self.zmax),
+				)
+			)
+			frequency = rng.uniform(0, max_frequency)
+
+			dipoles.append(
+				OscillatingDipole(numpy.array([px, py, pz]), s_pts, frequency)
+			)
+		return OscillatingDipoleArrangement(dipoles)
+
+	def get_monte_carlo_dipole_inputs(
+		self,
+		monte_carlo_n: int,
+		max_frequency: float,
+		rng_to_use: numpy.random.Generator = None,
+	) -> numpy.ndarray:
+
+		rng: numpy.random.Generator
+		if rng_to_use is None:
+			rng = self.rng
+		else:
+			rng = rng_to_use
+
+		shape = (monte_carlo_n, self.n)
+		theta = 2 * numpy.pi * rng.random(shape)
+		phi = numpy.arccos(2 * rng.random(shape) - 1)
+		px = self.pfixed * numpy.cos(theta) * numpy.sin(phi)
+		py = self.pfixed * numpy.sin(theta) * numpy.sin(phi)
+		pz = self.pfixed * numpy.cos(phi)
+
+		sx = rng.uniform(self.xmin, self.xmax, shape)
+		sy = rng.uniform(self.ymin, self.ymax, shape)
+		sz = rng.uniform(self.zmin, self.zmax, shape)
+
+		w = rng.uniform(1, max_frequency, shape)
+
+		return numpy.stack([px, py, pz, sx, sy, sz, w], axis=-1)

pdme/model/unrestricted_model.py

@@ -1,216 +0,0 @@
-import numpy
-from dataclasses import dataclass
-from typing import Sequence, Tuple
-import scipy.optimize
-from pdme.model.model import Model, Discretisation
-from pdme.measurement import (
-	DotMeasurement,
-	OscillatingDipoleArrangement,
-	OscillatingDipole,
-)
-
-
-class UnrestrictedModel(Model):
-	"""
-	Model of oscillating dipoles with no restrictions.
-	Additionally, each dipole is assumed to be orientated in the plus or minus z direction.
-
-	Parameters
-	----------
-	n : int
-	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:
-		self.xmin = xmin
-		self.xmax = xmax
-		self.ymin = ymin
-		self.ymax = ymax
-		self.zmin = zmin
-		self.zmax = zmax
-		self.max_p = max_p
-		self._n = n
-		self.rng = numpy.random.default_rng()
-
-	def __repr__(self) -> str:
-		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:
-		"""
-		Dipole is unconstrained in this model.
-		All seven degrees of freedom: (px, py, pz, sx, sy, sz, w).
-		"""
-		return 7
-
-	def n(self) -> int:
-		return self._n
-
-	def v_for_point_at_dot(self, dot: DotMeasurement, pt: numpy.ndarray) -> float:
-		p = pt[0:3]
-		s = pt[3:6]
-		w = pt[6]
-
-		diff = dot.r - s
-		alpha = p.dot(diff) / (numpy.linalg.norm(diff) ** 3)
-		b = (1 / numpy.pi) * (w / (w**2 + dot.f**2))
-		return alpha**2 * b
-
-	def jac_for_point_at_dot(
-		self, dot: DotMeasurement, pt: numpy.ndarray
-	) -> numpy.ndarray:
-		p = pt[0:3]
-		s = pt[3:6]
-		w = pt[6]
-
-		diff = dot.r - s
-		alpha = p.dot(diff) / (numpy.linalg.norm(diff) ** 3)
-		b = (1 / numpy.pi) * (w / (w**2 + dot.f**2))
-
-		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
-		)
-
-		f2 = dot.f**2
-		w2 = w**2
-
-		w_div = alpha**2 * (1 / numpy.pi) * ((f2 - w2) / ((f2 + w2) ** 2))
-
-		return numpy.concatenate((p_divs, r_divs, w_div), axis=None)
-
-	def get_dipoles(self, frequency: float) -> OscillatingDipoleArrangement:
-		theta = numpy.arccos(self.rng.uniform(-1, 1))
-		phi = self.rng.uniform(0, 2 * numpy.pi)
-		p = self.rng.uniform(0, self.max_p)
-		px = p * numpy.sin(theta) * numpy.cos(phi)
-		py = p * numpy.sin(theta) * numpy.sin(phi)
-		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),
-			)
-		)
-		return OscillatingDipoleArrangement(
-			[OscillatingDipole(numpy.array([px, py, pz]), s_pts, frequency)]
-		)
-
-
-@dataclass
-class UnrestrictedDiscretisation(Discretisation):
-	"""
-	Representation of a discretisation of a UnrestrictedModel.
-	Also captures a rough maximum value of dipole.
-
-	Parameters
-	----------
-	model : UnrestrictedModel
-	The parent model of the discretisation.
-
-	num_px: int
-	The number of partitions of the px.
-
-	num_py: int
-	The number of partitions of the py.
-
-	num_pz: int
-	The number of partitions of pz.
-
-	num_x : int
-	The number of partitions of the x axis.
-
-	num_y : int
-	The number of partitions of the y axis.
-
-	num_z : int
-	The number of partitions of the z axis.
-
-	max_p : int
-	The maximum p coordinate in any direction.
-	"""
-
-	model: UnrestrictedModel
-	num_px: int
-	num_py: int
-	num_pz: int
-	num_x: int
-	num_y: int
-	num_z: int
-
-	def __post_init__(self):
-		self.max_p = self.model.max_p
-		self.cell_count = self.num_x * self.num_y * self.num_z
-		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.z_step = (self.model.zmax - self.model.zmin) / self.num_z
-		self.px_step = 2 * self.max_p / self.num_px
-		self.py_step = 2 * self.max_p / self.num_py
-		self.pz_step = 2 * self.max_p / self.num_pz
-
-	def bounds(self, index: Tuple[float, ...]) -> Tuple:
-		pxi, pyi, pzi, xi, yi, zi = index
-
-		# For this model, a point is (px, py, pz, sx, sx, sy, w).
-		# We want to keep w unbounded, restrict sx, sy, sz, px and py based on step.
-		return (
-			[
-				pxi * self.px_step - self.max_p,
-				pyi * self.py_step - self.max_p,
-				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,
-				(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:
-		# 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
-
-	def solve_for_index(
-		self, dots: Sequence[DotMeasurement], index: Tuple[float, ...]
-	) -> scipy.optimize.OptimizeResult:
-		bounds = self.bounds(index)
-		px_mean = (bounds[0][0] + bounds[1][0]) / 2
-		py_mean = (bounds[0][1] + bounds[1][1]) / 2
-		pz_mean = (bounds[0][2] + bounds[1][2]) / 2
-		sx_mean = (bounds[0][3] + bounds[1][3]) / 2
-		sy_mean = (bounds[0][4] + bounds[1][4]) / 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, 0.1]
-			),
-			bounds=bounds,
-		)

pdme/util/__init__.py

@@ -1,3 +0,0 @@
-from pdme.util.normal_form import normalise_point_list
-
-__all__ = ["normalise_point_list"]

pdme/util/fast_nonlocal_spectrum.py

@@ -50,3 +50,51 @@ def fast_s_nonlocal(
 		_logger.debug(f"bses: {bses}")
 
 	return alphses1 * alphses2 * bses
+
+
+def fast_s_nonlocal_dipoleses(
+	dot_pair_inputs: numpy.ndarray, dipoleses: numpy.ndarray
+) -> numpy.ndarray:
+	"""
+	No error correction here baby.
+	"""
+	ps = dipoleses[:, :, 0:3]
+	ss = dipoleses[:, :, 3:6]
+	ws = dipoleses[:, :, 6]
+
+	_logger.debug(f"ps: {ps}")
+	_logger.debug(f"ss: {ss}")
+	_logger.debug(f"ws: {ws}")
+
+	r1s = dot_pair_inputs[:, 0, 0:3]
+	r2s = dot_pair_inputs[:, 1, 0:3]
+	f1s = dot_pair_inputs[:, 0, 3]
+	f2s = dot_pair_inputs[:, 1, 3]
+
+	if (f1s != f2s).all():
+		raise ValueError(f"Dot pair frequencies are inconsistent: {dot_pair_inputs}")
+
+	diffses1 = r1s[:, None] - ss[:, None, :]
+	diffses2 = r2s[:, None] - ss[:, None, :]
+	if _logger.isEnabledFor(logging.DEBUG):
+		_logger.debug(f"diffses1: {diffses1}")
+		_logger.debug(f"diffses2: {diffses2}")
+
+	norms1 = numpy.linalg.norm(diffses1, axis=3) ** 3
+	norms2 = numpy.linalg.norm(diffses2, axis=3) ** 3
+	if _logger.isEnabledFor(logging.DEBUG):
+		_logger.debug(f"norms1: {norms1}")
+		_logger.debug(f"norms2: {norms2}")
+
+	alphses1 = numpy.einsum("abcd,acd->abc", diffses1, ps) / norms1
+	alphses2 = numpy.einsum("abcd,acd->abc", diffses2, ps) / norms2
+	if _logger.isEnabledFor(logging.DEBUG):
+		_logger.debug(f"alphses1: {alphses1}")
+		_logger.debug(f"alphses2: {alphses2}")
+
+	bses = (1 / numpy.pi) * (ws[:, None, :] / (f1s[:, None] ** 2 + ws[:, None, :] ** 2))
+	if _logger.isEnabledFor(logging.DEBUG):
+		_logger.debug(f"bses: {bses}")
+
+	_logger.debug(f"Raw pair calc: [{alphses1 * alphses2 * bses}]")
+	return numpy.einsum("...j->...", alphses1 * alphses2 * bses)

pdme/util/fast_v_calc.py

@@ -10,6 +10,7 @@ def fast_vs_for_dipoles(
 ) -> numpy.ndarray:
 	"""
 	No error correction here baby.
+	Expects dot_inputs to be numpy array of [rx, ry, rz, f] entries, so a n by 4 where n is number of measurement points.
 	"""
 	ps = dipoles[:, 0:3]
 	ss = dipoles[:, 3:6]
@@ -35,6 +36,39 @@ def fast_vs_for_dipoles(
 	return ases * bses
 
 
+def fast_vs_for_dipoleses(
+	dot_inputs: numpy.ndarray, dipoleses: numpy.ndarray
+) -> numpy.ndarray:
+	"""
+	No error correction here baby.
+	Expects dot_inputs to be numpy array of [rx, ry, rz, f] entries, so a n by 4 where n is number of measurement points.
+
+	Dipoleses are expected to be array of arrays of arrays: list of sets of dipoles which are part of a single arrangement to be added together.
+	"""
+	ps = dipoleses[:, :, 0:3]
+	ss = dipoleses[:, :, 3:6]
+	ws = dipoleses[:, :, 6]
+
+	_logger.debug(f"ps: {ps}")
+	_logger.debug(f"ss: {ss}")
+	_logger.debug(f"ws: {ws}")
+
+	rs = dot_inputs[:, 0:3]
+	fs = dot_inputs[:, 3]
+
+	diffses = rs[:, None] - ss[:, None, :]
+
+	_logger.debug(f"diffses: {diffses}")
+	norms = numpy.linalg.norm(diffses, axis=3) ** 3
+	_logger.debug(f"norms: {norms}")
+	ases = (numpy.einsum("abcd,acd->abc", diffses, ps) / norms) ** 2
+	_logger.debug(f"ases: {ases}")
+
+	bses = (1 / numpy.pi) * (ws[:, None, :] / (fs[:, None] ** 2 + ws[:, None, :] ** 2))
+	_logger.debug(f"bses: {bses}")
+	return numpy.einsum("...j->...", ases * bses)
+
+
 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.

pdme/util/normal_form.py

@@ -1,45 +0,0 @@
-import numpy
-import operator
-import logging
-
-
-# flips px, py, pz
-SIGN_ARRAY_7 = numpy.array((-1, -1, -1, 1, 1, 1, 1))
-SIGN_ARRAY_4 = numpy.array((-1, 1, 1, 1))
-
-
-_logger = logging.getLogger(__name__)
-
-
-def flip_chunk_to_positive_px(pt: numpy.ndarray) -> numpy.ndarray:
-	if pt[0] > 0:
-		return pt
-	else:
-		# godawful hack.
-		if len(pt) == 7:
-			return SIGN_ARRAY_7 * pt
-		elif len(pt) == 4:
-			return SIGN_ARRAY_4 * pt
-		else:
-			_logger.warning(f"Could not normalise pt: {pt}. Returning as is...")
-			return pt
-
-
-def normalise_point_list(pts: numpy.ndarray, pt_length) -> numpy.ndarray:
-	chunked_pts = [
-		flip_chunk_to_positive_px(pts[i : i + pt_length])
-		for i in range(0, len(pts), pt_length)
-	]
-	range_to_length = list(range(pt_length))
-	rotated_range = (
-		range_to_length[pt_length - 1 :] + range_to_length[0 : pt_length - 1]
-	)
-	return numpy.concatenate(
-		sorted(
-			chunked_pts,
-			key=lambda x: tuple(
-				round(val, 3) for val in operator.itemgetter(*rotated_range)(x)
-			),
-		),
-		axis=None,
-	)

pyproject.toml

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "pdme"
-version = "0.6.2"
+version = "0.8.1"
 description = "Python dipole model evaluator"
 authors = ["Deepak <dmallubhotla+github@gmail.com>"]
 license = "GPL-3.0-only"

tests/model/test_fixed_dipole_moment_basic_solve.py

@@ -1,49 +0,0 @@
-from pdme.model import FixedDipoleModel
-from pdme.measurement import OscillatingDipole, OscillatingDipoleArrangement
-import logging
-import numpy
-import itertools
-import pytest
-
-
-@pytest.mark.skip(reason="No idea why this is failing, but it shouldn't!")
-def test_fixed_dipole_model_solve_basic():
-	# Initialise our dipole arrangement and create dot measurements along a square.
-	fixed_dipole_moment = numpy.array((2, 1, 2))
-	dipoles = OscillatingDipoleArrangement(
-		[OscillatingDipole(fixed_dipole_moment, (1, 2, 4), 6)]
-	)
-	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)
-
-	model = FixedDipoleModel(-3, 3, -3, 3, 3, 5, fixed_dipole_moment, 1)
-
-	# from the dipole, these are the unspecified variables in ((.5, 0, 2), (1, 2, 4), 8)
-	expected_solution = [1, 2, 4, 6]
-
-	result = model.solve(dots)
-	logging.info(result)
-	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,
-	)

tests/model/test_fixed_magnitude_basic_solve.py

@@ -1,74 +0,0 @@
-from pdme.model import FixedMagnitudeModel
-from pdme.measurement import OscillatingDipole, OscillatingDipoleArrangement
-import logging
-import numpy
-import itertools
-
-
-def test_fixed_magnitude_model_solve_basic():
-	# Initialise our dipole arrangement and create dot measurements along a square.
-	dipoles = OscillatingDipoleArrangement([OscillatingDipole((2, 0, 0), (1, 2, 4), 1)])
-	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)
-
-	model = FixedMagnitudeModel(-5, 5, -5, 5, -5, 5, 2, 1)
-
-	# from the dipole, these are the unspecified variables in ((0, 0, 2), (1, 2, 4), 1)
-	expected_solution = [numpy.pi / 2, 0, 1, 2, 4, 1]
-
-	result = model.solve(dots)
-	logging.info(result)
-	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,
-	)
-
-	solved_dipoles = model.solution_as_dipoles(result.normalised_x)
-
-	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(
-		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():
-	model = FixedMagnitudeModel(-5, 5, -5, 5, -5, 5, 2, 1)
-	assert repr(model) == "FixedMagnitudeModel(-5, 5, -5, 5, -5, 5, 2, 1)"

tests/model/test_fixed_z_basic_solve.py

@@ -1,54 +0,0 @@
-from pdme.model import FixedZPlaneModel
-from pdme.measurement import OscillatingDipole, OscillatingDipoleArrangement
-import logging
-import numpy
-import itertools
-import pytest
-
-
-def test_fixed_z_plane_model_solve_error_initial():
-
-	model = FixedZPlaneModel(4, -10, 10, -10, 10, 1)
-
-	with pytest.raises(ValueError):
-		model.solve([], initial_pt=[1, 2])
-
-
-def test_fixed_z_plane_model_solve_basic():
-	# Initialise our dipole arrangement and create dot measurements along a square.
-	dipoles = OscillatingDipoleArrangement([OscillatingDipole((0, 0, 2), (1, 2, 4), 1)])
-	dot_inputs = list(
-		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)
-
-	model = FixedZPlaneModel(4, -10, 10, -10, 10, 1)
-
-	# from the dipole, these are the unspecified variables in ((0, 0, pz), (sx, sy, 4), w)
-	expected_solution = [2, 1, 2, 1]
-
-	result = model.solve(dots)
-	logging.info(result)
-	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,
-	)
-
-	# Do it again with an initial point
-	result = model.solve(dots, initial_pt=[2, 2, 2, 2])
-	logging.info(result)
-	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,
-	)

tests/model/test_fixed_z_discretization.py

@@ -1,20 +0,0 @@
-from pdme.model.fixed_z_plane_model import FixedZPlaneModel, FixedZPlaneDiscretisation
-import numpy
-
-
-def test_fixed_z_plane_model_discretization():
-	model = FixedZPlaneModel(4, -10, 10, -10, 10, 1)
-	discretisation = FixedZPlaneDiscretisation(model, 1, 2, 5, 15)
-	# x: (-10, 0) and (0, 10)
-	# y: (-10, -6, -2, 2, 6, 10)
-	assert discretisation.cell_count == 10
-	assert discretisation.pz_step == 30
-	assert discretisation.x_step == 10
-	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(
-		list(discretisation.all_indices()), list(numpy.ndindex((1, 2, 5)))
-	)

tests/model/test_fixed_z_plane_model.py

@@ -1,46 +0,0 @@
-from pdme.model.fixed_z_plane_model import FixedZPlaneModel
-from pdme.measurement import DotMeasurement
-import logging
-import numpy
-
-
-def test_fixed_z_plane_model_repr():
-	model = FixedZPlaneModel(1, 2, 3, 4, 5, 6)
-	assert repr(model) == "FixedZPlaneModel(1, 2, 3, 4, 5, 6)"
-
-
-def test_fixed_z_plane_model_cost_and_jac_single():
-	model = FixedZPlaneModel(4, -10, 10, -10, 10, 1)
-	measured_v = 0.000191292  # from dipole with p=(0, 0, 2) at (1, 2, 4) with w = 1
-	dot = DotMeasurement(measured_v, (1, 2, 0), 5)
-	pt = [2, 2, 2, 2]
-
-	cost_function = model.costs([dot])
-
-	expected_cost = [0.0000946746]
-	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,
-	)
-
-	jac_function = model.jac([dot])
-
-	expected_jac = [
-		[0.0002859666151836802, -0.0001009293935942401, 0, 0.0001035396365320221]
-	]
-	actual_jac = jac_function(pt)
-
-	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,
-	)

tests/model/test_model_interface.py

@@ -1,35 +1,16 @@
-from pdme.model import Model
+from pdme.model import DipoleModel
-from pdme.measurement import DotMeasurement
 import pytest
 
 
-def test_model_interface_not_implemented_point_length():
+def test_model_interface_not_implemented_one_dipoles():
-	model = Model()
+	model = DipoleModel()
-	with pytest.raises(NotImplementedError):
-		model.point_length()
-
-
-def test_model_interface_not_implemented_point_n():
-	model = Model()
-	with pytest.raises(NotImplementedError):
-		model.n()
-
-
-def test_model_interface_not_implemented_cost():
-	model = Model()
-
-	model.point_length = lambda: 2
 
 	with pytest.raises(NotImplementedError):
-		cost = model.costs([DotMeasurement(0, [1, 2, 3], 4)])
+		model.get_dipoles(5)
-		cost([1, 2])
 
 
-def test_model_interface_not_implemented_jac():
+def test_model_interface_not_implemented_n_dipoles():
-	model = Model()
+	model = DipoleModel()
-
-	model.point_length = lambda: 2
 
 	with pytest.raises(NotImplementedError):
-		jac = model.jac([DotMeasurement(0, [1, 2, 3], 4)])
+		model.get_monte_carlo_dipole_inputs(5, 10)
-		jac([1, 2])

tests/model/test_multiple_dipole_fixed_magnitude_model.py

@@ -0,0 +1,194 @@
+from pdme.model import MultipleDipoleFixedMagnitudeModel
+import numpy
+import logging
+
+
+_logger = logging.getLogger(__name__)
+
+
+def test_repr_multiple_dipole_fixed_mag():
+	model = MultipleDipoleFixedMagnitudeModel(-10, 10, -5, 5, 2, 3, 10, 3.5)
+	assert (
+		repr(model)
+		== "MultipleDipoleFixedMagnitudeModel(-10, 10, -5, 5, 2, 3, 10, 3.5)"
+	), "Repr should be what I want."
+
+
+def test_multiple_dipole_fixed_mag_model_get_dipoles():
+
+	p_fixed = 10
+
+	model = MultipleDipoleFixedMagnitudeModel(-10, 10, -5, 5, 2, 3, p_fixed, 1)
+
+	dipole_arrangement = model.get_dipoles(5, numpy.random.default_rng(1234))
+	dipoles = dipole_arrangement.dipoles
+
+	assert len(dipoles) == 1, "Should have only had one dipole generated."
+	expected_p = numpy.array([-2.20191453, 2.06264523, 9.5339953])
+	expected_s = numpy.array([8.46492468, -2.38307576, 2.31909706])
+	expected_w = 0.5904561648332141
+
+	numpy.testing.assert_allclose(
+		dipoles[0].p, expected_p, err_msg="Random multiple dipole p wasn't as expected"
+	)
+	numpy.testing.assert_allclose(
+		dipoles[0].s, expected_s, err_msg="Random multiple dipole s wasn't as expected"
+	)
+	numpy.testing.assert_allclose(
+		dipoles[0].w, expected_w, err_msg="Random multiple dipole w wasn't as expected"
+	)
+	numpy.testing.assert_allclose(
+		numpy.linalg.norm(dipoles[0].p),
+		p_fixed,
+		err_msg="Should have had the expected dipole moment magnitude.",
+	)
+
+
+def test_multiple_dipole_fixed_mag_model_get_dipoles_multiple():
+
+	p_fixed = 10
+	dipole_count = 5
+
+	model = MultipleDipoleFixedMagnitudeModel(
+		-10, 10, -5, 5, 2, 3, p_fixed, dipole_count
+	)
+
+	dipole_arrangement = model.get_dipoles(20, numpy.random.default_rng(1234))
+	dipoles = dipole_arrangement.dipoles
+
+	assert (
+		len(dipoles) == dipole_count
+	), "Should have had multiple dipole based on count generated."
+
+
+def test_multiple_dipole_fixed_mag_model_get_dipoles_invariant():
+
+	x_min = -10
+	x_max = 10
+	y_min = -5
+	y_max = 5
+	z_min = 2
+	z_max = 3
+	p_fixed = 10
+	max_frequency = 5
+
+	model = MultipleDipoleFixedMagnitudeModel(
+		x_min, x_max, y_min, y_max, z_min, z_max, p_fixed, 1
+	)
+	model.rng = numpy.random.default_rng(1234)
+
+	dipole_arrangement = model.get_dipoles(5)
+	dipoles = dipole_arrangement.dipoles
+
+	assert len(dipoles) == 1, "Should have only had one dipole generated."
+	expected_p = numpy.array([-2.20191453, 2.06264523, 9.5339953])
+	expected_s = numpy.array([8.46492468, -2.38307576, 2.31909706])
+	expected_w = 0.5904561648332141
+
+	numpy.testing.assert_allclose(
+		dipoles[0].p, expected_p, err_msg="Random multiple dipole p wasn't as expected"
+	)
+	numpy.testing.assert_allclose(
+		dipoles[0].s, expected_s, err_msg="Random multiple dipole s wasn't as expected"
+	)
+	numpy.testing.assert_allclose(
+		dipoles[0].w, expected_w, err_msg="Random multiple dipole w wasn't as expected"
+	)
+	for i in range(10):
+		dipole_arrangement = model.get_dipoles(max_frequency)
+		dipoles = dipole_arrangement.dipoles
+
+		assert len(dipoles) == 1, "Should have only had one dipole generated."
+
+		min_s = numpy.array([x_min, y_min, z_min])
+		max_s = numpy.array([x_max, y_max, z_max])
+
+		numpy.testing.assert_equal(
+			numpy.logical_and(min_s < dipoles[0].s, max_s > dipoles[0].s),
+			True,
+			f"Dipole location [{dipoles[0].s}] should have been between min [{min_s}] and max [{max_s}] bounds.",
+		)
+		assert (
+			dipoles[0].w < max_frequency and dipoles[0].w > 0
+		), "Dipole frequency should have been between 0 and max."
+		numpy.testing.assert_allclose(
+			numpy.linalg.norm(dipoles[0].p),
+			p_fixed,
+			err_msg="Should have had the expected dipole moment magnitude.",
+		)
+
+
+def test_multiple_dipole_fixed_mag_model_get_n_dipoles():
+
+	x_min = -10
+	x_max = 10
+	y_min = -5
+	y_max = 5
+	z_min = 2
+	z_max = 3
+	p_fixed = 10
+	max_frequency = 5
+
+	model = MultipleDipoleFixedMagnitudeModel(
+		x_min, x_max, y_min, y_max, z_min, z_max, p_fixed, 1
+	)
+	model.rng = numpy.random.default_rng(1234)
+
+	dipole_array = model.get_monte_carlo_dipole_inputs(1, max_frequency)
+	expected_dipole_array = numpy.array(
+		[
+			[
+				[
+					9.60483896,
+					-1.41627817,
+					-2.3960853,
+					8.46492468,
+					-2.38307576,
+					2.31909706,
+					1.47236493,
+				]
+			]
+		]
+	)
+
+	numpy.testing.assert_allclose(
+		dipole_array,
+		expected_dipole_array,
+		err_msg="Should have had the expected output dipole array.",
+	)
+	numpy.testing.assert_allclose(
+		model.get_monte_carlo_dipole_inputs(
+			1, max_frequency, numpy.random.default_rng(1234)
+		),
+		expected_dipole_array,
+		err_msg="Should have had the expected output dipole array, even with explicitly passed rng.",
+	)
+
+
+def test_multiple_dipole_shape():
+
+	x_min = -10
+	x_max = 10
+	y_min = -5
+	y_max = 5
+	z_min = 2
+	z_max = 3
+	p_fixed = 10
+	max_frequency = 5
+	num_dipoles = 13
+	monte_carlo_n = 11
+
+	model = MultipleDipoleFixedMagnitudeModel(
+		x_min, x_max, y_min, y_max, z_min, z_max, p_fixed, num_dipoles
+	)
+	model.rng = numpy.random.default_rng(1234)
+
+	actual_shape = model.get_monte_carlo_dipole_inputs(
+		monte_carlo_n, max_frequency
+	).shape
+
+	numpy.testing.assert_equal(
+		actual_shape,
+		(monte_carlo_n, num_dipoles, 7),
+		err_msg="shape was wrong for monte carlo outputs",
+	)

tests/model/test_single_dipole_fixed_magnitude_model.py

@@ -0,0 +1,147 @@
+from pdme.model import SingleDipoleFixedMagnitudeModel
+import numpy
+import logging
+
+
+_logger = logging.getLogger(__name__)
+
+
+def test_repr_single_dipole_fixed_mag():
+	model = SingleDipoleFixedMagnitudeModel(-10, 10, -5, 5, 2, 3, 5)
+	assert (
+		repr(model) == "SingleDipoleFixedMagnitudeModel(-10, 10, -5, 5, 2, 3, 5)"
+	), "Repr should be what I want."
+
+
+def test_single_dipole_fixed_mag_model_get_dipoles():
+
+	p_fixed = 10
+
+	model = SingleDipoleFixedMagnitudeModel(-10, 10, -5, 5, 2, 3, p_fixed)
+
+	dipole_arrangement = model.get_dipoles(5, numpy.random.default_rng(1234))
+	dipoles = dipole_arrangement.dipoles
+
+	assert len(dipoles) == 1, "Should have only had one dipole generated."
+	expected_p = numpy.array([-2.20191453, 2.06264523, 9.5339953])
+	expected_s = numpy.array([8.46492468, -2.38307576, 2.31909706])
+	expected_w = 0.5904561648332141
+
+	numpy.testing.assert_allclose(
+		dipoles[0].p, expected_p, err_msg="Random single dipole p wasn't as expected"
+	)
+	numpy.testing.assert_allclose(
+		dipoles[0].s, expected_s, err_msg="Random single dipole s wasn't as expected"
+	)
+	numpy.testing.assert_allclose(
+		dipoles[0].w, expected_w, err_msg="Random single dipole w wasn't as expected"
+	)
+	numpy.testing.assert_allclose(
+		numpy.linalg.norm(dipoles[0].p),
+		p_fixed,
+		err_msg="Should have had the expected dipole moment magnitude.",
+	)
+
+
+def test_single_dipole_fixed_mag_model_get_dipoles_invariant():
+
+	x_min = -10
+	x_max = 10
+	y_min = -5
+	y_max = 5
+	z_min = 2
+	z_max = 3
+	p_fixed = 10
+	max_frequency = 5
+
+	model = SingleDipoleFixedMagnitudeModel(
+		x_min, x_max, y_min, y_max, z_min, z_max, p_fixed
+	)
+	model.rng = numpy.random.default_rng(1234)
+
+	dipole_arrangement = model.get_dipoles(5)
+	dipoles = dipole_arrangement.dipoles
+
+	assert len(dipoles) == 1, "Should have only had one dipole generated."
+	expected_p = numpy.array([-2.20191453, 2.06264523, 9.5339953])
+	expected_s = numpy.array([8.46492468, -2.38307576, 2.31909706])
+	expected_w = 0.5904561648332141
+
+	numpy.testing.assert_allclose(
+		dipoles[0].p, expected_p, err_msg="Random single dipole p wasn't as expected"
+	)
+	numpy.testing.assert_allclose(
+		dipoles[0].s, expected_s, err_msg="Random single dipole s wasn't as expected"
+	)
+	numpy.testing.assert_allclose(
+		dipoles[0].w, expected_w, err_msg="Random single dipole w wasn't as expected"
+	)
+	for i in range(10):
+		dipole_arrangement = model.get_dipoles(max_frequency)
+		dipoles = dipole_arrangement.dipoles
+
+		assert len(dipoles) == 1, "Should have only had one dipole generated."
+
+		min_s = numpy.array([x_min, y_min, z_min])
+		max_s = numpy.array([x_max, y_max, z_max])
+
+		numpy.testing.assert_equal(
+			numpy.logical_and(min_s < dipoles[0].s, max_s > dipoles[0].s),
+			True,
+			f"Dipole location [{dipoles[0].s}] should have been between min [{min_s}] and max [{max_s}] bounds.",
+		)
+		assert (
+			dipoles[0].w < max_frequency and dipoles[0].w > 0
+		), "Dipole frequency should have been between 0 and max."
+		numpy.testing.assert_allclose(
+			numpy.linalg.norm(dipoles[0].p),
+			p_fixed,
+			err_msg="Should have had the expected dipole moment magnitude.",
+		)
+
+
+def test_single_dipole_fixed_mag_model_get_n_dipoles():
+
+	x_min = -10
+	x_max = 10
+	y_min = -5
+	y_max = 5
+	z_min = 2
+	z_max = 3
+	p_fixed = 10
+	max_frequency = 5
+
+	model = SingleDipoleFixedMagnitudeModel(
+		x_min, x_max, y_min, y_max, z_min, z_max, p_fixed
+	)
+	model.rng = numpy.random.default_rng(1234)
+
+	dipole_array = model.get_monte_carlo_dipole_inputs(1, max_frequency)
+	expected_dipole_array = numpy.array(
+		[
+			[
+				[
+					9.60483896,
+					-1.41627817,
+					-2.3960853,
+					8.46492468,
+					-2.38307576,
+					2.31909706,
+					1.47236493,
+				]
+			]
+		]
+	)
+
+	numpy.testing.assert_allclose(
+		dipole_array,
+		expected_dipole_array,
+		err_msg="Should have had the expected output dipole array.",
+	)
+	numpy.testing.assert_allclose(
+		model.get_monte_carlo_dipole_inputs(
+			1, max_frequency, numpy.random.default_rng(1234)
+		),
+		expected_dipole_array,
+		err_msg="Should have had the expected output dipole array, even with explicitly passed rng.",
+	)

tests/model/test_unrestricted_basic_solve.py

@@ -1,46 +0,0 @@
-from pdme.model import UnrestrictedModel
-from pdme.measurement import OscillatingDipole, OscillatingDipoleArrangement
-import logging
-import numpy
-import itertools
-
-
-def test_unrestricted_model_solve_basic():
-	# Initialise our dipole arrangement and create dot measurements along a square.
-	dipoles = OscillatingDipoleArrangement(
-		[OscillatingDipole((0.2, 0, 2), (1, 2, 4), 1)]
-	)
-	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)
-
-	model = UnrestrictedModel(1, -1, 1, -1, 1, -1, 5, 1)
-
-	# from the dipole, these are the unspecified variables in ((0, 0, 2), (1, 2, 4), 1)
-	expected_solution = [0.2, 0, 2, 1, 2, 4, 1]
-
-	result = model.solve(dots)
-	logging.info(result)
-	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,
-	)

tests/model/test_unrestricted_discretization.py

@@ -1,23 +0,0 @@
-from pdme.model.unrestricted_model import UnrestrictedModel, UnrestrictedDiscretisation
-import numpy
-
-
-def test_unrestricted_model_discretization():
-	model = UnrestrictedModel(-10, 10, -10, 10, -10, 10, 15, 1)
-	discretisation = UnrestrictedDiscretisation(model, 1, 1, 2, 2, 5, 1)
-	# x: (-10, 0) and (0, 10)
-	# y: (-10, -6, -2, 2, 6, 10)
-	assert discretisation.cell_count == 10
-	assert discretisation.px_step == 30
-	assert discretisation.py_step == 30
-	assert discretisation.pz_step == 15
-	assert discretisation.x_step == 10
-	assert discretisation.y_step == 4
-	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(
-		list(discretisation.all_indices()), list(numpy.ndindex((1, 1, 2, 2, 5, 1)))
-	)

tests/model/test_unrestricted_model.py

@@ -1,54 +0,0 @@
-from pdme.model import UnrestrictedModel
-from pdme.measurement import DotMeasurement
-import logging
-import numpy
-
-
-def test_unrestricted_plane_model_repr():
-	model = UnrestrictedModel(1, 2, 3, 4, 5, 6, 7, 6)
-	assert repr(model) == "UnrestrictedModel(1, 2, 3, 4, 5, 6, 7, 6)"
-
-
-def test_unrestricted_model_cost_and_jac_single():
-	model = UnrestrictedModel(1, -1, 1, -1, 1, -1, 1, 1)
-	measured_v = 0.000191292  # from dipole with p=(0, 0, 2) at (1, 2, 4) with w = 1
-	dot = DotMeasurement(measured_v, (1, 2, 0), 5)
-	pt = numpy.array([0, 0, 2, 2, 2, 4, 2])
-
-	cost_function = model.costs([dot])
-
-	expected_cost = [0.0000946746]
-	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,
-	)
-
-	jac_function = model.jac([dot])
-
-	expected_jac = [
-		[
-			0.00007149165379592005,
-			0,
-			0.0002859666151836802,
-			-0.0001009293935942401,
-			0,
-			-0.0002607342667851202,
-			0.0001035396365320221,
-		]
-	]
-	actual_jac = jac_function(pt)
-
-	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,
-	)

tests/util/test_fast_nonlocal_spectrum_pairs.py

@@ -0,0 +1,55 @@
+import numpy
+import pdme.util.fast_nonlocal_spectrum
+import logging
+import pytest
+
+
+def test_fast_nonlocal_calc_multidipole():
+	d1 = [1, 2, 3, 4, 5, 6, 7]
+	d2 = [1, 2, 3, 4, 5, 6, 8]
+	d3 = [2, 5, 3, 4, -5, -6, 2]
+	d4 = [-3, 2, 1, 4, 5, 6, 10]
+
+	dipoleses = numpy.array([[d1, d2], [d3, d4]])
+
+	dot_pairs = numpy.array(
+		[[[-1, -2, -3, 11], [-1, 2, 5, 11]], [[-1, -2, -3, 6], [2, 4, 6, 6]]]
+	)
+	# expected_ij is for pair i, dipole j
+	expected_11 = 0.000021124454334947546213
+	expected_12 = 0.000022184755131682365135
+	expected_13 = 0.0000053860643617855849275
+	expected_14 = -0.0000023069501696755220764
+	expected_21 = 0.00022356021100884617796
+	expected_22 = 0.00021717277640859343002
+	expected_23 = 0.000017558321044891869169
+	expected_24 = -0.000034714318479634499683
+
+	expected = numpy.array(
+		[
+			[expected_11 + expected_12, expected_21 + expected_22],
+			[expected_13 + expected_14, expected_23 + 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.
+	pdme.util.fast_nonlocal_spectrum._logger.setLevel(logging.DEBUG)
+
+	numpy.testing.assert_allclose(
+		pdme.util.fast_nonlocal_spectrum.fast_s_nonlocal_dipoleses(
+			dot_pairs, dipoleses
+		),
+		expected,
+		err_msg="nonlocal voltages at dot aren't as expected for dipoleses.",
+	)
+
+
+def test_fast_nonlocal_frequency_check_multidipole():
+	d1 = [1, 2, 3, 4, 5, 6, 7]
+
+	dipoles = numpy.array([[d1]])
+
+	dot_pairs = numpy.array([[[-1, -2, -3, 11], [-1, 2, 5, 10]]])
+
+	with pytest.raises(ValueError):
+		pdme.util.fast_nonlocal_spectrum.fast_s_nonlocal_dipoleses(dot_pairs, dipoles)

tests/util/test_fast_v_calc.py

@@ -24,6 +24,78 @@ def test_fast_v_calc():
 	)
 
 
+def test_fast_v_calc_multidipoles():
+	d1 = [1, 2, 3, 4, 5, 6, 7]
+	d2 = [2, 5, 3, 4, -5, -6, 2]
+
+	dipoles = numpy.array([[d1, d2]])
+
+	dot_inputs = numpy.array([[-1, -1, -1, 11], [2, 3, 1, 5.5]])
+	# expected_ij is for dot i, dipole j
+	expected_11 = 0.00001421963287022476
+	expected_12 = 0.00001107180225755457
+	expected_21 = 0.000345021108583681380388722
+	expected_22 = 0.0000377061050587914705139781
+
+	expected = numpy.array([[expected_11 + expected_12, expected_21 + expected_22]])
+
+	numpy.testing.assert_allclose(
+		pdme.util.fast_v_calc.fast_vs_for_dipoleses(dot_inputs, dipoles),
+		expected,
+		err_msg="Voltages at dot aren't as expected for multidipole calc.",
+	)
+
+
+def test_fast_v_calc_big_multidipole():
+
+	dipoles = numpy.array(
+		[
+			[
+				[1, 1, 5, 6, 3, 1, 1],
+				[5, 3, 2, 13, 1, 1, 2],
+				[-5, -5, -3, -1, -3, 8, 3],
+			],
+			[
+				[-3, -1, -2, -2, -6, 3, 4],
+				[8, 0, 2, 0, 1, 5, 5],
+				[1, 4, -4, -1, -3, -5, 6],
+			],
+		]
+	)
+
+	dot_inputs = numpy.array(
+		[
+			[1, 1, 0, 1],
+			[2, 5, 6, 2],
+			[3, 1, 3, 3],
+			[0.5, 0.5, 0.5, 4],
+		]
+	)
+
+	expected = numpy.array(
+		[
+			[
+				0.0010151687742365581690202135,
+				0.00077627527320628609782627266,
+				0.00043313471258511003340648713,
+				0.000077184305988088453637005111,
+			],
+			[
+				0.000041099091967966890657097060,
+				0.0019377687238977568792327845,
+				0.0085903193415282984161225029,
+				0.00014557676715208209310911838,
+			],
+		]
+	)
+
+	numpy.testing.assert_allclose(
+		pdme.util.fast_v_calc.fast_vs_for_dipoleses(dot_inputs, dipoles),
+		expected,
+		err_msg="Voltages at dot aren't as expected for multidipole calc.",
+	)
+
+
 def test_between():
 	low = numpy.array([1, 2, 3])
 	high = numpy.array([6, 7, 8])