physics/pdme
2
0
Fork 0
Code Issues Pull Requests 8 Projects Releases Wiki Activity

Merge pull request 'mcmc' (#30) from mcmc into master
All checks were successful
gitea-physics/pdme/pipeline/head This commit looks good
Details

Browse Source 
Reviewed-on: #30
This commit is contained in:
Deepak Mallubhotla 2023-07-23 23:30:05 +00:00
commit dbf924490e
23 changed files with 976 additions and 56 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
do.sh
View File

@ -39,6 +39,15 @@ test() {
fi
}
updatesnap() {
echo "I am ${FUNCNAME[0]}ing"
if [[ "${DO_NIX_CUSTOM:=0}" -eq 1 ]]; then
pytest --snapshot-update
else
poetry run pytest --snapshot-update
fi
}
htmlcov() {
if [[ "${DO_NIX_CUSTOM:=0}" -eq 1 ]]; then
pytest --cov-report=html

53
pdme/model/log_spaced_random_choice_fixed_orientation_model.py
View File

@ -5,6 +5,11 @@ from pdme.measurement import (
OscillatingDipole,
OscillatingDipoleArrangement,
)
import logging
from typing import Optional
import pdme.subspace_simulation
_logger = logging.getLogger(__name__)
class LogSpacedRandomCountMultipleDipoleFixedMagnitudeFixedOrientationModel(
@ -138,3 +143,51 @@ class LogSpacedRandomCountMultipleDipoleFixedMagnitudeFixedOrientationModel(
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

72
pdme/model/log_spaced_random_choice_model.py
View File

@ -5,6 +5,8 @@ from pdme.measurement import (
OscillatingDipole,
OscillatingDipoleArrangement,
)
import pdme.subspace_simulation
from typing import Optional
class LogSpacedRandomCountMultipleDipoleFixedMagnitudeModel(DipoleModel):
@ -125,3 +127,73 @@ class LogSpacedRandomCountMultipleDipoleFixedMagnitudeModel(DipoleModel):
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

58
pdme/model/log_spaced_random_choice_xy_model.py
View File

@ -5,6 +5,8 @@ from pdme.measurement import (
OscillatingDipole,
OscillatingDipoleArrangement,
)
import pdme.subspace_simulation
from typing import Optional
class LogSpacedRandomCountMultipleDipoleFixedMagnitudeXYModel(DipoleModel):
@ -123,3 +125,59 @@ class LogSpacedRandomCountMultipleDipoleFixedMagnitudeXYModel(DipoleModel):
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

67
pdme/model/model.py
View File

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

55
pdme/subspace_simulation/__init__.py Normal file
View File

@ -0,0 +1,55 @@
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.
Utility function.
"""
return numpy.array(sorted(configuration, key=lambda l: l[6]))
__all__ = [
"DipoleStandardDeviation",
"MCMCStandardDeviation",
"sort_array_of_dipoles_by_frequency",
"proportional_cost",
"proportional_costs_vs_actual_measurement",
]

20
pdme/subspace_simulation/mcmc_costs.py Normal file
View File

@ -0,0 +1,20 @@
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, numpy.array([dipoles_to_test])
)
return proportional_cost(actual_measurement_array, vals)

90
poetry.lock generated
View File

@ -20,28 +20,6 @@ six = "*"
[package.extras]
test = ["astroid", "pytest"]
[[package]]
name = "atomicwrites"
version = "1.4.0"
description = "Atomic file writes."
category = "dev"
optional = false
python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*"
[[package]]
name = "attrs"
version = "21.4.0"
description = "Classes Without Boilerplate"
category = "dev"
optional = false
python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.*"
[package.extras]
dev = ["coverage[toml] (>=5.0.2)", "hypothesis", "pympler", "pytest (>=4.3.0)", "six", "mypy", "pytest-mypy-plugins", "zope.interface", "furo", "sphinx", "sphinx-notfound-page", "pre-commit", "cloudpickle"]
docs = ["furo", "sphinx", "zope.interface", "sphinx-notfound-page"]
tests = ["coverage[toml] (>=5.0.2)", "hypothesis", "pympler", "pytest (>=4.3.0)", "six", "mypy", "pytest-mypy-plugins", "zope.interface", "cloudpickle"]
tests_no_zope = ["coverage[toml] (>=5.0.2)", "hypothesis", "pympler", "pytest (>=4.3.0)", "six", "mypy", "pytest-mypy-plugins", "cloudpickle"]
[[package]]
name = "backcall"
version = "0.2.0"
@ -91,9 +69,17 @@ category = "dev"
optional = false
python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.*"
[[package]]
name = "colored"
version = "1.4.4"
description = "Simple library for color and formatting to terminal"
category = "dev"
optional = false
python-versions = "*"
[[package]]
name = "coverage"
version = "7.0.3"
version = "7.2.7"
description = "Code coverage measurement for Python"
category = "dev"
optional = false
@ -113,6 +99,17 @@ 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"
@ -349,14 +346,6 @@ python-versions = "*"
[package.extras]
tests = ["pytest"]
[[package]]
name = "py"
version = "1.11.0"
description = "library with cross-python path, ini-parsing, io, code, log facilities"
category = "dev"
optional = false
python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*, !=3.4.*"
[[package]]
name = "pycodestyle"
version = "2.8.0"
@ -394,32 +383,30 @@ diagrams = ["jinja2", "railroad-diagrams"]
[[package]]
name = "pytest"
version = "6.2.5"
version = "7.4.0"
description = "pytest: simple powerful testing with Python"
category = "dev"
optional = false
python-versions = ">=3.6"
python-versions = ">=3.7"
[package.dependencies]
atomicwrites = {version = ">=1.0", markers = "sys_platform == \"win32\""}
attrs = ">=19.2.0"
colorama = {version = "*", markers = "sys_platform == \"win32\""}
exceptiongroup = {version = ">=1.0.0rc8", markers = "python_version < \"3.11\""}
iniconfig = "*"
packaging = "*"
pluggy = ">=0.12,<2.0"
py = ">=1.8.2"
toml = "*"
tomli = {version = ">=1.0.0", markers = "python_version < \"3.11\""}
[package.extras]
testing = ["argcomplete", "hypothesis (>=3.56)", "mock", "nose", "requests", "xmlschema"]
testing = ["argcomplete", "attrs (>=19.2.0)", "hypothesis (>=3.56)", "mock", "nose", "pygments (>=2.7.2)", "requests", "setuptools", "xmlschema"]
[[package]]
name = "pytest-cov"
version = "3.0.0"
version = "4.1.0"
description = "Pytest plugin for measuring coverage."
category = "dev"
optional = false
python-versions = ">=3.6"
python-versions = ">=3.7"
[package.dependencies]
coverage = {version = ">=5.2.1", extras = ["toml"]}
@ -469,12 +456,16 @@ pure-eval = "*"
tests = ["pytest", "typeguard", "pygments", "littleutils", "cython"]
[[package]]
name = "toml"
version = "0.10.2"
description = "Python Library for Tom's Obvious, Minimal Language"
name = "syrupy"
version = "4.0.8"
description = "Pytest Snapshot Test Utility"
category = "dev"
optional = false
python-versions = ">=2.6, !=3.0.*, !=3.1.*, !=3.2.*"
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"
@ -513,20 +504,20 @@ python-versions = "*"
[metadata]
lock-version = "1.1"
python-versions = "^3.8,<3.10"
content-hash = "f5614947eb2f77e9cd8bd9dc026a924f4cb7394c3fef9cb175ab5afbe201ca73"
python-versions = ">=3.8.1,<3.10"
content-hash = "b5275c33449e8f85acbf9c0f6dfe1ec4e7296adfa16360d782b33534e1223638"
[metadata.files]
appnope = []
asttokens = []
atomicwrites = []
attrs = []
backcall = []
black = []
click = []
colorama = []
colored = []
coverage = []
decorator = []
exceptiongroup = []
executing = []
flake8 = []
iniconfig = []
@ -547,7 +538,6 @@ pluggy = []
prompt-toolkit = []
ptyprocess = []
pure-eval = []
py = []
pycodestyle = []
pyflakes = []
pygments = []
@ -557,7 +547,7 @@ pytest-cov = []
scipy = []
six = []
stack-data = []
toml = []
syrupy = []
tomli = []
traitlets = []
typing-extensions = []

5
pyproject.toml
View File

@ -7,17 +7,18 @@ license = "GPL-3.0-only"
readme = "README.md"
[tool.poetry.dependencies]
python = "^3.8,<3.10"
python = ">=3.8.1,<3.10"
numpy = "^1.22.3"
scipy = "~1.10"
[tool.poetry.dev-dependencies]
pytest = ">=6"
flake8 = "^4.0.0"
pytest-cov = "^3.0.0"
pytest-cov = "^4.1.0"
mypy = "^0.961"
ipython = "^8.2.0"
black = "^22.3.0"
syrupy = "^4.0.8"
[build-system]
requires = ["poetry-core>=1.0.0"]

55
tests/model/__snapshots__/test_log_spaced_fixed_orientation_fixedxy_orientation_mcmc.ambr Normal file
View File

@ -0,0 +1,55 @@
# serializer version: 1
# name: test_log_spaced_fixedxy_orientation_mcmc_basic
list([
tuple(
array([3984.46179656]),
array([[ 9.55610128, 2.94634152, 0. , 9.21529051, -2.46576127,
2.42481096, 9.19034554]]),
),
tuple(
array([8583.99087872]),
array([[ 9.99991539, 0.04113671, 0. , 8.71258954, -2.26599865,
2.60452102, 6.37042214]]),
),
tuple(
array([6215.6376616]),
array([[ 9.81950685, -1.89137124, 0. , 8.90637055, -2.48043039,
2.28444435, 8.84239221]]),
),
tuple(
array([424.73328466]),
array([[ 1.00028483, 9.94984574, 0. , 8.53064898, -2.59230757,
2.33774773, 8.6714416 ]]),
),
tuple(
array([300.92203808]),
array([[ 1.4003442 , 9.90146636, 0. , 8.05557992, -2.6753126 ,
2.65915755, 13.02021385]]),
),
tuple(
array([2400.01072771]),
array([[ 9.97761813, 0.66868263, 0. , 8.69171028, -2.73145011,
2.90140456, 19.94999593]]),
),
tuple(
array([5001.46205113]),
array([[ 9.93976109, -1.09596962, 0. , 8.95245025, -2.59409162,
2.90140456, 9.75535945]]),
),
tuple(
array([195.21980745]),
array([[ 0.20690762, 9.99785923, 0. , 9.59636585, -2.83240984,
2.90140456, 16.14771567]]),
),
tuple(
array([2698.2588445]),
array([[-9.68130127, -2.50447712, 0. , 8.94823619, -2.92889659,
2.77065328, 13.63173263]]),
),
tuple(
array([1193.69854739]),
array([[-6.16597091, -7.87278875, 0. , 9.62210721, -2.75993924,
2.77065328, 5.64553534]]),
),
])
# ---

55
tests/model/__snapshots__/test_log_spaced_fixed_orientation_free_orientation_mcmc.ambr Normal file
View File

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

55
tests/model/__snapshots__/test_log_spaced_fixed_orientation_mcmc.ambr Normal file
View File

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

31
tests/model/__snapshots__/test_log_spaced_random_count_multiple_dipole_fixed_orientation_fixed_magnitude_model.ambr Normal file
View File

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

90
tests/model/test_log_spaced_fixed_orientation_fixedxy_orientation_mcmc.py Normal file
View File

@ -0,0 +1,90 @@
from pdme.model import (
LogSpacedRandomCountMultipleDipoleFixedMagnitudeXYModel,
)
import pdme.inputs
import pdme.measurement.input_types
import pdme.subspace_simulation
import numpy
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_dot_measurements(dot_inputs)
actual_measurements_array = numpy.array([m.v for m in actual_measurements])
def cost_to_use(sample_dipoles: numpy.ndarray) -> numpy.ndarray:
return pdme.subspace_simulation.proportional_costs_vs_actual_measurement(
dot_input_array, actual_measurements_array, sample_dipoles
)
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)[0]
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, cost_function, 10, stdevs, rng_arg=numpy.random.default_rng(1515)
)
assert chain == snapshot

90
tests/model/test_log_spaced_fixed_orientation_free_orientation_mcmc.py Normal file
View File

@ -0,0 +1,90 @@
from pdme.model import (
LogSpacedRandomCountMultipleDipoleFixedMagnitudeModel,
)
import pdme.inputs
import pdme.measurement.input_types
import pdme.subspace_simulation
import numpy
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_dot_measurements(dot_inputs)
actual_measurements_array = numpy.array([m.v for m in actual_measurements])
def cost_to_use(sample_dipoles: numpy.ndarray) -> numpy.ndarray:
return pdme.subspace_simulation.proportional_costs_vs_actual_measurement(
dot_input_array, actual_measurements_array, sample_dipoles
)
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)[0]
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, cost_function, 10, stdevs, rng_arg=numpy.random.default_rng(1515)
)
assert chain == snapshot

98
tests/model/test_log_spaced_fixed_orientation_mcmc.py Normal file
View File

@ -0,0 +1,98 @@
from pdme.model import (
LogSpacedRandomCountMultipleDipoleFixedMagnitudeFixedOrientationModel,
)
import pdme.inputs
import pdme.measurement.input_types
import pdme.subspace_simulation
import numpy
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_dot_measurements(dot_inputs)
actual_measurements_array = numpy.array([m.v for m in actual_measurements])
def cost_to_use(sample_dipoles: numpy.ndarray) -> numpy.ndarray:
return pdme.subspace_simulation.proportional_costs_vs_actual_measurement(
dot_input_array, actual_measurements_array, sample_dipoles
)
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)[0]
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, cost_function, 10, stdevs, rng_arg=numpy.random.default_rng(1515)
)
assert chain == snapshot

8
tests/model/test_log_spaced_random_count_multiple_dipole_fixed_orientation_fixed_magnitude_model.py
View File

@ -115,7 +115,7 @@ def test_random_count_multiple_dipole_fixed_mag_model_get_dipoles_invariant():
)
def test_random_count_multiple_dipole_fixed_or_fixed_mag_model_get_n_dipoles():
def test_random_count_multiple_dipole_fixed_or_fixed_mag_model_get_n_dipoles(snapshot):
# TODO: this test is a bit garbage just calls things without testing.
x_min = -10
x_max = 10
@ -145,7 +145,9 @@ def test_random_count_multiple_dipole_fixed_or_fixed_mag_model_get_n_dipoles():
)
model.rng = numpy.random.default_rng(1234)
model.get_monte_carlo_dipole_inputs(1, max_frequency)
model.get_monte_carlo_dipole_inputs(
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

4
tests/subspace_simulation/__snapshots__/test_costs.ambr Normal file
View File

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

32
tests/subspace_simulation/__snapshots__/test_sort_dipoles.ambr Normal file
View File

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

22
tests/subspace_simulation/__snapshots__/test_stdevs.ambr Normal file
View File

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

10
tests/subspace_simulation/test_costs.py Normal file
View File

@ -0,0 +1,10 @@
import pdme.subspace_simulation
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

15
tests/subspace_simulation/test_sort_dipoles.py Normal file
View File

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

38
tests/subspace_simulation/test_stdevs.py Normal file
View File

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