From 55373cc6db1d8332d01dd2787fc93eb97b8cd1ad Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Thu, 9 May 2024 21:26:44 +0100 Subject: [PATCH 01/10] adds initial cuckoo implementation and corresponding tests, fixes the parameterisation integration tests --- .../notebooks/optimiser_calibration.ipynb | 218 ++++++++++++------ pybop/__init__.py | 1 + pybop/optimisers/_cuckoo.py | 201 ++++++++++++++++ tests/integration/test_parameterisations.py | 51 ++-- tests/unit/test_optimisation.py | 1 + 5 files changed, 386 insertions(+), 86 deletions(-) create mode 100644 pybop/optimisers/_cuckoo.py diff --git a/examples/notebooks/optimiser_calibration.ipynb b/examples/notebooks/optimiser_calibration.ipynb index beed72875..61ff509be 100644 --- a/examples/notebooks/optimiser_calibration.ipynb +++ b/examples/notebooks/optimiser_calibration.ipynb @@ -126,9 +126,24 @@ "outputs": [], "source": [ "parameter_set = pybop.ParameterSet.pybamm(\"Chen2020\")\n", + "parameter_set.update(\n", + " {\n", + " \"Negative electrode active material volume fraction\": 0.65,\n", + " \"Positive electrode active material volume fraction\": 0.51,\n", + " }\n", + ")\n", "model = pybop.lithium_ion.SPM(parameter_set=parameter_set)\n", - "t_eval = np.arange(0, 900, 3)\n", - "values = model.predict(t_eval=t_eval)" + "init_soc = 0.4\n", + "experiment = pybop.Experiment(\n", + " [\n", + " (\n", + " \"Discharge at 0.5C for 6 minutes (4 second period)\",\n", + " \"Charge at 0.5C for 6 minutes (4 second period)\",\n", + " ),\n", + " ]\n", + " * 2\n", + ")\n", + "values = model.predict(init_soc=init_soc, experiment=experiment)" ] }, { @@ -153,8 +168,10 @@ }, "outputs": [], "source": [ - "sigma = 0.001\n", - "corrupt_values = values[\"Voltage [V]\"].data + np.random.normal(0, sigma, len(t_eval))" + "sigma = 0.002\n", + "corrupt_values = values[\"Voltage [V]\"].data + np.random.normal(\n", + " 0, sigma, len(values[\"Voltage [V]\"].data)\n", + ")" ] }, { @@ -200,7 +217,7 @@ "source": [ "dataset = pybop.Dataset(\n", " {\n", - " \"Time [s]\": t_eval,\n", + " \"Time [s]\": values[\"Time [s]\"].data,\n", " \"Current function [A]\": values[\"Current [A]\"].data,\n", " \"Voltage [V]\": corrupt_values,\n", " }\n", @@ -235,13 +252,15 @@ "parameters = [\n", " pybop.Parameter(\n", " \"Negative electrode active material volume fraction\",\n", - " prior=pybop.Gaussian(0.7, 0.025),\n", - " bounds=[0.6, 0.9],\n", + " prior=pybop.Uniform(0.45, 0.7),\n", + " bounds=[0.4, 0.8],\n", + " true_value=0.65,\n", " ),\n", " pybop.Parameter(\n", " \"Positive electrode active material volume fraction\",\n", - " prior=pybop.Gaussian(0.6, 0.025),\n", - " bounds=[0.5, 0.8],\n", + " prior=pybop.Uniform(0.45, 0.7),\n", + " bounds=[0.4, 0.8],\n", + " true_value=0.51,\n", " ),\n", "]" ] @@ -279,7 +298,7 @@ } ], "source": [ - "problem = pybop.FittingProblem(model, parameters, dataset)\n", + "problem = pybop.FittingProblem(model, parameters, dataset, init_soc=init_soc)\n", "cost = pybop.SumSquaredError(problem)\n", "optim = pybop.Optimisation(cost, optimiser=pybop.GradientDescent, sigma0=0.2)\n", "optim.set_max_iterations(100)" @@ -308,26 +327,7 @@ }, "id": "-9OVt0EQ04qB" }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Error: Events ['Maximum voltage [V]'] are non-positive at initial conditions\n", - "Error: Events ['Maximum voltage [V]'] are non-positive at initial conditions\n", - "Error: Events ['Maximum voltage [V]'] are non-positive at initial conditions\n", - "Error: Events ['Maximum voltage [V]'] are non-positive at initial conditions\n", - "Error: Events ['Maximum voltage [V]'] are non-positive at initial conditions\n", - "Error: Events ['Maximum voltage [V]'] are non-positive at initial conditions\n", - "Error: Events ['Maximum voltage [V]'] are non-positive at initial conditions\n", - "Error: Events ['Maximum voltage [V]'] are non-positive at initial conditions\n", - "Error: Events ['Maximum voltage [V]'] are non-positive at initial conditions\n", - "Error: Events ['Maximum voltage [V]'] are non-positive at initial conditions\n", - "Error: Events ['Maximum voltage [V]'] are non-positive at initial conditions\n", - "Error: Events ['Maximum voltage [V]'] are non-positive at initial conditions\n" - ] - } - ], + "outputs": [], "source": [ "x, final_cost = optim.run()" ] @@ -363,7 +363,7 @@ { "data": { "text/plain": [ - "array([0.70742414, 0.58383355])" + "array([0.64609807, 0.51472958])" ] }, "execution_count": 9, @@ -397,7 +397,7 @@ { "data": { "image/svg+xml": [ - "02004006008003.753.83.853.93.9544.05ReferenceModelOptimised ComparisonTime / sVoltage / V" + "0500100015003.53.553.63.653.7ReferenceModelOptimised ComparisonTime / sVoltage / V" ] }, "metadata": {}, @@ -435,26 +435,30 @@ "text": [ "0.001\n", "NOTE: Boundaries ignored by Gradient Descent\n", - "0.0045\n", + "0.012285714285714285\n", + "NOTE: Boundaries ignored by Gradient Descent\n", + "0.023571428571428573\n", + "NOTE: Boundaries ignored by Gradient Descent\n", + "0.03485714285714286\n", "NOTE: Boundaries ignored by Gradient Descent\n", - "0.008\n", + "0.046142857142857145\n", "NOTE: Boundaries ignored by Gradient Descent\n", - "0.0115\n", + "0.05742857142857143\n", "NOTE: Boundaries ignored by Gradient Descent\n", - "0.015\n", + "0.06871428571428571\n", + "NOTE: Boundaries ignored by Gradient Descent\n", + "0.08\n", "NOTE: Boundaries ignored by Gradient Descent\n" ] } ], "source": [ - "sigmas = np.linspace(\n", - " 0.001, 0.015, 5\n", - ") # Change this to a smaller range for a quicker run\n", + "sigmas = np.linspace(0.001, 0.08, 8) # Change this to a smaller range for a quicker run\n", "xs = []\n", "optims = []\n", "for sigma in sigmas:\n", " print(sigma)\n", - " problem = pybop.FittingProblem(model, parameters, dataset)\n", + " problem = pybop.FittingProblem(model, parameters, dataset, init_soc=init_soc)\n", " cost = pybop.SumSquaredError(problem)\n", " optim = pybop.Optimisation(cost, optimiser=pybop.GradientDescent, sigma0=sigma)\n", " optim.set_max_iterations(100)\n", @@ -479,11 +483,14 @@ "name": "stdout", "output_type": "stream", "text": [ - "| Sigma: 0.001 | Num Iterations: 100 | Best Cost: 0.0013289907848209911 | Results: [0.69535773 0.67509662] |\n", - "| Sigma: 0.0045 | Num Iterations: 100 | Best Cost: 0.0007218197918308683 | Results: [0.71892626 0.67060898] |\n", - "| Sigma: 0.008 | Num Iterations: 100 | Best Cost: 0.0006371022763628136 | Results: [0.72396797 0.6696914 ] |\n", - "| Sigma: 0.0115 | Num Iterations: 18 | Best Cost: 0.0004608694532019237 | Results: [0.74070995 0.6667419 ] |\n", - "| Sigma: 0.015 | Num Iterations: 100 | Best Cost: 0.0007468897676990436 | Results: [0.71758655 0.67085529] |\n" + "| Sigma: 0.001 | Num Iterations: 100 | Best Cost: 0.008590687346571011 | Results: [0.58273999 0.64430015] |\n", + "| Sigma: 0.012285714285714285 | Num Iterations: 100 | Best Cost: 0.0017482878947612424 | Results: [0.62229759 0.5406604 ] |\n", + "| Sigma: 0.023571428571428573 | Num Iterations: 100 | Best Cost: 0.0013871420979637958 | Results: [0.63941964 0.52140605] |\n", + "| Sigma: 0.03485714285714286 | Num Iterations: 100 | Best Cost: 0.001571369568098984 | Results: [0.62907481 0.53267599] |\n", + "| Sigma: 0.046142857142857145 | Num Iterations: 28 | Best Cost: 0.0013533853388748253 | Results: [0.64673791 0.51409832] |\n", + "| Sigma: 0.05742857142857143 | Num Iterations: 25 | Best Cost: 0.0013584031053821507 | Results: [0.64390064 0.51673076] |\n", + "| Sigma: 0.06871428571428571 | Num Iterations: 74 | Best Cost: 0.0013568172573032275 | Results: [0.64444354 0.51631924] |\n", + "| Sigma: 0.08 | Num Iterations: 73 | Best Cost: 0.0013551215844470215 | Results: [0.64505654 0.51551585] |\n" ] } ], @@ -516,7 +523,34 @@ { "data": { "image/svg+xml": [ - "2040608010000.050.10.150.2Sigma: 0.001IterationCost" + "204060801000.00850.0090.00950.010.01050.0110.01150.012Sigma: 0.001IterationCost" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/svg+xml": [ + "0204060800.5840.5860.5880.590.5920.5940204060800.6440.6460.6480.650.6520.6540.6560.658Negative electrode active material volume fractionPositive electrode active material volume fractionParameter ConvergenceFunction CallFunction CallNegative electrode active material volume fractionPositive electrode active material volume fraction" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/svg+xml": [ + "2040608010000.0050.010.0150.020.0250.030.035Sigma: 0.012285714285714285IterationCost" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/svg+xml": [ + "0204060800.540.560.580.60.620204060800.520.5250.530.5350.540.5450.550.555Negative electrode active material volume fractionPositive electrode active material volume fractionParameter ConvergenceFunction CallFunction CallNegative electrode active material volume fractionPositive electrode active material volume fraction" ] }, "metadata": {}, @@ -525,7 +559,7 @@ { "data": { "image/svg+xml": [ - "0204060800.680.6820.6840.6860.6880.690.6920.6940.6960204060800.610.620.630.640.650.660.67Negative electrode active material volume fractionPositive electrode active material volume fractionParameter ConvergenceFunction CallFunction CallNegative electrode active material volume fractionPositive electrode active material volume fraction" + "2040608010000.020.040.060.080.1Sigma: 0.023571428571428573IterationCost" ] }, "metadata": {}, @@ -534,7 +568,7 @@ { "data": { "image/svg+xml": [ - "2040608010000.050.10.150.20.25Sigma: 0.0045IterationCost" + "0204060800.50.520.540.560.580.60.620.640204060800.450.460.470.480.490.50.510.520.530.54Negative electrode active material volume fractionPositive electrode active material volume fractionParameter ConvergenceFunction CallFunction CallNegative electrode active material volume fractionPositive electrode active material volume fraction" ] }, "metadata": {}, @@ -543,7 +577,7 @@ { "data": { "image/svg+xml": [ - "0204060800.70.7050.710.7150.720204060800.60.610.620.630.640.650.660.67Negative electrode active material volume fractionPositive electrode active material volume fractionParameter ConvergenceFunction CallFunction CallNegative electrode active material volume fractionPositive electrode active material volume fraction" + "204060801000.0050.010.0150.02Sigma: 0.03485714285714286IterationCost" ] }, "metadata": {}, @@ -552,7 +586,7 @@ { "data": { "image/svg+xml": [ - "2040608010000.050.10.150.20.25Sigma: 0.008IterationCost" + "0204060800.570.580.590.60.610.620.630204060800.540.560.580.60.620.640.660.680.7Negative electrode active material volume fractionPositive electrode active material volume fractionParameter ConvergenceFunction CallFunction CallNegative electrode active material volume fractionPositive electrode active material volume fraction" ] }, "metadata": {}, @@ -561,7 +595,7 @@ { "data": { "image/svg+xml": [ - "0204060800.70.7050.710.7150.720.7250204060800.60.610.620.630.640.650.660.67Negative electrode active material volume fractionPositive electrode active material volume fractionParameter ConvergenceFunction CallFunction CallNegative electrode active material volume fractionPositive electrode active material volume fraction" + "5101520250.0020.0040.0060.0080.010.0120.014Sigma: 0.046142857142857145IterationCost" ] }, "metadata": {}, @@ -570,7 +604,7 @@ { "data": { "image/svg+xml": [ - "5101500.010.020.030.04Sigma: 0.0115IterationCost" + "05101520250.650.660.670.680.6905101520250.520.530.540.550.56Negative electrode active material volume fractionPositive electrode active material volume fractionParameter ConvergenceFunction CallFunction CallNegative electrode active material volume fractionPositive electrode active material volume fraction" ] }, "metadata": {}, @@ -579,7 +613,7 @@ { "data": { "image/svg+xml": [ - "0510150.7350.7360.7370.7380.7390.740.7410510150.6350.640.6450.650.6550.660.665Negative electrode active material volume fractionPositive electrode active material volume fractionParameter ConvergenceFunction CallFunction CallNegative electrode active material volume fractionPositive electrode active material volume fraction" + "5101520250.001350.00140.001450.00150.001550.00160.00165Sigma: 0.05742857142857143IterationCost" ] }, "metadata": {}, @@ -588,7 +622,7 @@ { "data": { "image/svg+xml": [ - "2040608010000.10.20.30.40.5Sigma: 0.015IterationCost" + "051015200.6340.6360.6380.640.6420.644051015200.5150.51550.5160.51650.5170.51750.5180.51850.5190.5195Negative electrode active material volume fractionPositive electrode active material volume fractionParameter ConvergenceFunction CallFunction CallNegative electrode active material volume fractionPositive electrode active material volume fraction" ] }, "metadata": {}, @@ -597,7 +631,34 @@ { "data": { "image/svg+xml": [ - "0204060800.660.670.680.690.70.710.720204060800.580.60.620.640.660.680.70.720.740.76Negative electrode active material volume fractionPositive electrode active material volume fractionParameter ConvergenceFunction CallFunction CallNegative electrode active material volume fractionPositive electrode active material volume fraction" + "1020304050607000.0050.010.0150.020.0250.030.0350.04Sigma: 0.06871428571428571IterationCost" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/svg+xml": [ + "02040600.610.620.630.640.650.660.670.680.6902040600.520.540.560.580.60.620.64Negative electrode active material volume fractionPositive electrode active material volume fractionParameter ConvergenceFunction CallFunction CallNegative electrode active material volume fractionPositive electrode active material volume fraction" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/svg+xml": [ + "102030405060700.0020.0040.0060.0080.01Sigma: 0.08IterationCost" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/svg+xml": [ + "02040600.560.570.580.590.60.610.620.630.640.6502040600.520.530.540.550.560.57Negative electrode active material volume fractionPositive electrode active material volume fractionParameter ConvergenceFunction CallFunction CallNegative electrode active material volume fractionPositive electrode active material volume fraction" ] }, "metadata": {}, @@ -634,7 +695,34 @@ { "data": { "image/svg+xml": [ - "0.60.650.70.750.80.850.90.50.550.60.650.70.750.80.40.81.21.622.4Sigma: 0.001Negative electrode active material volume fractionPositive electrode active material volume fraction" + "0.40.50.60.70.80.40.450.50.550.60.650.70.750.80.10.20.30.4Sigma: 0.001Negative electrode active material volume fractionPositive electrode active material volume fraction" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/svg+xml": [ + "0.40.50.60.70.80.40.450.50.550.60.650.70.750.80.10.20.30.4Sigma: 0.012285714285714285Negative electrode active material volume fractionPositive electrode active material volume fraction" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/svg+xml": [ + "0.40.50.60.70.80.40.450.50.550.60.650.70.750.80.10.20.30.4Sigma: 0.023571428571428573Negative electrode active material volume fractionPositive electrode active material volume fraction" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/svg+xml": [ + "0.40.50.60.70.80.40.450.50.550.60.650.70.750.80.10.20.30.4Sigma: 0.03485714285714286Negative electrode active material volume fractionPositive electrode active material volume fraction" ] }, "metadata": {}, @@ -643,7 +731,7 @@ { "data": { "image/svg+xml": [ - "0.60.650.70.750.80.850.90.50.550.60.650.70.750.80.40.81.21.622.4Sigma: 0.0045Negative electrode active material volume fractionPositive electrode active material volume fraction" + "0.40.50.60.70.80.40.450.50.550.60.650.70.750.80.10.20.30.4Sigma: 0.046142857142857145Negative electrode active material volume fractionPositive electrode active material volume fraction" ] }, "metadata": {}, @@ -652,7 +740,7 @@ { "data": { "image/svg+xml": [ - "0.60.650.70.750.80.850.90.50.550.60.650.70.750.80.40.81.21.622.4Sigma: 0.008Negative electrode active material volume fractionPositive electrode active material volume fraction" + "0.40.50.60.70.80.40.450.50.550.60.650.70.750.80.10.20.30.4Sigma: 0.05742857142857143Negative electrode active material volume fractionPositive electrode active material volume fraction" ] }, "metadata": {}, @@ -661,7 +749,7 @@ { "data": { "image/svg+xml": [ - "0.60.650.70.750.80.850.90.50.550.60.650.70.750.80.40.81.21.622.4Sigma: 0.0115Negative electrode active material volume fractionPositive electrode active material volume fraction" + "0.40.50.60.70.80.40.450.50.550.60.650.70.750.80.10.20.30.4Sigma: 0.06871428571428571Negative electrode active material volume fractionPositive electrode active material volume fraction" ] }, "metadata": {}, @@ -670,7 +758,7 @@ { "data": { "image/svg+xml": [ - "0.60.650.70.750.80.850.90.50.550.60.650.70.750.80.40.81.21.622.4Sigma: 0.015Negative electrode active material volume fractionPositive electrode active material volume fraction" + "0.40.50.60.70.80.40.450.50.550.60.650.70.750.80.10.20.30.4Sigma: 0.08Negative electrode active material volume fractionPositive electrode active material volume fraction" ] }, "metadata": {}, @@ -679,7 +767,7 @@ ], "source": [ "# Plot the cost landscape with optimisation path and updated bounds\n", - "bounds = np.array([[0.6, 0.9], [0.5, 0.8]])\n", + "bounds = np.array([[0.4, 0.8], [0.4, 0.8]])\n", "for optim, sigma in zip(optims, sigmas):\n", " pybop.plot2d(optim, bounds=bounds, steps=10, title=f\"Sigma: {sigma}\")" ] @@ -690,12 +778,12 @@ "source": [ "### Updating the Learning Rate\n", "\n", - "Let's take `sigma0 = 0.0115` as the best learning rate for this problem and look at the time-series trajectories." + "Let's take `sigma0 = 0.08` as the best learning rate for this problem and look at the time-series trajectories." ] }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 16, "metadata": { "execution": { "iopub.execute_input": "2024-04-14T18:59:54.698068Z", @@ -715,7 +803,7 @@ { "data": { "image/svg+xml": [ - "02004006008003.83.853.93.9544.05ReferenceModelOptimised ComparisonTime / sVoltage / V" + "0500100015003.53.553.63.653.7ReferenceModelOptimised ComparisonTime / sVoltage / V" ] }, "metadata": {}, @@ -723,7 +811,7 @@ } ], "source": [ - "optim = pybop.Optimisation(cost, optimiser=pybop.GradientDescent, sigma0=0.0115)\n", + "optim = pybop.Optimisation(cost, optimiser=pybop.GradientDescent, sigma0=0.08)\n", "x, final_cost = optim.run()\n", "pybop.quick_plot(problem, parameter_values=x, title=\"Optimised Comparison\");" ] diff --git a/pybop/__init__.py b/pybop/__init__.py index c25abdc75..7c208f344 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -113,6 +113,7 @@ SNES, XNES, ) +from .optimisers._cuckoo import CuckooSearch # # Parameter classes diff --git a/pybop/optimisers/_cuckoo.py b/pybop/optimisers/_cuckoo.py new file mode 100644 index 000000000..b4a6f97b8 --- /dev/null +++ b/pybop/optimisers/_cuckoo.py @@ -0,0 +1,201 @@ +import numpy as np +import pints +from scipy.special import gamma + + +class CuckooSearch(pints.PopulationBasedOptimiser): + """ + Cuckoo Search (CS) optimization algorithm, inspired by the brood parasitism + of some cuckoo species. This algorithm was introduced by Yang and Deb in 2009. + + The algorithm uses a population of host nests (solutions), where each cuckoo + (new solution) tries to replace a worse nest in the population. The quality + or fitness of the nests is determined by the objective function. A fraction + of the worst nests is abandoned at each generation, and new ones are built + randomly. + + The pseudo-code for the Cuckoo Search is as follows: + + 1. Initialize population of n host nests + 2. While (t < max_generations): + a. Get a cuckoo randomly by Lévy flights + b. Evaluate its quality/fitness F + c. Choose a nest among n (say, j) randomly + d. If (F > fitness of j): + i. Replace j with the new solution + e. Abandon a fraction (pa) of the worst nests and build new ones + f. Keep the best solutions/nests + g. Rank the solutions and find the current best + 3. End While + + This implementation also uses a decreasing step size for the Lévy flights, calculated + as sigma = sigma0 / sqrt(iterations), where sigma0 is the initial step size and + iterations is the current iteration number. + + Parameters: + - pa: Probability of discovering alien eggs/solutions (abandoning rate) + + References: + - X. -S. Yang and Suash Deb, "Cuckoo Search via Lévy flights," + 2009 World Congress on Nature & Biologically Inspired Computing (NaBIC), + Coimbatore, India, 2009, pp. 210-214, https://doi.org/10.1109/NABIC.2009.5393690. + + - S. Walton, O. Hassan, K. Morgan, M.R. Brown, + Modified cuckoo search: A new gradient free optimisation algorithm, + Chaos, Solitons & Fractals, Volume 44, Issue 9, 2011, + Pages 710-718, ISSN 0960-0779, + https://doi.org/10.1016/j.chaos.2011.06.004. + """ + + def __init__(self, x0, sigma0=0.01, bounds=None, pa=0.25): + if bounds is None: + self.boundaries = None + elif not all( + np.isfinite(value) for sublist in bounds.values() for value in sublist + ): + raise ValueError( + "Either all bounds or no bounds must be set for Cuckoo Search." + ) + else: + self.boundaries = pints.RectangularBoundaries( + bounds["lower"], bounds["upper"] + ) + super().__init__(x0, sigma0, self.boundaries) + + # Problem dimensionality + self._dim = len(x0) + + # Population size and abandon rate + self._n = self._population_size + self._pa = pa + self.step_size = self._sigma0 + self.beta = 1.5 + + # Set states + self._running = False + self._ready_for_tell = False + + # Initialise nests + if self._boundaries is not None: + self._nests = np.random.uniform( + low=self._boundaries.lower(), + high=self._boundaries.upper(), + size=(self._n, self._dim), + ) + else: + self._nests = np.random.normal(self._x0, self._sigma0) + + self._fitness = np.full(self._n, np.inf) + + # Initialise best solutions + self._x_best = np.copy(x0) + self._f_best = np.inf + + # Set iteration count + self._iterations = 1 + + def ask(self): + """ + Returns a list of next points in the parameter-space + to evaluate from the optimiser. + """ + # Set flag to indicate that the optimiser is ready to receive replies + self._ready_for_tell = True + self._running = True + + # Generate new solutions (cuckoos) by Lévy flights + self.step_size = self._sigma0 / np.sqrt(self._iterations) + step = self.levy_flight(self.beta, self._dim) * self.step_size + self.cuckoos = self._nests + step + return self.clip_nests(self.cuckoos) + + def tell(self, replies): + """ + Receives a list of function values from the cost function from points + previously specified by `self.ask()`, and updates the optimiser state + accordingly. + """ + # Update iteration count + self._iterations += 1 + + # Compare cuckoos with current nests + for i in range(self._n): + f_new = replies[i] + if f_new < self._fitness[i]: + self._nests[i] = self.cuckoos[i] + self._fitness[i] = f_new + if f_new < self._f_best: + self._f_best = f_new + self._x_best = self.cuckoos[i] + + # Abandon some worse nests + n_abandon = int(self._pa * self._n) + worst_nests = np.argsort(self._fitness)[-n_abandon:] + for idx in worst_nests: + if self._boundaries is not None: + self._nests[idx] = np.random.uniform( + low=self._boundaries.lower(), + high=self._boundaries.upper(), + size=self._dim, + ) + else: + self._nests[idx] = np.random.normal(self._x0, self._sigma0) + + self._fitness[idx] = np.inf # reset fitness + + def levy_flight(self, alpha, size): + """ + Generate step sizes via the Mantegna's algorithm for Levy flights + """ + from numpy import pi, power, random, sin + + sigma_u = power( + (gamma(1 + alpha) * sin(pi * alpha / 2)) + / (gamma((1 + alpha) / 2) * alpha * power(2, (alpha - 1) / 2)), + 1 / alpha, + ) + sigma_v = 1 + + u = random.normal(0, sigma_u, size=size) + v = random.normal(0, sigma_v, size=size) + step = u / power(abs(v), 1 / alpha) + + return step + + def clip_nests(self, x): + """ + Clip the input array to the boundaries. + """ + return np.clip(x, self._boundaries.lower(), self._boundaries.upper()) + + def _suggested_population_size(self): + """ + Inherited from Pints:PopulationBasedOptimiser. + Returns a suggested population size, based on the + dimension of the parameter space. + """ + return 4 + int(3 * np.log(self._n_parameters)) + + def running(self): + """ + Returns ``True`` if the optimisation is in progress. + """ + return self._running + + def x_best(self): + """ + Returns the best parameter values found so far. + """ + return self._x_best + + def f_best(self): + """ + Returns the best score found so far. + """ + return self._f_best + + def name(self): + """ + Returns the name of the optimiser. + """ + return "Cuckoo Search" diff --git a/tests/integration/test_parameterisations.py b/tests/integration/test_parameterisations.py index 32bd40335..4c69bb97c 100644 --- a/tests/integration/test_parameterisations.py +++ b/tests/integration/test_parameterisations.py @@ -25,12 +25,12 @@ def parameters(self): return [ pybop.Parameter( "Negative electrode active material volume fraction", - prior=pybop.Gaussian(0.55, 0.05), + prior=pybop.Uniform(0.35, 0.75), bounds=[0.375, 0.75], ), pybop.Parameter( "Positive electrode active material volume fraction", - prior=pybop.Gaussian(0.55, 0.05), + prior=pybop.Uniform(0.35, 0.75), # no bounds ), ] @@ -84,6 +84,7 @@ def spm_costs(self, model, parameters, cost_class, init_soc): pybop.SciPyDifferentialEvolution, pybop.Adam, pybop.CMAES, + pybop.CuckooSearch, pybop.GradientDescent, pybop.IRPropMin, pybop.NelderMead, @@ -95,7 +96,11 @@ def spm_costs(self, model, parameters, cost_class, init_soc): @pytest.mark.integration def test_spm_optimisers(self, optimiser, spm_costs): # Some optimisers require a complete set of bounds - if optimiser in [pybop.SciPyDifferentialEvolution, pybop.PSO]: + if optimiser in [ + pybop.SciPyDifferentialEvolution, + pybop.PSO, + pybop.CuckooSearch, + ]: spm_costs.problem.parameters[1].set_bounds( [0.3, 0.8] ) # Large range to ensure IC within bounds @@ -107,29 +112,33 @@ def test_spm_optimisers(self, optimiser, spm_costs): spm_costs.bounds = bounds # Test each optimiser - parameterisation = pybop.Optimisation( - cost=spm_costs, optimiser=optimiser, sigma0=0.05 - ) - parameterisation.set_max_unchanged_iterations(iterations=35, threshold=1e-5) - parameterisation.set_max_iterations(125) - - initial_cost = parameterisation.cost(spm_costs.x0) - if optimiser in [pybop.GradientDescent]: if isinstance( spm_costs, (pybop.GaussianLogLikelihoodKnownSigma, pybop.MAP) ): - parameterisation.optimiser.set_learning_rate(1.8e-5) + parameterisation = pybop.Optimisation( + cost=spm_costs, optimiser=optimiser, sigma0=5e-5 + ) else: - parameterisation.optimiser.set_learning_rate(0.015) - x, final_cost = parameterisation.run() - + parameterisation = pybop.Optimisation( + cost=spm_costs, optimiser=optimiser, sigma0=0.02 + ) elif optimiser in [pybop.SciPyMinimize]: - parameterisation.cost.problem.model.allow_infeasible_solutions = False - x, final_cost = parameterisation.run() - + parameterisation = pybop.Optimisation( + cost=spm_costs, + optimiser=optimiser, + sigma0=0.05, + allow_infeasible_solutions=False, + ) else: - x, final_cost = parameterisation.run() + parameterisation = pybop.Optimisation( + cost=spm_costs, optimiser=optimiser, sigma0=0.05 + ) + + parameterisation.set_max_unchanged_iterations(iterations=35, threshold=1e-5) + parameterisation.set_max_iterations(125) + initial_cost = parameterisation.cost(spm_costs.x0) + x, final_cost = parameterisation.run() # Assertions assert initial_cost > final_cost @@ -248,8 +257,8 @@ def getdata(self, model, x, init_soc): experiment = pybop.Experiment( [ ( - "Discharge at 0.5C for 3 minutes (2 second period)", - "Charge at 0.5C for 3 minutes (2 second period)", + "Discharge at 0.5C for 6 minutes (4 second period)", + "Charge at 0.5C for 6 minutes (4 second period)", ), ] * 2 diff --git a/tests/unit/test_optimisation.py b/tests/unit/test_optimisation.py index 54674c954..5e63eacb6 100644 --- a/tests/unit/test_optimisation.py +++ b/tests/unit/test_optimisation.py @@ -76,6 +76,7 @@ def two_param_cost(self, model, two_parameters, dataset): (pybop.GradientDescent, "Gradient descent"), (pybop.Adam, "Adam"), (pybop.CMAES, "Covariance Matrix Adaptation Evolution Strategy (CMA-ES)"), + (pybop.CuckooSearch, "Cuckoo Search"), (pybop.SNES, "Seperable Natural Evolution Strategy (SNES)"), (pybop.XNES, "Exponential Natural Evolution Strategy (xNES)"), (pybop.PSO, "Particle Swarm Optimisation (PSO)"), From d0f56f3be9d58ba488691f9c5bde5c5758a4f91b Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Thu, 23 May 2024 11:03:26 +0100 Subject: [PATCH 02/10] feat: integrate pybop optimisers with base_pints structure --- pybop/__init__.py | 3 ++- pybop/optimisers/_cuckoo.py | 20 ++++------------ pybop/optimisers/base_pints_optimiser.py | 20 ++++++++-------- pybop/optimisers/pints_optimisers.py | 30 +++++++++++++++++++++++- 4 files changed, 45 insertions(+), 28 deletions(-) diff --git a/pybop/__init__.py b/pybop/__init__.py index 44dbde658..fab124921 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -95,6 +95,7 @@ # # Optimiser class # +from .optimisers._cuckoo import _CuckooSearch from .optimisers.base_optimiser import BaseOptimiser from .optimisers.base_pints_optimiser import BasePintsOptimiser from .optimisers.scipy_optimisers import ( @@ -111,8 +112,8 @@ PSO, SNES, XNES, + CuckooSearch, ) -from .optimisers._cuckoo import CuckooSearch from .optimisers.optimisation import Optimisation # diff --git a/pybop/optimisers/_cuckoo.py b/pybop/optimisers/_cuckoo.py index b4a6f97b8..467c273bd 100644 --- a/pybop/optimisers/_cuckoo.py +++ b/pybop/optimisers/_cuckoo.py @@ -1,9 +1,9 @@ import numpy as np -import pints +from pints import PopulationBasedOptimiser from scipy.special import gamma -class CuckooSearch(pints.PopulationBasedOptimiser): +class _CuckooSearch(PopulationBasedOptimiser): """ Cuckoo Search (CS) optimization algorithm, inspired by the brood parasitism of some cuckoo species. This algorithm was introduced by Yang and Deb in 2009. @@ -47,20 +47,8 @@ class CuckooSearch(pints.PopulationBasedOptimiser): https://doi.org/10.1016/j.chaos.2011.06.004. """ - def __init__(self, x0, sigma0=0.01, bounds=None, pa=0.25): - if bounds is None: - self.boundaries = None - elif not all( - np.isfinite(value) for sublist in bounds.values() for value in sublist - ): - raise ValueError( - "Either all bounds or no bounds must be set for Cuckoo Search." - ) - else: - self.boundaries = pints.RectangularBoundaries( - bounds["lower"], bounds["upper"] - ) - super().__init__(x0, sigma0, self.boundaries) + def __init__(self, x0, sigma0=0.01, boundaries=None, pa=0.25): + super().__init__(x0, sigma0, boundaries=boundaries) # Problem dimensionality self._dim = len(x0) diff --git a/pybop/optimisers/base_pints_optimiser.py b/pybop/optimisers/base_pints_optimiser.py index 543d32b64..5a781305d 100644 --- a/pybop/optimisers/base_pints_optimiser.py +++ b/pybop/optimisers/base_pints_optimiser.py @@ -1,7 +1,7 @@ import numpy as np import pints -from pybop import BaseOptimiser +from pybop import BaseOptimiser, _CuckooSearch class BasePintsOptimiser(BaseOptimiser): @@ -131,16 +131,16 @@ def _sanitise_inputs(self): ): print(f"NOTE: Boundaries ignored by {self.pints_optimiser}") self.bounds = None - elif issubclass(self.pints_optimiser, pints.PSO): - if not all( - np.isfinite(value) - for sublist in self.bounds.values() - for value in sublist - ): - raise ValueError( - "Either all bounds or no bounds must be set for Pints PSO." - ) else: + if issubclass(self.pints_optimiser, (pints.PSO, _CuckooSearch)): + if not all( + np.isfinite(value) + for sublist in self.bounds.values() + for value in sublist + ): + raise ValueError( + f"Either all bounds or no bounds must be set for {self.pints_optimiser.__name__}." + ) self._boundaries = pints.RectangularBoundaries( self.bounds["lower"], self.bounds["upper"] ) diff --git a/pybop/optimisers/pints_optimisers.py b/pybop/optimisers/pints_optimisers.py index e3d8ee31b..ab0cc0616 100644 --- a/pybop/optimisers/pints_optimisers.py +++ b/pybop/optimisers/pints_optimisers.py @@ -1,6 +1,6 @@ import pints -from pybop import BasePintsOptimiser +from pybop import BasePintsOptimiser, _CuckooSearch class GradientDescent(BasePintsOptimiser): @@ -233,3 +233,31 @@ def __init__(self, cost, **optimiser_kwargs): + "Please choose another optimiser." ) super().__init__(cost, pints.CMAES, **optimiser_kwargs) + + +class CuckooSearch(BasePintsOptimiser): + """ + Adapter for the Cuckoo Search optimiser in PyBOP. + + Cuckoo Search is a population-based optimization algorithm inspired by the brood parasitism of some cuckoo species. + It is designed to be simple, efficient, and robust, and is suitable for global optimization problems. + + Parameters + ---------- + **optimiser_kwargs : optional + Valid PyBOP option keys and their values, for example: + x0 : array_like + Initial + sigma0 : float + Initial step size. + bounds : dict + A dictionary with 'lower' and 'upper' keys containing arrays for lower and + upper bounds on the parameters. + + See Also + -------- + pybop.CuckooSearch : PyBOP implementation of Cuckoo Search algorithm. + """ + + def __init__(self, cost, **optimiser_kwargs): + super().__init__(cost, _CuckooSearch, **optimiser_kwargs) From 7e1dd31aa08b902395f4d733556b34bbe3a732d6 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Fri, 7 Jun 2024 09:34:34 +0100 Subject: [PATCH 03/10] refactor: add abandon_nest method, bugfix BaseOptimiser boundaries construction --- pybop/optimisers/_cuckoo.py | 26 +++++++++++++++--------- pybop/optimisers/base_pints_optimiser.py | 6 +++--- 2 files changed, 19 insertions(+), 13 deletions(-) diff --git a/pybop/optimisers/_cuckoo.py b/pybop/optimisers/_cuckoo.py index 467c273bd..3a7b8522f 100644 --- a/pybop/optimisers/_cuckoo.py +++ b/pybop/optimisers/_cuckoo.py @@ -120,15 +120,7 @@ def tell(self, replies): n_abandon = int(self._pa * self._n) worst_nests = np.argsort(self._fitness)[-n_abandon:] for idx in worst_nests: - if self._boundaries is not None: - self._nests[idx] = np.random.uniform( - low=self._boundaries.lower(), - high=self._boundaries.upper(), - size=self._dim, - ) - else: - self._nests[idx] = np.random.normal(self._x0, self._sigma0) - + self.abandon_nests(idx) self._fitness[idx] = np.inf # reset fitness def levy_flight(self, alpha, size): @@ -150,11 +142,25 @@ def levy_flight(self, alpha, size): return step + def abandon_nests(self, idx): + """ + Set the boundaries for the parameter space. + """ + if self._boundaries is not None: + self._nests[idx] = np.random.uniform( + low=self._boundaries.lower(), + high=self._boundaries.upper(), + ) + else: + self._nests[idx] = np.random.normal(self._x0, self._sigma0) + def clip_nests(self, x): """ Clip the input array to the boundaries. """ - return np.clip(x, self._boundaries.lower(), self._boundaries.upper()) + if self._boundaries is not None: + x = np.clip(x, self._boundaries.lower(), self._boundaries.upper()) + return x def _suggested_population_size(self): """ diff --git a/pybop/optimisers/base_pints_optimiser.py b/pybop/optimisers/base_pints_optimiser.py index d3fa741f1..bfea92274 100644 --- a/pybop/optimisers/base_pints_optimiser.py +++ b/pybop/optimisers/base_pints_optimiser.py @@ -150,9 +150,9 @@ def _sanitise_inputs(self): raise ValueError( "Either all bounds or no bounds must be set for {self.pints_optimiser.__name__}." ) - self._boundaries = PintsRectangularBoundaries( - self.bounds["lower"], self.bounds["upper"] - ) + self._boundaries = PintsRectangularBoundaries( + self.bounds["lower"], self.bounds["upper"] + ) def name(self): """ From 011f5f882e20300a7e21781726827187ccbea7fa Mon Sep 17 00:00:00 2001 From: Brady Planden <55357039+BradyPlanden@users.noreply.github.com> Date: Fri, 7 Jun 2024 10:23:10 +0100 Subject: [PATCH 04/10] Apply suggestions from code review --- pybop/optimisers/_cuckoo.py | 10 +++++----- pybop/optimisers/pints_optimisers.py | 6 +++--- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/pybop/optimisers/_cuckoo.py b/pybop/optimisers/_cuckoo.py index 3a7b8522f..e61ecada4 100644 --- a/pybop/optimisers/_cuckoo.py +++ b/pybop/optimisers/_cuckoo.py @@ -5,18 +5,18 @@ class _CuckooSearch(PopulationBasedOptimiser): """ - Cuckoo Search (CS) optimization algorithm, inspired by the brood parasitism + Cuckoo Search (CS) optimisation algorithm, inspired by the brood parasitism of some cuckoo species. This algorithm was introduced by Yang and Deb in 2009. The algorithm uses a population of host nests (solutions), where each cuckoo (new solution) tries to replace a worse nest in the population. The quality - or fitness of the nests is determined by the objective function. A fraction + or fitness of the nests is determined by the cost function. A fraction of the worst nests is abandoned at each generation, and new ones are built randomly. The pseudo-code for the Cuckoo Search is as follows: - 1. Initialize population of n host nests + 1. Initialise population of n host nests 2. While (t < max_generations): a. Get a cuckoo randomly by Lévy flights b. Evaluate its quality/fitness F @@ -144,7 +144,7 @@ def levy_flight(self, alpha, size): def abandon_nests(self, idx): """ - Set the boundaries for the parameter space. + Updates the nests to abandon the worst performers and reinitialise. """ if self._boundaries is not None: self._nests[idx] = np.random.uniform( @@ -156,7 +156,7 @@ def abandon_nests(self, idx): def clip_nests(self, x): """ - Clip the input array to the boundaries. + Clip the input array to the boundaries if available. """ if self._boundaries is not None: x = np.clip(x, self._boundaries.lower(), self._boundaries.upper()) diff --git a/pybop/optimisers/pints_optimisers.py b/pybop/optimisers/pints_optimisers.py index 1e17d7be3..1a470092f 100644 --- a/pybop/optimisers/pints_optimisers.py +++ b/pybop/optimisers/pints_optimisers.py @@ -246,15 +246,15 @@ class CuckooSearch(BasePintsOptimiser): """ Adapter for the Cuckoo Search optimiser in PyBOP. - Cuckoo Search is a population-based optimization algorithm inspired by the brood parasitism of some cuckoo species. - It is designed to be simple, efficient, and robust, and is suitable for global optimization problems. + Cuckoo Search is a population-based optimisation algorithm inspired by the brood parasitism of some cuckoo species. + It is designed to be simple, efficient, and robust, and is suitable for global optimisation problems. Parameters ---------- **optimiser_kwargs : optional Valid PyBOP option keys and their values, for example: x0 : array_like - Initial + Initial parameter values. sigma0 : float Initial step size. bounds : dict From 5bffaf85ed3bc2b9e6e65a1b424c9ba1f4e9f268 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 7 Jun 2024 09:24:08 +0000 Subject: [PATCH 05/10] style: pre-commit fixes --- pybop/optimisers/_cuckoo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pybop/optimisers/_cuckoo.py b/pybop/optimisers/_cuckoo.py index e61ecada4..91cd92dbf 100644 --- a/pybop/optimisers/_cuckoo.py +++ b/pybop/optimisers/_cuckoo.py @@ -144,7 +144,7 @@ def levy_flight(self, alpha, size): def abandon_nests(self, idx): """ - Updates the nests to abandon the worst performers and reinitialise. + Updates the nests to abandon the worst performers and reinitialise. """ if self._boundaries is not None: self._nests[idx] = np.random.uniform( From 0974d54e152dd00be1a71273d46ea9a3313faee2 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Fri, 7 Jun 2024 13:24:51 +0100 Subject: [PATCH 06/10] refactor: updts _CuckooSearch to CuckooSearchImpl --- pybop/__init__.py | 2 +- pybop/optimisers/_cuckoo.py | 2 +- pybop/optimisers/pints_optimisers.py | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/pybop/__init__.py b/pybop/__init__.py index 2983da396..95f538cf9 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -102,7 +102,7 @@ # # Optimiser class # -from .optimisers._cuckoo import _CuckooSearch +from .optimisers._cuckoo import CuckooSearchImpl from .optimisers.base_optimiser import BaseOptimiser from .optimisers.base_pints_optimiser import BasePintsOptimiser from .optimisers.scipy_optimisers import ( diff --git a/pybop/optimisers/_cuckoo.py b/pybop/optimisers/_cuckoo.py index 91cd92dbf..5481f3af7 100644 --- a/pybop/optimisers/_cuckoo.py +++ b/pybop/optimisers/_cuckoo.py @@ -3,7 +3,7 @@ from scipy.special import gamma -class _CuckooSearch(PopulationBasedOptimiser): +class CuckooSearchImpl(PopulationBasedOptimiser): """ Cuckoo Search (CS) optimisation algorithm, inspired by the brood parasitism of some cuckoo species. This algorithm was introduced by Yang and Deb in 2009. diff --git a/pybop/optimisers/pints_optimisers.py b/pybop/optimisers/pints_optimisers.py index 1a470092f..3942a4734 100644 --- a/pybop/optimisers/pints_optimisers.py +++ b/pybop/optimisers/pints_optimisers.py @@ -7,7 +7,7 @@ from pints import IRPropMin as PintsIRPropMin from pints import NelderMead as PintsNelderMead -from pybop import BasePintsOptimiser, _CuckooSearch +from pybop import BasePintsOptimiser, CuckooSearchImpl class GradientDescent(BasePintsOptimiser): @@ -267,4 +267,4 @@ class CuckooSearch(BasePintsOptimiser): """ def __init__(self, cost, **optimiser_kwargs): - super().__init__(cost, _CuckooSearch, **optimiser_kwargs) + super().__init__(cost, CuckooSearchImpl, **optimiser_kwargs) From 991a59ce2ca863db25db52213b90903377eb708d Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Fri, 5 Jul 2024 09:59:48 +0100 Subject: [PATCH 07/10] tests: increase coverage --- pybop/optimisers/_cuckoo.py | 6 +++--- tests/unit/test_optimisation.py | 15 ++++++++++++++- 2 files changed, 17 insertions(+), 4 deletions(-) diff --git a/pybop/optimisers/_cuckoo.py b/pybop/optimisers/_cuckoo.py index 5481f3af7..9dcb9b735 100644 --- a/pybop/optimisers/_cuckoo.py +++ b/pybop/optimisers/_cuckoo.py @@ -64,7 +64,7 @@ def __init__(self, x0, sigma0=0.01, boundaries=None, pa=0.25): self._ready_for_tell = False # Initialise nests - if self._boundaries is not None: + if self._boundaries: self._nests = np.random.uniform( low=self._boundaries.lower(), high=self._boundaries.upper(), @@ -146,7 +146,7 @@ def abandon_nests(self, idx): """ Updates the nests to abandon the worst performers and reinitialise. """ - if self._boundaries is not None: + if self._boundaries: self._nests[idx] = np.random.uniform( low=self._boundaries.lower(), high=self._boundaries.upper(), @@ -158,7 +158,7 @@ def clip_nests(self, x): """ Clip the input array to the boundaries if available. """ - if self._boundaries is not None: + if self._boundaries: x = np.clip(x, self._boundaries.lower(), self._boundaries.upper()) return x diff --git a/tests/unit/test_optimisation.py b/tests/unit/test_optimisation.py index e16fab574..921ec3f20 100644 --- a/tests/unit/test_optimisation.py +++ b/tests/unit/test_optimisation.py @@ -176,6 +176,7 @@ def test_optimiser_kwargs(self, cost, optimiser): ): warnings.simplefilter("always") optim = optimiser(cost=cost, unrecognised=10) + assert not optim.pints_optimiser.running() else: # Check bounds in list format and update tol bounds = [ @@ -250,7 +251,6 @@ def test_optimiser_kwargs(self, cost, optimiser): # Check defaults assert optim.pints_optimiser.n_hyper_parameters() == 5 - assert not optim.pints_optimiser.running() assert optim.pints_optimiser.x_guessed() == optim.pints_optimiser._x0 with pytest.raises(Exception): optim.pints_optimiser.tell([0.1]) @@ -264,6 +264,19 @@ def test_optimiser_kwargs(self, cost, optimiser): assert optim.x0 == x0_new assert optim.x0 != x0 + @pytest.mark.unit + def test_cuckoo_no_bounds(self, dataset, model): + parameter = pybop.Parameter( + "Negative electrode active material volume fraction", + prior=pybop.Gaussian(0.6, 0.2), + ) + + cost_no_bounds = pybop.SumSquaredError( + pybop.FittingProblem(model, parameter, dataset) + ) + optim = pybop.CuckooSearch(cost=cost_no_bounds, max_iterations=1) + optim.run() + @pytest.mark.unit def test_scipy_minimize_with_jac(self, cost): # Check a method that uses gradient information From aa94180c426e0a89a1761ac2e9f7226837e51ff2 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Fri, 5 Jul 2024 11:20:04 +0100 Subject: [PATCH 08/10] fix: cuckoo boundaries==None clipping, add missing cuckoo to unit optimisation tests, updt. integration tests --- pybop/costs/_likelihoods.py | 2 +- pybop/optimisers/_cuckoo.py | 4 +++- tests/integration/test_spm_parameterisations.py | 2 +- tests/unit/test_optimisation.py | 14 ++++---------- 4 files changed, 9 insertions(+), 13 deletions(-) diff --git a/pybop/costs/_likelihoods.py b/pybop/costs/_likelihoods.py index 896d0c0dd..ed05b1345 100644 --- a/pybop/costs/_likelihoods.py +++ b/pybop/costs/_likelihoods.py @@ -265,7 +265,7 @@ class MAP(BaseLikelihood): """ - def __init__(self, problem, likelihood, sigma0=None, gradient_step=1e-2): + def __init__(self, problem, likelihood, sigma0=None, gradient_step=1e-3): super(MAP, self).__init__(problem) self.sigma0 = sigma0 self.gradient_step = gradient_step diff --git a/pybop/optimisers/_cuckoo.py b/pybop/optimisers/_cuckoo.py index 9dcb9b735..eda3f5c90 100644 --- a/pybop/optimisers/_cuckoo.py +++ b/pybop/optimisers/_cuckoo.py @@ -71,7 +71,9 @@ def __init__(self, x0, sigma0=0.01, boundaries=None, pa=0.25): size=(self._n, self._dim), ) else: - self._nests = np.random.normal(self._x0, self._sigma0) + self._nests = np.random.normal( + self._x0, self._sigma0, size=(self._n, self._dim) + ) self._fitness = np.full(self._n, np.inf) diff --git a/tests/integration/test_spm_parameterisations.py b/tests/integration/test_spm_parameterisations.py index c183d0222..85ddd6007 100644 --- a/tests/integration/test_spm_parameterisations.py +++ b/tests/integration/test_spm_parameterisations.py @@ -109,7 +109,7 @@ def test_spm_optimisers(self, optimiser, spm_costs): ) # Set sigma0 and create optimiser - sigma0 = 0.01 if isinstance(spm_costs, pybop.GaussianLogLikelihood) else 0.05 + sigma0 = 0.006 if isinstance(spm_costs, pybop.MAP) else None optim = optimiser(sigma0=sigma0, **common_args) # Set max unchanged iterations for BasePintsOptimisers diff --git a/tests/unit/test_optimisation.py b/tests/unit/test_optimisation.py index 921ec3f20..8444159b9 100644 --- a/tests/unit/test_optimisation.py +++ b/tests/unit/test_optimisation.py @@ -127,6 +127,7 @@ def test_no_optimisation_parameters(self, model, dataset): pybop.PSO, pybop.IRPropMin, pybop.NelderMead, + pybop.CuckooSearch, ], ) @pytest.mark.unit @@ -265,17 +266,10 @@ def test_optimiser_kwargs(self, cost, optimiser): assert optim.x0 != x0 @pytest.mark.unit - def test_cuckoo_no_bounds(self, dataset, model): - parameter = pybop.Parameter( - "Negative electrode active material volume fraction", - prior=pybop.Gaussian(0.6, 0.2), - ) - - cost_no_bounds = pybop.SumSquaredError( - pybop.FittingProblem(model, parameter, dataset) - ) - optim = pybop.CuckooSearch(cost=cost_no_bounds, max_iterations=1) + def test_cuckoo_no_bounds(self, dataset, cost, model): + optim = pybop.CuckooSearch(cost=cost, bounds=None, max_iterations=1) optim.run() + assert optim.pints_optimiser._boundaries is None @pytest.mark.unit def test_scipy_minimize_with_jac(self, cost): From f3038f9da0bca5b5fbabad62109d893c095cd7cc Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Fri, 5 Jul 2024 13:33:06 +0100 Subject: [PATCH 09/10] fix: align self._iterations with codebase --- pybop/optimisers/_cuckoo.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pybop/optimisers/_cuckoo.py b/pybop/optimisers/_cuckoo.py index eda3f5c90..0b5907a60 100644 --- a/pybop/optimisers/_cuckoo.py +++ b/pybop/optimisers/_cuckoo.py @@ -82,7 +82,7 @@ def __init__(self, x0, sigma0=0.01, boundaries=None, pa=0.25): self._f_best = np.inf # Set iteration count - self._iterations = 1 + self._iterations = 0 def ask(self): """ @@ -94,7 +94,7 @@ def ask(self): self._running = True # Generate new solutions (cuckoos) by Lévy flights - self.step_size = self._sigma0 / np.sqrt(self._iterations) + self.step_size = self._sigma0 / max(1, np.sqrt(self._iterations)) step = self.levy_flight(self.beta, self._dim) * self.step_size self.cuckoos = self._nests + step return self.clip_nests(self.cuckoos) From 86a297a7336539e35ad149bbe2370fef20976600 Mon Sep 17 00:00:00 2001 From: Brady Planden Date: Fri, 5 Jul 2024 14:03:34 +0100 Subject: [PATCH 10/10] add changelog entry --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index e5e93849d..26a6f942a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- [#319](https://github.com/pybop-team/PyBOP/pull/319/) - Adds `CuckooSearch` optimiser with corresponding tests. - [#379](https://github.com/pybop-team/PyBOP/pull/379) - Adds model.simulateS1 to weekly benchmarks. - [#174](https://github.com/pybop-team/PyBOP/issues/174) - Adds new logo and updates Readme for accessibility. - [#316](https://github.com/pybop-team/PyBOP/pull/316) - Adds Adam with weight decay (AdamW) optimiser, adds depreciation warning for pints.Adam implementation.