pyduino.optimization.nelder_mead 

  1import numpy as np
  2import warnings
  3from typing import Union, List, Iterable
  4from . import Optimizer, linmap
  5
  6def sigmoid(x):
  7    return 1 / (1 + np.exp(-x))
  8def logit(x, inf=100):
  9    return np.where(x==0, -inf, np.where(x==1, inf, np.log(x) - np.log(1-x)))
 10def softmax(x):
 11    assert x.ndim == 2, "Softmax only implemented for 2D arrays"
 12    return np.exp(x) / np.sum(np.exp(x), axis=1, keepdims=True)
 13
 14class NelderMeadBounded(Optimizer):
 15    def __init__(self, population_size: int, ranges: List[float], rng_seed: int = 0, hypercube_radius = 100):
 16        """
 17        This Nelder-Mead algorithm assumes that the optimization function 'f' is to be minimized and time independent.
 18
 19        The space is assumed to be a hypercube, and the algorithm will map the hypercube to the unit hypercube, bounded by a sigmoid function.
 20
 21        ranges list of pairs:
 22            Maxima and minima that each parameter can cover.
 23            Example: [(0, 10), (0, 10)] for two parameters ranging from 0 to 10.
 24        population_size int:
 25            Number of individuals in the population.
 26        rng_seed int:
 27            Seed for the random number generator. Default is 0.
 28        """
 29        self.population_size = population_size
 30        self.ranges = np.array(ranges)
 31        self.rng_seed = np.random.default_rng(rng_seed)
 32
 33        # Derived attributes
 34        self.a = 1/hypercube_radius
 35        self.backward = lambda x: logit(linmap(self.ranges, np.array([[0, 1]] * len(self.ranges)))(x))/self.a
 36        self.forward = lambda x: linmap(np.array([[0, 1]] * len(self.ranges)), self.ranges)(sigmoid(self.a * x))
 37
 38        # Initialize the population (random position and initial momentum)
 39        self.population = self.rng_seed.random((self.population_size, len(self.ranges)))
 40        
 41        # Initialize y as vector of nans
 42        self.y = np.full(self.population_size, np.nan)
 43
 44    def view(self, x):
 45        """
 46        Maps the input from the domain to the codomain.
 47        """
 48        return self.forward(x)
 49    
 50    def view_g(self):
 51        """
 52        Maps the input from the domain to the codomain.
 53        """
 54        return self.forward(self.population)
 55
 56    def inverse_view(self, x):
 57        """
 58        Maps the input from the codomain to the domain.
 59        """
 60        return self.backward(x)
 61    
 62    def ask_oracle(self, X: np.ndarray) -> np.ndarray:
 63        return super().ask_oracle()
 64
 65    def init_oracle(self):
 66        return self.ask_oracle(self.view_g())
 67
 68    def step(self):
 69        """
 70        This function performs a single step of the Nelder-Mead algorithm.
 71        """
 72        
 73        # Sort the population by the value of the oracle
 74        y = self.ask_oracle(self.view_g())
 75        idx = np.argsort(y)
 76        self.population = self.population[idx]
 77
 78        # Compute the centroid of the population
 79        centroid = self.population[:-1].mean(axis=0)
 80
 81        # Reflect the worst point through the centroid
 82        reflected = centroid + (centroid - self.population[-1])
 83
 84        # Evaluate the reflected point
 85        
 86        y_reflected = self.ask_oracle(self.view(reflected.reshape(1,-1)))
 87
 88        # If the reflected point is better than the second worst, but not better than the best, then expand
 89        
 90        if y_reflected < y[-2] and y_reflected > y[0]:
 91            expanded = centroid + (reflected - centroid)
 92            y_expanded = self.ask_oracle(self.view(expanded.reshape(1,-1)))
 93            if y_expanded < y_reflected:
 94                self.population[-1] = expanded
 95            else:
 96                self.population[-1] = reflected
 97        # If the reflected point is better than the best, then expand
 98        elif y_reflected < y[0]:
 99            expanded = centroid + 2 * (reflected - centroid)
