Compare commits

...

15 Commits

Author SHA1 Message Date
885508e104
chore(release): 1.5.0
All checks were successful
gitea-physics/pdme/pipeline/head This commit looks good
gitea-physics/pdme/pipeline/tag This commit looks good
2024-05-16 23:12:43 -05:00
6193ecb9c9
feat: adds mcmc chain that returns number of repeats
All checks were successful
gitea-physics/pdme/pipeline/head This commit looks good
2024-05-16 23:12:15 -05:00
5ad442750e
chore(release): 1.4.0
All checks were successful
gitea-physics/pdme/pipeline/tag This commit looks good
gitea-physics/pdme/pipeline/head This commit looks good
2024-05-16 21:07:02 -05:00
9b1538b3c6
feat: adds relative squared diff calc utility method 2024-05-16 21:06:42 -05:00
7b277fdc85
chore(release): 1.3.0
All checks were successful
gitea-physics/pdme/pipeline/head This commit looks good
gitea-physics/pdme/pipeline/tag This commit looks good
2024-05-16 19:21:24 -05:00
e5fc1207a8
feat: adds utility function for sorting samples by frequency for subspace simulation 2024-05-16 19:20:55 -05:00
387a607e09
chore(release): 1.2.0
All checks were successful
gitea-physics/pdme/pipeline/head This commit looks good
gitea-physics/pdme/pipeline/tag This commit looks good
2024-05-02 20:56:39 -05:00
9e6d1df559
feat: adds pdme fast calc for e field xs 2024-05-02 20:56:15 -05:00
45031857f2
doc: adds more index comments 2024-05-02 20:26:21 -05:00
64181eeef2
chore: clean up pointless nb file
All checks were successful
gitea-physics/pdme/pipeline/head This commit looks good
2024-05-02 18:09:19 -05:00
d682d50554
chore(release): 1.1.0
All checks were successful
gitea-physics/pdme/pipeline/head This commit looks good
gitea-physics/pdme/pipeline/tag This commit looks good
2024-05-02 18:06:17 -05:00
e9e34162a3
feat: adds both electric potential and electric field x sources, makes some fast util tests use the slower explicit versions as double check 2024-05-02 18:05:37 -05:00
32a7812a43
just: adds snapshot just task 2024-05-02 18:01:22 -05:00
7c39475742
log: don't log warnings anymore for util calc
All checks were successful
gitea-physics/pdme/pipeline/head This commit looks good
2024-04-28 21:48:06 -05:00
1dd7569fde
doc: now set as maintained 2024-04-28 21:38:48 -05:00
28 changed files with 1064 additions and 5207 deletions

View File

