Compare commits

..

No commits in common. "master" and "0.5.4" have entirely different histories.

79 changed files with 6423 additions and 5740 deletions

View File

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

View File

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

4
.gitignore vendored
View File

@ -114,10 +114,6 @@ ENV/
env.bak/ env.bak/
venv.bak/ venv.bak/
#direnv
.envrc
.direnv
# Spyder project settings # Spyder project settings
.spyderproject .spyderproject
.spyproject .spyproject

View File

@ -2,226 +2,6 @@
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. 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) ### [0.5.4](https://gitea.deepak.science:2222/physics/pdme/compare/0.5.3...0.5.4) (2022-03-06)

View File

@ -5,7 +5,7 @@
[![Jenkins](https://img.shields.io/jenkins/build?jobUrl=https%3A%2F%2Fjenkins.deepak.science%2Fjob%2Fgitea-physics%2Fjob%2Fpdme%2Fjob%2Fmaster&style=flat-square)](https://jenkins.deepak.science/job/gitea-physics/job/pdme/job/master/) [![Jenkins](https://img.shields.io/jenkins/build?jobUrl=https%3A%2F%2Fjenkins.deepak.science%2Fjob%2Fgitea-physics%2Fjob%2Fpdme%2Fjob%2Fmaster&style=flat-square)](https://jenkins.deepak.science/job/gitea-physics/job/pdme/job/master/)
![Jenkins tests](https://img.shields.io/jenkins/tests?compact_message&jobUrl=https%3A%2F%2Fjenkins.deepak.science%2Fjob%2Fgitea-physics%2Fjob%2Fpdme%2Fjob%2Fmaster%2F&style=flat-square) ![Jenkins tests](https://img.shields.io/jenkins/tests?compact_message&jobUrl=https%3A%2F%2Fjenkins.deepak.science%2Fjob%2Fgitea-physics%2Fjob%2Fpdme%2Fjob%2Fmaster%2F&style=flat-square)
![Jenkins Coverage](https://img.shields.io/jenkins/coverage/cobertura?jobUrl=https%3A%2F%2Fjenkins.deepak.science%2Fjob%2Fgitea-physics%2Fjob%2Fpdme%2Fjob%2Fmaster%2F&style=flat-square) ![Jenkins Coverage](https://img.shields.io/jenkins/coverage/cobertura?jobUrl=https%3A%2F%2Fjenkins.deepak.science%2Fjob%2Fgitea-physics%2Fjob%2Fpdme%2Fjob%2Fmaster%2F&style=flat-square)
![Maintenance](https://img.shields.io/maintenance/yes/2024?style=flat-square) ![Maintenance](https://img.shields.io/maintenance/yes/2022?style=flat-square)
This repo has library code for evaluating dipole models. This repo has library code for evaluating dipole models.

5029
diagnosis1.nb Normal file

File diff suppressed because it is too large Load Diff

33
do.sh Normal file
View File

@ -0,0 +1,33 @@
#!/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
}
release() {
./scripts/release.sh
}
all() {
build && test
}
"$@" # <- execute the task
[ "$#" -gt 0 ] || printf "Usage:\n\t./do.sh %s\n" "($(compgen -A function | grep '^[^_]' | paste -sd '|' -))"

95
flake.lock generated
View File

@ -1,95 +0,0 @@
{
"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
}

View File

@ -1,55 +0,0 @@
{
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 ];
};
}));
}

View File

@ -1,64 +0,0 @@
# 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

View File

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

View File

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

View File

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

View File

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

View File