100            y_expanded = self.ask_oracle(self.view(expanded.reshape(1,-1)))
101            if y_expanded < y_reflected:
102                self.population[-1] = expanded
103            else:
104                self.population[-1] = reflected
105        # If the reflected point is worse than the second worst, then contract
106        elif y_reflected > y[-2]:
107            contracted = centroid + 0.5 * (self.population[-1] - centroid)
108            y_contracted = self.ask_oracle(self.view(contracted.reshape(1,-1)))
109            if y_contracted < y[-1]:
110                self.population[-1] = contracted
111            else:
112                for i in range(1, len(self.population)):
113                    self.population[i] = 0.5 * (self.population[i] + self.population[0])
114        # If the reflected point is worse than the worst, then shrink
115        elif y_reflected > y[-1]:
116            for i in range(1, len(self.population)):
117                self.population[i] = 0.5 * (self.population[i] + self.population[0])
118
119        self.y = y.copy()
120        return self.view_g()
121
122class NelderMeadConstant(Optimizer):
123    def __init__(self, population_size: int, ranges: Union[int, Iterable], rng_seed: int = 0, energy=10):
124        """
125        This Nelder-Mead algorithm assumes that the optimization function 'f' is to be minimized and time independent.
126
127        The parameters are constrained so that their sum always add up to `energy` though a softmax function.
128
129        Parameters:
130        - population_size (int): The size of the population.
131        - ranges (int or list): The number of parameters or a list of ranges for each parameter.
132        - rng_seed (int): The seed for the random number generator.
133        - energy (float): The energy parameter used in the optimization.
134
135        Note: The ranges parameter is merely a placeholder. The ranges are set to (0, energy) for all parameters.
136        """
137        self.population_size = population_size
138        if isinstance(ranges, Iterable):
139            self.ranges = np.array([(0, energy) for _ in range(len(ranges))])
140            warnings.warn("While using Nelder-MeadConstant, the ranges are set to (0, energy) for all parameters. The parameter `ranges` is merely a placeholder.")
141        elif isinstance(ranges, int):
142            self.ranges = np.array([(0, energy) for _ in range(ranges)])
143        self.rng_seed = np.random.default_rng(rng_seed)
144
145        # Initialize the population (random position and initial momentum)
146        self.population = self.rng_seed.random((self.population_size, len(self.ranges)))
147
148        # Derived attributes
149        self.energy = energy
150        self.e_summation = 1
151        # Initialize y as vector of nans
152        self.y = np.full(self.population_size, np.nan)
153    
154    def forward(self, x):
155        self.e_summation = np.sum(np.exp(x), axis=1, keepdims=True)
156        return self.energy * softmax(x)
157    
158    def backward(self, x):
159        """
160        Softmax is not injective. This is a pseudo-inverse.
161        """
162        return np.log(self.e_summation * x/self.energy)
163
164    def view(self, x):
165        """
166        Maps the input from the domain to the codomain.
167        """
168        return self.forward(x)
169    
170    def view_g(self):
171        """
172        Maps the input from the domain to the codomain.
173        """
174        return self.forward(self.population)
175
176    def inverse_view(self, x):
177        """
178        Maps the input from the codomain to the domain.
179        """
180        return self.backward(x)
181    
182    def ask_oracle(self, X: np.ndarray) -> np.ndarray:
183        return super().ask_oracle()
184
185    def init_oracle(self):
186        return self.ask_oracle(self.view_g())
187
188    def step(self):
189        """
190        This function performs a single step of the Nelder-Mead algorithm.
191        """
192        
193        # Sort the population by the value of the oracle
194        y = self.ask_oracle(self.view_g())
195        idx = np.argsort(y)
196        self.population = self.population[idx]
197
198        # Compute the centroid of the population
199        centroid = self.population[:-1].mean(axis=0)
200
201        # Reflect the worst point through the centroid
202        reflected = centroid + (centroid - self.population[-1])
203
204        # Evaluate the reflected point
205        
206        y_reflected = self.ask_oracle(self.view(reflected.reshape(1,-1)))
207
208        # If the reflected point is better than the second worst, but not better than the best, then expand
209        
210        if y_reflected < y[-2] and y_reflected > y[0]:
211            expanded = centroid + (reflected - centroid)
212            y_expanded = self.ask_oracle(self.view(expanded.reshape(1,-1)))
213            if y_expanded < y_reflected:
214                self.population[-1] = expanded
215            else:
216                self.population[-1] = reflected
217        # If the reflected point is better than the best, then expand
218        elif y_reflected < y[0]:
219            expanded = centroid + 2 * (reflected - centroid)
220            y_expanded = self.ask_oracle(self.view(expanded.reshape(1,-1)))
221            if y_expanded < y_reflected:
222                self.population[-1] = expanded
223            else:
224                self.population[-1] = reflected
225        # If the reflected point is worse than the second worst, then contract
226        elif y_reflected > y[-2]:
227            contracted = centroid + 0.5 * (self.population[-1] - centroid)
228            y_contracted = self.ask_oracle(self.view(contracted.reshape(1,-1)))
229            if y_contracted < y[-1]:
230                self.population[-1] = contracted
231            else:
232                for i in range(1, len(self.population)):
233                    self.population[i] = 0.5 * (self.population[i] + self.population[0])
234        # If the reflected point is worse than the worst, then shrink
235        elif y_reflected > y[-1]:
236            for i in range(1, len(self.population)):
237                self.population[i] = 0.5 * (self.population[i] + self.population[0])
238
239        self.y = y.copy()
240        return self.view_g()
def sigmoid(x): 
7def sigmoid(x):
8    return 1 / (1 + np.exp(-x))
def logit(x, inf=100): 
 9def logit(x, inf=100):
