[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 5517 - (download) (as text) (annotate)
Thu Sep 20 18:02:24 2007 UTC (12 years, 2 months ago) by dpvc
File size: 16513 byte(s)
Make adaptive parameters try again if they fail to be found the first
time (Gavin has found cases where the random points cause the
coefficient matrix to be singular, so a second try should resolve the
problem).

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

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9