[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 3517 - (download) (as text) (annotate)
Sat Aug 13 21:31:46 2005 UTC (14 years, 6 months ago) by dpvc
File size: 15897 byte(s)
Pass the canBeInUnion call on to the Formula's root node

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

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9