10    return np.where(x==0, -inf, np.where(x==1, inf, np.log(x) - np.log(1-x)))
def softmax(x): 
11def softmax(x):
12    assert x.ndim == 2, "Softmax only implemented for 2D arrays"
13    return np.exp(x) / np.sum(np.exp(x), axis=1, keepdims=True)
class NelderMeadBounded(pyduino.optimization.Optimizer): 
 15class NelderMeadBounded(Optimizer):
 16    def __init__(self, population_size: int, ranges: List[float], rng_seed: int = 0, hypercube_radius = 100):
 17        """
 18        This Nelder-Mead algorithm assumes that the optimization function 'f' is to be minimized and time independent.
 19
 20        The space is assumed to be a hypercube, and the algorithm will map the hypercube to the unit hypercube, bounded by a sigmoid function.
 21
 22        ranges list of pairs:
 23            Maxima and minima that each parameter can cover.
 24            Example: [(0, 10), (0, 10)] for two parameters ranging from 0 to 10.
 25        population_size int:
 26            Number of individuals in the population.
 27        rng_seed int:
 28            Seed for the random number generator. Default is 0.
 29        """
 30        self.population_size = population_size
 31        self.ranges = np.array(ranges)
 32        self.rng_seed = np.random.default_rng(rng_seed)
 33
 34        # Derived attributes
 35        self.a = 1/hypercube_radius
 36        self.backward = lambda x: logit(linmap(self.ranges, np.array([[0, 1]] * len(self.ranges)))(x))/self.a
 37        self.forward = lambda x: linmap(np.array([[0, 1]] * len(self.ranges)), self.ranges)(sigmoid(self.a * x))
 38
 39        # Initialize the population (random position and initial momentum)
 40        self.population = self.rng_seed.random((self.population_size, len(self.ranges)))
 41        
 42        # Initialize y as vector of nans
 43        self.y = np.full(self.population_size, np.nan)
 44
 45    def view(self, x):
 46        """
 47        Maps the input from the domain to the codomain.
 48        """
 49        return self.forward(x)
 50    
 51    def view_g(self):
 52        """
 53        Maps the input from the domain to the codomain.
 54        """
 55        return self.forward(self.population)
 56
 57    def inverse_view(self, x):
 58        """
 59        Maps the input from the codomain to the domain.
 60        """
 61        return self.backward(x)
 62    
 63    def ask_oracle(self, X: np.ndarray) -> np.ndarray:
 64        return super().ask_oracle()
 65
 66    def init_oracle(self):
 67        return self.ask_oracle(self.view_g())
 68
 69    def step(self):
 70        """
 71        This function performs a single step of the Nelder-Mead algorithm.
 72        """
 73        
 74        # Sort the population by the value of the oracle
 75        y = self.ask_oracle(self.view_g())
 76        idx = np.argsort(y)
 77        self.population = self.population[idx]
 78
 79        # Compute the centroid of the population
 80        centroid = self.population[:-1].mean(axis=0)
 81
 82        # Reflect the worst point through the centroid
 83        reflected = centroid + (centroid - self.population[-1])
 84
 85        # Evaluate the reflected point
 86        
 87        y_reflected = self.ask_oracle(self.view(reflected.reshape(1,-1)))
 88
 89        # If the reflected point is better than the second worst, but not better than the best, then expand
 90        
 91        if y_reflected < y[-2] and y_reflected > y[0]:
 92            expanded = centroid + (reflected - centroid)
 93            y_expanded = self.ask_oracle(self.view(expanded.reshape(1,-1)))
 94            if y_expanded < y_reflected:
 95                self.population[-1] = expanded
 96            else:
 97                self.population[-1] = reflected
 98        # If the reflected point is better than the best, then expand
 99        elif y_reflected < y[0]:
