[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 3716 - (download) (as text) (annotate)
Sun Oct 16 03:37:17 2005 UTC (14 years, 1 month ago) by dpvc
File size: 16143 byte(s)
In the past, when Value objects were inserted into strings, they would
automatically include parentheses so that if you had $f equal to 1+x
and $g equal to 1-x, then Formula("$f/$g") would mean (1+x)/(1-x)
rather than 1+(x/1)-x, which is what would happen as a straing string
substitution.

The problem is that this would also happen for real numbers, vectors,
and everything else, even when it wasn't necessary.  So if $x=Real(3),
then "Let x = $x" would be "Let x = (3)".

I have changed the behavior of the string concatenation for Value
objects so that parentheses are only added in a few cases: for
Formulas, Complex numbers, and Unions.  This makes the other Value
objects work more like regular variables in strings, but might cause
some problems with strings that are used as formulas.  For example, if
$a = Real(-3), then "x + 2 $a" will become "x + 2 -3", or "x-1" rather
than the expected "x - 6".  (The old approach would have made it "x +
2 (-3)" which would have worked properly).  For the most part, it is
easier to use something like "x + 2*$a" or even "x" + 2*$a in this
case, so the extra trouble of having to avoid parentheses when you
really meant to substitute the value into a string didn't seem worth
it.

    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     foreach $i (0..scalar(@{$lvalues})-1) {
  209       my $tol = $tolerance;
  210       my ($lv,$rv,$av) = ($lvalues->[$i]->value,$rvalues->[$i]->value,$avalues->[$i]->value);
  211       $tol *= abs($lv) if $isRelative && abs($lv) > $zeroLevel;
  212       return $rv <=> $av unless abs($rv - $av) < $tol;
  213     }
  214     return 0;
  215   }
  216 
  217   #
  218   #  Look through the two lists of values to see if they are equal.
  219   #  If not, return the comparison of the first unequal value
  220   #    (not good for < and >, but OK for ==).
  221   #
  222   my $domainError = 0;
  223   foreach $i (0..scalar(@{$lvalues})-1) {
  224     if (ref($lvalues->[$i]) eq 'UNDEF' ^ ref($rvalues->[$i]) eq 'UNDEF') {$domainError = 1; next}
  225     $cmp = $lvalues->[$i] <=> $rvalues->[$i];
  226     return $cmp if $cmp;
  227   }
  228   $l->{domainMismatch} = $domainError;  # return this value
  229 }
  230 
  231 #
  232 #  Create the value list from a given set of test points
  233 #
  234 sub createPointValues {
  235   my $self = shift;
  236   my $points = shift || $self->{test_points} || $self->createRandomPoints;
  237   my $showError = shift; my $cacheResults = shift;
  238   my @vars   = $self->{context}->variables->variables;
  239   my @params = $self->{context}->variables->parameters;
  240   my @zeros  = (0) x scalar(@params);
  241   my $f = $self->{f}; $f = $self->{f} = $self->perlFunction(undef,[@vars,@params]) unless $f;
  242   my $checkUndef = scalar(@params) == 0 && (shift || $self->getFlag('checkUndefinedPoints',0));
  243 
  244   my $values = []; my $v;
  245   foreach my $p (@{$points}) {
  246     $v = eval {&$f(@{$p},@zeros)};
  247     if (!defined($v) && !$checkUndef) {
  248       return unless $showError;
  249       Value::Error("Can't evaluate formula on test point (%s)",join(',',@{$p}));
  250     }
  251     push @{$values}, (defined($v)? Value::makeValue($v): $UNDEF);
  252   }
  253   if ($cacheResults) {
  254     $self->{test_points} = $points;
  255     $self->{test_values} = $values;
  256   }
  257   return $values;
  258 }
  259 
  260 #
  261 #  Create the adapted value list for the test points
  262 #
  263 sub createAdaptedValues {
  264   my $self = shift;
  265   my $points = shift || $self->{test_points} || $self->createRandomPoints;
  266   my $showError = shift;
  267   my @vars   = $self->{context}->variables->variables;
  268   my @params = $self->{context}->variables->parameters;
  269   my $f = $self->{f}; $f = $self->{f} = $self->perlFunction(undef,[@vars,@params]) unless $f;
  270 
  271   my $values = []; my $v;
  272   my @adapt = @{$self->{parameters}};
  273   foreach my $p (@{$points}) {
  274     $v = eval {&$f(@{$p},@adapt)};
  275     if (!defined($v)) {
  276       return unless $showError;
  277       Value::Error("Can't evaluate formula on test point (%s) with parameters (%s)",
  278        join(',',@{$p}),join(',',@adapt));
  279     }
  280     push @{$values}, Value::makeValue($v);
  281   }
  282   $self->{test_adapt} = $values;
  283 }
  284 
  285 #
  286 #  Create a list of random points, making sure that the function
  287 #  is defined at the given points.  Error if we can't find enough.
  288 #
  289 sub createRandomPoints {
  290   my $self = shift;
  291   my ($num_points,$include) = @_; my $cacheResults = !defined($num_points);
  292   $num_points = int($self->getFlag('num_points',5)) unless defined($num_points);
  293   $num_points = 1 if $num_points < 1;
  294 
  295   my @vars   = $self->{context}->variables->variables;
  296   my @params = $self->{context}->variables->parameters;
  297   my @limits = $self->getVariableLimits(@vars);
  298   my @make   = $self->getVariableTypes(@vars);
  299   my @zeros  = (0) x scalar(@params);
  300   my $f = $self->{f}; $f = $self->{f} = $self->perlFunction(undef,[@vars,@params]) unless $f;
  301   my $seedRandom = $self->{context}->flag('random_seed')? 'PGseedRandom' : 'seedRandom';
  302   my $getRandom  = $self->{context}->flag('random_seed')? 'PGgetRandom'  : 'getRandom';
  303   my $checkUndef = scalar(@params) == 0 && $self->getFlag('checkUndefinedPoints',0);
  304   my $max_undef  = $self->getFlag('max_undefined',$num_points);
  305 
  306   $self->$seedRandom;
  307   my $points = []; my $values = []; my $num_undef = 0;
  308   if ($include) {
  309     push(@{$points},@{$include});
  310     push(@{$values},@{$self->createPointValues($include,1,$cacheResults,$self->{checkundefinedPoints})});
  311   }
  312   my (@P,@p,$v,$i); my $k = 0;
  313   while (scalar(@{$points}) < $num_points+$num_undef && $k < 10) {
  314     @P = (); $i = 0;
  315     foreach my $limit (@limits) {
  316       @p = (); foreach my $I (@{$limit}) {push @p, $self->$getRandom(@{$I})}
  317       push @P, $make[$i++]->make(@p);
  318     }
  319     $v = eval {&$f(@P,@zeros)};
  320     if (!defined($v)) {
  321       if ($checkUndef && $num_undef < $max_undef) {
  322   push @{$points}, [@P];
  323   push @{$values}, $UNDEF;
  324   $num_undef++;
  325       }
  326       $k++;
  327     } else {
  328       push @{$points}, [@P];
  329       push @{$values}, Value::makeValue($v);
  330       $k = 0; # reset count when we find a point
  331     }
  332   }
  333 
  334   Value::Error("Can't generate enough valid points for comparison") if $k;
  335   return ($points,$values) unless $cacheResults;
  336   $self->{test_values} = $values;
  337   $self->{test_points} = $points;
  338 }
  339 
  340 #
  341 #  Get the array of variable limits
  342 #
  343 sub getVariableLimits {
  344   my $self = shift;
  345   my $userlimits = $self->{limits};
  346   if (defined($userlimits)) {
  347     $userlimits = [[[-$userlimits,$userlimits]]] unless ref($userlimits) eq 'ARRAY';
  348     $userlimits = [$userlimits] unless ref($userlimits->[0]) eq 'ARRAY';
  349     $userlimits = [$userlimits] if scalar(@_) == 1 && ref($userlimits->[0][0]) ne 'ARRAY';
  350     foreach my $I (@{$userlimits}) {$I = [$I] unless ref($I->[0]) eq 'ARRAY'};
  351   }
  352   $userlimits = [] unless $userlimits; my @limits;
  353   my $default;  $default = $userlimits->[0][0] if defined($userlimits->[0]);
  354   $default = $default || $self->{context}{flags}{limits} || [-2,2];
  355   my $granularity = $self->getFlag('granularity',1000);
  356   my $resolution = $self->getFlag('resolution');
  357   my $i = 0;
  358   foreach my $x (@_) {
  359     my $def = $self->{context}{variables}{$x};
  360     my $limit = $userlimits->[$i++] || $def->{limits} || [];
  361     $limit = [$limit] if defined($limit->[0]) && ref($limit->[0]) ne 'ARRAY';
  362     push(@{$limit},$limit->[0] || $default) while (scalar(@{$limit}) < $def->{type}{length});
  363     pop(@{$limit}) while (scalar(@{$limit}) > $def->{type}{length});
  364     push @limits, $self->addGranularity($limit,$def,$granularity,$resolution);
  365   }
  366   return @limits;
  367 }
  368 
  369 #
  370 #  Add the granularity to the limit intervals
  371 #
  372 sub addGranularity {
  373   my $self = shift; my $limit = shift; my $def = shift;
  374   my $granularity = shift; my $resolution = shift;
  375   $granularity = $def->{granularity} || $granularity;
  376   $resolution = $def->{resolution} || $resolution;
  377   foreach my $I (@{$limit}) {
  378     my ($a,$b,$n) = @{$I}; $b = -$a unless defined $b;
  379     $I = [$a,$b,($n || $resolution || abs($b-$a)/$granularity)];
  380   }
  381   return $limit;
  382 }
  383 
  384 #
  385 #  Get the routines to make the coordinates of the points
  386 #
  387 sub getVariableTypes {
  388   my $self = shift;
  389   my @make;
  390   foreach my $x (@_) {
  391     my $type = $self->{context}{variables}{$x}{type};
  392     if ($type->{name} eq 'Number') {
  393       push @make,($type->{length} == 1)? 'Value::Formula::number': 'Value::Complex';
  394     } else {
  395       push @make, "Value::$type->{name}";
  396     }
  397   }
  398   return @make;
  399 }
  400 
  401 #
  402 #  Fake object for making reals (rather than use overhead of Value::Real)
  403 #
  404 sub Value::Formula::number::make {shift; shift}
  405 
  406 #
  407 #  Find adaptive parameters, if any
  408 #
  409 sub AdaptParameters {
  410   my $l = shift; my $r = shift;
  411   my @params = @_; my $d = scalar(@params);
  412   return 0 if $d == 0; return 0 unless $l->usesOneOf(@params);
  413   $l->Error("Adaptive parameters can only be used for real-valued functions")
  414     unless $l->{tree}->isRealNumber;
  415   #
  416   #  Get coefficient matrix of adaptive parameters
  417   #  and value vector for linear system
  418   #
  419   my ($p,$v) = $l->createRandomPoints($d);
  420   my @P = (0) x $d; my ($f,$F) = ($l->{f},$r->{f});
  421   my @A = (); my @b = ();
  422   foreach my $i (0..$d-1) {
  423     my @a = (); my @p = @{$p->[$i]};
  424     foreach my $j (0..$d-1) {
  425       $P[$j] = 1; push(@a,&$f(@p,@P)-$v->[$i]);
  426       $P[$j] = 0;
  427     }
  428     push @A, [@a]; push @b, [&$F(@p,@P)-$v->[$i]];
  429   }
  430   #
  431   #  Use MatrixReal1.pm to solve system of linear equations
  432   #
  433   my $M = MatrixReal1->new($d,$d); $M->[0] = \@A;
  434   my $B = MatrixReal1->new($d,1);  $B->[0] = \@b;
  435   ($M,$B) = $M->normalize($B);
  436   $M = $M->decompose_LR;
  437   if (($d,$B,$M) = $M->solve_LR($B)) {
  438     if ($d == 0) {
  439       #
  440       #  Get parameter values and recompute the points using them
  441       #
  442       my @a; my $i = 0; my $max = Value::Real->new($l->getFlag('max_adapt',1E8));
  443       foreach my $row (@{$B->[0]}) {
  444   if (abs($row->[0]) > $max) {
  445     $l->Error("Constant of integration is too large: %s\n(maximum allowed is %s)",
  446         $row->[0]->string,$max->string) if ($params[$i] eq 'C0');
  447     $l->Error("Adaptive constant is too large: %s = %s\n(maximum allowed is %s)",
  448         $params[$i],$row->[0]->string,$max->string);
  449   }
  450   push @a, $row->[0]; $i++;
  451       }
  452       $l->{parameters} = [@a];
  453       $l->createAdaptedValues;
  454       return 1;
  455     }
  456   }
  457   $l->Error("Can't solve for adaptive parameters");
  458 }
  459 
  460 sub usesOneOf {
  461   my $self = shift;
  462   foreach my $x (@_) {return 1 if $self->{variables}{$x}}
  463   return 0;
  464 }
  465 
  466 ##
  467 ##  debugging routine
  468 ##
  469 #sub main::Format {
  470 #  my $v = scalar(@_) > 1? [@_]: shift;
  471 #  $v = [%{$v}] if ref($v) eq 'HASH';
  472 #  return $v unless ref($v) eq 'ARRAY';
  473 #  my @V; foreach my $x (@{$v}) {push @V, main::Format($x)}
  474 #  return '['.join(",",@V).']';
  475 #}
  476 
  477 #
  478 #  Random number generator  (replaced by Value::WeBWorK.pm)
  479 #
  480 sub seedRandom {srand}
  481 sub getRandom {
  482   my $self = shift;
  483   my ($m,$M,$n) = @_; $n = 1 unless $n;
  484   return $m + $n*int(rand()*(int(($M-$m)/$n)+1));
  485 }
  486 
  487 ############################################
  488 #
  489 #  Check if the value of a formula is constant
  490 #    (could use shift->{tree}{isConstant}, but I don't trust it)
  491 #
  492 sub isConstant {
  493   my @vars = (%{shift->{variables}});
  494   return scalar(@vars) == 0;
  495 }
  496 
  497 ############################################
  498 #
  499 #  Provide output formats
  500 #
  501 sub stringify {
  502   my $self = shift;
  503   return $self->TeX if $$Value::context->flag('StringifyAsTeX');
  504   $self->string;
  505 }
  506 
  507 ###########################################################################
  508 
  509 1;

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9