1#    This file is part of DEAP.
2#
3#    DEAP is free software: you can redistribute it and/or modify
4#    it under the terms of the GNU Lesser General Public License as
5#    published by the Free Software Foundation, either version 3 of
6#    the License, or (at your option) any later version.
7#
8#    DEAP is distributed in the hope that it will be useful,
9#    but WITHOUT ANY WARRANTY; without even the implied warranty of
10#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11#    GNU Lesser General Public License for more details.
12#
13#    You should have received a copy of the GNU Lesser General Public
14#    License along with DEAP. If not, see <http://www.gnu.org/licenses/>.
15
16#    Special thanks to Nikolaus Hansen for providing major part of
17#    this code. The CMA-ES algorithm is provided in many other languages
18#    and advanced versions at http://www.lri.fr/~hansen/cmaesintro.html.
19
20"""A module that provides support for the Covariance Matrix Adaptation
21Evolution Strategy.
22"""
23import copy
24from math import sqrt, log, exp
25import numpy
26
27import tools
28
29
30class Strategy(object):
31    """
32    A strategy that will keep track of the basic parameters of the CMA-ES
33    algorithm ([Hansen2001]_).
34
35    :param centroid: An iterable object that indicates where to start the
36                     evolution.
37    :param sigma: The initial standard deviation of the distribution.
38    :param parameter: One or more parameter to pass to the strategy as
39                      described in the following table, optional.
40
41    +----------------+---------------------------+----------------------------+
42    | Parameter      | Default                   | Details                    |
43    +================+===========================+============================+
44    | ``lambda_``    | ``int(4 + 3 * log(N))``   | Number of children to      |
45    |                |                           | produce at each generation,|
46    |                |                           | ``N`` is the individual's  |
47    |                |                           | size (integer).            |
48    +----------------+---------------------------+----------------------------+
49    | ``mu``         | ``int(lambda_ / 2)``      | The number of parents to   |
50    |                |                           | keep from the              |
51    |                |                           | lambda children (integer). |
52    +----------------+---------------------------+----------------------------+
53    | ``cmatrix``    | ``identity(N)``           | The initial covariance     |
54    |                |                           | matrix of the distribution |
55    |                |                           | that will be sampled.      |
56    +----------------+---------------------------+----------------------------+
57    | ``weights``    | ``"superlinear"``         | Decrease speed, can be     |
58    |                |                           | ``"superlinear"``,         |
59    |                |                           | ``"linear"`` or            |
60    |                |                           | ``"equal"``.               |
61    +----------------+---------------------------+----------------------------+
62    | ``cs``         | ``(mueff + 2) /           | Cumulation constant for    |
63    |                | (N + mueff + 3)``         | step-size.                 |
64    +----------------+---------------------------+----------------------------+
65    | ``damps``      | ``1 + 2 * max(0, sqrt((   | Damping for step-size.     |
66    |                | mueff - 1) / (N + 1)) - 1)|                            |
67    |                | + cs``                    |                            |
68    +----------------+---------------------------+----------------------------+
69    | ``ccum``       | ``4 / (N + 4)``           | Cumulation constant for    |
70    |                |                           | covariance matrix.         |
71    +----------------+---------------------------+----------------------------+
72    | ``ccov1``      | ``2 / ((N + 1.3)^2 +      | Learning rate for rank-one |
73    |                | mueff)``                  | update.                    |
74    +----------------+---------------------------+----------------------------+
75    | ``ccovmu``     | ``2 * (mueff - 2 + 1 /    | Learning rate for rank-mu  |
76    |                | mueff) / ((N + 2)^2 +     | update.                    |
77    |                | mueff)``                  |                            |
78    +----------------+---------------------------+----------------------------+
79
80    .. [Hansen2001] Hansen and Ostermeier, 2001. Completely Derandomized
81       Self-Adaptation in Evolution Strategies. *Evolutionary Computation*
82
83    """
84    def __init__(self, centroid, sigma, **kargs):
85        self.params = kargs
86
87        # Create a centroid as a numpy array
88        self.centroid = numpy.array(centroid)
89
90        self.dim = len(self.centroid)
91        self.sigma = sigma
92        self.pc = numpy.zeros(self.dim)
93        self.ps = numpy.zeros(self.dim)
94        self.chiN = sqrt(self.dim) * (1 - 1. / (4. * self.dim) +
95                                      1. / (21. * self.dim ** 2))
96
97        self.C = self.params.get("cmatrix", numpy.identity(self.dim))
98        self.diagD, self.B = numpy.linalg.eigh(self.C)
99
100        indx = numpy.argsort(self.diagD)
101        self.diagD = self.diagD[indx] ** 0.5
102        self.B = self.B[:, indx]
103        self.BD = self.B * self.diagD
104
105        self.cond = self.diagD[indx[-1]] / self.diagD[indx[0]]
106
107        self.lambda_ = self.params.get("lambda_", int(4 + 3 * log(self.dim)))
108        self.update_count = 0
109        self.computeParams(self.params)
110
111    def generate(self, ind_init):
112        """Generate a population of :math:`\lambda` individuals of type
113        *ind_init* from the current strategy.
114
115        :param ind_init: A function object that is able to initialize an
116                         individual from a list.
117        :returns: A list of individuals.
118        """
119        arz = numpy.random.standard_normal((self.lambda_, self.dim))
120        arz = self.centroid + self.sigma * numpy.dot(arz, self.BD.T)
121        return map(ind_init, arz)
122
123    def update(self, population):
124        """Update the current covariance matrix strategy from the
125        *population*.
126
127        :param population: A list of individuals from which to update the
128                           parameters.
129        """
130        population.sort(key=lambda ind: ind.fitness, reverse=True)
131
132        old_centroid = self.centroid
133        self.centroid = numpy.dot(self.weights, population[0:self.mu])
134
135        c_diff = self.centroid - old_centroid
136
137        # Cumulation : update evolution path
138        self.ps = (1 - self.cs) * self.ps \
139            + sqrt(self.cs * (2 - self.cs) * self.mueff) / self.sigma \
140            * numpy.dot(self.B, (1. / self.diagD) *
141                        numpy.dot(self.B.T, c_diff))
142
143        hsig = float((numpy.linalg.norm(self.ps) /
144                      sqrt(1. - (1. - self.cs) ** (2. * (self.update_count + 1.))) / self.chiN <
145                      (1.4 + 2. / (self.dim + 1.))))
146
147        self.update_count += 1
148
149        self.pc = (1 - self.cc) * self.pc + hsig \
150            * sqrt(self.cc * (2 - self.cc) * self.mueff) / self.sigma \
151            * c_diff
152
153        # Update covariance matrix
154        artmp = population[0:self.mu] - old_centroid
155        self.C = (1 - self.ccov1 - self.ccovmu + (1 - hsig) *
156                  self.ccov1 * self.cc * (2 - self.cc)) * self.C \
157            + self.ccov1 * numpy.outer(self.pc, self.pc) \
158            + self.ccovmu * numpy.dot((self.weights * artmp.T), artmp) \
159            / self.sigma ** 2
160
161        self.sigma *= numpy.exp((numpy.linalg.norm(self.ps) / self.chiN - 1.) *
162                                self.cs / self.damps)
163
164        self.diagD, self.B = numpy.linalg.eigh(self.C)
165        indx = numpy.argsort(self.diagD)
166
167        self.cond = self.diagD[indx[-1]] / self.diagD[indx[0]]
168
169        self.diagD = self.diagD[indx] ** 0.5
170        self.B = self.B[:, indx]
171        self.BD = self.B * self.diagD
172
173    def computeParams(self, params):
174        """Computes the parameters depending on :math:`\lambda`. It needs to
175        be called again if :math:`\lambda` changes during evolution.
176
177        :param params: A dictionary of the manually set parameters.
178        """
179        self.mu = params.get("mu", int(self.lambda_ / 2))
180        rweights = params.get("weights", "superlinear")
181        if rweights == "superlinear":
182            self.weights = log(self.mu + 0.5) - \
183                numpy.log(numpy.arange(1, self.mu + 1))
184        elif rweights == "linear":
185            self.weights = self.mu + 0.5 - numpy.arange(1, self.mu + 1)
186        elif rweights == "equal":
187            self.weights = numpy.ones(self.mu)
188        else:
189            raise RuntimeError("Unknown weights : %s" % rweights)
190
191        self.weights /= sum(self.weights)
192        self.mueff = 1. / sum(self.weights ** 2)
193
194        self.cc = params.get("ccum", 4. / (self.dim + 4.))
195        self.cs = params.get("cs", (self.mueff + 2.) /
196                             (self.dim + self.mueff + 3.))
197        self.ccov1 = params.get("ccov1", 2. / ((self.dim + 1.3) ** 2 +
198                                               self.mueff))
199        self.ccovmu = params.get("ccovmu", 2. * (self.mueff - 2. +
200                                                 1. / self.mueff) /
201                                 ((self.dim + 2.) ** 2 + self.mueff))
202        self.ccovmu = min(1 - self.ccov1, self.ccovmu)
203        self.damps = 1. + 2. * max(0, sqrt((self.mueff - 1.) /
204                                           (self.dim + 1.)) - 1.) + self.cs
205        self.damps = params.get("damps", self.damps)
206
207
208class StrategyOnePlusLambda(object):
209    """
210    A CMA-ES strategy that uses the :math:`1 + \lambda` paradigm ([Igel2007]_).
211
212    :param parent: An iterable object that indicates where to start the
213                   evolution. The parent requires a fitness attribute.
214    :param sigma: The initial standard deviation of the distribution.
215    :param lambda_: Number of offspring to produce from the parent.
216                    (optional, defaults to 1)
217    :param parameter: One or more parameter to pass to the strategy as
218                      described in the following table. (optional)
219
220    Other parameters can be provided as described in the next table
221
222    +----------------+---------------------------+----------------------------+
223    | Parameter      | Default                   | Details                    |
224    +================+===========================+============================+
225    | ``d``          | ``1.0 + N / (2.0 *        | Damping for step-size.     |
226    |                | lambda_)``                |                            |
227    +----------------+---------------------------+----------------------------+
228    | ``ptarg``      | ``1.0 / (5 + sqrt(lambda_)| Taget success rate.        |
229    |                | / 2.0)``                  |                            |
230    +----------------+---------------------------+----------------------------+
231    | ``cp``         | ``ptarg * lambda_ / (2.0 +| Step size learning rate.   |
232    |                | ptarg * lambda_)``        |                            |
233    +----------------+---------------------------+----------------------------+
234    | ``cc``         | ``2.0 / (N + 2.0)``       | Cumulation time horizon.   |
235    +----------------+---------------------------+----------------------------+
236    | ``ccov``       | ``2.0 / (N**2 + 6.0)``    | Covariance matrix learning |
237    |                |                           | rate.                      |
238    +----------------+---------------------------+----------------------------+
239    | ``pthresh``    | ``0.44``                  | Threshold success rate.    |
240    +----------------+---------------------------+----------------------------+
241
242    .. [Igel2007] Igel, Hansen, Roth, 2007. Covariance matrix adaptation for
243    multi-objective optimization. *Evolutionary Computation* Spring;15(1):1-28
244
245    """
246    def __init__(self, parent, sigma, **kargs):
247        self.parent = parent
248        self.sigma = sigma
249        self.dim = len(self.parent)
250
251        self.C = numpy.identity(self.dim)
252        self.A = numpy.identity(self.dim)
253
254        self.pc = numpy.zeros(self.dim)
255
256        self.computeParams(kargs)
257        self.psucc = self.ptarg
258
259    def computeParams(self, params):
260        """Computes the parameters depending on :math:`\lambda`. It needs to
261        be called again if :math:`\lambda` changes during evolution.
262
263        :param params: A dictionary of the manually set parameters.
264        """
265        # Selection :
266        self.lambda_ = params.get("lambda_", 1)
267
268        # Step size control :
269        self.d = params.get("d", 1.0 + self.dim / (2.0 * self.lambda_))
270        self.ptarg = params.get("ptarg", 1.0 / (5 + sqrt(self.lambda_) / 2.0))
271        self.cp = params.get("cp", self.ptarg * self.lambda_ / (2 + self.ptarg * self.lambda_))
272
273        # Covariance matrix adaptation
274        self.cc = params.get("cc", 2.0 / (self.dim + 2.0))
275        self.ccov = params.get("ccov", 2.0 / (self.dim ** 2 + 6.0))
276        self.pthresh = params.get("pthresh", 0.44)
277
278    def generate(self, ind_init):
279        """Generate a population of :math:`\lambda` individuals of type
280        *ind_init* from the current strategy.
281
282        :param ind_init: A function object that is able to initialize an
283                         individual from a list.
284        :returns: A list of individuals.
285        """
286        # self.y = numpy.dot(self.A, numpy.random.standard_normal(self.dim))
287        arz = numpy.random.standard_normal((self.lambda_, self.dim))
288        arz = self.parent + self.sigma * numpy.dot(arz, self.A.T)
289        return map(ind_init, arz)
290
291    def update(self, population):
292        """Update the current covariance matrix strategy from the
293        *population*.
294
295        :param population: A list of individuals from which to update the
296                           parameters.
297        """
298        population.sort(key=lambda ind: ind.fitness, reverse=True)
299        lambda_succ = sum(self.parent.fitness <= ind.fitness for ind in population)
300        p_succ = float(lambda_succ) / self.lambda_
301        self.psucc = (1 - self.cp) * self.psucc + self.cp * p_succ
302
303        if self.parent.fitness <= population[0].fitness:
304            x_step = (population[0] - numpy.array(self.parent)) / self.sigma
305            self.parent = copy.deepcopy(population[0])
306            if self.psucc < self.pthresh:
307                self.pc = (1 - self.cc) * self.pc + sqrt(self.cc * (2 - self.cc)) * x_step
308                self.C = (1 - self.ccov) * self.C + self.ccov * numpy.outer(self.pc, self.pc)
309            else:
310                self.pc = (1 - self.cc) * self.pc
311                self.C = (1 - self.ccov) * self.C + self.ccov * (numpy.outer(self.pc, self.pc) + self.cc * (2 - self.cc) * self.C)
312
313        self.sigma = self.sigma * exp(1.0 / self.d * (self.psucc - self.ptarg) / (1.0 - self.ptarg))
314
315        # We use Cholesky since for now we have no use of eigen decomposition
316        # Basically, Cholesky returns a matrix A as C = A*A.T
317        # Eigen decomposition returns two matrix B and D^2 as C = B*D^2*B.T = B*D*D*B.T
318        # So A == B*D
319        # To compute the new individual we need to multiply each vector z by A
320        # as y = centroid + sigma * A*z
321        # So the Cholesky is more straightforward as we don't need to compute
322        # the squareroot of D^2, and multiply B and D in order to get A, we directly get A.
323        # This can't be done (without cost) with the standard CMA-ES as the eigen decomposition is used
324        # to compute covariance matrix inverse in the step-size evolutionary path computation.
325        self.A = numpy.linalg.cholesky(self.C)
326
327
328class StrategyMultiObjective(object):
329    """Multiobjective CMA-ES strategy based on the paper [Voss2010]_. It
330    is used similarly as the standard CMA-ES strategy with a generate-update
331    scheme.
332
333    :param population: An initial population of individual.
334    :param sigma: The initial step size of the complete system.
335    :param mu: The number of parents to use in the evolution. When not
336               provided it defaults to the length of *population*. (optional)
337    :param lambda_: The number of offspring to produce at each generation.
338                    (optional, defaults to 1)
339    :param indicator: The indicator function to use. (optional, default to
340                      :func:`~deap.tools.hypervolume`)
341
342    Other parameters can be provided as described in the next table
343
344    +----------------+---------------------------+----------------------------+
345    | Parameter      | Default                   | Details                    |
346    +================+===========================+============================+
347    | ``d``          | ``1.0 + N / 2.0``         | Damping for step-size.     |
348    +----------------+---------------------------+----------------------------+
349    | ``ptarg``      | ``1.0 / (5 + 1.0 / 2.0)`` | Taget success rate.        |
350    +----------------+---------------------------+----------------------------+
351    | ``cp``         | ``ptarg / (2.0 + ptarg)`` | Step size learning rate.   |
352    +----------------+---------------------------+----------------------------+
353    | ``cc``         | ``2.0 / (N + 2.0)``       | Cumulation time horizon.   |
354    +----------------+---------------------------+----------------------------+
355    | ``ccov``       | ``2.0 / (N**2 + 6.0)``    | Covariance matrix learning |
356    |                |                           | rate.                      |
357    +----------------+---------------------------+----------------------------+
358    | ``pthresh``    | ``0.44``                  | Threshold success rate.    |
359    +----------------+---------------------------+----------------------------+
360
361    .. [Voss2010] Voss, Hansen, Igel, "Improved Step Size Adaptation
362       for the MO-CMA-ES", 2010.
363
364    """
365    def __init__(self, population, sigma, **params):
366        self.parents = population
367        self.dim = len(self.parents[0])
368
369        # Selection
370        self.mu = params.get("mu", len(self.parents))
371        self.lambda_ = params.get("lambda_", 1)
372
373        # Step size control
374        self.d = params.get("d", 1.0 + self.dim / 2.0)
375        self.ptarg = params.get("ptarg", 1.0 / (5.0 + 0.5))
376        self.cp = params.get("cp", self.ptarg / (2.0 + self.ptarg))
377
378        # Covariance matrix adaptation
379        self.cc = params.get("cc", 2.0 / (self.dim + 2.0))
380        self.ccov = params.get("ccov", 2.0 / (self.dim ** 2 + 6.0))
381        self.pthresh = params.get("pthresh", 0.44)
382
383        # Internal parameters associated to the mu parent
384        self.sigmas = [sigma] * len(population)
385        # Lower Cholesky matrix (Sampling matrix)
386        self.A = [numpy.identity(self.dim) for _ in range(len(population))]
387        # Inverse Cholesky matrix (Used in the update of A)
388        self.invCholesky = [numpy.identity(self.dim) for _ in range(len(population))]
389        self.pc = [numpy.zeros(self.dim) for _ in range(len(population))]
390        self.psucc = [self.ptarg] * len(population)
391
392        self.indicator = params.get("indicator", tools.hypervolume)
393
394    def generate(self, ind_init):
395        """Generate a population of :math:`\lambda` individuals of type
396        *ind_init* from the current strategy.
397
398        :param ind_init: A function object that is able to initialize an
399                         individual from a list.
400        :returns: A list of individuals with a private attribute :attr:`_ps`.
401                  This last attribute is essential to the update function, it
402                  indicates that the individual is an offspring and the index
403                  of its parent.
404        """
405        arz = numpy.random.randn(self.lambda_, self.dim)
406        individuals = list()
407
408        # Make sure every parent has a parent tag and index
409        for i, p in enumerate(self.parents):
410            p._ps = "p", i
411
412        # Each parent produce an offspring
413        if self.lambda_ == self.mu:
414            for i in range(self.lambda_):
415                # print "Z", list(arz[i])
416                individuals.append(ind_init(self.parents[i] + self.sigmas[i] * numpy.dot(self.A[i], arz[i])))
417                individuals[-1]._ps = "o", i
418
419        # Parents producing an offspring are chosen at random from the first front
420        else:
421            ndom = tools.sortLogNondominated(self.parents, len(self.parents), first_front_only=True)
422            for i in range(self.lambda_):
423                j = numpy.random.randint(0, len(ndom))
424                _, p_idx = ndom[j]._ps
425                individuals.append(ind_init(self.parents[p_idx] + self.sigmas[p_idx] * numpy.dot(self.A[p_idx], arz[i])))
426                individuals[-1]._ps = "o", p_idx
427
428        return individuals
429
430    def _select(self, candidates):
431        if len(candidates) <= self.mu:
432            return candidates, []
433
434        pareto_fronts = tools.sortLogNondominated(candidates, len(candidates))
435
436        chosen = list()
437        mid_front = None
438        not_chosen = list()
439
440        # Fill the next population (chosen) with the fronts until there is not enouch space
441        # When an entire front does not fit in the space left we rely on the hypervolume
442        # for this front
443        # The remaining fronts are explicitely not chosen
444        full = False
445        for front in pareto_fronts:
446            if len(chosen) + len(front) <= self.mu and not full:
447                chosen += front
448            elif mid_front is None and len(chosen) < self.mu:
449                mid_front = front
450                # With this front, we selected enough individuals
451                full = True
452            else:
453                not_chosen += front
454
455        # Separate the mid front to accept only k individuals
456        k = self.mu - len(chosen)
457        if k > 0:
458            # reference point is chosen in the complete population
459            # as the worst in each dimension +1
460            ref = numpy.array([ind.fitness.wvalues for ind in candidates]) * -1
461            ref = numpy.max(ref, axis=0) + 1
462
463            for _ in range(len(mid_front) - k):
464                idx = self.indicator(mid_front, ref=ref)
465                not_chosen.append(mid_front.pop(idx))
466
467            chosen += mid_front
468
469        return chosen, not_chosen
470
471    def _rankOneUpdate(self, invCholesky, A, alpha, beta, v):
472        w = numpy.dot(invCholesky, v)
473
474        # Under this threshold, the update is mostly noise
475        if w.max() > 1e-20:
476            w_inv = numpy.dot(w, invCholesky)
477            norm_w2 = numpy.sum(w ** 2)
478            a = sqrt(alpha)
479            root = numpy.sqrt(1 + beta / alpha * norm_w2)
480            b = a / norm_w2 * (root - 1)
481
482            A = a * A + b * numpy.outer(v, w)
483            invCholesky = 1.0 / a * invCholesky - b / (a ** 2 + a * b * norm_w2) * numpy.outer(w, w_inv)
484
485        return invCholesky, A
486
487    def update(self, population):
488        """Update the current covariance matrix strategies from the
489        *population*.
490
491        :param population: A list of individuals from which to update the
492                           parameters.
493        """
494        chosen, not_chosen = self._select(population + self.parents)
495
496        cp, cc, ccov = self.cp, self.cc, self.ccov
497        d, ptarg, pthresh = self.d, self.ptarg, self.pthresh
498
499        # Make copies for chosen offspring only
500        last_steps = [self.sigmas[ind._ps[1]] if ind._ps[0] == "o" else None for ind in chosen]
501        sigmas = [self.sigmas[ind._ps[1]] if ind._ps[0] == "o" else None for ind in chosen]
502        invCholesky = [self.invCholesky[ind._ps[1]].copy() if ind._ps[0] == "o" else None for ind in chosen]
503        A = [self.A[ind._ps[1]].copy() if ind._ps[0] == "o" else None for ind in chosen]
504        pc = [self.pc[ind._ps[1]].copy() if ind._ps[0] == "o" else None for ind in chosen]
505        psucc = [self.psucc[ind._ps[1]] if ind._ps[0] == "o" else None for ind in chosen]
506
507        # Update the internal parameters for successful offspring
508        for i, ind in enumerate(chosen):
509            t, p_idx = ind._ps
510
511            # Only the offspring update the parameter set
512            if t == "o":
513                # Update (Success = 1 since it is chosen)
514                psucc[i] = (1.0 - cp) * psucc[i] + cp
515                sigmas[i] = sigmas[i] * exp((psucc[i] - ptarg) / (d * (1.0 - ptarg)))
516
517                if psucc[i] < pthresh:
518                    xp = numpy.array(ind)
519                    x = numpy.array(self.parents[p_idx])
520                    pc[i] = (1.0 - cc) * pc[i] + sqrt(cc * (2.0 - cc)) * (xp - x) / last_steps[i]
521                    invCholesky[i], A[i] = self._rankOneUpdate(invCholesky[i], A[i], 1 - ccov, ccov, pc[i])
522                else:
523                    pc[i] = (1.0 - cc) * pc[i]
524                    pc_weight = cc * (2.0 - cc)
525                    invCholesky[i], A[i] = self._rankOneUpdate(invCholesky[i], A[i], 1 - ccov + pc_weight, ccov, pc[i])
526
527                self.psucc[p_idx] = (1.0 - cp) * self.psucc[p_idx] + cp
528                self.sigmas[p_idx] = self.sigmas[p_idx] * exp((self.psucc[p_idx] - ptarg) / (d * (1.0 - ptarg)))
529
530        # It is unnecessary to update the entire parameter set for not chosen individuals
531        # Their parameters will not make it to the next generation
532        for ind in not_chosen:
533            t, p_idx = ind._ps
534
535            # Only the offspring update the parameter set
536            if t == "o":
537                self.psucc[p_idx] = (1.0 - cp) * self.psucc[p_idx]
538                self.sigmas[p_idx] = self.sigmas[p_idx] * exp((self.psucc[p_idx] - ptarg) / (d * (1.0 - ptarg)))
539
540        # Make a copy of the internal parameters
541        # The parameter is in the temporary variable for offspring and in the original one for parents
542        self.parents = chosen
543        self.sigmas = [sigmas[i] if ind._ps[0] == "o" else self.sigmas[ind._ps[1]] for i, ind in enumerate(chosen)]
544        self.invCholesky = [invCholesky[i] if ind._ps[0] == "o" else self.invCholesky[ind._ps[1]] for i, ind in enumerate(chosen)]
545        self.A = [A[i] if ind._ps[0] == "o" else self.A[ind._ps[1]] for i, ind in enumerate(chosen)]
546        self.pc = [pc[i] if ind._ps[0] == "o" else self.pc[ind._ps[1]] for i, ind in enumerate(chosen)]
547        self.psucc = [psucc[i] if ind._ps[0] == "o" else self.psucc[ind._ps[1]] for i, ind in enumerate(chosen)]
548