diff --git a/pdme/calculations/__init__.py b/pdme/calculations/__init__.py new file mode 100644 index 0000000..a03bbc2 --- /dev/null +++ b/pdme/calculations/__init__.py @@ -0,0 +1,46 @@ +""" +This module is a canonical source of the accurate expressions we want to use for calculating our noise. +No reference to class or anything, just a straight set of math functions. +""" +import numpy + + +def telegraph_beta(f: float, w: float) -> float: + """ + This function represents the frequency component of analytic telegraph noise. + + We're assuming we care about the one-sided PSD where we are ignoring negative frequencies. + This matches with experimental data from say Connors et al., and I think is better than keeping with one-sided. + Note that this means that it will only be comparable then with time series data assuming one-sided! + + Don't bikeshed yet, if we care about two-sided things for any reason down the line divide this by two or just change it then. + """ + return 2 * w / ((numpy.pi * f) ** 2 + w**2) + + +def electric_potential(p: numpy.ndarray, s: numpy.ndarray, r: numpy.ndarray) -> float: + """ + Gives the electric potential of a defect with dipole moment p, located at position s, + as measured from position r. + + p, s, r, are numpy arrays of length 3 + """ + diff = r - s + return (p.dot(diff) / (numpy.linalg.norm(diff) ** 3)).item() + + +def electric_field( + p: numpy.ndarray, s: numpy.ndarray, r: numpy.ndarray +) -> numpy.ndarray: + """ + Gives the electric field of a defect with dipole moment p, located at position s, + as measured from position r. + + p, s, r, are numpy arrays of length 3 + + Returns an array of length 3, ideally. + """ + diff = r - s + norm_diff = numpy.linalg.norm(diff) + + return ((3 * (p.dot(diff) * diff) / (norm_diff**2)) - p) / (norm_diff**3) diff --git a/pdme/measurement/oscillating_dipole.py b/pdme/measurement/oscillating_dipole.py index 5f356e3..741e398 100644 --- a/pdme/measurement/oscillating_dipole.py +++ b/pdme/measurement/oscillating_dipole.py @@ -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 diff --git a/pdme/util/fast_nonlocal_spectrum.py b/pdme/util/fast_nonlocal_spectrum.py index 39160f9..6636584 100644 --- a/pdme/util/fast_nonlocal_spectrum.py +++ b/pdme/util/fast_nonlocal_spectrum.py @@ -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) @@ -218,7 +252,7 @@ def fast_s_spin_qubit_tarucha_nonlocal_dipoleses( _logger.debug(f"alphses1: {alphses1}") _logger.debug(f"alphses2: {alphses2}") - bses = (1 / numpy.pi) * (ws[:, None, :] / (f1s[:, None] ** 2 + ws[:, None, :] ** 2)) + bses = 2 * ws[:, None, :] / ((numpy.pi * f1s[:, None]) ** 2 + ws[:, None, :] ** 2) if _logger.isEnabledFor(logging.DEBUG): _logger.debug(f"bses: {bses}") diff --git a/pdme/util/fast_v_calc.py b/pdme/util/fast_v_calc.py index 833db57..0f8351c 100644 --- a/pdme/util/fast_v_calc.py +++ b/pdme/util/fast_v_calc.py @@ -31,7 +31,8 @@ def fast_vs_for_dipoles( 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)) + bses = 2 * ws[:, None] / ((numpy.pi * fs) ** 2 + ws[:, None] ** 2) + _logger.debug(f"bses: {bses}") return ases * bses @@ -64,7 +65,7 @@ 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) @@ -92,6 +93,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 +109,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) diff --git a/tests/calculations/__snapshots__/test_calculations.ambr b/tests/calculations/__snapshots__/test_calculations.ambr new file mode 100644 index 0000000..5a3ba84 --- /dev/null +++ b/tests/calculations/__snapshots__/test_calculations.ambr @@ -0,0 +1,81 @@ +# serializer version: 1 +# name: test_multiple_electric_field_alphas + list([ + list([ + 0.009579434215333742, + 0.007714411624737791, + 0.0035604976729559034, + ]), + list([ + 0.005036717012002481, + 0.006259221141129298, + -0.009682232702684382, + ]), + list([ + 0.01503516326014651, + 0.012028130608117207, + 0.009021097956087907, + ]), + list([ + 0.0033871215792458027, + 0.0038182097802407235, + -0.005604146612933966, + ]), + list([ + 0.007409666906433089, + 0.006778734778841714, + 0.003467602973242188, + ]), + list([ + 0.004730939083250055, + 0.0046546336141653774, + -0.011766303332857395, + ]), + list([ + -0.010766266772985707, + -0.012631289363581657, + -0.010851040527103707, + ]), + list([ + -0.008410828408392494, + -0.020391368873835285, + 0.009877833363344673, + ]), + list([ + -0.015035163260146511, + -0.018042195912175815, + -0.021049228564205113, + ]), + list([ + -0.005850482727788205, + -0.012193637685284892, + 0.0054809785555068454, + ]), + list([ + -0.007682229585552562, + -0.010584517372472879, + -0.009352937859414517, + ]), + list([ + -0.009767100042838822, + -0.020831393060117175, + 0.014024945217763873, + ]), + ]) +# --- +# name: test_multiple_electric_potential_alphas + list([ + 0.05035560994609065, + -0.03369221379873504, + 0.07216878364870323, + -0.024633611485424027, + 0.04496779459769236, + -0.03108684810509794, + -0.08019597139562586, + 0.08293468011996319, + -0.10825317547305485, + 0.060044427995721066, + -0.06935710692186449, + 0.08733923991432278, + ]) +# --- diff --git a/tests/calculations/test_calculations.py b/tests/calculations/test_calculations.py new file mode 100644 index 0000000..6a7267f --- /dev/null +++ b/tests/calculations/test_calculations.py @@ -0,0 +1,102 @@ +import pdme.calculations +import pytest +import numpy +import numpy.testing + + +# generated in mathematica to compare here +beta_test_data = [ + [-2, -2, 0.8072976151], + [-1, -2, 0.008105366193], + [0, -2, 0.00008105691406], + [1, -2, 8.105694659e-7], + [2, -2, 8.105694691e-9], + [-2, -1, 5.768008783], + [-1, -1, 0.08072976151], + [0, -1, 0.0008105366193], + [1, -1, 8.105691406e-6], + [2, -1, 8.105694659e-8], + [-2, 0, 1.951840272], + [-1, 0, 0.5768008783], + [0, 0, 0.008072976151], + [1, 0, 0.00008105366193], + [2, 0, 8.105691406e-7], + [-2, 1, 0.1999506642], + [-1, 1, 0.1951840272], + [0, 1, 0.05768008783], + [1, 1, 0.0008072976151], + [2, 1, 8.105366193e-6], + [-2, 2, 0.01999995065], + [-1, 2, 0.01999506642], + [0, 2, 0.01951840272], + [1, 2, 0.005768008783], + [2, 2, 0.00008072976151], +] + + +@pytest.mark.parametrize("f, w, expected", beta_test_data) +def test_calculations_beta_func(f, w, expected): + # there's nothing special about the 5 * 10^ f and 10^w passing in logs + # this was just to get a variety of orders of magnitude in results. + actual = pdme.calculations.telegraph_beta(5 * 10**f, 10**w) + numpy.testing.assert_allclose(actual, expected, atol=0, rtol=1e-8) + + +def test_multiple_electric_potential_alphas(snapshot): + """ + Manually compare these with mathematica stuff because manually including a list is a bit annoying. + Basically just visually compare the snapshot values to make sure they're actually correct. + """ + dipole_ps = [ + numpy.array([1, 2, 3]), + numpy.array([-4, -3, -2]), + ] + dipole_ss = [ + numpy.array([0, 0, 1]), + numpy.array([1, 1, 1]), + numpy.array([0, -0.5, 0.5]), + ] + + test_rs = [ + numpy.array([5, 5, 5]), + numpy.array([-4, -6, 2]), + ] + + actuals = [ + pdme.calculations.electric_potential(p, s, r) + for p in dipole_ps + for s in dipole_ss + for r in test_rs + ] + + assert actuals == snapshot + + +def test_multiple_electric_field_alphas(snapshot): + """ + Manually compare these with mathematica stuff because manually including a list is a bit annoying. + Basically just visually compare the snapshot values to make sure they're actually correct. + """ + dipole_ps = [ + numpy.array([1, 2, 3]), + numpy.array([-4, -3, -2]), + ] + dipole_ss = [ + numpy.array([0, 0, 1]), + numpy.array([1, 1, 1]), + numpy.array([0, -0.5, 0.5]), + ] + + test_rs = [ + numpy.array([5, 5, 5]), + numpy.array([-4, -6, 2]), + ] + + actuals = [ + pdme.calculations.electric_field(p, s, r).tolist() + for p in dipole_ps + for s in dipole_ss + for r in test_rs + ] + + assert actuals == snapshot diff --git a/tests/measurement/test_dipole_noise_measurement.py b/tests/measurement/test_dipole_noise_measurement.py index 3e0ff13..0c85fdb 100644 --- a/tests/measurement/test_dipole_noise_measurement.py +++ b/tests/measurement/test_dipole_noise_measurement.py @@ -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." diff --git a/tests/model/__snapshots__/test_log_spaced_fixed_orientation_fixedxy_orientation_mcmc.ambr b/tests/model/__snapshots__/test_log_spaced_fixed_orientation_fixedxy_orientation_mcmc.ambr index 5f34e56..5bd485c 100644 --- a/tests/model/__snapshots__/test_log_spaced_fixed_orientation_fixedxy_orientation_mcmc.ambr +++ b/tests/model/__snapshots__/test_log_spaced_fixed_orientation_fixedxy_orientation_mcmc.ambr @@ -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]]), ), diff --git a/tests/model/__snapshots__/test_log_spaced_fixed_orientation_free_orientation_mcmc.ambr b/tests/model/__snapshots__/test_log_spaced_fixed_orientation_free_orientation_mcmc.ambr index e6030e7..d039e5d 100644 --- a/tests/model/__snapshots__/test_log_spaced_fixed_orientation_free_orientation_mcmc.ambr +++ b/tests/model/__snapshots__/test_log_spaced_fixed_orientation_free_orientation_mcmc.ambr @@ -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]]), ), diff --git a/tests/model/__snapshots__/test_log_spaced_fixed_orientation_mcmc.ambr b/tests/model/__snapshots__/test_log_spaced_fixed_orientation_mcmc.ambr index 4b01374..09c6659 100644 --- a/tests/model/__snapshots__/test_log_spaced_fixed_orientation_mcmc.ambr +++ b/tests/model/__snapshots__/test_log_spaced_fixed_orientation_mcmc.ambr @@ -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]]), ), diff --git a/tests/model/test_log_spaced_fixed_orientation_fixedxy_orientation_mcmc.py b/tests/model/test_log_spaced_fixed_orientation_fixedxy_orientation_mcmc.py index b0eb6bc..7d518fa 100644 --- a/tests/model/test_log_spaced_fixed_orientation_fixedxy_orientation_mcmc.py +++ b/tests/model/test_log_spaced_fixed_orientation_fixedxy_orientation_mcmc.py @@ -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: diff --git a/tests/model/test_log_spaced_fixed_orientation_free_orientation_mcmc.py b/tests/model/test_log_spaced_fixed_orientation_free_orientation_mcmc.py index ea84ac0..18178a7 100644 --- a/tests/model/test_log_spaced_fixed_orientation_free_orientation_mcmc.py +++ b/tests/model/test_log_spaced_fixed_orientation_free_orientation_mcmc.py @@ -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: diff --git a/tests/model/test_log_spaced_fixed_orientation_mcmc.py b/tests/model/test_log_spaced_fixed_orientation_mcmc.py index 843937f..269685f 100644 --- a/tests/model/test_log_spaced_fixed_orientation_mcmc.py +++ b/tests/model/test_log_spaced_fixed_orientation_mcmc.py @@ -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: diff --git a/tests/util/__snapshots__/test_fast_nonlocal_spectrum_pairs.ambr b/tests/util/__snapshots__/test_fast_nonlocal_spectrum_pairs.ambr index 114fa2d..8daefbe 100644 --- a/tests/util/__snapshots__/test_fast_nonlocal_spectrum_pairs.ambr +++ b/tests/util/__snapshots__/test_fast_nonlocal_spectrum_pairs.ambr @@ -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, ]), ]) # --- diff --git a/tests/util/test_fast_nonlocal_spectrum.py b/tests/util/test_fast_nonlocal_spectrum.py index 2f1a6c1..a671c18 100644 --- a/tests/util/test_fast_nonlocal_spectrum.py +++ b/tests/util/test_fast_nonlocal_spectrum.py @@ -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 ] ) diff --git a/tests/util/test_fast_nonlocal_spectrum_pairs.py b/tests/util/test_fast_nonlocal_spectrum_pairs.py index 5c71d9d..19dbc37 100644 --- a/tests/util/test_fast_nonlocal_spectrum_pairs.py +++ b/tests/util/test_fast_nonlocal_spectrum_pairs.py @@ -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 ] ) diff --git a/tests/util/test_fast_v_calc.py b/tests/util/test_fast_v_calc.py index b5c44d9..b14b891 100644 --- a/tests/util/test_fast_v_calc.py +++ b/tests/util/test_fast_v_calc.py @@ -1,6 +1,21 @@ 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 test_fast_v_calc(): d1 = [1, 2, 3, 4, 5, 6, 7] @@ -9,11 +24,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 +46,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 +63,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,25 +87,21 @@ def test_fast_v_calc_big_multidipole(): ] ) - expected = numpy.array( + expected = [ [ - [ - 0.0010151687742365581690202135, - 0.00077627527320628609782627266, - 0.00043313471258511003340648713, - 0.000077184305988088453637005111, - ], - [ - 0.000041099091967966890657097060, - 0.0019377687238977568792327845, - 0.0085903193415282984161225029, - 0.00014557676715208209310911838, - ], + 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, dipoles), + 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.", ) @@ -107,28 +118,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.", - )