chore(release): 1.5.0

CHANGELOG.md

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

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

pyproject.toml

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

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)