* Bug #13015 fixed - Genetic_algorithms: proper computation of Efficiency
[scilab.git] / scilab / modules / genetic_algorithms / macros / optim_nsga.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 // This file must be used under the terms of the CeCILL.
5 // This source file is licensed as described in the file COPYING, which
6 // you should have received as part of this distribution.  The terms
7 // are also available at
8 // http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
9
10 function [pop_opt, fobj_pop_opt, pop_init, fobj_pop_init] = optim_nsga(ga_f, pop_size, nb_generation, p_mut, p_cross, Log, param, sigma, pow)
11
12     [nargout, nargin] = argn();
13
14     if ~isdef("param","local") then
15         param = [];
16     end
17
18     [codage_func,err]    = get_param(param,"codage_func",coding_ga_identity);
19     [init_func,err]      = get_param(param,"init_func",init_ga_default);
20     [crossover_func,err] = get_param(param,"crossover_func",crossover_ga_default);
21     [mutation_func,err]  = get_param(param,"mutation_func",mutation_ga_default);
22     [selection_func,err] = get_param(param,"selection_func",selection_ga_elitist);
23     [nb_couples,err]     = get_param(param,"nb_couples",100);
24     [pressure,err]       = get_param(param,"pressure",0.05);
25
26     if ~isdef("ga_f","local") then
27         error(gettext("optim_moga: ga_f is mandatory"));
28     else
29         if typeof(ga_f)=="list" then
30             deff("y=_ga_f(x)","y=ga_f(1)(x, ga_f(2:$))");
31         else
32             deff("y=_ga_f(x)","y=ga_f(x)");
33         end
34     end
35
36     if ~isdef("pop_size","local") then
37         pop_size = 100;
38     end
39     if ~isdef("nb_generation","local") then
40         nb_generation = 10;
41     end
42     if ~isdef("p_mut","local") then
43         p_mut = 0.1;
44     end
45     if ~isdef("p_cross","local") then
46         p_cross = 0.1;
47     end
48     if ~isdef("Log","local") then
49         Log = %F;
50     end
51     if ~isdef("sigma","local") then
52         sigma = 0.01;
53     end
54     if ~isdef("pow","local") then
55         pow = 2;
56     end
57
58     // Initialization of the population
59     if (Log) then
60         printf(gettext("%s: Initialization of the population\n"),"optim_nsga");
61     end
62
63     Pop = list();
64     Pop = init_func(pop_size,param);
65
66
67     if (nargout>=3) then
68         pop_init = Pop;
69     end
70
71     // Code the individuals
72     Pop = codage_func(Pop,"code",param);
73
74     for i=1:length(Pop)
75         MO_FObj_Pop(i,:) = _ga_f(Pop(i));
76     end
77
78     // Compute the domination rank
79     for i=1:size(MO_FObj_Pop,1)
80         Index = 0;
81         for j=1:size(MO_FObj_Pop,1)
82             Index = Index + double(and(MO_FObj_Pop(i,:)<=MO_FObj_Pop(j,:)) & or(MO_FObj_Pop(i,:)<MO_FObj_Pop(j,:)));
83         end
84         FObj_Pop(i) = - (Index + 1);
85     end
86
87     if (nargout==4) then
88         fobj_pop_init = MO_FObj_Pop;
89     end
90
91     // The genetic algorithm
92     for i=1:nb_generation
93         if (Log) then
94             printf(gettext("%s: iteration %d / %d"), "optim_nsga", i, nb_generation);
95         end
96
97         // Computation of the niching penality
98         for j=1:size(MO_FObj_Pop,1)
99             Niching(j) = 0;
100             for k=1:size(MO_FObj_Pop,1)
101                 Distance = sqrt(sum((MO_FObj_Pop(j,:) - MO_FObj_Pop(k,:)).^2));
102                 if Distance < sigma then
103                     Niching(j) = Niching(j) + (1 - Distance / sigma)^pow;
104                 end
105             end
106         end
107
108         // Apply niching penality to fobj
109         FObj_Pop = FObj_Pop ./ Niching;
110
111         FObj_Pop_Max = max(FObj_Pop);
112         FObj_Pop_Min = min(FObj_Pop);
113
114         // Normalization of the efficiency
115         Efficiency = (1 - pressure) * (FObj_Pop_Max - FObj_Pop) / max([FObj_Pop_Max - FObj_Pop_Min %eps]) + pressure;
116
117         //
118         // Selection
119         //
120         Indiv1 = list();
121         Indiv2 = list();
122         Wheel = cumsum(Efficiency);
123
124         for j=1:nb_couples
125             // Selection of the first individual in the couple
126             Shoot = grand(1,1,"def")*Wheel($);
127             Index = find(Wheel<=Shoot);
128             if length(Index)>1 then Index = Index($); end;
129             if isempty(Index)  then Index = 1; end;
130             Indiv1(j)           = Pop(Index);
131             MO_FObj_Indiv1(j,:) = MO_FObj_Pop(Index,:);
132             // Selection of the second individual in the couple
133             Shoot = grand(1,1,"def")*Wheel($);
134             Index = find(Wheel<=Shoot);
135             if length(Index)>1 then Index = Index($); end;
136             if isempty(Index)  then Index = 1; end;
137             Indiv2(j)           = Pop(Index);
138             MO_FObj_Indiv2(j,:) = MO_FObj_Pop(Index,:);
139         end
140         //
141         // Crossover
142         //
143         for j=1:nb_couples
144             if (p_cross>grand(1,1,"def")) then
145                 [x1, x2]  = crossover_func(Indiv1(j), Indiv2(j),param);
146                 Indiv1(j) = x1;
147                 Indiv2(j) = x2;
148                 ToCompute_I1(j) = %T;
149                 ToCompute_I2(j) = %T;
150             else
151                 ToCompute_I1(j) = %F;
152                 ToCompute_I2(j) = %F;
153             end
154         end
155         //
156         // Mutation
157         //
158         for j=1:nb_couples
159             if (p_mut>grand(1,1,"def")) then
160                 x1 = mutation_func(Indiv1(j),param);
161                 Indiv1(j) = x1;
162                 ToCompute_I1(j) = %T;
163             end
164             if (p_mut>grand(1,1,"def")) then
165                 x2 = mutation_func(Indiv2(j),param);
166                 Indiv2(j) = x2;
167                 ToCompute_I2(j) = %T;
168             end
169         end
170         //
171         // Computation of the objective functions
172         //
173         for j=1:length(Indiv1)
174             if ToCompute_I1(j) then MO_FObj_Indiv1(j,:) = _ga_f(Indiv1(j)); end
175             if ToCompute_I2(j) then MO_FObj_Indiv2(j,:) = _ga_f(Indiv2(j)); end
176         end
177
178         // Reinit ToCompute lists
179         ToCompute_I1 = ToCompute_I1 & %F;
180         ToCompute_I2 = ToCompute_I2 & %F;
181
182         // Compute the domination rank
183         for j=1:size(MO_FObj_Indiv1,1)
184             // We compute the rank for Indiv1
185             Index1 = 0; Index2 = 0; Index3 = 0;
186             for k=1:size(MO_FObj_Indiv1,1)
187                 Index1 = Index1 + double(and(MO_FObj_Indiv1(j,:)<=MO_FObj_Indiv1(k,:)) & or(MO_FObj_Indiv1(j,:)<MO_FObj_Indiv1(k,:)));
188                 Index2 = Index2 + double(and(MO_FObj_Indiv1(j,:)<=MO_FObj_Indiv2(k,:)) & or(MO_FObj_Indiv1(j,:)<MO_FObj_Indiv2(k,:)));
189             end
190             for k=1:size(MO_FObj_Pop,1)
191                 Index3 = Index3 + and(MO_FObj_Indiv1(j,:)<=MO_FObj_Pop(k,:)) & or(MO_FObj_Indiv1(j,:)<MO_FObj_Pop(k,:));
192             end
193             FObj_Indiv1(j) = - (Index1 + Index2 + Index3 + 1);
194
195             // We compute the rank for Indiv2
196             Index1 = 0; Index2 = 0; Index3 = 0;
197             for k=1:size(MO_FObj_Indiv1,1)
198                 Index1 = Index1 + double(and(MO_FObj_Indiv2(j,:)<=MO_FObj_Indiv1(k,:)) & or(MO_FObj_Indiv2(j,:)<MO_FObj_Indiv1(k,:)));
199                 Index2 = Index2 + double(and(MO_FObj_Indiv2(j,:)<=MO_FObj_Indiv2(k,:)) & or(MO_FObj_Indiv2(j,:)<MO_FObj_Indiv2(k,:)));
200             end
201             for k=1:size(MO_FObj_Pop,1)
202                 Index3 = Index3 + double(and(MO_FObj_Indiv2(j,:)<=MO_FObj_Pop(k,:)) & or(MO_FObj_Indiv2(j,:)<MO_FObj_Pop(k,:)));
203             end
204             FObj_Indiv2(j) = - (Index1 + Index2 + Index3 + 1);
205         end
206
207         // We compute the rank for Pop
208         for j=1:size(MO_FObj_Pop,1)
209             Index1 = 0; Index2 = 0; Index3 = 0;
210             for k=1:size(MO_FObj_Indiv1,1)
211                 Index1 = Index1 + double(and(MO_FObj_Pop(j,:)<=MO_FObj_Indiv1(k,:)) & or(MO_FObj_Pop(j,:)<MO_FObj_Indiv1(k,:)));
212                 Index2 = Index2 + double(and(MO_FObj_Pop(j,:)<=MO_FObj_Indiv2(k,:)) & or(MO_FObj_Pop(j,:)<MO_FObj_Indiv2(k,:)));
213             end
214             for k=1:size(FObj_Pop,1)
215                 Index3 = Index3 + double(and(MO_FObj_Pop(j,:)<=MO_FObj_Pop(k,:)) & or(MO_FObj_Pop(j,:)<MO_FObj_Pop(k,:)));
216             end
217             FObj_Pop(j) = - (Index1 + Index2 + Index3 + 1);
218         end
219
220         //
221         // Recombination
222         //
223
224         [Pop,FObj_Pop,Efficiency,MO_FObj_Pop] = selection_func(Pop,Indiv1,Indiv2,FObj_Pop,FObj_Indiv1,FObj_Indiv2, ...
225         MO_FObj_Pop,MO_FObj_Indiv1,MO_FObj_Indiv2,param);
226         if (Log) then
227             printf(gettext(" - min / max value found = %f / %f\n"), min(FObj_Pop), max(FObj_Pop));
228         end
229     end
230
231     pop_opt      = codage_func(Pop,"decode",param);
232     fobj_pop_opt = MO_FObj_Pop;
233 endfunction