100            expanded = centroid + 2 * (reflected - centroid)
101            y_expanded = self.ask_oracle(self.view(expanded.reshape(1,-1)))
102            if y_expanded < y_reflected:
103                self.population[-1] = expanded
104            else:
105                self.population[-1] = reflected
106        # If the reflected point is worse than the second worst, then contract
107        elif y_reflected > y[-2]:
108            contracted = centroid + 0.5 * (self.population[-1] - centroid)
109            y_contracted = self.ask_oracle(self.view(contracted.reshape(1,-1)))
110            if y_contracted < y[-1]:
111                self.population[-1] = contracted
112            else:
113                for i in range(1, len(self.population)):
114                    self.population[i] = 0.5 * (self.population[i] + self.population[0])
115        # If the reflected point is worse than the worst, then shrink
116        elif y_reflected > y[-1]:
117            for i in range(1, len(self.population)):
118                self.population[i] = 0.5 * (self.population[i] + self.population[0])
119
120        self.y = y.copy()
121        return self.view_g()

Abstract class for optimization algorithms

NelderMeadBounded( population_size: int, ranges: List[float], rng_seed: int = 0, hypercube_radius=100) 
16    def __init__(self, population_size: int, ranges: List[float], rng_seed: int = 0, hypercube_radius = 100):
17        """
18        This Nelder-Mead algorithm assumes that the optimization function 'f' is to be minimized and time independent.
19
20        The space is assumed to be a hypercube, and the algorithm will map the hypercube to the unit hypercube, bounded by a sigmoid function.
21
22        ranges list of pairs:
23            Maxima and minima that each parameter can cover.
24            Example: [(0, 10), (0, 10)] for two parameters ranging from 0 to 10.
25        population_size int:
26            Number of individuals in the population.
27        rng_seed int:
28            Seed for the random number generator. Default is 0.
29        """
30        self.population_size = population_size
31        self.ranges = np.array(ranges)
32        self.rng_seed = np.random.default_rng(rng_seed)
33
34        # Derived attributes
35        self.a = 1/hypercube_radius
36        self.backward = lambda x: logit(linmap(self.ranges, np.array([[0, 1]] * len(self.ranges)))(x))/self.a
37        self.forward = lambda x: linmap(np.array([[0, 1]] * len(self.ranges)), self.ranges)(sigmoid(self.a * x))
38
39        # Initialize the population (random position and initial momentum)
40        self.population = self.rng_seed.random((self.population_size, len(self.ranges)))
41        
42        # Initialize y as vector of nans
43        self.y = np.full(self.population_size, np.nan)

This Nelder-Mead algorithm assumes that the optimization function 'f' is to be minimized and time independent.

The space is assumed to be a hypercube, and the algorithm will map the hypercube to the unit hypercube, bounded by a sigmoid function.

ranges list of pairs: Maxima and minima that each parameter can cover. Example: [(0, 10), (0, 10)] for two parameters ranging from 0 to 10. population_size int: Number of individuals in the population. rng_seed int: Seed for the random number generator. Default is 0.

population_size
ranges
rng_seed
a
backward
forward
population
y
def view(self, x): 
45    def view(self, x):
46        """
47        Maps the input from the domain to the codomain.
48        """
49        return self.forward(x)

Maps the input from the domain to the codomain.

def view_g(self): 
51    def view_g(self):
52        """
53        Maps the input from the domain to the codomain.
54        """
55        return self.forward(self.population)

Maps the input from the domain to the codomain.

def inverse_view(self, x): 
57    def inverse_view(self, x):
58        """
59        Maps the input from the codomain to the domain.
60        """
61        return self.backward(x)

Maps the input from the codomain to the domain.

