chore(release): 1.5.0

CHANGELOG.md

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

pdme/model/model.py

@@ -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)

pdme/subspace_simulation/__init__.py

@@ -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",

pdme/subspace_simulation/mcmc_costs.py

@@ -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)

pdme/util/fast_nonlocal_spectrum.py

@@ -252,10 +252,15 @@ def fast_s_spin_qubit_tarucha_nonlocal_dipoleses(
 		_logger.debug(f"alphses1: {alphses1}")
 		_logger.debug(f"alphses2: {alphses2}")
 
+	# ws has shape (A, j), so it becomes (A, 1, j) in numerator with the new axis
+	# f1s has shape (m), so we get in the denominator (m, 1) + (A, 1, j)
+	# This becomes (A, m, j)
 	bses = 2 * ws[:, None, :] / ((numpy.pi * f1s[:, None]) ** 2 + ws[:, None, :] ** 2)
 	if _logger.isEnabledFor(logging.DEBUG):
 		_logger.debug(f"bses: {bses}")
 
+	# alphas have (A, 1, j), betas have (A, m, j)
+	# Final result is (A, m, j)
 	_logger.debug(f"Raw pair calc: [{alphses1 * alphses2 * bses}]")
 	return numpy.einsum("...j->...", alphses1 * alphses2 * bses)
 

pdme/util/fast_v_calc.py

@@ -12,28 +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}")
 
+	# [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
 
 
@@ -70,6 +86,121 @@ def fast_vs_for_dipoleses(
 	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:

pyproject.toml

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

tests/subspace_simulation/__snapshots__/test_sort_dipoles.ambr

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

tests/subspace_simulation/test_costs.py

@@ -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)

tests/subspace_simulation/test_sort_dipoles.py

@@ -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

tests/util/test_fast_v_calc.py

@@ -17,6 +17,15 @@ def s_potential_from_arrays(
 	return dipole.s_electric_potential_at_position(r, f)
 
 
+def s_electric_field_x_from_arrays(
+	dipole_array: numpy.ndarray, dotf_array: numpy.ndarray
+) -> float:
+	dipole = dipole_from_array(dipole_array)
+	r = dotf_array[0:3]
+	f = dotf_array[3]
+	return dipole.s_electric_fieldx_at_position(r, f)
+
+
 def test_fast_v_calc():
 	d1 = [1, 2, 3, 4, 5, 6, 7]
 	d2 = [2, 5, 3, 4, -5, -6, 2]
@@ -107,6 +116,96 @@ def test_fast_v_calc_big_multidipole():
 	)
 
 
+def test_fast_electric_field_x_calc():
+	d1 = [1, 2, 3, 4, 5, 6, 7]
+	d2 = [2, 5, 3, 4, -5, -6, 2]
+
+	dipoles = numpy.array([d1, d2])
+
+	dot_inputs = numpy.array([[-1, -1, -1, 11], [2, 3, 1, 5.5]])
+
+	expected_11 = s_electric_field_x_from_arrays(dipoles[0], dot_inputs[0])
+	expected_12 = s_electric_field_x_from_arrays(dipoles[1], dot_inputs[0])
+	expected_21 = s_electric_field_x_from_arrays(dipoles[0], dot_inputs[1])
+	expected_22 = s_electric_field_x_from_arrays(dipoles[1], dot_inputs[1])
+
+	expected = numpy.array([[expected_11, expected_21], [expected_12, expected_22]])
+
+	numpy.testing.assert_allclose(
+		pdme.util.fast_v_calc.fast_efieldxs_for_dipoles(dot_inputs, dipoles),
+		expected,
+		err_msg="E x fast calc at dot aren't as expected.",
+	)
+
+
+def test_fast_electric_field_x_calc_multidipoles():
+	d1 = [1, 2, 3, 4, 5, 6, 7]
+	d2 = [2, 5, 3, 4, -5, -6, 2]
+
+	dipoles = numpy.array([[d1, d2]])
+
+	dot_inputs = numpy.array([[-1, -1, -1, 11], [2, 3, 1, 5.5]])
+
+	expected_11 = s_electric_field_x_from_arrays(dipoles[0][0], dot_inputs[0])
+	expected_12 = s_electric_field_x_from_arrays(dipoles[0][1], dot_inputs[0])
+	expected_21 = s_electric_field_x_from_arrays(dipoles[0][0], dot_inputs[1])
+	expected_22 = s_electric_field_x_from_arrays(dipoles[0][1], dot_inputs[1])
+
+	expected = numpy.array([[expected_11 + expected_12, expected_21 + expected_22]])
+
+	numpy.testing.assert_allclose(
+		pdme.util.fast_v_calc.fast_efieldxs_for_dipoleses(dot_inputs, dipoles),
+		expected,
+		err_msg="E x fast calc at dot aren't as expected for multidipole calc.",
+	)
+
+
+def test_fast_electric_field_x_calc_big_multidipole():
+
+	dipoleses = numpy.array(
+		[
+			[
+				[1, 1, 5, 6, 3, 1, 1],
+				[5, 3, 2, 13, 1, 1, 2],
+				[-5, -5, -3, -1, -3, 8, 3],
+			],
+			[
+				[-3, -1, -2, -2, -6, 3, 4],
+				[8, 0, 2, 0, 1, 5, 5],
+				[1, 4, -4, -1, -3, -5, 6],
+			],
+		]
+	)
+
+	dot_inputs = numpy.array(
+		[
+			[1, 1, 0, 1],
+			[2, 5, 6, 2],
+			[3, 1, 3, 3],
+			[0.5, 0.5, 0.5, 4],
+		]
+	)
+
+	expected = [
+		[
+			sum(
+				[
+					s_electric_field_x_from_arrays(dipole_array, dot_input)
+					for dipole_array in dipole_config
+				]
+			)
+			for dot_input in dot_inputs
+		]
+		for dipole_config in dipoleses
+	]
+
+	numpy.testing.assert_allclose(
+		pdme.util.fast_v_calc.fast_efieldxs_for_dipoleses(dot_inputs, dipoleses),
+		expected,
+		err_msg="E x fast calc at dot aren't as expected for multidipole calc.",
+	)
+
+
 def test_between():
 	low = numpy.array([1, 2, 3])
 	high = numpy.array([6, 7, 8])