← Back to branch summary

~ubuntu-branches/ubuntu/utopic/mricron/utopic

« back to all changes in this revision

Viewing changes to npm/dmath/simopt.pas

  • Committer: Bazaar Package Importer
  • Author(s): Michael Hanke
  • Date: 2010-07-29 22:07:43 UTC
  • Revision ID: james.westby@ubuntu.com-20100729220743-q621ts2zj806gu0n
Tags: upstream-0.20100725.1~dfsg.1
Import upstream version 0.20100725.1~dfsg.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
npm/dmath/simopt.pas
 
1
{ **********************************************************************
 
2
  *                           Unit SIMOPT.PAS                          *
 
3
  *                             Version 1.0                            *
 
4
  *                     (c) J. Debord, August 2000                     *
 
5
  **********************************************************************
 
6
  This unit implements simulated annealing for function minimization
 
7
  **********************************************************************
 
8
  Reference: Program SIMANN.FOR by Bill Goffe
 
9
  (http://www.netlib.org/simann)
 
10
  ********************************************************************** }
 
11
 
 
12
unit SimOpt;
 
13
 
 
14
interface
 
15
 
 
16
uses
 
17
  FMath, Matrices, Optim, Stat;
 
18
 
 
19
const
 
20
  SA_Nt      : Integer = 5;    { Number of loops at constant temperature }
 
21
  SA_Ns      : Integer = 15;   { Number of loops before step adjustment }
 
22
  SA_Rt      : Float   = 0.9;  { Temperature reduction factor }
 
23
  SA_NCycles : Integer = 1;    { Number of cycles }
 
24
 
 
25
function SimAnn(Func           : TFuncNVar;
 
26
                X, Xmin, Xmax  : PVector;
 
27
                Lbound, Ubound : Integer;
 
28
                MaxIter        : Integer;
 
29
                Tol            : Float;
 
30
                var F_min      : Float) : Integer;
 
31
{ ----------------------------------------------------------------------
 
32
  Minimization of a function of several variables by simulated annealing
 
33
  ----------------------------------------------------------------------
 
34
  Input parameters : Func    = objective function to be minimized
 
35
                     X       = initial minimum coordinates
 
36
                     Xmin    = minimum value of X
 
37
                     Xmax    = maximum value of X
 
38
                     Lbound,
 
39
                     Ubound  = indices of first and last variables
 
40
                     MaxIter = max number of annealing steps
 
41
                     Tol     = required precision
 
42
  ----------------------------------------------------------------------
 
43
  Output parameter : X       = refined minimum coordinates
 
44
                     F_min   = function value at minimum
 
45
  ----------------------------------------------------------------------
 
46
  Possible results : OPT_OK
 
47
                     OPT_NON_CONV
 
48
  ---------------------------------------------------------------------- }
 
49
 
 
50
implementation
 
51
 
 
52
var
 
53
  LogFile : Text;  { Stores the result of each minimization step }
 
54
 
 
55
  procedure CreateLogFile;
 
56
  begin
 
57
    Assign(LogFile, LogFileName);
 
58
    Rewrite(LogFile);
 
59
  end;
 
60
 
 
61
  function InitTemp(Func           : TFuncNVar;
 
62
                    X, Xmin, Range : PVector;
 
63
                    Lbound, Ubound : Integer) : Float;
 
64
{ ----------------------------------------------------------------------
 
65
  Computes the initial temperature so that the probability
 
66
  of accepting an increase of the function is about 0.5
 
67
  ---------------------------------------------------------------------- }
 
68
  const
 
69
    N_EVAL = 50;  { Number of function evaluations }
 
70
  var
 
71
    T : Float;         { Temperature }
 
72
    F, F1 : Float;     { Function values }
 
73
    DeltaF : PVector;  { Function increases }
 
74
    N_inc : Integer;   { Number of function increases }
 
75
    I : Integer;       { Index of function evaluation }
 
76
    K : Integer;       { Index of parameter }
 
77
  begin
 
78
    DimVector(DeltaF, N_EVAL);
 
79
 
 
80
    T := 0.0;
 
81
    N_inc := 0;
 
82
    F := Func(X);
 
83
 
 
84
    { Compute N_EVAL function values, changing each parameter in turn }
 
85
    K := Lbound;
 
86
    for I := 1 to N_EVAL do
 
87
      begin
 
88
        X^[K] := Xmin^[K] + RanMar * Range^[K];
 
89
        F1 := Func(X);
 
90
        if F1 > F then
 
91
          begin
 
92
            Inc(N_inc);
 
93
            DeltaF^[N_inc] := F1 - F;
 
94
          end;
 
95
        F := F1;
 
96
        Inc(K);
 
97
        if K > Ubound then K := Lbound;
 
98
      end;
 
99
 
 
100
    { The median M of these N_eval values has a probability of 1/2.
 
101
      From Boltzmann's formula: Exp(-M/T) = 1/2 ==> T = M / Ln(2) }
 
102
    T := Median(DeltaF, 1, N_inc) / LN2;
 
103
    if T = 0.0 then T := 1.0;
 
104
    InitTemp := T;
 
105
 
 
106
    DelVector(DeltaF, N_EVAL);
 
107
  end;
 
108
 
 
109
  function ParamConv(X, Step        : PVector;
 
110
                     Lbound, Ubound : Integer;
 
111
                     Tol            : Float) : Boolean;
 
112
{ ----------------------------------------------------------------------
 
113
  Checks for convergence on parameters
 
114
  ---------------------------------------------------------------------- }
 
115
  var
 
116
    I : Integer;
 
117
    Conv : Boolean;
 
118
  begin
 
119
    I := Lbound;
 
120
    Conv := True;
 
121
    repeat
 
122
      Conv := Conv and (Step^[I] < FMax(Tol, Tol * Abs(X^[I])));
 
123
      Inc(I);
 
124
    until (Conv = False) or (I > Ubound);
 
125
    ParamConv := Conv;
 
126
  end;
 
127
 
 
128
  function Accept(DeltaF, T        : Float;
 
129
                  var N_inc, N_acc : Integer) : Boolean;
 
130
{ ----------------------------------------------------------------------
 
131
  Checks if a variation DeltaF of the function at temperature T is
 
132
  acceptable. Updates the counters N_inc (number of increases of the
 
133
  function) and N_acc (number of accepted increases).
 
134
  ---------------------------------------------------------------------- }
 
135
  begin
 
136
    if DeltaF < 0.0 then
 
137
      Accept := True
 
138
    else
 
139
      begin
 
140
        Inc(N_inc);
 
141
        if Expo(- DeltaF / T) > RanMar then
 
142
          begin
 
143
            Accept := True;
 
144
            Inc(N_acc);
 
145
          end
 
146
        else
 
147
          Accept := False;
 
148
      end;
 
149
  end;
 
150
 
 
151
  function SimAnnCycle(Func           : TFuncNVar;
 
152
                       X, Xmin, Xmax  : PVector;
 
153
                       Lbound, Ubound : Integer;
 
154
                       MaxIter        : Integer;
 
155
                       Tol            : Float;
 
156
                       var LogFile    : Text;
 
157
                       var F_min      : Float) : Integer;
 
158
{ ----------------------------------------------------------------------
 
159
  Performs one cycle of simulated annealing
 
160
  ---------------------------------------------------------------------- }
 
161
  const
 
162
    N_FACT = 2.0;  { Factor for step reduction }
 
163
  var
 
164
    I, Iter, J, K, N_inc, N_acc : Integer;
 
165
    F, F1, DeltaF, Ratio, T, OldX : Float;
 
166
    Range, Step, Xopt : PVector;
 
167
    Nacc : PIntVector;
 
168
  begin
 
169
    DimVector(Step, Ubound);
 
170
    DimVector(Xopt, Ubound);
 
171
    DimVector(Range, Ubound);
 
172
    DimIntVector(Nacc, Ubound);
 
173
 
 
174
    { Determine parameter range, step and optimum }
 
175
    for K := Lbound to Ubound do
 
176
      begin
 
177
        Range^[K] := Xmax^[K] - Xmin^[K];
 
178
        Step^[K] := 0.5 * Range^[K];
 
179
        Xopt^[K] := X^[K];
 
180
      end;
 
181
 
 
182
    { Initialize function values }
 
183
    F := Func(X);
 
184
    F_min := F;
 
185
 
 
186
    { Initialize temperature and iteration count }
 
187
    T := InitTemp(Func, X, Xmin, Range, Lbound, Ubound);
 
188
    Iter := 0;
 
189
 
 
190
    repeat
 
191
      { Perform SA_Nt evaluations at constant temperature }
 
192
      N_inc := 0; N_acc := 0;
 
193
      for I := 1 to SA_Nt do
 
194
        begin
 
195
          for J := 1 to SA_Ns do
 
196
            for K := Lbound to Ubound do
 
197
              begin
 
198
                { Save current parameter value }
 
199
                OldX := X^[K];
 
200
 
 
201
                { Pick new value, keeping it within Range }
 
202
                X^[K] := X^[K] + (2.0 * RanMar - 1.0) * Step^[K];
 
203
                if (X^[K] < Xmin^[K]) or (X^[K] > Xmax^[K]) then
 
204
                  X^[K] := Xmin^[K] + RanMar * Range^[K];
 
205
 
 
206
                { Compute new function value }
 
207
                F1 := Func(X);
 
208
                DeltaF := F1 - F;
 
209
 
 
210
                { Check for acceptance }
 
211
                if Accept(DeltaF, T, N_inc, N_acc) then
 
212
                  begin
 
213
                    Inc(Nacc^[K]);
 
214
                    F := F1;
 
215
                  end
 
216
                else
 
217
                  { Restore parameter value }
 
218
                  X^[K] := OldX;
 
219
 
 
220
                { Update minimum if necessary }
 
221
                if F < F_min then
 
222
                  begin
 
223
                    Xopt^[K] := X^[K];
 
224
                    F_min := F;
 
225
                  end;
 
226
              end;
 
227
 
 
228
          { Ajust step length to maintain an acceptance
 
229
            ratio of about 50% for each parameter }
 
230
          for K := Lbound to Ubound do
 
231
            begin
 
232
              Ratio := Int(Nacc^[K]) / Int(SA_Ns);
 
233
              if Ratio > 0.6 then
 
234
                begin
 
235
                  { Increase step length, keeping it within Range }
 
236
                  Step^[K] := Step^[K] * (1.0 + ((Ratio - 0.6) / 0.4) * N_FACT);
 
237
                  if Step^[K] > Range^[K] then Step^[K] := Range^[K];
 
238
                end
 
239
              else if Ratio < 0.4 then
 
240
                { Reduce step length }
 
241
                Step^[K] := Step^[K] / (1.0 + ((0.4 - Ratio) / 0.4) * N_FACT);
 
242
 
 
243
              { Restore counter }
 
244
              Nacc^[K] := 0;
 
245
            end;
 
246
        end;
 
247
 
 
248
      if WriteLogFile then
 
249
        WriteLn(LogFile, Iter:4, '   ', T:12, '   ', F:12, N_inc:6, N_acc:6);
 
250
 
 
251
      { Update temperature and iteration count }
 
252
      T := T * SA_Rt;
 
253
      Inc(Iter);
 
254
    until ParamConv(Xopt, Step, Lbound, Ubound, Tol) or (Iter > MaxIter);
 
255
 
 
256
    for K := Lbound to Ubound do
 
257
      X^[K] := Xopt^[K];
 
258
 
 
259
    DelVector(Step, Ubound);
 
260
    DelVector(Xopt, Ubound);
 
261
    DelVector(Range, Ubound);
 
262
    DelIntVector(Nacc, Ubound);
 
263
 
 
264
    if Iter > MaxIter then
 
265
      SimAnnCycle := OPT_NON_CONV
 
266
    else
 
267
      SimAnnCycle := OPT_OK;
 
268
  end;
 
269
 
 
270
  function SimAnn(Func           : TFuncNVar;
 
271
                  X, Xmin, Xmax  : PVector;
 
272
                  Lbound, Ubound : Integer;
 
273
                  MaxIter        : Integer;
 
274
                  Tol            : Float;
 
275
                  var F_min      : Float) : Integer;
 
276
  var
 
277
    Cycle, ErrCode : Integer;
 
278
  begin
 
279
    if WriteLogFile then
 
280
      CreateLogFile;
 
281
 
 
282
    { Initialize the Marsaglia random number generator
 
283
      using the standard Pascal generator }
 
284
    Randomize;
 
285
    RMarIn(System.Random(10000), System.Random(10000));
 
286
 
 
287
    Cycle := 1;
 
288
    repeat
 
289
      if WriteLogFile then
 
290
        begin
 
291
          WriteLn(LogFile, 'Simulated annealing: Cycle ', Cycle);
 
292
          WriteLn(LogFile);
 
293
          WriteLn(LogFile, 'Iter         T              F        Inc   Acc');
 
294
        end;
 
295
 
 
296
      ErrCode := SimAnnCycle(Func, X, Xmin, Xmax, Lbound, Ubound,
 
297
                             MaxIter, Tol, LogFile, F_min);
 
298
 
 
299
      Inc(Cycle);
 
300
    until (Cycle > SA_NCycles) or (ErrCode <> OPT_OK);
 
301
 
 
302
    if WriteLogFile then
 
303
      Close(LogFile);
 
304
 
 
305
    SimAnn := ErrCode;
 
306
  end;
 
307
 
 
308
end.