def ask_oracle(self, X: numpy.ndarray) -> numpy.ndarray: 
63    def ask_oracle(self, X: np.ndarray) -> np.ndarray:
64        return super().ask_oracle()
def init_oracle(self): 
66    def init_oracle(self):
67        return self.ask_oracle(self.view_g())
def step(self): 
 69    def step(self):
 70        """
 71        This function performs a single step of the Nelder-Mead algorithm.
 72        """
 73        
 74        # Sort the population by the value of the oracle
 75        y = self.ask_oracle(self.view_g())
 76        idx = np.argsort(y)
 77        self.population = self.population[idx]
 78
 79        # Compute the centroid of the population
 80        centroid = self.population[:-1].mean(axis=0)
 81
 82        # Reflect the worst point through the centroid
 83        reflected = centroid + (centroid - self.population[-1])
 84
 85        # Evaluate the reflected point
 86        
 87        y_reflected = self.ask_oracle(self.view(reflected.reshape(1,-1)))
 88
 89        # If the reflected point is better than the second worst, but not better than the best, then expand
 90        
 91        if y_reflected < y[-2] and y_reflected > y[0]:
 92            expanded = centroid + (reflected - centroid)
 93            y_expanded = self.ask_oracle(self.view(expanded.reshape(1,-1)))
 94            if y_expanded < y_reflected:
 95                self.population[-1] = expanded
 96            else:
 97                self.population[-1] = reflected
 98        # If the reflected point is better than the best, then expand
 99        elif y_reflected < y[0]:
100            expanded = centroid + 2 * (reflected - centroid)
101            y_expanded = self.ask_oracle(self.view(expanded.reshape(1,-1)))
102            if y_expanded < y_reflected:
103                self.population[-1] = expanded
104            else:
105                self.population[-1] = reflected
106        # If the reflected point is worse than the second worst, then contract
107        elif y_reflected > y[-2]:
108            contracted = centroid + 0.5 * (self.population[-1] - centroid)
109            y_contracted = self.ask_oracle(self.view(contracted.reshape(1,-1)))
110            if y_contracted < y[-1]:
111                self.population[-1] = contracted
112            else:
113                for i in range(1, len(self.population)):
114                    self.population[i] = 0.5 * (self.population[i] + self.population[0])
115        # If the reflected point is worse than the worst, then shrink
116        elif y_reflected > y[-1]:
117            for i in range(1, len(self.population)):
118                self.population[i] = 0.5 * (self.population[i] + self.population[0])
119
120        self.y = y.copy()
121        return self.view_g()

This function performs a single step of the Nelder-Mead algorithm.

