[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 5990 - (download) (as text) (annotate)
Thu Feb 5 14:57:48 2009 UTC (10 years, 10 months ago) by dpvc
File size: 17211 byte(s)
Fix inheritance to copy the tree rather than just the pointer.  This
makes sure that the nodes point to the proper equation.

In the routine for adapting to parameters, check that the coefficient
matrix isn't singular, and try new random points if it is.

    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;
  228   $self = $self->SUPER::inherit(@_);
  229   $self->{tree} = $self->{tree}->copy($self);
  230   return $self;
  231 }
  232 
  233 #
  234 #  Don't inherit test values or adapted values, or other temporary items
  235 #
  236 sub noinherit {
  237   my $self = shift;
  238   ($self->SUPER::noinherit(@_),"test_values","test_adapt","tree","string","variables",
  239     "f","stack","ref","tokens","values","space","domainMismatch");
  240 }
  241 
  242 
  243 #
  244 #  Create the value list from a given set of test points
  245 #
  246 sub createPointValues {
  247   my $self = shift; my $context = $self->context;
  248   my $points = shift || $self->{test_points} || $self->createRandomPoints;
  249   my $showError = shift; my $cacheResults = shift;
  250   my @vars   = $context->variables->variables;
  251   my @params = $context->variables->parameters;
  252   my @zeros  = (0) x scalar(@params);
  253   my $f = $self->{f}; $f = $self->{f} = $self->perlFunction(undef,[@vars,@params]) unless $f;
  254   my $checkUndef = scalar(@params) == 0 && (shift || $self->getFlag('checkUndefinedPoints',0));
  255 
  256   my $values = []; my $v;
  257   foreach my $p (@{$points}) {
  258     $v = eval {&$f(@{$p},@zeros)};
  259     if (!defined($v) && !$checkUndef) {
  260       return unless $showError;
  261       Value::Error("Can't evaluate formula on test point (%s)",join(',',@{$p}));
  262     }
  263     if (defined($v)) {
  264       $v = Value::makeValue($v,context=>$context)->with(equation=>$self);
  265       $v->transferFlags("equation");
  266       push @{$values}, $v;
  267     } else {
  268       push @{$values}, $UNDEF;
  269     }
  270   }
  271   if ($cacheResults) {
  272     $self->{test_points} = $points;
  273     $self->{test_values} = $values;
  274   }
  275   return $values;
  276 }
  277 
  278 #
  279 #  Create the adapted value list for the test points
  280 #
  281 sub createAdaptedValues {
  282   my $self = shift; my $context = $self->context;
  283   my $points = shift || $self->{test_points} || $self->createRandomPoints;
  284   my $showError = shift;
  285   my @vars   = $context->variables->variables;
  286   my @params = $context->variables->parameters;
  287   my $f = $self->{f}; $f = $self->{f} = $self->perlFunction(undef,[@vars,@params]) unless $f;
  288 
  289   my $values = []; my $v;
  290   my @adapt = @{$self->{parameters}};
  291   foreach my $p (@{$points}) {
  292     $v = eval {&$f(@{$p},@adapt)};
  293     if (!defined($v)) {
  294       return unless $showError;
  295       Value::Error("Can't evaluate formula on test point (%s) with parameters (%s)",
  296        join(',',@{$p}),join(',',@adapt));
  297     }
  298     $v = Value::makeValue($v,context=>$context)->with(equation=>$self);
  299     $v->transferFlags("equation");
  300     push @{$values}, $v;
  301   }
  302   $self->{test_adapt} = $values;
  303 }
  304 
  305 #
  306 #  Create a list of random points, making sure that the function
  307 #  is defined at the given points.  Error if we can't find enough.
  308 #
  309 sub createRandomPoints {
  310   my $self = shift; my $context = $self->context;
  311   my ($num_points,$include) = @_; my $cacheResults = !defined($num_points);
  312   $num_points = int($self->getFlag('num_points',5)) unless defined($num_points);
  313   $num_points = 1 if $num_points < 1;
  314 
  315   my @vars   = $context->variables->variables;
  316   my @params = $context->variables->parameters;
  317   my @limits = $self->getVariableLimits(@vars);
  318   my @make   = $self->getVariableTypes(@vars);
  319   my @zeros  = (0) x scalar(@params);
  320   my $f = $self->{f}; $f = $self->{f} = $self->perlFunction(undef,[@vars,@params]) unless $f;
  321   my $seedRandom = $context->flag('random_seed')? 'PGseedRandom' : 'seedRandom';
  322   my $getRandom  = $context->flag('random_seed')? 'PGgetRandom'  : 'getRandom';
  323   my $checkUndef = scalar(@params) == 0 && $self->getFlag('checkUndefinedPoints',0);
  324   my $max_undef  = $self->getFlag('max_undefined',$num_points);
  325 
  326   $self->$seedRandom;
  327   my $points = []; my $values = []; my $num_undef = 0;
  328   if ($include) {
  329     push(@{$points},@{$include});
  330     push(@{$values},@{$self->createPointValues($include,1,$cacheResults,$self->{checkundefinedPoints})});
  331   }
  332   my (@P,@p,$v,$i); my $k = 0;
  333   while (scalar(@{$points}) < $num_points+$num_undef && $k < 10) {
  334     @P = (); $i = 0;
  335     foreach my $limit (@limits) {
  336       @p = (); foreach my $I (@{$limit})
  337         {push @p, $context->Package("Real")->make($context,$self->$getRandom(@{$I}))}
  338       push @P, $make[$i++]->make($context,@p);
  339     }
  340     $v = eval {&$f(@P,@zeros)};
  341     if (!defined($v)) {
  342       if ($checkUndef && $num_undef < $max_undef) {
  343   push @{$points}, [@P];
  344   push @{$values}, $UNDEF;
  345   $num_undef++;
  346       }
  347       $k++;
  348     } else {
  349       $v = Value::makeValue($v,context=>$context)->with(equation=>$self);
  350       $v->transferFlags("equation");
  351       push @{$points}, [@P];
  352       push @{$values}, $v;
  353       $k = 0; # reset count when we find a point
  354     }
  355   }
  356 
  357   if ($k) {
  358     my $error = "Can't generate enough valid points for comparison";
  359     $error .= ':<div style="margin-left:1em">'.($context->{error}{message} || $@).'</div>'
  360       if ($self->getFlag('showTestPointErrors'));
  361     $error =~ s/ (in \S+ )?at line \d+.*//s;
  362     Value::Error($error);
  363   }
  364 
  365   return ($points,$values) unless $cacheResults;
  366   $self->{test_values} = $values;
  367   $self->{test_points} = $points;
  368 }
  369 
  370 #
  371 #  Get the array of variable limits
  372 #
  373 sub getVariableLimits {
  374   my $self = shift;
  375   my $userlimits = $self->{limits};
  376   if (defined($userlimits)) {
  377     $userlimits = [[[-$userlimits,$userlimits]]] unless ref($userlimits) eq 'ARRAY';
  378     $userlimits = [$userlimits] unless ref($userlimits->[0]) eq 'ARRAY';
  379     $userlimits = [$userlimits] if scalar(@_) == 1 && ref($userlimits->[0][0]) ne 'ARRAY';
  380     foreach my $I (@{$userlimits}) {$I = [$I] unless ref($I->[0]) eq 'ARRAY'};
  381   }
  382   $userlimits = [] unless $userlimits; my @limits;
  383   my $default;  $default = $userlimits->[0][0] if defined($userlimits->[0]);
  384   $default = $default || $self->getFlag('limits',[-2,2]);
  385   my $granularity = $self->getFlag('granularity',1000);
  386   my $resolution = $self->getFlag('resolution');
  387   my $i = 0;
  388   foreach my $x (@_) {
  389     my $def = $self->{context}{variables}{$x};
  390     my $limit = $userlimits->[$i++] || $def->{limits} || [];
  391     $limit = [$limit] if defined($limit->[0]) && ref($limit->[0]) ne 'ARRAY';
  392     push(@{$limit},$limit->[0] || $default) while (scalar(@{$limit}) < $def->{type}{length});
  393     pop(@{$limit}) while (scalar(@{$limit}) > $def->{type}{length});
  394     push @limits, $self->addGranularity($limit,$def,$granularity,$resolution);
  395   }
  396   return @limits;
  397 }
  398 
  399 #
  400 #  Add the granularity to the limit intervals
  401 #
  402 sub addGranularity {
  403   my $self = shift; my $limit = shift; my $def = shift;
  404   my $granularity = shift; my $resolution = shift;
  405   $granularity = $def->{granularity} || $granularity;
  406   $resolution = $def->{resolution} || $resolution;
  407   foreach my $I (@{$limit}) {
  408     my ($a,$b,$n) = @{$I}; $b = -$a unless defined $b;
  409     $I = [$a,$b,($n || $resolution || abs($b-$a)/$granularity)];
  410   }
  411   return $limit;
  412 }
  413 
  414 #
  415 #  Get the routines to make the coordinates of the points
  416 #
  417 sub getVariableTypes {
  418   my $self = shift;
  419   my @make;
  420   foreach my $x (@_) {
  421     my $type = $self->{context}{variables}{$x}{type};
  422     if ($type->{name} eq 'Number') {
  423       push @make,($type->{length} == 1)? 'Value::Formula::number': $self->Package("Complex");
  424     } else {
  425       push @make, $self->Package($type->{name});
  426     }
  427   }
  428   return @make;
  429 }
  430 
  431 #
  432 #  Fake object for making reals (rather than use overhead of Value::Real)
  433 #
  434 sub Value::Formula::number::make {shift; shift; shift->value}
  435 
  436 #
  437 #  Find adaptive parameters, if any
  438 #
  439 sub AdaptParameters {
  440   my $l = shift; my $r = shift;
  441   my @params = @_; my $d = scalar(@params); my $D;
  442   return 0 if $d == 0; return 0 unless $l->usesOneOf(@params);
  443   $l->Error("Adaptive parameters can only be used for real-valued formulas")
  444     unless $l->{tree}->isRealNumber;
  445 
  446   #
  447   #  Try up to three times (the random points might not work the first time)
  448   #
  449   foreach my $attempt (1..3) {
  450     #
  451     #  Get coefficient matrix of adaptive parameters
  452     #  and value vector for linear system
  453     #
  454     my ($p,$v) = $l->createRandomPoints($d);
  455     my @P = (0) x $d; my ($f,$F) = ($l->{f},$r->{f});
  456     my @A = (); my @b = ();
  457     foreach my $i (0..$d-1) {
  458       my @a = (); my @p = @{$p->[$i]};
  459       foreach my $j (0..$d-1) {
  460   $P[$j] = 1;
  461   my $y = eval {&$f(@p,@P)};
  462   $l->Error(["Can't evaluate correct answer at adapted point (%s)",join(",",@$p,@P)])
  463          unless defined $y;
  464   push(@a,($y-$v->[$i])->value);
  465   $P[$j] = 0;
  466       }
  467       my $y = eval {&$F(@p,@P)}; return unless defined $y;
  468       push @A, [@a]; push @b, [($y-$v->[$i])->value];
  469     }
  470     #
  471     #  Use MatrixReal1.pm to solve system of linear equations
  472     #
  473     my $M = MatrixReal1->new($d,$d); $M->[0] = \@A;
  474     my $B = MatrixReal1->new($d,1);  $B->[0] = \@b;
  475     ($M,$B) = $M->normalize($B);
  476     $M = $M->decompose_LR;
  477     if (abs($M->det_LR) > 1E-6) {
  478       if (($D,$B,$M) = $M->solve_LR($B)) {
  479   if ($D == 0) {
  480     #
  481     #  Get parameter values and recompute the points using them
  482     #
  483     my @a; my $i = 0; my $max = $l->getFlag('max_adapt',1E8);
  484     foreach my $row (@{$B->[0]}) {
  485       if (abs($row->[0]) > $max) {
  486         $max = Value::makeValue($max); $row->[0] = Value::makeValue($row->[0]);
  487         $l->Error(["Constant of integration is too large: %s\n(maximum allowed is %s)",
  488        $row->[0]->string,$max->string]) if $params[$i] eq 'C0' or $params[$i] eq 'n00';
  489         $l->Error(["Adaptive constant is too large: %s = %s\n(maximum allowed is %s)",
  490        $params[$i],$row->[0]->string,$max->string]);
  491       }
  492       push @a, $row->[0]; $i++;
  493     }
  494     my $context = $l->context;
  495     foreach my $i (0..$#a) {$context->{variables}{$params[$i]}{value} = $a[$i]}
  496     $l->{parameters} = [@a];
  497     $l->createAdaptedValues;
  498     return 1;
  499   }
  500       }
  501     }
  502   }
  503   $l->Error("Can't solve for adaptive parameters");
  504 }
  505 
  506 ##
  507 ##  debugging routine
  508 ##
  509 #sub main::Format {
  510 #  my $v = scalar(@_) > 1? [@_]: shift;
  511 #  $v = [%{$v}] if ref($v) eq 'HASH';
  512 #  return $v unless ref($v) eq 'ARRAY';
  513 #  my @V; foreach my $x (@{$v}) {push @V, main::Format($x)}
  514 #  return '['.join(",",@V).']';
  515 #}
  516 
  517 #
  518 #  Random number generator  (replaced by Value::WeBWorK.pm)
  519 #
  520 sub seedRandom {srand}
  521 sub getRandom {
  522   my $self = shift;
  523   my ($m,$M,$n) = @_; $n = 1 unless $n;
  524   return $m + $n*int(rand()*(int(($M-$m)/$n)+1));
  525 }
  526 
  527 ############################################
  528 #
  529 #  Check if the value of a formula is constant
  530 #    (could use shift->{tree}{isConstant}, but I don't trust it)
  531 #
  532 sub isConstant {
  533   my @vars = (%{shift->{variables}});
  534   return scalar(@vars) == 0;
  535 }
  536 
  537 #
  538 #  Check if the Formula includes one of the named variables
  539 #
  540 sub usesOneOf {
  541   my $self = shift;
  542   foreach my $x (@_) {return 1 if $self->{variables}{$x}}
  543   return 0;
  544 }
  545 
  546 ###########################################################################
  547 
  548 1;

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9