doc:

Reviewed-on: #30

.editorconfig

@@ -3,5 +3,5 @@ root = true
 # Unix-style newlines with a newline ending every file
 [*]
 end_of_line = lf
-insert_final_newline = true
+# insert_final_newline = true
 indent_style = tab

.flake8

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

.gitignore

@@ -114,6 +114,10 @@ ENV/
 env.bak/
 venv.bak/
 
+#direnv
+.envrc
+.direnv
+
 # Spyder project settings
 .spyderproject
 .spyproject
@@ -136,4 +140,7 @@ dmypy.json
 .pytype/
 
 # Cython debug symbols
-cython_debug/
+cython_debug/
+
+.csv
+.idea

.versionrc

@@ -0,0 +1,10 @@
+{
+	"bumpFiles": [
+		{
+			"filename": "pyproject.toml",
+			"updater": "scripts/standard-version/pyproject-updater.js"
+		}
+	],
+	"sign": true,
+	"tag-prefix": ""
+}

CHANGELOG.md

@@ -0,0 +1,249 @@
+# Changelog
+
+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.
+
+## [1.5.0](https://gitea.deepak.science:2222/physics/pdme/compare/1.4.0...1.5.0) (2024-05-17)
+
+
+### Features
+
+* adds mcmc chain that returns number of repeats ([6193ecb](https://gitea.deepak.science:2222/physics/pdme/commit/6193ecb9c9f7a21d24e860987a7107549a4b2fa7))
+
+## [1.4.0](https://gitea.deepak.science:2222/physics/pdme/compare/1.3.0...1.4.0) (2024-05-17)
+
+
+### Features
+
+* adds relative squared diff calc utility method ([9b1538b](https://gitea.deepak.science:2222/physics/pdme/commit/9b1538b3c63bfaf2a779bb109cd160a8d7887195))
+
+## [1.3.0](https://gitea.deepak.science:2222/physics/pdme/compare/1.2.0...1.3.0) (2024-05-17)
+
+
+### Features
+
+* adds utility function for sorting samples by frequency for subspace simulation ([e5fc120](https://gitea.deepak.science:2222/physics/pdme/commit/e5fc1207a8b7d5b67208ad825907baa442eec648))
+
+## [1.2.0](https://gitea.deepak.science:2222/physics/pdme/compare/1.1.0...1.2.0) (2024-05-03)
+
+
+### Features
+
+* adds pdme fast calc for e field xs ([9e6d1df](https://gitea.deepak.science:2222/physics/pdme/commit/9e6d1df559e58998851a1c2bf24fcc46d8c1b148))
+
+## [1.1.0](https://gitea.deepak.science:2222/physics/pdme/compare/1.0.0...1.1.0) (2024-05-02)
+
+
+### Features
+
+* adds both electric potential and electric field x sources, makes some fast util tests use the slower explicit versions as double check ([e9e3416](https://gitea.deepak.science:2222/physics/pdme/commit/e9e34162a3b84faad5c18ddeda327c2f7f5ac5aa))
+
+## [1.0.0](https://gitea.deepak.science:2222/physics/pdme/compare/0.9.3...1.0.0) (2024-04-29)
+
+
+### Bug Fixes
+
+* fixes the broken implementation of the tarucha frequency calculation ([631ba13](https://gitea.deepak.science:2222/physics/pdme/commit/631ba13c791c71ad8922d39a13b780a40eac2391))
+
+### [0.9.3](https://gitea.deepak.science:2222/physics/pdme/compare/0.9.2...0.9.3) (2024-02-26)
+
+
+### Features
+
+* adds util func for calculating arg using sign instead of complex arithmetic ([3ebe2bb](https://gitea.deepak.science:2222/physics/pdme/commit/3ebe2bb82430d677680383c42a1c269df83d99cd))
+
+
+### Bug Fixes
+
+* fixes stupid cost shape issue ([ed9dd2c](https://gitea.deepak.science:2222/physics/pdme/commit/ed9dd2c94f88a08c36f581f05b26a87a6b780d5b))
+
+### [0.9.2](https://gitea.deepak.science:2222/physics/pdme/compare/0.9.1...0.9.2) (2023-07-24)
+
+
+### Bug Fixes
+
+* update tests but for git also don't wrap costs ([50f98ed](https://gitea.deepak.science:2222/physics/pdme/commit/50f98ed89b2a05cd47c41958036dd50bc872e07c))
+
+### [0.9.1](https://gitea.deepak.science:2222/physics/pdme/compare/0.9.0...0.9.1) (2023-07-24)
+
+
+### Bug Fixes
+
+* fixes some of the shape mangling of our mcmc code ([e01d0e1](https://gitea.deepak.science:2222/physics/pdme/commit/e01d0e14a9bcd6d7e8fe9449ce562dbf1b8fd25c))
+
+## [0.9.0](https://gitea.deepak.science:2222/physics/pdme/compare/0.8.9...0.9.0) (2023-07-24)
+
+
+### ⚠ BREAKING CHANGES
+
+* separates threshold cost and the seed_cost in mcmc
+
+### Features
+
+* separates threshold cost and the seed_cost in mcmc ([ca710e3](https://gitea.deepak.science:2222/physics/pdme/commit/ca710e359fd0cfbb620a3574a2fa4fab1be2b52a))
+
+### [0.8.9](https://gitea.deepak.science:2222/physics/pdme/compare/0.8.8...0.8.9) (2023-07-23)
+
+
+### Features
+
+* adds a bunch of mcmc generation code for log spaced models, yay ([f280448](https://gitea.deepak.science:2222/physics/pdme/commit/f280448cfe2fcf5bdc5ac2317ee52b27523bb49d))
+* adds utility functions for dealing with markov chain monte carlo ([feb0a5f](https://gitea.deepak.science:2222/physics/pdme/commit/feb0a5f6453dcb5e71a07c7749cd579dab15171c))
+
+### [0.8.8](https://gitea.deepak.science:2222/physics/pdme/compare/0.8.7...0.8.8) (2023-04-09)
+
+
+### Features
+
+* adds fast calc that allows for variable temp ([36454d5](https://gitea.deepak.science:2222/physics/pdme/commit/36454d5044d93b6b178e016b84dd59a5ebaf15e2))
+
+### [0.8.7](https://gitea.deepak.science:2222/physics/pdme/compare/0.8.6...0.8.7) (2022-09-17)
+
+
+### Features
+
+* adds xy model for convenience to pdme ([e894c89](https://gitea.deepak.science:2222/physics/pdme/commit/e894c897029c05a1d4754e7930ae9ba2be7a1cfd))
+* moves xy model up to model package ([2a1ae3b](https://gitea.deepak.science:2222/physics/pdme/commit/2a1ae3b1a7f7e10469b7fd2930fee0b338f0c03f))
+
+
+### Bug Fixes
+
+* Correctly generates monte carlo version of xy model dipoles ([5acf0ac](https://gitea.deepak.science:2222/physics/pdme/commit/5acf0ac347382705674bb596440d27cba3730bac))
+
+### [0.8.6](https://gitea.deepak.science:2222/physics/pdme/compare/0.8.5...0.8.6) (2022-06-13)
+
+
+### Features
+
+* makes library build system poetry-core ([ed0c6e2](https://gitea.deepak.science:2222/physics/pdme/commit/ed0c6e2858f45fec9b9a673d9b5bc98605e73508))
+
+
+### Bug Fixes
+
+* pyproject build system now core as well ([6277e84](https://gitea.deepak.science:2222/physics/pdme/commit/6277e843d5a7df8eba1878db490d2fae4052af57))
+
+### [0.8.5](https://gitea.deepak.science:2222/physics/pdme/compare/0.8.4...0.8.5) (2022-06-04)
+
+
+### Features
+
+* adds fixedorientation model ([2cdf46a](https://gitea.deepak.science:2222/physics/pdme/commit/2cdf46afa1492f86b4c403e4c5013cefc89b21d6))
+
+### [0.8.4](https://gitea.deepak.science:2222/physics/pdme/compare/0.8.3...0.8.4) (2022-05-26)
+
+
+### Bug Fixes
+
+* probability of occupancy fixed ([01f905a](https://gitea.deepak.science:2222/physics/pdme/commit/01f905a237a2423f5637ee6a0f43e0937c55d2ea))
+
+### [0.8.3](https://gitea.deepak.science:2222/physics/pdme/compare/0.8.2...0.8.3) (2022-05-22)
+
+
+### Features
+
+* Adds log spaced in frequency space model ([7a9fa1b](https://gitea.deepak.science:2222/physics/pdme/commit/7a9fa1ba04a12586ef09c89e35e706817012faab))
+
+### [0.8.2](https://gitea.deepak.science:2222/physics/pdme/compare/0.8.1...0.8.2) (2022-05-07)
+
+
+### Features
+
+* Adds Random count dipole model with binomial type distribution for dipole number ([27cbb36](https://gitea.deepak.science:2222/physics/pdme/commit/27cbb364a6d9625abca6145dafedf0df5743f816))
+
+### [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)
+
+
+### Features
+
+* adds methods for converting pdme objects into flatter numpy arrays ([5c6c4c7](https://gitea.deepak.science:2222/physics/pdme/commit/5c6c4c79d1b55d8d6bf70e2e54226d282d831f66))
+* Allows you to pass in rng to generate monte carlo dipoles ([ab244ed](https://gitea.deepak.science:2222/physics/pdme/commit/ab244ed91dede1e3c29fe94ea96089f4c15aa1fb))
+
+### [0.6.1](https://gitea.deepak.science:2222/physics/pdme/compare/0.6.0...0.6.1) (2022-03-28)
+
+
+### Bug Fixes
+
+* Swaps high and lows for ranges if needed to make negatives behave nicer ([8a60967](https://gitea.deepak.science:2222/physics/pdme/commit/8a60967ddfa18e49460b05edc7ec575f06db868f))
+
+## [0.6.0](https://gitea.deepak.science:2222/physics/pdme/compare/0.5.4...0.6.0) (2022-03-27)
+
+
+### ⚠ BREAKING CHANGES
+
+* adds pair inputs to array function
+
+### Features
+
+* adds calc for pairs and pair ranges ([e72489f](https://gitea.deepak.science:2222/physics/pdme/commit/e72489f0cb6586756180c2021649f4eb019c77fc))
+* adds fast nonlocal noise Sij calculator ([5fbff2a](https://gitea.deepak.science:2222/physics/pdme/commit/5fbff2a5c0f128ffc5b749788c8bdcd6f88de234))
+* adds methods for better creation of input lists, including pairs ([852836b](https://gitea.deepak.science:2222/physics/pdme/commit/852836b924c8ffefd6b28a8ea1bc1ca47b77756e))
+* adds pair inputs to array function ([a7508b8](https://gitea.deepak.science:2222/physics/pdme/commit/a7508b8906923a579dfa05edc22963d2f0102be9))
+
+### [0.5.4](https://gitea.deepak.science:2222/physics/pdme/compare/0.5.3...0.5.4) (2022-03-06)
+
+
+### Bug Fixes
+
+* lets you get model from discretisation for fixed magnitude models in type hintable way ([548dcfe](https://gitea.deepak.science:2222/physics/pdme/commit/548dcfebfc7cd6354ad708b3ab0e664881c58de4))
+
+### [0.5.3](https://gitea.deepak.science:2222/physics/pdme/compare/0.5.2...0.5.3) (2022-03-06)
+
+
+### Features
+
+* Adds between utility function ([a782b49](https://gitea.deepak.science:2222/physics/pdme/commit/a782b4957f2aebfb0bf9a3c7d4166277a5450b76))
+* Adds fast calc expressions ([0855911](https://gitea.deepak.science:2222/physics/pdme/commit/08559110be197b7f4a685d5e44e248f914853a3c))
+* Adds numpy make more dipoles at once and converters ([c70a971](https://gitea.deepak.science:2222/physics/pdme/commit/c70a971ce09bda236a84cfcff309d9a9cf2092d9))
+* Calculates range measurements. ([2d773df](https://gitea.deepak.science:2222/physics/pdme/commit/2d773df7259378ba9285776407b2cb5ae380d809))
+
+### [0.5.2](https://gitea.deepak.science:2222/physics/pdme/compare/0.5.1...0.5.2) (2022-03-06)
+
+### [0.5.1](https://gitea.deepak.science:2222/physics/pdme/compare/0.5.0...0.5.1) (2022-03-05)
+
+
+### Bug Fixes
+
+* Fixes typo ([31762de](https://gitea.deepak.science:2222/physics/pdme/commit/31762de69785b9c25704605b132c52036ad6dad3))

Jenkinsfile

@@ -4,7 +4,7 @@ pipeline {
 		label 'pdme'  // all your pods will be named with this prefix, followed by a unique id
 		idleMinutes 5  // how long the pod will live after no jobs have run on it
 		yamlFile 'jenkins/ci-agent-pod.yaml'  // path to the pod definition relative to the root of our project
-		defaultContainer 'python'  // define a default container if more than a few stages use it, will default to jnlp container
+		defaultContainer 'poetry'  // define a default container if more than a few stages use it, will default to jnlp container
 	  }
 	}
 
@@ -12,36 +12,30 @@ pipeline {
 		parallelsAlwaysFailFast()
 	}
 
-	environment {
-		POETRY_HOME="/opt/poetry"
-		POETRY_VERSION="1.1.12"
-	}
-
 	stages {
 		stage('Build') {
 			steps {
 				echo 'Building...'
 				sh 'python --version'
-				sh 'curl -sSL https://raw.githubusercontent.com/python-poetry/poetry/master/get-poetry.py | python'
-				sh '${POETRY_HOME}/bin/poetry --version'
-				sh '${POETRY_HOME}/bin/poetry install'
+				sh 'poetry --version'
+				sh 'poetry install'
 			}
 		}
 		stage('Test') {
 			parallel{
 				stage('pytest') {
 					steps {
-						sh '${POETRY_HOME}/bin/poetry run pytest'
+						sh 'poetry run pytest'
 					}
 				}
 				stage('lint') {
 					steps {
-						sh '${POETRY_HOME}/bin/poetry run flake8'
+						sh 'poetry run flake8 pdme tests'
 					}
 				}
 				stage('mypy') {
 					steps {
-						sh '${POETRY_HOME}/bin/poetry run mypy pdme'
+						sh 'poetry run mypy pdme'
 					}
 				}
 			}
@@ -57,7 +51,7 @@ pipeline {
 			}
 			steps {
 				echo 'Deploying...'
-				sh '${POETRY_HOME}/bin/poetry publish -u ${PYPI_USR} -p ${PYPI_PSW} --build'
+				sh 'poetry publish -u ${PYPI_USR} -p ${PYPI_PSW} --build'
 			}
 		}
 

README.md

@@ -1 +1,23 @@
 # pdme - the python dipole model evaluator
+
+[![Conventional Commits](https://img.shields.io/badge/Conventional%20Commits-1.0.0-green.svg?style=flat-square)](https://conventionalcommits.org)
+[![PyPI](https://img.shields.io/pypi/v/pdme?style=flat-square)](https://pypi.org/project/pdme/)
+![Maintenance](https://img.shields.io/maintenance/yes/2024?style=flat-square)
+
+This library includes a variety of utility functions for calculating charge noise in qubits from TLSs for the purposes of Bayesian analysis.
+This was a major part of my PhD thesis on characterizing charge noise.
+New code here should be mostly low-level numpy calculations or common types that can be shared elsewhere.
+
+## Related libraries
+
+Some other libraries I wrote related to this project include:
+- [deepdog](https://gitea.deepak.science/physics/deepdog): performs Monte Carlo simulations with two underlying algorithm options for direct MC and subset simulation
+- [kalpaa](https://gitea.deepak.science/physics/kalpa): high level wrapper for pdme, deepdog and tantri, implementing retry logic, data transformation and summarising. the place to start to actually run analyses
+- [tantri](https://gitea.deepak.science/physics/tantri): generates mock telegraph noise and utility functions for binning and transforming measured or simulated power spectral densities (PSDs)
+- [dreader](https://gitea.deepak.science/physics/dreader): utility library for parsing output files
+
+## Getting started
+
+`poetry install` to start locally
+
+Commit using [Conventional Commits](https://www.conventionalcommits.org/en/v1.0.0/), and when commits are on master, release with `doo release`.

do.sh

@@ -1,29 +0,0 @@
-#!/usr/bin/env bash
-# Do - The Simplest Build Tool on Earth.
-# Documentation and examples see https://github.com/8gears/do
-
-set -Eeuo pipefail # -e "Automatic exit from bash shell script on error"  -u "Treat unset variables and parameters as errors"
-
-build() {
-   echo "I am ${FUNCNAME[0]}ing"
-   poetry build
-}
-
-test() {
-   echo "I am ${FUNCNAME[0]}ing"
-   poetry run flake8 pdme tests
-   poetry run mypy pdme
-   poetry run pytest
-}
-
-htmlcov() {
-	poetry run pytest --cov-report=html
-}
-
-all() {
-   build && test
-}
-
-"$@" # <- execute the task
-
-[ "$#" -gt 0 ] || printf "Usage:\n\t./do.sh %s\n" "($(compgen -A function | grep '^[^_]' | paste -sd '|' -))"

flake.lock

@@ -0,0 +1,95 @@
+{
+  "nodes": {
+    "flake-utils": {
+      "locked": {
+        "lastModified": 1648297722,
+        "narHash": "sha256-W+qlPsiZd8F3XkzXOzAoR+mpFqzm3ekQkJNa+PIh1BQ=",
+        "owner": "numtide",
+        "repo": "flake-utils",
+        "rev": "0f8662f1319ad6abf89b3380dd2722369fc51ade",
+        "type": "github"
+      },
+      "original": {
+        "owner": "numtide",
+        "repo": "flake-utils",
+        "rev": "0f8662f1319ad6abf89b3380dd2722369fc51ade",
+        "type": "github"
+      }
+    },
+    "flake-utils_2": {
+      "locked": {
+        "lastModified": 1653893745,
+        "narHash": "sha256-0jntwV3Z8//YwuOjzhV2sgJJPt+HY6KhU7VZUL0fKZQ=",
+        "owner": "numtide",
+        "repo": "flake-utils",
+        "rev": "1ed9fb1935d260de5fe1c2f7ee0ebaae17ed2fa1",
+        "type": "github"
+      },
+      "original": {
+        "owner": "numtide",
+        "repo": "flake-utils",
+        "type": "github"
+      }
+    },
+    "nixpkgs": {
+      "locked": {
+        "lastModified": 1648854265,
+        "narHash": "sha256-e/RlfodBOMr2SH9diDPYMraTWvhOWSSsXDQikHFdUvM=",
+        "owner": "NixOS",
+        "repo": "nixpkgs",
+        "rev": "e194871435cad8ffb1d64b64fb7df3b2b8a10088",
+        "type": "github"
+      },
+      "original": {
+        "owner": "NixOS",
+        "repo": "nixpkgs",
+        "rev": "e194871435cad8ffb1d64b64fb7df3b2b8a10088",
+        "type": "github"
+      }
+    },
+    "nixpkgs_2": {
+      "locked": {
+        "lastModified": 1655043425,
+        "narHash": "sha256-A+oT+aQGhW5lXy8H0cqBLsYtgcnT5glmGOXWQDcGw6I=",
+        "owner": "NixOS",
+        "repo": "nixpkgs",
+        "rev": "914ef51ffa88d9b386c71bdc88bffc5273c08ada",
+        "type": "github"
+      },
+      "original": {
+        "owner": "NixOS",
+        "repo": "nixpkgs",
+        "type": "github"
+      }
+    },
+    "poetry2nix": {
+      "inputs": {
+        "flake-utils": "flake-utils_2",
+        "nixpkgs": "nixpkgs_2"
+      },
+      "locked": {
+        "lastModified": 1654921554,
+        "narHash": "sha256-hkfMdQAHSwLWlg0sBVvgrQdIiBP45U1/ktmFpY4g2Mo=",
+        "owner": "nix-community",
+        "repo": "poetry2nix",
+        "rev": "7b71679fa7df00e1678fc3f1d1d4f5f372341b63",
+        "type": "github"
+      },
+      "original": {
+        "owner": "nix-community",
+        "repo": "poetry2nix",
+        "rev": "7b71679fa7df00e1678fc3f1d1d4f5f372341b63",
+        "type": "github"
+      }
+    },
+    "root": {
+      "inputs": {
+        "flake-utils": "flake-utils",
+        "nixpkgs": "nixpkgs",
+        "poetry2nix": "poetry2nix"
+      }
+    }
+  },
+  "root": "root",
+  "version": 7
+}

flake.nix

@@ -0,0 +1,55 @@
+{
+  description = "Application packaged using poetry2nix";
+
+  inputs.flake-utils.url = "github:numtide/flake-utils?rev=0f8662f1319ad6abf89b3380dd2722369fc51ade";
+  inputs.nixpkgs.url = "github:NixOS/nixpkgs?rev=e194871435cad8ffb1d64b64fb7df3b2b8a10088";
+  inputs.poetry2nix.url = "github:nix-community/poetry2nix?rev=7b71679fa7df00e1678fc3f1d1d4f5f372341b63";
+
+  outputs = { self, nixpkgs, flake-utils, poetry2nix }:
+    {
+      # Nixpkgs overlay providing the application
+      overlay = nixpkgs.lib.composeManyExtensions [
+        poetry2nix.overlay
+        (final: prev: {
+          # The application
+          pdme = prev.poetry2nix.mkPoetryApplication {
+            overrides = [
+              prev.poetry2nix.defaultPoetryOverrides
+            ];
+            projectDir = ./.;
+          };
+          pdmeEnv = prev.poetry2nix.mkPoetryEnv {
+            overrides = [
+              prev.poetry2nix.defaultPoetryOverrides
+            ];
+            projectDir = ./.;
+          };
+        })
+      ];
+    } // (flake-utils.lib.eachDefaultSystem (system:
+      let
+        pkgs = import nixpkgs {
+          inherit system;
+          overlays = [ self.overlay ];
+        };
+      in
+      {
+        apps = {
+          pdme = pkgs.pdme;
+        };
+
+        defaultApp = pkgs.pdme;
+        devShell = pkgs.mkShell {
+          buildInputs = [
+            pkgs.poetry
+            pkgs.pdmeEnv
+            pkgs.pdme
+          ];
+          shellHook = ''
+            export DO_NIX_CUSTOM=1
+          '';
+          packages = [ pkgs.nodejs-16_x pkgs.gnupg ];
+        };
+
+      }));
+}

jenkins/ci-agent-pod.yaml

@@ -1,9 +1,11 @@
 apiVersion: v1
 kind: Pod
 spec:
+  imagePullSecrets:
+    - name: regcreds
   containers:  # list of containers that you want present for your build, you can define a default container in the Jenkinsfile
-    - name: python
-      image: python:3.8
+    - name: poetry
+      image: ghcr.io/dmallubhotla/poetry-image:1
       command: ["tail", "-f", "/dev/null"]  # this or any command that is bascially a noop is required, this is so that you don't overwrite the entrypoint of the base container
       imagePullPolicy: Always # use cache or pull image for agent
       resources:  # limits the resources your build contaienr

justfile

@@ -0,0 +1,64 @@
+
+# execute default build
+default: build
+
+# builds the python module using poetry
+build:
+	echo "building..."
+	poetry build
+
+# print a message displaying whether nix is being used
+checknix:
+	#!/usr/bin/env bash
+	set -euxo pipefail
+	if [[ "${DO_NIX_CUSTOM:=0}" -eq 1 ]]; then
+		echo "In an interactive nix env."
+	else
+		echo "Using poetry as runner, no nix detected."
+	fi
+
+# run all tests
+test: fmt
+	#!/usr/bin/env bash
+	set -euxo pipefail
+	
+	if [[ "${DO_NIX_CUSTOM:=0}" -eq 1 ]]; then
+		echo "testing, using nix..."
+		flake8 pdme tests
+		mypy pdme
+		pytest
+	else
+		echo "testing..."
+		poetry run flake8 pdme tests
+		poetry run mypy pdme
+		poetry run pytest
+	fi
+
+# update all test snapshots, use if snapshots are out of date
+update-snapshots:
+		#!/usr/bin/env bash
+		set -euxo pipefail
+		if [[ "${DO_NIX_CUSTOM:=0}" -eq 1 ]]; then
+				pytest --snapshot-update
+		else
+				poetry run pytest --snapshot-update
+		fi
+
+# format code
+fmt:
+	#!/usr/bin/env bash
+	set -euxo pipefail
+	if [[ "${DO_NIX_CUSTOM:=0}" -eq 1 ]]; then
+	      black .
+	else
+		poetry run black .
+	fi
+	find pdme -type f -name "*.py" -exec sed -i -e 's/    /\t/g' {} \;
+	find tests -type f -name "*.py" -exec sed -i -e 's/    /\t/g' {} \;
+
+# release the app, checking that our working tree is clean and ready for release
+release:
+	./scripts/release.sh
+
+htmlcov:
+	poetry run pytest --cov-report=html

pdme/calculations/__init__.py

@@ -0,0 +1,46 @@
+"""
+This module is a canonical source of the accurate expressions we want to use for calculating our noise.
+No reference to class or anything, just a straight set of math functions.
+"""
+import numpy
+
+
+def telegraph_beta(f: float, w: float) -> float:
+	"""
+	This function represents the frequency component of analytic telegraph noise.
+
+	We're assuming we care about the one-sided PSD where we are ignoring negative frequencies.
+	This matches with experimental data from say Connors et al., and I think is better than keeping with one-sided.
+	Note that this means that it will only be comparable then with time series data assuming one-sided!
+
+	Don't bikeshed yet, if we care about two-sided things for any reason down the line divide this by two or just change it then.
+	"""
+	return 2 * w / ((numpy.pi * f) ** 2 + w**2)
+
+
+def electric_potential(p: numpy.ndarray, s: numpy.ndarray, r: numpy.ndarray) -> float:
+	"""
+	Gives the electric potential of a defect with dipole moment p, located at position s,
+	as measured from position r.
+
+	p, s, r, are numpy arrays of length 3
+	"""
+	diff = r - s
+	return (p.dot(diff) / (numpy.linalg.norm(diff) ** 3)).item()
+
+
+def electric_field(
+	p: numpy.ndarray, s: numpy.ndarray, r: numpy.ndarray
+) -> numpy.ndarray:
+	"""
+	Gives the electric field of a defect with dipole moment p, located at position s,
+	as measured from position r.
+
+	p, s, r, are numpy arrays of length 3
+
+	Returns an array of length 3, ideally.
+	"""
+	diff = r - s
+	norm_diff = numpy.linalg.norm(diff)
+
+	return ((3 * (p.dot(diff) * diff) / (norm_diff**2)) - p) / (norm_diff**3)

pdme/inputs/__init__.py

@@ -0,0 +1,7 @@
+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"]

pdme/inputs/dot_inputs.py

@@ -0,0 +1,21 @@
+import numpy
+import numpy.typing
+import itertools
+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]]:
+	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]]:
+	all_pairs = itertools.combinations(dots, 2)
+	return list(
+		itertools.chain(
+			*[[(dot1, dot2, f) for f in frequency_range] for (dot1, dot2) in all_pairs]
+		)
+	)

pdme/measurement/__init__.py

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

pdme/measurement/dot_measure.py

@@ -0,0 +1,59 @@
+from dataclasses import dataclass
+import numpy
+import numpy.typing
+
+
+@dataclass
+class DotMeasurement:
+	"""
+	Representation of a dot measuring oscillating dipoles.
+
+	Parameters
+	----------
+	v : float
+	The voltage measured at the dot.
+
+	r : numpy.ndarray
+	The position of the dot.
+
+	f : float
+	The measurement frequency.
+	"""
+
+	v: float
+	r: numpy.ndarray
+	f: float
+
+	def __post_init__(self) -> None:
+		self.r = numpy.array(self.r)
+
+
+@dataclass
+class DotRangeMeasurement:
+	"""
+	Representation of a dot measuring oscillating dipoles.
+
+	Parameters
+	----------
+	v_low : float
+	The lower range of voltage measured at the dot.
+
+	v_high : float
+	The upper range of voltage measured at the dot.
+
+	r : numpy.ndarray
+	The position of the dot.
+
+	f : float
+	The measurement frequency.
+	"""
+
+	v_low: float
+	v_high: float
+	r: numpy.ndarray
+	f: float
+
+	def __post_init__(self) -> None:
+		self.r = numpy.array(self.r)
+		if self.v_low > self.v_high:
+			self.v_low, self.v_high = self.v_high, self.v_low

pdme/measurement/dot_pair_measure.py

@@ -0,0 +1,69 @@
+from dataclasses import dataclass
+import numpy
+import numpy.typing
+
+
+@dataclass
+class DotPairMeasurement:
+	"""
+	Representation of a dot measuring oscillating dipoles.
+
+	Parameters
+	----------
+	v : float
+	The voltage measured at the dot.
+
+	r1 : numpy.ndarray
+	The position of the first dot.
+
+	r2 : numpy.ndarray
+	The position of the second dot.
+
+	f : float
+	The measurement frequency.
+	"""
+
+	v: float
+	r1: numpy.ndarray
+	r2: numpy.ndarray
+	f: float
+
+	def __post_init__(self) -> None:
+		self.r1 = numpy.array(self.r1)
+		self.r2 = numpy.array(self.r2)
+
+
+@dataclass
+class DotPairRangeMeasurement:
+	"""
+	Representation of a dot measuring oscillating dipoles.
+
+	Parameters
+	----------
+	v_low : float
+	The lower range of voltage measured at the dot.
+
+	v_high : float
+	The upper range of voltage measured at the dot.
+
+	r1 : numpy.ndarray
+	The position of the first dot.
+
+	r2 : numpy.ndarray
+	The position of the second dot.
+
+	f : float
+	The measurement frequency.
+	"""
+
+	v_low: float
+	v_high: float
+	r1: numpy.ndarray
+	r2: numpy.ndarray
+	f: float
+
+	def __post_init__(self) -> None:
+		self.r1 = numpy.array(self.r1)
+		self.r2 = numpy.array(self.r2)
+		if self.v_low > self.v_high:
+			self.v_low, self.v_high = self.v_high, self.v_low

pdme/measurement/input_types.py

@@ -0,0 +1,36 @@
+import numpy.typing
+from typing import Tuple, Sequence, Union
+from pdme.measurement.dot_measure import DotRangeMeasurement
+from pdme.measurement.dot_pair_measure import DotPairRangeMeasurement
+
+
+DotInput = Tuple[numpy.typing.ArrayLike, float]
+DotPairInput = Tuple[numpy.typing.ArrayLike, numpy.typing.ArrayLike, float]
+
+
+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]
+	)
+
+
+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
+		]
+	)
+
+
+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]
+	highs = [measurement.v_high for measurement in dot_range_measurements]
+	return (numpy.array(lows), numpy.array(highs))

pdme/measurement/oscillating_dipole.py

@@ -0,0 +1,267 @@
+from dataclasses import dataclass
+import numpy
+import numpy.typing
+from typing import Sequence, List
+from pdme.measurement.dot_measure import DotMeasurement, DotRangeMeasurement
+from pdme.measurement.dot_pair_measure import (
+	DotPairMeasurement,
+	DotPairRangeMeasurement,
+)
+import pdme.calculations
+from pdme.measurement.input_types import DotInput, DotPairInput
+
+
+@dataclass
+class OscillatingDipole:
+	"""
+	Representation of an oscillating dipole, either known or guessed.
+
+	Parameters
+	----------
+	p : numpy.ndarray
+	The oscillating dipole moment, with overall sign arbitrary.
+
+	s : numpy.ndarray
+	The position of the dipole.
+
+	w : float
+	The oscillation frequency.
+	"""
+
+	p: numpy.ndarray
+	s: numpy.ndarray
+	w: float
+
+	def __post_init__(self) -> None:
+		"""
+		Coerce the inputs into numpy arrays.
+		"""
+		self.p = numpy.array(self.p)
+		self.s = numpy.array(self.s)
+
+	def _alpha_electric_potential(self, r: numpy.ndarray) -> float:
+		"""
+		Returns the electric potential of this dipole at position r.
+		"""
+		return pdme.calculations.electric_potential(self.p, self.s, r)
+
+	def _alpha_electric_field(self, r: numpy.ndarray) -> numpy.ndarray:
+		"""
+		Returns the electric field of this dipole at position r.
+		"""
+		return pdme.calculations.electric_field(self.p, self.s, r)
+
+	def _b(self, f: float) -> float:
+		return pdme.calculations.telegraph_beta(f, self.w)
+
+	def s_electric_potential_at_position(self, r: numpy.ndarray, f: float) -> float:
+		"""
+		Returns the noise potential at a point r, at some frequency f.
+
+		Specifically for electric potential!
+
+		Parameters
+		----------
+		r : numpy.ndarray
+		The position of the dot.
+
+		f : float
+		The dot frequency to sample.
+		"""
+		return (self._alpha_electric_potential(r)) ** 2 * self._b(f)
+
+	def s_electric_potential_for_dot_pair(
+		self, r1: numpy.ndarray, r2: numpy.ndarray, f: float
+	) -> float:
+		"""
+		This is specifically the analytic cpsd for electric potential noise.
+		This should be deprecated
+		"""
+		return (
+			self._alpha_electric_potential(r1)
+			* self._alpha_electric_potential(r2)
+			* self._b(f)
+		)
+
+	def s_electric_fieldx_at_position(self, r: numpy.ndarray, f: float) -> float:
+		"""
+		Returns the noise potential at a point r, at some frequency f.
+
+		Specifically for electric potential!
+
+		Parameters
+		----------
+		r : numpy.ndarray
+		The position of the dot.
+
+		f : float
+		The dot frequency to sample.
+		"""
+		return (self._alpha_electric_field(r)[0]) ** 2 * self._b(f)
+
+	def s_electric_fieldx_for_dot_pair(
+		self, r1: numpy.ndarray, r2: numpy.ndarray, f: float
+	) -> float:
+		"""
+		This is specifically the analytic cpsd for electric potential noise.
+		This should be deprecated
+		"""
+		return (
+			self._alpha_electric_field(r1)[0]
+			* self._alpha_electric_field(r2)[0]
+			* self._b(f)
+		)
+
+	def to_flat_array(self) -> numpy.ndarray:
+		return numpy.concatenate([self.p, self.s, numpy.array([self.w])])
+
+
+class OscillatingDipoleArrangement:
+	"""
+	A collection of oscillating dipoles, which we are interested in being able to characterise.
+
+	Parameters
+	--------
+	dipoles : Sequence[OscillatingDipole]
+	"""
+
+	def __init__(self, dipoles: Sequence[OscillatingDipole]):
+		self.dipoles = dipoles
+
+	def get_potential_dot_measurement(self, dot_input: DotInput) -> DotMeasurement:
+		r = numpy.array(dot_input[0])
+		f = dot_input[1]
+		return DotMeasurement(
+			sum(
+				[
+					dipole.s_electric_potential_at_position(r, f)
+					for dipole in self.dipoles
+				]
+			),
+			r,
+			f,
+		)
+
+	def get_potential_dot_pair_measurement(
+		self, dot_pair_input: DotPairInput
+	) -> DotPairMeasurement:
+		r1 = numpy.array(dot_pair_input[0])
+		r2 = numpy.array(dot_pair_input[1])
+		f = dot_pair_input[2]
+		return DotPairMeasurement(
+			sum(
+				[
+					dipole.s_electric_potential_for_dot_pair(r1, r2, f)
+					for dipole in self.dipoles
+				]
+			),
+			r1,
+			r2,
+			f,
+		)
+
+	def get_potential_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.
+		"""
+		return [
+			self.get_potential_dot_measurement(dot_input) for dot_input in dot_inputs
+		]
+
+	def get_potential_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.
+		"""
+		return [
+			self.get_potential_dot_pair_measurement(dot_pair_input)
+			for dot_pair_input in dot_pair_inputs
+		]
+
+	def get_percent_range_potential_dot_measurement(
+		self, dot_input: DotInput, low_percent: float, high_percent: float
+	) -> DotRangeMeasurement:
+		r = numpy.array(dot_input[0])
+		f = dot_input[1]
+		return DotRangeMeasurement(
+			low_percent
+			* sum(
+				[
+					dipole.s_electric_potential_at_position(r, f)
+					for dipole in self.dipoles
+				]
+			),
+			high_percent
+			* sum(
+				[
+					dipole.s_electric_potential_at_position(r, f)
+					for dipole in self.dipoles
+				]
+			),
+			r,
+			f,
+		)
+
+	def get_percent_range_potential_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.
+		"""
+		return [
+			self.get_percent_range_potential_dot_measurement(
+				dot_input, low_percent, high_percent
+			)
+			for dot_input in dot_inputs
+		]
+
+	def get_percent_range_potential_dot_pair_measurement(
+		self, pair_input: DotPairInput, low_percent: float, high_percent: float
+	) -> DotPairRangeMeasurement:
+		r1 = numpy.array(pair_input[0])
+		r2 = numpy.array(pair_input[1])
+		f = pair_input[2]
+		return DotPairRangeMeasurement(
+			low_percent
+			* sum(
+				[
+					dipole.s_electric_potential_for_dot_pair(r1, r2, f)
+					for dipole in self.dipoles
+				]
+			),
+			high_percent
+			* sum(
+				[
+					dipole.s_electric_potential_for_dot_pair(r1, r2, f)
+					for dipole in self.dipoles
+				]
+			),
+			r1,
+			r2,
+			f,
+		)
+
+	def get_percent_range_potential_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.
+		"""
+		return [
+			self.get_percent_range_potential_dot_pair_measurement(
+				pair_input, low_percent, high_percent
+			)
+			for pair_input in pair_inputs
+		]
+
+	def to_numpy_array(self) -> numpy.ndarray:
+		"""
+		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])

pdme/meta.py

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

pdme/model/__init__.py

@@ -0,0 +1,30 @@
+from pdme.model.model import DipoleModel
+from pdme.model.fixed_magnitude_model import SingleDipoleFixedMagnitudeModel
+from pdme.model.multidipole_fixed_magnitude_model import (
+	MultipleDipoleFixedMagnitudeModel,
+)
+from pdme.model.random_count_multidipole_fixed_magnitude_model import (
+	RandomCountMultipleDipoleFixedMagnitudeModel,
+)
+
+from pdme.model.log_spaced_random_choice_model import (
+	LogSpacedRandomCountMultipleDipoleFixedMagnitudeModel,
+)
+
+from pdme.model.log_spaced_random_choice_xy_model import (
+	LogSpacedRandomCountMultipleDipoleFixedMagnitudeXYModel,
+)
+
+from pdme.model.log_spaced_random_choice_fixed_orientation_model import (
+	LogSpacedRandomCountMultipleDipoleFixedMagnitudeFixedOrientationModel,
+)
+
+__all__ = [
+	"DipoleModel",
+	"SingleDipoleFixedMagnitudeModel",
+	"MultipleDipoleFixedMagnitudeModel",
+	"RandomCountMultipleDipoleFixedMagnitudeModel",
+	"LogSpacedRandomCountMultipleDipoleFixedMagnitudeModel",
+	"LogSpacedRandomCountMultipleDipoleFixedMagnitudeXYModel",
+	"LogSpacedRandomCountMultipleDipoleFixedMagnitudeFixedOrientationModel",
+]

pdme/model/fixed_magnitude_model.py

@@ -0,0 +1,91 @@
+import numpy
+import numpy.random
+from pdme.model.model import DipoleModel
+from pdme.measurement import (
+	OscillatingDipole,
+	OscillatingDipoleArrangement,
+)
+
+
+class SingleDipoleFixedMagnitudeModel(DipoleModel):
+	"""
+	Model of single oscillating dipole with a fixed magnitude, but free rotation.
+
+	Parameters
+	----------
+	pfixed : float
+	The fixed dipole magnitude.
+	"""
+
+	def __init__(
+		self,
+		xmin: float,
+		xmax: float,
+		ymin: float,
+		ymax: float,
+		zmin: float,
+		zmax: float,
+		pfixed: float,
+	) -> 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()
+
+	def __repr__(self) -> str:
+		return f"SingleDipoleFixedMagnitudeModel({self.xmin}, {self.xmax}, {self.ymin}, {self.ymax}, {self.zmin}, {self.zmax}, {self.pfixed})"
+
+	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
+
+		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)
+		return OscillatingDipoleArrangement(
+			[OscillatingDipole(numpy.array([px, py, pz]), s_pts, frequency)]
+		)
+
+	def get_monte_carlo_dipole_inputs(
+		self, 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 = (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, 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/log_spaced_random_choice_fixed_orientation_model.py

@@ -0,0 +1,195 @@
+import numpy
+import numpy.random
+from pdme.model.model import DipoleModel
+from pdme.measurement import (
+	OscillatingDipole,
+	OscillatingDipoleArrangement,
+)
+import logging
+from typing import Optional
+import pdme.subspace_simulation
+
+_logger = logging.getLogger(__name__)
+
+
+class LogSpacedRandomCountMultipleDipoleFixedMagnitudeFixedOrientationModel(
+	DipoleModel
+):
+	"""
+	Model of multiple oscillating dipoles with a fixed magnitude and fixed rotation. Spaced log uniformly in relaxation time.
+
+	Parameters
+	----------
+
+	wexp_min: log-10 lower bound for dipole frequency
+	wexp_min: log-10 upper bound for dipole frequency
+
+	pfixed : float
+	The fixed dipole magnitude.
+
+	thetafixed: float
+	The fixed theta (polar angle).
+	Should be between 0 and pi.
+
+	phifixed: float
+	The fixed phi (azimuthal angle).
+	Should be between 0 and 2 pi.
+
+	n_max : int
+	The maximum number of dipoles.
+
+	prob_occupancy : float
+	The probability of dipole occupancy
+	"""
+
+	def __init__(
+		self,
+		xmin: float,
+		xmax: float,
+		ymin: float,
+		ymax: float,
+		zmin: float,
+		zmax: float,
+		wexp_min: float,
+		wexp_max: float,
+		pfixed: float,
+		thetafixed: float,
+		phifixed: float,
+		n_max: int,
+		prob_occupancy: float,
+	) -> None:
+		self.xmin = xmin
+		self.xmax = xmax
+		self.ymin = ymin
+		self.ymax = ymax
+		self.zmin = zmin
+		self.zmax = zmax
+		self.wexp_min = wexp_min
+		self.wexp_max = wexp_max
+		self.pfixed = pfixed
+		self.thetafixed = thetafixed
+		self.phifixed = phifixed
+		self.rng = numpy.random.default_rng()
+		self.n_max = n_max
+
+		px = self.pfixed * numpy.sin(self.thetafixed) * numpy.cos(self.phifixed)
+		py = self.pfixed * numpy.sin(self.thetafixed) * numpy.sin(self.phifixed)
+		pz = self.pfixed * numpy.cos(self.thetafixed)
+
+		self.moment_fixed = numpy.array([px, py, pz])
+		if prob_occupancy >= 1 or prob_occupancy <= 0:
+			raise ValueError(
+				f"The probability of a dipole site occupancy must be between 0 and 1, got {prob_occupancy}"
+			)
+		self.prob_occupancy = prob_occupancy
+
+	def __repr__(self) -> str:
+		return f"LogSpacedRandomCountMultipleDipoleFixedMagnitudeFixedOrientationModel({self.xmin}, {self.xmax}, {self.ymin}, {self.ymax}, {self.zmin}, {self.zmax}, {self.wexp_min}, {self.wexp_max}, {self.pfixed}, {self.thetafixed}, {self.phifixed}, {self.n_max}, {self.prob_occupancy})"
+
+	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 = []
+
+		n = rng.binomial(self.n_max, self.prob_occupancy)
+
+		for i in range(n):
+			s_pts = numpy.array(
+				(
+					rng.uniform(self.xmin, self.xmax),
+					rng.uniform(self.ymin, self.ymax),
+					rng.uniform(self.zmin, self.zmax),
+				)
+			)
+			frequency = 10 ** rng.uniform(self.wexp_min, self.wexp_max)
+
+			dipoles.append(OscillatingDipole(self.moment_fixed, s_pts, frequency))
+		return OscillatingDipoleArrangement(dipoles)
+
+	def get_monte_carlo_dipole_inputs(
+		self,
+		monte_carlo_n: int,
+		_: 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_max)
+
+		p_mask = rng.binomial(1, self.prob_occupancy, shape)
+
+		# dipoles = numpy.einsum("ij,k->ijk", p_mask, self.moment_fixed)
+		# Is there a better way to create the final array? probably! can create a flatter guy then reshape.
+		# this is easier to reason about.
+		p_magnitude = self.pfixed * p_mask
+
+		px = p_magnitude * numpy.sin(self.thetafixed) * numpy.cos(self.phifixed)
+		py = p_magnitude * numpy.sin(self.thetafixed) * numpy.sin(self.phifixed)
+		pz = p_magnitude * numpy.cos(self.thetafixed)
+
+		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 = 10 ** rng.uniform(self.wexp_min, self.wexp_max, shape)
+
+		return numpy.stack([px, py, pz, sx, sy, sz, w], axis=-1)
+
+	def markov_chain_monte_carlo_proposal(
+		self,
+		dipole: numpy.ndarray,
+		stdev: pdme.subspace_simulation.DipoleStandardDeviation,
+		rng_arg: Optional[numpy.random.Generator] = None,
+	) -> numpy.ndarray:
+		if rng_arg is None:
+			rng_to_use = self.rng
+		else:
+			rng_to_use = rng_arg
+
+		px = dipole[0]
+		py = dipole[1]
+		pz = dipole[2]
+		# won't change p for this model of fixed dipole moment.
+
+		rx = dipole[3]
+		ry = dipole[4]
+		rz = dipole[5]
+
+		tentative_rx = rx + stdev.rx_step * rng_to_use.uniform(-1, 1)
+		if tentative_rx < self.xmin or tentative_rx > self.xmax:
+			tentative_rx = rx
+
+		tentative_ry = ry + stdev.ry_step * rng_to_use.uniform(-1, 1)
+		if tentative_ry < self.ymin or tentative_ry > self.ymax:
+			tentative_ry = ry
+		tentative_rz = rz + stdev.rz_step * rng_to_use.uniform(-1, 1)
+		if tentative_rz < self.zmin or tentative_rz > self.zmax:
+			tentative_rz = rz
+
+		w = dipole[6]
+		tentative_w = numpy.exp(
+			numpy.log(w) + (stdev.w_log_step * rng_to_use.uniform(-1, 1))
+		)
+		tentative_dip = numpy.array(
+			[
+				px,
+				py,
+				pz,
+				tentative_rx,
+				tentative_ry,
+				tentative_rz,
+				tentative_w,
+			]
+		)
+		return tentative_dip

pdme/model/log_spaced_random_choice_model.py

@@ -0,0 +1,199 @@
+import numpy
+import numpy.random
+from pdme.model.model import DipoleModel
+from pdme.measurement import (
+	OscillatingDipole,
+	OscillatingDipoleArrangement,
+)
+import pdme.subspace_simulation
+from typing import Optional
+
+
+class LogSpacedRandomCountMultipleDipoleFixedMagnitudeModel(DipoleModel):
+	"""
+	Model of multiple oscillating dipoles with a fixed magnitude, but free rotation. Spaced logarithmically.
+
+	Parameters
+	----------
+
+	wexp_min: log-10 lower bound for dipole frequency
+	wexp_min: log-10 upper bound for dipole frequency
+
+	pfixed : float
+	The fixed dipole magnitude.
+
+	n_max : int
+	The maximum number of dipoles.
+
+	prob_occupancy : float
+	The probability of dipole occupancy
+	"""
+
+	def __init__(
+		self,
+		xmin: float,
+		xmax: float,
+		ymin: float,
+		ymax: float,
+		zmin: float,
+		zmax: float,
+		wexp_min: float,
+		wexp_max: float,
+		pfixed: float,
+		n_max: int,
+		prob_occupancy: float,
+	) -> None:
+		self.xmin = xmin
+		self.xmax = xmax
+		self.ymin = ymin
+		self.ymax = ymax
+		self.zmin = zmin
+		self.zmax = zmax
+		self.wexp_min = wexp_min
+		self.wexp_max = wexp_max
+		self.pfixed = pfixed
+		self.rng = numpy.random.default_rng()
+		self.n_max = n_max
+		if prob_occupancy >= 1 or prob_occupancy <= 0:
+			raise ValueError(
+				f"The probability of a dipole site occupancy must be between 0 and 1, got {prob_occupancy}"
+			)
+		self.prob_occupancy = prob_occupancy
+
+	def __repr__(self) -> str:
+		return f"LogSpacedRandomCountMultipleDipoleFixedMagnitudeModel({self.xmin}, {self.xmax}, {self.ymin}, {self.ymax}, {self.zmin}, {self.zmax}, {self.wexp_min}, {self.wexp_max}, {self.pfixed}, {self.n_max}, {self.prob_occupancy})"
+
+	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 = []
+
+		n = rng.binomial(self.n_max, self.prob_occupancy)
+
+		for i in range(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 = 10 ** rng.uniform(self.wexp_min, self.wexp_max)
+
+			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,
+		_: 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_max)
+		theta = 2 * numpy.pi * rng.random(shape)
+		phi = numpy.arccos(2 * rng.random(shape) - 1)
+
+		p_mask = rng.binomial(1, self.prob_occupancy, shape)
+		p_magnitude = self.pfixed * p_mask
+
+		px = p_magnitude * numpy.cos(theta) * numpy.sin(phi)
+		py = p_magnitude * numpy.sin(theta) * numpy.sin(phi)
+		pz = p_magnitude * 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 = 10 ** rng.uniform(self.wexp_min, self.wexp_max, shape)
+
+		return numpy.stack([px, py, pz, sx, sy, sz, w], axis=-1)
+
+	def markov_chain_monte_carlo_proposal(
+		self,
+		dipole: numpy.ndarray,
+		stdev: pdme.subspace_simulation.DipoleStandardDeviation,
+		rng_arg: Optional[numpy.random.Generator] = None,
+	) -> numpy.ndarray:
+		if rng_arg is None:
+			rng_to_use = self.rng
+		else:
+			rng_to_use = rng_arg
+
+		px = dipole[0]
+		py = dipole[1]
+		pz = dipole[2]
+		theta = numpy.arccos(pz / self.pfixed)
+		phi = numpy.arctan2(py, px)
+
+		# need to step phi, theta, rx, ry, rz, w
+		# then p^\ast is 1/(2 phi_step) and Delta = phi_step(2 * {0, 1} - 1)
+		delta_phi = stdev.p_phi_step * rng_to_use.uniform(-1, 1)
+		tentative_phi = phi + delta_phi
+
+		# theta
+		delta_theta = stdev.p_theta_step * rng_to_use.uniform(-1, 1)
+		r = (numpy.sin(theta + delta_theta)) / (numpy.sin(theta))
+		if r > rng_to_use.uniform(0, 1):
+			tentative_theta = theta + delta_theta
+		else:
+			tentative_theta = theta
+
+		tentative_px = (
+			self.pfixed * numpy.sin(tentative_theta) * numpy.cos(tentative_phi)
+		)
+		tentative_py = (
+			self.pfixed * numpy.sin(tentative_theta) * numpy.sin(tentative_phi)
+		)
+		tentative_pz = self.pfixed * numpy.cos(tentative_theta)
+
+		rx = dipole[3]
+		ry = dipole[4]
+		rz = dipole[5]
+
+		tentative_rx = rx + stdev.rx_step * rng_to_use.uniform(-1, 1)
+		if tentative_rx < self.xmin or tentative_rx > self.xmax:
+			tentative_rx = rx
+
+		tentative_ry = ry + stdev.ry_step * rng_to_use.uniform(-1, 1)
+		if tentative_ry < self.ymin or tentative_ry > self.ymax:
+			tentative_ry = ry
+		tentative_rz = rz + stdev.rz_step * rng_to_use.uniform(-1, 1)
+		if tentative_rz < self.zmin or tentative_rz > self.zmax:
+			tentative_rz = rz
+
+		w = dipole[6]
+		tentative_w = numpy.exp(
+			numpy.log(w) + (stdev.w_log_step * rng_to_use.uniform(-1, 1))
+		)
+		tentative_dip = numpy.array(
+			[
+				tentative_px,
+				tentative_py,
+				tentative_pz,
+				tentative_rx,
+				tentative_ry,
+				tentative_rz,
+				tentative_w,
+			]
+		)
+		return tentative_dip

pdme/model/log_spaced_random_choice_xy_model.py

@@ -0,0 +1,183 @@
+import numpy
+import numpy.random
+from pdme.model.model import DipoleModel
+from pdme.measurement import (
+	OscillatingDipole,
+	OscillatingDipoleArrangement,
+)
+import pdme.subspace_simulation
+from typing import Optional
+
+
+class LogSpacedRandomCountMultipleDipoleFixedMagnitudeXYModel(DipoleModel):
+	"""
+	Model of multiple oscillating dipoles with a fixed magnitude, but free rotation in XY plane. Spaced logarithmically.
+
+	Parameters
+	----------
+
+	wexp_min: log-10 lower bound for dipole frequency
+	wexp_min: log-10 upper bound for dipole frequency
+
+	pfixed : float
+	The fixed dipole magnitude.
+
+	n_max : int
+	The maximum number of dipoles.
+
+	prob_occupancy : float
+	The probability of dipole occupancy
+	"""
+
+	def __init__(
+		self,
+		xmin: float,
+		xmax: float,
+		ymin: float,
+		ymax: float,
+		zmin: float,
+		zmax: float,
+		wexp_min: float,
+		wexp_max: float,
+		pfixed: float,
+		n_max: int,
+		prob_occupancy: float,
+	) -> None:
+		self.xmin = xmin
+		self.xmax = xmax
+		self.ymin = ymin
+		self.ymax = ymax
+		self.zmin = zmin
+		self.zmax = zmax
+		self.wexp_min = wexp_min
+		self.wexp_max = wexp_max
+		self.pfixed = pfixed
+		self.rng = numpy.random.default_rng()
+		self.n_max = n_max
+		if prob_occupancy >= 1 or prob_occupancy <= 0:
+			raise ValueError(
+				f"The probability of a dipole site occupancy must be between 0 and 1, got {prob_occupancy}"
+			)
+		self.prob_occupancy = prob_occupancy
+
+	def __repr__(self) -> str:
+		return f"LogSpacedRandomCountMultipleDipoleFixedMagnitudeXYModel({self.xmin}, {self.xmax}, {self.ymin}, {self.ymax}, {self.zmin}, {self.zmax}, {self.wexp_min}, {self.wexp_max}, {self.pfixed}, {self.n_max}, {self.prob_occupancy})"
+
+	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 = []
+
+		n = rng.binomial(self.n_max, self.prob_occupancy)
+
+		for i in range(n):
+			phi = rng.uniform(0, 2 * numpy.pi)
+			px = self.pfixed * numpy.cos(phi)
+			py = self.pfixed * numpy.sin(phi)
+			pz = 0
+			s_pts = numpy.array(
+				(
+					rng.uniform(self.xmin, self.xmax),
+					rng.uniform(self.ymin, self.ymax),
+					rng.uniform(self.zmin, self.zmax),
+				)
+			)
+			frequency = 10 ** rng.uniform(self.wexp_min, self.wexp_max)
+
+			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,
+		_: 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_max)
+		phi = 2 * numpy.pi * rng.random(shape)
+
+		p_mask = rng.binomial(1, self.prob_occupancy, shape)
+		p_magnitude = self.pfixed * p_mask
+
+		px = p_magnitude * numpy.cos(phi)
+		py = p_magnitude * numpy.sin(phi)
+		pz = p_magnitude * 0
+
+		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 = 10 ** rng.uniform(self.wexp_min, self.wexp_max, shape)
+
+		return numpy.stack([px, py, pz, sx, sy, sz, w], axis=-1)
+
+	def markov_chain_monte_carlo_proposal(
+		self,
+		dipole: numpy.ndarray,
+		stdev: pdme.subspace_simulation.DipoleStandardDeviation,
+		rng_arg: Optional[numpy.random.Generator] = None,
+	) -> numpy.ndarray:
+		if rng_arg is None:
+			rng_to_use = self.rng
+		else:
+			rng_to_use = rng_arg
+
+		px = dipole[0]
+		py = dipole[1]
+		pz = dipole[2]
+		phi = numpy.arctan2(py, px)
+
+		# need to step phi, rx, ry, rz, w
+		# then p^\ast is 1/(2 phi_step) and Delta = phi_step(2 * {0, 1} - 1)
+		delta_phi = stdev.p_phi_step * rng_to_use.uniform(-1, 1)
+		tentative_phi = phi + delta_phi
+
+		tentative_px = self.pfixed * numpy.cos(tentative_phi)
+		tentative_py = self.pfixed * numpy.sin(tentative_phi)
+
+		rx = dipole[3]
+		ry = dipole[4]
+		rz = dipole[5]
+
+		tentative_rx = rx + stdev.rx_step * rng_to_use.uniform(-1, 1)
+		if tentative_rx < self.xmin or tentative_rx > self.xmax:
+			tentative_rx = rx
+
+		tentative_ry = ry + stdev.ry_step * rng_to_use.uniform(-1, 1)
+		if tentative_ry < self.ymin or tentative_ry > self.ymax:
+			tentative_ry = ry
+		tentative_rz = rz + stdev.rz_step * rng_to_use.uniform(-1, 1)
+		if tentative_rz < self.zmin or tentative_rz > self.zmax:
+			tentative_rz = rz
+
+		w = dipole[6]
+		tentative_w = numpy.exp(
+			numpy.log(w) + (stdev.w_log_step * rng_to_use.uniform(-1, 1))
+		)
+		tentative_dip = numpy.array(
+			[
+				tentative_px,
+				tentative_py,
+				pz,
+				tentative_rx,
+				tentative_ry,
+				tentative_rz,
+				tentative_w,
+			]
+		)
+		return tentative_dip

pdme/model/model.py

@@ -0,0 +1,163 @@
+import numpy
+import numpy.random
+from pdme.measurement import (
+	OscillatingDipoleArrangement,
+)
+import logging
+import pdme.subspace_simulation
+from typing import List, Tuple, Optional
+
+_logger = logging.getLogger(__name__)
+
+
+class DipoleModel:
+	"""
+	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 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 get_monte_carlo_dipole_inputs(
+		self, n: int, max_frequency: float, rng: numpy.random.Generator = None
+	) -> numpy.ndarray:
+		"""
+		For a given DipoleModel, gets a set of dipole collections as a monte_carlo_n x dipole_count x 7 numpy array.
+		"""
+		raise NotImplementedError
+
+	def markov_chain_monte_carlo_proposal(
+		self,
+		dipole: numpy.ndarray,
+		stdev: pdme.subspace_simulation.DipoleStandardDeviation,
+		rng_arg: Optional[numpy.random.Generator] = None,
+	) -> numpy.ndarray:
+		raise NotImplementedError
+
+	def get_mcmc_chain(
+		self,
+		seed,
+		cost_function,
+		chain_length,
+		threshold_cost: float,
+		stdevs: pdme.subspace_simulation.MCMCStandardDeviation,
+		initial_cost: Optional[float] = None,
+		rng_arg: Optional[numpy.random.Generator] = None,
+	) -> List[Tuple[float, numpy.ndarray]]:
+		"""
+		performs constrained markov chain monte carlo starting on seed parameter.
+		The cost function given is used as a constrained to condition the chain;
+		a new state is only accepted if cost_function(state) < cost_function(previous_state).
+		The stdevs passed in are the stdevs we're expected to use.
+
+		Because we're using this for subspace simulation where our proposal function is not too important, we're in good shape.
+		Note that for our adaptive stdevs to work, there's an unwritten contract that we sort each dipole in the state by frequency (increasing).
+
+		The seed is a list of dipoles, and each chain state is a list of dipoles as well.
+
+		initial_cost is a performance guy that lets you pre-populate the initial cost used to define the condition.
+		Probably premature optimisation.
+
+		Returns a chain of [ (cost: float, state: dipole_ndarray ) ] format.
+		"""
+		_logger.debug(
+			f"Starting Markov Chain Monte Carlo with seed: {seed} for chain length {chain_length} and provided stdevs {stdevs}"
+		)
+		chain: List[Tuple[float, numpy.ndarray]] = []
+		if initial_cost is None:
+			current_cost = cost_function(numpy.array([seed]))
+		else:
+			current_cost = initial_cost
+		current = seed
+		for i in range(chain_length):
+			dips = []
+			for dipole_index, dipole in enumerate(current):
+				_logger.debug(dipole_index)
+				_logger.debug(dipole)
+				stdev = stdevs[dipole_index]
+				tentative_dip = self.markov_chain_monte_carlo_proposal(
+					dipole, stdev, rng_arg
+				)
+
+				dips.append(tentative_dip)
+			dips_array = pdme.subspace_simulation.sort_array_of_dipoles_by_frequency(
+				dips
+			)
+			tentative_cost = cost_function(numpy.array([dips_array]))[0]
+			if tentative_cost < threshold_cost:
+				chain.append((numpy.squeeze(tentative_cost).item(), dips_array))
+				current = dips_array
+				current_cost = tentative_cost
+			else:
+				chain.append((numpy.squeeze(current_cost).item(), current))
+		return chain
+
+	def get_repeat_counting_mcmc_chain(
+		self,
+		seed,
+		cost_function,
+		chain_length,
+		threshold_cost: float,
+		stdevs: pdme.subspace_simulation.MCMCStandardDeviation,
+		initial_cost: Optional[float] = None,
+		rng_arg: Optional[numpy.random.Generator] = None,
+	) -> Tuple[int, List[Tuple[float, numpy.ndarray]]]:
+		"""
+		performs constrained markov chain monte carlo starting on seed parameter.
+		The cost function given is used as a constrained to condition the chain;
+		a new state is only accepted if cost_function(state) < cost_function(previous_state).
+		The stdevs passed in are the stdevs we're expected to use.
+
+		Because we're using this for subspace simulation where our proposal function is not too important, we're in good shape.
+		Note that for our adaptive stdevs to work, there's an unwritten contract that we sort each dipole in the state by frequency (increasing).
+
+		The seed is a list of dipoles, and each chain state is a list of dipoles as well.
+
+		initial_cost is a performance guy that lets you pre-populate the initial cost used to define the condition.
+		Probably premature optimisation.
+
+		Chain has type of [ (cost: float, state: dipole_ndarray ) ] format,
+		returning (repeat_count, chain) to keep track of number of repeats
+		"""
+		_logger.debug(
+			f"Starting Markov Chain Monte Carlo with seed: {seed} for chain length {chain_length} and provided stdevs {stdevs}"
+		)
+		chain: List[Tuple[float, numpy.ndarray]] = []
+		if initial_cost is None:
+			current_cost = cost_function(numpy.array([seed]))
+		else:
+			current_cost = initial_cost
+		current = seed
+		repeat_event_count = 0
+		for _ in range(chain_length):
+			dips = []
+			for dipole_index, dipole in enumerate(current):
+				_logger.debug(dipole_index)
+				_logger.debug(dipole)
+				stdev = stdevs[dipole_index]
+				tentative_dip = self.markov_chain_monte_carlo_proposal(
+					dipole, stdev, rng_arg
+				)
+
+				dips.append(tentative_dip)
+			dips_array = pdme.subspace_simulation.sort_array_of_dipoles_by_frequency(
+				dips
+			)
+			tentative_cost = cost_function(numpy.array([dips_array]))[0]
+			if tentative_cost < threshold_cost:
+				chain.append((numpy.squeeze(tentative_cost).item(), dips_array))
+				current = dips_array
+				current_cost = tentative_cost
+			else:
+				# repeating a sample, increase count
+				repeat_event_count += 1
+				chain.append((numpy.squeeze(current_cost).item(), current))
+		return (repeat_event_count, chain)

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/random_count_multidipole_fixed_magnitude_model.py

@@ -0,0 +1,119 @@
+import numpy
+import numpy.random
+from pdme.model.model import DipoleModel
+from pdme.measurement import (
+	OscillatingDipole,
+	OscillatingDipoleArrangement,
+)
+
+
+class RandomCountMultipleDipoleFixedMagnitudeModel(DipoleModel):
+	"""
+	Model of multiple oscillating dipoles with a fixed magnitude, but free rotation.
+
+	Parameters
+	----------
+	pfixed : float
+	The fixed dipole magnitude.
+
+	n_max : int
+	The maximum number of dipoles.
+
+	prob_occupancy : float
+	The probability of dipole occupancy
+	"""
+
+	def __init__(
+		self,
+		xmin: float,
+		xmax: float,
+		ymin: float,
+		ymax: float,
+		zmin: float,
+		zmax: float,
+		pfixed: float,
+		n_max: int,
+		prob_occupancy: float,
+	) -> 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_max = n_max
+		if prob_occupancy >= 1 or prob_occupancy <= 0:
+			raise ValueError(
+				f"The probability of a dipole site occupancy must be between 0 and 1, got {prob_occupancy}"
+			)
+		self.prob_occupancy = prob_occupancy
+
+	def __repr__(self) -> str:
+		return f"RandomCountMultipleDipoleFixedMagnitudeModel({self.xmin}, {self.xmax}, {self.ymin}, {self.ymax}, {self.zmin}, {self.zmax}, {self.pfixed}, {self.n_max}, {self.prob_occupancy})"
+
+	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 = []
+
+		n = rng.binomial(self.n_max, self.prob_occupancy)
+
+		for i in range(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_max)
+		theta = 2 * numpy.pi * rng.random(shape)
+		phi = numpy.arccos(2 * rng.random(shape) - 1)
+
+		p_mask = rng.binomial(1, self.prob_occupancy, shape)
+		p_magnitude = self.pfixed * p_mask
+
+		px = p_magnitude * numpy.cos(theta) * numpy.sin(phi)
+		py = p_magnitude * numpy.sin(theta) * numpy.sin(phi)
+		pz = p_magnitude * 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/subspace_simulation/__init__.py

@@ -0,0 +1,74 @@
+from dataclasses import dataclass
+from typing import Sequence
+import numpy
+from pdme.subspace_simulation.mcmc_costs import (
+	proportional_cost,
+	proportional_costs_vs_actual_measurement,
+)
+
+
+@dataclass
+class DipoleStandardDeviation:
+	"""
+	contains the dipole standard deviation to be used in porposals for markov chain monte carlo
+	"""
+
+	p_phi_step: float
+	p_theta_step: float
+	rx_step: float
+	ry_step: float
+	rz_step: float
+	w_log_step: float
+
+
+class MCMCStandardDeviation:
+	"""
+	wrapper for multiple standard deviations, allows for flexible length stuff
+	"""
+
+	def __init__(self, stdevs: Sequence[DipoleStandardDeviation]):
+		self.stdevs = stdevs
+		if len(stdevs) < 1:
+			raise ValueError(f"Got stdevs: {stdevs}, must have length > 1")
+
+	def __getitem__(self, key):
+		newkey = key % len(self.stdevs)
+		return self.stdevs[newkey]
+
+
+def sort_array_of_dipoles_by_frequency(configuration) -> numpy.ndarray:
+	"""
+	Say we have a situation of 2 dipoles, and we've created 8 samples. Then we'll have an (8, 2, 7) numpy array.
+	For each of the 8 samples, we want the 2 dipoles to be in order of frequency.
+
+	This just sorts each sample, the 2x7 array.
+
+	Utility function.
+	"""
+	return numpy.array(sorted(configuration, key=lambda l: l[6]))
+
+
+def sort_array_of_dipoleses_by_frequency(configurations) -> numpy.ndarray:
+	"""
+	Say we have a situation of 2 dipoles, and we've created 8 samples. Then we'll have an (8, 2, 7) numpy array.
+	For each of the 8 samples, we want the 2 dipoles to be in order of frequency.
+
+	This is the wrapper that sorts everything.
+
+	Utility function.
+	"""
+	return numpy.array(
+		[
+			sort_array_of_dipoles_by_frequency(configuration)
+			for configuration in configurations
+		]
+	)
+
+
+__all__ = [
+	"DipoleStandardDeviation",
+	"MCMCStandardDeviation",
+	"sort_array_of_dipoles_by_frequency",
+	"proportional_cost",
+	"proportional_costs_vs_actual_measurement",
+]

pdme/subspace_simulation/mcmc_costs.py

@@ -0,0 +1,30 @@
+import numpy
+import numpy.typing
+import pdme.util.fast_v_calc
+
+
+def proportional_cost(a: numpy.ndarray, b: numpy.ndarray) -> numpy.ndarray:
+	tops = numpy.max(b / a, axis=-1)
+	bottoms = numpy.max(a / b, axis=-1)
+	return numpy.maximum(tops, bottoms)
+
+
+def proportional_costs_vs_actual_measurement(
+	dot_inputs_array: numpy.ndarray,
+	actual_measurement_array: numpy.ndarray,
+	dipoles_to_test: numpy.ndarray,
+) -> numpy.ndarray:
+	vals = pdme.util.fast_v_calc.fast_vs_for_dipoleses(
+		dot_inputs_array, dipoles_to_test
+	)
+	return proportional_cost(actual_measurement_array, vals)
+
+
+def relative_square_diffs(
+	approx: numpy.ndarray, target: numpy.ndarray
+) -> numpy.ndarray:
+	# Assume that both approx and target are arrays of length m
+	# Approx can broadcast if additional indexes to the left
+	# diffs.shape = [ m ]
+	diffs = (approx - target) ** 2 / (target**2)
+	return diffs.sum(axis=-1)

pdme/util/fast_nonlocal_spectrum.py

@@ -0,0 +1,273 @@
+import numpy
+import logging
+
+
+_logger = logging.getLogger(__name__)
+
+
+def fast_s_nonlocal(
+	dot_pair_inputs: numpy.ndarray, dipoles: numpy.ndarray
+) -> numpy.ndarray:
+	"""
+	No error correction here baby.
+	"""
+	ps = dipoles[:, 0:3]
+	ss = dipoles[:, 3:6]
+	ws = dipoles[:, 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 - ss[:, None]
+	diffses2 = r2s - ss[:, None]
+	if _logger.isEnabledFor(logging.DEBUG):
+		_logger.debug(f"diffses1: {diffses1}")
+		_logger.debug(f"diffses2: {diffses2}")
+
+	norms1 = numpy.linalg.norm(diffses1, axis=2) ** 3
+	norms2 = numpy.linalg.norm(diffses2, axis=2) ** 3
+	if _logger.isEnabledFor(logging.DEBUG):
+		_logger.debug(f"norms1: {norms1}")
+		_logger.debug(f"norms2: {norms2}")
+
+	alphses1 = numpy.einsum("...ji, ...i", diffses1, ps) / norms1
+	alphses2 = numpy.einsum("...ji, ...i", diffses2, ps) / norms2
+	if _logger.isEnabledFor(logging.DEBUG):
+		_logger.debug(f"alphses1: {alphses1}")
+		_logger.debug(f"alphses2: {alphses2}")
+
+	bses = 2 * ws[:, None] / ((numpy.pi * f1s) ** 2 + ws[:, None] ** 2)
+	if _logger.isEnabledFor(logging.DEBUG):
+		_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.
+	"""
+
+	# We're going to annotate the indices on this class.
+	# Let's define some indices:
+	# A -> index of dipoleses configurations
+	# measurement_index -> if we have 100 frequencies for example, indexes which one of them it is
+	# j -> within a particular configuration, indexes dipole j
+	# If we need to use numbers, let's use A -> 2, j -> 10, measurement_index -> 9 for consistency with
+	# my other notes
+	# cart -> {x, y, z} is a cartesian axis
+
+	# ps, ss have shape [A, j, cart]
+	ps = dipoleses[:, :, 0:3]
+	ss = dipoleses[:, :, 3:6]
+	# ws shape [A, j]
+	ws = dipoleses[:, :, 6]
+
+	_logger.debug(f"ps: {ps}")
+	_logger.debug(f"ss: {ss}")
+	_logger.debug(f"ws: {ws}")
+
+	# rs have shape [meas_idx, {}, cart], where the inner index goes away leaving
+	# [meas, cart]
+	r1s = dot_pair_inputs[:, 0, 0:3]
+	r2s = dot_pair_inputs[:, 1, 0:3]
+	# fs have index [meas_idx], this makes sense
+	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}")
+
+	# r1s have shape [meas, cart], adding the none makes it
+	# r1s[:, None].shape = [meas, 1, cart]
+	# ss[:, None, :].shape = [A, 1, j, cart]
+	# subtracting broadcasts by matching from the right to the left, giving a final shape of
+	# diffses.shape [A, meas, j, cart]
+	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}")
+
+	# norming on the cartesian axis, which is axis 3 as seen above
+	# norms.shape [A, meas, j]
+	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}")
+
+	# diffses shape [A, meas, j, cart]
+	# ps shape [A, j, cart]
+	# so we're summing over the d axis, the cartesian one.
+	# final shape of numerator is [A, meas, j]
+	# denom shape is [A, meas, j]
+	# final shape stayes [A, meas, j]
+	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}")
+
+	# ws shape [A, j], so numerator has shape [A, 1, j]
+	# f1s shape is [meas], so first term of denom is [meas, 1]
+	# ws[:, None, :].shape [A, 1, j] so breadcasting the sum in denom gives
+	# denom.shape [A meas, j]
+	# so now overall shape is [A, meas, j]
+	bses = 2 * ws[:, None, :] / ((numpy.pi * f1s[:, None]) ** 2 + ws[:, None, :] ** 2)
+	if _logger.isEnabledFor(logging.DEBUG):
+		_logger.debug(f"bses: {bses}")
+
+	# so our output shape is [A, meas, j]
+	_logger.debug(f"Raw pair calc: [{alphses1 * alphses2 * bses}]")
+	return numpy.einsum("...j->...", alphses1 * alphses2 * bses)
+
+
+def fast_s_spin_qubit_tarucha_nonlocal_dipoleses(
+	dot_pair_inputs: numpy.ndarray, dipoleses: numpy.ndarray
+) -> numpy.ndarray:
+	"""
+	No error correction here baby.
+	"""
+
+	# We're going to annotate the indices on this class.
+	# Let's define some indices:
+	# A -> index of dipoleses configurations
+	# j -> within a particular configuration, indexes dipole j
+	# measurement_index -> if we have 100 frequencies for example, indexes which one of them it is
+	# If we need to use numbers, let's use A -> 2, j -> 10, measurement_index -> 9 for consistency with
+	# my other notes
+
+	# axes are [dipole_config_idx A, dipole_idx j, {px, py, pz}3]
+	ps = dipoleses[:, :, 0:3]
+	# axes are [dipole_config_idx A, dipole_idx j, {sx, sy, sz}3]
+	ss = dipoleses[:, :, 3:6]
+	# axes are [dipole_config_idx A, dipole_idx j, w], last axis is just 1
+	ws = dipoleses[:, :, 6]
+
+	_logger.debug(f"ps: {ps}")
+	_logger.debug(f"ss: {ss}")
+	_logger.debug(f"ws: {ws}")
+
+	# dot_index is either 0 or 1 for dot1 or dot2
+	# hopefully this adhoc grammar is making sense, with the explicit labelling of the values of the last axis in cartesian space
+	# axes are [measurement_idx, {dot_index}, {rx, ry, rz}] where the inner {dot_index} is gone
+	# [measurement_idx, cartesian3]
+	r1s = dot_pair_inputs[:, 0, 0:3]
+	r2s = dot_pair_inputs[:, 1, 0:3]
+	# axes are [measurement_idx]
+	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}")
+
+	# first operation!
+	# r1s has shape [measurement_idx, rxs]
+	# None inserts an extra axis so the r1s[:, None] has shape
+	# [measurement_idx, 1]([rxs]) with the last rxs hidden
+	#
+	# ss has shape [ A, j, {sx, sy, sz}3], so second term has shape [A, 1, j]([sxs])
+	# these broadcast from right to left
+	# [	 measurement_idx, 1, rxs]
+	# [A,	  1,			   j, sxs]
+	# resulting in [A, measurement_idx, j, cart3] sxs rxs are both cart3
+	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}")
+
+	# norms takes out axis 3, the last one, giving [A, measurement_idx, j]
+	norms1 = numpy.linalg.norm(diffses1, axis=3)
+	norms2 = numpy.linalg.norm(diffses2, axis=3)
+	if _logger.isEnabledFor(logging.DEBUG):
+		_logger.debug(f"norms1: {norms1}")
+		_logger.debug(f"norms2: {norms2}")
+
+	# _logger.info(f"norms1: {norms1}")
+	# _logger.info(f"norms1 shape: {norms1.shape}")
+	#
+	# diffses1 (A, measurement_idx, j, xs)
+	# ps:  (A, j, px)
+	# result is (A, measurement_idx, j)
+	# intermediate_dot_prod = numpy.einsum("abcd,acd->abc", diffses1, ps)
+	# _logger.info(f"dot product shape: {intermediate_dot_prod.shape}")
+
+	# transpose makes it (j, measurement_idx, A)
+	# transp_intermediate_dot_prod = numpy.transpose(numpy.einsum("abcd,acd->abc", diffses1, ps) / (norms1**3))
+
+	# transpose of diffses has shape (xs, j, measurement_idx, A)
+	# numpy.transpose(diffses1)
+	# _logger.info(f"dot product shape: {transp_intermediate_dot_prod.shape}")
+
+	# inner transpose is (j, measurement_idx, A) * (xs, j, measurement_idx, A)
+	# next transpose puts it back to (A, measurement_idx, j, xs)
+	# p_dot_r_times_r_term = 3 * numpy.transpose(numpy.transpose(numpy.einsum("abcd,acd->abc", diffses1, ps) / (norms1**3)) * numpy.transpose(diffses1))
+	# _logger.info(f"p_dot_r_times_r_term: {p_dot_r_times_r_term.shape}")
+
+	# only x axis puts us at (A, measurement_idx, j)
+	# p_dot_r_times_r_term_x_only = p_dot_r_times_r_term[:, :, :, 0]
+	# _logger.info(f"p_dot_r_times_r_term_x_only.shape: {p_dot_r_times_r_term_x_only.shape}")
+
+	# now to complete the numerator we subtract the ps, which are (A, j, px):
+	# slicing off the end gives us (A, j), so we newaxis to get (A, 1, j)
+	# _logger.info(ps[:, numpy.newaxis, :, 0].shape)
+	alphses1 = (
+		(
+			3
+			* numpy.transpose(
+				numpy.transpose(
+					numpy.einsum("abcd,acd->abc", diffses1, ps) / (norms1**2)
+				)
+				* numpy.transpose(diffses1)
+			)[:, :, :, 0]
+		)
+		- ps[:, numpy.newaxis, :, 0]
+	) / (norms1**3)
+	alphses2 = (
+		(
+			3
+			* numpy.transpose(
+				numpy.transpose(
+					numpy.einsum("abcd,acd->abc", diffses2, ps) / (norms2**2)
+				)
+				* numpy.transpose(diffses2)
+			)[:, :, :, 0]
+		)
+		- ps[:, numpy.newaxis, :, 0]
+	) / (norms2**3)
+	if _logger.isEnabledFor(logging.DEBUG):
+		_logger.debug(f"alphses1: {alphses1}")
+		_logger.debug(f"alphses2: {alphses2}")
+
+	# ws has shape (A, j), so it becomes (A, 1, j) in numerator with the new axis
+	# f1s has shape (m), so we get in the denominator (m, 1) + (A, 1, j)
+	# This becomes (A, m, j)
+	bses = 2 * ws[:, None, :] / ((numpy.pi * f1s[:, None]) ** 2 + ws[:, None, :] ** 2)
+	if _logger.isEnabledFor(logging.DEBUG):
+		_logger.debug(f"bses: {bses}")
+
+	# alphas have (A, 1, j), betas have (A, m, j)
+	# Final result is (A, m, j)
+	_logger.debug(f"Raw pair calc: [{alphses1 * alphses2 * bses}]")
+	return numpy.einsum("...j->...", alphses1 * alphses2 * bses)
+
+
+def signarg(x, **kwargs):
+	"""
+	uses numpy.sign to implement Arg for real numbers only. Should return pi for negative inputs, 0 for positive.
+	Passes through args to numpy.sign
+	"""
+	return numpy.pi * (numpy.sign(x, **kwargs) - 1) / (-2)

pdme/util/fast_v_calc.py

@@ -0,0 +1,252 @@
+import numpy
+import logging
+
+
+_logger = logging.getLogger(__name__)
+
+
+def fast_vs_for_dipoles(
+	dot_inputs: numpy.ndarray, dipoles: 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.
+	"""
+	# name indexes:
+	# A: dipole config index
+	# cart: cartesian index
+	# m: measurement index
+
+	# [A, cart]
+	ps = dipoles[:, 0:3]
+	ss = dipoles[:, 3:6]
+	# [A]
+	ws = dipoles[:, 6]
+
+	_logger.debug(f"ps: {ps}")
+	_logger.debug(f"ss: {ss}")
+	_logger.debug(f"ws: {ws}")
+
+	# [m, cart]
+	rs = dot_inputs[:, 0:3]
+	# [m]
+	fs = dot_inputs[:, 3]
+
+	# [m, cart] - [A, 1, cart]
+	# [A, m, cart]
+	diffses = rs - ss[:, None]
+
+	_logger.debug(f"diffses: {diffses}")
+	# [A, m]
+	norms = numpy.linalg.norm(diffses, axis=2) ** 3
+	_logger.debug(f"norms: {norms}")
+	# [A, m, cart] [A, cart] -> [A, m]
+	ases = (numpy.einsum("...ji, ...i", diffses, ps) / norms) ** 2
+	_logger.debug(f"ases: {ases}")
+
+	# [A, 1], denom [m] + [A, 1] -> [A, m]
+	# [A, m]
+	bses = 2 * ws[:, None] / ((numpy.pi * fs) ** 2 + ws[:, None] ** 2)
+
+	_logger.debug(f"bses: {bses}")
+	# returns shape [A, m]
+	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 = 2 * ws[:, None, :] / ((numpy.pi * fs[:, None]) ** 2 + ws[:, None, :] ** 2)
+	_logger.debug(f"bses: {bses}")
+	return numpy.einsum("...j->...", ases * bses)
+
+
+def fast_efieldxs_for_dipoles(
+	dot_inputs: numpy.ndarray, dipoles: 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.
+	"""
+	# name indexes:
+	# A: dipole config index
+	# j: dipole index within a config
+	# cart: cartesian index
+	# m: measurement index
+
+	# Indexes [A, cart]
+	ps = dipoles[:, 0:3]
+	ss = dipoles[:, 3:6]
+	# Indexes [A]
+	ws = dipoles[:, 6]
+
+	# Indexes [m, cart]
+	rs = dot_inputs[:, 0:3]
+	# Indexes [m]
+	fs = dot_inputs[:, 3]
+
+	# Indexes [m, cart] - [A, 1, cart]
+	# Broadcasting from right
+	# diffses.indexes [A, m, cart]
+	diffses = rs - ss[:, None, :]
+
+	# [A, m, cart][2] = cart
+	# norms takes out axis 2, the last one, giving [A, m]
+	norms = numpy.linalg.norm(diffses, axis=2)
+
+	# long story but this ends up becoming (A, 1, j)
+	# Some evidence is looking at ps term, which has shape (A, 1, j, cart=0) becoming (A, 1, j)
+
+	# [A, m, cart] einsum [A, cart] explicitly gives [A, m]
+	dot_products = numpy.einsum("amc,ac->am", diffses, ps) / (norms**2)
+
+	# [m, A] * [cart, m, A] -> [cart, m, A], transpose that and you get [A, m, cart]
+	projections = numpy.transpose(
+		numpy.transpose(dot_products) * numpy.transpose(diffses)
+	)
+
+	# numerator [A, m, cart] - [A, 1, cart] -> [A, m, cart][:, :, 0] -> [A, m]
+	alphas = (3 * projections - ps[:, numpy.newaxis])[:, :, 0] / norms**3
+
+	# [A, m]
+	ases = alphas**2
+
+	# [A, 1], denom [m] + [A, 1] -> [A, m]
+	# [A, m]
+	bses = 2 * ws[:, None] / ((numpy.pi * fs) ** 2 + ws[:, None] ** 2)
+
+	# return shape [A, m, j]
+	return ases * bses
+
+
+def fast_efieldxs_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.
+	"""
+	# name indexes:
+	# A: dipole config index
+	# j: dipole index within a config
+	# cart: cartesian index
+	# m: measurement index
+
+	# Indexes [A, j, cart]
+	ps = dipoleses[:, :, 0:3]
+	ss = dipoleses[:, :, 3:6]
+	# Indexes [A, j]
+	ws = dipoleses[:, :, 6]
+
+	# Indexes [m, cart]
+	rs = dot_inputs[:, 0:3]
+	# Indexes [m]
+	fs = dot_inputs[:, 3]
+
+	# Indexes [m, 1, cart] - [A, 1, j, cart]
+	# Broadcasting from right
+	# diffses.indexes [A, m, j, cart]
+	diffses = rs[:, None] - ss[:, None, :]
+
+	# norms takes out axis 3, the last one, giving [A, m, j]
+	norms = numpy.linalg.norm(diffses, axis=3)
+
+	# long story but this ends up becoming (A, 1, j)
+	# Some evidence is looking at ps term, which has shape (A, 1, j, cart=0) becoming (A, 1, j)
+	alphas = (
+		(
+			3
+			* numpy.transpose(
+				numpy.transpose(
+					numpy.einsum("abcd,acd->abc", diffses, ps) / (norms**2)
+				)
+				* numpy.transpose(diffses)
+			)[:, :, :, 0]
+		)
+		- ps[:, numpy.newaxis, :, 0]
+	) / (norms**3)
+	ases = alphas**2
+	# bses.shape [A, m, j)]
+	bses = 2 * ws[:, None, :] / ((numpy.pi * fs[:, None]) ** 2 + ws[:, None, :] ** 2)
+
+	# return shape [A, m, j]
+	return numpy.einsum("...j->...", ases * bses)
+
+
+def fast_vs_for_asymmetric_dipoleses(
+	dot_inputs: numpy.ndarray, dipoleses: numpy.ndarray, temp: 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.
+	Within each dipole, the expected format is [px, py, pz, sx, sy, sz, e1, e2, w]
+	The passed in w is expected to be half the actual. This is bad, but here for historical reasons to be changed later.
+	"""
+	raw_ps = dipoleses[:, :, 0:3]
+	ss = dipoleses[:, :, 3:6]
+	e1s = dipoleses[:, :, 6]
+	e2s = dipoleses[:, :, 7]
+	raw_ws = dipoleses[:, :, 8]
+
+	rs = dot_inputs[:, 0:3]
+	fs = dot_inputs[:, 3]
+
+	diffses = rs[:, None] - ss[:, None, :]
+
+	_logger.warning(
+		"This method is very likely to be broken, and should not be used without more thought"
+	)
+	w1s = numpy.exp(-e1s / temp) * raw_ws
+	w2s = numpy.exp(-e2s / temp) * raw_ws
+
+	mag_prefactor = 4 * ((w1s * w2s) / ((w1s + w2s) ** 2))
+	ws = w1s + w2s
+
+	# some annoying broadcast thing here?
+	ps = (raw_ps.T * mag_prefactor.T).T
+
+	norms = numpy.linalg.norm(diffses, axis=3) ** 3
+
+	ases = (numpy.einsum("abcd,acd->abc", diffses, ps) / norms) ** 2
+
+	bses = ws[:, None, :] / ((numpy.pi * fs[:, None]) ** 2 + ws[:, None, :] ** 2)
+
+	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.
+	"""
+	return numpy.all(numpy.logical_and(low < a, high > a), axis=1)

poetry.lock

@@ -1,24 +1,65 @@
 [[package]]
-name = "atomicwrites"
-version = "1.4.0"
-description = "Atomic file writes."
+name = "appnope"
+version = "0.1.2"
+description = "Disable App Nap on macOS >= 10.9"
 category = "dev"
 optional = false
-python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*"
+python-versions = "*"
 
 [[package]]
-name = "attrs"
-version = "21.4.0"
-description = "Classes Without Boilerplate"
+name = "asttokens"
+version = "2.0.5"
+description = "Annotate AST trees with source code positions"
 category = "dev"
 optional = false
-python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.*"
+python-versions = "*"
+
+[package.dependencies]
+six = "*"
 
 [package.extras]
-dev = ["coverage[toml] (>=5.0.2)", "hypothesis", "pympler", "pytest (>=4.3.0)", "six", "mypy", "pytest-mypy-plugins", "zope.interface", "furo", "sphinx", "sphinx-notfound-page", "pre-commit", "cloudpickle"]
-docs = ["furo", "sphinx", "zope.interface", "sphinx-notfound-page"]
-tests = ["coverage[toml] (>=5.0.2)", "hypothesis", "pympler", "pytest (>=4.3.0)", "six", "mypy", "pytest-mypy-plugins", "zope.interface", "cloudpickle"]
-tests_no_zope = ["coverage[toml] (>=5.0.2)", "hypothesis", "pympler", "pytest (>=4.3.0)", "six", "mypy", "pytest-mypy-plugins", "cloudpickle"]
+test = ["astroid", "pytest"]
+
+[[package]]
+name = "backcall"
+version = "0.2.0"
+description = "Specifications for callback functions passed in to an API"
+category = "dev"
+optional = false
+python-versions = "*"
+
+[[package]]
+name = "black"
+version = "22.3.0"
+description = "The uncompromising code formatter."
+category = "dev"
+optional = false
+python-versions = ">=3.6.2"
+
+[package.dependencies]
+click = ">=8.0.0"
+mypy-extensions = ">=0.4.3"
+pathspec = ">=0.9.0"
+platformdirs = ">=2"
+tomli = {version = ">=1.1.0", markers = "python_version < \"3.11\""}
+typing-extensions = {version = ">=3.10.0.0", markers = "python_version < \"3.10\""}
+
+[package.extras]
+colorama = ["colorama (>=0.4.3)"]
+d = ["aiohttp (>=3.7.4)"]
+jupyter = ["ipython (>=7.8.0)", "tokenize-rt (>=3.2.0)"]
+uvloop = ["uvloop (>=0.15.2)"]
+
+[[package]]
+name = "click"
+version = "8.1.2"
+description = "Composable command line interface toolkit"
+category = "dev"
+optional = false
+python-versions = ">=3.7"
+
+[package.dependencies]
+colorama = {version = "*", markers = "platform_system == \"Windows\""}
 
 [[package]]
 name = "colorama"
@@ -28,29 +69,67 @@ category = "dev"
 optional = false
 python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.*"
 
+[[package]]
+name = "colored"
+version = "1.4.4"
+description = "Simple library for color and formatting to terminal"
+category = "dev"
+optional = false
+python-versions = "*"
+
 [[package]]
 name = "coverage"
-version = "6.2"
+version = "7.2.7"
 description = "Code coverage measurement for Python"
 category = "dev"
 optional = false
-python-versions = ">=3.6"
+python-versions = ">=3.7"
+
+[package.dependencies]
+tomli = {version = "*", optional = true, markers = "python_full_version <= \"3.11.0a6\" and extra == \"toml\""}
 
 [package.extras]
 toml = ["tomli"]
 
+[[package]]
+name = "decorator"
+version = "5.1.1"
+description = "Decorators for Humans"
+category = "dev"
+optional = false
+python-versions = ">=3.5"
+
+[[package]]
+name = "exceptiongroup"
+version = "1.1.2"
+description = "Backport of PEP 654 (exception groups)"
+category = "dev"
+optional = false
+python-versions = ">=3.7"
+
+[package.extras]
+test = ["pytest (>=6)"]
+
+[[package]]
+name = "executing"
+version = "0.8.3"
+description = "Get the currently executing AST node of a frame, and other information"
+category = "dev"
+optional = false
+python-versions = "*"
+
 [[package]]
 name = "flake8"
-version = "3.9.2"
+version = "4.0.1"
 description = "the modular source code checker: pep8 pyflakes and co"
 category = "dev"
 optional = false
-python-versions = "!=3.0.*,!=3.1.*,!=3.2.*,!=3.3.*,!=3.4.*,>=2.7"
+python-versions = ">=3.6"
 
 [package.dependencies]
 mccabe = ">=0.6.0,<0.7.0"
-pycodestyle = ">=2.7.0,<2.8.0"
-pyflakes = ">=2.3.0,<2.4.0"
+pycodestyle = ">=2.8.0,<2.9.0"
+pyflakes = ">=2.4.0,<2.5.0"
 
 [[package]]
 name = "iniconfig"
@@ -60,6 +139,67 @@ category = "dev"
 optional = false
 python-versions = "*"
 
+[[package]]
+name = "ipython"
+version = "8.2.0"
+description = "IPython: Productive Interactive Computing"
+category = "dev"
+optional = false
+python-versions = ">=3.8"
+
+[package.dependencies]
+appnope = {version = "*", markers = "sys_platform == \"darwin\""}
+backcall = "*"
+colorama = {version = "*", markers = "sys_platform == \"win32\""}
+decorator = "*"
+jedi = ">=0.16"
+matplotlib-inline = "*"
+pexpect = {version = ">4.3", markers = "sys_platform != \"win32\""}
+pickleshare = "*"
+prompt-toolkit = ">=2.0.0,<3.0.0 || >3.0.0,<3.0.1 || >3.0.1,<3.1.0"
+pygments = ">=2.4.0"
+stack-data = "*"
+traitlets = ">=5"
+
+[package.extras]
+all = ["black", "Sphinx (>=1.3)", "ipykernel", "nbconvert", "nbformat", "ipywidgets", "notebook", "ipyparallel", "qtconsole", "pytest (<7.1)", "pytest-asyncio", "testpath", "curio", "matplotlib (!=3.2.0)", "numpy (>=1.19)", "pandas", "trio"]
+black = ["black"]
+doc = ["Sphinx (>=1.3)"]
+kernel = ["ipykernel"]
+nbconvert = ["nbconvert"]
+nbformat = ["nbformat"]
+notebook = ["ipywidgets", "notebook"]
+parallel = ["ipyparallel"]
+qtconsole = ["qtconsole"]
+test = ["pytest (<7.1)", "pytest-asyncio", "testpath"]
+test_extra = ["pytest (<7.1)", "pytest-asyncio", "testpath", "curio", "matplotlib (!=3.2.0)", "nbformat", "numpy (>=1.19)", "pandas", "trio"]
+
+[[package]]
+name = "jedi"
+version = "0.18.1"
+description = "An autocompletion tool for Python that can be used for text editors."
+category = "dev"
+optional = false
+python-versions = ">=3.6"
+
+[package.dependencies]
+parso = ">=0.8.0,<0.9.0"
+
+[package.extras]
+qa = ["flake8 (==3.8.3)", "mypy (==0.782)"]
+testing = ["Django (<3.1)", "colorama", "docopt", "pytest (<7.0.0)"]
+
+[[package]]
+name = "matplotlib-inline"
+version = "0.1.3"
+description = "Inline Matplotlib backend for Jupyter"
+category = "dev"
+optional = false
+python-versions = ">=3.5"
+
+[package.dependencies]
+traitlets = "*"
+
 [[package]]
 name = "mccabe"
 version = "0.6.1"
@@ -70,20 +210,21 @@ python-versions = "*"
 
 [[package]]
 name = "mypy"
-version = "0.910"
+version = "0.961"
 description = "Optional static typing for Python"
 category = "dev"
 optional = false
-python-versions = ">=3.5"
+python-versions = ">=3.6"
 
 [package.dependencies]
-mypy-extensions = ">=0.4.3,<0.5.0"
-toml = "*"
-typing-extensions = ">=3.7.4"
+mypy-extensions = ">=0.4.3"
+tomli = {version = ">=1.1.0", markers = "python_version < \"3.11\""}
+typing-extensions = ">=3.10"
 
 [package.extras]
 dmypy = ["psutil (>=4.0)"]
-python2 = ["typed-ast (>=1.4.0,<1.5.0)"]
+python2 = ["typed-ast (>=1.4.0,<2)"]
+reports = ["lxml"]
 
 [[package]]
 name = "mypy-extensions"
@@ -95,7 +236,7 @@ python-versions = "*"
 
 [[package]]
 name = "numpy"
-version = "1.22.0"
+version = "1.22.3"
 description = "NumPy is the fundamental package for array computing with Python."
 category = "main"
 optional = false
@@ -112,6 +253,57 @@ python-versions = ">=3.6"
 [package.dependencies]
 pyparsing = ">=2.0.2,<3.0.5 || >3.0.5"
 
+[[package]]
+name = "parso"
+version = "0.8.3"
+description = "A Python Parser"
+category = "dev"
+optional = false
+python-versions = ">=3.6"
+
+[package.extras]
+qa = ["flake8 (==3.8.3)", "mypy (==0.782)"]
+testing = ["docopt", "pytest (<6.0.0)"]
+
+[[package]]
+name = "pathspec"
+version = "0.9.0"
+description = "Utility library for gitignore style pattern matching of file paths."
+category = "dev"
+optional = false
+python-versions = "!=3.0.*,!=3.1.*,!=3.2.*,!=3.3.*,!=3.4.*,>=2.7"
+
+[[package]]
+name = "pexpect"
+version = "4.8.0"
+description = "Pexpect allows easy control of interactive console applications."
+category = "dev"
+optional = false
+python-versions = "*"
+
+[package.dependencies]
+ptyprocess = ">=0.5"
+
+[[package]]
+name = "pickleshare"
+version = "0.7.5"
+description = "Tiny 'shelve'-like database with concurrency support"
+category = "dev"
+optional = false
+python-versions = "*"
+
+[[package]]
+name = "platformdirs"
+version = "2.5.1"
+description = "A small Python module for determining appropriate platform-specific dirs, e.g. a \"user data dir\"."
+category = "dev"
+optional = false
+python-versions = ">=3.7"
+
+[package.extras]
+docs = ["Sphinx (>=4)", "furo (>=2021.7.5b38)", "proselint (>=0.10.2)", "sphinx-autodoc-typehints (>=1.12)"]
+test = ["appdirs (==1.4.4)", "pytest (>=6)", "pytest-cov (>=2.7)", "pytest-mock (>=3.6)"]
+
 [[package]]
 name = "pluggy"
 version = "1.0.0"
@@ -125,28 +317,58 @@ dev = ["pre-commit", "tox"]
 testing = ["pytest", "pytest-benchmark"]
 
 [[package]]
-name = "py"
-version = "1.11.0"
-description = "library with cross-python path, ini-parsing, io, code, log facilities"
+name = "prompt-toolkit"
+version = "3.0.28"
+description = "Library for building powerful interactive command lines in Python"
+category = "dev"
+optional = false
+python-versions = ">=3.6.2"
+
+[package.dependencies]
+wcwidth = "*"
+
+[[package]]
+name = "ptyprocess"
+version = "0.7.0"
+description = "Run a subprocess in a pseudo terminal"
+category = "dev"
+optional = false
+python-versions = "*"
+
+[[package]]
+name = "pure-eval"
+version = "0.2.2"
+description = "Safely evaluate AST nodes without side effects"
+category = "dev"
+optional = false
+python-versions = "*"
+
+[package.extras]
+tests = ["pytest"]
+
+[[package]]
+name = "pycodestyle"
+version = "2.8.0"
+description = "Python style guide checker"
 category = "dev"
 optional = false
 python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.*"
 
 [[package]]
-name = "pycodestyle"
-version = "2.7.0"
-description = "Python style guide checker"
+name = "pyflakes"
+version = "2.4.0"
+description = "passive checker of Python programs"
 category = "dev"
 optional = false
 python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*"
 
 [[package]]
-name = "pyflakes"
-version = "2.3.1"
-description = "passive checker of Python programs"
+name = "pygments"
+version = "2.11.2"
+description = "Pygments is a syntax highlighting package written in Python."
 category = "dev"
 optional = false
-python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*"
+python-versions = ">=3.5"
 
 [[package]]
 name = "pyparsing"
@@ -161,48 +383,108 @@ diagrams = ["jinja2", "railroad-diagrams"]
 
 [[package]]
 name = "pytest"
-version = "6.2.5"
+version = "7.4.0"
 description = "pytest: simple powerful testing with Python"
 category = "dev"
 optional = false
-python-versions = ">=3.6"
+python-versions = ">=3.7"
 
 [package.dependencies]
-atomicwrites = {version = ">=1.0", markers = "sys_platform == \"win32\""}
-attrs = ">=19.2.0"
 colorama = {version = "*", markers = "sys_platform == \"win32\""}
+exceptiongroup = {version = ">=1.0.0rc8", markers = "python_version < \"3.11\""}
 iniconfig = "*"
 packaging = "*"
 pluggy = ">=0.12,<2.0"
-py = ">=1.8.2"
-toml = "*"
+tomli = {version = ">=1.0.0", markers = "python_version < \"3.11\""}
 
 [package.extras]
-testing = ["argcomplete", "hypothesis (>=3.56)", "mock", "nose", "requests", "xmlschema"]
+testing = ["argcomplete", "attrs (>=19.2.0)", "hypothesis (>=3.56)", "mock", "nose", "pygments (>=2.7.2)", "requests", "setuptools", "xmlschema"]
 
 [[package]]
 name = "pytest-cov"
-version = "2.12.1"
+version = "4.1.0"
 description = "Pytest plugin for measuring coverage."
 category = "dev"
 optional = false
-python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.*"
+python-versions = ">=3.7"
 
 [package.dependencies]
-coverage = ">=5.2.1"
+coverage = {version = ">=5.2.1", extras = ["toml"]}
 pytest = ">=4.6"
-toml = "*"
 
 [package.extras]
 testing = ["fields", "hunter", "process-tests", "six", "pytest-xdist", "virtualenv"]
 
 [[package]]
-name = "toml"
-version = "0.10.2"
-description = "Python Library for Tom's Obvious, Minimal Language"
+name = "scipy"
+version = "1.10.0"
+description = "Fundamental algorithms for scientific computing in Python"
+category = "main"
+optional = false
+python-versions = "<3.12,>=3.8"
+
+[package.dependencies]
+numpy = ">=1.19.5,<1.27.0"
+
+[package.extras]
+test = ["pytest", "pytest-cov", "pytest-timeout", "pytest-xdist", "asv", "mpmath", "gmpy2", "threadpoolctl", "scikit-umfpack", "pooch"]
+doc = ["sphinx (!=4.1.0)", "pydata-sphinx-theme (==0.9.0)", "sphinx-design (>=0.2.0)", "matplotlib (>2)", "numpydoc"]
+dev = ["mypy", "typing-extensions", "pycodestyle", "flake8", "rich-click", "click", "doit (>=0.36.0)", "pydevtool"]
+
+[[package]]
+name = "six"
+version = "1.16.0"
+description = "Python 2 and 3 compatibility utilities"
 category = "dev"
 optional = false
-python-versions = ">=2.6, !=3.0.*, !=3.1.*, !=3.2.*"
+python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*"
+
+[[package]]
+name = "stack-data"
+version = "0.2.0"
+description = "Extract data from python stack frames and tracebacks for informative displays"
+category = "dev"
+optional = false
+python-versions = "*"
+
+[package.dependencies]
+asttokens = "*"
+executing = "*"
+pure-eval = "*"
+
+[package.extras]
+tests = ["pytest", "typeguard", "pygments", "littleutils", "cython"]
+
+[[package]]
+name = "syrupy"
+version = "4.0.8"
+description = "Pytest Snapshot Test Utility"
+category = "dev"
+optional = false
+python-versions = ">=3.8.1,<4"
+
+[package.dependencies]
+colored = ">=1.3.92,<2.0.0"
+pytest = ">=7.0.0,<8.0.0"
+
+[[package]]
+name = "tomli"
+version = "2.0.1"
+description = "A lil' TOML parser"
+category = "dev"
+optional = false
+python-versions = ">=3.7"
+
+[[package]]
+name = "traitlets"
+version = "5.1.1"
+description = "Traitlets Python configuration system"
+category = "dev"
+optional = false
+python-versions = ">=3.7"
+
+[package.extras]
+test = ["pytest"]
 
 [[package]]
 name = "typing-extensions"
@@ -212,175 +494,61 @@ category = "dev"
 optional = false
 python-versions = ">=3.6"
 
+[[package]]
+name = "wcwidth"
+version = "0.2.5"
+description = "Measures the displayed width of unicode strings in a terminal"
+category = "dev"
+optional = false
+python-versions = "*"
+
 [metadata]
 lock-version = "1.1"
-python-versions = "^3.8,<3.10"
-content-hash = "d9790d575f5c79aca211a7f2a4b3539f35ea13d0ed72d7fb796e2f36126738eb"
+python-versions = ">=3.8.1,<3.10"
+content-hash = "b5275c33449e8f85acbf9c0f6dfe1ec4e7296adfa16360d782b33534e1223638"
 
 [metadata.files]
-atomicwrites = [
-    {file = "atomicwrites-1.4.0-py2.py3-none-any.whl", hash = "sha256:6d1784dea7c0c8d4a5172b6c620f40b6e4cbfdf96d783691f2e1302a7b88e197"},
-    {file = "atomicwrites-1.4.0.tar.gz", hash = "sha256:ae70396ad1a434f9c7046fd2dd196fc04b12f9e91ffb859164193be8b6168a7a"},
-]
-attrs = [
-    {file = "attrs-21.4.0-py2.py3-none-any.whl", hash = "sha256:2d27e3784d7a565d36ab851fe94887c5eccd6a463168875832a1be79c82828b4"},
-    {file = "attrs-21.4.0.tar.gz", hash = "sha256:626ba8234211db98e869df76230a137c4c40a12d72445c45d5f5b716f076e2fd"},
-]
-colorama = [
-    {file = "colorama-0.4.4-py2.py3-none-any.whl", hash = "sha256:9f47eda37229f68eee03b24b9748937c7dc3868f906e8ba69fbcbdd3bc5dc3e2"},
-    {file = "colorama-0.4.4.tar.gz", hash = "sha256:5941b2b48a20143d2267e95b1c2a7603ce057ee39fd88e7329b0c292aa16869b"},
-]
-coverage = [
-    {file = "coverage-6.2-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:6dbc1536e105adda7a6312c778f15aaabe583b0e9a0b0a324990334fd458c94b"},
-    {file = "coverage-6.2-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:174cf9b4bef0db2e8244f82059a5a72bd47e1d40e71c68ab055425172b16b7d0"},
-    {file = "coverage-6.2-cp310-cp310-manylinux_2_5_i686.manylinux1_i686.manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:92b8c845527eae547a2a6617d336adc56394050c3ed8a6918683646328fbb6da"},
-    {file = "coverage-6.2-cp310-cp310-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:c7912d1526299cb04c88288e148c6c87c0df600eca76efd99d84396cfe00ef1d"},
-    {file = "coverage-6.2-cp310-cp310-musllinux_1_1_aarch64.whl", hash = "sha256:d5d2033d5db1d58ae2d62f095e1aefb6988af65b4b12cb8987af409587cc0739"},
-    {file = "coverage-6.2-cp310-cp310-musllinux_1_1_i686.whl", hash = "sha256:3feac4084291642165c3a0d9eaebedf19ffa505016c4d3db15bfe235718d4971"},
-    {file = "coverage-6.2-cp310-cp310-musllinux_1_1_x86_64.whl", hash = "sha256:276651978c94a8c5672ea60a2656e95a3cce2a3f31e9fb2d5ebd4c215d095840"},
-    {file = "coverage-6.2-cp310-cp310-win32.whl", hash = "sha256:f506af4f27def639ba45789fa6fde45f9a217da0be05f8910458e4557eed020c"},
-    {file = "coverage-6.2-cp310-cp310-win_amd64.whl", hash = "sha256:3f7c17209eef285c86f819ff04a6d4cbee9b33ef05cbcaae4c0b4e8e06b3ec8f"},
-    {file = "coverage-6.2-cp311-cp311-macosx_10_14_x86_64.whl", hash = "sha256:13362889b2d46e8d9f97c421539c97c963e34031ab0cb89e8ca83a10cc71ac76"},
-    {file = "coverage-6.2-cp311-cp311-manylinux_2_5_x86_64.manylinux1_x86_64.whl", hash = "sha256:22e60a3ca5acba37d1d4a2ee66e051f5b0e1b9ac950b5b0cf4aa5366eda41d47"},
-    {file = "coverage-6.2-cp311-cp311-win_amd64.whl", hash = "sha256:b637c57fdb8be84e91fac60d9325a66a5981f8086c954ea2772efe28425eaf64"},
-    {file = "coverage-6.2-cp36-cp36m-macosx_10_9_x86_64.whl", hash = "sha256:f467bbb837691ab5a8ca359199d3429a11a01e6dfb3d9dcc676dc035ca93c0a9"},
-    {file = "coverage-6.2-cp36-cp36m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:2641f803ee9f95b1f387f3e8f3bf28d83d9b69a39e9911e5bfee832bea75240d"},
-    {file = "coverage-6.2-cp36-cp36m-manylinux_2_5_i686.manylinux1_i686.manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:1219d760ccfafc03c0822ae2e06e3b1248a8e6d1a70928966bafc6838d3c9e48"},
-    {file = "coverage-6.2-cp36-cp36m-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:9a2b5b52be0a8626fcbffd7e689781bf8c2ac01613e77feda93d96184949a98e"},
-    {file = "coverage-6.2-cp36-cp36m-musllinux_1_1_aarch64.whl", hash = "sha256:8e2c35a4c1f269704e90888e56f794e2d9c0262fb0c1b1c8c4ee44d9b9e77b5d"},
-    {file = "coverage-6.2-cp36-cp36m-musllinux_1_1_i686.whl", hash = "sha256:5d6b09c972ce9200264c35a1d53d43ca55ef61836d9ec60f0d44273a31aa9f17"},
-    {file = "coverage-6.2-cp36-cp36m-musllinux_1_1_x86_64.whl", hash = "sha256:e3db840a4dee542e37e09f30859f1612da90e1c5239a6a2498c473183a50e781"},
-    {file = "coverage-6.2-cp36-cp36m-win32.whl", hash = "sha256:4e547122ca2d244f7c090fe3f4b5a5861255ff66b7ab6d98f44a0222aaf8671a"},
-    {file = "coverage-6.2-cp36-cp36m-win_amd64.whl", hash = "sha256:01774a2c2c729619760320270e42cd9e797427ecfddd32c2a7b639cdc481f3c0"},
-    {file = "coverage-6.2-cp37-cp37m-macosx_10_9_x86_64.whl", hash = "sha256:fb8b8ee99b3fffe4fd86f4c81b35a6bf7e4462cba019997af2fe679365db0c49"},
-    {file = "coverage-6.2-cp37-cp37m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:619346d57c7126ae49ac95b11b0dc8e36c1dd49d148477461bb66c8cf13bb521"},
-    {file = "coverage-6.2-cp37-cp37m-manylinux_2_5_i686.manylinux1_i686.manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:0a7726f74ff63f41e95ed3a89fef002916c828bb5fcae83b505b49d81a066884"},
-    {file = "coverage-6.2-cp37-cp37m-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:cfd9386c1d6f13b37e05a91a8583e802f8059bebfccde61a418c5808dea6bbfa"},
-    {file = "coverage-6.2-cp37-cp37m-musllinux_1_1_aarch64.whl", hash = "sha256:17e6c11038d4ed6e8af1407d9e89a2904d573be29d51515f14262d7f10ef0a64"},
-    {file = "coverage-6.2-cp37-cp37m-musllinux_1_1_i686.whl", hash = "sha256:c254b03032d5a06de049ce8bca8338a5185f07fb76600afff3c161e053d88617"},
-    {file = "coverage-6.2-cp37-cp37m-musllinux_1_1_x86_64.whl", hash = "sha256:dca38a21e4423f3edb821292e97cec7ad38086f84313462098568baedf4331f8"},
-    {file = "coverage-6.2-cp37-cp37m-win32.whl", hash = "sha256:600617008aa82032ddeace2535626d1bc212dfff32b43989539deda63b3f36e4"},
-    {file = "coverage-6.2-cp37-cp37m-win_amd64.whl", hash = "sha256:bf154ba7ee2fd613eb541c2bc03d3d9ac667080a737449d1a3fb342740eb1a74"},
-    {file = "coverage-6.2-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:f9afb5b746781fc2abce26193d1c817b7eb0e11459510fba65d2bd77fe161d9e"},
-    {file = "coverage-6.2-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:edcada2e24ed68f019175c2b2af2a8b481d3d084798b8c20d15d34f5c733fa58"},
-    {file = "coverage-6.2-cp38-cp38-manylinux_2_5_i686.manylinux1_i686.manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:a9c8c4283e17690ff1a7427123ffb428ad6a52ed720d550e299e8291e33184dc"},
-    {file = "coverage-6.2-cp38-cp38-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:f614fc9956d76d8a88a88bb41ddc12709caa755666f580af3a688899721efecd"},
-    {file = "coverage-6.2-cp38-cp38-musllinux_1_1_aarch64.whl", hash = "sha256:9365ed5cce5d0cf2c10afc6add145c5037d3148585b8ae0e77cc1efdd6aa2953"},
-    {file = "coverage-6.2-cp38-cp38-musllinux_1_1_i686.whl", hash = "sha256:8bdfe9ff3a4ea37d17f172ac0dff1e1c383aec17a636b9b35906babc9f0f5475"},
-    {file = "coverage-6.2-cp38-cp38-musllinux_1_1_x86_64.whl", hash = "sha256:63c424e6f5b4ab1cf1e23a43b12f542b0ec2e54f99ec9f11b75382152981df57"},
-    {file = "coverage-6.2-cp38-cp38-win32.whl", hash = "sha256:49dbff64961bc9bdd2289a2bda6a3a5a331964ba5497f694e2cbd540d656dc1c"},
-    {file = "coverage-6.2-cp38-cp38-win_amd64.whl", hash = "sha256:9a29311bd6429be317c1f3fe4bc06c4c5ee45e2fa61b2a19d4d1d6111cb94af2"},
-    {file = "coverage-6.2-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:03b20e52b7d31be571c9c06b74746746d4eb82fc260e594dc662ed48145e9efd"},
-    {file = "coverage-6.2-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:215f8afcc02a24c2d9a10d3790b21054b58d71f4b3c6f055d4bb1b15cecce685"},
-    {file = "coverage-6.2-cp39-cp39-manylinux_2_5_i686.manylinux1_i686.manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:a4bdeb0a52d1d04123b41d90a4390b096f3ef38eee35e11f0b22c2d031222c6c"},
-    {file = "coverage-6.2-cp39-cp39-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:c332d8f8d448ded473b97fefe4a0983265af21917d8b0cdcb8bb06b2afe632c3"},
-    {file = "coverage-6.2-cp39-cp39-win32.whl", hash = "sha256:6e1394d24d5938e561fbeaa0cd3d356207579c28bd1792f25a068743f2d5b282"},
-    {file = "coverage-6.2-cp39-cp39-win_amd64.whl", hash = "sha256:86f2e78b1eff847609b1ca8050c9e1fa3bd44ce755b2ec30e70f2d3ba3844644"},
-    {file = "coverage-6.2-pp36.pp37.pp38-none-any.whl", hash = "sha256:5829192582c0ec8ca4a2532407bc14c2f338d9878a10442f5d03804a95fac9de"},
-    {file = "coverage-6.2.tar.gz", hash = "sha256:e2cad8093172b7d1595b4ad66f24270808658e11acf43a8f95b41276162eb5b8"},
-]
-flake8 = [
-    {file = "flake8-3.9.2-py2.py3-none-any.whl", hash = "sha256:bf8fd333346d844f616e8d47905ef3a3384edae6b4e9beb0c5101e25e3110907"},
-    {file = "flake8-3.9.2.tar.gz", hash = "sha256:07528381786f2a6237b061f6e96610a4167b226cb926e2aa2b6b1d78057c576b"},
-]
-iniconfig = [
-    {file = "iniconfig-1.1.1-py2.py3-none-any.whl", hash = "sha256:011e24c64b7f47f6ebd835bb12a743f2fbe9a26d4cecaa7f53bc4f35ee9da8b3"},
-    {file = "iniconfig-1.1.1.tar.gz", hash = "sha256:bc3af051d7d14b2ee5ef9969666def0cd1a000e121eaea580d4a313df4b37f32"},
-]
-mccabe = [
-    {file = "mccabe-0.6.1-py2.py3-none-any.whl", hash = "sha256:ab8a6258860da4b6677da4bd2fe5dc2c659cff31b3ee4f7f5d64e79735b80d42"},
-    {file = "mccabe-0.6.1.tar.gz", hash = "sha256:dd8d182285a0fe56bace7f45b5e7d1a6ebcbf524e8f3bd87eb0f125271b8831f"},
-]
-mypy = [
-    {file = "mypy-0.910-cp35-cp35m-macosx_10_9_x86_64.whl", hash = "sha256:a155d80ea6cee511a3694b108c4494a39f42de11ee4e61e72bc424c490e46457"},
-    {file = "mypy-0.910-cp35-cp35m-manylinux1_x86_64.whl", hash = "sha256:b94e4b785e304a04ea0828759172a15add27088520dc7e49ceade7834275bedb"},
-    {file = "mypy-0.910-cp35-cp35m-manylinux2010_x86_64.whl", hash = "sha256:088cd9c7904b4ad80bec811053272986611b84221835e079be5bcad029e79dd9"},
-    {file = "mypy-0.910-cp35-cp35m-win_amd64.whl", hash = "sha256:adaeee09bfde366d2c13fe6093a7df5df83c9a2ba98638c7d76b010694db760e"},
-    {file = "mypy-0.910-cp36-cp36m-macosx_10_9_x86_64.whl", hash = "sha256:ecd2c3fe726758037234c93df7e98deb257fd15c24c9180dacf1ef829da5f921"},
-    {file = "mypy-0.910-cp36-cp36m-manylinux1_x86_64.whl", hash = "sha256:d9dd839eb0dc1bbe866a288ba3c1afc33a202015d2ad83b31e875b5905a079b6"},
-    {file = "mypy-0.910-cp36-cp36m-manylinux2010_x86_64.whl", hash = "sha256:3e382b29f8e0ccf19a2df2b29a167591245df90c0b5a2542249873b5c1d78212"},
-    {file = "mypy-0.910-cp36-cp36m-win_amd64.whl", hash = "sha256:53fd2eb27a8ee2892614370896956af2ff61254c275aaee4c230ae771cadd885"},
-    {file = "mypy-0.910-cp37-cp37m-macosx_10_9_x86_64.whl", hash = "sha256:b6fb13123aeef4a3abbcfd7e71773ff3ff1526a7d3dc538f3929a49b42be03f0"},
-    {file = "mypy-0.910-cp37-cp37m-manylinux1_x86_64.whl", hash = "sha256:e4dab234478e3bd3ce83bac4193b2ecd9cf94e720ddd95ce69840273bf44f6de"},
-    {file = "mypy-0.910-cp37-cp37m-manylinux2010_x86_64.whl", hash = "sha256:7df1ead20c81371ccd6091fa3e2878559b5c4d4caadaf1a484cf88d93ca06703"},
-    {file = "mypy-0.910-cp37-cp37m-win_amd64.whl", hash = "sha256:0aadfb2d3935988ec3815952e44058a3100499f5be5b28c34ac9d79f002a4a9a"},
-    {file = "mypy-0.910-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:ec4e0cd079db280b6bdabdc807047ff3e199f334050db5cbb91ba3e959a67504"},
-    {file = "mypy-0.910-cp38-cp38-manylinux1_x86_64.whl", hash = "sha256:119bed3832d961f3a880787bf621634ba042cb8dc850a7429f643508eeac97b9"},
-    {file = "mypy-0.910-cp38-cp38-manylinux2010_x86_64.whl", hash = "sha256:866c41f28cee548475f146aa4d39a51cf3b6a84246969f3759cb3e9c742fc072"},
-    {file = "mypy-0.910-cp38-cp38-win_amd64.whl", hash = "sha256:ceb6e0a6e27fb364fb3853389607cf7eb3a126ad335790fa1e14ed02fba50811"},
-    {file = "mypy-0.910-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:1a85e280d4d217150ce8cb1a6dddffd14e753a4e0c3cf90baabb32cefa41b59e"},
-    {file = "mypy-0.910-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:42c266ced41b65ed40a282c575705325fa7991af370036d3f134518336636f5b"},
-    {file = "mypy-0.910-cp39-cp39-manylinux1_x86_64.whl", hash = "sha256:3c4b8ca36877fc75339253721f69603a9c7fdb5d4d5a95a1a1b899d8b86a4de2"},
-    {file = "mypy-0.910-cp39-cp39-manylinux2010_x86_64.whl", hash = "sha256:c0df2d30ed496a08de5daed2a9ea807d07c21ae0ab23acf541ab88c24b26ab97"},
-    {file = "mypy-0.910-cp39-cp39-win_amd64.whl", hash = "sha256:c6c2602dffb74867498f86e6129fd52a2770c48b7cd3ece77ada4fa38f94eba8"},
-    {file = "mypy-0.910-py3-none-any.whl", hash = "sha256:ef565033fa5a958e62796867b1df10c40263ea9ded87164d67572834e57a174d"},
-    {file = "mypy-0.910.tar.gz", hash = "sha256:704098302473cb31a218f1775a873b376b30b4c18229421e9e9dc8916fd16150"},
-]
-mypy-extensions = [
-    {file = "mypy_extensions-0.4.3-py2.py3-none-any.whl", hash = "sha256:090fedd75945a69ae91ce1303b5824f428daf5a028d2f6ab8a299250a846f15d"},
-    {file = "mypy_extensions-0.4.3.tar.gz", hash = "sha256:2d82818f5bb3e369420cb3c4060a7970edba416647068eb4c5343488a6c604a8"},
-]
-numpy = [
-    {file = "numpy-1.22.0-cp310-cp310-macosx_10_9_universal2.whl", hash = "sha256:3d22662b4b10112c545c91a0741f2436f8ca979ab3d69d03d19322aa970f9695"},
-    {file = "numpy-1.22.0-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:11a1f3816ea82eed4178102c56281782690ab5993251fdfd75039aad4d20385f"},
-    {file = "numpy-1.22.0-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:5dc65644f75a4c2970f21394ad8bea1a844104f0fe01f278631be1c7eae27226"},
-    {file = "numpy-1.22.0-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:42c16cec1c8cf2728f1d539bd55aaa9d6bb48a7de2f41eb944697293ef65a559"},
-    {file = "numpy-1.22.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:a97e82c39d9856fe7d4f9b86d8a1e66eff99cf3a8b7ba48202f659703d27c46f"},
-    {file = "numpy-1.22.0-cp310-cp310-win_amd64.whl", hash = "sha256:e41e8951749c4b5c9a2dc5fdbc1a4eec6ab2a140fdae9b460b0f557eed870f4d"},
-    {file = "numpy-1.22.0-cp38-cp38-macosx_10_9_universal2.whl", hash = "sha256:bece0a4a49e60e472a6d1f70ac6cdea00f9ab80ff01132f96bd970cdd8a9e5a9"},
-    {file = "numpy-1.22.0-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:818b9be7900e8dc23e013a92779135623476f44a0de58b40c32a15368c01d471"},
-    {file = "numpy-1.22.0-cp38-cp38-macosx_11_0_arm64.whl", hash = "sha256:47ee7a839f5885bc0c63a74aabb91f6f40d7d7b639253768c4199b37aede7982"},
-    {file = "numpy-1.22.0-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:a024181d7aef0004d76fb3bce2a4c9f2e67a609a9e2a6ff2571d30e9976aa383"},
-    {file = "numpy-1.22.0-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:f71d57cc8645f14816ae249407d309be250ad8de93ef61d9709b45a0ddf4050c"},
-    {file = "numpy-1.22.0-cp38-cp38-win32.whl", hash = "sha256:283d9de87c0133ef98f93dfc09fad3fb382f2a15580de75c02b5bb36a5a159a5"},
-    {file = "numpy-1.22.0-cp38-cp38-win_amd64.whl", hash = "sha256:2762331de395739c91f1abb88041f94a080cb1143aeec791b3b223976228af3f"},
-    {file = "numpy-1.22.0-cp39-cp39-macosx_10_9_universal2.whl", hash = "sha256:76ba7c40e80f9dc815c5e896330700fd6e20814e69da9c1267d65a4d051080f1"},
-    {file = "numpy-1.22.0-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:0cfe07133fd00b27edee5e6385e333e9eeb010607e8a46e1cd673f05f8596595"},
-    {file = "numpy-1.22.0-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:6ed0d073a9c54ac40c41a9c2d53fcc3d4d4ed607670b9e7b0de1ba13b4cbfe6f"},
-    {file = "numpy-1.22.0-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:41388e32e40b41dd56eb37fcaa7488b2b47b0adf77c66154d6b89622c110dfe9"},
-    {file = "numpy-1.22.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:b55b953a1bdb465f4dc181758570d321db4ac23005f90ffd2b434cc6609a63dd"},
-    {file = "numpy-1.22.0-cp39-cp39-win32.whl", hash = "sha256:5a311ee4d983c487a0ab546708edbdd759393a3dc9cd30305170149fedd23c88"},
-    {file = "numpy-1.22.0-cp39-cp39-win_amd64.whl", hash = "sha256:a97a954a8c2f046d3817c2bce16e3c7e9a9c2afffaf0400f5c16df5172a67c9c"},
-    {file = "numpy-1.22.0-pp38-pypy38_pp73-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:bb02929b0d6bfab4c48a79bd805bd7419114606947ec8284476167415171f55b"},
-    {file = "numpy-1.22.0.zip", hash = "sha256:a955e4128ac36797aaffd49ab44ec74a71c11d6938df83b1285492d277db5397"},
-]
-packaging = [
-    {file = "packaging-21.3-py3-none-any.whl", hash = "sha256:ef103e05f519cdc783ae24ea4e2e0f508a9c99b2d4969652eed6a2e1ea5bd522"},
-    {file = "packaging-21.3.tar.gz", hash = "sha256:dd47c42927d89ab911e606518907cc2d3a1f38bbd026385970643f9c5b8ecfeb"},
-]
-pluggy = [
-    {file = "pluggy-1.0.0-py2.py3-none-any.whl", hash = "sha256:74134bbf457f031a36d68416e1509f34bd5ccc019f0bcc952c7b909d06b37bd3"},
-    {file = "pluggy-1.0.0.tar.gz", hash = "sha256:4224373bacce55f955a878bf9cfa763c1e360858e330072059e10bad68531159"},
-]
-py = [
-    {file = "py-1.11.0-py2.py3-none-any.whl", hash = "sha256:607c53218732647dff4acdfcd50cb62615cedf612e72d1724fb1a0cc6405b378"},
-    {file = "py-1.11.0.tar.gz", hash = "sha256:51c75c4126074b472f746a24399ad32f6053d1b34b68d2fa41e558e6f4a98719"},
-]
-pycodestyle = [
-    {file = "pycodestyle-2.7.0-py2.py3-none-any.whl", hash = "sha256:514f76d918fcc0b55c6680472f0a37970994e07bbb80725808c17089be302068"},
-    {file = "pycodestyle-2.7.0.tar.gz", hash = "sha256:c389c1d06bf7904078ca03399a4816f974a1d590090fecea0c63ec26ebaf1cef"},
-]
-pyflakes = [
-    {file = "pyflakes-2.3.1-py2.py3-none-any.whl", hash = "sha256:7893783d01b8a89811dd72d7dfd4d84ff098e5eed95cfa8905b22bbffe52efc3"},
-    {file = "pyflakes-2.3.1.tar.gz", hash = "sha256:f5bc8ecabc05bb9d291eb5203d6810b49040f6ff446a756326104746cc00c1db"},
-]
-pyparsing = [
-    {file = "pyparsing-3.0.6-py3-none-any.whl", hash = "sha256:04ff808a5b90911829c55c4e26f75fa5ca8a2f5f36aa3a51f68e27033341d3e4"},
-    {file = "pyparsing-3.0.6.tar.gz", hash = "sha256:d9bdec0013ef1eb5a84ab39a3b3868911598afa494f5faa038647101504e2b81"},
-]
-pytest = [
-    {file = "pytest-6.2.5-py3-none-any.whl", hash = "sha256:7310f8d27bc79ced999e760ca304d69f6ba6c6649c0b60fb0e04a4a77cacc134"},
-    {file = "pytest-6.2.5.tar.gz", hash = "sha256:131b36680866a76e6781d13f101efb86cf674ebb9762eb70d3082b6f29889e89"},
-]
-pytest-cov = [
-    {file = "pytest-cov-2.12.1.tar.gz", hash = "sha256:261ceeb8c227b726249b376b8526b600f38667ee314f910353fa318caa01f4d7"},
-    {file = "pytest_cov-2.12.1-py2.py3-none-any.whl", hash = "sha256:261bb9e47e65bd099c89c3edf92972865210c36813f80ede5277dceb77a4a62a"},
-]
-toml = [
-    {file = "toml-0.10.2-py2.py3-none-any.whl", hash = "sha256:806143ae5bfb6a3c6e736a764057db0e6a0e05e338b5630894a5f779cabb4f9b"},
-    {file = "toml-0.10.2.tar.gz", hash = "sha256:b3bda1d108d5dd99f4a20d24d9c348e91c4db7ab1b749200bded2f839ccbe68f"},
-]
-typing-extensions = [
-    {file = "typing_extensions-4.0.1-py3-none-any.whl", hash = "sha256:7f001e5ac290a0c0401508864c7ec868be4e701886d5b573a9528ed3973d9d3b"},
-    {file = "typing_extensions-4.0.1.tar.gz", hash = "sha256:4ca091dea149f945ec56afb48dae714f21e8692ef22a395223bcd328961b6a0e"},
-]
+appnope = []
+asttokens = []
+backcall = []
+black = []
+click = []
+colorama = []
+colored = []
+coverage = []
+decorator = []
+exceptiongroup = []
+executing = []
+flake8 = []
+iniconfig = []
+ipython = []
+jedi = []
+matplotlib-inline = []
+mccabe = []
+mypy = []
+mypy-extensions = []
+numpy = []
+packaging = []
+parso = []
+pathspec = []
+pexpect = []
+pickleshare = []
+platformdirs = []
+pluggy = []
+prompt-toolkit = []
+ptyprocess = []
+pure-eval = []
+pycodestyle = []
+pyflakes = []
+pygments = []
+pyparsing = []
+pytest = []
+pytest-cov = []
+scipy = []
+six = []
+stack-data = []
+syrupy = []
+tomli = []
+traitlets = []
+typing-extensions = []
+wcwidth = []

pyproject.toml

@@ -1,29 +1,35 @@
 [tool.poetry]
 name = "pdme"
-version = "0.0.1"
+version = "1.5.0"
 description = "Python dipole model evaluator"
 authors = ["Deepak <dmallubhotla+github@gmail.com>"]
 license = "GPL-3.0-only"
 readme = "README.md"
 
 [tool.poetry.dependencies]
-python = "^3.8,<3.10"
-numpy = "^1.21.1"
+python = ">=3.8.1,<3.10"
+numpy = "^1.22.3"
+scipy = "~1.10"
 
 [tool.poetry.dev-dependencies]
 pytest = ">=6"
-flake8 = "^3.8.4"
-pytest-cov = "^2.10.1"
-mypy = "^0.910"
+flake8 = "^4.0.0"
+pytest-cov = "^4.1.0"
+mypy = "^0.961"
+ipython = "^8.2.0"
+black = "^22.3.0"
+syrupy = "^4.0.8"
 
 [build-system]
-requires = ["poetry>=0.12"]
-build-backend = "poetry.masonry.api"
+requires = ["poetry-core>=1.0.0"]
+build-backend = "poetry.core.masonry.api"
 
 [tool.pytest.ini_options]
 testpaths = ["tests"]
-addopts = "--junitxml pytest.xml --cov pdme --cov-report=xml:coverage.xml --cov-fail-under=90 --cov-report=html"
+addopts = "--junitxml pytest.xml --cov pdme --cov-report=xml:coverage.xml --cov-fail-under=50 --cov-report=html"
 junit_family = "xunit1"
+log_format = "%(asctime)s | %(levelname)s | %(pathname)s:%(lineno)d | %(message)s"
+log_level = "WARNING"
 
 [tool.mypy]
 plugins = "numpy.typing.mypy_plugin"

renovate.json

@@ -0,0 +1,3 @@
+{
+	"$schema": "https://docs.renovatebot.com/renovate-schema.json"
+}

scripts/patch.sh

@@ -1,29 +0,0 @@
-#!/usr/bin/env bash
-set -Eeuo pipefail
-
-if [ -z "$(git status --porcelain)" ]; then 
-	# Working directory clean
-	branch_name=$(git symbolic-ref -q HEAD)
-	branch_name=${branch_name##refs/heads/}
-	branch_name=${branch_name:-HEAD}
-	
-	poetry version patch
-	version=`sed 's/version = "\([0-9]*.[0-9]*.[0-9]*\)"/\1/p' -n <pyproject.toml`
-	read -p "Create commit for version $version? " -n 1 -r
-	echo    # (optional) move to a new line
-	if [[ $REPLY =~ ^[Yy]$ ]]
-	then
-		# do dangerous stuff
-		echo "Creating a new patch"
-		git add pyproject.toml
-		git commit -m "Created version $version"
-		git tag -a "$version" -m "patch.sh created version $version"
-		git push --tags
-	else
-		echo "Surrendering, clean up by reverting pyproject.toml..."
-		exit 2
-	fi
-else 
-	echo "Can't create patch version, working tree unclean..."
-	exit 1
-fi

scripts/release.sh

@@ -0,0 +1,45 @@
+#!/usr/bin/env bash
+set -Eeuo pipefail
+
+if [ -z "$(git status --porcelain)" ]; then
+	branch_name=$(git symbolic-ref -q HEAD)
+	branch_name=${branch_name##refs/heads/}
+	branch_name=${branch_name:-HEAD}
+	if [ $branch_name != "master" ]; then
+		echo "The current branch is not master!"
+		echo "I'd feel uncomfortable releasing from here..."
+		exit 3
+	fi
+
+	release_needed=false
+	if \
+		{ git log "$( git describe --tags --abbrev=0 )..HEAD" --format='%s' | cut -d: -f1 | sort -u | sed -e 's/([^)]*)//' | grep -q -i -E '^feat|fix|perf|refactor|revert$' ; } || \
+		{ git log "$( git describe --tags --abbrev=0 )..HEAD" --format='%s' | cut -d: -f1 | sort -u | sed -e 's/([^)]*)//' | grep -q -E '\!$' ; } || \
+		{ git log "$( git describe --tags --abbrev=0 )..HEAD" --format='%b' | grep -q -E '^BREAKING CHANGE:' ; }
+	then
+		release_needed=true
+	fi
+	
+	if ! [ "$release_needed" = true ]; then
+		echo "No release needed..."
+		exit 0
+	fi
+
+	# Working directory clean
+	echo "Doing a dry run..."
+	npx standard-version --dry-run
+	read -p "Does that look good? [y/N] " -n 1 -r
+	echo    # (optional) move to a new line
+	if [[ $REPLY =~ ^[Yy]$ ]]
+	then
+		# do dangerous stuff
+		npx standard-version
+		git push --follow-tags origin master
+	else
+		echo "okay, never mind then..."
+		exit 2
+	fi
+else 
+	echo "Can't create release, working tree unclean..."
+	exit 1
+fi

scripts/standard-version/pyproject-updater.js

@@ -0,0 +1,11 @@
+const pattern = /(\[tool\.poetry\]\nname = "pdme"\nversion = ")(?<vers>\d+\.\d+\.\d+)(")/mg;
+
+module.exports.readVersion = function (contents) {
+	const result = pattern.exec(contents);
+	return result.groups.vers;
+}
+
+module.exports.writeVersion = function (contents, version) {
+	const newContents = contents.replace(pattern, `$1${version}$3`);
+	return newContents;
+}

tests/calculations/__snapshots__/test_calculations.ambr

@@ -0,0 +1,81 @@
+# serializer version: 1
+# name: test_multiple_electric_field_alphas
+  list([
+    list([
+      0.009579434215333742,
+      0.007714411624737791,
+      0.0035604976729559034,
+    ]),
+    list([
+      0.005036717012002481,
+      0.006259221141129298,
+      -0.009682232702684382,
+    ]),
+    list([
+      0.01503516326014651,
+      0.012028130608117207,
+      0.009021097956087907,
+    ]),
+    list([
+      0.0033871215792458027,
+      0.0038182097802407235,
+      -0.005604146612933966,
+    ]),
+    list([
+      0.007409666906433089,
+      0.006778734778841714,
+      0.003467602973242188,
+    ]),
+    list([
+      0.004730939083250055,
+      0.0046546336141653774,
+      -0.011766303332857395,
+    ]),
+    list([
+      -0.010766266772985707,
+      -0.012631289363581657,
+      -0.010851040527103707,
+    ]),
+    list([
+      -0.008410828408392494,
+      -0.020391368873835285,
+      0.009877833363344673,
+    ]),
+    list([
+      -0.015035163260146511,
+      -0.018042195912175815,
+      -0.021049228564205113,
+    ]),
+    list([
+      -0.005850482727788205,
+      -0.012193637685284892,
+      0.0054809785555068454,
+    ]),
+    list([
+      -0.007682229585552562,
+      -0.010584517372472879,
+      -0.009352937859414517,
+    ]),
+    list([
+      -0.009767100042838822,
+      -0.020831393060117175,
+      0.014024945217763873,
+    ]),
+  ])
+# ---
+# name: test_multiple_electric_potential_alphas
+  list([
+    0.05035560994609065,
+    -0.03369221379873504,
+    0.07216878364870323,
+    -0.024633611485424027,
+    0.04496779459769236,
+    -0.03108684810509794,
+    -0.08019597139562586,
+    0.08293468011996319,
+    -0.10825317547305485,
+    0.060044427995721066,
+    -0.06935710692186449,
+    0.08733923991432278,
+  ])
+# ---

tests/calculations/test_calculations.py

@@ -0,0 +1,102 @@
+import pdme.calculations
+import pytest
+import numpy
+import numpy.testing
+
+
+# generated in mathematica to compare here
+beta_test_data = [
+	[-2, -2, 0.8072976151],
+	[-1, -2, 0.008105366193],
+	[0, -2, 0.00008105691406],
+	[1, -2, 8.105694659e-7],
+	[2, -2, 8.105694691e-9],
+	[-2, -1, 5.768008783],
+	[-1, -1, 0.08072976151],
+	[0, -1, 0.0008105366193],
+	[1, -1, 8.105691406e-6],
+	[2, -1, 8.105694659e-8],
+	[-2, 0, 1.951840272],
+	[-1, 0, 0.5768008783],
+	[0, 0, 0.008072976151],
+	[1, 0, 0.00008105366193],
+	[2, 0, 8.105691406e-7],
+	[-2, 1, 0.1999506642],
+	[-1, 1, 0.1951840272],
+	[0, 1, 0.05768008783],
+	[1, 1, 0.0008072976151],
+	[2, 1, 8.105366193e-6],
+	[-2, 2, 0.01999995065],
+	[-1, 2, 0.01999506642],
+	[0, 2, 0.01951840272],
+	[1, 2, 0.005768008783],
+	[2, 2, 0.00008072976151],
+]
+
+
+@pytest.mark.parametrize("f, w, expected", beta_test_data)
+def test_calculations_beta_func(f, w, expected):
+	# there's nothing special about the 5 * 10^ f and 10^w passing in logs
+	# this was just to get a variety of orders of magnitude in results.
+	actual = pdme.calculations.telegraph_beta(5 * 10**f, 10**w)
+	numpy.testing.assert_allclose(actual, expected, atol=0, rtol=1e-8)
+
+
+def test_multiple_electric_potential_alphas(snapshot):
+	"""
+	Manually compare these with mathematica stuff because manually including a list is a bit annoying.
+	Basically just visually compare the snapshot values to make sure they're actually correct.
+	"""
+	dipole_ps = [
+		numpy.array([1, 2, 3]),
+		numpy.array([-4, -3, -2]),
+	]
+	dipole_ss = [
+		numpy.array([0, 0, 1]),
+		numpy.array([1, 1, 1]),
+		numpy.array([0, -0.5, 0.5]),
+	]
+
+	test_rs = [
+		numpy.array([5, 5, 5]),
+		numpy.array([-4, -6, 2]),
+	]
+
+	actuals = [
+		pdme.calculations.electric_potential(p, s, r)
+		for p in dipole_ps
+		for s in dipole_ss
+		for r in test_rs
+	]
+
+	assert actuals == snapshot
+
+
+def test_multiple_electric_field_alphas(snapshot):
+	"""
+	Manually compare these with mathematica stuff because manually including a list is a bit annoying.
+	Basically just visually compare the snapshot values to make sure they're actually correct.
+	"""
+	dipole_ps = [
+		numpy.array([1, 2, 3]),
+		numpy.array([-4, -3, -2]),
+	]
+	dipole_ss = [
+		numpy.array([0, 0, 1]),
+		numpy.array([1, 1, 1]),
+		numpy.array([0, -0.5, 0.5]),
+	]
+
+	test_rs = [
+		numpy.array([5, 5, 5]),
+		numpy.array([-4, -6, 2]),
+	]
+
+	actuals = [
+		pdme.calculations.electric_field(p, s, r).tolist()
+		for p in dipole_ps
+		for s in dipole_ss
+		for r in test_rs
+	]
+
+	assert actuals == snapshot

tests/inputs/test_dot_inputs.py

@@ -0,0 +1,38 @@
+import pdme.inputs
+
+
+def test_inputs_with_frequency_range():
+	dot1 = [1, 2, 3]
+	dot2 = [2, 4, 6]
+
+	frequencies = [5, 7, 9]
+
+	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),
+	]
+
+	actual = pdme.inputs.inputs_with_frequency_range([dot1, dot2], frequencies)
+
+	assert sorted(actual) == sorted(expected), "Did not actually match dot inputs!"
+
+
+def test_input_pairs_with_frequency_range():
+	dot1 = [1, 2, 3]
+	dot2 = [2, 4, 6]
+
+	frequencies = [5, 7, 9]
+
+	expected = [
+		([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)
+
+	assert sorted(actual) == sorted(expected), "Did not actually match dot inputs!"

tests/measurement/test_dipole_noise_measurement.py

@@ -0,0 +1,174 @@
+import numpy
+import pdme.measurement
+
+
+def test_static_dipole_electric_potential():
+	d1 = pdme.measurement.OscillatingDipole(
+		numpy.array((1, 2, 3)), numpy.array([4, 5, 6]), 7
+	)
+	d2 = pdme.measurement.OscillatingDipole(
+		numpy.array([2, 5, 3]), numpy.array([4, -5, -6]), 2
+	)
+	dipoles = pdme.measurement.OscillatingDipoleArrangement([d1, d2])
+
+	dot_position1 = numpy.array([-1, -1, -1])
+	dot_frequency1 = 11
+
+	expected_v1 = 0.00001221710876727626
+	expected_v2 = 7.257229625870065e-6
+
+	numpy.testing.assert_allclose(
+		d1.s_electric_potential_at_position(dot_position1, dot_frequency1),
+		expected_v1,
+		err_msg="Voltage at dot isn't as expected.",
+	)
+
+	dot_measurements = dipoles.get_potential_dot_measurements(
+		[(dot_position1, dot_frequency1)]
+	)
+	assert len(dot_measurements) == 1, "Should have only had one dot measurement."
+	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(
+		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_static_dipole_electric_field_x():
+	d1 = pdme.measurement.OscillatingDipole(
+		numpy.array((1, 2, 3)), numpy.array([4, 5, 6]), 7
+	)
+	d2 = pdme.measurement.OscillatingDipole(
+		numpy.array([2, 5, 3]), numpy.array([4, -5, -6]), 2
+	)
+	# dipoles = pdme.measurement.OscillatingDipoleArrangement([d1, d2])
+
+	dot_position1 = numpy.array([-1, -1, -1])
+	dot_frequency1 = 11
+
+	# these should be true for electric field
+	expected_v1 = 1.479556451978925e-7
+	expected_v2 = 6.852024308908262e-7
+
+	numpy.testing.assert_allclose(
+		d1.s_electric_fieldx_at_position(dot_position1, dot_frequency1),
+		expected_v1,
+		err_msg="Fieldx noise at dot isn't as expected.",
+	)
+
+	numpy.testing.assert_allclose(
+		d2.s_electric_fieldx_at_position(dot_position1, dot_frequency1),
+		expected_v2,
+		err_msg="Fieldx at dot isn't as expected.",
+	)
+
+
+def test_dipole_dot_pair():
+	d1 = pdme.measurement.OscillatingDipole(
+		numpy.array([1, 2, 3]), numpy.array([4, 5, 6]), 7
+	)
+	dipoles = pdme.measurement.OscillatingDipoleArrangement([d1])
+
+	dot_position1 = numpy.array([-1, -1, -1])
+	dot_position2 = numpy.array([1, 2, 3])
+	dot_frequency = 8
+	expected_sij = 0.00008692058236162933
+
+	numpy.testing.assert_allclose(
+		d1.s_electric_potential_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_potential_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():
+	d1 = pdme.measurement.OscillatingDipole(
+		numpy.array([1, 2, 3]), numpy.array([4, 5, 6]), 7
+	)
+	dipoles = pdme.measurement.OscillatingDipoleArrangement([d1])
+
+	dot_position1 = numpy.array([-1, -1, -1])
+	dot_position2 = numpy.array([1, 2, 3])
+	dot_frequency = 8
+	expected_sij = 0.00008692058236162933
+
+	actuals = dipoles.get_percent_range_potential_dot_pair_measurements(
+		[(dot_position1, dot_position2, dot_frequency)], 0.5, 1.5
+	)
+	assert len(actuals) == 1, "should have only been one pair"
+	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.",
+	)
+
+
+def test_range_dipole_measurements():
+	d1 = pdme.measurement.OscillatingDipole(
+		numpy.array([1, 2, 3]), numpy.array([4, 5, 6]), 7
+	)
+	d2 = pdme.measurement.OscillatingDipole(
+		numpy.array([2, 5, 3]), numpy.array([4, -5, -6]), 2
+	)
+	dipoles = pdme.measurement.OscillatingDipoleArrangement([d1, d2])
+
+	dot_position1 = numpy.array([-1, -1, -1])
+	dot_frequency1 = 11
+	expected_v1 = 0.00001221710876727626
+	expected_v2 = 7.257229625870065e-6
+
+	numpy.testing.assert_allclose(
+		d1.s_electric_potential_at_position(dot_position1, dot_frequency1),
+		expected_v1,
+		err_msg="Voltage at dot isn't as expected.",
+	)
+
+	range_dot_measurements = dipoles.get_percent_range_potential_dot_measurements(
+		[(dot_position1, dot_frequency1)], 0.5, 1.5
+	)
+	assert len(range_dot_measurements) == 1, "Should have only had one dot measurement."
+	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(
+		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.",
+	)

tests/measurement/test_dipole_to_array.py

@@ -0,0 +1,24 @@
+import numpy
+import pdme.measurement
+
+
+def test_dipole_to_array():
+	d1 = pdme.measurement.OscillatingDipole((1, 2, 3), (4, 5, 6), 7)
+	d2 = pdme.measurement.OscillatingDipole((1.5, 2.5, 3.5), (4.5, 5.5, 6.5), 7.5)
+
+	expected1 = numpy.array([1, 2, 3, 4, 5, 6, 7])
+	numpy.testing.assert_array_equal(
+		expected1, d1.to_flat_array(), err_msg="Didn't convert 1 correctly to array"
+	)
+
+	expected2 = numpy.array([1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5])
+	numpy.testing.assert_array_equal(
+		expected2, d2.to_flat_array(), err_msg="Didn't convert 2 correctly to array"
+	)
+
+	arrangement = pdme.measurement.OscillatingDipoleArrangement([d1, d2])
+	numpy.testing.assert_array_equal(
+		[expected1, expected2],
+		arrangement.to_numpy_array(),
+		err_msg="Didn't convert multiple dipoles right",
+	)

tests/measurement/test_measurement_range_swap.py

@@ -0,0 +1,75 @@
+from pdme.measurement import DotRangeMeasurement
+from pdme.measurement import DotPairRangeMeasurement
+import numpy
+
+
+def test_swap_high_low():
+	actual_high = 2
+	actual_low = 1
+	m1 = DotRangeMeasurement(actual_high, actual_low, 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(
+		[m2.v_low, m2.v_high],
+		[actual_low, actual_high],
+		err_msg="Highs were wrong without swap",
+	)
+
+
+def test_swap_high_low_negative():
+	actual_high = -1
+	actual_low = -2
+	m1 = DotRangeMeasurement(actual_high, actual_low, 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(
+		[m2.v_low, m2.v_high],
+		[actual_low, actual_high],
+		err_msg="Highs were wrong without swap, negative",
+	)
+
+
+def test_swap_high_low_pair():
+	actual_high = 2
+	actual_low = 1
+	m1 = DotPairRangeMeasurement(actual_high, actual_low, 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(
+		[m2.v_low, m2.v_high],
+		[actual_low, actual_high],
+		err_msg="Highs were wrong without swap",
+	)
+
+
+def test_swap_high_low_pair_negative():
+	actual_high = -1
+	actual_low = -2
+	m1 = DotPairRangeMeasurement(actual_high, actual_low, 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(
+		[m2.v_low, m2.v_high],
+		[actual_low, actual_high],
+		err_msg="Highs were wrong without swap, negative",
+	)

tests/measurement/test_measurement_to_arrays.py

@@ -0,0 +1,57 @@
+import numpy
+import pdme.measurement.input_types
+from pdme.measurement.dot_measure import DotRangeMeasurement
+from pdme.measurement.dot_pair_measure import DotPairRangeMeasurement
+
+
+def test_inputs_to_array():
+	i1 = ([1, 2, 3], 5)
+	i2 = ([-1, 4, -2], 9)
+
+	actual = pdme.measurement.input_types.dot_inputs_to_array([i1, i2])
+	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"
+	)
+
+
+def test_pair_inputs_to_array():
+	i1 = ([1, 2, 3], [-1, 4, -2], 5)
+	i2 = ([-1, 4, -2], [6, 7, 8], 9)
+
+	actual = pdme.measurement.input_types.dot_pair_inputs_to_array([i1, i2])
+	expected = numpy.array(
+		[
+			[[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"
+	)
+
+
+def test_ranges_to_array():
+	m1 = DotRangeMeasurement(1, 2, 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])
+	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")
+
+
+def test_pair_ranges_to_array():
+	m1 = DotPairRangeMeasurement(1, 2, 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])
+	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")

tests/model/__snapshots__/test_log_spaced_fixed_orientation_fixedxy_orientation_mcmc.ambr

@@ -0,0 +1,55 @@
+# serializer version: 1
+# name: test_log_spaced_fixedxy_orientation_mcmc_basic
+  list([
+    tuple(
+      3984.4199552664,
+      array([[ 9.55610128,  2.94634152,  0.        ,  9.21529051, -2.46576127,
+           2.42481096,  9.19034554]]),
+    ),
+    tuple(
+      8583.8032728084,
+      array([[ 9.99991539,  0.04113671,  0.        ,  8.71258954, -2.26599865,
+           2.60452102,  6.37042214]]),
+    ),
+    tuple(
+      6215.5671525452,
+      array([[ 9.81950685, -1.89137124,  0.        ,  8.90637055, -2.48043039,
+           2.28444435,  8.84239221]]),
+    ),
+    tuple(
+      424.7282747232,
+      array([[ 1.00028483,  9.94984574,  0.        ,  8.53064898, -2.59230757,
+           2.33774773,  8.6714416 ]]),
+    ),
+    tuple(
+      300.920463669,
+      array([[ 1.4003442 ,  9.90146636,  0.        ,  8.05557992, -2.6753126 ,
+           2.65915755, 13.02021385]]),
+    ),
+    tuple(
+      2400.0053792245,
+      array([[ 9.97761813,  0.66868263,  0.        ,  8.69171028, -2.73145011,
+           2.90140456, 19.94999593]]),
+    ),
+    tuple(
+      5001.4154377966,
+      array([[ 9.93976109, -1.09596962,  0.        ,  8.95245025, -2.59409162,
+           2.90140456,  9.75535945]]),
+    ),
+    tuple(
+      195.2191433934,
+      array([[ 0.20690762,  9.99785923,  0.        ,  9.59636585, -2.83240984,
+           2.90140456, 16.14771567]]),
+    ),
+    tuple(
+      2698.2459654601,
+      array([[-9.68130127, -2.50447712,  0.        ,  8.94823619, -2.92889659,
+           2.77065328, 13.63173263]]),
+    ),
+    tuple(
+      1193.6653292625,
+      array([[-6.16597091, -7.87278875,  0.        ,  9.62210721, -2.75993924,
+           2.77065328,  5.64553534]]),
+    ),
+  ])
+# ---

tests/model/__snapshots__/test_log_spaced_fixed_orientation_free_orientation_mcmc.ambr

@@ -0,0 +1,55 @@
+# serializer version: 1
+# name: test_log_spaced_free_orientation_mcmc_basic
+  list([
+    tuple(
+      3167.6603875494,
+      array([[ 9.60483896, -1.41627817, -2.3960853 , -4.76615152, -1.80902942,
+           2.11809123, 16.17452242]]),
+    ),
+    tuple(
+      3167.6603875494,
+      array([[ 9.60483896, -1.41627817, -2.3960853 , -4.76615152, -1.80902942,
+           2.11809123, 16.17452242]]),
+    ),
+    tuple(
+      3167.6603875494,
+      array([[ 9.60483896, -1.41627817, -2.3960853 , -4.76615152, -1.80902942,
+           2.11809123, 16.17452242]]),
+    ),
+    tuple(
+      2320.6090682509,
+      array([[ 4.1660069 , -8.11557337,  4.0965663 , -4.35968351, -1.97945216,
+           2.43615641, 12.92143144]]),
+    ),
+    tuple(
+      2320.6090682509,
+      array([[ 4.1660069 , -8.11557337,  4.0965663 , -4.35968351, -1.97945216,
+           2.43615641, 12.92143144]]),
+    ),
+    tuple(
+      2320.6090682509,
+      array([[ 4.1660069 , -8.11557337,  4.0965663 , -4.35968351, -1.97945216,
+           2.43615641, 12.92143144]]),
+    ),
+    tuple(
+      2248.0760851141,
+      array([[-1.71755535, -5.59925137,  8.10545419, -4.03306318, -1.81098441,
+           2.77407111, 32.28020575]]),
+    ),
+    tuple(
+      1663.3101472167,
+      array([[-5.16785855,  2.7558756 ,  8.10545419, -3.34620897, -1.74763642,
+           2.42770463, 52.98214008]]),
+    ),
+    tuple(
+      1329.2703365134,
+      array([[ -1.39600464,   9.69718343,  -2.00394725,  -2.59147366,
+           -1.91246681,   2.07361175, 123.01833742]]),
+    ),
+    tuple(
+      355.7695549121,
+      array([[  9.76047401,   0.84696075,  -2.00394725,  -3.04310053,
+           -1.99338573,   2.1185589 , 271.35743739]]),
+    ),
+  ])
+# ---

tests/model/__snapshots__/test_log_spaced_fixed_orientation_mcmc.ambr

@@ -0,0 +1,55 @@
+# serializer version: 1
+# name: test_log_spaced_fixed_orientation_mcmc_basic
+  list([
+    tuple(
+      50.5632126772,
+      array([[ 0.        ,  0.        , 10.        , -2.3960853 ,  4.23246234,
+           2.26169242, 39.39900844]]),
+    ),
+    tuple(
+      50.5632126772,
+      array([[ 0.        ,  0.        , 10.        , -2.3960853 ,  4.23246234,
+           2.26169242, 39.39900844]]),
+    ),
+    tuple(
+      47.4038652151,
+      array([[ 0.        ,  0.        , 10.        , -2.03666518,  4.14084039,
+           2.21309317, 47.82371559]]),
+    ),
+    tuple(
+      47.4038652151,
+      array([[ 0.        ,  0.        , 10.        , -2.03666518,  4.14084039,
+           2.21309317, 47.82371559]]),
+    ),
+    tuple(
+      47.4038652151,
+      array([[ 0.        ,  0.        , 10.        , -2.03666518,  4.14084039,
+           2.21309317, 47.82371559]]),
+    ),
+    tuple(
+      47.4038652151,
+      array([[ 0.        ,  0.        , 10.        , -2.03666518,  4.14084039,
+           2.21309317, 47.82371559]]),
+    ),
+    tuple(
+      22.9304785991,
+      array([[ 0.        ,  0.        , 10.        , -1.63019717,  3.97041764,
+           2.53115835, 38.2051999 ]]),
+    ),
+    tuple(
+      28.8090658953,
+      array([[ 0.        ,  0.        , 10.        , -1.14570315,  4.07709911,
+           2.48697441, 49.58615195]]),
+    ),
+    tuple(
+      28.8090658953,
+      array([[ 0.        ,  0.        , 10.        , -1.14570315,  4.07709911,
+           2.48697441, 49.58615195]]),
+    ),
+    tuple(
+      40.9699102596,
+      array([[ 0.        ,  0.        , 10.        , -0.50178755,  3.83878089,
+           2.93560796, 82.07827571]]),
+    ),
+  ])
+# ---

tests/model/__snapshots__/test_log_spaced_random_count_multiple_dipole_fixed_orientation_fixed_magnitude_model.ambr

@@ -0,0 +1,31 @@
+# serializer version: 1
+# name: test_random_count_multiple_dipole_fixed_or_fixed_mag_model_get_n_dipoles
+  list([
+    list([
+      list([
+        10.0,
+        0.0,
+        6.123233995736766e-16,
+        -2.3960852996076447,
+        4.232462337639554,
+        2.2616924238635443,
+        39.399008444891905,
+      ]),
+    ]),
+  ])
+# ---
+# name: test_random_count_multiple_dipole_fixed_or_fixed_mag_model_get_n_dipoles.1
+  list([
+    list([
+      list([
+        10.0,
+        0.0,
+        6.123233995736766e-16,
+        -2.3960852996076447,
+        4.232462337639554,
+        2.2616924238635443,
+        39.399008444891905,
+      ]),
+    ]),
+  ])
+# ---

tests/model/test_log_spaced_fixed_orientation_fixedxy_orientation_mcmc.py

@@ -0,0 +1,100 @@
+from pdme.model import (
+	LogSpacedRandomCountMultipleDipoleFixedMagnitudeXYModel,
+)
+import pdme.inputs
+import pdme.measurement.input_types
+import pdme.subspace_simulation
+import numpy
+import logging
+
+_logger = logging.getLogger(__name__)
+
+SEED_TO_USE = 42
+
+
+def get_cost_function():
+	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 = LogSpacedRandomCountMultipleDipoleFixedMagnitudeXYModel(
+		x_min,
+		x_max,
+		y_min,
+		y_max,
+		z_min,
+		z_max,
+		0,
+		max_frequency,
+		p_fixed,
+		1,
+		0.5,
+	)
+	model.rng = numpy.random.default_rng(SEED_TO_USE)
+
+	freqs = [0.01, 0.1, 1, 10, 100]
+	dot_positions = [[-1.5, 0, 0], [-0.5, 0, 0], [0.5, 0, 0], [1.5, 0, 0]]
+	dot_inputs = pdme.inputs.inputs_with_frequency_range(dot_positions, freqs)
+	dot_input_array = pdme.measurement.input_types.dot_inputs_to_array(dot_inputs)
+
+	actual_dipoles = model.get_dipoles(0, numpy.random.default_rng(SEED_TO_USE))
+	actual_measurements = actual_dipoles.get_potential_dot_measurements(dot_inputs)
+	actual_measurements_array = numpy.array([m.v for m in actual_measurements])
+
+	def cost_to_use(sample_dipoleses: numpy.ndarray) -> numpy.ndarray:
+		return pdme.subspace_simulation.proportional_costs_vs_actual_measurement(
+			dot_input_array, actual_measurements_array, sample_dipoleses
+		)
+
+	return cost_to_use
+
+
+def test_log_spaced_fixedxy_orientation_mcmc_basic(snapshot):
+
+	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 = LogSpacedRandomCountMultipleDipoleFixedMagnitudeXYModel(
+		x_min,
+		x_max,
+		y_min,
+		y_max,
+		z_min,
+		z_max,
+		0,
+		max_frequency,
+		p_fixed,
+		1,
+		0.5,
+	)
+	model.rng = numpy.random.default_rng(1234)
+
+	seed = model.get_monte_carlo_dipole_inputs(1, -1)
+
+	cost_function = get_cost_function()
+	stdev = pdme.subspace_simulation.DipoleStandardDeviation(2, 2, 1, 0.25, 0.5, 1)
+	stdevs = pdme.subspace_simulation.MCMCStandardDeviation([stdev])
+
+	chain = model.get_mcmc_chain(
+		seed[0],
+		cost_function,
+		10,
+		cost_function(seed)[0],
+		stdevs,
+		rng_arg=numpy.random.default_rng(1515),
+	)
+
+	chain_rounded = [(round(cost, 10), dipoles) for (cost, dipoles) in chain]
+
+	assert chain_rounded == snapshot

tests/model/test_log_spaced_fixed_orientation_free_orientation_mcmc.py

@@ -0,0 +1,100 @@
+from pdme.model import (
+	LogSpacedRandomCountMultipleDipoleFixedMagnitudeModel,
+)
+import pdme.inputs
+import pdme.measurement.input_types
+import pdme.subspace_simulation
+import numpy
+import logging
+
+_logger = logging.getLogger(__name__)
+
+SEED_TO_USE = 42
+
+
+def get_cost_function():
+	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 = LogSpacedRandomCountMultipleDipoleFixedMagnitudeModel(
+		x_min,
+		x_max,
+		y_min,
+		y_max,
+		z_min,
+		z_max,
+		0,
+		max_frequency,
+		p_fixed,
+		1,
+		0.5,
+	)
+	model.rng = numpy.random.default_rng(SEED_TO_USE)
+
+	freqs = [0.01, 0.1, 1, 10, 100]
+	dot_positions = [[-1.5, 0, 0], [-0.5, 0, 0], [0.5, 0, 0], [1.5, 0, 0]]
+	dot_inputs = pdme.inputs.inputs_with_frequency_range(dot_positions, freqs)
+	dot_input_array = pdme.measurement.input_types.dot_inputs_to_array(dot_inputs)
+
+	actual_dipoles = model.get_dipoles(0, numpy.random.default_rng(SEED_TO_USE))
+	actual_measurements = actual_dipoles.get_potential_dot_measurements(dot_inputs)
+	actual_measurements_array = numpy.array([m.v for m in actual_measurements])
+
+	def cost_to_use(sample_dipoleses: numpy.ndarray) -> numpy.ndarray:
+		return pdme.subspace_simulation.proportional_costs_vs_actual_measurement(
+			dot_input_array, actual_measurements_array, sample_dipoleses
+		)
+
+	return cost_to_use
+
+
+def test_log_spaced_free_orientation_mcmc_basic(snapshot):
+
+	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 = LogSpacedRandomCountMultipleDipoleFixedMagnitudeModel(
+		x_min,
+		x_max,
+		y_min,
+		y_max,
+		z_min,
+		z_max,
+		0,
+		max_frequency,
+		p_fixed,
+		1,
+		0.5,
+	)
+	model.rng = numpy.random.default_rng(1234)
+
+	seed = model.get_monte_carlo_dipole_inputs(1, -1)
+
+	cost_function = get_cost_function()
+	stdev = pdme.subspace_simulation.DipoleStandardDeviation(2, 2, 1, 0.25, 0.5, 1)
+	stdevs = pdme.subspace_simulation.MCMCStandardDeviation([stdev])
+
+	chain = model.get_mcmc_chain(
+		seed[0],
+		cost_function,
+		10,
+		cost_function(seed)[0],
+		stdevs,
+		rng_arg=numpy.random.default_rng(1515),
+	)
+
+	chain_rounded = [(round(cost, 10), dipoles) for (cost, dipoles) in chain]
+
+	assert chain_rounded == snapshot

tests/model/test_log_spaced_fixed_orientation_mcmc.py

@@ -0,0 +1,108 @@
+from pdme.model import (
+	LogSpacedRandomCountMultipleDipoleFixedMagnitudeFixedOrientationModel,
+)
+import pdme.inputs
+import pdme.measurement.input_types
+import pdme.subspace_simulation
+import numpy
+import logging
+
+_logger = logging.getLogger(__name__)
+
+SEED_TO_USE = 42
+
+
+def get_cost_function():
+	x_min = -10
+	x_max = 10
+	y_min = -5
+	y_max = 5
+	z_min = 2
+	z_max = 3
+	p_fixed = 10
+	theta = 0
+	phi = 0
+	max_frequency = 5
+
+	model = LogSpacedRandomCountMultipleDipoleFixedMagnitudeFixedOrientationModel(
+		x_min,
+		x_max,
+		y_min,
+		y_max,
+		z_min,
+		z_max,
+		0,
+		max_frequency,
+		p_fixed,
+		theta,
+		phi,
+		1,
+		0.5,
+	)
+	model.rng = numpy.random.default_rng(SEED_TO_USE)
+
+	freqs = [0.01, 0.1, 1, 10, 100]
+	dot_positions = [[-1.5, 0, 0], [-0.5, 0, 0], [0.5, 0, 0], [1.5, 0, 0]]
+	dot_inputs = pdme.inputs.inputs_with_frequency_range(dot_positions, freqs)
+	dot_input_array = pdme.measurement.input_types.dot_inputs_to_array(dot_inputs)
+
+	actual_dipoles = model.get_dipoles(0, numpy.random.default_rng(SEED_TO_USE))
+	actual_measurements = actual_dipoles.get_potential_dot_measurements(dot_inputs)
+	actual_measurements_array = numpy.array([m.v for m in actual_measurements])
+
+	def cost_to_use(sample_dipoleses: numpy.ndarray) -> numpy.ndarray:
+		return pdme.subspace_simulation.proportional_costs_vs_actual_measurement(
+			dot_input_array, actual_measurements_array, sample_dipoleses
+		)
+
+	return cost_to_use
+
+
+def test_log_spaced_fixed_orientation_mcmc_basic(snapshot):
+
+	x_min = -10
+	x_max = 10
+	y_min = -5
+	y_max = 5
+	z_min = 2
+	z_max = 3
+	p_fixed = 10
+	theta = 0
+	phi = 0
+	max_frequency = 5
+
+	model = LogSpacedRandomCountMultipleDipoleFixedMagnitudeFixedOrientationModel(
+		x_min,
+		x_max,
+		y_min,
+		y_max,
+		z_min,
+		z_max,
+		0,
+		max_frequency,
+		p_fixed,
+		theta,
+		phi,
+		1,
+		0.5,
+	)
+	model.rng = numpy.random.default_rng(1234)
+
+	seed = model.get_monte_carlo_dipole_inputs(1, -1)
+
+	cost_function = get_cost_function()
+	stdev = pdme.subspace_simulation.DipoleStandardDeviation(2, 2, 1, 0.25, 0.5, 1)
+	stdevs = pdme.subspace_simulation.MCMCStandardDeviation([stdev])
+
+	chain = model.get_mcmc_chain(
+		seed[0],
+		cost_function,
+		10,
+		cost_function(seed)[0],
+		stdevs,
+		rng_arg=numpy.random.default_rng(1515),
+	)
+
+	chain_rounded = [(round(cost, 10), dipoles) for (cost, dipoles) in chain]
+
+	assert chain_rounded == snapshot

tests/model/test_log_spaced_random_count_multiple_dipole_fixed_magnitude_model.py

@@ -0,0 +1,229 @@
+from pdme.model import LogSpacedRandomCountMultipleDipoleFixedMagnitudeModel
+import numpy
+import logging
+import pytest
+
+
+_logger = logging.getLogger(__name__)
+
+
+def test_random_count_multiple_dipole_wrong_probability():
+	with pytest.raises(ValueError):
+		LogSpacedRandomCountMultipleDipoleFixedMagnitudeModel(
+			-10, 10, -5, 5, 2, 3, 1, 2, 10, 5, 2
+		)
+
+
+def test_repr_random_count_multiple_dipole_fixed_mag():
+	model = LogSpacedRandomCountMultipleDipoleFixedMagnitudeModel(
+		-10, 10, -5, 5, 2, 3, 1, 2, 10, 5, 0.5
+	)
+	assert (
+		repr(model)
+		== "LogSpacedRandomCountMultipleDipoleFixedMagnitudeModel(-10, 10, -5, 5, 2, 3, 1, 2, 10, 5, 0.5)"
+	), "Repr should be what I want."
+
+
+def test_random_count_multiple_dipole_fixed_mag_model_get_dipoles():
+
+	p_fixed = 10
+
+	model = LogSpacedRandomCountMultipleDipoleFixedMagnitudeModel(
+		-10, 10, -5, 5, 2, 3, 0, 5, p_fixed, 1, 0.5
+	)
+
+	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([8.60141814, -4.50270821, -2.3960853])
+	expected_s = numpy.array([-4.76615152, -1.80902942, 2.11809123])
+	expected_w = 10**1.2088314662639255
+
+	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_random_count_multiple_dipole_fixed_mag_model_get_dipoles_multiple():
+
+	p_fixed = 10
+	dipole_count = 5
+
+	model = LogSpacedRandomCountMultipleDipoleFixedMagnitudeModel(
+		-10, 10, -5, 5, 2, 3, 0, 5, p_fixed, dipole_count, 0.5
+	)
+
+	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_random_count_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 = LogSpacedRandomCountMultipleDipoleFixedMagnitudeModel(
+		x_min, x_max, y_min, y_max, z_min, z_max, 0, max_frequency, p_fixed, 1, 0.5
+	)
+	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([8.60141814, -4.50270821, -2.3960853])
+	expected_s = numpy.array([-4.76615152, -1.80902942, 2.11809123])
+	expected_w = 10**1.2088314662639255
+
+	numpy.testing.assert_allclose(
+		dipoles[0].p,
+		expected_p,
+		err_msg="Random multiple dipole p wasn't as expected",
+		rtol=1e-06,
+	)
+	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) in (
+			0,
+			1,
+		), "Should have either zero or one dipole generated."
+
+		if len(dipoles) > 0:
+			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 < 10 ** max_frequency and dipoles[0].w > 10**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_random_count_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 = LogSpacedRandomCountMultipleDipoleFixedMagnitudeModel(
+		x_min, x_max, y_min, y_max, z_min, z_max, 0, max_frequency, p_fixed, 1, 0.5
+	)
+	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,
+					-4.76615152,
+					-1.80902942,
+					2.11809123,
+					10**1.2088314662639255,
+				]
+			]
+		]
+	)
+
+	numpy.testing.assert_allclose(
+		dipole_array,
+		expected_dipole_array,
+		err_msg="Should have had the expected output dipole array.",
+		rtol=1e-6,
+	)
+	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_random_count_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 = LogSpacedRandomCountMultipleDipoleFixedMagnitudeModel(
+		x_min,
+		x_max,
+		y_min,
+		y_max,
+		z_min,
+		z_max,
+		0,
+		max_frequency,
+		p_fixed,
+		num_dipoles,
+		0.5,
+	)
+	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_log_spaced_random_count_multiple_dipole_fixed_magnitude_xy_model.py

@@ -0,0 +1,206 @@
+from pdme.model import LogSpacedRandomCountMultipleDipoleFixedMagnitudeXYModel
+import numpy
+import logging
+import pytest
+
+
+_logger = logging.getLogger(__name__)
+
+
+def test_random_count_multiple_dipole_xy_wrong_probability():
+	with pytest.raises(ValueError):
+		LogSpacedRandomCountMultipleDipoleFixedMagnitudeXYModel(
+			-10, 10, -5, 5, 2, 3, 1, 2, 10, 5, 2
+		)
+
+
+def test_repr_random_count_multiple_dipole_fixed_mag_xy():
+	model = LogSpacedRandomCountMultipleDipoleFixedMagnitudeXYModel(
+		-10, 10, -5, 5, 2, 3, 1, 2, 10, 5, 0.5
+	)
+	assert (
+		repr(model)
+		== "LogSpacedRandomCountMultipleDipoleFixedMagnitudeXYModel(-10, 10, -5, 5, 2, 3, 1, 2, 10, 5, 0.5)"
+	), "Repr should be what I want."
+
+
+def test_random_count_multiple_dipole_fixed_mag_model_get_dipoles_multiple_xy():
+
+	p_fixed = 10
+	dipole_count = 5
+
+	model = LogSpacedRandomCountMultipleDipoleFixedMagnitudeXYModel(
+		-10, 10, -5, 5, 2, 3, 0, 5, p_fixed, dipole_count, 0.5
+	)
+
+	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_random_count_multiple_dipole_fixed_mag_model_get_dipoles_invariant_xy():
+
+	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 = LogSpacedRandomCountMultipleDipoleFixedMagnitudeXYModel(
+		x_min, x_max, y_min, y_max, z_min, z_max, 0, max_frequency, p_fixed, 1, 0.5
+	)
+	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."
+	for i in range(10):
+		dipole_arrangement = model.get_dipoles(max_frequency)
+		dipoles = dipole_arrangement.dipoles
+
+		assert len(dipoles) in (
+			0,
+			1,
+		), "Should have either zero or one dipole generated."
+
+		if len(dipoles) > 0:
+			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 < 10 ** max_frequency and dipoles[0].w > 10**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.",
+			)
+
+			_logger.warning(dipoles[0].p)
+			numpy.testing.assert_allclose(
+				dipoles[0].p[2],
+				0,
+				err_msg="Should have had zero z magnitude.",
+			)
+
+
+def test_random_count_multiple_dipole_fixed_mag_model_get_dipoles_invariant_monte_carlo_xy():
+
+	x_min = -10
+	x_max = 10
+	y_min = -5
+	y_max = 5
+	z_min = 2
+	z_max = 3
+	p_fixed = 10
+	max_frequency = 5
+	monte_carlo_n = 20
+
+	model = LogSpacedRandomCountMultipleDipoleFixedMagnitudeXYModel(
+		x_min, x_max, y_min, y_max, z_min, z_max, 0, max_frequency, p_fixed, 1, 0.5
+	)
+	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."
+	for i in range(10):
+		dipoles = model.get_monte_carlo_dipole_inputs(monte_carlo_n, max_frequency)
+
+		min_s = numpy.array([x_min, y_min, z_min])
+		max_s = numpy.array([x_max, y_max, z_max])
+
+		_logger.warning(dipoles)
+		_logger.warning(dipoles[:, 0, 0:3])
+		_logger.warning(dipoles[:, 0, 3:6])
+		_logger.warning(dipoles[:, 0, 6])
+
+		ps = dipoles[:, 0, 0:3]
+		ss = dipoles[:, 0, 3:6]
+		ws = dipoles[:, 0, 6]
+
+		numpy.testing.assert_equal(
+			numpy.logical_and(min_s < ss, max_s > ss).all(),
+			True,
+			f"Dipole location [{ss}] should have been between min [{min_s}] and max [{max_s}] bounds.",
+		)
+		assert (ws < 10**max_frequency).all() and (
+			ws > 10**0
+		).all(), "Dipole frequency should have been between 0 and max."
+
+		norms = numpy.linalg.norm(ps, axis=1)
+		filtered_norms = norms[norms > 0]
+		numpy.testing.assert_allclose(
+			filtered_norms,
+			p_fixed,
+			err_msg="Should have had the expected dipole moment magnitude.",
+		)
+
+		numpy.testing.assert_allclose(
+			ps[:, 2],
+			0,
+			err_msg="Should have had zero z magnitude.",
+		)
+
+
+def test_random_count_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 = LogSpacedRandomCountMultipleDipoleFixedMagnitudeXYModel(
+		x_min,
+		x_max,
+		y_min,
+		y_max,
+		z_min,
+		z_max,
+		0,
+		max_frequency,
+		p_fixed,
+		num_dipoles,
+		0.5,
+	)
+	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",
+	)
+
+	actual_shape = model.get_monte_carlo_dipole_inputs(
+		monte_carlo_n, max_frequency, rng_to_use=numpy.random.default_rng(1515)
+	).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_log_spaced_random_count_multiple_dipole_fixed_orientation_fixed_magnitude_model.py

@@ -0,0 +1,152 @@
+from pdme.model import (
+	LogSpacedRandomCountMultipleDipoleFixedMagnitudeFixedOrientationModel,
+)
+
+import numpy
+import logging
+import pytest
+
+
+_logger = logging.getLogger(__name__)
+
+
+def test_random_count_fixedorientation_multiple_dipole_wrong_probability():
+	with pytest.raises(ValueError):
+		LogSpacedRandomCountMultipleDipoleFixedMagnitudeFixedOrientationModel(
+			-10, 10, -5, 5, 2, 3, 1, 2, 10, 0, 0, 5, 2
+		)
+
+
+def test_repr_random_count_multiple_dipole_fixed_orientation_mag():
+	model = LogSpacedRandomCountMultipleDipoleFixedMagnitudeFixedOrientationModel(
+		-10, 10, -5, 5, 2, 3, 1, 2, 10, 0, 1, 5, 0.5
+	)
+	assert (
+		repr(model)
+		== "LogSpacedRandomCountMultipleDipoleFixedMagnitudeFixedOrientationModel(-10, 10, -5, 5, 2, 3, 1, 2, 10, 0, 1, 5, 0.5)"
+	), "Repr should be same as instantiation."
+
+
+def test_random_count_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
+	theta = 0
+	phi = 0
+	max_frequency = 5
+
+	model = LogSpacedRandomCountMultipleDipoleFixedMagnitudeFixedOrientationModel(
+		x_min,
+		x_max,
+		y_min,
+		y_max,
+		z_min,
+		z_max,
+		0,
+		max_frequency,
+		p_fixed,
+		theta,
+		phi,
+		1,
+		0.5,
+	)
+	model.rng = numpy.random.default_rng(1234)
+
+	for i in range(10):
+		dipole_arrangement = model.get_dipoles(max_frequency)
+		dipoles = dipole_arrangement.dipoles
+
+		assert len(dipoles) in (
+			0,
+			1,
+		), "Should have either zero or one dipole generated."
+
+		if len(dipoles) > 0:
+			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 < 10 ** max_frequency and dipoles[0].w > 10**0
+			), "Dipole frequency should have been between 0 and max."
+
+			numpy.testing.assert_allclose(
+				dipoles[0].p,
+				numpy.array([0, 0, p_fixed]),
+				err_msg="Should have had the expected dipole moment.",
+			)
+
+	custom_rng = numpy.random.default_rng(1234)
+	for i in range(10):
+		dipole_arrangement = model.get_dipoles(max_frequency, custom_rng)
+		dipoles = dipole_arrangement.dipoles
+
+		assert len(dipoles) in (
+			0,
+			1,
+		), "Should have either zero or one dipole generated."
+
+		if len(dipoles) > 0:
+			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 < 10 ** max_frequency and dipoles[0].w > 10**0
+			), "Dipole frequency should have been between 0 and max."
+
+			numpy.testing.assert_allclose(
+				dipoles[0].p,
+				numpy.array([0, 0, p_fixed]),
+				err_msg="Should have had the expected dipole moment.",
+			)
+
+
+def test_random_count_multiple_dipole_fixed_or_fixed_mag_model_get_n_dipoles(snapshot):
+	x_min = -10
+	x_max = 10
+	y_min = -5
+	y_max = 5
+	z_min = 2
+	z_max = 3
+	p_fixed = 10
+	theta = numpy.pi / 2
+	phi = 0
+	max_frequency = 5
+
+	model = LogSpacedRandomCountMultipleDipoleFixedMagnitudeFixedOrientationModel(
+		x_min,
+		x_max,
+		y_min,
+		y_max,
+		z_min,
+		z_max,
+		0,
+		max_frequency,
+		p_fixed,
+		theta,
+		phi,
+		1,
+		0.5,
+	)
+	model.rng = numpy.random.default_rng(1234)
+
+	actual_inputs = model.get_monte_carlo_dipole_inputs(1, max_frequency)
+	assert actual_inputs.tolist() == snapshot
+	actual_monte_carlo_inputs = model.get_monte_carlo_dipole_inputs(
+		1, max_frequency, numpy.random.default_rng(1234)
+	)
+	assert actual_monte_carlo_inputs.tolist() == snapshot

tests/model/test_model_interface.py

@@ -0,0 +1,16 @@
+from pdme.model import DipoleModel
+import pytest
+
+
+def test_model_interface_not_implemented_one_dipoles():
+	model = DipoleModel()
+
+	with pytest.raises(NotImplementedError):
+		model.get_dipoles(5)
+
+
+def test_model_interface_not_implemented_n_dipoles():
+	model = DipoleModel()
+
+	with pytest.raises(NotImplementedError):
+		model.get_monte_carlo_dipole_inputs(5, 10)

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_random_count_multiple_dipole_fixed_magnitude_model.py

@@ -0,0 +1,217 @@
+from pdme.model import RandomCountMultipleDipoleFixedMagnitudeModel
+import numpy
+import logging
+import pytest
+
+
+_logger = logging.getLogger(__name__)
+
+
+def test_random_count_multiple_dipole_wrong_probability():
+	with pytest.raises(ValueError):
+		RandomCountMultipleDipoleFixedMagnitudeModel(-10, 10, -5, 5, 2, 3, 10, 5, 2)
+
+
+def test_repr_random_count_multiple_dipole_fixed_mag():
+	model = RandomCountMultipleDipoleFixedMagnitudeModel(
+		-10, 10, -5, 5, 2, 3, 10, 5, 0.5
+	)
+	assert (
+		repr(model)
+		== "RandomCountMultipleDipoleFixedMagnitudeModel(-10, 10, -5, 5, 2, 3, 10, 5, 0.5)"
+	), "Repr should be what I want."
+
+
+def test_random_count_multiple_dipole_fixed_mag_model_get_dipoles():
+
+	p_fixed = 10
+
+	model = RandomCountMultipleDipoleFixedMagnitudeModel(
+		-10, 10, -5, 5, 2, 3, p_fixed, 1, 0.5
+	)
+
+	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([8.60141814, -4.50270821, -2.3960853])
+	expected_s = numpy.array([-4.76615152, -1.80902942, 2.11809123])
+	expected_w = 1.2088314662639255
+
+	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_random_count_multiple_dipole_fixed_mag_model_get_dipoles_multiple():
+
+	p_fixed = 10
+	dipole_count = 5
+
+	model = RandomCountMultipleDipoleFixedMagnitudeModel(
+		-10, 10, -5, 5, 2, 3, p_fixed, dipole_count, 0.5
+	)
+
+	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_random_count_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 = RandomCountMultipleDipoleFixedMagnitudeModel(
+		x_min, x_max, y_min, y_max, z_min, z_max, p_fixed, 1, 0.5
+	)
+	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([8.60141814, -4.50270821, -2.3960853])
+	expected_s = numpy.array([-4.76615152, -1.80902942, 2.11809123])
+	expected_w = 1.2088314662639255
+
+	numpy.testing.assert_allclose(
+		dipoles[0].p,
+		expected_p,
+		err_msg="Random multiple dipole p wasn't as expected",
+		rtol=1e-06,
+	)
+	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) in (
+			0,
+			1,
+		), "Should have either zero or one dipole generated."
+
+		if len(dipoles) > 0:
+			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_random_count_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 = RandomCountMultipleDipoleFixedMagnitudeModel(
+		x_min, x_max, y_min, y_max, z_min, z_max, p_fixed, 1, 0.5
+	)
+	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,
+					-4.76615152,
+					-1.80902942,
+					2.11809123,
+					1.96706517,
+				]
+			]
+		]
+	)
+
+	numpy.testing.assert_allclose(
+		dipole_array,
+		expected_dipole_array,
+		err_msg="Should have had the expected output dipole array.",
+		rtol=1e-6,
+	)
+	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_random_count_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 = RandomCountMultipleDipoleFixedMagnitudeModel(
+		x_min, x_max, y_min, y_max, z_min, z_max, p_fixed, num_dipoles, 0.5
+	)
+	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/subspace_simulation/__snapshots__/test_costs.ambr

@@ -0,0 +1,4 @@
+# serializer version: 1
+# name: test_proportional_costs
+  7000.0
+# ---

tests/subspace_simulation/__snapshots__/test_sort_dipoles.ambr

@@ -0,0 +1,94 @@
+# serializer version: 1
+# name: test_sort_dipoles_by_freq
+  list([
+    list([
+      100.0,
+      200.0,
+      300.0,
+      400.0,
+      500.0,
+      600.0,
+      0.07,
+    ]),
+    list([
+      1.0,
+      2.0,
+      3.0,
+      4.0,
+      5.0,
+      6.0,
+      7.0,
+    ]),
+    list([
+      10.0,
+      200.0,
+      30.0,
+      41.0,
+      315.0,
+      0.31,
+      100.0,
+    ]),
+  ])
+# ---
+# name: test_sort_dipoleses_by_freq
+  list([
+    list([
+      list([
+        100.0,
+        200.0,
+        300.0,
+        400.0,
+        500.0,
+        600.0,
+        0.07,
+      ]),
+      list([
+        1.0,
+        2.0,
+        3.0,
+        4.0,
+        5.0,
+        6.0,
+        7.0,
+      ]),
+      list([
+        10.0,
+        200.0,
+        30.0,
+        41.0,
+        315.0,
+        0.31,
+        100.0,
+      ]),
+    ]),
+    list([
+      list([
+        1.0,
+        1.0,
+        1.0,
+        1.0,
+        1.0,
+        1.0,
+        100.0,
+      ]),
+      list([
+        22.0,
+        22.2,
+        2.2,
+        222.0,
+        22.0,
+        2.0,
+        200.0,
+      ]),
+      list([
+        33.0,
+        33.3,
+        33.0,
+        3.3,
+        0.33,
+        0.3,
+        300.0,
+      ]),
+    ]),
+  ])
+# ---

tests/subspace_simulation/__snapshots__/test_stdevs.ambr

@@ -0,0 +1,22 @@
+# serializer version: 1
+# name: test_return_four
+  DipoleStandardDeviation(p_phi_step=1, p_theta_step=2, rx_step=3, ry_step=4, rz_step=5, w_log_step=6)
+# ---
+# name: test_return_four.1
+  DipoleStandardDeviation(p_phi_step=10, p_theta_step=20, rx_step=30, ry_step=40, rz_step=50, w_log_step=60)
+# ---
+# name: test_return_four.2
+  DipoleStandardDeviation(p_phi_step=0.1, p_theta_step=0.2, rx_step=0.3, ry_step=0.4, rz_step=0.5, w_log_step=0.6)
+# ---
+# name: test_return_four.3
+  DipoleStandardDeviation(p_phi_step=1, p_theta_step=2, rx_step=3, ry_step=4, rz_step=5, w_log_step=6)
+# ---
+# name: test_return_four.4
+  DipoleStandardDeviation(p_phi_step=10, p_theta_step=20, rx_step=30, ry_step=40, rz_step=50, w_log_step=60)
+# ---
+# name: test_return_four.5
+  DipoleStandardDeviation(p_phi_step=0.1, p_theta_step=0.2, rx_step=0.3, ry_step=0.4, rz_step=0.5, w_log_step=0.6)
+# ---
+# name: test_return_one
+  DipoleStandardDeviation(p_phi_step=1, p_theta_step=2, rx_step=3, ry_step=4, rz_step=5, w_log_step=6)
+# ---

tests/subspace_simulation/test_costs.py

@@ -0,0 +1,31 @@
+import pdme.subspace_simulation
+import pdme.subspace_simulation.mcmc_costs
+import numpy
+
+
+def test_proportional_costs(snapshot):
+	a = numpy.array([2, 4, 5, 6, 7, 8, 10])
+	b = numpy.array([51, 13, 1, 31, 0.001, 3, 1])
+
+	actual_result = pdme.subspace_simulation.proportional_cost(a, b).tolist()
+	assert actual_result == snapshot
+
+
+def test_squared_costs_manual():
+	target = numpy.array([100, 400, 900])
+	approx1 = numpy.array([0, 400, 800])
+	approx2 = numpy.array([200, 400, 600])
+
+	expected1 = 1.0123456790123457
+	expected2 = 1.1111111111111111
+
+	actual1 = pdme.subspace_simulation.mcmc_costs.relative_square_diffs(approx1, target)
+	assert actual1 == expected1
+
+	actual2 = pdme.subspace_simulation.mcmc_costs.relative_square_diffs(approx2, target)
+	assert actual2 == expected2
+
+	combined_actual = pdme.subspace_simulation.mcmc_costs.relative_square_diffs(
+		numpy.array([approx1, approx2]), target
+	)
+	numpy.testing.assert_allclose(combined_actual, [expected1, expected2], rtol=1e-14)

tests/subspace_simulation/test_sort_dipoles.py

@@ -0,0 +1,40 @@
+import numpy
+import pdme.subspace_simulation
+
+
+def test_sort_dipoles_by_freq(snapshot):
+	orig = numpy.array(
+		[
+			[1, 2, 3, 4, 5, 6, 7],
+			[100, 200, 300, 400, 500, 600, 0.07],
+			[10, 200, 30, 41, 315, 0.31, 100],
+		]
+	)
+
+	actual_sorted = pdme.subspace_simulation.sort_array_of_dipoles_by_frequency(orig)
+	assert actual_sorted.tolist() == snapshot
+
+
+def test_sort_dipoleses_by_freq(snapshot):
+	sample_1 = numpy.array(
+		[
+			[1, 2, 3, 4, 5, 6, 7],
+			[100, 200, 300, 400, 500, 600, 0.07],
+			[10, 200, 30, 41, 315, 0.31, 100],
+		]
+	)
+
+	sample_2 = numpy.array(
+		[
+			[1, 1, 1, 1, 1, 1, 100],
+			[33, 33.3, 33, 3.3, 0.33, 0.3, 300],
+			[22, 22.2, 2.2, 222, 22, 2, 200],
+		]
+	)
+
+	original_samples = numpy.array([sample_1, sample_2])
+
+	actual_sorted = pdme.subspace_simulation.sort_array_of_dipoleses_by_frequency(
+		original_samples
+	)
+	assert actual_sorted.tolist() == snapshot

tests/subspace_simulation/test_stdevs.py

@@ -0,0 +1,38 @@
+import pytest
+import pdme.subspace_simulation
+
+
+def test_empty():
+	with pytest.raises(ValueError):
+		pdme.subspace_simulation.MCMCStandardDeviation([])
+
+
+def test_return_one(snapshot):
+	stdev = pdme.subspace_simulation.DipoleStandardDeviation(
+		1,
+		2,
+		3,
+		4,
+		5,
+		6,
+	)
+	stdevs = pdme.subspace_simulation.MCMCStandardDeviation([stdev])
+
+	assert stdevs[3] == snapshot
+	assert stdevs[3] == stdev
+
+
+def test_return_four(snapshot):
+	stdev_list = [
+		pdme.subspace_simulation.DipoleStandardDeviation(1, 2, 3, 4, 5, 6),
+		pdme.subspace_simulation.DipoleStandardDeviation(10, 20, 30, 40, 50, 60),
+		pdme.subspace_simulation.DipoleStandardDeviation(0.1, 0.2, 0.3, 0.4, 0.5, 0.6),
+	]
+	stdevs = pdme.subspace_simulation.MCMCStandardDeviation(stdev_list)
+
+	assert stdevs[0] == snapshot
+	assert stdevs[1] == snapshot
+	assert stdevs[2] == snapshot
+	assert stdevs[3] == snapshot
+	assert stdevs[4] == snapshot
+	assert stdevs[5] == snapshot

tests/test_pdme.py

@@ -3,5 +3,4 @@ import pdme
 
 
 def test_version():
-	assert __version__ == '0.0.1'
 	assert pdme.get_version() == __version__

tests/util/__snapshots__/test_fast_nonlocal_spectrum.ambr

@@ -0,0 +1,20 @@
+# serializer version: 1
+# name: test_arg
+  list([
+    list([
+      -0.0,
+      -0.0,
+      -0.0,
+    ]),
+    list([
+      3.141592653589793,
+      -0.0,
+      -0.0,
+    ]),
+    list([
+      -0.0,
+      -0.0,
+      3.141592653589793,
+    ]),
+  ])
+# ---

tests/util/__snapshots__/test_fast_nonlocal_spectrum_pairs.ambr

@@ -0,0 +1,25 @@
+# serializer version: 1
+# name: test_fast_nonlocal_calc_multidipole_phase_snapshot
+  list([
+    list([
+      -0.0,
+      -0.0,
+    ]),
+    list([
+      -0.0,
+      3.141592653589793,
+    ]),
+  ])
+# ---
+# name: test_fast_spin_qubit_frequency_tarucha_calc
+  list([
+    list([
+      1.1471308805198364,
+      0.042486328908142086,
+    ]),
+    list([
+      0.0008292459978410822,
+      0.0006214312613006971,
+    ]),
+  ])
+# ---

tests/util/test_fast_nonlocal_spectrum.py

@@ -0,0 +1,68 @@
+import numpy
+import pdme.util.fast_nonlocal_spectrum
+import pdme.measurement
+import logging
+import pytest
+
+
+def dipole_from_array(arr: numpy.ndarray) -> pdme.measurement.OscillatingDipole:
+	return pdme.measurement.OscillatingDipole(arr[0:3], arr[3:6], arr[6])
+
+
+def s_potential_from_arrays(
+	dipole_array: numpy.ndarray, dotf_pair_array: numpy.ndarray
+) -> float:
+	dipole = dipole_from_array(dipole_array)
+	r1 = dotf_pair_array[0][0:3]
+	f1 = dotf_pair_array[0][3]
+	r2 = dotf_pair_array[1][0:3]
+	return dipole.s_electric_potential_for_dot_pair(r1, r2, f1)
+
+
+def test_fast_nonlocal_calc():
+	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]
+
+	dipoles = 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 = numpy.array(
+		[
+			[s_potential_from_arrays(dipole_array, dot_pair) for dot_pair in dot_pairs]
+			for dipole_array in dipoles
+		]
+	)
+
+	# 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(dot_pairs, dipoles),
+		expected,
+		err_msg="nonlocal voltages at dot aren't as expected.",
+	)
+
+
+def test_fast_nonlocal_frequency_check():
+	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(dot_pairs, dipoles)
+
+
+def test_arg(snapshot):
+
+	test_input = numpy.array([[1, 2, 3], [-1, 1, 3], [3, 5, -1]])
+
+	actual_result = pdme.util.fast_nonlocal_spectrum.signarg(test_input)
+	assert actual_result.tolist() == snapshot

tests/util/test_fast_nonlocal_spectrum_pairs.py

@@ -0,0 +1,113 @@
+import numpy
+import pdme.util.fast_nonlocal_spectrum
+import pdme.measurement
+import logging
+import pytest
+
+_logger = logging.getLogger(__name__)
+
+
+def dipole_from_array(arr: numpy.ndarray) -> pdme.measurement.OscillatingDipole:
+	return pdme.measurement.OscillatingDipole(arr[0:3], arr[3:6], arr[6])
+
+
+def s_potential_from_arrays(
+	dipole_array: numpy.ndarray, dotf_pair_array: numpy.ndarray
+) -> float:
+	dipole = dipole_from_array(dipole_array)
+	r1 = dotf_pair_array[0][0:3]
+	f1 = dotf_pair_array[0][3]
+	r2 = dotf_pair_array[1][0:3]
+	return dipole.s_electric_potential_for_dot_pair(r1, r2, f1)
+
+
+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 = numpy.array(
+		[
+			[
+				sum(
+					[
+						s_potential_from_arrays(dipole_array, dot_pair)
+						for dipole_array in dipoles
+					]
+				)
+				for dot_pair in dot_pairs
+			]
+			for dipoles in dipoleses
+		]
+	)
+
+	# 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)
+
+
+def test_fast_nonlocal_calc_multidipole_phase_snapshot(snapshot):
+	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]]]
+	)
+
+	# 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)
+
+	actual_phases = pdme.util.fast_nonlocal_spectrum.signarg(
+		pdme.util.fast_nonlocal_spectrum.fast_s_nonlocal_dipoleses(dot_pairs, dipoleses)
+	)
+	assert actual_phases.tolist() == snapshot
+
+
+def test_fast_spin_qubit_frequency_tarucha_calc(snapshot):
+	d1 = [1, 2, 3, 0, 0, 0, 5]
+	d2 = [6, 7, 8, 5, 4, 3, 8]
+	dipoleses = numpy.array([[d1], [d2]])
+
+	dot_pairs = numpy.array(
+		[[[1, 0, 0, 1], [1, 0, 0, 1]], [[1, 0, 0, 1], [3, 0, 0, 1]]]
+	)
+
+	actual = (
+		pdme.util.fast_nonlocal_spectrum.fast_s_spin_qubit_tarucha_nonlocal_dipoleses(
+			dot_pairs, dipoleses
+		)
+	)
+
+	pdme.util.fast_nonlocal_spectrum._logger.setLevel(logging.DEBUG)
+	_logger.info(actual)
+	assert actual.tolist() == snapshot

tests/util/test_fast_v_calc.py

@@ -0,0 +1,219 @@
+import numpy
+import pdme.util.fast_v_calc
+
+import pdme.measurement
+
+
+def dipole_from_array(arr: numpy.ndarray) -> pdme.measurement.OscillatingDipole:
+	return pdme.measurement.OscillatingDipole(arr[0:3], arr[3:6], arr[6])
+
+
+def s_potential_from_arrays(
+	dipole_array: numpy.ndarray, dotf_array: numpy.ndarray
+) -> float:
+	dipole = dipole_from_array(dipole_array)
+	r = dotf_array[0:3]
+	f = dotf_array[3]
+	return dipole.s_electric_potential_at_position(r, f)
+
+
+def s_electric_field_x_from_arrays(
+	dipole_array: numpy.ndarray, dotf_array: numpy.ndarray
+) -> float:
+	dipole = dipole_from_array(dipole_array)
+	r = dotf_array[0:3]
+	f = dotf_array[3]
+	return dipole.s_electric_fieldx_at_position(r, f)
+
+
+def test_fast_v_calc():
+	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_11 = s_potential_from_arrays(dipoles[0], dot_inputs[0])
+	expected_12 = s_potential_from_arrays(dipoles[1], dot_inputs[0])
+	expected_21 = s_potential_from_arrays(dipoles[0], dot_inputs[1])
+	expected_22 = s_potential_from_arrays(dipoles[1], dot_inputs[1])
+
+	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.",
+	)
+
+
+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_11 = s_potential_from_arrays(dipoles[0][0], dot_inputs[0])
+	expected_12 = s_potential_from_arrays(dipoles[0][1], dot_inputs[0])
+	expected_21 = s_potential_from_arrays(dipoles[0][0], dot_inputs[1])
+	expected_22 = s_potential_from_arrays(dipoles[0][1], dot_inputs[1])
+
+	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():
+
+	dipoleses = 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 = [
+		[
+			sum(
+				[
+					s_potential_from_arrays(dipole_array, dot_input)
+					for dipole_array in dipole_config
+				]
+			)
+			for dot_input in dot_inputs
+		]
+		for dipole_config in dipoleses
+	]
+
+	numpy.testing.assert_allclose(
+		pdme.util.fast_v_calc.fast_vs_for_dipoleses(dot_inputs, dipoleses),
+		expected,
+		err_msg="Voltages at dot aren't as expected for multidipole calc.",
+	)
+
+
+def test_fast_electric_field_x_calc():
+	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_11 = s_electric_field_x_from_arrays(dipoles[0], dot_inputs[0])
+	expected_12 = s_electric_field_x_from_arrays(dipoles[1], dot_inputs[0])
+	expected_21 = s_electric_field_x_from_arrays(dipoles[0], dot_inputs[1])
+	expected_22 = s_electric_field_x_from_arrays(dipoles[1], dot_inputs[1])
+
+	expected = numpy.array([[expected_11, expected_21], [expected_12, expected_22]])
+
+	numpy.testing.assert_allclose(
+		pdme.util.fast_v_calc.fast_efieldxs_for_dipoles(dot_inputs, dipoles),
+		expected,
+		err_msg="E x fast calc at dot aren't as expected.",
+	)
+
+
+def test_fast_electric_field_x_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_11 = s_electric_field_x_from_arrays(dipoles[0][0], dot_inputs[0])
+	expected_12 = s_electric_field_x_from_arrays(dipoles[0][1], dot_inputs[0])
+	expected_21 = s_electric_field_x_from_arrays(dipoles[0][0], dot_inputs[1])
+	expected_22 = s_electric_field_x_from_arrays(dipoles[0][1], dot_inputs[1])
+
+	expected = numpy.array([[expected_11 + expected_12, expected_21 + expected_22]])
+
+	numpy.testing.assert_allclose(
+		pdme.util.fast_v_calc.fast_efieldxs_for_dipoleses(dot_inputs, dipoles),
+		expected,
+		err_msg="E x fast calc at dot aren't as expected for multidipole calc.",
+	)
+
+
+def test_fast_electric_field_x_calc_big_multidipole():
+
+	dipoleses = 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 = [
+		[
+			sum(
+				[
+					s_electric_field_x_from_arrays(dipole_array, dot_input)
+					for dipole_array in dipole_config
+				]
+			)
+			for dot_input in dot_inputs
+		]
+		for dipole_config in dipoleses
+	]
+
+	numpy.testing.assert_allclose(
+		pdme.util.fast_v_calc.fast_efieldxs_for_dipoleses(dot_inputs, dipoleses),
+		expected,
+		err_msg="E x fast calc at dot aren't as expected for multidipole calc.",
+	)
+
+
+def test_between():
+	low = numpy.array([1, 2, 3])
+	high = numpy.array([6, 7, 8])
+
+	# 	  FALSE		FALSE	   TRUE
+	a = [[0, 1, 2], [0, 9, 5], [4, 5, 6]]
+
+	actual = pdme.util.fast_v_calc.between(a, low, high)
+	expected = numpy.array([False, False, True])
+
+	numpy.testing.assert_array_equal(actual, expected, err_msg="Between calc wrong")