Xcos tests: fix blocks_set.tst after 59f9727d
[scilab.git] / scilab / modules / optimization / macros / genetic / optim_nsga2.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2008 - Yann COLLETTE <yann.collette@renault.com>
3 //
4 // Copyright (C) 2012 - 2016 - Scilab Enterprises
5 //
6 // This file is hereby licensed under the terms of the GNU GPL v2.0,
7 // pursuant to article 5.3.4 of the CeCILL v.2.1.
8 // This file was originally licensed under the terms of the CeCILL v2.1,
9 // and continues to be available under such terms.
10 // For more information, see the COPYING file which you should have received
11 // along with this program.
12
13 function [pop_opt, fobj_pop_opt, pop_init, fobj_pop_init] = optim_nsga2(ga_f, pop_size, nb_generation, p_mut, p_cross, Log, param)
14
15     [nargout, nargin] = argn();
16
17     if ~isdef("param", "local") then
18         param = [];
19     end
20
21     [codage_func, err] = get_param(param, "codage_func", coding_ga_identity);
22     [init_func, err] = get_param(param, "init_func", init_ga_default);
23     [crossover_func, err] = get_param(param, "crossover_func", crossover_ga_default);
24     [mutation_func, err] = get_param(param, "mutation_func", mutation_ga_default);
25     [nb_couples, err] = get_param(param, "nb_couples", 100);
26     [output_func, err] = get_param(param, 'output_func', output_nsga2_default);
27
28     if ~isdef("ga_f", "local") then
29         error(gettext("optim_nsga2: ga_f is mandatory"));
30     else
31         if typeof(ga_f)=="list" then
32             deff("y=_ga_f(x)", "y=ga_f(1)(x, ga_f(2:$))");
33         else
34             deff("y=_ga_f(x)", "y=ga_f(x)");
35         end
36     end
37
38     if ~isdef("pop_size", "local") then
39         pop_size = 100;
40     end
41     if ~isdef("nb_generation", "local") then
42         nb_generation = 10;
43     end
44     if ~isdef("p_mut", "local") then
45         p_mut = 0.1;
46     end
47     if ~isdef("p_cross", "local") then
48         p_cross = 0.1;
49     end
50     if ~isdef("Log", "local") then
51         Log = %F;
52     end
53
54     // Initialization of the population
55     if (Log) then
56         printf(gettext("%s: Initialization of the population\n"), "optim_nsga2");
57     end
58
59     Pop = init_func(pop_size, param);
60
61     if (nargout>=3) then
62         pop_init = Pop;
63     end
64
65     // Code the individuals
66     Pop = codage_func(Pop, "code", param);
67
68     for i=1:length(Pop)
69         FObj_Pop(i, :) = _ga_f(Pop(i));
70     end
71
72     // Compute the domination rank
73     Rank=DominationRank(FObj_Pop);
74
75     // Compute the crowding distance
76     MO_FObj_Pop = FObj_Pop;
77     Index    = 1:size(MO_FObj_Pop, 1);
78     Crowdist = zeros(size(MO_FObj_Pop, 1), 1);
79     for i=1:size(FObj_Pop, 2)
80         [tmp, Index_List] = gsort(MO_FObj_Pop(:, i));
81         MO_FObj_Pop       = MO_FObj_Pop(Index_List, :);
82         Index             = Index(Index_List);
83         Crowdist(Index_List(1)) = %inf;
84         Crowdist(Index_List($)) = %inf;
85         _Max = max(MO_FObj_Pop(:, i));
86         _Min = min(MO_FObj_Pop(:, i));
87         for j=2:size(MO_FObj_Pop, 1)-1
88             Crowdist(Index(j)) = Crowdist(Index(j)) - (MO_FObj_Pop(j+1, i) - MO_FObj_Pop(j-1, i)) / (_Max - _Min);
89         end
90     end
91
92     if (nargout==4) then
93         fobj_pop_init = FObj_Pop;
94     end
95
96     // The genetic algorithm
97     for It=1:nb_generation
98         //
99         // Selection
100         //
101         Indiv1 = list();
102         Indiv2 = list();
103         for j=1:nb_couples
104             // Selection of 2 individuals via binary tournament selection to fill Indiv1
105             Index1 = ceil((size(FObj_Pop, 1) - 1)*grand(1, 1, "def")+1);
106             Index2 = ceil((size(FObj_Pop, 1) - 1)*grand(1, 1, "def")+1);
107             if (Rank(Index1)<Rank(Index2)) | ((Rank(Index1)==Rank(Index2)) & (Crowdist(Index1)>Crowdist(Index2))) then
108                 Indiv1(j)        = Pop(Index1);
109                 FObj_Indiv1(j, :) = MO_FObj_Pop(Index1, :);
110             else
111                 Indiv1(j)        = Pop(Index2);
112                 FObj_Indiv1(j, :) = MO_FObj_Pop(Index2, :);
113             end
114             // Selection of 2 individuals via binary tournament selection to fill Indiv2
115             Index1 = ceil((size(FObj_Pop, 1) - 1)*grand(1, 1, "def")+1);
116             Index2 = ceil((size(FObj_Pop, 1) - 1)*grand(1, 1, "def")+1);
117             if (Rank(Index1)<Rank(Index2)) | ((Rank(Index1)==Rank(Index2)) & (Crowdist(Index1)>Crowdist(Index2))) then
118                 Indiv2(j)        = Pop(Index1);
119                 FObj_Indiv2(j, :) = MO_FObj_Pop(Index1, :);
120             else
121                 Indiv2(j)        = Pop(Index2);
122                 FObj_Indiv2(j, :) = MO_FObj_Pop(Index2, :);
123             end
124         end
125         //
126         // Crossover
127         //
128         for j=1:nb_couples
129             if (p_cross>grand(1, 1, "def")) then
130                 [x1, x2] = crossover_func(Indiv1(j), Indiv2(j), param);
131                 Indiv1(j) = x1;
132                 Indiv2(j) = x2;
133                 ToCompute_I1(j) = %T;
134                 ToCompute_I2(j) = %T;
135             else
136                 ToCompute_I1(j) = %F;
137                 ToCompute_I2(j) = %F;
138             end
139         end
140         //
141         // Mutation
142         //
143         for j=1:nb_couples
144             if (p_mut>grand(1, 1, "def")) then
145                 x1 = mutation_func(Indiv1(j), param);
146                 Indiv1(j) = x1;
147                 ToCompute_I1(j) = %T;
148             end
149             if (p_mut>grand(1, 1, "def")) then
150                 x2 = mutation_func(Indiv2(j), param);
151                 Indiv2(j) = x2;
152                 ToCompute_I2(j) = %T;
153             end
154         end
155         //
156         // Computation of the objective functions
157         //
158         for j=1:length(Indiv1)
159             if ToCompute_I1(j) then FObj_Indiv1(j, :) = _ga_f(Indiv1(j)); end
160             if ToCompute_I2(j) then FObj_Indiv2(j, :) = _ga_f(Indiv2(j)); end
161         end
162
163         // Reinit ToCompute lists
164         ToCompute_I1 = ToCompute_I1 & %F;
165         ToCompute_I2 = ToCompute_I2 & %F;
166
167         // We merge all the individuals in one list ...
168         All_Pop  = lstcat(Pop, Indiv1, Indiv2);
169         All_FObj = [FObj_Pop' FObj_Indiv1' FObj_Indiv2']';
170
171         // Compute the domination rank on all the population
172         Rank=DominationRank(All_FObj);
173
174         // Compute the crowding distance
175         MO_All_FObj = All_FObj;
176         Index    = 1:size(MO_All_FObj, 1);
177         Crowdist = zeros(size(MO_All_FObj, 1), 1);
178         for k=1:size(MO_All_FObj, 2)
179             [tmp, Index_List] = gsort(MO_All_FObj(:, k));
180             MO_All_FObj = MO_All_FObj(Index_List, :);
181             Index = Index(Index_List);
182             Crowdist(Index_List(1)) = %inf;
183             Crowdist(Index_List($)) = %inf;
184             _Max = max(MO_All_FObj(:, k));
185             _Min = min(MO_All_FObj(:, k));
186             for j=2:size(MO_All_FObj, 1)-1
187                 Crowdist(Index(j)) = Crowdist(Index(j)) - (MO_All_FObj(j+1, k) - MO_All_FObj(j-1, k)) / (_Max - _Min);
188             end
189         end
190         //
191         // Recombination
192         //
193         // We rank all the individual wrt to the partial order
194         for k=1:size(All_FObj, 1)-1
195             for j=k+1:size(All_FObj, 1)
196                 if (Rank(j)<Rank(k)) | ((Rank(j)==Rank(k)) & (Crowdist(j)>Crowdist(k))) then
197                     tmp           = Rank(k);
198                     Rank(k)       = Rank(j);
199                     Rank(j)       = tmp;
200                     tmp           = Crowdist(k);
201                     Crowdist(k)   = Crowdist(j);
202                     Crowdist(j)   = tmp;
203                     tmp           = All_Pop(k);
204                     All_Pop(k)    = All_Pop(j);
205                     All_Pop(j)    = tmp;
206                     tmp           = All_FObj(k, :);
207                     All_FObj(k, :) = All_FObj(j, :);
208                     All_FObj(j, :) = tmp;
209                 end
210             end
211         end
212         // Extraction and selection of the phenotype
213         FObj_Pop = All_FObj(1:pop_size, :);
214         // Extraction and selection of the genotype
215         Pop = list(All_Pop(1:pop_size));
216         // Extraction of the ranks and Crow distance
217         Rank     = Rank(1:pop_size);
218         Crowdist = Crowdist(1:pop_size);
219
220         if (Log) then
221             stop = output_func(i, nb_generation, Pop, FObj_Pop, param);
222             if stop then
223                 break
224             end
225         end
226     end
227
228     pop_opt      = codage_func(Pop, "decode", param);
229     fobj_pop_opt = FObj_Pop;
230 endfunction