@ -2,6 +2,41 @@
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)

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 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)
![Maintenance](https://img.shields.io/maintenance/yes/2023?style=flat-square)
![Maintenance](https://img.shields.io/maintenance/yes/2024?style=flat-square)
This repo has library code for evaluating dipole models.

File diff suppressed because it is too large Load Diff

View File

@ -34,6 +34,16 @@ test: fmt
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

View File

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

View File

@ -7,6 +7,7 @@ from pdme.measurement.dot_pair_measure import (
DotPairMeasurement,
DotPairRangeMeasurement,
)
import pdme.calculations
from pdme.measurement.input_types import DotInput, DotPairInput
@ -38,10 +39,27 @@ class OscillatingDipole:
self.p = numpy.array(self.p)
self.s = numpy.array(self.s)
def s_at_position(self, r: numpy.ndarray, f: float) -> float:
def _alpha_electric_potential(self, r: numpy.ndarray) -> float:
"""
Returns the electric potential of this dipole at position r.
"""
return pdme.calculations.electric_potential(self.p, self.s, r)
def _alpha_electric_field(self, r: numpy.ndarray) -> numpy.ndarray:
"""
Returns the electric field of this dipole at position r.
"""
return pdme.calculations.electric_field(self.p, self.s, r)
def _b(self, f: float) -> float:
return pdme.calculations.telegraph_beta(f, self.w)
def s_electric_potential_at_position(self, r: numpy.ndarray, f: float) -> float:
"""
Returns the noise potential at a point r, at some frequency f.
Specifically for electric potential!
Parameters
----------
r : numpy.ndarray
@ -50,17 +68,49 @@ class OscillatingDipole:
f : float
The dot frequency to sample.
"""
return (self._alpha(r)) ** 2 * self._b(f)
return (self._alpha_electric_potential(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 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 _b(self, f: float) -> float:
return (1 / numpy.pi) * (self.w / (f**2 + self.w**2))
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.
def s_for_dot_pair(self, r1: numpy.ndarray, r2: numpy.ndarray, f: float) -> float:
return self._alpha(r1) * self._alpha(r2) * self._b(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])])
@ -78,69 +128,97 @@ class OscillatingDipoleArrangement:
def __init__(self, dipoles: Sequence[OscillatingDipole]):
self.dipoles = dipoles
def get_dot_measurement(self, dot_input: DotInput) -> DotMeasurement:
def get_potential_dot_measurement(self, dot_input: DotInput) -> DotMeasurement:
r = numpy.array(dot_input[0])
f = dot_input[1]
return DotMeasurement(
sum([dipole.s_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_dot_pair_measurement(
def get_potential_dot_pair_measurement(
self, dot_pair_input: DotPairInput
) -> DotPairMeasurement:
r1 = numpy.array(dot_pair_input[0])
r2 = numpy.array(dot_pair_input[1])
f = dot_pair_input[2]
return DotPairMeasurement(
sum([dipole.s_for_dot_pair(r1, r2, f) for dipole in self.dipoles]),
sum(
[
dipole.s_electric_potential_for_dot_pair(r1, r2, f)
for dipole in self.dipoles
]
),
r1,
r2,
f,
)
def get_dot_measurements(
def get_potential_dot_measurements(
self, dot_inputs: Sequence[DotInput]
) -> List[DotMeasurement]:
"""
For a series of points, each with three coordinates and a frequency, return a list of the corresponding DotMeasurements.
"""
return [self.get_dot_measurement(dot_input) for dot_input in dot_inputs]
return [
self.get_potential_dot_measurement(dot_input) for dot_input in dot_inputs
]
def get_dot_pair_measurements(
def get_potential_dot_pair_measurements(
self, dot_pair_inputs: Sequence[DotPairInput]
) -> List[DotPairMeasurement]:
"""
For a series of pairs of points, each with three coordinates and a frequency, return a list of the corresponding DotPairMeasurements.
"""
return [
self.get_dot_pair_measurement(dot_pair_input)
self.get_potential_dot_pair_measurement(dot_pair_input)
for dot_pair_input in dot_pair_inputs
]
def get_percent_range_dot_measurement(
def get_percent_range_potential_dot_measurement(
self, dot_input: DotInput, low_percent: float, high_percent: float
) -> DotRangeMeasurement:
r = numpy.array(dot_input[0])
f = dot_input[1]
return DotRangeMeasurement(
low_percent * sum([dipole.s_at_position(r, f) for dipole in self.dipoles]),
high_percent * sum([dipole.s_at_position(r, f) for dipole in self.dipoles]),
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_dot_measurements(
def get_percent_range_potential_dot_measurements(
self, dot_inputs: Sequence[DotInput], low_percent: float, high_percent: float
) -> List[DotRangeMeasurement]:
"""
For a series of pairs of points, each with three coordinates and a frequency, and also a lower error range and upper error range, return a list of the corresponding DotPairRangeMeasurements.
"""
return [
self.get_percent_range_dot_measurement(dot_input, low_percent, high_percent)
self.get_percent_range_potential_dot_measurement(
dot_input, low_percent, high_percent
)
for dot_input in dot_inputs
]
def get_percent_range_dot_pair_measurement(
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])
@ -148,15 +226,25 @@ class OscillatingDipoleArrangement:
f = pair_input[2]
return DotPairRangeMeasurement(
low_percent
* sum([dipole.s_for_dot_pair(r1, r2, f) for dipole in self.dipoles]),
* sum(
[
dipole.s_electric_potential_for_dot_pair(r1, r2, f)
for dipole in self.dipoles
]
),
high_percent
* sum([dipole.s_for_dot_pair(r1, r2, f) for dipole in self.dipoles]),
* sum(
[
dipole.s_electric_potential_for_dot_pair(r1, r2, f)
for dipole in self.dipoles
]
),
r1,
r2,
f,
)
def get_percent_range_dot_pair_measurements(
def get_percent_range_potential_dot_pair_measurements(
self,
pair_inputs: Sequence[DotPairInput],
low_percent: float,
@ -166,7 +254,7 @@ class OscillatingDipoleArrangement:
For a series of pairs of points, each with three coordinates and a frequency, and also a lower error range and upper error range, return a list of the corresponding DotPairRangeMeasurements.
"""
return [
self.get_percent_range_dot_pair_measurement(
self.get_percent_range_potential_dot_pair_measurement(
pair_input, low_percent, high_percent
)
for pair_input in pair_inputs

View File

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

View File

@ -41,11 +41,30 @@ 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",

View File

@ -18,3 +18,13 @@ def proportional_costs_vs_actual_measurement(
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

@ -45,7 +45,7 @@ def fast_s_nonlocal(
_logger.debug(f"alphses1: {alphses1}")
_logger.debug(f"alphses2: {alphses2}")
bses = (1 / numpy.pi) * (ws[:, None] / (f1s**2 + ws[:, None] ** 2))
bses = 2 * ws[:, None] / ((numpy.pi * f1s) ** 2 + ws[:, None] ** 2)
if _logger.isEnabledFor(logging.DEBUG):
_logger.debug(f"bses: {bses}")
@ -58,44 +58,78 @@ def fast_s_nonlocal_dipoleses(
"""
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}")
bses = (1 / numpy.pi) * (ws[:, None, :] / (f1s[:, None] ** 2 + ws[:, None, :] ** 2))
# 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)
@ -152,15 +186,15 @@ def fast_s_spin_qubit_tarucha_nonlocal_dipoleses(
diffses1 = r1s[:, None] - ss[:, None, :]
diffses2 = r2s[:, None] - ss[:, None, :]
if _logger.isEnabledFor(logging.DEBUG):
_logger.warning(f"diffses1: {diffses1}")
_logger.warning(f"diffses2: {diffses2}")
_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.warning(f"norms1: {norms1}")
_logger.warning(f"norms2: {norms2}")
_logger.debug(f"norms1: {norms1}")
_logger.debug(f"norms2: {norms2}")
# _logger.info(f"norms1: {norms1}")
# _logger.info(f"norms1 shape: {norms1.shape}")
@ -215,14 +249,19 @@ def fast_s_spin_qubit_tarucha_nonlocal_dipoleses(
- ps[:, numpy.newaxis, :, 0]
) / (norms2**3)
if _logger.isEnabledFor(logging.DEBUG):
_logger.warning(f"alphses1: {alphses1}")
_logger.warning(f"alphses2: {alphses2}")
_logger.debug(f"alphses1: {alphses1}")
_logger.debug(f"alphses2: {alphses2}")
bses = (1 / numpy.pi) * (ws[:, None, :] / (f1s[:, None] ** 2 + ws[:, None, :] ** 2))
# 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.warning(f"bses: {bses}")
_logger.debug(f"bses: {bses}")
_logger.warning(f"Raw pair calc: [{alphses1 * alphses2 * 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)

View File

@ -12,27 +12,44 @@ def fast_vs_for_dipoles(
No error correction here baby.
Expects dot_inputs to be numpy array of [rx, ry, rz, f] entries, so a n by 4 where n is number of measurement points.
"""
# name indexes:
# A: dipole config index
# cart: cartesian index
# m: measurement index
# [A, cart]
ps = dipoles[:, 0:3]
ss = dipoles[:, 3:6]
# [A]
ws = dipoles[:, 6]
_logger.debug(f"ps: {ps}")
_logger.debug(f"ss: {ss}")
_logger.debug(f"ws: {ws}")
# [m, cart]
rs = dot_inputs[:, 0:3]
# [m]
fs = dot_inputs[:, 3]
# [m, cart] - [A, 1, cart]
# [A, m, cart]
diffses = rs - ss[:, None]
_logger.debug(f"diffses: {diffses}")
# [A, m]
norms = numpy.linalg.norm(diffses, axis=2) ** 3
_logger.debug(f"norms: {norms}")
# [A, m, cart] [A, cart] -> [A, m]
ases = (numpy.einsum("...ji, ...i", diffses, ps) / norms) ** 2
_logger.debug(f"ases: {ases}")
bses = (1 / numpy.pi) * (ws[:, None] / (fs**2 + ws[:, None] ** 2))
# [A, 1], denom [m] + [A, 1] -> [A, m]
# [A, m]
bses = 2 * ws[:, None] / ((numpy.pi * fs) ** 2 + ws[:, None] ** 2)
_logger.debug(f"bses: {bses}")
# returns shape [A, m]
return ases * bses
@ -64,11 +81,126 @@ def fast_vs_for_dipoleses(
ases = (numpy.einsum("abcd,acd->abc", diffses, ps) / norms) ** 2
_logger.debug(f"ases: {ases}")
bses = (1 / numpy.pi) * (ws[:, None, :] / (fs[:, None] ** 2 + ws[:, None, :] ** 2))
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:
@ -92,6 +224,9 @@ def fast_vs_for_asymmetric_dipoleses(
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
@ -105,7 +240,7 @@ def fast_vs_for_asymmetric_dipoleses(
ases = (numpy.einsum("abcd,acd->abc", diffses, ps) / norms) ** 2
bses = (1 / numpy.pi) * (ws[:, None, :] / (fs[:, None] ** 2 + ws[:, None, :] ** 2))
bses = ws[:, None, :] / ((numpy.pi * fs[:, None]) ** 2 + ws[:, None, :] ** 2)
return numpy.einsum("...j->...", ases * bses)

View File

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

View File

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

View File

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

View File

@ -2,23 +2,30 @@ import numpy
import pdme.measurement
def test_static_dipole():
d1 = pdme.measurement.OscillatingDipole((1, 2, 3), (4, 5, 6), 7)
d2 = pdme.measurement.OscillatingDipole((2, 5, 3), (4, -5, -6), 2)
def test_static_dipole_electric_potential():
d1 = pdme.measurement.OscillatingDipole(
numpy.array((1, 2, 3)), numpy.array([4, 5, 6]), 7
)
d2 = pdme.measurement.OscillatingDipole(
numpy.array([2, 5, 3]), numpy.array([4, -5, -6]), 2
)
dipoles = pdme.measurement.OscillatingDipoleArrangement([d1, d2])
dot_position1 = (-1, -1, -1)
dot_position1 = numpy.array([-1, -1, -1])
dot_frequency1 = 11
expected_v1 = 0.00001421963287022476
expected_v2 = 0.00001107180225755457
expected_v1 = 0.00001221710876727626
expected_v2 = 7.257229625870065e-6
numpy.testing.assert_allclose(
d1.s_at_position(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_dot_measurements([(dot_position1, dot_frequency1)])
dot_measurements = dipoles.get_potential_dot_measurements(
[(dot_position1, dot_frequency1)]
)
assert len(dot_measurements) == 1, "Should have only had one dot measurement."
measurement = dot_measurements[0]
numpy.testing.assert_allclose(
@ -38,24 +45,57 @@ def test_static_dipole():
)
def test_dipole_dot_pair():
d1 = pdme.measurement.OscillatingDipole((1, 2, 3), (4, 5, 6), 7)
dipoles = pdme.measurement.OscillatingDipoleArrangement([d1])
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 = (-1, -1, -1)
dot_position2 = (1, 2, 3)
dot_frequency = 8
expected_sij = 0.000083328037100902801698
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_for_dot_pair(dot_position1, dot_position2, dot_frequency),
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_dot_pair_measurements(
for m in dipoles.get_potential_dot_pair_measurements(
[(dot_position1, dot_position2, dot_frequency)]
)
],
@ -65,15 +105,17 @@ def test_dipole_dot_pair():
def test_range_pairs():
d1 = pdme.measurement.OscillatingDipole((1, 2, 3), (4, 5, 6), 7)
d1 = pdme.measurement.OscillatingDipole(
numpy.array([1, 2, 3]), numpy.array([4, 5, 6]), 7
)
dipoles = pdme.measurement.OscillatingDipoleArrangement([d1])
dot_position1 = (-1, -1, -1)
dot_position2 = (1, 2, 3)
dot_position1 = numpy.array([-1, -1, -1])
dot_position2 = numpy.array([1, 2, 3])
dot_frequency = 8
expected_sij = 0.000083328037100902801698
expected_sij = 0.00008692058236162933
actuals = dipoles.get_percent_range_dot_pair_measurements(
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"
@ -86,22 +128,26 @@ def test_range_pairs():
def test_range_dipole_measurements():
d1 = pdme.measurement.OscillatingDipole((1, 2, 3), (4, 5, 6), 7)
d2 = pdme.measurement.OscillatingDipole((2, 5, 3), (4, -5, -6), 2)
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 = (-1, -1, -1)
dot_position1 = numpy.array([-1, -1, -1])
dot_frequency1 = 11
expected_v1 = 0.00001421963287022476
expected_v2 = 0.00001107180225755457
expected_v1 = 0.00001221710876727626
expected_v2 = 7.257229625870065e-6
numpy.testing.assert_allclose(
d1.s_at_position(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.",
)
range_dot_measurements = dipoles.get_percent_range_dot_measurements(
range_dot_measurements = dipoles.get_percent_range_potential_dot_measurements(
[(dot_position1, dot_frequency1)], 0.5, 1.5
)
assert len(range_dot_measurements) == 1, "Should have only had one dot measurement."

View File

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

View File

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

View File

@ -2,52 +2,52 @@
# name: test_log_spaced_fixed_orientation_mcmc_basic
list([
tuple(
50.5683119299,
50.5632126772,
array([[ 0. , 0. , 10. , -2.3960853 , 4.23246234,
2.26169242, 39.39900844]]),
),
tuple(
50.5683119299,
50.5632126772,
array([[ 0. , 0. , 10. , -2.3960853 , 4.23246234,
2.26169242, 39.39900844]]),
),
tuple(
47.408654554,
47.4038652151,
array([[ 0. , 0. , 10. , -2.03666518, 4.14084039,
2.21309317, 47.82371559]]),
),
tuple(
47.408654554,
47.4038652151,
array([[ 0. , 0. , 10. , -2.03666518, 4.14084039,
2.21309317, 47.82371559]]),
),
tuple(
47.408654554,
47.4038652151,
array([[ 0. , 0. , 10. , -2.03666518, 4.14084039,
2.21309317, 47.82371559]]),
),
tuple(
47.408654554,
47.4038652151,
array([[ 0. , 0. , 10. , -2.03666518, 4.14084039,
2.21309317, 47.82371559]]),
),
tuple(
22.9327902847,
22.9304785991,
array([[ 0. , 0. , 10. , -1.63019717, 3.97041764,
2.53115835, 38.2051999 ]]),
),
tuple(
28.8119773322,
28.8090658953,
array([[ 0. , 0. , 10. , -1.14570315, 4.07709911,
2.48697441, 49.58615195]]),
),
tuple(
28.8119773322,
28.8090658953,
array([[ 0. , 0. , 10. , -1.14570315, 4.07709911,
2.48697441, 49.58615195]]),
),
tuple(
40.9740600543,
40.9699102596,
array([[ 0. , 0. , 10. , -0.50178755, 3.83878089,
2.93560796, 82.07827571]]),
),

View File

@ -43,7 +43,7 @@ def get_cost_function():
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 = 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:

View File

@ -43,7 +43,7 @@ def get_cost_function():
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 = 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:

View File

@ -47,7 +47,7 @@ def get_cost_function():
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 = 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:

View File

@ -30,3 +30,65 @@
]),
])
# ---
# 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,4 +1,5 @@
import pdme.subspace_simulation
import pdme.subspace_simulation.mcmc_costs
import numpy
@ -8,3 +9,23 @@ def test_proportional_costs(snapshot):
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

@ -13,3 +13,28 @@ def test_sort_dipoles_by_freq(snapshot):
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

@ -14,12 +14,12 @@
# name: test_fast_spin_qubit_frequency_tarucha_calc
list([
list([
0.24485375860291592,
0.009068657726033923,
1.1471308805198364,
0.042486328908142086,
]),
list([
0.0001499877991005419,
0.00011239982757520363,
0.0008292459978410822,
0.0006214312613006971,
]),
])
# ---

View File

@ -1,9 +1,24 @@
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]
@ -16,21 +31,11 @@ def test_fast_nonlocal_calc():
[[[-1, -2, -3, 11], [-1, 2, 5, 11]], [[-1, -2, -3, 6], [2, 4, 6, 6]]]
)
# expected_ij is for pair i, dipole j
expected_11 = 0.000021124454334947546213
expected_12 = 0.000022184755131682365135
expected_13 = 0.0000053860643617855849275
expected_14 = -0.0000023069501696755220764
expected_21 = 0.00022356021100884617796
expected_22 = 0.00021717277640859343002
expected_23 = 0.000017558321044891869169
expected_24 = -0.000034714318479634499683
expected = numpy.array(
[
[expected_11, expected_21],
[expected_12, expected_22],
[expected_13, expected_23],
[expected_14, expected_24],
[s_potential_from_arrays(dipole_array, dot_pair) for dot_pair in dot_pairs]
for dipole_array in dipoles
]
)

View File

@ -1,11 +1,26 @@
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]
@ -18,19 +33,19 @@ def test_fast_nonlocal_calc_multidipole():
[[[-1, -2, -3, 11], [-1, 2, 5, 11]], [[-1, -2, -3, 6], [2, 4, 6, 6]]]
)
# expected_ij is for pair i, dipole j
expected_11 = 0.000021124454334947546213
expected_12 = 0.000022184755131682365135
expected_13 = 0.0000053860643617855849275
expected_14 = -0.0000023069501696755220764
expected_21 = 0.00022356021100884617796
expected_22 = 0.00021717277640859343002
expected_23 = 0.000017558321044891869169
expected_24 = -0.000034714318479634499683
expected = numpy.array(
[
[expected_11 + expected_12, expected_21 + expected_22],
[expected_13 + expected_14, expected_23 + expected_24],
[
sum(
[
s_potential_from_arrays(dipole_array, dot_pair)
for dipole_array in dipoles
]
)
for dot_pair in dot_pairs
]
for dipoles in dipoleses
]
)

View File

@ -1,6 +1,30 @@
import numpy
import pdme.util.fast_v_calc
import pdme.measurement
def dipole_from_array(arr: numpy.ndarray) -> pdme.measurement.OscillatingDipole:
return pdme.measurement.OscillatingDipole(arr[0:3], arr[3:6], arr[6])
def s_potential_from_arrays(
dipole_array: numpy.ndarray, dotf_array: numpy.ndarray
) -> float:
dipole = dipole_from_array(dipole_array)
r = dotf_array[0:3]
f = dotf_array[3]
return dipole.s_electric_potential_at_position(r, f)
def s_electric_field_x_from_arrays(
dipole_array: numpy.ndarray, dotf_array: numpy.ndarray
) -> float:
dipole = dipole_from_array(dipole_array)
r = dotf_array[0:3]
f = dotf_array[3]
return dipole.s_electric_fieldx_at_position(r, f)
def test_fast_v_calc():
d1 = [1, 2, 3, 4, 5, 6, 7]
@ -9,11 +33,11 @@ def test_fast_v_calc():
dipoles = numpy.array([d1, d2])
dot_inputs = numpy.array([[-1, -1, -1, 11], [2, 3, 1, 5.5]])
# expected_ij is for dot i, dipole j
expected_11 = 0.00001421963287022476
expected_12 = 0.00001107180225755457
expected_21 = 0.000345021108583681380388722
expected_22 = 0.0000377061050587914705139781
expected_11 = s_potential_from_arrays(dipoles[0], dot_inputs[0])
expected_12 = s_potential_from_arrays(dipoles[1], dot_inputs[0])
expected_21 = s_potential_from_arrays(dipoles[0], dot_inputs[1])
expected_22 = s_potential_from_arrays(dipoles[1], dot_inputs[1])
expected = numpy.array([[expected_11, expected_21], [expected_12, expected_22]])
@ -31,11 +55,11 @@ def test_fast_v_calc_multidipoles():
dipoles = numpy.array([[d1, d2]])
dot_inputs = numpy.array([[-1, -1, -1, 11], [2, 3, 1, 5.5]])
# expected_ij is for dot i, dipole j
expected_11 = 0.00001421963287022476
expected_12 = 0.00001107180225755457
expected_21 = 0.000345021108583681380388722
expected_22 = 0.0000377061050587914705139781
expected_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]])
@ -48,7 +72,7 @@ def test_fast_v_calc_multidipoles():
def test_fast_v_calc_big_multidipole():
dipoles = numpy.array(
dipoleses = numpy.array(
[
[
[1, 1, 5, 6, 3, 1, 1],
@ -72,27 +96,113 @@ def test_fast_v_calc_big_multidipole():
]
)
expected = numpy.array(
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(
[
[
0.0010151687742365581690202135,
0.00077627527320628609782627266,
0.00043313471258511003340648713,
0.000077184305988088453637005111,
[1, 1, 5, 6, 3, 1, 1],
[5, 3, 2, 13, 1, 1, 2],
[-5, -5, -3, -1, -3, 8, 3],
],
[
0.000041099091967966890657097060,
0.0019377687238977568792327845,
0.0085903193415282984161225029,
0.00014557676715208209310911838,
[-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_vs_for_dipoleses(dot_inputs, dipoles),
pdme.util.fast_v_calc.fast_efieldxs_for_dipoleses(dot_inputs, dipoleses),
expected,
err_msg="Voltages at dot aren't as expected for multidipole calc.",
err_msg="E x fast calc at dot aren't as expected for multidipole calc.",
)
@ -107,28 +217,3 @@ def test_between():
expected = numpy.array([False, False, True])
numpy.testing.assert_array_equal(actual, expected, err_msg="Between calc wrong")
def test_fast_v_calc_asymmetric_multidipoles_but_symmetric():
# expected format is [px, py, pz, sx, sy, sz, e1, e2, w]
d1 = [1, 2, 3, 4, 5, 6, 1, 1, 7 / 2]
d2 = [2, 5, 3, 4, -5, -6, 2, 2, 2 / 2]
dipoles = numpy.array([[d1, d2]])
dot_inputs = numpy.array([[-1, -1, -1, 11], [2, 3, 1, 5.5]])
# expected_ij is for dot i, dipole j
expected_11 = 0.00001421963287022476
expected_12 = 0.00001107180225755457
expected_21 = 0.000345021108583681380388722
expected_22 = 0.0000377061050587914705139781
expected = numpy.array([[expected_11 + expected_12, expected_21 + expected_22]])
numpy.testing.assert_allclose(
pdme.util.fast_v_calc.fast_vs_for_asymmetric_dipoleses(
dot_inputs, dipoles, 1e10
),
expected,
err_msg="Voltages at dot aren't as expected for multidipole calc.",
)