class NelderMeadConstant(pyduino.optimization.Optimizer): 
123class NelderMeadConstant(Optimizer):
124    def __init__(self, population_size: int, ranges: Union[int, Iterable], rng_seed: int = 0, energy=10):
125        """
126        This Nelder-Mead algorithm assumes that the optimization function 'f' is to be minimized and time independent.
127
128        The parameters are constrained so that their sum always add up to `energy` though a softmax function.
129
130        Parameters:
131        - population_size (int): The size of the population.
132        - ranges (int or list): The number of parameters or a list of ranges for each parameter.
133        - rng_seed (int): The seed for the random number generator.
134        - energy (float): The energy parameter used in the optimization.
135
136        Note: The ranges parameter is merely a placeholder. The ranges are set to (0, energy) for all parameters.
137        """
138        self.population_size = population_size
139        if isinstance(ranges, Iterable):
140            self.ranges = np.array([(0, energy) for _ in range(len(ranges))])
141            warnings.warn("While using Nelder-MeadConstant, the ranges are set to (0, energy) for all parameters. The parameter `ranges` is merely a placeholder.")
142        elif isinstance(ranges, int):
143            self.ranges = np.array([(0, energy) for _ in range(ranges)])
144        self.rng_seed = np.random.default_rng(rng_seed)
145
146        # Initialize the population (random position and initial momentum)
147        self.population = self.rng_seed.random((self.population_size, len(self.ranges)))
148
149        # Derived attributes
150        self.energy = energy
151        self.e_summation = 1
152        # Initialize y as vector of nans
153        self.y = np.full(self.population_size, np.nan)
154    
155    def forward(self, x):
156        self.e_summation = np.sum(np.exp(x), axis=1, keepdims=True)
157        return self.energy * softmax(x)
158    
159    def backward(self, x):
160        """
161        Softmax is not injective. This is a pseudo-inverse.
162        """
163        return np.log(self.e_summation * x/self.energy)
164
165    def view(self, x):
166        """
167        Maps the input from the domain to the codomain.
168        """
169        return self.forward(x)
170    
171    def view_g(self):
172        """
173        Maps the input from the domain to the codomain.
174        """
175        return self.forward(self.population)
176
177    def inverse_view(self, x):
178        """
179        Maps the input from the codomain to the domain.
180        """
181        return self.backward(x)
182    
183    def ask_oracle(self, X: np.ndarray) -> np.ndarray:
184        return super().ask_oracle()
185
186    def init_oracle(self):
187        return self.ask_oracle(self.view_g())
188
189    def step(self):
190        """
191        This function performs a single step of the Nelder-Mead algorithm.
192        """
193        
194        # Sort the population by the value of the oracle
195        y = self.ask_oracle(self.view_g())
196        idx = np.argsort(y)
197        self.population = self.population[idx]
198
199        # Compute the centroid of the population
200        centroid = self.population[:-1].mean(axis=0)
201
202        # Reflect the worst point through the centroid
203        reflected = centroid + (centroid - self.population[-1])
204
205        # Evaluate the reflected point
206        
207        y_reflected = self.ask_oracle(self.view(reflected.reshape(1,-1)))
208
209        # If the reflected point is better than the second worst, but not better than the best, then expand
210        
211        if y_reflected < y[-2] and y_reflected > y[0]:
212            expanded = centroid + (reflected - centroid)
213            y_expanded = self.ask_oracle(self.view(expanded.reshape(1,-1)))
214            if y_expanded < y_reflected:
215                self.population[-1] = expanded
216            else:
217                self.population[-1] = reflected
218        # If the reflected point is better than the best, then expand
219        elif y_reflected < y[0]:
220            expanded = centroid + 2 * (reflected - centroid)
221            y_expanded = self.ask_oracle(self.view(expanded.reshape(1,-1)))
222            if y_expanded < y_reflected:
223                self.population[-1] = expanded
224            else:
225                self.population[-1] = reflected
226        # If the reflected point is worse than the second worst, then contract
227        elif y_reflected > y[-2]:
228            contracted = centroid + 0.5 * (self.population[-1] - centroid)
229            y_contracted = self.ask_oracle(self.view(contracted.reshape(1,-1)))
230            if y_contracted < y[-1]:
231                self.population[-1] = contracted
232            else:
233                for i in range(1, len(self.population)):
234                    self.population[i] = 0.5 * (self.population[i] + self.population[0])
235        # If the reflected point is worse than the worst, then shrink
236        elif y_reflected > y[-1]:
237            for i in range(1, len(self.population)):
238                self.population[i] = 0.5 * (self.population[i] + self.population[0])
239
240        self.y = y.copy()
241        return self.view_g()

Abstract class for optimization algorithms

NelderMeadConstant( population_size: int, ranges: Union[int, Iterable], rng_seed: int = 0, energy=10) 
124    def __init__(self, population_size: int, ranges: Union[int, Iterable], rng_seed: int = 0, energy=10):
125        """
126        This Nelder-Mead algorithm assumes that the optimization function 'f' is to be minimized and time independent.
127
128        The parameters are constrained so that their sum always add up to `energy` though a softmax function.
129
130        Parameters:
131        - population_size (int): The size of the population.
132        - ranges (int or list): The number of parameters or a list of ranges for each parameter.
133        - rng_seed (int): The seed for the random number generator.
134        - energy (float): The energy parameter used in the optimization.
135
136        Note: The ranges parameter is merely a placeholder. The ranges are set to (0, energy) for all parameters.
137        """
138        self.population_size = population_size
139        if isinstance(ranges, Iterable):
140            self.ranges = np.array([(0, energy) for _ in range(len(ranges))])
141            warnings.warn("While using Nelder-MeadConstant, the ranges are set to (0, energy) for all parameters. The parameter `ranges` is merely a placeholder.")
142        elif isinstance(ranges, int):
143            self.ranges = np.array([(0, energy) for _ in range(ranges)])
144        self.rng_seed = np.random.default_rng(rng_seed)
145
146        # Initialize the population (random position and initial momentum)
147        self.population = self.rng_seed.random((self.population_size, len(self.ranges)))
148
149        # Derived attributes
150        self.energy = energy
151        self.e_summation = 1
152        # Initialize y as vector of nans
153        self.y = np.full(self.population_size, np.nan)

