From 170845165a6dbd8b15fedfc668c2afb8217fec7d Mon Sep 17 00:00:00 2001 From: Prithvi-Rao Date: Fri, 18 Oct 2024 21:48:55 +0530 Subject: [PATCH 01/22] Added Simulated Annealing --- src/navigation.md | 1 + src/num_methods/simulated_annealing.md | 159 +++++++++++++++++++++++++ 2 files changed, 160 insertions(+) create mode 100644 src/num_methods/simulated_annealing.md diff --git a/src/navigation.md b/src/navigation.md index 29d6a94cb..6b7caef53 100644 --- a/src/navigation.md +++ b/src/navigation.md @@ -110,6 +110,7 @@ search: - [Binary Search](num_methods/binary_search.md) - [Ternary Search](num_methods/ternary_search.md) - [Newton's method for finding roots](num_methods/roots_newton.md) + - [Simulated Annealing](num_methods/simulated_annealing.md) - Integration - [Integration by Simpson's formula](num_methods/simpson-integration.md) - Geometry diff --git a/src/num_methods/simulated_annealing.md b/src/num_methods/simulated_annealing.md new file mode 100644 index 000000000..42879fe99 --- /dev/null +++ b/src/num_methods/simulated_annealing.md @@ -0,0 +1,159 @@ +--- +tags: + - Original +--- + +# Simulated Annealing + +**Simulated Annealing (SA)** is a randomized algorithm, which approximates the global optimum of a function. It's called a randomized algorithm, because it employs a certain amount of randomness in its search and thus its output can vary for the same input. + +## The Problem + +We are given a function $ E(s) $, which calculates the potential of the state $ s $. We are tasked with finding the state $ s_{best} $ at which $ E(s) $ is minimized. SA is suited for problems where the states are discrete and $ E(s) $ has multiple local minima. We'll take the example of the [Travelling Salesman Problem (TSP)](https://en.wikipedia.org/wiki/Travelling_salesman_problem). + +### Travelling Salesman Problem (TSP) + +You are given a set of nodes in 2 dimensional space. Each node is characterised by its $ x $ and $ y $ coordinates. Your task is to find the an ordering of the nodes, which will minimise the distance to be travelled when visiting these nodes in that order. + +### State + +State space is the collection of all possible values that can be taken by the independent variables. +A State is a unique point in the state space of the problem. In the case of TSP, all possible paths that we can take to visit all the nodes is the state space, and any single one of these paths can be considered as a State. + +### Neighbouring State + +It is a State in the State space which is close to the previous State. This usually means that we can obtain the neighbouring State from the original State using a simple transform. In the case of the Travelling salesman problem, a neighbouring state is obtained by randomly choosing 2 nodes, and swapping their positions in the current state. + +### E(s) + +$ E(s) $ is the function which needs to be minimised (or maximised). It maps every state to a real number. In the case of TSP, $ E(s) $ returns the distance of travelling one full circle in the order of nodes in the state. + +## Approach + +We start of with a random state $ s $. In every step, we choose a neighbouring state $ s_{next} $ of the current state $ s $. If $ E(s_{next}) < E(s) $, then we update $ s = s_{next} $. Otherwise, we use a probability acceptance function $ P(E(s),E(s_{next}),T) $ which decides whether we should move to $ s_{next} $ or stay at s. T here is the temperature, which is initially set to a high value and decays slowly with every step. The higher the temperature, the more likely it is to move to s_next. +At the same time we also keep a track of the best state $ s_{best} $ across all iterations. Proceed till convergence or time runs out. + + +## Probability Acceptance Function +$$ +P(E,E_{next},T) = + \begin{cases} + \text{True} &\quad\text{if } \exp{\frac{-(E_{next}-E)}{T}} \ge random(0,1) \\ + \text{False} &\quad\text{otherwise}\\ + \end{cases} + +$$ +This function takes in the current state, the next state and the Temperature , returning a boolean value, which tells our search whether it should move to $ s_{next} $ or stay at s. Note that for $E_{next} < E$ , this function will always return True. + +```cpp +bool P(double E,double E_next,double T){ + double e = 2.71828; + double prob = (double)rand()/RAND_MAX; // Generate a random number between 0 and 1 + if(pow(e,-(E_next-E)/T) > prob) return true; + else return false; +} +``` +## Code Template + +```cpp +class state{ + public: + state(){ + // Generate the initial state + } + state next(){ + state next; + next = s; + // Make changes to the state "next" and then return it + return next; + } + double E(){ + // implement the cost function here + }; +}; + +int main(){ + double T = 1000; // Initial temperature + double u = 0.99; // decay rate + state s = state(); + state best = s; + double E = s.E(); + double E_next; + double E_best = E; + while (T > 1){ + state next = s.next(); + E_next = next.E(); + if(P(E,E_next,T)){ + s = next; + if(E_next < E_best){ + best = s; + E_best = E_next; + } + } + E = E_next; + T *= u; + } + cout << E_best << "\n"; + return 0; +} +``` +## How to use: +Fill in the state class functions as appropriate. If you are trying to find a global maxima and not a minima, ensure that the $ E() $ function returns negative of the function you are maximising and finally print out $ -E_{best} $. Set the below parameters as per your need. + +### Parameters +- T : Temperature. Set it to a higher value if you want the search to run for a longer time +- u : Decay. Decides the rate of cooling. A slower cooling rate (larger value of u) usually gives better results. Ensure $u < 1$. + +The number of iterations the loop will run for is given by the expression +$$ +N = \lceil -\log_{u}{T} \rceil +$$ + +To see the effect of decay rate on solution results, run simulated annealing for decay rates 0.95 , 0.97 and 0.99 and see the difference. + +### Example State class for TSP +```cpp +class state{ + public: + vector> points; + + state(){ // Initial random order of points + points = {{0,0},{2,2},{0,2},{2,0},{0,1},{1,2},{2,1},{1,0}}; + } + state next(){ // picks 2 random indices and swaps them + state s_next; + s_next.points = points; + int a = ((rand()*points.size())/RAND_MAX); + int b = ((rand()*points.size())/RAND_MAX); + pair t = s_next.points[a]; + s_next.points[a] = s_next.points[b]; + s_next.points[b] = t; + return s_next; + } + + double euclidean(pair a, pair b){ // return euclidean distance between 2 points + return pow(pow((a.first-b.first),2)+pow((a.second-b.second),2),0.5); + } + double E(){ // calculates the round cost of travelling one full circle. + double dist = 0; + bool first = true; + int n = points.size(); + for(int i = 0;i < n; i++){ + dist += euclidean(points[i],points[(i+1)%n]); + } + return dist; + }; +}; +``` + +## Extra Modifications to the Algorithm: + +- Add a time based exit condition to the while loop to prevent TLE +- You can replace the e value in the Probability Acceptance function to any real number > 1. For a given $ E_{next} - E > 0 $, a higher e value reduces the chance of accepting that state and a smaller e value, increases it. + + +## Problems + +- https://usaco.org/index.php?page=viewproblem2&cpid=698 +- https://codeforces.com/contest/1556/problem/H +- https://atcoder.jp/contests/intro-heuristics/tasks/intro_heuristics_a \ No newline at end of file From 1eb4b14898ec98310b15c19d9aff20915f9f66f9 Mon Sep 17 00:00:00 2001 From: Prithvi-Rao Date: Sat, 19 Oct 2024 11:11:24 +0530 Subject: [PATCH 02/22] Formatting --- README.md | 1 + src/num_methods/simulated_annealing.md | 36 ++++++++++++-------------- 2 files changed, 18 insertions(+), 19 deletions(-) diff --git a/README.md b/README.md index 31ce219bc..87fce368e 100644 --- a/README.md +++ b/README.md @@ -29,6 +29,7 @@ Compiled pages are published at [https://cp-algorithms.com/](https://cp-algorith ### New articles +- (19 October 2024) [Simulated Annealing](https://cp-algorithms.com/num_methods/simulated_annealing.html) - (12 July 2024) [Manhattan distance](https://cp-algorithms.com/geometry/manhattan-distance.html) - (8 June 2024) [Knapsack Problem](https://cp-algorithms.com/dynamic_programming/knapsack.html) - (28 January 2024) [Introduction to Dynamic Programming](https://cp-algorithms.com/dynamic_programming/intro-to-dp.html) diff --git a/src/num_methods/simulated_annealing.md b/src/num_methods/simulated_annealing.md index 42879fe99..33dee6449 100644 --- a/src/num_methods/simulated_annealing.md +++ b/src/num_methods/simulated_annealing.md @@ -9,11 +9,11 @@ tags: ## The Problem -We are given a function $ E(s) $, which calculates the potential of the state $ s $. We are tasked with finding the state $ s_{best} $ at which $ E(s) $ is minimized. SA is suited for problems where the states are discrete and $ E(s) $ has multiple local minima. We'll take the example of the [Travelling Salesman Problem (TSP)](https://en.wikipedia.org/wiki/Travelling_salesman_problem). +We are given a function $E(s)$, which calculates the potential of the state $s$. We are tasked with finding the state $s_{best}$ at which $E(s)$ is minimized. SA is suited for problems where the states are discrete and $E(s)$ has multiple local minima. We'll take the example of the [Travelling Salesman Problem (TSP)](https://en.wikipedia.org/wiki/Travelling_salesman_problem). ### Travelling Salesman Problem (TSP) -You are given a set of nodes in 2 dimensional space. Each node is characterised by its $ x $ and $ y $ coordinates. Your task is to find the an ordering of the nodes, which will minimise the distance to be travelled when visiting these nodes in that order. +You are given a set of nodes in 2 dimensional space. Each node is characterised by its $x$ and $y$ coordinates. Your task is to find the an ordering of the nodes, which will minimise the distance to be travelled when visiting these nodes in that order. ### State @@ -26,12 +26,12 @@ It is a State in the State space which is close to the previous State. This usua ### E(s) -$ E(s) $ is the function which needs to be minimised (or maximised). It maps every state to a real number. In the case of TSP, $ E(s) $ returns the distance of travelling one full circle in the order of nodes in the state. +$E(s)$ is the function which needs to be minimised (or maximised). It maps every state to a real number. In the case of TSP, $E(s)$ returns the distance of travelling one full circle in the order of nodes in the state. ## Approach -We start of with a random state $ s $. In every step, we choose a neighbouring state $ s_{next} $ of the current state $ s $. If $ E(s_{next}) < E(s) $, then we update $ s = s_{next} $. Otherwise, we use a probability acceptance function $ P(E(s),E(s_{next}),T) $ which decides whether we should move to $ s_{next} $ or stay at s. T here is the temperature, which is initially set to a high value and decays slowly with every step. The higher the temperature, the more likely it is to move to s_next. -At the same time we also keep a track of the best state $ s_{best} $ across all iterations. Proceed till convergence or time runs out. +We start of with a random state $s$. In every step, we choose a neighbouring state $s_{next}$ of the current state $s$. If $E(s_{next}) < E(s)$, then we update $s = s_{next}$. Otherwise, we use a probability acceptance function $P(E(s),E(s_{next}),T)$ which decides whether we should move to $s_{next}$ or stay at s. T here is the temperature, which is initially set to a high value and decays slowly with every step. The higher the temperature, the more likely it is to move to s_next. +At the same time we also keep a track of the best state $s_{best}$ across all iterations. Proceed till convergence or time runs out. ## Probability Acceptance Function @@ -41,9 +41,8 @@ P(E,E_{next},T) = \text{True} &\quad\text{if } \exp{\frac{-(E_{next}-E)}{T}} \ge random(0,1) \\ \text{False} &\quad\text{otherwise}\\ \end{cases} - $$ -This function takes in the current state, the next state and the Temperature , returning a boolean value, which tells our search whether it should move to $ s_{next} $ or stay at s. Note that for $E_{next} < E$ , this function will always return True. +This function takes in the current state, the next state and the Temperature , returning a boolean value, which tells our search whether it should move to $s_{next}$ or stay at s. Note that for $E_{next} < E$ , this function will always return True. ```cpp bool P(double E,double E_next,double T){ @@ -72,13 +71,12 @@ class state{ }; }; -int main(){ - double T = 1000; // Initial temperature - double u = 0.99; // decay rate +pair simAnneal(){ state s = state(); state best = s; - double E = s.E(); - double E_next; + double T = 1000; // Initial temperature + double u = 0.99; // decay rate + double E = s.E(),E_next; double E_best = E; while (T > 1){ state next = s.next(); @@ -93,12 +91,12 @@ int main(){ E = E_next; T *= u; } - cout << E_best << "\n"; - return 0; + return {E_best,best}; } + ``` ## How to use: -Fill in the state class functions as appropriate. If you are trying to find a global maxima and not a minima, ensure that the $ E() $ function returns negative of the function you are maximising and finally print out $ -E_{best} $. Set the below parameters as per your need. +Fill in the state class functions as appropriate. If you are trying to find a global maxima and not a minima, ensure that the $E()$ function returns negative of the function you are maximising and finally print out $-E_{best}$. Set the below parameters as per your need. ### Parameters - T : Temperature. Set it to a higher value if you want the search to run for a longer time @@ -149,11 +147,11 @@ class state{ ## Extra Modifications to the Algorithm: - Add a time based exit condition to the while loop to prevent TLE -- You can replace the e value in the Probability Acceptance function to any real number > 1. For a given $ E_{next} - E > 0 $, a higher e value reduces the chance of accepting that state and a smaller e value, increases it. +- You can replace the e value in the Probability Acceptance function to any real number > 1. For a given $E_{next} - E > 0$, a higher e value reduces the chance of accepting that state and a smaller e value, increases it. ## Problems -- https://usaco.org/index.php?page=viewproblem2&cpid=698 -- https://codeforces.com/contest/1556/problem/H -- https://atcoder.jp/contests/intro-heuristics/tasks/intro_heuristics_a \ No newline at end of file +- [USACO Jan 2017 - Subsequence Reversal](https://usaco.org/index.php?page=viewproblem2&cpid=698) +- [Deltix Summer 2021 - DIY Tree](https://codeforces.com/contest/1556/problem/H) +- [AtCoder Contest Scheduling](https://atcoder.jp/contests/intro-heuristics/tasks/intro_heuristics_a) \ No newline at end of file From ebd3cdf18838fc6ca551bf646a647a4e3b885e28 Mon Sep 17 00:00:00 2001 From: Prithvi-Rao Date: Sat, 19 Oct 2024 11:18:54 +0530 Subject: [PATCH 03/22] fixed compile issue --- src/num_methods/simulated_annealing.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/num_methods/simulated_annealing.md b/src/num_methods/simulated_annealing.md index 33dee6449..6fca91cd8 100644 --- a/src/num_methods/simulated_annealing.md +++ b/src/num_methods/simulated_annealing.md @@ -116,7 +116,7 @@ class state{ vector> points; state(){ // Initial random order of points - points = {{0,0},{2,2},{0,2},{2,0},{0,1},{1,2},{2,1},{1,0}}; + points = {}; // Fill in some points to start with, or generate them randomly } state next(){ // picks 2 random indices and swaps them state s_next; From 9c81f7917ad11673fbab4be323b08043de8ea45a Mon Sep 17 00:00:00 2001 From: Prithvi-Rao Date: Tue, 22 Oct 2024 00:37:59 +0530 Subject: [PATCH 04/22] Elaborated temp --- src/num_methods/simulated_annealing.md | 47 ++++++++++++++++---------- 1 file changed, 30 insertions(+), 17 deletions(-) diff --git a/src/num_methods/simulated_annealing.md b/src/num_methods/simulated_annealing.md index 6fca91cd8..1af1d49cd 100644 --- a/src/num_methods/simulated_annealing.md +++ b/src/num_methods/simulated_annealing.md @@ -9,7 +9,7 @@ tags: ## The Problem -We are given a function $E(s)$, which calculates the potential of the state $s$. We are tasked with finding the state $s_{best}$ at which $E(s)$ is minimized. SA is suited for problems where the states are discrete and $E(s)$ has multiple local minima. We'll take the example of the [Travelling Salesman Problem (TSP)](https://en.wikipedia.org/wiki/Travelling_salesman_problem). +We are given a function $E(s)$, which calculates the potential of the state $s$. We are tasked with finding the state $s_{best}$ at which $E(s)$ is minimized. **SA** is suited for problems where the states are discrete and $E(s)$ has multiple local minima. We'll take the example of the [Travelling Salesman Problem (TSP)](https://en.wikipedia.org/wiki/Travelling_salesman_problem). ### Travelling Salesman Problem (TSP) @@ -20,29 +20,43 @@ You are given a set of nodes in 2 dimensional space. Each node is characterised State space is the collection of all possible values that can be taken by the independent variables. A State is a unique point in the state space of the problem. In the case of TSP, all possible paths that we can take to visit all the nodes is the state space, and any single one of these paths can be considered as a State. -### Neighbouring State +### Neighbouring state -It is a State in the State space which is close to the previous State. This usually means that we can obtain the neighbouring State from the original State using a simple transform. In the case of the Travelling salesman problem, a neighbouring state is obtained by randomly choosing 2 nodes, and swapping their positions in the current state. +It is a state in the state space which is close to the previous state. This usually means that we can obtain the neighbouring state from the original state using a simple transform. In the case of the Travelling Salesman Problem, a neighbouring state is obtained by randomly choosing 2 nodes, and swapping their positions in the current state. -### E(s) +### The Energy Function E(s) $E(s)$ is the function which needs to be minimised (or maximised). It maps every state to a real number. In the case of TSP, $E(s)$ returns the distance of travelling one full circle in the order of nodes in the state. -## Approach +## The Approach -We start of with a random state $s$. In every step, we choose a neighbouring state $s_{next}$ of the current state $s$. If $E(s_{next}) < E(s)$, then we update $s = s_{next}$. Otherwise, we use a probability acceptance function $P(E(s),E(s_{next}),T)$ which decides whether we should move to $s_{next}$ or stay at s. T here is the temperature, which is initially set to a high value and decays slowly with every step. The higher the temperature, the more likely it is to move to s_next. -At the same time we also keep a track of the best state $s_{best}$ across all iterations. Proceed till convergence or time runs out. +We start of with a random state $s$. In every step, we choose a neighbouring state $s_{next}$ of the current state $s$. If $E(s_{next}) < E(s)$, then we update $s = s_{next}$. Otherwise, we use a probability acceptance function $P(E(s),E(s_{next}),T)$ which decides whether we should move to $s_{next}$ or stay at $s$. T here is the temperature, which is initially set to a high value and decays slowly with every step. The higher the temperature, the more likely it is to move to $s_{next}$. +At the same time we also keep a track of the best state $s_{best}$ across all iterations. Proceeding till convergence or time runs out. +### How does this work? + +This algorithm is called simulated annealing because we are simulating the process of annealing, wherein a material is heated up and allowed to cool, in order to allow the atoms inside to rearrange themselves in an arrangement with minimal internal energy, which in turn causes the material to have different properties. The state is the arrangement of atoms and the internal energy is the function being minimised. We can think of the original state of the atoms, as a local minima for its internal energy. To make the material rearrange its atoms, we need to motivate it to go across a region where its internal energy is not minimised in order to reach the global minima. This motivation is given by heating the material to a higher temperature. + +Simulated annealing, literally simulates this process. We start off with some random state (material) and set a high temperature (heat it up). Now, the algorithm is ready to accept states which have a higher energy than the current state, as it is motivated by the high value of $T$. This prevents the algorithm from getting stuck inside local minimas and move towards the global minima. As time progresses, the algorithm cools down and refuses the states with higher energy and moves into the closest minima it has found. + + +
+ +
+A visual representation of simulated annealing, searching for the maxima of this function with multiple local maxima.. +
+This gif by [Kingpin13](https://commons.wikimedia.org/wiki/User:Kingpin13) is distributed under CC0 1.0 license. +
## Probability Acceptance Function -$$ -P(E,E_{next},T) = + +$P(E,E_{next},T) = \begin{cases} \text{True} &\quad\text{if } \exp{\frac{-(E_{next}-E)}{T}} \ge random(0,1) \\ \text{False} &\quad\text{otherwise}\\ - \end{cases} -$$ -This function takes in the current state, the next state and the Temperature , returning a boolean value, which tells our search whether it should move to $s_{next}$ or stay at s. Note that for $E_{next} < E$ , this function will always return True. + \end{cases}$ + +This function takes in the current state, the next state and the Temperature , returning a boolean value, which tells our search whether it should move to $s_{next}$ or stay at $s$. Note that for $E_{next} < E$ , this function will always return True. ```cpp bool P(double E,double E_next,double T){ @@ -99,13 +113,12 @@ pair simAnneal(){ Fill in the state class functions as appropriate. If you are trying to find a global maxima and not a minima, ensure that the $E()$ function returns negative of the function you are maximising and finally print out $-E_{best}$. Set the below parameters as per your need. ### Parameters -- T : Temperature. Set it to a higher value if you want the search to run for a longer time -- u : Decay. Decides the rate of cooling. A slower cooling rate (larger value of u) usually gives better results. Ensure $u < 1$. +- $T$ : Temperature. Set it to a higher value if you want the search to run for a longer time +- $u$ : Decay. Decides the rate of cooling. A slower cooling rate (larger value of u) usually gives better results. Ensure $u < 1$. The number of iterations the loop will run for is given by the expression -$$ -N = \lceil -\log_{u}{T} \rceil -$$ + +$N = \lceil -\log_{u}{T} \rceil$ To see the effect of decay rate on solution results, run simulated annealing for decay rates 0.95 , 0.97 and 0.99 and see the difference. From 1962436bf3f4c970a7d6c265263186afdf721e12 Mon Sep 17 00:00:00 2001 From: Prithvi-Rao Date: Tue, 29 Oct 2024 09:53:44 +0530 Subject: [PATCH 05/22] updated code --- src/num_methods/simulated_annealing.md | 96 +++++++++++++++++--------- 1 file changed, 62 insertions(+), 34 deletions(-) diff --git a/src/num_methods/simulated_annealing.md b/src/num_methods/simulated_annealing.md index 1af1d49cd..0e650625d 100644 --- a/src/num_methods/simulated_annealing.md +++ b/src/num_methods/simulated_annealing.md @@ -7,7 +7,7 @@ tags: **Simulated Annealing (SA)** is a randomized algorithm, which approximates the global optimum of a function. It's called a randomized algorithm, because it employs a certain amount of randomness in its search and thus its output can vary for the same input. -## The Problem +## The problem We are given a function $E(s)$, which calculates the potential of the state $s$. We are tasked with finding the state $s_{best}$ at which $E(s)$ is minimized. **SA** is suited for problems where the states are discrete and $E(s)$ has multiple local minima. We'll take the example of the [Travelling Salesman Problem (TSP)](https://en.wikipedia.org/wiki/Travelling_salesman_problem). @@ -18,17 +18,17 @@ You are given a set of nodes in 2 dimensional space. Each node is characterised ### State State space is the collection of all possible values that can be taken by the independent variables. -A State is a unique point in the state space of the problem. In the case of TSP, all possible paths that we can take to visit all the nodes is the state space, and any single one of these paths can be considered as a State. +A state is a unique point in the state space of the problem. In the case of TSP, all possible paths that we can take to visit all the nodes is the state space, and any single one of these paths can be considered as a state. ### Neighbouring state It is a state in the state space which is close to the previous state. This usually means that we can obtain the neighbouring state from the original state using a simple transform. In the case of the Travelling Salesman Problem, a neighbouring state is obtained by randomly choosing 2 nodes, and swapping their positions in the current state. -### The Energy Function E(s) +### The energy function E(s) $E(s)$ is the function which needs to be minimised (or maximised). It maps every state to a real number. In the case of TSP, $E(s)$ returns the distance of travelling one full circle in the order of nodes in the state. -## The Approach +## The approach We start of with a random state $s$. In every step, we choose a neighbouring state $s_{next}$ of the current state $s$. If $E(s_{next}) < E(s)$, then we update $s = s_{next}$. Otherwise, we use a probability acceptance function $P(E(s),E(s_{next}),T)$ which decides whether we should move to $s_{next}$ or stay at $s$. T here is the temperature, which is initially set to a high value and decays slowly with every step. The higher the temperature, the more likely it is to move to $s_{next}$. At the same time we also keep a track of the best state $s_{best}$ across all iterations. Proceeding till convergence or time runs out. @@ -48,7 +48,7 @@ Simulated annealing, literally simulates this process. We start off with some ra This gif by [Kingpin13](https://commons.wikimedia.org/wiki/User:Kingpin13) is distributed under CC0 1.0 license. -## Probability Acceptance Function +## Probability Acceptance Function(PAF) $P(E,E_{next},T) = \begin{cases} @@ -60,9 +60,8 @@ This function takes in the current state, the next state and the Temperature , r ```cpp bool P(double E,double E_next,double T){ - double e = 2.71828; double prob = (double)rand()/RAND_MAX; // Generate a random number between 0 and 1 - if(pow(e,-(E_next-E)/T) > prob) return true; + if(exp(-(E_next-E)/T) > prob) return true; else return false; } ``` @@ -74,23 +73,28 @@ class state{ state(){ // Generate the initial state } + state(state& s){ + // Implement the copy constructor + } state next(){ - state next; - next = s; - // Make changes to the state "next" and then return it - return next; + state s_next = state(*this); + + // Modify s_next to the neighbouring state + + return s_next; } double E(){ - // implement the cost function here + // Implement the cost function here }; }; -pair simAnneal(){ +pair simAnneal(){ state s = state(); - state best = s; + state best = state(s); double T = 1000; // Initial temperature double u = 0.99; // decay rate - double E = s.E(),E_next; + double E = s.E(); + double E_next; double E_best = E; while (T > 1){ state next = s.next(); @@ -113,39 +117,42 @@ pair simAnneal(){ Fill in the state class functions as appropriate. If you are trying to find a global maxima and not a minima, ensure that the $E()$ function returns negative of the function you are maximising and finally print out $-E_{best}$. Set the below parameters as per your need. ### Parameters -- $T$ : Temperature. Set it to a higher value if you want the search to run for a longer time -- $u$ : Decay. Decides the rate of cooling. A slower cooling rate (larger value of u) usually gives better results. Ensure $u < 1$. +- $T$ : Temperature. Set it to a higher value if you want the search to run for a longer time. +- $u$ : Decay. Decides the rate of cooling. A slower cooling rate (larger value of u) usually gives better results, at the cost of running for a longer time. Ensure $u < 1$. The number of iterations the loop will run for is given by the expression $N = \lceil -\log_{u}{T} \rceil$ -To see the effect of decay rate on solution results, run simulated annealing for decay rates 0.95 , 0.97 and 0.99 and see the difference. +Tips for choosing $T$ and $u$ : If there are many local minimas and a wide state space, set $u = 0.999$, for a slow cooling rate, which will allow the algorithm to explore more possibilities. On the other hand, if the state space is narrower, $u = 0.99$ should suffice. If you are not sure, play it safe by setting $u = 0.998$ or higher. Calculate the time complexity of a single iteration of the algorithm, and use this to approximate a value of $N$ which will prevent TLE, then use the below formula to obtain $T$. + +$T = u^{-N}$ -### Example State class for TSP +### Example implementation for TSP ```cpp + class state{ public: vector> points; - state(){ // Initial random order of points - points = {}; // Fill in some points to start with, or generate them randomly + state(){ + points = { {0,0},{2,2},{0,2},{2,0},{0,1},{1,2},{2,1},{1,0} }; + } + state(state& s){ + points = s.points; } - state next(){ // picks 2 random indices and swaps them - state s_next; - s_next.points = points; - int a = ((rand()*points.size())/RAND_MAX); - int b = ((rand()*points.size())/RAND_MAX); - pair t = s_next.points[a]; - s_next.points[a] = s_next.points[b]; - s_next.points[b] = t; + state next(){ + state s_next = state(*this); + int a = (rand()%points.size()); + int b = (rand()%points.size()); + s_next.points[a].swap(s_next.points[b]); return s_next; } - double euclidean(pair a, pair b){ // return euclidean distance between 2 points + double euclidean(pair a, pair b){ return pow(pow((a.first-b.first),2)+pow((a.second-b.second),2),0.5); } - double E(){ // calculates the round cost of travelling one full circle. + double E(){ double dist = 0; bool first = true; int n = points.size(); @@ -155,13 +162,34 @@ class state{ return dist; }; }; + + +int main(){ + pair res; + res = simAnneal(); + double E_best = res.first; + state best = res.second; + cout << "Lenght of shortest path found : " << E_best << "\n"; + cout << "Order of points in shortest path : \n"; + for(auto x: best.points){ + cout << x.first << " " << x.second << "\n"; + } +} ``` -## Extra Modifications to the Algorithm: +## Further modifications to the algorithm: - Add a time based exit condition to the while loop to prevent TLE -- You can replace the e value in the Probability Acceptance function to any real number > 1. For a given $E_{next} - E > 0$, a higher e value reduces the chance of accepting that state and a smaller e value, increases it. - +- The Probability acceptance function given above, prefers accepting states which are lower in energy because of the $E_{next} - E$ factor in the numerator of the exponent. You can simply remove this factor, to make the PAF independent of the difference in energies. +- The effect of the difference in energies, $E_{next} - E$ on the PAF can be increased/decreased by increasing/decreasing the base of the exponent as shown below: +```cpp +bool P(double E,double E_next,double T){ + double e = 2; // set e to any real number greater than 1 + double prob = (double)rand()/RAND_MAX; // Generate a random number between 0 and 1 + if(pow(e,-(E_next-E)/T) > prob) return true; + else return false; +} +``` ## Problems From 7b46a269250c13aa7ae3c63a2c75d1e984ffa871 Mon Sep 17 00:00:00 2001 From: Prithvi-Rao Date: Tue, 29 Oct 2024 10:38:33 +0530 Subject: [PATCH 06/22] fixed braces error --- src/num_methods/simulated_annealing.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/num_methods/simulated_annealing.md b/src/num_methods/simulated_annealing.md index 0e650625d..5c06de057 100644 --- a/src/num_methods/simulated_annealing.md +++ b/src/num_methods/simulated_annealing.md @@ -136,7 +136,7 @@ class state{ vector> points; state(){ - points = { {0,0},{2,2},{0,2},{2,0},{0,1},{1,2},{2,1},{1,0} }; + {% raw %}points = {{0,0},{2,2},{0,2},{2,0},{0,1},{1,2},{2,1},{1,0}};{% endraw %} } state(state& s){ points = s.points; @@ -181,7 +181,7 @@ int main(){ - Add a time based exit condition to the while loop to prevent TLE - The Probability acceptance function given above, prefers accepting states which are lower in energy because of the $E_{next} - E$ factor in the numerator of the exponent. You can simply remove this factor, to make the PAF independent of the difference in energies. -- The effect of the difference in energies, $E_{next} - E$ on the PAF can be increased/decreased by increasing/decreasing the base of the exponent as shown below: +- The effect of the difference in energies, $E_{next} - E$, on the PAF can be increased/decreased by increasing/decreasing the base of the exponent as shown below: ```cpp bool P(double E,double E_next,double T){ double e = 2; // set e to any real number greater than 1 From bab84cda941be4061f3d139252b8cdc370ed6002 Mon Sep 17 00:00:00 2001 From: Prithvi-Rao Date: Tue, 5 Nov 2024 09:16:10 +0530 Subject: [PATCH 07/22] fixes --- src/num_methods/simulated_annealing.md | 150 +++++++++++++------------ 1 file changed, 78 insertions(+), 72 deletions(-) diff --git a/src/num_methods/simulated_annealing.md b/src/num_methods/simulated_annealing.md index 5c06de057..abd035950 100644 --- a/src/num_methods/simulated_annealing.md +++ b/src/num_methods/simulated_annealing.md @@ -9,35 +9,33 @@ tags: ## The problem -We are given a function $E(s)$, which calculates the potential of the state $s$. We are tasked with finding the state $s_{best}$ at which $E(s)$ is minimized. **SA** is suited for problems where the states are discrete and $E(s)$ has multiple local minima. We'll take the example of the [Travelling Salesman Problem (TSP)](https://en.wikipedia.org/wiki/Travelling_salesman_problem). +We are given a function $E(s)$, which calculates the energy of the state $s$. We are tasked with finding the state $s_{best}$ at which $E(s)$ is minimized. **SA** is suited for problems where the states are discrete and $E(s)$ has multiple local minima. We'll take the example of the [Travelling Salesman Problem (TSP)](https://en.wikipedia.org/wiki/Travelling_salesman_problem). ### Travelling Salesman Problem (TSP) -You are given a set of nodes in 2 dimensional space. Each node is characterised by its $x$ and $y$ coordinates. Your task is to find the an ordering of the nodes, which will minimise the distance to be travelled when visiting these nodes in that order. +You are given a set of nodes in 2 dimensional space. Each node is characterised by its $x$ and $y$ coordinates. Your task is to find an ordering of the nodes, which will minimise the distance to be travelled when visiting these nodes in that order. -### State - -State space is the collection of all possible values that can be taken by the independent variables. -A state is a unique point in the state space of the problem. In the case of TSP, all possible paths that we can take to visit all the nodes is the state space, and any single one of these paths can be considered as a state. - -### Neighbouring state +## Motivation +Annealing is a metallurgical process , wherein a material is heated up and allowed to cool, in order to allow the atoms inside to rearrange themselves in an arrangement with minimal internal energy, which in turn causes the material to have different properties. The state is the arrangement of atoms and the internal energy is the function being minimised. We can think of the original state of the atoms, as a local minima for its internal energy. To make the material rearrange its atoms, we need to motivate it to go across a region where its internal energy is not minimised in order to reach the global minima. This motivation is given by heating the material to a higher temperature. -It is a state in the state space which is close to the previous state. This usually means that we can obtain the neighbouring state from the original state using a simple transform. In the case of the Travelling Salesman Problem, a neighbouring state is obtained by randomly choosing 2 nodes, and swapping their positions in the current state. +Simulated annealing, literally, simulates this process. We start off with some random state (material) and set a high temperature (heat it up). Now, the algorithm is ready to accept states which have a higher energy than the current state, as it is motivated by the high temperature. This prevents the algorithm from getting stuck inside local minimas and move towards the global minima. As time progresses, the algorithm cools down and refuses the states with higher energy and moves into the closest minima it has found. ### The energy function E(s) $E(s)$ is the function which needs to be minimised (or maximised). It maps every state to a real number. In the case of TSP, $E(s)$ returns the distance of travelling one full circle in the order of nodes in the state. -## The approach +### State + +The state space is the domain of the energy function, E(s), and a state is any element which belongs to the state space. In the case of TSP, all possible paths that we can take to visit all the nodes is the state space, and any single one of these paths can be considered as a state. -We start of with a random state $s$. In every step, we choose a neighbouring state $s_{next}$ of the current state $s$. If $E(s_{next}) < E(s)$, then we update $s = s_{next}$. Otherwise, we use a probability acceptance function $P(E(s),E(s_{next}),T)$ which decides whether we should move to $s_{next}$ or stay at $s$. T here is the temperature, which is initially set to a high value and decays slowly with every step. The higher the temperature, the more likely it is to move to $s_{next}$. -At the same time we also keep a track of the best state $s_{best}$ across all iterations. Proceeding till convergence or time runs out. +### Neighbouring state -### How does this work? +It is a state in the state space which is close to the previous state. This usually means that we can obtain the neighbouring state from the original state using a simple transform. In the case of the Travelling Salesman Problem, a neighbouring state is obtained by randomly choosing 2 nodes, and swapping their positions in the current state. -This algorithm is called simulated annealing because we are simulating the process of annealing, wherein a material is heated up and allowed to cool, in order to allow the atoms inside to rearrange themselves in an arrangement with minimal internal energy, which in turn causes the material to have different properties. The state is the arrangement of atoms and the internal energy is the function being minimised. We can think of the original state of the atoms, as a local minima for its internal energy. To make the material rearrange its atoms, we need to motivate it to go across a region where its internal energy is not minimised in order to reach the global minima. This motivation is given by heating the material to a higher temperature. +## Algorithm -Simulated annealing, literally simulates this process. We start off with some random state (material) and set a high temperature (heat it up). Now, the algorithm is ready to accept states which have a higher energy than the current state, as it is motivated by the high value of $T$. This prevents the algorithm from getting stuck inside local minimas and move towards the global minima. As time progresses, the algorithm cools down and refuses the states with higher energy and moves into the closest minima it has found. +We start with a random state $s$. In every step, we choose a neighbouring state $s_{next}$ of the current state $s$. If $E(s_{next}) < E(s)$, then we update $s = s_{next}$. Otherwise, we use a probability acceptance function $P(E(s),E(s_{next}),T)$ which decides whether we should move to $s_{next}$ or stay at $s$. T here is the temperature, which is initially set to a high value and decays slowly with every step. The higher the temperature, the more likely it is to move to $s_{next}$. +At the same time we also keep a track of the best state $s_{best}$ across all iterations. Proceeding till convergence or time runs out.
@@ -45,63 +43,66 @@ Simulated annealing, literally simulates this process. We start off with some ra
A visual representation of simulated annealing, searching for the maxima of this function with multiple local maxima..
-This gif by [Kingpin13](https://commons.wikimedia.org/wiki/User:Kingpin13) is distributed under CC0 1.0 license. -
+This animation by [Kingpin13](https://commons.wikimedia.org/wiki/User:Kingpin13) is distributed under CC0 1.0 license. + +### Temperature($T$) and decay($u$) + +The temperature of the system quantifies the willingness of the algorithm to accept a state with a higher energy. The decay is a constant which quantifies the "cooling rate" of the algorithm. A slow cooling rate (larger $u$) is known to give better results. ## Probability Acceptance Function(PAF) $P(E,E_{next},T) = \begin{cases} - \text{True} &\quad\text{if } \exp{\frac{-(E_{next}-E)}{T}} \ge random(0,1) \\ + \text{True} &\quad\text{if } \mathcal{U}_{[0,1]} \le \exp(-\frac{E_{next}-E}{T}) \\ \text{False} &\quad\text{otherwise}\\ \end{cases}$ -This function takes in the current state, the next state and the Temperature , returning a boolean value, which tells our search whether it should move to $s_{next}$ or stay at $s$. Note that for $E_{next} < E$ , this function will always return True. +Here, $\mathcal{U}_{[0,1]}$ is a continuous uniform random value on $[0,1]$. This function takes in the current state, the next state and the temperature, returning a boolean value, which tells our search whether it should move to $s_{next}$ or stay at $s$. Note that for $E_{next} < E$ , this function will always return True, otherwise it can still make the move with probability $\exp(-\frac{E_{next}-E}{T})$, which corresponds to the [Gibbs measure](https://en.wikipedia.org/wiki/Gibbs_measure). ```cpp -bool P(double E,double E_next,double T){ - double prob = (double)rand()/RAND_MAX; // Generate a random number between 0 and 1 - if(exp(-(E_next-E)/T) > prob) return true; - else return false; +bool P(double E,double E_next,double T,mt19937 rng){ + double prob = exp(-(E_next-E)/T); + if(prob > 1) return true; + else{ + bernoulli_distribution d(prob); + return d(rng); + } } ``` ## Code Template ```cpp -class state{ +class state { public: - state(){ + state() { // Generate the initial state } - state(state& s){ - // Implement the copy constructor - } - state next(){ - state s_next = state(*this); - - // Modify s_next to the neighbouring state - + state next() { + state s_next; + // Modify s_next to a random neighboring state return s_next; } - double E(){ - // Implement the cost function here + double E() { + // Implement the energy function here }; }; -pair simAnneal(){ + +pair simAnneal() { state s = state(); - state best = state(s); - double T = 1000; // Initial temperature - double u = 0.99; // decay rate + state best = s; + double T = 10000; // Initial temperature + double u = 0.995; // decay rate double E = s.E(); double E_next; double E_best = E; - while (T > 1){ + mt19937 rng(chrono::steady_clock::now().time_since_epoch().count()); + while (T > 1) { state next = s.next(); E_next = next.E(); - if(P(E,E_next,T)){ + if (P(E, E_next, T, rng)) { s = next; - if(E_next < E_best){ + if (E_next < E_best) { best = s; E_best = E_next; } @@ -109,15 +110,15 @@ pair simAnneal(){ E = E_next; T *= u; } - return {E_best,best}; + return {E_best, best}; } ``` ## How to use: -Fill in the state class functions as appropriate. If you are trying to find a global maxima and not a minima, ensure that the $E()$ function returns negative of the function you are maximising and finally print out $-E_{best}$. Set the below parameters as per your need. +Fill in the state class functions as appropriate. If you are trying to find a global maxima and not a minima, ensure that the $E()$ function returns negative of the function you are maximizing and print $-E_{best}$ in the end. Set the below parameters as per your need. ### Parameters -- $T$ : Temperature. Set it to a higher value if you want the search to run for a longer time. +- $T$ : Initial temperature. Set it to a higher value if you want the search to run for a longer time. - $u$ : Decay. Decides the rate of cooling. A slower cooling rate (larger value of u) usually gives better results, at the cost of running for a longer time. Ensure $u < 1$. The number of iterations the loop will run for is given by the expression @@ -131,47 +132,47 @@ $T = u^{-N}$ ### Example implementation for TSP ```cpp -class state{ +class state { public: - vector> points; - - state(){ - {% raw %}points = {{0,0},{2,2},{0,2},{2,0},{0,1},{1,2},{2,1},{1,0}};{% endraw %} - } - state(state& s){ - points = s.points; + vector> points; + std::mt19937 mt{ static_cast( + std::chrono::steady_clock::now().time_since_epoch().count() + ) }; + state() { + points = {%raw%} {{0,0},{2,2},{0,2},{2,0},{0,1},{1,2},{2,1},{1,0}} {%endraw%}; } - state next(){ - state s_next = state(*this); - int a = (rand()%points.size()); - int b = (rand()%points.size()); + state next() { + state s_next; + s_next.points = points; + uniform_int_distribution<> choose(0, points.size()-1); + int a = choose(mt); + int b = choose(mt); s_next.points[a].swap(s_next.points[b]); return s_next; } - double euclidean(pair a, pair b){ - return pow(pow((a.first-b.first),2)+pow((a.second-b.second),2),0.5); + double euclidean(pair a, pair b) { + return hypot(a.first - b.first, a.second - b.second); } - double E(){ + + double E() { double dist = 0; bool first = true; int n = points.size(); - for(int i = 0;i < n; i++){ - dist += euclidean(points[i],points[(i+1)%n]); - } + for (int i = 0;i < n; i++) + dist += euclidean(points[i], points[(i+1)%n]); return dist; }; }; - -int main(){ - pair res; +int main() { + pair res; res = simAnneal(); double E_best = res.first; state best = res.second; cout << "Lenght of shortest path found : " << E_best << "\n"; cout << "Order of points in shortest path : \n"; - for(auto x: best.points){ + for(auto x: best.points) { cout << x.first << " " << x.second << "\n"; } } @@ -180,14 +181,19 @@ int main(){ ## Further modifications to the algorithm: - Add a time based exit condition to the while loop to prevent TLE +- The decay implemented above is an exponential decay. You can always replace this with a decay function as per your needs. - The Probability acceptance function given above, prefers accepting states which are lower in energy because of the $E_{next} - E$ factor in the numerator of the exponent. You can simply remove this factor, to make the PAF independent of the difference in energies. - The effect of the difference in energies, $E_{next} - E$, on the PAF can be increased/decreased by increasing/decreasing the base of the exponent as shown below: ```cpp -bool P(double E,double E_next,double T){ - double e = 2; // set e to any real number greater than 1 - double prob = (double)rand()/RAND_MAX; // Generate a random number between 0 and 1 - if(pow(e,-(E_next-E)/T) > prob) return true; - else return false; +bool P(double E, double E_next, double T, mt19937 rng) { + e = 2 // set e to any real number greater than 1 + double prob = exp(-(E_next-E)/T); + if (prob > 1) + return true; + else { + bernoulli_distribution d(prob); + return d(rng); + } } ``` From 7c1254ba49761d3eaedb7e7387c3489a1e5647de Mon Sep 17 00:00:00 2001 From: Oleksandr Kulkov Date: Mon, 24 Mar 2025 21:26:18 +0100 Subject: [PATCH 08/22] Empty commit for GitHub Actions --- src/num_methods/simulated_annealing.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/num_methods/simulated_annealing.md b/src/num_methods/simulated_annealing.md index abd035950..a5f6b466a 100644 --- a/src/num_methods/simulated_annealing.md +++ b/src/num_methods/simulated_annealing.md @@ -201,4 +201,4 @@ bool P(double E, double E_next, double T, mt19937 rng) { - [USACO Jan 2017 - Subsequence Reversal](https://usaco.org/index.php?page=viewproblem2&cpid=698) - [Deltix Summer 2021 - DIY Tree](https://codeforces.com/contest/1556/problem/H) -- [AtCoder Contest Scheduling](https://atcoder.jp/contests/intro-heuristics/tasks/intro_heuristics_a) \ No newline at end of file +- [AtCoder Contest Scheduling](https://atcoder.jp/contests/intro-heuristics/tasks/intro_heuristics_a) From 9b8c5b638057a351ed79439d9ebd73002190fd81 Mon Sep 17 00:00:00 2001 From: Oleksandr Kulkov Date: Mon, 24 Mar 2025 21:33:53 +0100 Subject: [PATCH 09/22] Close
tag --- src/num_methods/simulated_annealing.md | 1 + 1 file changed, 1 insertion(+) diff --git a/src/num_methods/simulated_annealing.md b/src/num_methods/simulated_annealing.md index a5f6b466a..f9bd1686d 100644 --- a/src/num_methods/simulated_annealing.md +++ b/src/num_methods/simulated_annealing.md @@ -44,6 +44,7 @@ At the same time we also keep a track of the best state $s_{best}$ across all it A visual representation of simulated annealing, searching for the maxima of this function with multiple local maxima..
This animation by [Kingpin13](https://commons.wikimedia.org/wiki/User:Kingpin13) is distributed under CC0 1.0 license. +
### Temperature($T$) and decay($u$) From bee2358e0ea051328e6340babeeb57d1bd719c29 Mon Sep 17 00:00:00 2001 From: Proxihox Date: Thu, 1 May 2025 14:57:31 +0530 Subject: [PATCH 10/22] final fixes --- README.md | 2 +- src/num_methods/simulated_annealing.md | 16 +++++++--------- 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/README.md b/README.md index 87fce368e..4bd8eacfc 100644 --- a/README.md +++ b/README.md @@ -29,7 +29,7 @@ Compiled pages are published at [https://cp-algorithms.com/](https://cp-algorith ### New articles -- (19 October 2024) [Simulated Annealing](https://cp-algorithms.com/num_methods/simulated_annealing.html) +- (1 May 2025) [Simulated Annealing](https://cp-algorithms.com/num_methods/simulated_annealing.html) - (12 July 2024) [Manhattan distance](https://cp-algorithms.com/geometry/manhattan-distance.html) - (8 June 2024) [Knapsack Problem](https://cp-algorithms.com/dynamic_programming/knapsack.html) - (28 January 2024) [Introduction to Dynamic Programming](https://cp-algorithms.com/dynamic_programming/intro-to-dp.html) diff --git a/src/num_methods/simulated_annealing.md b/src/num_methods/simulated_annealing.md index f9bd1686d..054f71fa4 100644 --- a/src/num_methods/simulated_annealing.md +++ b/src/num_methods/simulated_annealing.md @@ -16,7 +16,7 @@ We are given a function $E(s)$, which calculates the energy of the state $s$. We You are given a set of nodes in 2 dimensional space. Each node is characterised by its $x$ and $y$ coordinates. Your task is to find an ordering of the nodes, which will minimise the distance to be travelled when visiting these nodes in that order. ## Motivation -Annealing is a metallurgical process , wherein a material is heated up and allowed to cool, in order to allow the atoms inside to rearrange themselves in an arrangement with minimal internal energy, which in turn causes the material to have different properties. The state is the arrangement of atoms and the internal energy is the function being minimised. We can think of the original state of the atoms, as a local minima for its internal energy. To make the material rearrange its atoms, we need to motivate it to go across a region where its internal energy is not minimised in order to reach the global minima. This motivation is given by heating the material to a higher temperature. +Annealing is a metallurgical process, wherein a material is heated up and allowed to cool, in order to allow the atoms inside to rearrange themselves in an arrangement with minimal internal energy, which in turn causes the material to have different properties. The state is the arrangement of atoms and the internal energy is the function being minimised. We can think of the original state of the atoms, as a local minima for its internal energy. To make the material rearrange its atoms, we need to motivate it to go across a region where its internal energy is not minimised in order to reach the global minima. This motivation is given by heating the material to a higher temperature. Simulated annealing, literally, simulates this process. We start off with some random state (material) and set a high temperature (heat it up). Now, the algorithm is ready to accept states which have a higher energy than the current state, as it is motivated by the high temperature. This prevents the algorithm from getting stuck inside local minimas and move towards the global minima. As time progresses, the algorithm cools down and refuses the states with higher energy and moves into the closest minima it has found. @@ -26,7 +26,7 @@ $E(s)$ is the function which needs to be minimised (or maximised). It maps every ### State -The state space is the domain of the energy function, E(s), and a state is any element which belongs to the state space. In the case of TSP, all possible paths that we can take to visit all the nodes is the state space, and any single one of these paths can be considered as a state. +The state space is the domain of the energy function, $E(s)$, and a state is any element which belongs to the state space. In the case of TSP, all possible paths that we can take to visit all the nodes is the state space, and any single one of these paths can be considered as a state. ### Neighbouring state @@ -41,12 +41,11 @@ At the same time we also keep a track of the best state $s_{best}$ across all it

-A visual representation of simulated annealing, searching for the maxima of this function with multiple local maxima.. +A visual representation of simulated annealing, searching for the maxima of this function with multiple local maxima.
-This animation by [Kingpin13](https://commons.wikimedia.org/wiki/User:Kingpin13) is distributed under CC0 1.0 license.
-### Temperature($T$) and decay($u$) +### Temperature(T) and decay(u) The temperature of the system quantifies the willingness of the algorithm to accept a state with a higher energy. The decay is a constant which quantifies the "cooling rate" of the algorithm. A slow cooling rate (larger $u$) is known to give better results. @@ -107,8 +106,8 @@ pair simAnneal() { best = s; E_best = E_next; } + E = E_next; } - E = E_next; T *= u; } return {E_best, best}; @@ -158,7 +157,6 @@ class state { double E() { double dist = 0; - bool first = true; int n = points.size(); for (int i = 0;i < n; i++) dist += euclidean(points[i], points[(i+1)%n]); @@ -187,8 +185,8 @@ int main() { - The effect of the difference in energies, $E_{next} - E$, on the PAF can be increased/decreased by increasing/decreasing the base of the exponent as shown below: ```cpp bool P(double E, double E_next, double T, mt19937 rng) { - e = 2 // set e to any real number greater than 1 - double prob = exp(-(E_next-E)/T); + double e = 2 // set e to any real number greater than 1 + double prob = pow(e,-(E_next-E)/T); if (prob > 1) return true; else { From 46e4efa5180bfa059694dec086269806832105a7 Mon Sep 17 00:00:00 2001 From: Shashank Sahu <52148284+bit-shashank@users.noreply.github.com> Date: Sat, 3 May 2025 00:07:17 +0530 Subject: [PATCH 11/22] Typo fix in graph/fixed_length_paths.md MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Replaced It is obvious that the constructed adjacency matrix if the answer to the problem for the case   $k = 1$ . It contains the number of paths of length   $1$  between each pair of vertices. To It is obvious that the constructed adjacency matrix is the answer to the problem for the case   $k = 1$ . It contains the number of paths of length   $1$  between each pair of vertices. --- src/graph/fixed_length_paths.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/graph/fixed_length_paths.md b/src/graph/fixed_length_paths.md index e5ecf76bc..9e1c1efad 100644 --- a/src/graph/fixed_length_paths.md +++ b/src/graph/fixed_length_paths.md @@ -21,7 +21,7 @@ The following algorithm works also in the case of multiple edges: if some pair of vertices $(i, j)$ is connected with $m$ edges, then we can record this in the adjacency matrix by setting $G[i][j] = m$. Also the algorithm works if the graph contains loops (a loop is an edge that connect a vertex with itself). -It is obvious that the constructed adjacency matrix if the answer to the problem for the case $k = 1$. +It is obvious that the constructed adjacency matrix is the answer to the problem for the case $k = 1$. It contains the number of paths of length $1$ between each pair of vertices. We will build the solution iteratively: From 4611aee0bdc668908f5def6a51228b6a0ed51619 Mon Sep 17 00:00:00 2001 From: Michael Hayter Date: Wed, 7 May 2025 17:57:39 -0400 Subject: [PATCH 12/22] Update CONTRIBUTING.md --- CONTRIBUTING.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index e1d3f2be1..fb6e7afee 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -19,7 +19,7 @@ Follow these steps to start contributing: 4. **Preview your changes** using the [preview page](preview.md) to ensure they look correct. 5. **Commit your changes** by clicking the _Propose changes_ button. 6. **Create a Pull Request (PR)** by clicking _Compare & pull request_. -7. **Review process**: Someone from the core team will review your changes. This may take a few hours to a few days. +7. **Review process**: Someone from the core team will review your changes. This may take a few days to a few weeks. ### Making Larger Changes From 27e90b5f3b84e58577f411d56a1e1582deb62850 Mon Sep 17 00:00:00 2001 From: t0wbo2t <52655804+t0wbo2t@users.noreply.github.com> Date: Thu, 15 May 2025 10:18:15 +0530 Subject: [PATCH 13/22] Update factorization.md [Update Powersmooth Definition] The definition of powersmooth was a bit confusing so I added a formal definition inspired by a number theory book. --- src/algebra/factorization.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algebra/factorization.md b/src/algebra/factorization.md index 11bf4049c..51f206482 100644 --- a/src/algebra/factorization.md +++ b/src/algebra/factorization.md @@ -160,7 +160,7 @@ By looking at the squares $a^2$ modulo a fixed small number, it can be observed ## Pollard's $p - 1$ method { data-toc-label="Pollard's method" } It is very likely that at least one factor of a number is $B$**-powersmooth** for small $B$. -$B$-powersmooth means that every prime power $d^k$ that divides $p-1$ is at most $B$. +$B$-powersmooth means that every prime power $d^k$ that divides $p-1$ is at most $B$. Formally, let $\mathrm{B} \geqslant 1$ and $n \geqslant 1$ with prime factorization $n = \prod {p_i}^{e_i},$ then $n$ is $\mathrm{B}$-powersmooth if, for all $i,$ ${p_i}^{e_i} \leqslant \mathrm{B}$. E.g. the prime factorization of $4817191$ is $1303 \cdot 3697$. And the factors are $31$-powersmooth and $16$-powersmooth respectably, because $1303 - 1 = 2 \cdot 3 \cdot 7 \cdot 31$ and $3697 - 1 = 2^4 \cdot 3 \cdot 7 \cdot 11$. In 1974 John Pollard invented a method to extracts $B$-powersmooth factors from a composite number. From 6ec2fe6634b6e2bd4d879d08dfdf3bbcf87d0d6b Mon Sep 17 00:00:00 2001 From: 100daysummer <138024460+100daysummer@users.noreply.github.com> Date: Wed, 21 May 2025 12:00:17 +0300 Subject: [PATCH 14/22] Update longest_increasing_subsequence.md Small improvement of the wording of a question --- src/sequences/longest_increasing_subsequence.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sequences/longest_increasing_subsequence.md b/src/sequences/longest_increasing_subsequence.md index a4c8f0fe4..dca4f020b 100644 --- a/src/sequences/longest_increasing_subsequence.md +++ b/src/sequences/longest_increasing_subsequence.md @@ -174,7 +174,7 @@ We will again gradually process the numbers, first $a[0]$, then $a[1]$, etc, and $$ When we process $a[i]$, we can ask ourselves. -What have the conditions to be, that we write the current number $a[i]$ into the $d[0 \dots n]$ array? +Under what conditions should we write the current number $a[i]$ into the $d[0 \dots n]$ array? We set $d[l] = a[i]$, if there is a longest increasing sequence of length $l$ that ends in $a[i]$, and there is no longest increasing sequence of length $l$ that ends in a smaller number. Similar to the previous approach, if we remove the number $a[i]$ from the longest increasing sequence of length $l$, we get another longest increasing sequence of length $l -1$. From d1cfad25305a0047c822116f4369bdb33d157364 Mon Sep 17 00:00:00 2001 From: Oleksandr Kulkov Date: Wed, 21 May 2025 13:11:00 +0200 Subject: [PATCH 15/22] semicolon after e --- src/num_methods/simulated_annealing.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/num_methods/simulated_annealing.md b/src/num_methods/simulated_annealing.md index 054f71fa4..bd0dd6ed2 100644 --- a/src/num_methods/simulated_annealing.md +++ b/src/num_methods/simulated_annealing.md @@ -185,7 +185,7 @@ int main() { - The effect of the difference in energies, $E_{next} - E$, on the PAF can be increased/decreased by increasing/decreasing the base of the exponent as shown below: ```cpp bool P(double E, double E_next, double T, mt19937 rng) { - double e = 2 // set e to any real number greater than 1 + double e = 2; // set e to any real number greater than 1 double prob = pow(e,-(E_next-E)/T); if (prob > 1) return true; From a82cd128a8f3040c8e2658c317f853f9d0ec2d8d Mon Sep 17 00:00:00 2001 From: Oleksandr Kulkov Date: Wed, 21 May 2025 13:11:36 +0200 Subject: [PATCH 16/22] update date --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 4bd8eacfc..513e55f32 100644 --- a/README.md +++ b/README.md @@ -29,7 +29,7 @@ Compiled pages are published at [https://cp-algorithms.com/](https://cp-algorith ### New articles -- (1 May 2025) [Simulated Annealing](https://cp-algorithms.com/num_methods/simulated_annealing.html) +- (21 May 2025) [Simulated Annealing](https://cp-algorithms.com/num_methods/simulated_annealing.html) - (12 July 2024) [Manhattan distance](https://cp-algorithms.com/geometry/manhattan-distance.html) - (8 June 2024) [Knapsack Problem](https://cp-algorithms.com/dynamic_programming/knapsack.html) - (28 January 2024) [Introduction to Dynamic Programming](https://cp-algorithms.com/dynamic_programming/intro-to-dp.html) From 422fb19bdd7c39a56d8fb51caefa419251b37fab Mon Sep 17 00:00:00 2001 From: t0wbo2t <52655804+t0wbo2t@users.noreply.github.com> Date: Wed, 21 May 2025 17:33:44 +0530 Subject: [PATCH 17/22] Update factorization.md [Update powersmooth definition with suggestions] Commas are now outside $...$ and definiton is for (p - 1) for consistency. --- src/algebra/factorization.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algebra/factorization.md b/src/algebra/factorization.md index 51f206482..997a11be2 100644 --- a/src/algebra/factorization.md +++ b/src/algebra/factorization.md @@ -160,7 +160,7 @@ By looking at the squares $a^2$ modulo a fixed small number, it can be observed ## Pollard's $p - 1$ method { data-toc-label="Pollard's method" } It is very likely that at least one factor of a number is $B$**-powersmooth** for small $B$. -$B$-powersmooth means that every prime power $d^k$ that divides $p-1$ is at most $B$. Formally, let $\mathrm{B} \geqslant 1$ and $n \geqslant 1$ with prime factorization $n = \prod {p_i}^{e_i},$ then $n$ is $\mathrm{B}$-powersmooth if, for all $i,$ ${p_i}^{e_i} \leqslant \mathrm{B}$. +$B$-powersmooth means that every prime power $d^k$ that divides $p-1$ is at most $B$. Formally, let $\mathrm{B} \geqslant 1$ and let $p$ be a prime such that $(p - 1) \geqslant 1$. Suppose the prime factorization of $(p - 1)$ is $(p - 1) = \prod {q_i}^{e_i}$, where each $q_i$ is a prime and $e_i \geqslant 1$ then $(p - 1)$ is $\mathrm{B}$-powersmooth if, for all $i$, ${q_i}^{e_i} \leqslant \mathrm{B}$. E.g. the prime factorization of $4817191$ is $1303 \cdot 3697$. And the factors are $31$-powersmooth and $16$-powersmooth respectably, because $1303 - 1 = 2 \cdot 3 \cdot 7 \cdot 31$ and $3697 - 1 = 2^4 \cdot 3 \cdot 7 \cdot 11$. In 1974 John Pollard invented a method to extracts $B$-powersmooth factors from a composite number. From 1a7db31f720455b647a929596ff70936875f7532 Mon Sep 17 00:00:00 2001 From: t0wbo2t <52655804+t0wbo2t@users.noreply.github.com> Date: Thu, 22 May 2025 09:24:40 +0530 Subject: [PATCH 18/22] Update factorization.md [Update Pollard's (p - 1) Method] Updated materials related to powersmoothness. Corrected some minor mistakes. --- src/algebra/factorization.md | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/algebra/factorization.md b/src/algebra/factorization.md index 997a11be2..84bdd356d 100644 --- a/src/algebra/factorization.md +++ b/src/algebra/factorization.md @@ -159,11 +159,10 @@ By looking at the squares $a^2$ modulo a fixed small number, it can be observed ## Pollard's $p - 1$ method { data-toc-label="Pollard's method" } -It is very likely that at least one factor of a number is $B$**-powersmooth** for small $B$. -$B$-powersmooth means that every prime power $d^k$ that divides $p-1$ is at most $B$. Formally, let $\mathrm{B} \geqslant 1$ and let $p$ be a prime such that $(p - 1) \geqslant 1$. Suppose the prime factorization of $(p - 1)$ is $(p - 1) = \prod {q_i}^{e_i}$, where each $q_i$ is a prime and $e_i \geqslant 1$ then $(p - 1)$ is $\mathrm{B}$-powersmooth if, for all $i$, ${q_i}^{e_i} \leqslant \mathrm{B}$. +It is very likely that a number $n$ has at least one prime factor $p$ such that $p - 1$ is $\mathrm{B}$**-powersmooth** for small $\mathrm{B}$. An integer $m$ is said to be $\mathrm{B}$-powersmooth if every prime power dividing $m$ is at most $\mathrm{B}$. Formally, let $\mathrm{B} \geqslant 1$ and let $m$ be any positive integer. Suppose the prime factorization of $m$ is $m = \prod {q_i}^{e_i}$, where each $q_i$ is a prime and $e_i \geqslant 1$. Then $m$ is $\mathrm{B}$-powersmooth if, for all $i$, ${q_i}^{e_i} \leqslant \mathrm{B}$. E.g. the prime factorization of $4817191$ is $1303 \cdot 3697$. -And the factors are $31$-powersmooth and $16$-powersmooth respectably, because $1303 - 1 = 2 \cdot 3 \cdot 7 \cdot 31$ and $3697 - 1 = 2^4 \cdot 3 \cdot 7 \cdot 11$. -In 1974 John Pollard invented a method to extracts $B$-powersmooth factors from a composite number. +And the values, $1303 - 1$ and $3697 - 1$, are $31$-powersmooth and $16$-powersmooth respectively, because $1303 - 1 = 2 \cdot 3 \cdot 7 \cdot 31$ and $3697 - 1 = 2^4 \cdot 3 \cdot 7 \cdot 11$. +In 1974 John Pollard invented a method to extracts $\mathrm{B}$-powersmooth factors from a composite number. The idea comes from [Fermat's little theorem](phi-function.md#application). Let a factorization of $n$ be $n = p \cdot q$. @@ -180,7 +179,7 @@ This means that $a^M - 1 = p \cdot r$, and because of that also $p ~|~ \gcd(a^M Therefore, if $p - 1$ for a factor $p$ of $n$ divides $M$, we can extract a factor using [Euclid's algorithm](euclid-algorithm.md). -It is clear, that the smallest $M$ that is a multiple of every $B$-powersmooth number is $\text{lcm}(1,~2~,3~,4~,~\dots,~B)$. +It is clear, that the smallest $M$ that is a multiple of every $\mathrm{B}$-powersmooth number is $\text{lcm}(1,~2~,3~,4~,~\dots,~B)$. Or alternatively: $$M = \prod_{\text{prime } q \le B} q^{\lfloor \log_q B \rfloor}$$ @@ -189,11 +188,11 @@ Notice, if $p-1$ divides $M$ for all prime factors $p$ of $n$, then $\gcd(a^M - In this case we don't receive a factor. Therefore, we will try to perform the $\gcd$ multiple times, while we compute $M$. -Some composite numbers don't have $B$-powersmooth factors for small $B$. +Some composite numbers don't have $\mathrm{B}$-powersmooth factors for small $\mathrm{B}$. For example, the factors of the composite number $100~000~000~000~000~493 = 763~013 \cdot 131~059~365~961$ are $190~753$-powersmooth and $1~092~161~383$-powersmooth. We will have to choose $B >= 190~753$ to factorize the number. -In the following implementation we start with $B = 10$ and increase $B$ after each each iteration. +In the following implementation we start with $\mathrm{B} = 10$ and increase $\mathrm{B}$ after each each iteration. ```{.cpp file=factorization_p_minus_1} long long pollards_p_minus_1(long long n) { From bb7e13757ec14e867b3cd4de3da8966d87417270 Mon Sep 17 00:00:00 2001 From: Oleksandr Kulkov Date: Sat, 24 May 2025 20:58:13 +0200 Subject: [PATCH 19/22] fix #1372 --- src/string/manacher.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/string/manacher.md b/src/string/manacher.md index 0c8bd5928..26589b83c 100644 --- a/src/string/manacher.md +++ b/src/string/manacher.md @@ -147,7 +147,7 @@ vector manacher_odd(string s) { vector p(n + 2); int l = 0, r = 1; for(int i = 1; i <= n; i++) { - p[i] = max(0, min(r - i, p[l + (r - i)])); + p[i] = min(r - i, p[l + (r - i)]); while(s[i - p[i]] == s[i + p[i]]) { p[i]++; } From 2597e0558304678a7fa7f92c66a19f9f7913c974 Mon Sep 17 00:00:00 2001 From: t0wbo2t <52655804+t0wbo2t@users.noreply.github.com> Date: Thu, 29 May 2025 11:19:43 +0530 Subject: [PATCH 20/22] Update src/algebra/factorization.md Co-authored-by: Oleksandr Kulkov --- src/algebra/factorization.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/algebra/factorization.md b/src/algebra/factorization.md index 84bdd356d..bc606607f 100644 --- a/src/algebra/factorization.md +++ b/src/algebra/factorization.md @@ -188,8 +188,8 @@ Notice, if $p-1$ divides $M$ for all prime factors $p$ of $n$, then $\gcd(a^M - In this case we don't receive a factor. Therefore, we will try to perform the $\gcd$ multiple times, while we compute $M$. -Some composite numbers don't have $\mathrm{B}$-powersmooth factors for small $\mathrm{B}$. -For example, the factors of the composite number $100~000~000~000~000~493 = 763~013 \cdot 131~059~365~961$ are $190~753$-powersmooth and $1~092~161~383$-powersmooth. +Some composite numbers don't have factors $p$ s.t. $p-1$ is $\mathrm{B}$-powersmooth for small $\mathrm{B}$. +For example, for the composite number $100~000~000~000~000~493 = 763~013 \cdot 131~059~365~961$, values $p-1$ are $190~753$-powersmooth and $1~092~161~383$-powersmooth correspondingly. We will have to choose $B >= 190~753$ to factorize the number. In the following implementation we start with $\mathrm{B} = 10$ and increase $\mathrm{B}$ after each each iteration. From 54fec62526c6454ecd43b38544c6a56880031544 Mon Sep 17 00:00:00 2001 From: t0wbo2t <52655804+t0wbo2t@users.noreply.github.com> Date: Thu, 29 May 2025 11:19:54 +0530 Subject: [PATCH 21/22] Update src/algebra/factorization.md Co-authored-by: Oleksandr Kulkov --- src/algebra/factorization.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algebra/factorization.md b/src/algebra/factorization.md index bc606607f..9d9ab7ed7 100644 --- a/src/algebra/factorization.md +++ b/src/algebra/factorization.md @@ -162,7 +162,7 @@ By looking at the squares $a^2$ modulo a fixed small number, it can be observed It is very likely that a number $n$ has at least one prime factor $p$ such that $p - 1$ is $\mathrm{B}$**-powersmooth** for small $\mathrm{B}$. An integer $m$ is said to be $\mathrm{B}$-powersmooth if every prime power dividing $m$ is at most $\mathrm{B}$. Formally, let $\mathrm{B} \geqslant 1$ and let $m$ be any positive integer. Suppose the prime factorization of $m$ is $m = \prod {q_i}^{e_i}$, where each $q_i$ is a prime and $e_i \geqslant 1$. Then $m$ is $\mathrm{B}$-powersmooth if, for all $i$, ${q_i}^{e_i} \leqslant \mathrm{B}$. E.g. the prime factorization of $4817191$ is $1303 \cdot 3697$. And the values, $1303 - 1$ and $3697 - 1$, are $31$-powersmooth and $16$-powersmooth respectively, because $1303 - 1 = 2 \cdot 3 \cdot 7 \cdot 31$ and $3697 - 1 = 2^4 \cdot 3 \cdot 7 \cdot 11$. -In 1974 John Pollard invented a method to extracts $\mathrm{B}$-powersmooth factors from a composite number. +In 1974 John Pollard invented a method to extract factors $p$, s.t. $p-1$ is $\mathrm{B}$-powersmooth, from a composite number. The idea comes from [Fermat's little theorem](phi-function.md#application). Let a factorization of $n$ be $n = p \cdot q$. From 53639d500284ae160ca2e9f968957c360b6f2040 Mon Sep 17 00:00:00 2001 From: t0wbo2t <52655804+t0wbo2t@users.noreply.github.com> Date: Thu, 29 May 2025 11:20:00 +0530 Subject: [PATCH 22/22] Update src/algebra/factorization.md Co-authored-by: Oleksandr Kulkov --- src/algebra/factorization.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algebra/factorization.md b/src/algebra/factorization.md index 9d9ab7ed7..14715605f 100644 --- a/src/algebra/factorization.md +++ b/src/algebra/factorization.md @@ -190,7 +190,7 @@ Therefore, we will try to perform the $\gcd$ multiple times, while we compute $M Some composite numbers don't have factors $p$ s.t. $p-1$ is $\mathrm{B}$-powersmooth for small $\mathrm{B}$. For example, for the composite number $100~000~000~000~000~493 = 763~013 \cdot 131~059~365~961$, values $p-1$ are $190~753$-powersmooth and $1~092~161~383$-powersmooth correspondingly. -We will have to choose $B >= 190~753$ to factorize the number. +We will have to choose $B \geq 190~753$ to factorize the number. In the following implementation we start with $\mathrm{B} = 10$ and increase $\mathrm{B}$ after each each iteration.