[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 5516 Revision 5517
179 # 179 #
180 # Handle adaptive parameters: 180 # Handle adaptive parameters:
181 # Get the tolerances, and check each adapted point relative 181 # Get the tolerances, and check each adapted point relative
182 # to the ORIGINAL correct answer. (This will have to be 182 # to the ORIGINAL correct answer. (This will have to be
183 # fixed if we ever do adaptive parameters for non-real formulas) 183 # fixed if we ever do adaptive parameters for non-real formulas)
184 #
185 # FIXME: it doesn't make sense to apply the ORIGINAL value's
186 # tolerance, and causes problems when the values
187 # differ in magnitude by much. Gavin has found several
188 # situations where this is a problem.
184 # 189 #
185 if ($l->AdaptParameters($r,$self->{context}->variables->parameters)) { 190 if ($l->AdaptParameters($r,$self->{context}->variables->parameters)) {
186 my $avalues = $l->{test_adapt}; 191 my $avalues = $l->{test_adapt};
187 my $tolerance = $self->getFlag('tolerance',1E-4); 192 my $tolerance = $self->getFlag('tolerance',1E-4);
188 my $isRelative = $self->getFlag('tolType','relative') eq 'relative'; 193 my $isRelative = $self->getFlag('tolType','relative') eq 'relative';
410# 415#
411# Find adaptive parameters, if any 416# Find adaptive parameters, if any
412# 417#
413sub AdaptParameters { 418sub AdaptParameters {
414 my $l = shift; my $r = shift; 419 my $l = shift; my $r = shift;
415 my @params = @_; my $d = scalar(@params); 420 my @params = @_; my $d = scalar(@params); my $D;
416 return 0 if $d == 0; return 0 unless $l->usesOneOf(@params); 421 return 0 if $d == 0; return 0 unless $l->usesOneOf(@params);
417 $l->Error("Adaptive parameters can only be used for real-valued formulas") 422 $l->Error("Adaptive parameters can only be used for real-valued formulas")
418 unless $l->{tree}->isRealNumber; 423 unless $l->{tree}->isRealNumber;
424
425 #
426 # Try up to three times (the random points might not work the first time)
427 #
428 foreach my $attempt (1..3) {
419 # 429 #
420 # Get coefficient matrix of adaptive parameters 430 # Get coefficient matrix of adaptive parameters
421 # and value vector for linear system 431 # and value vector for linear system
422 # 432 #
423 my ($p,$v) = $l->createRandomPoints($d); 433 my ($p,$v) = $l->createRandomPoints($d);
424 my @P = (0) x $d; my ($f,$F) = ($l->{f},$r->{f}); 434 my @P = (0) x $d; my ($f,$F) = ($l->{f},$r->{f});
425 my @A = (); my @b = (); 435 my @A = (); my @b = ();
426 foreach my $i (0..$d-1) { 436 foreach my $i (0..$d-1) {
427 my @a = (); my @p = @{$p->[$i]}; 437 my @a = (); my @p = @{$p->[$i]};
428 foreach my $j (0..$d-1) { 438 foreach my $j (0..$d-1) {
429 $P[$j] = 1; push(@a,(&$f(@p,@P)-$v->[$i])->value); 439 $P[$j] = 1; push(@a,(&$f(@p,@P)-$v->[$i])->value);
430 $P[$j] = 0; 440 $P[$j] = 0;
431 }
432 push @A, [@a]; push @b, [(&$F(@p,@P)-$v->[$i])->value];
433 }
434 #
435 # Use MatrixReal1.pm to solve system of linear equations
436 #
437 my $M = MatrixReal1->new($d,$d); $M->[0] = \@A;
438 my $B = MatrixReal1->new($d,1); $B->[0] = \@b;
439 ($M,$B) = $M->normalize($B);
440 $M = $M->decompose_LR;
441 if (($d,$B,$M) = $M->solve_LR($B)) {
442 if ($d == 0) {
443 #
444 # Get parameter values and recompute the points using them
445 #
446 my @a; my $i = 0; my $max = $l->getFlag('max_adapt',1E8);
447 foreach my $row (@{$B->[0]}) {
448 if (abs($row->[0]) > $max) {
449 $max = Value::makeValue($max); $row->[0] = Value::makeValue($row->[0]);
450 $l->Error(["Constant of integration is too large: %s\n(maximum allowed is %s)",
451 $row->[0]->string,$max->string]) if $params[$i] eq 'C0';
452 $l->Error(["Adaptive constant is too large: %s = %s\n(maximum allowed is %s)",
453 $params[$i],$row->[0]->string,$max->string]);
454 }
455 push @a, $row->[0]; $i++;
456 } 441 }
442 push @A, [@a]; push @b, [(&$F(@p,@P)-$v->[$i])->value];
443 }
444 #
445 # Use MatrixReal1.pm to solve system of linear equations
446 #
447 my $M = MatrixReal1->new($d,$d); $M->[0] = \@A;
448 my $B = MatrixReal1->new($d,1); $B->[0] = \@b;
449 ($M,$B) = $M->normalize($B);
450 $M = $M->decompose_LR;
451 if (($D,$B,$M) = $M->solve_LR($B)) {
452 if ($D == 0) {
453 #
454 # Get parameter values and recompute the points using them
455 #
456 my @a; my $i = 0; my $max = $l->getFlag('max_adapt',1E8);
457 foreach my $row (@{$B->[0]}) {
458 if (abs($row->[0]) > $max) {
459 $max = Value::makeValue($max); $row->[0] = Value::makeValue($row->[0]);
460 $l->Error(["Constant of integration is too large: %s\n(maximum allowed is %s)",
461 $row->[0]->string,$max->string]) if $params[$i] eq 'C0' or $params[$i] eq 'n00';
462 $l->Error(["Adaptive constant is too large: %s = %s\n(maximum allowed is %s)",
463 $params[$i],$row->[0]->string,$max->string]);
464 }
465 push @a, $row->[0]; $i++;
466 }
457 my $context = $l->context; 467 my $context = $l->context;
458 foreach my $i (0..$#a) {$context->{variables}{$params[$i]}{value} = $a[$i]} 468 foreach my $i (0..$#a) {$context->{variables}{$params[$i]}{value} = $a[$i]}
459 $l->{parameters} = [@a]; 469 $l->{parameters} = [@a];
460 $l->createAdaptedValues; 470 $l->createAdaptedValues;
461 return 1; 471 return 1;
472 }
462 } 473 }
463 } 474 }
464 $l->Error("Can't solve for adaptive parameters"); 475 $l->Error("Can't solve for adaptive parameters");
465} 476}
466 477

Legend:
Removed from v.5516  
changed lines
  Added in v.5517

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9