[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 5509 - (download) (as text) (annotate)
Sat Sep 15 00:56:51 2007 UTC (12 years, 3 months ago) by dpvc
File size: 16026 byte(s)
Formula objects and Context objects contain reference loops, which
prevent them from being freed properly by perl when they are no longer
needed.  This is a source of an important memory leak in WeBWorK.  The
problem has been fixed by using Scalar::Util::weaken for these
recursive references, so these objects can be freed properly when they
go out of scope.  This should cause an improvement in the memory usage
of the httpd child processes.

    1 ###########################################################################
    2 #
    3 #  Implements the Formula class.
    4 #
    5 package Value::Formula;
    6 my $pkg = 'Value::Formula';
    7 
    8 use strict;
    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   if ($l->AdaptParameters($r,$self->{context}->variables->parameters)) {
  186     my $avalues = $l->{test_adapt};
  187     my $tolerance    = $self->getFlag('tolerance',1E-4);
  188     my $isRelative   = $self->getFlag('tolType','relative') eq 'relative';
  189     my $zeroLevel    = $self->getFlag('zeroLevel',1E-14);
  190     my $zeroLevelTol = $self->getFlag('zeroLevelTol',1E-12);
  191     foreach $i (0..scalar(@{$lvalues})-1) {
  192       my $tol = $tolerance;
  193       my ($lv,$rv,$av) = ($lvalues->[$i]->value,$rvalues->[$i]->value,$avalues->[$i]->value);
  194       if ($isRelative) {
  195   if (abs($lv) <= $zeroLevel) {$tol = $zeroLevelTol}
  196                          else {$tol *= abs($lv)}
  197       }
  198       return $rv <=> $av unless abs($rv - $av) < $tol;
  199     }
  200     return 0;
  201   }
  202 
  203   #
  204   #  Look through the two lists of values to see if they are equal.
  205   #  If not, return the comparison of the first unequal value
  206   #    (not good for < and >, but OK for ==).
  207   #
  208   my $domainError = 0;
  209   foreach $i (0..scalar(@{$lvalues})-1) {
  210     if (ref($lvalues->[$i]) eq 'UNDEF' ^ ref($rvalues->[$i]) eq 'UNDEF') {$domainError = 1; next}
  211     $cmp = $lvalues->[$i] <=> $rvalues->[$i];
  212     return $cmp if $cmp;
  213   }
  214   $l->{domainMismatch} = $domainError;  # return this value
  215 }
  216 
  217 #
  218 #  Create the value list from a given set of test points
  219 #
  220 sub createPointValues {
  221   my $self = shift; my $context = $self->context;
  222   my $points = shift || $self->{test_points} || $self->createRandomPoints;
  223   my $showError = shift; my $cacheResults = shift;
  224   my @vars   = $context->variables->variables;
  225   my @params = $context->variables->parameters;
  226   my @zeros  = (0) x scalar(@params);
  227   my $f = $self->{f}; $f = $self->{f} = $self->perlFunction(undef,[@vars,@params]) unless $f;
  228   my $checkUndef = scalar(@params) == 0 && (shift || $self->getFlag('checkUndefinedPoints',0));
  229 
  230   my $values = []; my $v;
  231   foreach my $p (@{$points}) {
  232     $v = eval {&$f(@{$p},@zeros)};
  233     if (!defined($v) && !$checkUndef) {
  234       return unless $showError;
  235       Value::Error("Can't evaluate formula on test point (%s)",join(',',@{$p}));
  236     }
  237     if (defined($v)) {
  238       $v = Value::makeValue($v,context=>$context)->with(equation=>$self);
  239       $v->transferFlags("equation");
  240       push @{$values}, $v;
  241     } else {
  242       push @{$values}, $UNDEF;
  243     }
  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; my $context = $self->context;
  257   my $points = shift || $self->{test_points} || $self->createRandomPoints;
  258   my $showError = shift;
  259   my @vars   = $context->variables->variables;
  260   my @params = $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     $v = Value::makeValue($v,context=>$context)->with(equation=>$self);
  273     $v->transferFlags("equation");
  274     push @{$values}, $v;
  275   }
  276   $self->{test_adapt} = $values;
  277 }
  278 
  279 #
  280 #  Create a list of random points, making sure that the function
  281 #  is defined at the given points.  Error if we can't find enough.
  282 #
  283 sub createRandomPoints {
  284   my $self = shift; my $context = $self->context;
  285   my ($num_points,$include) = @_; my $cacheResults = !defined($num_points);
  286   $num_points = int($self->getFlag('num_points',5)) unless defined($num_points);
  287   $num_points = 1 if $num_points < 1;
  288 
  289   my @vars   = $context->variables->variables;
  290   my @params = $context->variables->parameters;
  291   my @limits = $self->getVariableLimits(@vars);
  292   my @make   = $self->getVariableTypes(@vars);
  293   my @zeros  = (0) x scalar(@params);
  294   my $f = $self->{f}; $f = $self->{f} = $self->perlFunction(undef,[@vars,@params]) unless $f;
  295   my $seedRandom = $context->flag('random_seed')? 'PGseedRandom' : 'seedRandom';
  296   my $getRandom  = $context->flag('random_seed')? 'PGgetRandom'  : 'getRandom';
  297   my $checkUndef = scalar(@params) == 0 && $self->getFlag('checkUndefinedPoints',0);
  298   my $max_undef  = $self->getFlag('max_undefined',$num_points);
  299 
  300   $self->$seedRandom;
  301   my $points = []; my $values = []; my $num_undef = 0;
  302   if ($include) {
  303     push(@{$points},@{$include});
  304     push(@{$values},@{$self->createPointValues($include,1,$cacheResults,$self->{checkundefinedPoints})});
  305   }
  306   my (@P,@p,$v,$i); my $k = 0;
  307   while (scalar(@{$points}) < $num_points+$num_undef && $k < 10) {
  308     @P = (); $i = 0;
  309     foreach my $limit (@limits) {
  310       @p = (); foreach my $I (@{$limit})
  311         {push @p, $context->Package("Real")->make($context,$self->$getRandom(@{$I}))}
  312       push @P, $make[$i++]->make($context,@p);
  313     }
  314     $v = eval {&$f(@P,@zeros)};
  315     if (!defined($v)) {
  316       if ($checkUndef && $num_undef < $max_undef) {
  317   push @{$points}, [@P];
  318   push @{$values}, $UNDEF;
  319   $num_undef++;
  320       }
  321       $k++;
  322     } else {
  323       $v = Value::makeValue($v,context=>$context)->with(equation=>$self);
  324       $v->transferFlags("equation");
  325       push @{$points}, [@P];
  326       push @{$values}, $v;
  327       $k = 0; # reset count when we find a point
  328     }
  329   }
  330 
  331   if ($k) {
  332     my $error = "Can't generate enough valid points for comparison";
  333     $error .= ':<div style="margin-left:1em">'.($context->{error}{message} || $@).'</div>'
  334       if ($self->getFlag('showTestPointErrors'));
  335     $error =~ s/ (in \S+ )?at line \d+.*//s;
  336     Value::Error($error);
  337   }
  338 
  339   return ($points,$values) unless $cacheResults;
  340   $self->{test_values} = $values;
  341   $self->{test_points} = $points;
  342 }
  343 
  344 #
  345 #  Get the array of variable limits
  346 #
  347 sub getVariableLimits {
  348   my $self = shift;
  349   my $userlimits = $self->{limits};
  350   if (defined($userlimits)) {
  351     $userlimits = [[[-$userlimits,$userlimits]]] unless ref($userlimits) eq 'ARRAY';
  352     $userlimits = [$userlimits] unless ref($userlimits->[0]) eq 'ARRAY';
  353     $userlimits = [$userlimits] if scalar(@_) == 1 && ref($userlimits->[0][0]) ne 'ARRAY';
  354     foreach my $I (@{$userlimits}) {$I = [$I] unless ref($I->[0]) eq 'ARRAY'};
  355   }
  356   $userlimits = [] unless $userlimits; my @limits;
  357   my $default;  $default = $userlimits->[0][0] if defined($userlimits->[0]);
  358   $default = $default || $self->getFlag('limits',[-2,2]);
  359   my $granularity = $self->getFlag('granularity',1000);
  360   my $resolution = $self->getFlag('resolution');
  361   my $i = 0;
  362   foreach my $x (@_) {
  363     my $def = $self->{context}{variables}{$x};
  364     my $limit = $userlimits->[$i++] || $def->{limits} || [];
  365     $limit = [$limit] if defined($limit->[0]) && ref($limit->[0]) ne 'ARRAY';
  366     push(@{$limit},$limit->[0] || $default) while (scalar(@{$limit}) < $def->{type}{length});
  367     pop(@{$limit}) while (scalar(@{$limit}) > $def->{type}{length});
  368     push @limits, $self->addGranularity($limit,$def,$granularity,$resolution);
  369   }
  370   return @limits;
  371 }
  372 
  373 #
  374 #  Add the granularity to the limit intervals
  375 #
  376 sub addGranularity {
  377   my $self = shift; my $limit = shift; my $def = shift;
  378   my $granularity = shift; my $resolution = shift;
  379   $granularity = $def->{granularity} || $granularity;
  380   $resolution = $def->{resolution} || $resolution;
  381   foreach my $I (@{$limit}) {
  382     my ($a,$b,$n) = @{$I}; $b = -$a unless defined $b;
  383     $I = [$a,$b,($n || $resolution || abs($b-$a)/$granularity)];
  384   }
  385   return $limit;
  386 }
  387 
  388 #
  389 #  Get the routines to make the coordinates of the points
  390 #
  391 sub getVariableTypes {
  392   my $self = shift;
  393   my @make;
  394   foreach my $x (@_) {
  395     my $type = $self->{context}{variables}{$x}{type};
  396     if ($type->{name} eq 'Number') {
  397       push @make,($type->{length} == 1)? 'Value::Formula::number': $self->Package("Complex");
  398     } else {
  399       push @make, $self->Package($type->{name});
  400     }
  401   }
  402   return @make;
  403 }
  404 
  405 #
  406 #  Fake object for making reals (rather than use overhead of Value::Real)
  407 #
  408 sub Value::Formula::number::make {shift; shift; shift->value}
  409 
  410 #
  411 #  Find adaptive parameters, if any
  412 #
  413 sub AdaptParameters {
  414   my $l = shift; my $r = shift;
  415   my @params = @_; my $d = scalar(@params);
  416   return 0 if $d == 0; return 0 unless $l->usesOneOf(@params);
  417   $l->Error("Adaptive parameters can only be used for real-valued formulas")
  418     unless $l->{tree}->isRealNumber;
  419   #
  420   #  Get coefficient matrix of adaptive parameters
  421   #  and value vector for linear system
  422   #
  423   my ($p,$v) = $l->createRandomPoints($d);
  424   my @P = (0) x $d; my ($f,$F) = ($l->{f},$r->{f});
  425   my @A = (); my @b = ();
  426   foreach my $i (0..$d-1) {
  427     my @a = (); my @p = @{$p->[$i]};
  428     foreach my $j (0..$d-1) {
  429       $P[$j] = 1; push(@a,(&$f(@p,@P)-$v->[$i])->value);
  430       $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       }
  457       my $context = $l->context;
  458       foreach my $i (0..$#a) {$context->{variables}{$params[$i]}{value} = $a[$i]}
  459       $l->{parameters} = [@a];
  460       $l->createAdaptedValues;
  461       return 1;
  462     }
  463   }
  464   $l->Error("Can't solve for adaptive parameters");
  465 }
  466 
  467 ##
  468 ##  debugging routine
  469 ##
  470 #sub main::Format {
  471 #  my $v = scalar(@_) > 1? [@_]: shift;
  472 #  $v = [%{$v}] if ref($v) eq 'HASH';
  473 #  return $v unless ref($v) eq 'ARRAY';
  474 #  my @V; foreach my $x (@{$v}) {push @V, main::Format($x)}
  475 #  return '['.join(",",@V).']';
  476 #}
  477 
  478 #
  479 #  Random number generator  (replaced by Value::WeBWorK.pm)
  480 #
  481 sub seedRandom {srand}
  482 sub getRandom {
  483   my $self = shift;
  484   my ($m,$M,$n) = @_; $n = 1 unless $n;
  485   return $m + $n*int(rand()*(int(($M-$m)/$n)+1));
  486 }
  487 
  488 ############################################
  489 #
  490 #  Check if the value of a formula is constant
  491 #    (could use shift->{tree}{isConstant}, but I don't trust it)
  492 #
  493 sub isConstant {
  494   my @vars = (%{shift->{variables}});
  495   return scalar(@vars) == 0;
  496 }
  497 
  498 #
  499 #  Check if the Formula includes one of the named variables
  500 #
  501 sub usesOneOf {
  502   my $self = shift;
  503   foreach my $x (@_) {return 1 if $self->{variables}{$x}}
  504   return 0;
  505 }
  506 
  507 ###########################################################################
  508 
  509 1;

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9