1 /*
2 * hillclimb.cpp
3 * scorealign
4 *
5 * Created by Roger Dannenberg on 10/20/07.
6 * Copyright 2007 __MyCompanyName__. All rights reserved.
7 *
8 * Hillclimb is an abstract class for optimization. It models problems where
9 * you have a vector of parameters (stored as an array), a corresponding set
10 * of step sizes, and a non-linear function. The function is a virtual
11 * member function that subclasses must implement.
12 *
13 * The optimization algorithm is as follows:
14 * An initial set of parameters and step sizes is given.
15 *
16 * Estimate the partial derivatives with respect to each parameter value
17 * by taking a step along that dimension (use step sizes to determine
18 * how far to go) and calling the evaluate virtual function.
19 * Find the parameter that causes the maximum absolute change. If the
20 * change is positive for that parameter, take the step along that
21 * dimension. If the change is negative, take a negative step along that
22 * dimension.
23 *
24 * Repeat the previous paragraph as long as the result of evaluate is
25 * increasing. When it stops, you are at the top of a hill, a local
26 * maximum.
27 */
28
29
30 #include "stdio.h"
31 #include <stdlib.h>
32 #include "sautils.h"
33 #include "hillclimb.h"
34
35 #define HC_VERBOSE 0
36 #define V if (HC_VERBOSE)
37
~Hillclimb()38 Hillclimb::~Hillclimb()
39 {
40 if (parameters) FREE(parameters);
41 if (step_size) FREE(step_size);
42 if (min_param) FREE(min_param);
43 if (max_param) FREE(max_param);
44 }
45
46
setup(int n_)47 void Hillclimb::setup(int n_) {
48 n = n_;
49 parameters = ALLOC(double, n);
50 step_size = ALLOC(double, n);
51 min_param = ALLOC(double, n);
52 max_param = ALLOC(double, n);
53 }
54
55
set_parameters(double * p,double * ss,double * min_,double * max_,int plen)56 void Hillclimb::set_parameters(double *p, double *ss,
57 double *min_, double *max_, int plen)
58 {
59 parameters = p;
60 step_size = ss;
61 min_param = min_;
62 max_param = max_;
63 n = plen;
64 }
65
66 /* this optimize assumes that the surface is smooth enought that if the
67 * function decreases when parameter[i] increases, then the function will
68 * increase when parameter[i] decreases. The alternative version does more
69 * evaluation, but checks in both directions to find the best overall move.
70
71 double Hillclimb::optimize()
72 {
73 double best = evaluate();
74 while (true) {
75 printf("best %g ", best);
76 // eval partial derivatives
77 int i;
78 // variables to search for max partial derivative
79 double max = 0; // max of |dy| so far
80 int max_i; // index where max was found
81 int max_sign = 1; // sign of dy
82 double max_y; // value of evaluate() at max_i
83 // now search over all parameters for max change
84 for (i = 0; i < n; i++) {
85 int sign = 1; // sign of derivative in the +step direction
86 int step_direction = 1; // how to undo parameter variation
87 parameters[i] += step_size[i];
88 if (parameters[i] > max_param[i]) {
89 // try stepping in the other direction
90 parameters[i] -= step_size[i] * 2;
91 sign = -1;
92 step_direction = -1;
93 }
94
95 double y = evaluate();
96 // restore parameter i
97 parameters[i] -= step_size[i] * step_direction;
98
99 double dy = y - best;
100 if (dy < 0) {
101 dy = -dy;
102 sign = -sign;
103 }
104 // is this the best yet and legal move?
105 double proposal = parameters[i] + step_size[i] * sign;
106 if (dy > max && proposal <= max_param[i] &&
107 proposal >= min_param[i]) {
108 max = dy;
109 max_i = i;
110 max_y = y;
111 max_sign = sign;
112 }
113 }
114 // best move is parameter max_i in max_sign direction
115 parameters[max_i] += step_size[max_i] * max_sign;
116 printf("moved %d to %g", max_i, parameters[max_i]);
117 // what's the value now? put it in max_y
118 if (max_sign == -1) max_y = evaluate();
119 printf(" to get %g (vs. best %g)\n", max_y, best);
120 // otherwise, max_y already has the new value
121 if (max_y <= best) { // no improvement, we're done
122 parameters[max_i] -= step_size[max_i] * max_sign;
123 printf("\nCompleted hillclimbing, best %g\n", best);
124 return best;
125 }
126 // improvement because max_y higher than best:
127 best = max_y;
128 }
129 }
130 */
131
optimize(Report_fn_ptr report,void * cookie)132 double Hillclimb::optimize(Report_fn_ptr report, void *cookie)
133 {
134 double best = evaluate();
135 int iterations = 0;
136 while (true) {
137 (*report)(cookie, iterations, best);
138 V printf("best %g ", best);
139 // eval partial derivatives
140 int i;
141 // variables to search for max partial derivative
142 double max_y = best; // max of evaluate() so far
143 int max_i = 0; // index where best max was found
144 // the good parameter value for max_i:
145 double max_parameter = parameters[0];
146 // now search over all parameters for best improvement
147 for (i = 0; i < n; i++) {
148 V printf("optimize at %d param %g ", i, parameters[i]);
149 double save_param = parameters[i];
150 parameters[i] = save_param + step_size[i];
151 if (parameters[i] <= max_param[i]) {
152 double y = evaluate();
153 V printf("up->%g ", y);
154 if (y > max_y) {
155 V printf("NEW MAX! ");
156 max_y = y;
157 max_i = i;
158 max_parameter = parameters[i];
159 }
160 }
161 parameters[i] = save_param - step_size[i];
162 if (parameters[i] >= min_param[i]) {
163 double y = evaluate();
164 V printf("dn->%g ", y);
165 if (y > max_y) {
166 V printf("NEW MAX! ");
167 max_y = y;
168 max_i = i;
169 max_parameter = parameters[i];
170 }
171 }
172 parameters[i] = save_param;
173 V printf("\n");
174 }
175 iterations++; // for debugging, reporting
176 if (max_y <= best) { // no improvement, we're done
177 V printf("\nCompleted hillclimbing, best %g\n", best);
178 (*report)(cookie, iterations, best);
179 return best;
180 }
181 // improvement because max_y higher than best:
182 parameters[max_i] = max_parameter;
183 best = max_y;
184 }
185 }
186
187
188