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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6223 - (download) (as text) (annotate)
Wed Apr 14 14:26:12 2010 UTC (9 years, 10 months ago) by dpvc
File size: 17347 byte(s)
Make sure inherited formulas get the right variables

    1 ###########################################################################
    2 #
    3 #  Implements the Formula class.
    4 #
    5 package Value::Formula;
    6 my $pkg = 'Value::Formula';
    7 
    8 use strict; no strict "refs";
    9 our @ISA = qw(Parser Value);
   10 
   11 my $UNDEF = bless {}, "UNDEF"; # used for undefined points
   12 
   13 
   14 #
   15 #  Call Parser to make the new Formula
   16 #
   17 sub new {
   18   my $self = shift;
   19   my $f = $self->SUPER::new(@_);
   20   foreach my $id ('open','close') {$f->{$id} = $f->{tree}{$id}}
   21   return $f;
   22 }
   23 
   24 #
   25 #  Create a new Formula with no string
   26 #    (we'll fill in its tree by hand)
   27 #
   28 sub blank {shift->SUPER::new(@_)}
   29 
   30 #
   31 #  with() changes tree element as well
   32 #    as the formula itself.
   33 #
   34 sub with {
   35   my $self = shift; my %hash = @_;
   36   $self = $self->SUPER::with(@_);
   37   $self->{tree} = $self->{tree}->copy($self); # make a new copy pointing to the new equation.
   38   foreach my $id (keys(%hash)) {$self->{tree}{$id} = $hash{$id}}
   39   return $self;
   40 }
   41 
   42 #
   43 #  Get the type from the tree
   44 #
   45 sub typeRef {(shift)->{tree}->typeRef}
   46 sub length {(shift)->{tree}->typeRef->{length}}
   47 
   48 sub isZero {(shift)->{tree}{isZero}}
   49 sub isOne {(shift)->{tree}{isOne}}
   50 
   51 sub isSetOfReals {(shift)->{tree}->isSetOfReals}
   52 sub canBeInUnion {(shift)->{tree}->canBeInUnion}
   53 
   54 sub transferFlags {}
   55 
   56 ############################################
   57 #
   58 #  Create a BOP from two operands
   59 #
   60 #  Get the context and variables from the left and right operands
   61 #    if they are formulas
   62 #  Make them into Value objects if they aren't already.
   63 #  Convert '+' to union for intervals or unions.
   64 #  Make a new BOP with the two operands.
   65 #  Record the variables.
   66 #  Evaluate the formula if it is constant.
   67 #
   68 sub bop {
   69   my $bop = shift;
   70   my ($self,$l,$r) = Value::checkOpOrder(@_);
   71   my $class = ref($self) || $self;
   72   my $call = $self->context->{method}{$bop};
   73   my $formula = $self->blank($self->context);
   74   if (ref($r) eq $class || ref($r) eq $pkg) {
   75     $formula->{context} = $r->{context};
   76     $r = $r->{tree}->copy($formula);
   77   } else {
   78     $r = $self->new($r)->{tree}->copy($formula);
   79   }
   80   if (ref($l) eq $class || ref($l) eq $pkg) {
   81     $formula->{context} = $l->{context};
   82     $l = $l->{tree}->copy($formula);
   83   } else {
   84     $l = $self->new($l)->{tree}->copy($formula);
   85   }
   86   $bop = 'U' if $bop eq '+' &&
   87     ($l->type =~ m/Interval|Set|Union/ || $r->type =~ m/Interval|Set|Union/);
   88   $formula->{tree} = $formula->Item("BOP")->new($formula,$bop,$l,$r);
   89   $formula->{variables} = $formula->{tree}->getVariables;
   90   return $formula;
   91 }
   92 
   93 sub add   {bop('+',@_)}
   94 sub sub   {bop('-',@_)}
   95 sub mult  {bop('*',@_)}
   96 sub div   {bop('/',@_)}
   97 sub power {bop('**',@_)}
   98 sub cross {bop('><',@_)}
   99 
  100 #
  101 #  Make dot work for vector operands
  102 #
  103 sub _dot   {
  104   my ($l,$r,$flag) = @_;
  105   if ($l->promotePrecedence($r)) {return $r->_dot($l,!$flag)}
  106   return bop('.',@_) if ($l->type eq 'Vector' || $l->{isVector}) &&
  107      Value::isValue($r) && ($r->type eq 'Vector' || $r->{isVector});
  108   $l->SUPER::_dot($r,$flag);
  109 }
  110 
  111 sub pdot {'('.(shift->stringify).')'}
  112 
  113 #
  114 #  Call the Parser::Function call function
  115 #
  116 sub call {
  117   my $self = shift; my $name = shift;
  118   Parser::Function->call($name,$self);
  119 }
  120 
  121 ############################################
  122 #
  123 #  Form the negation of a formula
  124 #
  125 sub neg {
  126   my $self = shift;
  127   my $formula = $self->blank($self->context);
  128   $formula->{variables} = $self->{variables};
  129   $formula->{tree} = $formula->Item("UOP")->new($formula,'u-',$self->{tree}->copy($formula));
  130   return $formula;
  131 }
  132 
  133 #
  134 #  Form the function atan2 function call on two operands
  135 #
  136 sub atan2 {
  137   my ($self,$l,$r) = Value::checkOpOrderWithPromote(@_);
  138   Parser::Function->call('atan2',$l,$r);
  139 }
  140 
  141 #
  142 #  Other overloaded functions
  143 #
  144 sub sin  {shift->call('sin',@_)}
  145 sub cos  {shift->call('cos',@_)}
  146 sub abs  {shift->call('abs',@_)}
  147 sub exp  {shift->call('exp',@_)}
  148 sub log  {shift->call('log',@_)}
  149 sub sqrt {shift->call('sqrt',@_)}
  150 
  151 sub twiddle {shift->call('conj',@_)}
  152 
  153 ############################################
  154 #
  155 #  Compare two functions for equality
  156 #
  157 sub compare {
  158   my ($l,$r) = @_; my $self = $l;
  159   my $context = $self->context;
  160   $r = $context->Package("Formula")->new($context,$r) unless Value::isFormula($r);
  161   Value::Error("Formulas from different contexts can't be compared")
  162     unless $l->{context} == $r->{context};
  163 
  164   #
  165   #  Get the test points and evaluate the functions at those points
  166   #
  167   ##  FIXME: Check given points for consistency
  168   ##  FIXME: make arrays if only a single value is given
  169   ##  FIXME: insert additional values if vars in use in formula aren't all the vars in the context
  170   my $points  = $l->{test_points} || $l->createRandomPoints(undef,$l->{test_at});
  171   my $lvalues = $l->{test_values} || $l->createPointValues($points,1,1);
  172   my $rvalues = $r->createPointValues($points,0,1,$l->{checkUndefinedPoints});
  173   #
  174   # Note: $l is bigger if $r can't be evaluated at one of the points
  175   return 1 unless $rvalues;
  176 
  177   my ($i, $cmp);
  178 
  179   #
  180   #  Handle adaptive parameters:
  181   #    Get the tolerances, and check each adapted point relative
  182   #    to the ORIGINAL correct answer.  (This will have to be
  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.
  189   #
  190   if ($l->AdaptParameters($r,$self->{context}->variables->parameters)) {
  191     my $avalues      = $l->{test_adapt};
  192     my $tolerance    = $self->getFlag('tolerance',1E-4);
  193     my $isRelative   = $self->getFlag('tolType','relative') eq 'relative';
  194     my $zeroLevel    = $self->getFlag('zeroLevel',1E-14);
  195     my $zeroLevelTol = $self->getFlag('zeroLevelTol',1E-12);
  196     foreach $i (0..scalar(@{$lvalues})-1) {
  197       my $tol = $tolerance;
  198       my ($lv,$rv,$av) = ($lvalues->[$i]->value,$rvalues->[$i]->value,$avalues->[$i]->value);
  199       if ($isRelative) {
  200   if (abs($lv) <= $zeroLevel) {$tol = $zeroLevelTol}
  201                          else {$tol *= abs($lv)}
  202       }
  203       return $rv <=> $av unless abs($rv - $av) < $tol;
  204     }
  205     return 0;
  206   }
  207 
  208   #
  209   #  Look through the two lists of values to see if they are equal.
  210   #  If not, return the comparison of the first unequal value
  211   #    (not good for < and >, but OK for ==).
  212   #
  213   my $domainError = 0;
  214   foreach $i (0..scalar(@{$lvalues})-1) {
  215     if (ref($lvalues->[$i]) eq 'UNDEF' ^ ref($rvalues->[$i]) eq 'UNDEF') {$domainError = 1; next}
  216     $cmp = $lvalues->[$i] <=> $rvalues->[$i];
  217     return $cmp if $cmp;
  218   }
  219   $l->{domainMismatch} = $domainError;  # return this value
  220 }
  221 
  222 #
  223 #  Inherit should make sure the tree is copied
  224 #  (so it's nodes point to the correct equation, for one thing)
  225 #
  226 sub inherit {
  227   my $self = shift; my $tree = $self->{tree};
  228   $self = $self->SUPER::inherit(@_);
  229   $self->{tree} = $tree->copy($self);
  230   $self->{variables} = $tree->getVariables;
  231   return $self;
  232 }
  233 
  234 #
  235 #  Don't inherit test values or adapted values, or other temporary items
  236 #
  237 sub noinherit {
  238   my $self = shift;
  239   ($self->SUPER::noinherit(@_),"test_values","test_adapt","tree","string","variables",
  240     "f","stack","ref","tokens","values","space","domainMismatch");
  241 }
  242 
  243 
  244 #
  245 #  Create the value list from a given set of test points
  246 #
  247 sub createPointValues {
  248   my $self = shift; my $context = $self->context;
  249   my $points = shift || $self->{test_points} || $self->createRandomPoints;
  250   my $showError = shift; my $cacheResults = shift;
  251   my @vars   = $context->variables->variables;
  252   my @params = $context->variables->parameters;
  253   my @zeros  = (0) x scalar(@params);
  254   my $f = $self->{f}; $f = $self->{f} = $self->perlFunction(undef,[@vars,@params]) unless $f;
  255   my $checkUndef = scalar(@params) == 0 && (shift || $self->getFlag('checkUndefinedPoints',0));
  256 
  257   my $values = []; my $v;
  258   foreach my $p (@{$points}) {
  259     $v = eval {&$f(@{$p},@zeros)};
  260     if (!defined($v) && !$checkUndef) {
  261       return unless $showError;
  262       Value::Error("Can't evaluate formula on test point (%s)",join(',',@{$p}));
  263     }
  264     if (defined($v)) {
  265       $v = Value::makeValue($v,context=>$context)->with(equation=>$self);
  266       $v->transferFlags("equation");
  267       push @{$values}, $v;
  268     } else {
  269       push @{$values}, $UNDEF;
  270     }
  271   }
  272   if ($cacheResults) {
  273     $self->{test_points} = $points;
  274     $self->{test_values} = $values;
  275   }
  276   return $values;
  277 }
  278 
  279 #
  280 #  Create the adapted value list for the test points
  281 #
  282 sub createAdaptedValues {
  283   my $self = shift; my $context = $self->context;
  284   my $points = shift || $self->{test_points} || $self->createRandomPoints;
  285   my $showError = shift;
  286   my @vars   = $context->variables->variables;
  287   my @params = $context->variables->parameters;
  288   my $f = $self->{f}; $f = $self->{f} = $self->perlFunction(undef,[@vars,@params]) unless $f;
  289 
  290   my $values = []; my $v;
  291   my @adapt = @{$self->{parameters}};
  292   foreach my $p (@{$points}) {
  293     $v = eval {&$f(@{$p},@adapt)};
  294     if (!defined($v)) {
  295       return unless $showError;
  296       Value::Error("Can't evaluate formula on test point (%s) with parameters (%s)",
  297        join(',',@{$p}),join(',',@adapt));
  298     }
  299     $v = Value::makeValue($v,context=>$context)->with(equation=>$self);
  300     $v->transferFlags("equation");
  301     push @{$values}, $v;
  302   }
  303   $self->{test_adapt} = $values;
  304 }
  305 
  306 #
  307 #  Create a list of random points, making sure that the function
  308 #  is defined at the given points.  Error if we can't find enough.
  309 #
  310 sub createRandomPoints {
  311   my $self = shift; my $context = $self->context;
  312   my ($num_points,$include,$noErrors) = @_; my $cacheResults = !defined($num_points);
  313   $num_points = int($self->getFlag('num_points',5)) unless defined($num_points);
  314   $num_points = 1 if $num_points < 1;
  315 
  316   my @vars   = $context->variables->variables;
  317   my @params = $context->variables->parameters;
  318   my @limits = $self->getVariableLimits(@vars);
  319   my @make   = $self->getVariableTypes(@vars);
  320   my @zeros  = (0) x scalar(@params);
  321   my $f = $self->{f}; $f = $self->{f} = $self->perlFunction(undef,[@vars,@params]) unless $f;
  322   my $seedRandom = $context->flag('random_seed')? 'PGseedRandom' : 'seedRandom';
  323   my $getRandom  = $context->flag('random_seed')? 'PGgetRandom'  : 'getRandom';
  324   my $checkUndef = scalar(@params) == 0 && $self->getFlag('checkUndefinedPoints',0);
  325   my $max_undef  = $self->getFlag('max_undefined',$num_points);
  326 
  327   $self->$seedRandom;
  328   my $points = []; my $values = []; my $num_undef = 0;
  329   if ($include) {
  330     push(@{$points},@{$include});
  331     push(@{$values},@{$self->createPointValues($include,1,$cacheResults,$self->{checkundefinedPoints})});
  332   }
  333   my (@P,@p,$v,$i); my $k = 0;
  334   while (scalar(@{$points}) < $num_points+$num_undef && $k < 10) {
  335     @P = (); $i = 0;
  336     foreach my $limit (@limits) {
  337       @p = (); foreach my $I (@{$limit})
  338         {push @p, $context->Package("Real")->make($context,$self->$getRandom(@{$I}))}
  339       push @P, $make[$i++]->make($context,@p);
  340     }
  341     $v = eval {&$f(@P,@zeros)};
  342     if (!defined($v)) {
  343       if ($checkUndef && $num_undef < $max_undef) {
  344   push @{$points}, [@P];
  345   push @{$values}, $UNDEF;
  346   $num_undef++;
  347       }
  348       $k++;
  349     } else {
  350       $v = Value::makeValue($v,context=>$context)->with(equation=>$self);
  351       $v->transferFlags("equation");
  352       push @{$points}, [@P];
  353       push @{$values}, $v;
  354       $k = 0; # reset count when we find a point
  355     }
  356   }
  357 
  358   if ($k) {
  359     my $error = "Can't generate enough valid points for comparison";
  360     $error .= ':<div style="margin-left:1em">'.($context->{error}{message} || $@).'</div>'
  361       if ($self->getFlag('showTestPointErrors'));
  362     $error =~ s! (in \S+ )?at line \d+.*</div>!</div>!s;
  363     Value::Error($error) unless $noErrors;
  364     return ($points,$values,1);
  365   }
  366 
  367   return ($points,$values) unless $cacheResults;
  368   $self->{test_values} = $values;
  369   $self->{test_points} = $points;
  370 }
  371 
  372 #
  373 #  Get the array of variable limits
  374 #
  375 sub getVariableLimits {
  376   my $self = shift;
  377   my $userlimits = $self->{limits};
  378   if (defined($userlimits)) {
  379     $userlimits = [[[-$userlimits,$userlimits]]] unless ref($userlimits) eq 'ARRAY';
  380     $userlimits = [$userlimits] unless ref($userlimits->[0]) eq 'ARRAY';
  381     $userlimits = [$userlimits] if scalar(@_) == 1 && ref($userlimits->[0][0]) ne 'ARRAY';
  382     foreach my $I (@{$userlimits}) {$I = [$I] unless ref($I->[0]) eq 'ARRAY'};
  383   }
  384   $userlimits = [] unless $userlimits; my @limits;
  385   my $default;  $default = $userlimits->[0][0] if defined($userlimits->[0]);
  386   $default = $default || $self->getFlag('limits',[-2,2]);
  387   my $granularity = $self->getFlag('granularity',1000);
  388   my $resolution = $self->getFlag('resolution');
  389   my $i = 0;
  390   foreach my $x (@_) {
  391     my $def = $self->{context}{variables}{$x};
  392     my $limit = $userlimits->[$i++] || $def->{limits} || [];
  393     $limit = [$limit] if defined($limit->[0]) && ref($limit->[0]) ne 'ARRAY';
  394     push(@{$limit},$limit->[0] || $default) while (scalar(@{$limit}) < $def->{type}{length});
  395     pop(@{$limit}) while (scalar(@{$limit}) > $def->{type}{length});
  396     push @limits, $self->addGranularity($limit,$def,$granularity,$resolution);
  397   }
  398   return @limits;
  399 }
  400 
  401 #
  402 #  Add the granularity to the limit intervals
  403 #
  404 sub addGranularity {
  405   my $self = shift; my $limit = shift; my $def = shift;
  406   my $granularity = shift; my $resolution = shift;
  407   $granularity = $def->{granularity} || $granularity;
  408   $resolution = $def->{resolution} || $resolution;
  409   foreach my $I (@{$limit}) {
  410     my ($a,$b,$n) = @{$I}; $b = -$a unless defined $b;
  411     $I = [$a,$b,($n || $resolution || abs($b-$a)/$granularity)];
  412   }
  413   return $limit;
  414 }
  415 
  416 #
  417 #  Get the routines to make the coordinates of the points
  418 #
  419 sub getVariableTypes {
  420   my $self = shift;
  421   my @make;
  422   foreach my $x (@_) {
  423     my $type = $self->{context}{variables}{$x}{type};
  424     if ($type->{name} eq 'Number') {
  425       push @make,($type->{length} == 1)? 'Value::Formula::number': $self->Package("Complex");
  426     } else {
  427       push @make, $self->Package($type->{name});
  428     }
  429   }
  430   return @make;
  431 }
  432 
  433 #
  434 #  Fake object for making reals (rather than use overhead of Value::Real)
  435 #
  436 sub Value::Formula::number::make {shift; shift; shift->value}
  437 
  438 #
  439 #  Find adaptive parameters, if any
  440 #
  441 sub AdaptParameters {
  442   my $l = shift; my $r = shift;
  443   my @params = @_; my $d = scalar(@params); my $D;
  444   return 0 if $d == 0; return 0 unless $l->usesOneOf(@params);
  445   $l->Error("Adaptive parameters can only be used for real-valued formulas")
  446     unless $l->{tree}->isRealNumber;
  447 
  448   #
  449   #  Try up to three times (the random points might not work the first time)
  450   #
  451   foreach my $attempt (1..3) {
  452     #
  453     #  Get coefficient matrix of adaptive parameters
  454     #  and value vector for linear system
  455     #
  456     my ($p,$v) = $l->createRandomPoints($d);
  457     my @P = (0) x $d; my ($f,$F) = ($l->{f},$r->{f});
  458     my @A = (); my @b = ();
  459     foreach my $i (0..$d-1) {
  460       my @a = (); my @p = @{$p->[$i]};
  461       foreach my $j (0..$d-1) {
  462   $P[$j] = 1;
  463   my $y = eval {&$f(@p,@P)};
  464   $l->Error(["Can't evaluate correct answer at adapted point (%s)",join(",",@$p,@P)])
  465          unless defined $y;
  466   push(@a,($y-$v->[$i])->value);
  467   $P[$j] = 0;
  468       }
  469       my $y = eval {&$F(@p,@P)}; return unless defined $y;
  470       push @A, [@a]; push @b, [($y-$v->[$i])->value];
  471     }
  472     #
  473     #  Use MatrixReal1.pm to solve system of linear equations
  474     #
  475     my $M = MatrixReal1->new($d,$d); $M->[0] = \@A;
  476     my $B = MatrixReal1->new($d,1);  $B->[0] = \@b;
  477     ($M,$B) = $M->normalize($B);
  478     $M = $M->decompose_LR;
  479     if (abs($M->det_LR) > 1E-6) {
  480       if (($D,$B,$M) = $M->solve_LR($B)) {
  481     if ($D == 0) {
  482       #
  483       #  Get parameter values and recompute the points using them
  484       #
  485       my @a; my $i = 0; my $max = $l->getFlag('max_adapt',1E8);
  486       foreach my $row (@{$B->[0]}) {
  487       if (abs($row->[0]) > $max) {
  488         $max = Value::makeValue($max); $row->[0] = Value::makeValue($row->[0]);
  489         $l->Error(["Constant of integration is too large: %s\n(maximum allowed is %s)",
  490          $row->[0]->string,$max->string]) if $params[$i] eq 'C0' or $params[$i] eq 'n00';
  491         $l->Error(["Adaptive constant is too large: %s = %s\n(maximum allowed is %s)",
  492          $params[$i],$row->[0]->string,$max->string]);
  493       }
  494       push @a, $row->[0]; $i++;
  495       }
  496       my $context = $l->context;
  497       foreach my $i (0..$#a) {$context->{variables}{$params[$i]}{value} = $a[$i]}
  498       $l->{parameters} = [@a];
  499       $l->createAdaptedValues;
  500       return 1;
  501     }
  502       }
  503     }
  504   }
  505   $l->Error("Can't solve for adaptive parameters");
  506 }
  507 
  508 ##
  509 ##  debugging routine
  510 ##
  511 #sub main::Format {
  512 #  my $v = scalar(@_) > 1? [@_]: shift;
  513 #  $v = [%{$v}] if ref($v) eq 'HASH';
  514 #  return $v unless ref($v) eq 'ARRAY';
  515 #  my @V; foreach my $x (@{$v}) {push @V, main::Format($x)}
  516 #  return '['.join(",",@V).']';
  517 #}
  518 
  519 #
  520 #  Random number generator  (replaced by Value::WeBWorK.pm)
  521 #
  522 sub seedRandom {srand}
  523 sub getRandom {
  524   my $self = shift;
  525   my ($m,$M,$n) = @_; $n = 1 unless $n;
  526   return $m + $n*int(rand()*(int(($M-$m)/$n)+1));
  527 }
  528 
  529 ############################################
  530 #
  531 #  Check if the value of a formula is constant
  532 #    (could use shift->{tree}{isConstant}, but I don't trust it)
  533 #
  534 sub isConstant {
  535   my @vars = (%{shift->{variables}});
  536   return scalar(@vars) == 0;
  537 }
  538 
  539 #
  540 #  Check if the Formula includes one of the named variables
  541 #
  542 sub usesOneOf {
  543   my $self = shift;
  544   foreach my $x (@_) {return 1 if $self->{variables}{$x}}
  545   return 0;
  546 }
  547 
  548 ###########################################################################
  549 
  550 1;

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9