[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 4822 - (download) (as text) (annotate)
Sat Mar 3 00:11:28 2007 UTC (12 years, 9 months ago) by dpvc
File size: 16554 byte(s)
Use zeroLevelTol values when checking answers that have adaptive parameters.

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

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9