[system] / trunk / pg / lib / Value / Formula.pm Repository: Repository Listing bbplugincoursesdistsnplrochestersystemwww

# Diff of /trunk/pg/lib/Value/Formula.pm

Revision 2686 Revision 2687
150 # 150 #
151 # Note: \$l is bigger if \$r can't be evaluated at one of the points 151 # Note: \$l is bigger if \$r can't be evaluated at one of the points
152 return 1 unless \$rvalues; 152 return 1 unless \$rvalues;
153 153
154 # 154 #
155 # Handle parameters
156 #
157 \$lvalues = \$l->{test_values}
159
160 #
155 # Look through the two lists to see if they are equal. 161 # Look through the two lists to see if they are equal.
156 # If not, return the comparison of the first unequal value 162 # If not, return the comparison of the first unequal value
157 # (not good for < and >, but OK for ==). 163 # (not good for < and >, but OK for ==).
158 # 164 #
159 my (\$i, \$cmp); 165 my (\$i, \$cmp);
169# 175#
170sub createPointValues { 176sub createPointValues {
171 my \$self = shift; 177 my \$self = shift;
172 my \$points = shift || \$self->{test_points} || \$self->createRandomPoints; 178 my \$points = shift || \$self->{test_points} || \$self->createRandomPoints;
173 my \$showError = shift; 179 my \$showError = shift;
174 my \$f = \$self->{f}; 180 my @vars = \$self->{context}->variables->variables;
175 \$f = \$self->{f} = \$self->perlFunction(undef,[\$self->{context}->variables->names]) 181 my @params = \$self->{context}->variables->parameters;
176 unless \$f; 182 my @zeros = @{\$self->{parameters} || [split('',"0" x scalar(@params))]};
183 my \$f = \$self->{f}; \$f = \$self->{f} = \$self->perlFunction(undef,[@vars,@params]) unless \$f;
177 184
178 my \$values = []; my \$v; 185 my \$values = []; my \$v;
179 foreach my \$p (@{\$points}) { 186 foreach my \$p (@{\$points}) {
180 \$v = eval {&\$f(@{\$p})}; 187 \$v = eval {&\$f(@{\$p},@zeros)};
181 if (!defined(\$v)) { 188 if (!defined(\$v)) {
182 return unless \$showError; 189 return unless \$showError;
183 Value::Error("Can't evaluate formula on test point (".join(',',@{\$p}).")"); 190 Value::Error("Can't evaluate formula on test point (".join(',',@{\$p}).")");
184 } 191 }
185 push @{\$values}, Value::makeValue(\$v); 192 push @{\$values}, Value::makeValue(\$v);
197 my \$self = shift; 204 my \$self = shift;
198 my \$num_points = @_[0]; 205 my \$num_points = @_[0];
199 \$num_points = int(\$self->getFlag('num_points',5)) unless defined(\$num_points); 206 \$num_points = int(\$self->getFlag('num_points',5)) unless defined(\$num_points);
200 \$num_points = 1 if \$num_points < 1; 207 \$num_points = 1 if \$num_points < 1;
201 208
202 my @vars = \$self->{context}->variables->names; 209 my @vars = \$self->{context}->variables->variables;
210 my @params = \$self->{context}->variables->parameters;
203 my @limits = \$self->getVariableLimits(@vars); 211 my @limits = \$self->getVariableLimits(@vars);
204 my @make = \$self->getVariableTypes(@vars); 212 my @make = \$self->getVariableTypes(@vars);
213 my @zeros = split('',"0" x scalar(@params));
205 my \$f = \$self->{f}; \$f = \$self->{f} = \$self->perlFunction(undef,[@vars]) unless \$f; 214 my \$f = \$self->{f}; \$f = \$self->{f} = \$self->perlFunction(undef,[@vars,@params]) unless \$f;
206 my \$seedRandom = \$self->{context}->flag('random_seed')? 'PGseedRandom' : 'seedRandom'; 215 my \$seedRandom = \$self->{context}->flag('random_seed')? 'PGseedRandom' : 'seedRandom';
207 my \$getRandom = \$self->{context}->flag('random_seed')? 'PGgetRandom' : 'getRandom'; 216 my \$getRandom = \$self->{context}->flag('random_seed')? 'PGgetRandom' : 'getRandom';
208 217
209 \$self->\$seedRandom; 218 \$self->\$seedRandom;
210 my \$points = []; my \$values = []; 219 my \$points = []; my \$values = [];
213 @P = (); \$i = 0; 222 @P = (); \$i = 0;
214 foreach my \$limit (@limits) { 223 foreach my \$limit (@limits) {
215 @p = (); foreach my \$I (@{\$limit}) {push @p, \$self->\$getRandom(@{\$I})} 224 @p = (); foreach my \$I (@{\$limit}) {push @p, \$self->\$getRandom(@{\$I})}
216 push @P, \$make[\$i++]->make(@p); 225 push @P, \$make[\$i++]->make(@p);
217 } 226 }
218 \$v = eval {&\$f(@P)}; 227 \$v = eval {&\$f(@P,@zeros)};
219 if (!defined(\$v)) {\$k++} else { 228 if (!defined(\$v)) {\$k++} else {
220 push @{\$points}, [@P]; 229 push @{\$points}, [@P];
221 push @{\$values}, Value::makeValue(\$v); 230 push @{\$values}, Value::makeValue(\$v);
222 \$k = 0; # reset count when we find a point 231 \$k = 0; # reset count when we find a point
223 } 232 }
293# 302#
294# Fake object for making reals (rather than use overhead of Value::Real) 303# Fake object for making reals (rather than use overhead of Value::Real)
295# 304#
296sub Value::Formula::number::make {shift; shift} 305sub Value::Formula::number::make {shift; shift}
297 306
307#
308# Find adaptive parameters, if any
309#
311 my \$l = shift; my \$r = shift;
312 my @params = @_; my \$d = scalar(@params);
313 return 0 if \$d == 0;
314 \$l->Error("Adaptive parameters can only be used for real-valued functions")
315 unless \$l->{tree}->isRealNumber;
316 #
317 # Get coefficient matrix of adaptive parameters
318 # and value vector for linear system
319 #
320 my (\$p,\$v) = \$l->createRandomPoints(\$d,1);
321 my @P = split('',"0" x \$d); my (\$f,\$F) = (\$l->{f},\$r->{f});
322 my @A = (); my @b = ();
323 foreach my \$i (0..\$d-1) {
324 my @a = (); my @p = @{\$p->[\$i]};
325 foreach my \$j (0..\$d-1) {
326 \$P[\$j] = 1; push(@a,&\$f(@p,@P)-\$v->[\$i]);
327 \$P[\$j] = 0;
328 }
329 push @A, [@a]; push @b, [&\$F(@p,@P)-\$v->[\$i]];
330 }
331 #
332 # Use MatrixReal1.pm to solve system of linear equations
333 #
334 my \$M = MatrixReal1->new(\$d,\$d); \$M->[0] = \@A;
335 my \$B = MatrixReal1->new(\$d,1); \$B->[0] = \@b;
336 (\$M,\$B) = \$M->normalize(\$B);
337 \$M = \$M->decompose_LR;
338 if ((\$d,\$B,\$M) = \$M->solve_LR(\$B)) {
339 if (\$d == 0) {
340 #
341 # Get values and recompute the points using them
342 #
343 my @a; foreach my \$r (@{\$B->[0]}) {push @a, \$r->[0]}
344 \$l->{parameters} = [@a];
345 \$l->createPointValues;
346 return 1;
347 }
348 }
349 \$l->Error("Can't solve for adaptive parameters");
350}
351
298## 352##
299## debugging routine 353## debugging routine
300## 354##
301#sub main::Format { 355#sub main::Format {
302# my \$v = scalar(@_) > 1? [@_]: shift; 356# my \$v = scalar(@_) > 1? [@_]: shift;

Legend:
 Removed from v.2686 changed lines Added in v.2687