-
Notifications
You must be signed in to change notification settings - Fork 2
/
cma.py
5981 lines (5181 loc) · 261 KB
/
cma.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#!/usr/bin/env python
"""Module cma implements the CMA-ES, Covariance Matrix Adaptation Evolution
Strategy, a stochastic optimizer for robust non-linear non-convex
derivative-free function minimization for Python versions 2.6 and 2.7
(for Python 2.5 class SolutionDict would need to be re-implemented, because
it depends on collections.MutableMapping, since version 0.91.01).
CMA-ES searches for a minimizer (a solution x in R**n) of an
objective function f (cost function), such that f(x) is
minimal. Regarding f, only function values for candidate solutions
need to be available, gradients are not necessary. Even less
restrictive, only a passably reliable ranking of the candidate
solutions in each iteration is necessary, the function values
itself do not matter. Some termination criteria however depend
on actual f-values.
Two interfaces are provided:
- function `fmin(func, x0, sigma0,...)`
runs a complete minimization
of the objective function func with CMA-ES.
- class `CMAEvolutionStrategy`
allows for minimization such that the
control of the iteration loop remains with the user.
Used packages:
- unavoidable: `numpy` (see `barecmaes2.py` if `numpy` is not
available),
- avoidable with small changes: `time`, `sys`
- optional: `matplotlib.pylab` (for `plot` etc., highly
recommended), `pprint` (pretty print), `pickle` (in class
`Sections`), `doctest`, `inspect`, `pygsl` (never by default)
Testing
-------
The code can be tested on a given system. Typing::
python cma.py --test --quiet
or in the Python shell ``ipython -pylab``::
run cma.py --test --quiet
runs ``doctest.testmod(cma)`` showing only exceptions (and not the
tests that fail due to small differences in the output) and should
run without complaints in about under two minutes. On some systems,
the pop up windows must be closed manually to continue and finish
the test.
Install
-------
The code can be installed by::
python cma.py --install
where the ``setup`` function from package ``distutils.core`` is used.
Example
-------
::
import cma
help(cma) # "this" help message, use cma? in ipython
help(cma.fmin)
help(cma.CMAEvolutionStrategy)
help(cma.Options)
cma.Options('tol') # display 'tolerance' termination options
cma.Options('verb') # display verbosity options
res = cma.fmin(cma.Fcts.tablet, 15 * [1], 1)
res[0] # best evaluated solution
res[5] # mean solution, presumably better with noise
:See: `fmin()`, `Options`, `CMAEvolutionStrategy`
:Author: Nikolaus Hansen, 2008-2012
:License: GPL 2 and 3
"""
from __future__ import division # future is >= 3.0, this code has been used with 2.6 & 2.7
from __future__ import with_statement # only necessary for python 2.5 and not in heavy use
# from __future__ import collections.MutableMapping # does not exist in future, otherwise 2.5 would work
# from __future__ import print_function # for cross-checking, available from python 2.6
__version__ = "0.91.02 $Revision: 3168 $"
# $Date: 2012-03-09 18:35:03 +0100 (Fri, 09 Mar 2012) $
# bash: svn propset svn:keywords 'Date Revision' cma.py
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, version 2 or 3.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# for testing:
# pyflakes cma.py # finds bugs by static analysis
# pychecker --limit 60 cma.py # also executes, gives 60 warnings (all checked)
# python cma.py -t -quiet # executes implemented tests based on doctest
# to create a html documentation file:
# pydoc -w cma # edit the header (remove local pointers)
# epydoc cma.py # comes close to javadoc but does not find the
# # links of function references etc
# doxygen needs @package cma as first line in the module docstring
# some things like class attributes are not interpreted correctly
# sphinx: doc style of doc.python.org, could not make it work
# TODO: make those options that are only used in fmin an error in init of CMA, but still Options() should
# work as input to CMA.
# TODO: add a default logger in CMAEvolutionStrategy, see fmin() and optimize() first
# tell() should probably not add data, but optimize() should handle even an after_iteration_handler.
# TODO: CMAEvolutionStrategy(ones(10), 1).optimize(cma.fcts.elli) # should work like fmin
# one problem: the data logger is not default and seemingly cannot be attached in one line
# TODO: check combination of boundary handling and transformation: penalty must be computed
# on gp.pheno(x_geno, bounds=None), but without bounds, check/remove usage of .geno everywhere
# TODO: check whether all new solutions are put into self.sent_solutions
# TODO: separate initialize==reset_state from __init__
# TODO: introduce Zpos == diffC which makes the code more consistent and the active update "exact"
# TODO: split tell into a variable transformation part and the "pure" functionality
# usecase: es.tell_geno(X, [func(es.pheno(x)) for x in X])
# genotypic repair is not part of tell_geno
# TODO: read settable "options" from a (properties) file, see myproperties.py
#
# typical parameters in scipy.optimize: disp, xtol, ftol, maxiter, maxfun, callback=None
# maxfev, diag (A sequency of N positive entries that serve as
# scale factors for the variables.)
# full_output -- non-zero to return all optional outputs.
# If xtol < 0.0, xtol is set to sqrt(machine_precision)
# 'infot -- a dictionary of optional outputs with the keys:
# 'nfev': the number of function calls...
#
# see eg fmin_powell
# typical returns
# x, f, dictionary d
# (xopt, {fopt, gopt, Hopt, func_calls, grad_calls, warnflag}, <allvecs>)
#
# TODO: keep best ten solutions
# TODO: implement constraints handling
# TODO: option full_output -- non-zero to return all optional outputs.
# TODO: extend function unitdoctest, or use unittest?
# TODO: implement equal-fitness termination, covered by stagnation?
# TODO: apply style guide: no capitalizations!?
# TODO: check and test dispdata()
# TODO: eigh(): thorough testing would not hurt
#
# TODO (later): implement readSignals from a file like properties file (to be called after tell())
import sys, time # not really essential
import collections, numpy as np # arange, cos, size, eye, inf, dot, floor, outer, zeros, linalg.eigh, sort, argsort, random, ones,...
from numpy import inf, array, dot, exp, log, sqrt, sum # to access the built-in sum fct: __builtins__.sum or del sum removes the imported sum and recovers the shadowed
try:
import matplotlib.pylab as pylab # also: use ipython -pylab
show = pylab.show
savefig = pylab.savefig # we would like to be able to use cma.savefig() etc
closefig = pylab.close
except:
pylab = None
print(' Could not import matplotlib.pylab, therefore ``cma.plot()`` etc. is not available')
def show():
pass
__docformat__ = "reStructuredText" # this hides some comments entirely?
sys.py3kwarning = True # TODO: out-comment from version 2.6
# why not package math?
# TODO: check scitools.easyviz and how big the adaptation would be
# changes:
# 12/07/21: convert value True for noisehandling into 1 making the output compatible
# 12/01/30: class Solution and more old stuff removed r3101
# 12/01/29: class Solution is depreciated, GenoPheno and SolutionDict do the job (v0.91.00, r3100)
# 12/01/06: CMA_eigenmethod option now takes a function (integer still works)
# 11/09/30: flat fitness termination checks also history length
# 11/09/30: elitist option (using method clip_or_fit_solutions)
# 11/09/xx: method clip_or_fit_solutions for check_points option for all sorts of
# injected or modified solutions and even reliable adaptive encoding
# 11/08/19: fixed: scaling and typical_x type clashes 1 vs array(1) vs ones(dim) vs dim * [1]
# 11/07/25: fixed: fmin wrote first and last line even with verb_log==0
# fixed: method settableOptionsList, also renamed to versatileOptions
# default seed depends on time now
# 11/07/xx (0.9.92): added: active CMA, selective mirrored sampling, noise/uncertainty handling
# fixed: output argument ordering in fmin, print now only used as function
# removed: parallel option in fmin
# 11/07/01: another try to get rid of the memory leak by replacing self.unrepaired = self[:]
# 11/07/01: major clean-up and reworking of abstract base classes and of the documentation,
# also the return value of fmin changed and attribute stop is now a method.
# 11/04/22: bug-fix: option fixed_variables in combination with scaling
# 11/04/21: stopdict is not a copy anymore
# 11/04/15: option fixed_variables implemented
# 11/03/23: bug-fix boundary update was computed even without boundaries
# 11/03/12: bug-fix of variable annotation in plots
# 11/02/05: work around a memory leak in numpy
# 11/02/05: plotting routines improved
# 10/10/17: cleaning up, now version 0.9.30
# 10/10/17: bug-fix: return values of fmin now use phenotyp (relevant
# if input scaling_of_variables is given)
# 08/10/01: option evalparallel introduced,
# bug-fix for scaling being a vector
# 08/09/26: option CMAseparable becomes CMA_diagonal
# 08/10/18: some names change, test functions go into a class
# 08/10/24: more refactorizing
# 10/03/09: upper bound exp(min(1,...)) for step-size control
# TODO: this would define the visible interface
# __all__ = ['fmin', 'CMAEvolutionStrategy', 'plot', ...]
#
# emptysets = ('', (), [], {}) # array([]) does not work but also np.size(.) == 0
# "x in emptysets" cannot be well replaced by "not x"
# which is also True for array([]) and None, but also for 0 and False, and False for NaN
use_sent_solutions = True # 5-30% CPU slower, particularly for large lambda, will be mandatory soon
#____________________________________________________________
#____________________________________________________________
#
def unitdoctest():
"""is used to describe test cases and might in future become helpful
as an experimental tutorial as well. The main testing feature at the
moment is by doctest with ``cma._test()`` or conveniently by
``python cma.py --test``. Unfortunately, depending on the
system, the results will slightly differ and many "failed" test cases
might be reported. This is prevented with the --quiet option.
A simple first overall test:
>>> import cma
>>> res = cma.fmin(cma.fcts.elli, 3*[1], 1, CMA_diagonal=2, seed=1, verb_time=0)
(3_w,7)-CMA-ES (mu_w=2.3,w_1=58%) in dimension 3 (seed=1)
Covariance matrix is diagonal for 2 iterations (1/ccov=7.0)
Iterat #Fevals function value axis ratio sigma minstd maxstd min:sec
1 7 1.453161670768570e+04 1.2e+00 1.08e+00 1e+00 1e+00
2 14 3.281197961927601e+04 1.3e+00 1.22e+00 1e+00 2e+00
3 21 1.082851071704020e+04 1.3e+00 1.24e+00 1e+00 2e+00
100 700 8.544042012075362e+00 1.4e+02 3.18e-01 1e-03 2e-01
200 1400 5.691152415221861e-12 1.0e+03 3.82e-05 1e-09 1e-06
220 1540 3.890107746209078e-15 9.5e+02 4.56e-06 8e-11 7e-08
termination on tolfun : 1e-11
final/bestever f-value = 3.89010774621e-15 2.52273602735e-15
mean solution: [ -4.63614606e-08 -3.42761465e-10 1.59957987e-11]
std deviation: [ 6.96066282e-08 2.28704425e-09 7.63875911e-11]
Test on the Rosenbrock function with 3 restarts. The first trial only
finds the local optimum, which happens in about 20% of the cases.
>>> import cma
>>> res = cma.fmin(cma.fcts.rosen, 4*[-1],1, ftarget=1e-6, restarts=3, verb_time=0, verb_disp=500, seed=3)
(4_w,8)-CMA-ES (mu_w=2.6,w_1=52%) in dimension 4 (seed=3)
Iterat #Fevals function value axis ratio sigma minstd maxstd min:sec
1 8 4.875315645656848e+01 1.0e+00 8.43e-01 8e-01 8e-01
2 16 1.662319948123120e+02 1.1e+00 7.67e-01 7e-01 8e-01
3 24 6.747063604799602e+01 1.2e+00 7.08e-01 6e-01 7e-01
184 1472 3.701428610430019e+00 4.3e+01 9.41e-07 3e-08 5e-08
termination on tolfun : 1e-11
final/bestever f-value = 3.70142861043 3.70142861043
mean solution: [-0.77565922 0.61309336 0.38206284 0.14597202]
std deviation: [ 2.54211502e-08 3.88803698e-08 4.74481641e-08 3.64398108e-08]
(8_w,16)-CMA-ES (mu_w=4.8,w_1=32%) in dimension 4 (seed=4)
Iterat #Fevals function value axis ratio sigma minstd maxstd min:sec
1 1489 2.011376859371495e+02 1.0e+00 8.90e-01 8e-01 9e-01
2 1505 4.157106647905128e+01 1.1e+00 8.02e-01 7e-01 7e-01
3 1521 3.548184889359060e+01 1.1e+00 1.02e+00 8e-01 1e+00
111 3249 6.831867555502181e-07 5.1e+01 2.62e-02 2e-04 2e-03
termination on ftarget : 1e-06
final/bestever f-value = 6.8318675555e-07 1.18576673231e-07
mean solution: [ 0.99997004 0.99993938 0.99984868 0.99969505]
std deviation: [ 0.00018973 0.00038006 0.00076479 0.00151402]
>>> assert res[1] <= 1e-6
Notice the different termination conditions. Termination on the target
function value ftarget prevents further restarts.
Test of scaling_of_variables option
>>> import cma
>>> opts = cma.Options()
>>> opts['seed'] = 456
>>> opts['verb_disp'] = 0
>>> opts['CMA_active'] = 1
>>> # rescaling of third variable: for searching in roughly
>>> # x0 plus/minus 1e3*sigma0 (instead of plus/minus sigma0)
>>> opts.scaling_of_variables = [1, 1, 1e3, 1]
>>> res = cma.fmin(cma.fcts.rosen, 4 * [0.1], 0.1, **opts)
termination on tolfun : 1e-11
final/bestever f-value = 2.68096173031e-14 1.09714829146e-14
mean solution: [ 1.00000001 1.00000002 1.00000004 1.00000007]
std deviation: [ 3.00466854e-08 5.88400826e-08 1.18482371e-07 2.34837383e-07]
The printed std deviations reflect the actual true value (not the one
in the internal representation which would be different).
:See: cma.main(), cma._test()
"""
pass
#____________________________________________________________
#____________________________________________________________
#
class BlancClass(object):
"""blanc container class for having a collection of attributes"""
#_____________________________________________________________________
#_____________________________________________________________________
#
class DerivedDictBase(collections.MutableMapping):
"""for conveniently adding features to a dictionary. The actual
dictionary is in ``self.data``. Copy-paste
and modify setitem, getitem, and delitem, if necessary"""
def __init__(self, *args, **kwargs):
# collections.MutableMapping.__init__(self)
super(DerivedDictBase, self).__init__()
# super(SolutionDict, self).__init__() # the same
self.data = dict(*args, **kwargs)
def __len__(self):
return len(self.data)
def __contains__(self, value):
return value in self.data
def __iter__(self):
return iter(self.data)
def __setitem__(self, key, value):
"""defines self[key] = value"""
self.data[key] = value
def __getitem__(self, key):
"""defines self[key]"""
return self.data[key]
def __delitem__(self, key):
del self.data[key]
class SolutionDict(DerivedDictBase):
"""dictionary with computation of an hash key for the inserted solutions and
a stack of previously inserted same solutions.
Each entry is meant to store additional information related to the solution.
>>> import cma, numpy as np
>>> d = cma.SolutionDict()
>>> x = np.array([1,2,4])
>>> d[x] = {'x': x, 'iteration': 1}
>>> d.get(x) == (d[x] if d.key(x) in d.keys() else None)
The last line is always true.
TODO: data_with_same_key behaves like a stack (see setitem and delitem), but rather should behave like a queue?!
A queue is less consistent with the operation self[key] = ..., if self.data_with_same_key[key] is not empty.
"""
def __init__(self, *args, **kwargs):
DerivedDictBase.__init__(self, *args, **kwargs)
self.data_with_same_key = {}
def key(self, x):
try:
return tuple(x)
except TypeError:
return x
def __setitem__(self, key, value):
"""defines self[key] = value"""
key = self.key(key)
if key in self.data_with_same_key:
self.data_with_same_key[key] += [self.data[key]]
elif key in self.data:
self.data_with_same_key[key] = [self.data[key]]
self.data[key] = value
def __getitem__(self, key):
"""defines self[key]"""
return self.data[self.key(key)]
def __delitem__(self, key):
"""remove only most current key-entry"""
key = self.key(key)
if key in self.data_with_same_key:
if len(self.data_with_same_key[key]) == 1:
self.data[key] = self.data_with_same_key.pop(key)[0]
else:
self.data[key] = self.data_with_same_key[key].pop(-1)
else:
del self.data[key]
def truncate(self, max_len, min_iter):
if len(self) > max_len:
for k in self.keys():
if self[k]['iteration'] < min_iter:
del self[k] # only deletes one item with k as key, should delete all?
class SolutionDictOld(dict):
"""depreciated, SolutionDict should do, to be removed after SolutionDict
has been successfully applied.
dictionary with computation of an hash key for the inserted solutions and
stack of previously inserted same solutions.
Each entry is meant to store additional information related to the solution.
Methods ``pop`` and ``get`` are modified accordingly.
d = SolutionDict()
x = array([1,2,4])
d.insert(x, {'x': x, 'iteration': 1})
d.get(x) == d[d.key(x)] if d.key(x) in d.keys() else d.get(x) is None
TODO: not yet tested
TODO: behaves like a stack (see _pop_derived), but rather should behave like a queue?!
A queue is less consistent with the operation self[key] = ..., if self.more[key] is not empty.
"""
def __init__(self):
self.more = {} # previously inserted same solutions
self._pop_base = self.pop
self.pop = self._pop_derived
self._get_base = self.get
self.get = self._get_derived
def key(self, x):
"""compute the hash key of ``x``"""
return tuple(x)
def insert(self, x, datadict):
key = self.key(x)
if key in self.more:
self.more[key] += [self[key]]
elif key in self:
self.more[key] = [self[key]]
self[key] = datadict
def _get_derived(self, x, default=None):
return self._get_base(self.key(x), default)
def _pop_derived(self, x):
key = self.key(x)
res = self[key]
if key in self.more:
if len(self.more[key]) == 1:
self[key] = self.more.pop(key)[0]
else:
self[key] = self.more[key].pop(-1)
return res
class BestSolution(object):
"""container to keep track of the best solution seen"""
def __init__(self, x=None, f=np.inf, evals=None):
"""initialize the best solution with `x`, `f`, and `evals`.
Better solutions have smaller `f`-values.
"""
self.x = x
self.x_geno = None
self.f = f if f is not None and f is not np.nan else np.inf
self.evals = evals
self.evalsall = evals
self.last = BlancClass()
self.last.x = x
self.last.f = f
def update(self, arx, xarchive=None, arf=None, evals=None):
"""checks for better solutions in list `arx`, based on the smallest
corresponding value in `arf`, alternatively, `update` may be called
with a `BestSolution` instance like ``update(another_best_solution)``
in which case the better solution becomes the current best.
`xarchive` is used to retrieve the genotype of a solution.
"""
if arf is not None: # find failsave minimum
minidx = np.nanargmin(arf)
if minidx is np.nan:
return
minarf = arf[minidx]
# minarf = reduce(lambda x, y: y if y and y is not np.nan and y < x else x, arf, np.inf)
if type(arx) == BestSolution:
self.evalsall = max((self.evalsall, arx.evalsall))
if arx.f is not None and arx.f < np.inf:
self.update([arx.x], xarchive, [arx.f], arx.evals)
return self
elif minarf < np.inf and (minarf < self.f or self.f is None):
self.x, self.f = arx[minidx], arf[minidx]
self.x_geno = xarchive[self.x]['geno'] if xarchive is not None else None
self.evals = None if not evals else evals - len(arf) + minidx+1
self.evalsall = evals
elif evals:
self.evalsall = evals
self.last.x = arx[minidx]
self.last.f = minarf
def get(self):
"""return ``(x, f, evals)`` """
return self.x, self.f, self.evals, self.x_geno
#____________________________________________________________
#____________________________________________________________
#
class BoundPenalty(object):
"""Computes the boundary penalty. Must be updated each iteration,
using the `update` method.
Details
-------
The penalty computes like ``sum(w[i] * (x[i]-xfeas[i])**2)``,
where `xfeas` is the closest feasible (in-bounds) solution from `x`.
The weight `w[i]` should be updated during each iteration using
the update method.
This class uses `GenoPheno.into_bounds` in method `update` to access
domain boundary values and repair. This inconsistency might be
removed in future.
"""
def __init__(self, bounds=None):
"""Argument bounds can be `None` or ``bounds[0]`` and ``bounds[1]``
are lower and upper domain boundaries, each is either `None` or
a scalar or a list or array of appropriate size.
"""
##
# bounds attribute reminds the domain boundary values
self.bounds = bounds
self.gamma = 1 # a very crude assumption
self.weights_initialized = False # gamma becomes a vector after initialization
self.hist = [] # delta-f history
def has_bounds(self):
"""return True, if any variable is bounded"""
bounds = self.bounds
if bounds in (None, [None, None]):
return False
for i in xrange(bounds[0]):
if bounds[0][i] is not None and bounds[0][i] > -np.inf:
return True
for i in xrange(bounds[1]):
if bounds[1][i] is not None and bounds[1][i] < np.inf:
return True
return False
def repair(self, x, bounds=None, copy=False, copy_always=False):
"""sets out-of-bounds components of ``x`` on the bounds.
Arguments
---------
`bounds`
can be `None`, in which case the "default" bounds are used,
or ``[lb, ub]``, where `lb` and `ub`
represent lower and upper domain bounds respectively that
can be `None` or a scalar or a list or array of length ``len(self)``
code is more or less copy-paste from Solution.repair, but never tested
"""
# TODO (old data): CPU(N,lam,iter=20,200,100): 3.3s of 8s for two bounds, 1.8s of 6.5s for one bound
# TODO: test whether np.max([bounds[0], x], axis=0) etc is speed relevant
if bounds is None:
bounds = self.bounds
if copy_always:
x_out = array(x, copy=True)
if bounds not in (None, [None, None], (None, None)): # solely for effiency
x_out = array(x, copy=True) if copy and not copy_always else x
if bounds[0] is not None:
if np.isscalar(bounds[0]):
for i in xrange(len(x)):
x_out[i] = max([bounds[0], x[i]])
else:
for i in xrange(len(x)):
if bounds[0][i] is not None:
x_out[i] = max([bounds[0][i], x[i]])
if bounds[1] is not None:
if np.isscalar(bounds[1]):
for i in xrange(len(x)):
x_out[i] = min([bounds[1], x[i]])
else:
for i in xrange(len(x)):
if bounds[1][i] is not None:
x_out[i] = min([bounds[1][i], x[i]])
return x_out # convenience return
#____________________________________________________________
#
def __call__(self, x, archive, gp):
"""returns the boundary violation penalty for `x` ,where `x` is a
single solution or a list or array of solutions.
If `bounds` is not `None`, the values in `bounds` are used, see `__init__`"""
if x in (None, (), []):
return x
if gp.bounds in (None, [None, None], (None, None)):
return 0.0 if np.isscalar(x[0]) else [0.0] * len(x) # no penalty
x_is_single_vector = np.isscalar(x[0])
x = [x] if x_is_single_vector else x
pen = []
for xi in x:
# CAVE: this does not work with already repaired values!!
# CPU(N,lam,iter=20,200,100)?: 3s of 10s, array(xi): 1s (check again)
# remark: one deep copy can be prevented by xold = xi first
xpheno = gp.pheno(archive[xi]['geno'])
xinbounds = gp.into_bounds(xpheno)
fac = 1 # exp(0.1 * (log(self.scal) - np.mean(self.scal)))
pen.append(sum(self.gamma * ((xinbounds - xpheno) / fac)**2) / len(xi))
return pen[0] if x_is_single_vector else pen
#____________________________________________________________
#
def feasible_ratio(self, solutions):
"""counts for each coordinate the number of feasible values in
``solutions`` and returns an array of length ``len(solutions[0])``
with the ratios.
`solutions` is a list or array of repaired `Solution` instances
"""
count = np.zeros(len(solutions[0]))
for x in solutions:
count += x.unrepaired == x
return count / float(len(solutions))
#____________________________________________________________
#
def update(self, function_values, es, bounds=None):
"""updates the weights for computing a boundary penalty.
Arguments
---------
`function_values`
all function values of recent population of solutions
`es`
`CMAEvolutionStrategy` object instance, in particular the
method `into_bounds` of the attribute `gp` of type `GenoPheno`
is used.
`bounds`
not (yet) in use other than for ``bounds == [None, None]`` nothing
is updated.
Reference: Hansen et al 2009, A Method for Handling Uncertainty...
IEEE TEC, with addendum at http://www.lri.fr/~hansen/TEC2009online.pdf
"""
if bounds is None:
bounds = self.bounds
if bounds is None or (bounds[0] is None and bounds[1] is None): # no bounds ==> no penalty
return self # len(function_values) * [0.0] # case without voilations
N = es.N
### prepare
# compute varis = sigma**2 * C_ii
varis = es.sigma**2 * array(N * [es.C] if np.isscalar(es.C) else ( # scalar case
es.C if np.isscalar(es.C[0]) else # diagonal matrix case
[es.C[i][i] for i in xrange(N)])) # full matrix case
dmean = (es.mean - es.gp.into_bounds(es.mean)) / varis**0.5
### Store/update a history of delta fitness value
fvals = sorted(function_values)
l = 1 + len(fvals)
val = fvals[3*l // 4] - fvals[l // 4] # exact interquartile range apart interpolation
val = val / np.mean(varis) # new: val is normalized with sigma of the same iteration
# insert val in history
if np.isfinite(val) and val > 0:
self.hist.insert(0, val)
elif val == inf and len(self.hist) > 1:
self.hist.insert(0, max(self.hist))
else:
pass # ignore 0 or nan values
if len(self.hist) > 20 + (3*N) / es.popsize:
self.hist.pop()
### prepare
dfit = np.median(self.hist) # median interquartile range
damp = min(1, es.sp.mueff/10./N)
### set/update weights
# Throw initialization error
if len(self.hist) == 0:
raise _Error('wrongful initialization, no feasible solution sampled. ' +
'Reasons can be mistakenly set bounds (lower bound not smaller than upper bound) or a too large initial sigma0 or... ' +
'See description of argument func in help(cma.fmin) or an example handling infeasible solutions in help(cma.CMAEvolutionStrategy). ')
# initialize weights
if (dmean.any() and (not self.weights_initialized or es.countiter == 2)): # TODO
self.gamma = array(N * [2*dfit])
self.weights_initialized = True
# update weights gamma
if self.weights_initialized:
edist = array(abs(dmean) - 3 * max(1, N**0.5/es.sp.mueff))
if 1 < 3: # this is better, around a factor of two
# increase single weights possibly with a faster rate than they can decrease
# value unit of edst is std dev, 3==random walk of 9 steps
self.gamma *= exp((edist>0) * np.tanh(edist/3) / 2.)**damp
# decrease all weights up to the same level to avoid single extremely small weights
# use a constant factor for pseudo-keeping invariance
self.gamma[self.gamma > 5 * dfit] *= exp(-1./3)**damp
# self.gamma[idx] *= exp(5*dfit/self.gamma[idx] - 1)**(damp/3)
elif 1 < 3 and (edist>0).any(): # previous method
# CAVE: min was max in TEC 2009
self.gamma[edist>0] *= 1.1**min(1, es.sp.mueff/10./N)
# max fails on cigtab(N=12,bounds=[0.1,None]):
# self.gamma[edist>0] *= 1.1**max(1, es.sp.mueff/10./N) # this was a bug!?
# self.gamma *= exp((edist>0) * np.tanh(edist))**min(1, es.sp.mueff/10./N)
else: # alternative version, but not better
solutions = es.pop # this has not been checked
r = self.feasible_ratio(solutions) # has to be the averaged over N iterations
self.gamma *= exp(np.max([N*[0], 0.3 - r], axis=0))**min(1, es.sp.mueff/10/N)
es.more_to_write = self.gamma if self.weights_initialized else np.ones(N)
### return penalty
# es.more_to_write = self.gamma if not np.isscalar(self.gamma) else N*[1]
return self # bound penalty values
#____________________________________________________________
#____________________________________________________________
#
class GenoPhenoBase(object):
"""depreciated, abstract base class for genotyp-phenotype transformation,
to be implemented.
See (and rather use) option ``transformation`` of ``fmin`` or ``CMAEvolutionStrategy``.
Example
-------
::
import cma
class Mygpt(cma.GenoPhenoBase):
def pheno(self, x):
return x # identity for the time being
gpt = Mygpt()
optim = cma.CMAEvolutionStrategy(...)
while not optim.stop():
X = optim.ask()
f = [func(gpt.pheno(x)) for x in X]
optim.tell(X, f)
In case of a repair, we might pass the repaired solution into `tell()`
(with check_points being True).
TODO: check usecases in `CMAEvolutionStrategy` and implement option GenoPhenoBase
"""
def pheno(self, x):
raise NotImplementedError()
return x
#____________________________________________________________
#____________________________________________________________
#
class GenoPheno(object):
"""Genotype-phenotype transformation.
Method `pheno` provides the transformation from geno- to phenotype,
that is from the internal representation to the representation used
in the objective function. Method `geno` provides the "inverse" pheno-
to genotype transformation. The geno-phenotype transformation comprises,
in this order:
- insert fixed variables (with the phenotypic and therefore quite
possibly "wrong" values)
- affine linear transformation (scaling and shift)
- user-defined transformation
- projection into feasible domain (boundaries)
- assign fixed variables their original phenotypic value
By default all transformations are the identity. The boundary
transformation is only applied, if the boundaries are given as argument to
the method `pheno` or `geno` respectively.
``geno`` is not really necessary and might disappear in future.
"""
def __init__(self, dim, scaling=None, typical_x=None, bounds=None, fixed_values=None, tf=None):
"""return `GenoPheno` instance with fixed dimension `dim`.
Keyword Arguments
-----------------
`scaling`
the diagonal of a scaling transformation matrix, multipliers
in the genotyp-phenotyp transformation, see `typical_x`
`typical_x`
``pheno = scaling*geno + typical_x``
`bounds` (obsolete, might disappear)
list with two elements,
lower and upper bounds both can be a scalar or a "vector"
of length dim or `None`. Without effect, as `bounds` must
be given as argument to `pheno()`.
`fixed_values`
a dictionary of variable indices and values, like ``{0:2.0, 2:1.1}``,
that a not subject to change, negative indices are dropped
(they act like incommenting the index), values are phenotypic
values.
`tf`
list of two user-defined transformation functions, or `None`.
``tf[0]`` is a function that transforms the internal representation
as used by the optimizer into a solution as used by the
objective function. ``tf[1]`` does the back-transformation.
For example ::
tf_0 = lambda x: [xi**2 for xi in x]
tf_1 = lambda x: [abs(xi)**0.5 fox xi in x]
or "equivalently" without the `lambda` construct ::
def tf_0(x):
return [xi**2 for xi in x]
def tf_1(x):
return [abs(xi)**0.5 fox xi in x]
``tf=[tf_0, tf_1]`` is a reasonable way to guaranty that only positive
values are used in the objective function.
Details
-------
If ``tf_1`` is ommitted, the initial x-value must be given as genotype (as the
phenotype-genotype transformation is unknown) and injection of solutions
might lead to unexpected results.
"""
self.N = dim
self.bounds = bounds
self.fixed_values = fixed_values
if tf is not None:
self.tf_pheno = tf[0]
self.tf_geno = tf[1] # TODO: should not necessarily be needed
# r = np.random.randn(dim)
# assert all(tf[0](tf[1](r)) - r < 1e-7)
# r = np.random.randn(dim)
# assert all(tf[0](tf[1](r)) - r > -1e-7)
print("WARNING in class GenoPheno: user defined transformations have not been tested thoroughly")
else:
self.tf_geno = None
self.tf_pheno = None
if fixed_values:
if type(fixed_values) is not dict:
raise _Error("fixed_values must be a dictionary {index:value,...}")
if max(fixed_values.keys()) >= dim:
raise _Error("max(fixed_values.keys()) = " + str(max(fixed_values.keys())) +
" >= dim=N=" + str(dim) + " is not a feasible index")
# convenience commenting functionality: drop negative keys
for k in fixed_values.keys():
if k < 0:
fixed_values.pop(k)
if bounds:
if len(bounds) != 2:
raise _Error('len(bounds) must be 2 for lower and upper bounds')
for i in (0,1):
if bounds[i] is not None:
bounds[i] = array(dim * [bounds[i]] if np.isscalar(bounds[i]) else
[b for b in bounds[i]])
def vec_is_default(vec, default_val=0):
"""None or [None] are also recognized as default"""
try:
if len(vec) == 1:
vec = vec[0] # [None] becomes None and is always default
else:
return False
except TypeError:
pass # vec is a scalar
if vec is None or vec == array(None) or vec == default_val:
return True
return False
self.scales = array(scaling)
if vec_is_default(self.scales, 1):
self.scales = 1 # CAVE: 1 is not array(1)
elif self.scales.shape is not () and len(self.scales) != self.N:
raise _Error('len(scales) == ' + str(len(self.scales)) +
' does not match dimension N == ' + str(self.N))
self.typical_x = array(typical_x)
if vec_is_default(self.typical_x, 0):
self.typical_x = 0
elif self.typical_x.shape is not () and len(self.typical_x) != self.N:
raise _Error('len(typical_x) == ' + str(len(self.typical_x)) +
' does not match dimension N == ' + str(self.N))
if (self.scales is 1 and
self.typical_x is 0 and
self.bounds in (None, [None, None]) and
self.fixed_values is None and
self.tf_pheno is None):
self.isidentity = True
else:
self.isidentity = False
def into_bounds(self, y, bounds=None, copy_never=False, copy_always=False):
"""Argument `y` is a phenotypic vector,
return `y` put into boundaries, as a copy iff ``y != into_bounds(y)``.
Note: this code is duplicated in `Solution.repair` and might
disappear in future.
"""
bounds = bounds if bounds is not None else self.bounds
if bounds in (None, [None, None]):
return y if not copy_always else array(y, copy=True)
if bounds[0] is not None:
if len(bounds[0]) not in (1, len(y)):
raise ValueError('len(bounds[0]) = ' + str(len(bounds[0])) +
' and len of initial solution (' + str(len(y)) + ') disagree')
if copy_never: # is rather slower
for i in xrange(len(y)):
y[i] = max(bounds[0][i], y[i])
else:
y = np.max([bounds[0], y], axis=0)
if bounds[1] is not None:
if len(bounds[1]) not in (1, len(y)):
raise ValueError('len(bounds[1]) = ' + str(len(bounds[1])) +
' and initial solution (' + str(len(y)) + ') disagree')
if copy_never:
for i in xrange(len(y)):
y[i] = min(bounds[1][i], y[i])
else:
y = np.min([bounds[1], y], axis=0)
return y
def pheno(self, x, bounds=None, copy=True, copy_always=False):
"""maps the genotypic input argument into the phenotypic space,
boundaries are only applied if argument ``bounds is not None``, see
help for class `GenoPheno`
"""
if copy_always and not copy:
raise ValueError('arguments copy_always=' + str(copy_always) +
' and copy=' + str(copy) + ' have inconsistent values')
if self.isidentity and bounds in (None, [None, None], (None, None)):
return x if not copy_always else array(x, copy=copy_always)
if self.fixed_values is None:
y = array(x, copy=copy) # make a copy, in case
else: # expand with fixed values
y = list(x) # is a copy
for i in sorted(self.fixed_values.keys()):
y.insert(i, self.fixed_values[i])
y = array(y, copy=False)
if self.scales is not 1: # just for efficiency
y *= self.scales
if self.typical_x is not 0:
y += self.typical_x
if self.tf_pheno is not None:
y = array(self.tf_pheno(y), copy=False)
if bounds is not None:
y = self.into_bounds(y, bounds)
if self.fixed_values is not None:
for i, k in self.fixed_values.items():
y[i] = k
return y
def geno(self, y, bounds=None, copy=True, copy_always=False, archive=None):
"""maps the phenotypic input argument into the genotypic space.
If `bounds` are given, first `y` is projected into the feasible
domain. In this case ``copy==False`` leads to a copy.
by default a copy is made only to prevent to modify ``y``
method geno is only needed if external solutions are injected
(geno(initial_solution) is depreciated and will disappear)
TODO: arg copy=True should become copy_never=False
"""
if archive is not None and bounds is not None:
try:
return archive[y]['geno']
except:
pass
x = array(y, copy=(copy and not self.isidentity) or copy_always)
# bounds = self.bounds if bounds is None else bounds
if bounds is not None: # map phenotyp into bounds first
x = self.into_bounds(x, bounds)
if self.isidentity:
return x
# user-defined transformation
if self.tf_geno is not None:
x = array(self.tf_geno(x), copy=False)
# affine-linear transformation: shift and scaling
if self.typical_x is not 0:
x -= self.typical_x
if self.scales is not 1: # just for efficiency
x /= self.scales
# kick out fixed_values
if self.fixed_values is not None:
# keeping the transformed values does not help much
# therefore it is omitted
if 1 < 3:
keys = sorted(self.fixed_values.keys())
x = array([x[i] for i in range(len(x)) if i not in keys], copy=False)
else: # TODO: is this more efficient?
x = list(x)
for key in sorted(self.fixed_values.keys(), reverse=True):