[system] / trunk / pg / lib / Value / Formula.pm Repository:
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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}
158 if $l->AdaptParameters($r,$self->{context}->variables->parameters);
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#
310sub AdaptParameters {
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

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9