@ -4,22 +4,19 @@ import numpy.typing
@dataclass @dataclass
class DotMeasurement: class DotMeasurement():
""" '''
Representation of a dot measuring oscillating dipoles. Representation of a dot measuring oscillating dipoles.
Parameters Parameters
---------- ----------
v : float v : float
The voltage measured at the dot. The voltage measured at the dot.
r : numpy.ndarray r : numpy.ndarray
The position of the dot. The position of the dot.
f : float f : float
The measurement frequency. The measurement frequency.
""" '''
v: float v: float
r: numpy.ndarray r: numpy.ndarray
f: float f: float
@ -29,25 +26,21 @@ class DotMeasurement:
@dataclass @dataclass
class DotRangeMeasurement: class DotRangeMeasurement():
""" '''
Representation of a dot measuring oscillating dipoles. Representation of a dot measuring oscillating dipoles.
Parameters Parameters
---------- ----------
v_low : float v_low : float
The lower range of voltage measured at the dot. The lower range of voltage measured at the dot.
v_high : float v_high : float
The upper range of voltage measured at the dot. The upper range of voltage measured at the dot.
r : numpy.ndarray r : numpy.ndarray
The position of the dot. The position of the dot.
f : float f : float
The measurement frequency. The measurement frequency.
""" '''
v_low: float v_low: float
v_high: float v_high: float
r: numpy.ndarray r: numpy.ndarray
@ -55,5 +48,3 @@ class DotRangeMeasurement:
def __post_init__(self) -> None: def __post_init__(self) -> None:
self.r = numpy.array(self.r) 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

View File

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

View File

@ -1,36 +0,0 @@
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))

View File

@ -1,267 +1,98 @@
from dataclasses import dataclass from dataclasses import dataclass
import numpy import numpy
import numpy.typing import numpy.typing
from typing import Sequence, List from typing import Sequence, List, Tuple
from pdme.measurement.dot_measure import DotMeasurement, DotRangeMeasurement from pdme.measurement.dot_measure import DotMeasurement, DotRangeMeasurement
from pdme.measurement.dot_pair_measure import (
DotPairMeasurement,
DotPairRangeMeasurement, DotInput = Tuple[numpy.typing.ArrayLike, float]
)
import pdme.calculations
from pdme.measurement.input_types import DotInput, DotPairInput
@dataclass @dataclass
class OscillatingDipole: class OscillatingDipole():
""" '''
Representation of an oscillating dipole, either known or guessed. Representation of an oscillating dipole, either known or guessed.
Parameters Parameters
---------- ----------
p : numpy.ndarray p : numpy.ndarray
The oscillating dipole moment, with overall sign arbitrary. The oscillating dipole moment, with overall sign arbitrary.
s : numpy.ndarray s : numpy.ndarray
The position of the dipole. The position of the dipole.
w : float w : float
The oscillation frequency. The oscillation frequency.
""" '''
p: numpy.ndarray p: numpy.ndarray
s: numpy.ndarray s: numpy.ndarray
w: float w: float
def __post_init__(self) -> None: def __post_init__(self) -> None:
""" '''
Coerce the inputs into numpy arrays. Coerce the inputs into numpy arrays.
""" '''
self.p = numpy.array(self.p) self.p = numpy.array(self.p)
self.s = numpy.array(self.s) self.s = numpy.array(self.s)
def _alpha_electric_potential(self, r: numpy.ndarray) -> float: def s_at_position(self, r: numpy.ndarray, f: float) -> float:
""" '''
Returns the electric potential of this dipole at position r. Returns the noise potential at a point r, at some frequency f.
"""
return pdme.calculations.electric_potential(self.p, self.s, r)
def _alpha_electric_field(self, r: numpy.ndarray) -> numpy.ndarray: Parameters
""" ----------
Returns the electric field of this dipole at position r. r : numpy.ndarray
""" The position of the dot.
return pdme.calculations.electric_field(self.p, self.s, r) f : float
The dot frequency to sample.
'''
return (self._alpha(r))**2 * self._b(f)
def _alpha(self, r: numpy.ndarray) -> float:
diff = r - self.s
return self.p.dot(diff) / (numpy.linalg.norm(diff)**3)
def _b(self, f: float) -> float: def _b(self, f: float) -> float:
return pdme.calculations.telegraph_beta(f, self.w) return (1 / numpy.pi) * (self.w / (f**2 + self.w**2))
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: 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_range_measurements_low_high_arrays(dot_range_measurements: Sequence[DotRangeMeasurement]) -> 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))
class OscillatingDipoleArrangement():
'''
A collection of oscillating dipoles, which we are interested in being able to characterise. A collection of oscillating dipoles, which we are interested in being able to characterise.
Parameters Parameters
-------- --------
dipoles : Sequence[OscillatingDipole] dipoles : Sequence[OscillatingDipole]
""" '''
def __init__(self, dipoles: Sequence[OscillatingDipole]): def __init__(self, dipoles: Sequence[OscillatingDipole]):
self.dipoles = dipoles self.dipoles = dipoles
def get_potential_dot_measurement(self, dot_input: DotInput) -> DotMeasurement: def get_dot_measurement(self, dot_input: DotInput) -> DotMeasurement:
r = numpy.array(dot_input[0]) r = numpy.array(dot_input[0])
f = dot_input[1] f = dot_input[1]
return DotMeasurement( return DotMeasurement(sum([dipole.s_at_position(r, f) for dipole in self.dipoles]), r, f)
sum(
[
dipole.s_electric_potential_at_position(r, f)
for dipole in self.dipoles
]
),
r,
f,
)
def get_potential_dot_pair_measurement( def get_dot_measurements(self, dot_inputs: Sequence[DotInput]) -> List[DotMeasurement]:
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. For a series of points, each with three coordinates and a frequency, return a list of the corresponding DotMeasurements.
""" '''
return [ return [self.get_dot_measurement(dot_input) for dot_input in dot_inputs]
self.get_potential_dot_measurement(dot_input) for dot_input in dot_inputs
]
def get_potential_dot_pair_measurements( def get_percent_range_dot_measurement(self, dot_input: DotInput, low_percent: float, high_percent: float) -> DotRangeMeasurement:
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]) r = numpy.array(dot_input[0])
f = dot_input[1] f = dot_input[1]
return DotRangeMeasurement( return DotRangeMeasurement(low_percent * sum([dipole.s_at_position(r, f) for dipole in self.dipoles]), high_percent * sum([dipole.s_at_position(r, f) for dipole in self.dipoles]), r, f)
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( def get_percent_range_dot_measurements(self, dot_inputs: Sequence[DotInput], low_percent: float, high_percent: float) -> List[DotRangeMeasurement]:
self, dot_inputs: Sequence[DotInput], low_percent: float, high_percent: float '''
) -> List[DotRangeMeasurement]: For a series 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 DotRangeMeasurements.
""" '''
For a series of pairs of points, each with three coordinates and a frequency, and also a lower error range and upper error range, return a list of the corresponding DotPairRangeMeasurements. return [self.get_percent_range_dot_measurement(dot_input, low_percent, high_percent) for dot_input in dot_inputs]
"""
return [
self.get_percent_range_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])

View File

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

View File

@ -1,30 +1,7 @@
from pdme.model.model import DipoleModel from pdme.model.model import Model, Discretisation
from pdme.model.fixed_magnitude_model import SingleDipoleFixedMagnitudeModel from pdme.model.fixed_z_plane_model import FixedZPlaneModel
from pdme.model.multidipole_fixed_magnitude_model import ( from pdme.model.unrestricted_model import UnrestrictedModel
MultipleDipoleFixedMagnitudeModel, from pdme.model.fixed_dipole_model import FixedDipoleModel
) from pdme.model.fixed_magnitude_model import FixedMagnitudeModel
from pdme.model.random_count_multidipole_fixed_magnitude_model import (
RandomCountMultipleDipoleFixedMagnitudeModel,
)
from pdme.model.log_spaced_random_choice_model import ( __all__ = ["Model", "Discretisation", "FixedZPlaneModel", "UnrestrictedModel", "FixedDipoleModel", "FixedMagnitudeModel"]
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",
]

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

@ -1,163 +1,119 @@
import numpy import numpy
import numpy.random import scipy.optimize
from pdme.measurement import ( from typing import Callable, Sequence, Tuple, List
OscillatingDipoleArrangement, from pdme.measurement import DotMeasurement, OscillatingDipoleArrangement, OscillatingDipole
) import pdme.util
import logging import logging
import pdme.subspace_simulation
from typing import List, Tuple, Optional
_logger = logging.getLogger(__name__) _logger = logging.getLogger(__name__)
class DipoleModel: class Model():
""" """
Interface for models based on dipoles. Interface for models.
Some concepts are kept specific for dipole-based models, even though other types of models could be useful later on.
""" """
def get_dipoles( def point_length(self) -> int:
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 raise NotImplementedError
def get_monte_carlo_dipole_inputs( def n(self) -> int:
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 raise NotImplementedError
def markov_chain_monte_carlo_proposal( def v_for_point_at_dot(self, dot: DotMeasurement, pt: numpy.ndarray) -> float:
self,
dipole: numpy.ndarray,
stdev: pdme.subspace_simulation.DipoleStandardDeviation,
rng_arg: Optional[numpy.random.Generator] = None,
) -> numpy.ndarray:
raise NotImplementedError raise NotImplementedError
def get_mcmc_chain( def get_dipoles(self, frequency: float) -> OscillatingDipoleArrangement:
self, raise NotImplementedError
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. def get_n_single_dipoles(self, n: int, max_frequency: float) -> numpy.ndarray:
Note that for our adaptive stdevs to work, there's an unwritten contract that we sort each dipole in the state by frequency (increasing). raise NotImplementedError
The seed is a list of dipoles, and each chain state is a list of dipoles as well. def solution_single_dipole(self, pt: numpy.ndarray) -> OscillatingDipole:
raise NotImplementedError
initial_cost is a performance guy that lets you pre-populate the initial cost used to define the condition. def solution_as_dipoles(self, pts: numpy.ndarray) -> List[OscillatingDipole]:
Probably premature optimisation. pt_length = self.point_length()
chunked_pts = [pts[i: i + pt_length] for i in range(0, len(pts), pt_length)]
return [self.solution_single_dipole(pt) for pt in chunked_pts]
Returns a chain of [ (cost: float, state: dipole_ndarray ) ] format. def cost_for_dot(self, dot: DotMeasurement, pts: numpy.ndarray) -> float:
""" # creates numpy.ndarrays in groups of self.point_length().
_logger.debug( # Will throw problems for irregular points, but that's okay for now.
f"Starting Markov Chain Monte Carlo with seed: {seed} for chain length {chain_length} and provided stdevs {stdevs}" pt_length = self.point_length()
) chunked_pts = [pts[i: i + pt_length] for i in range(0, len(pts), pt_length)]
chain: List[Tuple[float, numpy.ndarray]] = [] return sum(self.v_for_point_at_dot(dot, pt) for pt in chunked_pts) - dot.v
if initial_cost is None:
current_cost = cost_function(numpy.array([seed])) def costs(self, dots: Sequence[DotMeasurement]) -> Callable[[numpy.ndarray], numpy.ndarray]:
'''
Returns a function that returns the cost for the given list of DotMeasurements for a particular model-dependent phase space point.
Default implementation assumes a single dot cost from which to build the list.
Parameters
----------
dots: A list of dot measurements to use to find the cost functions.
Returns
----------
Returns the model's cost function.
'''
_logger.debug(f"Constructing costs for dots: {dots}")
def costs_to_return(pts: numpy.ndarray) -> numpy.ndarray:
return numpy.array([self.cost_for_dot(dot, pts) for dot in dots])
return costs_to_return
def jac_for_point_at_dot(self, dot: DotMeasurement, pt: numpy.ndarray) -> numpy.ndarray:
raise NotImplementedError
def jac_for_dot(self, dot: DotMeasurement, pts: numpy.ndarray) -> numpy.ndarray:
# creates numpy.ndarrays in groups of self.point_length().
# Will throw problems for irregular points, but that's okay for now.
pt_length = self.point_length()
chunked_pts = [pts[i: i + pt_length] for i in range(0, len(pts), pt_length)]
return numpy.append([], [self.jac_for_point_at_dot(dot, pt) for pt in chunked_pts])
def jac(self, dots: Sequence[DotMeasurement]) -> Callable[[numpy.ndarray], numpy.ndarray]:
'''
Returns a function that returns the cost function's Jacobian for the given list of DotMeasurements for a particular model-dependent phase space point.
Default implementation assumes a single dot jacobian from which to build the list.
Parameters
----------
dots: A list of dot measurements to use to find the cost functions and their Jacobian.
Returns
----------
Returns the model's cost function's Jacobian.
'''
def jac_to_return(pts: numpy.ndarray) -> numpy.ndarray:
return numpy.array([self.jac_for_dot(dot, pts) for dot in dots])
return jac_to_return
def solve(self, dots: Sequence[DotMeasurement], initial_pt: numpy.ndarray = None, bounds=(-numpy.inf, numpy.inf)) -> scipy.optimize.OptimizeResult:
if initial_pt is None:
initial = numpy.tile(.1, self.n() * self.point_length())
else: else:
current_cost = initial_cost if len(initial_pt) != self.point_length():
current = seed raise ValueError(f"The initial point {initial_pt} does not have the model's expected length: {self.point_length()}")
for i in range(chain_length): initial = numpy.tile(initial_pt, self.n())
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) result = scipy.optimize.least_squares(self.costs(dots), initial, jac=self.jac(dots), ftol=1e-15, gtol=3e-16, xtol=None, bounds=bounds)
dips_array = pdme.subspace_simulation.sort_array_of_dipoles_by_frequency( result.normalised_x = pdme.util.normalise_point_list(result.x, self.point_length())
dips return result
)
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. class Discretisation():
Note that for our adaptive stdevs to work, there's an unwritten contract that we sort each dipole in the state by frequency (increasing). def bounds(self, index: Tuple[float, ...]) -> Tuple:
raise NotImplementedError
The seed is a list of dipoles, and each chain state is a list of dipoles as well. def all_indices(self) -> numpy.ndindex:
raise NotImplementedError
initial_cost is a performance guy that lets you pre-populate the initial cost used to define the condition. def solve_for_index(self, dots: Sequence[DotMeasurement], index: Tuple) -> scipy.optimize.OptimizeResult:
Probably premature optimisation. raise NotImplementedError
Chain has type of [ (cost: float, state: dipole_ndarray ) ] format, def get_model(self) -> Model:
returning (repeat_count, chain) to keep track of number of repeats raise NotImplementedError
"""
_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)

View File

@ -1,104 +0,0 @@
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)

View File

@ -1,119 +0,0 @@
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)

View File

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

View File

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

View File

@ -1,30 +0,0 @@
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)

View File

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

View File

@ -1,273 +0,0 @@
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)

View File

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

32
pdme/util/normal_form.py Normal file
View File

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

628
poetry.lock generated
View File

@ -1,65 +1,24 @@
[[package]] [[package]]
name = "appnope" name = "atomicwrites"
version = "0.1.2" version = "1.4.0"
description = "Disable App Nap on macOS >= 10.9" description = "Atomic file writes."
category = "dev" category = "dev"
optional = false optional = false
python-versions = "*" python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*"
[[package]] [[package]]
name = "asttokens" name = "attrs"
version = "2.0.5" version = "21.4.0"
description = "Annotate AST trees with source code positions" description = "Classes Without Boilerplate"
category = "dev" category = "dev"
optional = false optional = false
python-versions = "*" python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.*"
[package.dependencies]
six = "*"
[package.extras] [package.extras]
test = ["astroid", "pytest"] 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"]
[[package]] tests = ["coverage[toml] (>=5.0.2)", "hypothesis", "pympler", "pytest (>=4.3.0)", "six", "mypy", "pytest-mypy-plugins", "zope.interface", "cloudpickle"]
name = "backcall" tests_no_zope = ["coverage[toml] (>=5.0.2)", "hypothesis", "pympler", "pytest (>=4.3.0)", "six", "mypy", "pytest-mypy-plugins", "cloudpickle"]
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]] [[package]]
name = "colorama" name = "colorama"
@ -69,67 +28,29 @@ category = "dev"
optional = false optional = false
python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.*" 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]] [[package]]
name = "coverage" name = "coverage"
version = "7.2.7" version = "6.2"
description = "Code coverage measurement for Python" description = "Code coverage measurement for Python"
category = "dev" category = "dev"
optional = false optional = false
python-versions = ">=3.7" python-versions = ">=3.6"
[package.dependencies]
tomli = {version = "*", optional = true, markers = "python_full_version <= \"3.11.0a6\" and extra == \"toml\""}
[package.extras] [package.extras]
toml = ["tomli"] 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]] [[package]]
name = "flake8" name = "flake8"
version = "4.0.1" version = "3.9.2"
description = "the modular source code checker: pep8 pyflakes and co" description = "the modular source code checker: pep8 pyflakes and co"
category = "dev" category = "dev"
optional = false optional = false
python-versions = ">=3.6" python-versions = "!=3.0.*,!=3.1.*,!=3.2.*,!=3.3.*,!=3.4.*,>=2.7"
[package.dependencies] [package.dependencies]
mccabe = ">=0.6.0,<0.7.0" mccabe = ">=0.6.0,<0.7.0"
pycodestyle = ">=2.8.0,<2.9.0" pycodestyle = ">=2.7.0,<2.8.0"
pyflakes = ">=2.4.0,<2.5.0" pyflakes = ">=2.3.0,<2.4.0"
[[package]] [[package]]
name = "iniconfig" name = "iniconfig"
@ -139,67 +60,6 @@ category = "dev"
optional = false optional = false
python-versions = "*" 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]] [[package]]
name = "mccabe" name = "mccabe"
version = "0.6.1" version = "0.6.1"
@ -210,21 +70,20 @@ python-versions = "*"
[[package]] [[package]]
name = "mypy" name = "mypy"
version = "0.961" version = "0.910"
description = "Optional static typing for Python" description = "Optional static typing for Python"
category = "dev" category = "dev"
optional = false optional = false
python-versions = ">=3.6" python-versions = ">=3.5"
[package.dependencies] [package.dependencies]
mypy-extensions = ">=0.4.3" mypy-extensions = ">=0.4.3,<0.5.0"
tomli = {version = ">=1.1.0", markers = "python_version < \"3.11\""} toml = "*"
typing-extensions = ">=3.10" typing-extensions = ">=3.7.4"
[package.extras] [package.extras]
dmypy = ["psutil (>=4.0)"] dmypy = ["psutil (>=4.0)"]
python2 = ["typed-ast (>=1.4.0,<2)"] python2 = ["typed-ast (>=1.4.0,<1.5.0)"]
reports = ["lxml"]
[[package]] [[package]]
name = "mypy-extensions" name = "mypy-extensions"
@ -236,7 +95,7 @@ python-versions = "*"
[[package]] [[package]]
name = "numpy" name = "numpy"
version = "1.22.3" version = "1.22.0"
description = "NumPy is the fundamental package for array computing with Python." description = "NumPy is the fundamental package for array computing with Python."
category = "main" category = "main"
optional = false optional = false
@ -253,57 +112,6 @@ python-versions = ">=3.6"
[package.dependencies] [package.dependencies]
pyparsing = ">=2.0.2,<3.0.5 || >3.0.5" 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]] [[package]]
name = "pluggy" name = "pluggy"
version = "1.0.0" version = "1.0.0"
@ -317,58 +125,28 @@ dev = ["pre-commit", "tox"]
testing = ["pytest", "pytest-benchmark"] testing = ["pytest", "pytest-benchmark"]
[[package]] [[package]]
name = "prompt-toolkit" name = "py"
version = "3.0.28" version = "1.11.0"
description = "Library for building powerful interactive command lines in Python" description = "library with cross-python path, ini-parsing, io, code, log facilities"
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" category = "dev"
optional = false optional = false
python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.*" python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.*"
[[package]] [[package]]
name = "pyflakes" name = "pycodestyle"
version = "2.4.0" version = "2.7.0"
description = "passive checker of Python programs" description = "Python style guide checker"
category = "dev" category = "dev"
optional = false optional = false
python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*" python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*"
[[package]] [[package]]
name = "pygments" name = "pyflakes"
version = "2.11.2" version = "2.3.1"
description = "Pygments is a syntax highlighting package written in Python." description = "passive checker of Python programs"
category = "dev" category = "dev"
optional = false optional = false
python-versions = ">=3.5" python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*"
[[package]] [[package]]
name = "pyparsing" name = "pyparsing"
@ -383,108 +161,59 @@ diagrams = ["jinja2", "railroad-diagrams"]
[[package]] [[package]]
name = "pytest" name = "pytest"
version = "7.4.0" version = "6.2.5"
description = "pytest: simple powerful testing with Python" description = "pytest: simple powerful testing with Python"
category = "dev" category = "dev"
optional = false optional = false
python-versions = ">=3.7" python-versions = ">=3.6"
[package.dependencies] [package.dependencies]
atomicwrites = {version = ">=1.0", markers = "sys_platform == \"win32\""}
attrs = ">=19.2.0"
colorama = {version = "*", markers = "sys_platform == \"win32\""} colorama = {version = "*", markers = "sys_platform == \"win32\""}
exceptiongroup = {version = ">=1.0.0rc8", markers = "python_version < \"3.11\""}
iniconfig = "*" iniconfig = "*"
packaging = "*" packaging = "*"
pluggy = ">=0.12,<2.0" pluggy = ">=0.12,<2.0"
tomli = {version = ">=1.0.0", markers = "python_version < \"3.11\""} py = ">=1.8.2"
toml = "*"
[package.extras] [package.extras]
testing = ["argcomplete", "attrs (>=19.2.0)", "hypothesis (>=3.56)", "mock", "nose", "pygments (>=2.7.2)", "requests", "setuptools", "xmlschema"] testing = ["argcomplete", "hypothesis (>=3.56)", "mock", "nose", "requests", "xmlschema"]
[[package]] [[package]]
name = "pytest-cov" name = "pytest-cov"
version = "4.1.0" version = "2.12.1"
description = "Pytest plugin for measuring coverage." description = "Pytest plugin for measuring coverage."
category = "dev" category = "dev"
optional = false optional = false
python-versions = ">=3.7" python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.*"
[package.dependencies] [package.dependencies]
coverage = {version = ">=5.2.1", extras = ["toml"]} coverage = ">=5.2.1"
pytest = ">=4.6" pytest = ">=4.6"
toml = "*"
[package.extras] [package.extras]
testing = ["fields", "hunter", "process-tests", "six", "pytest-xdist", "virtualenv"] testing = ["fields", "hunter", "process-tests", "six", "pytest-xdist", "virtualenv"]
[[package]] [[package]]
name = "scipy" name = "scipy"
version = "1.10.0" version = "1.5.4"
description = "Fundamental algorithms for scientific computing in Python" description = "SciPy: Scientific Library for Python"
category = "main" category = "main"
optional = false optional = false
python-versions = "<3.12,>=3.8" python-versions = ">=3.6"
[package.dependencies] [package.dependencies]
numpy = ">=1.19.5,<1.27.0" numpy = ">=1.14.5"
[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]] [[package]]
name = "six" name = "toml"
version = "1.16.0" version = "0.10.2"
description = "Python 2 and 3 compatibility utilities" description = "Python Library for Tom's Obvious, Minimal Language"
category = "dev" category = "dev"
optional = false optional = false
python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*" python-versions = ">=2.6, !=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]] [[package]]
name = "typing-extensions" name = "typing-extensions"
@ -494,61 +223,202 @@ category = "dev"
optional = false optional = false
python-versions = ">=3.6" 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] [metadata]
lock-version = "1.1" lock-version = "1.1"
python-versions = ">=3.8.1,<3.10" python-versions = "^3.8,<3.10"
content-hash = "b5275c33449e8f85acbf9c0f6dfe1ec4e7296adfa16360d782b33534e1223638" content-hash = "1f00423d851dd44448cfd16f1b00f1672664d766f8add03a4c927f01f39cea80"
[metadata.files] [metadata.files]
appnope = [] atomicwrites = [
asttokens = [] {file = "atomicwrites-1.4.0-py2.py3-none-any.whl", hash = "sha256:6d1784dea7c0c8d4a5172b6c620f40b6e4cbfdf96d783691f2e1302a7b88e197"},
backcall = [] {file = "atomicwrites-1.4.0.tar.gz", hash = "sha256:ae70396ad1a434f9c7046fd2dd196fc04b12f9e91ffb859164193be8b6168a7a"},
black = [] ]
click = [] attrs = [
colorama = [] {file = "attrs-21.4.0-py2.py3-none-any.whl", hash = "sha256:2d27e3784d7a565d36ab851fe94887c5eccd6a463168875832a1be79c82828b4"},
colored = [] {file = "attrs-21.4.0.tar.gz", hash = "sha256:626ba8234211db98e869df76230a137c4c40a12d72445c45d5f5b716f076e2fd"},
coverage = [] ]
decorator = [] colorama = [
exceptiongroup = [] {file = "colorama-0.4.4-py2.py3-none-any.whl", hash = "sha256:9f47eda37229f68eee03b24b9748937c7dc3868f906e8ba69fbcbdd3bc5dc3e2"},
executing = [] {file = "colorama-0.4.4.tar.gz", hash = "sha256:5941b2b48a20143d2267e95b1c2a7603ce057ee39fd88e7329b0c292aa16869b"},
flake8 = [] ]
iniconfig = [] coverage = [
ipython = [] {file = "coverage-6.2-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:6dbc1536e105adda7a6312c778f15aaabe583b0e9a0b0a324990334fd458c94b"},
jedi = [] {file = "coverage-6.2-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:174cf9b4bef0db2e8244f82059a5a72bd47e1d40e71c68ab055425172b16b7d0"},
matplotlib-inline = [] {file = "coverage-6.2-cp310-cp310-manylinux_2_5_i686.manylinux1_i686.manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:92b8c845527eae547a2a6617d336adc56394050c3ed8a6918683646328fbb6da"},
mccabe = [] {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"},
mypy = [] {file = "coverage-6.2-cp310-cp310-musllinux_1_1_aarch64.whl", hash = "sha256:d5d2033d5db1d58ae2d62f095e1aefb6988af65b4b12cb8987af409587cc0739"},
mypy-extensions = [] {file = "coverage-6.2-cp310-cp310-musllinux_1_1_i686.whl", hash = "sha256:3feac4084291642165c3a0d9eaebedf19ffa505016c4d3db15bfe235718d4971"},
numpy = [] {file = "coverage-6.2-cp310-cp310-musllinux_1_1_x86_64.whl", hash = "sha256:276651978c94a8c5672ea60a2656e95a3cce2a3f31e9fb2d5ebd4c215d095840"},
packaging = [] {file = "coverage-6.2-cp310-cp310-win32.whl", hash = "sha256:f506af4f27def639ba45789fa6fde45f9a217da0be05f8910458e4557eed020c"},
parso = [] {file = "coverage-6.2-cp310-cp310-win_amd64.whl", hash = "sha256:3f7c17209eef285c86f819ff04a6d4cbee9b33ef05cbcaae4c0b4e8e06b3ec8f"},
pathspec = [] {file = "coverage-6.2-cp311-cp311-macosx_10_14_x86_64.whl", hash = "sha256:13362889b2d46e8d9f97c421539c97c963e34031ab0cb89e8ca83a10cc71ac76"},
pexpect = [] {file = "coverage-6.2-cp311-cp311-manylinux_2_5_x86_64.manylinux1_x86_64.whl", hash = "sha256:22e60a3ca5acba37d1d4a2ee66e051f5b0e1b9ac950b5b0cf4aa5366eda41d47"},
pickleshare = [] {file = "coverage-6.2-cp311-cp311-win_amd64.whl", hash = "sha256:b637c57fdb8be84e91fac60d9325a66a5981f8086c954ea2772efe28425eaf64"},
platformdirs = [] {file = "coverage-6.2-cp36-cp36m-macosx_10_9_x86_64.whl", hash = "sha256:f467bbb837691ab5a8ca359199d3429a11a01e6dfb3d9dcc676dc035ca93c0a9"},
pluggy = [] {file = "coverage-6.2-cp36-cp36m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:2641f803ee9f95b1f387f3e8f3bf28d83d9b69a39e9911e5bfee832bea75240d"},
prompt-toolkit = [] {file = "coverage-6.2-cp36-cp36m-manylinux_2_5_i686.manylinux1_i686.manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:1219d760ccfafc03c0822ae2e06e3b1248a8e6d1a70928966bafc6838d3c9e48"},
ptyprocess = [] {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"},
pure-eval = [] {file = "coverage-6.2-cp36-cp36m-musllinux_1_1_aarch64.whl", hash = "sha256:8e2c35a4c1f269704e90888e56f794e2d9c0262fb0c1b1c8c4ee44d9b9e77b5d"},
pycodestyle = [] {file = "coverage-6.2-cp36-cp36m-musllinux_1_1_i686.whl", hash = "sha256:5d6b09c972ce9200264c35a1d53d43ca55ef61836d9ec60f0d44273a31aa9f17"},
pyflakes = [] {file = "coverage-6.2-cp36-cp36m-musllinux_1_1_x86_64.whl", hash = "sha256:e3db840a4dee542e37e09f30859f1612da90e1c5239a6a2498c473183a50e781"},
pygments = [] {file = "coverage-6.2-cp36-cp36m-win32.whl", hash = "sha256:4e547122ca2d244f7c090fe3f4b5a5861255ff66b7ab6d98f44a0222aaf8671a"},
pyparsing = [] {file = "coverage-6.2-cp36-cp36m-win_amd64.whl", hash = "sha256:01774a2c2c729619760320270e42cd9e797427ecfddd32c2a7b639cdc481f3c0"},
pytest = [] {file = "coverage-6.2-cp37-cp37m-macosx_10_9_x86_64.whl", hash = "sha256:fb8b8ee99b3fffe4fd86f4c81b35a6bf7e4462cba019997af2fe679365db0c49"},
pytest-cov = [] {file = "coverage-6.2-cp37-cp37m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:619346d57c7126ae49ac95b11b0dc8e36c1dd49d148477461bb66c8cf13bb521"},
scipy = [] {file = "coverage-6.2-cp37-cp37m-manylinux_2_5_i686.manylinux1_i686.manylinux_2_12_i686.manylinux2010_i686.whl", hash = "sha256:0a7726f74ff63f41e95ed3a89fef002916c828bb5fcae83b505b49d81a066884"},
six = [] {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"},
stack-data = [] {file = "coverage-6.2-cp37-cp37m-musllinux_1_1_aarch64.whl", hash = "sha256:17e6c11038d4ed6e8af1407d9e89a2904d573be29d51515f14262d7f10ef0a64"},
syrupy = [] {file = "coverage-6.2-cp37-cp37m-musllinux_1_1_i686.whl", hash = "sha256:c254b03032d5a06de049ce8bca8338a5185f07fb76600afff3c161e053d88617"},
tomli = [] {file = "coverage-6.2-cp37-cp37m-musllinux_1_1_x86_64.whl", hash = "sha256:dca38a21e4423f3edb821292e97cec7ad38086f84313462098568baedf4331f8"},
traitlets = [] {file = "coverage-6.2-cp37-cp37m-win32.whl", hash = "sha256:600617008aa82032ddeace2535626d1bc212dfff32b43989539deda63b3f36e4"},
typing-extensions = [] {file = "coverage-6.2-cp37-cp37m-win_amd64.whl", hash = "sha256:bf154ba7ee2fd613eb541c2bc03d3d9ac667080a737449d1a3fb342740eb1a74"},
wcwidth = [] {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"},
]
scipy = [
{file = "scipy-1.5.4-cp36-cp36m-macosx_10_9_x86_64.whl", hash = "sha256:4f12d13ffbc16e988fa40809cbbd7a8b45bc05ff6ea0ba8e3e41f6f4db3a9e47"},
{file = "scipy-1.5.4-cp36-cp36m-manylinux1_i686.whl", hash = "sha256:a254b98dbcc744c723a838c03b74a8a34c0558c9ac5c86d5561703362231107d"},
{file = "scipy-1.5.4-cp36-cp36m-manylinux1_x86_64.whl", hash = "sha256:368c0f69f93186309e1b4beb8e26d51dd6f5010b79264c0f1e9ca00cd92ea8c9"},
{file = "scipy-1.5.4-cp36-cp36m-manylinux2014_aarch64.whl", hash = "sha256:4598cf03136067000855d6b44d7a1f4f46994164bcd450fb2c3d481afc25dd06"},
{file = "scipy-1.5.4-cp36-cp36m-win32.whl", hash = "sha256:e98d49a5717369d8241d6cf33ecb0ca72deee392414118198a8e5b4c35c56340"},
{file = "scipy-1.5.4-cp36-cp36m-win_amd64.whl", hash = "sha256:65923bc3809524e46fb7eb4d6346552cbb6a1ffc41be748535aa502a2e3d3389"},
{file = "scipy-1.5.4-cp37-cp37m-macosx_10_9_x86_64.whl", hash = "sha256:9ad4fcddcbf5dc67619379782e6aeef41218a79e17979aaed01ed099876c0e62"},
{file = "scipy-1.5.4-cp37-cp37m-manylinux1_i686.whl", hash = "sha256:f87b39f4d69cf7d7529d7b1098cb712033b17ea7714aed831b95628f483fd012"},
{file = "scipy-1.5.4-cp37-cp37m-manylinux1_x86_64.whl", hash = "sha256:25b241034215247481f53355e05f9e25462682b13bd9191359075682adcd9554"},
{file = "scipy-1.5.4-cp37-cp37m-manylinux2014_aarch64.whl", hash = "sha256:fa789583fc94a7689b45834453fec095245c7e69c58561dc159b5d5277057e4c"},
{file = "scipy-1.5.4-cp37-cp37m-win32.whl", hash = "sha256:d6d25c41a009e3c6b7e757338948d0076ee1dd1770d1c09ec131f11946883c54"},
{file = "scipy-1.5.4-cp37-cp37m-win_amd64.whl", hash = "sha256:2c872de0c69ed20fb1a9b9cf6f77298b04a26f0b8720a5457be08be254366c6e"},
{file = "scipy-1.5.4-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:e360cb2299028d0b0d0f65a5c5e51fc16a335f1603aa2357c25766c8dab56938"},
{file = "scipy-1.5.4-cp38-cp38-manylinux1_i686.whl", hash = "sha256:3397c129b479846d7eaa18f999369a24322d008fac0782e7828fa567358c36ce"},
{file = "scipy-1.5.4-cp38-cp38-manylinux1_x86_64.whl", hash = "sha256:168c45c0c32e23f613db7c9e4e780bc61982d71dcd406ead746c7c7c2f2004ce"},
{file = "scipy-1.5.4-cp38-cp38-manylinux2014_aarch64.whl", hash = "sha256:213bc59191da2f479984ad4ec39406bf949a99aba70e9237b916ce7547b6ef42"},
{file = "scipy-1.5.4-cp38-cp38-win32.whl", hash = "sha256:634568a3018bc16a83cda28d4f7aed0d803dd5618facb36e977e53b2df868443"},
{file = "scipy-1.5.4-cp38-cp38-win_amd64.whl", hash = "sha256:b03c4338d6d3d299e8ca494194c0ae4f611548da59e3c038813f1a43976cb437"},
{file = "scipy-1.5.4-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:3d5db5d815370c28d938cf9b0809dade4acf7aba57eaf7ef733bfedc9b2474c4"},
{file = "scipy-1.5.4-cp39-cp39-manylinux1_i686.whl", hash = "sha256:6b0ceb23560f46dd236a8ad4378fc40bad1783e997604ba845e131d6c680963e"},
{file = "scipy-1.5.4-cp39-cp39-manylinux1_x86_64.whl", hash = "sha256:ed572470af2438b526ea574ff8f05e7f39b44ac37f712105e57fc4d53a6fb660"},
{file = "scipy-1.5.4-cp39-cp39-manylinux2014_aarch64.whl", hash = "sha256:8c8d6ca19c8497344b810b0b0344f8375af5f6bb9c98bd42e33f747417ab3f57"},
{file = "scipy-1.5.4-cp39-cp39-win32.whl", hash = "sha256:d84cadd7d7998433334c99fa55bcba0d8b4aeff0edb123b2a1dfcface538e474"},
{file = "scipy-1.5.4-cp39-cp39-win_amd64.whl", hash = "sha256:cc1f78ebc982cd0602c9a7615d878396bec94908db67d4ecddca864d049112f2"},
{file = "scipy-1.5.4.tar.gz", hash = "sha256:4a453d5e5689de62e5d38edf40af3f17560bfd63c9c5bd228c18c1f99afa155b"},
]
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"},
]

View File

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

View File

@ -1,4 +1,4 @@
const pattern = /(\[tool\.poetry\]\nname = "pdme"\nversion = ")(?<vers>\d+\.\d+\.\d+)(")/mg; const pattern = /(\[tool\.poetry\]\nname = "pdme"\nversion = ")(?<vers>\d+\.\d+\.\d)(")/mg;
module.exports.readVersion = function (contents) { module.exports.readVersion = function (contents) {
const result = pattern.exec(contents); const result = pattern.exec(contents);

View File

@ -1,81 +0,0 @@
# 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,
])
# ---

View File

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

View File

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

View File

@ -2,173 +2,42 @@ import numpy
import pdme.measurement import pdme.measurement
def test_static_dipole_electric_potential(): def test_static_dipole():
d1 = pdme.measurement.OscillatingDipole( d1 = pdme.measurement.OscillatingDipole((1, 2, 3), (4, 5, 6), 7)
numpy.array((1, 2, 3)), numpy.array([4, 5, 6]), 7 d2 = pdme.measurement.OscillatingDipole((2, 5, 3), (4, -5, -6), 2)
)
d2 = pdme.measurement.OscillatingDipole(
numpy.array([2, 5, 3]), numpy.array([4, -5, -6]), 2
)
dipoles = pdme.measurement.OscillatingDipoleArrangement([d1, d2]) dipoles = pdme.measurement.OscillatingDipoleArrangement([d1, d2])
dot_position1 = numpy.array([-1, -1, -1]) dot_position1 = (-1, -1, -1)
dot_frequency1 = 11 dot_frequency1 = 11
expected_v1 = 0.00001421963287022476
expected_v2 = 0.00001107180225755457
expected_v1 = 0.00001221710876727626 numpy.testing.assert_allclose(d1.s_at_position(dot_position1, dot_frequency1), expected_v1, err_msg="Voltage at dot isn't as expected.")
expected_v2 = 7.257229625870065e-6
numpy.testing.assert_allclose( dot_measurements = dipoles.get_dot_measurements([(dot_position1, dot_frequency1)])
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." assert len(dot_measurements) == 1, "Should have only had one dot measurement."
measurement = dot_measurements[0] measurement = dot_measurements[0]
numpy.testing.assert_allclose( numpy.testing.assert_allclose(measurement.r, dot_position1, err_msg="Dot position should have been passed over")
measurement.r, numpy.testing.assert_allclose(measurement.f, dot_frequency1, err_msg="Dot frequency should have been passed over")
dot_position1, numpy.testing.assert_allclose(measurement.v, expected_v1 + expected_v2, err_msg="Voltage at dot isn't as expected.")
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(): def test_range_dipole_measurements():
d1 = pdme.measurement.OscillatingDipole( d1 = pdme.measurement.OscillatingDipole((1, 2, 3), (4, 5, 6), 7)
numpy.array([1, 2, 3]), numpy.array([4, 5, 6]), 7 d2 = pdme.measurement.OscillatingDipole((2, 5, 3), (4, -5, -6), 2)
)
d2 = pdme.measurement.OscillatingDipole(
numpy.array([2, 5, 3]), numpy.array([4, -5, -6]), 2
)
dipoles = pdme.measurement.OscillatingDipoleArrangement([d1, d2]) dipoles = pdme.measurement.OscillatingDipoleArrangement([d1, d2])
dot_position1 = numpy.array([-1, -1, -1]) dot_position1 = (-1, -1, -1)
dot_frequency1 = 11 dot_frequency1 = 11
expected_v1 = 0.00001221710876727626 expected_v1 = 0.00001421963287022476
expected_v2 = 7.257229625870065e-6 expected_v2 = 0.00001107180225755457
numpy.testing.assert_allclose( numpy.testing.assert_allclose(d1.s_at_position(dot_position1, dot_frequency1), expected_v1, err_msg="Voltage at dot isn't as expected.")
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( range_dot_measurements = dipoles.get_percent_range_dot_measurements([(dot_position1, dot_frequency1)], 0.5, 1.5)
[(dot_position1, dot_frequency1)], 0.5, 1.5
)
assert len(range_dot_measurements) == 1, "Should have only had one dot measurement." assert len(range_dot_measurements) == 1, "Should have only had one dot measurement."
range_measurement = range_dot_measurements[0] range_measurement = range_dot_measurements[0]
numpy.testing.assert_allclose( numpy.testing.assert_allclose(range_measurement.r, dot_position1, err_msg="Dot position should have been passed over")
range_measurement.r, numpy.testing.assert_allclose(range_measurement.f, dot_frequency1, err_msg="Dot frequency should have been passed over")
dot_position1, numpy.testing.assert_allclose(range_measurement.v_low, (expected_v1 + expected_v2) / 2, err_msg="Lower voltage at dot isn't as expected.")
err_msg="Dot position should have been passed over", numpy.testing.assert_allclose(range_measurement.v_high, (expected_v1 + expected_v2) * 3 / 2, err_msg="Lower oltage at dot isn't as expected.")
)
numpy.testing.assert_allclose(
range_measurement.f,
dot_frequency1,
err_msg="Dot frequency should have been passed over",
)
numpy.testing.assert_allclose(
range_measurement.v_low,
(expected_v1 + expected_v2) / 2,
err_msg="Lower voltage at dot isn't as expected.",
)
numpy.testing.assert_allclose(
range_measurement.v_high,
(expected_v1 + expected_v2) * 3 / 2,
err_msg="Lower oltage at dot isn't as expected.",
)

View File

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

View File

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

View File

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

View File

@ -1,55 +0,0 @@
# 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]]),
),
])
# ---

View File

@ -1,55 +0,0 @@
# 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]]),
),
])
# ---

View File

@ -1,55 +0,0 @@
# 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]]),
),
])
# ---

View File

@ -1,31 +0,0 @@
# 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,
]),
]),
])
# ---

View File

@ -0,0 +1,27 @@
from pdme.model import FixedDipoleModel
from pdme.measurement import OscillatingDipole, OscillatingDipoleArrangement
import logging
import numpy
import itertools
import pytest
@pytest.mark.skip(reason="No idea why this is failing, but it shouldn't!")
def test_fixed_dipole_model_solve_basic():
# Initialise our dipole arrangement and create dot measurements along a square.
fixed_dipole_moment = numpy.array((2, 1, 2))
dipoles = OscillatingDipoleArrangement([OscillatingDipole(fixed_dipole_moment, (1, 2, 4), 6)])
dot_inputs = list(itertools.chain.from_iterable(
(([1, 2, 1], f), ([1, 1, -2], f), ([1.5, 2, 0.1], f), ([1.5, 1, -0.2], f), ([2, 1, 0], f), ([2, 2, 0], f), ([0, 2, -.1], f), ([0, 1, 0.04], f), ([2, 0, 2], f), ([1, 0, 0], f)) for f in numpy.arange(1, 10, 1)
))
dots = dipoles.get_dot_measurements(dot_inputs)
model = FixedDipoleModel(-3, 3, -3, 3, 3, 5, fixed_dipole_moment, 1)
# from the dipole, these are the unspecified variables in ((.5, 0, 2), (1, 2, 4), 8)
expected_solution = [1, 2, 4, 6]
result = model.solve(dots)
logging.info(result)
assert result.success
numpy.testing.assert_allclose(result.normalised_x, expected_solution, err_msg="Even well specified problem solution was wrong.", rtol=1e-6, atol=1e-11)

View File

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

View File

@ -0,0 +1,39 @@
from pdme.model import FixedZPlaneModel
from pdme.measurement import OscillatingDipole, OscillatingDipoleArrangement
import logging
import numpy
import itertools
import pytest
def test_fixed_z_plane_model_solve_error_initial():
model = FixedZPlaneModel(4, -10, 10, -10, 10, 1)
with pytest.raises(ValueError):
model.solve([], initial_pt=[1, 2])
def test_fixed_z_plane_model_solve_basic():
# Initialise our dipole arrangement and create dot measurements along a square.
dipoles = OscillatingDipoleArrangement([OscillatingDipole((0, 0, 2), (1, 2, 4), 1)])
dot_inputs = list(itertools.chain.from_iterable(
(([1, 2, 0], f), ([1, 1, 0], f), ([2, 1, 0], f), ([2, 2, 0], f)) for f in numpy.arange(1, 10, 2)
))
dots = dipoles.get_dot_measurements(dot_inputs)
model = FixedZPlaneModel(4, -10, 10, -10, 10, 1)
# from the dipole, these are the unspecified variables in ((0, 0, pz), (sx, sy, 4), w)
expected_solution = [2, 1, 2, 1]
result = model.solve(dots)
logging.info(result)
assert result.success
numpy.testing.assert_allclose(result.x, expected_solution, err_msg="Even well specified problem solution was wrong.", rtol=1e-6, atol=1e-11)
# Do it again with an initial point
result = model.solve(dots, initial_pt=[2, 2, 2, 2])
logging.info(result)
assert result.success
numpy.testing.assert_allclose(result.x, expected_solution, err_msg="Even well specified problem solution was wrong.", rtol=1e-6, atol=1e-11)

View File

@ -0,0 +1,15 @@
from pdme.model.fixed_z_plane_model import FixedZPlaneModel, FixedZPlaneDiscretisation
import numpy
def test_fixed_z_plane_model_discretization():
model = FixedZPlaneModel(4, -10, 10, -10, 10, 1)
discretisation = FixedZPlaneDiscretisation(model, 1, 2, 5, 15)
# x: (-10, 0) and (0, 10)
# y: (-10, -6, -2, 2, 6, 10)
assert discretisation.cell_count == 10
assert discretisation.pz_step == 30
assert discretisation.x_step == 10
assert discretisation.y_step == 4
numpy.testing.assert_allclose(discretisation.bounds((0, 0, 0)), ((-15, -10, -10, -numpy.inf), (15, 0, -6, numpy.inf)))
numpy.testing.assert_allclose(list(discretisation.all_indices()), list(numpy.ndindex((1, 2, 5))))

View File

@ -0,0 +1,32 @@
from pdme.model.fixed_z_plane_model import FixedZPlaneModel
from pdme.measurement import DotMeasurement
import logging
import numpy
def test_fixed_z_plane_model_repr():
model = FixedZPlaneModel(1, 2, 3, 4, 5, 6)
assert repr(model) == "FixedZPlaneModel(1, 2, 3, 4, 5, 6)"
def test_fixed_z_plane_model_cost_and_jac_single():
model = FixedZPlaneModel(4, -10, 10, -10, 10, 1)
measured_v = 0.000191292 # from dipole with p=(0, 0, 2) at (1, 2, 4) with w = 1
dot = DotMeasurement(measured_v, (1, 2, 0), 5)
pt = [2, 2, 2, 2]
cost_function = model.costs([dot])
expected_cost = [0.0000946746]
actual_cost = cost_function(pt)
numpy.testing.assert_allclose(actual_cost, expected_cost, err_msg="Cost wasn't as expected.", rtol=1e-6, atol=1e-11)
jac_function = model.jac([dot])
expected_jac = [[0.0002859666151836802, -0.0001009293935942401, 0, 0.0001035396365320221]]
actual_jac = jac_function(pt)
logging.warning(actual_jac)
numpy.testing.assert_allclose(actual_jac, expected_jac, err_msg="Jac wasn't as expected.", rtol=1e-6, atol=1e-11)

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

@ -1,147 +0,0 @@
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.",
)

View File

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

View File

@ -0,0 +1,18 @@
from pdme.model.unrestricted_model import UnrestrictedModel, UnrestrictedDiscretisation
import numpy
def test_unrestricted_model_discretization():
model = UnrestrictedModel(-10, 10, -10, 10, -10, 10, 15, 1)
discretisation = UnrestrictedDiscretisation(model, 1, 1, 2, 2, 5, 1)
# x: (-10, 0) and (0, 10)
# y: (-10, -6, -2, 2, 6, 10)
assert discretisation.cell_count == 10
assert discretisation.px_step == 30
assert discretisation.py_step == 30
assert discretisation.pz_step == 15
assert discretisation.x_step == 10
assert discretisation.y_step == 4
assert discretisation.z_step == 20
numpy.testing.assert_allclose(discretisation.bounds((0, 0, 0, 0, 0, 0)), ((-15, -15, -15, -10, -10, -10, -numpy.inf), (15, 15, 0, 0, -6, 10, numpy.inf)))
numpy.testing.assert_allclose(list(discretisation.all_indices()), list(numpy.ndindex((1, 1, 2, 2, 5, 1))))

View File

@ -0,0 +1,38 @@
from pdme.model import UnrestrictedModel
from pdme.measurement import DotMeasurement
import logging
import numpy
def test_unrestricted_plane_model_repr():
model = UnrestrictedModel(1, 2, 3, 4, 5, 6, 7, 6)
assert repr(model) == "UnrestrictedModel(1, 2, 3, 4, 5, 6, 7, 6)"
def test_unrestricted_model_cost_and_jac_single():
model = UnrestrictedModel(1, -1, 1, -1, 1, -1, 1, 1)
measured_v = 0.000191292 # from dipole with p=(0, 0, 2) at (1, 2, 4) with w = 1
dot = DotMeasurement(measured_v, (1, 2, 0), 5)
pt = numpy.array([0, 0, 2, 2, 2, 4, 2])
cost_function = model.costs([dot])
expected_cost = [0.0000946746]
actual_cost = cost_function(pt)
numpy.testing.assert_allclose(actual_cost, expected_cost, err_msg="Cost wasn't as expected.", rtol=1e-6, atol=1e-11)
jac_function = model.jac([dot])
expected_jac = [
[
0.00007149165379592005, 0, 0.0002859666151836802,
-0.0001009293935942401, 0, -0.0002607342667851202,
0.0001035396365320221
]
]
actual_jac = jac_function(pt)
logging.warning(actual_jac)
numpy.testing.assert_allclose(actual_jac, expected_jac, err_msg="Jac wasn't as expected.", rtol=1e-6, atol=1e-11)

View File

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

View File

@ -1,94 +0,0 @@
# 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,
]),
]),
])
# ---

View File

@ -1,22 +0,0 @@
# 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)
# ---

View File

@ -1,31 +0,0 @@
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)

View File

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

View File

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

View File

@ -1,20 +0,0 @@
# 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,
]),
])
# ---

View File

@ -1,25 +0,0 @@
# 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,
]),
])
# ---

View File

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

View File

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

View File

@ -1,30 +1,6 @@
import numpy import numpy
import pdme.util.fast_v_calc 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(): def test_fast_v_calc():
d1 = [1, 2, 3, 4, 5, 6, 7] d1 = [1, 2, 3, 4, 5, 6, 7]
@ -33,184 +9,22 @@ def test_fast_v_calc():
dipoles = numpy.array([d1, d2]) dipoles = numpy.array([d1, d2])
dot_inputs = numpy.array([[-1, -1, -1, 11], [2, 3, 1, 5.5]]) dot_inputs = numpy.array([[-1, -1, -1, 11], [2, 3, 1, 5.5]])
# expected_ij is for dot i, dipole j
expected_11 = s_potential_from_arrays(dipoles[0], dot_inputs[0]) expected_11 = 0.00001421963287022476
expected_12 = s_potential_from_arrays(dipoles[1], dot_inputs[0]) expected_12 = 0.00001107180225755457
expected_21 = s_potential_from_arrays(dipoles[0], dot_inputs[1]) expected_21 = 0.000345021108583681380388722
expected_22 = s_potential_from_arrays(dipoles[1], dot_inputs[1]) expected_22 = 0.0000377061050587914705139781
expected = numpy.array([[expected_11, expected_21], [expected_12, expected_22]]) expected = numpy.array([[expected_11, expected_21], [expected_12, expected_22]])
numpy.testing.assert_allclose( 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.")
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(): def test_between():
low = numpy.array([1, 2, 3]) low = numpy.array([1, 2, 3])
high = numpy.array([6, 7, 8]) high = numpy.array([6, 7, 8])
# FALSE FALSE TRUE # FALSE FALSE TRUE
a = [[0, 1, 2], [0, 9, 5], [4, 5, 6]] a = [[0, 1, 2], [0, 9, 5], [4, 5, 6]]
actual = pdme.util.fast_v_calc.between(a, low, high) actual = pdme.util.fast_v_calc.between(a, low, high)