This Nelder-Mead algorithm assumes that the optimization function 'f' is to be minimized and time independent.

The parameters are constrained so that their sum always add up to energy though a softmax function.

Parameters:

  • population_size (int): The size of the population.
  • ranges (int or list): The number of parameters or a list of ranges for each parameter.
  • rng_seed (int): The seed for the random number generator.
  • energy (float): The energy parameter used in the optimization.

Note: The ranges parameter is merely a placeholder. The ranges are set to (0, energy) for all parameters.

population_size
rng_seed
population
energy
e_summation
y
def forward(self, x): 
155    def forward(self, x):
156        self.e_summation = np.sum(np.exp(x), axis=1, keepdims=True)
157        return self.energy * softmax(x)
def backward(self, x): 
159    def backward(self, x):
160        """
161        Softmax is not injective. This is a pseudo-inverse.
162        """
163        return np.log(self.e_summation * x/self.energy)

Softmax is not injective. This is a pseudo-inverse.

def view(self, x): 
165    def view(self, x):
166        """
167        Maps the input from the domain to the codomain.
168        """
169        return self.forward(x)

Maps the input from the domain to the codomain.

def view_g(self): 
171    def view_g(self):
172        """
173        Maps the input from the domain to the codomain.
174        """
175        return self.forward(self.population)

Maps the input from the domain to the codomain.

def inverse_view(self, x): 
177    def inverse_view(self, x):
178        """
179        Maps the input from the codomain to the domain.
180        """
181        return self.backward(x)

Maps the input from the codomain to the domain.

def ask_oracle(self, X: numpy.ndarray) -> numpy.ndarray: 
183    def ask_oracle(self, X: np.ndarray) -> np.ndarray:
184        return super().ask_oracle()
def init_oracle(self): 
186    def init_oracle(self):
187        return self.ask_oracle(self.view_g())
def step(self): 
189    def step(self):
190        """
191        This function performs a single step of the Nelder-Mead algorithm.
192        """
193        
194        # Sort the population by the value of the oracle
195        y = self.ask_oracle(self.view_g())
196        idx = np.argsort(y)
197        self.population = self.population[idx]
198
199        # Compute the centroid of the population
200        centroid = self.population[:-1].mean(axis=0)
201
202        # Reflect the worst point through the centroid
203        reflected = centroid + (centroid - self.population[-1])
204
205        # Evaluate the reflected point
206        
207        y_reflected = self.ask_oracle(self.view(reflected.reshape(1,-1)))
208
209        # If the reflected point is better than the second worst, but not better than the best, then expand
210        
211        if y_reflected < y[-2] and y_reflected > y[0]:
212            expanded = centroid + (reflected - centroid)
213            y_expanded = self.ask_oracle(self.view(expanded.reshape(1,-1)))
214            if y_expanded < y_reflected:
215                self.population[-1] = expanded
216            else:
217                self.population[-1] = reflected
218        # If the reflected point is better than the best, then expand
219        elif y_reflected < y[0]:
220            expanded = centroid + 2 * (reflected - centroid)
221            y_expanded = self.ask_oracle(self.view(expanded.reshape(1,-1)))
222            if y_expanded < y_reflected:
223                self.population[-1] = expanded
224            else:
225                self.population[-1] = reflected
226        # If the reflected point is worse than the second worst, then contract
227        elif y_reflected > y[-2]:
228            contracted = centroid + 0.5 * (self.population[-1] - centroid)
229            y_contracted = self.ask_oracle(self.view(contracted.reshape(1,-1)))
230            if y_contracted < y[-1]:
231                self.population[-1] = contracted
232            else:
233                for i in range(1, len(self.population)):
234                    self.population[i] = 0.5 * (self.population[i] + self.population[0])
235        # If the reflected point is worse than the worst, then shrink
236        elif y_reflected > y[-1]:
237            for i in range(1, len(self.population)):
238                self.population[i] = 0.5 * (self.population[i] + self.population[0])
239
240        self.y = y.copy()
241        return self.view_g()

This function performs a single step of the Nelder-Mead algorithm.