[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 2688 - (download) (as text) (annotate)
Fri Aug 27 19:42:59 2004 UTC (15 years, 3 months ago) by dpvc
File size: 12780 byte(s)
Added in functMaxConstantOfIntegration support in the adaptive
parameter answer checking.  Also added a flag to the formula answer
checker that allows checking a formula up to a constant (for
integration problems).  This really just adds a new parameter to the
context and adds that to the formula, so there is nothing deep about
this.

	ANS(Formula("2x")->cmp(upToConstant=>1));

Finally, don't look for adaptive parameters if they aren't actually
used in the professor's formula (even if they are defined).

    1 ###########################################################################
    2 #
    3 #  Implements the Formula class.
    4 #
    5 package Value::Formula;
    6 my $pkg = 'Value::Formula';
    7 
    8 use strict;
    9 use vars qw(@ISA);
   10 @ISA = qw(Parser Value);
   11 
   12 use overload
   13        '+'    => \&add,
   14        '-'    => \&sub,
   15        '*'    => \&mult,
   16        '/'    => \&div,
   17        '**'   => \&power,
   18        '.'    => \&dot,
   19        'x'    => \&cross,
   20        '<=>'  => \&compare,
   21        'cmp'  => \&Value::cmp,
   22        '~'    => sub {Parser::Function->call('conj',$_[0])},
   23        'neg'  => sub {$_[0]->neg},
   24        'sin'  => sub {Parser::Function->call('sin',$_[0])},
   25        'cos'  => sub {Parser::Function->call('cos',$_[0])},
   26        'exp'  => sub {Parser::Function->call('exp',$_[0])},
   27        'abs'  => sub {Parser::Function->call('abs',$_[0])},
   28        'log'  => sub {Parser::Function->call('log',$_[0])},
   29        'sqrt' => sub {Parser::Function->call('sqrt',$_[0])},
   30       'atan2' => \&atan2,
   31    'nomethod' => \&Value::nomethod,
   32          '""' => \&Value::stringify;
   33 
   34 #
   35 #  Call Parser to make the new item
   36 #
   37 sub new {shift; $pkg->SUPER::new(@_)}
   38 
   39 #
   40 #  Create the new parser with no string
   41 #    (we'll fill in its tree by hand)
   42 #
   43 sub blank {$pkg->SUPER::new('')}
   44 
   45 #
   46 #  Get the type from the tree
   47 #
   48 sub typeRef {(shift)->{tree}->typeRef}
   49 
   50 ############################################
   51 #
   52 #  Create a BOP from two operands
   53 #
   54 #  Get the context and variables from the left and right operands
   55 #    if they are formulas
   56 #  Make them into Value objects if they aren't already.
   57 #  Convert '+' to union for intervals or unions.
   58 #  Make a new BOP with the two operands.
   59 #  Record the variables.
   60 #  Evaluate the formula if it is constant.
   61 #
   62 sub bop {
   63   my ($l,$r,$flag,$bop) = @_;
   64   if ($l->promotePrecedence($r)) {return $r->add($l,!$flag)}
   65   if ($flag) {my $tmp = $l; $l = $r; $r = $tmp}
   66   my $formula = $pkg->blank; my $parser = $formula->{context}{parser};
   67   my $vars = {};
   68   if (ref($r) eq $pkg) {
   69     $formula->{context} = $r->{context};
   70     $vars = {%{$vars},%{$r->{variables}}};
   71     $r = $r->{tree}->copy($formula);
   72   }
   73   if (ref($l) eq $pkg) {
   74     $formula->{context} = $l->{context};
   75     $vars = {%{$vars},%{$l->{variables}}};
   76     $l = $l->{tree}->copy($formula);
   77   }
   78   $l = $pkg->new($l) if (!ref($l) && Value::getType($formula,$l) eq "unknown");
   79   $r = $pkg->new($r) if (!ref($r) && Value::getType($formula,$r) eq "unknown");
   80   $l = $parser->{Value}->new($formula,$l) unless ref($l) =~ m/^Parser::/;
   81   $r = $parser->{Value}->new($formula,$r) unless ref($r) =~ m/^Parser::/;
   82   $bop = 'U' if $bop eq '+' &&
   83     ($l->type =~ m/Interval|Union/ || $r->type =~ m/Interval|Union/);
   84   $formula->{tree} = $parser->{BOP}->new($formula,$bop,$l,$r);
   85   $formula->{variables} = {%{$vars}};
   86   return $formula->eval if scalar(%{$vars}) == 0;
   87   return $formula;
   88 }
   89 
   90 sub add   {bop(@_,'+')}
   91 sub sub   {bop(@_,'-')}
   92 sub mult  {bop(@_,'*')}
   93 sub div   {bop(@_,'/')}
   94 sub power {bop(@_,'**')}
   95 sub cross {bop(@_,'><')}
   96 
   97 #
   98 #  Make dot work for vector operands
   99 #
  100 sub dot   {
  101   my ($l,$r,$flag) = @_;
  102   if ($l->promotePrecedence($r)) {return $r->compare($l,!$flag)}
  103   return bop(@_,'.') if $l->type eq 'Vector' &&
  104      Value::isValue($r) && $r->type eq 'Vector';
  105   Value::_dot(@_);
  106 }
  107 
  108 ############################################
  109 #
  110 #  Form the negation of a formula
  111 #
  112 sub neg {
  113   my $self = shift;
  114   my $formula = $self->blank;
  115   $formula->{context} = $self->{context};
  116   $formula->{variables} = $self->{variables};
  117   $formula->{tree} = $formula->{context}{parser}{UOP}->new($formula,'u-',$self->{tree}->copy($formula));
  118   return $formula->eval if scalar(%{$formula->{variables}}) == 0;
  119   return $formula;
  120 }
  121 
  122 #
  123 #  Form the function atan2 function call on two operands
  124 #
  125 sub atan2 {
  126   my ($l,$r,$flag) = @_;
  127   if ($l->promotePrecedence($r)) {return $r->compare($l,!$flag)}
  128   if ($flag) {my $tmp = $l; $l = $r; $r = $tmp}
  129   Parser::Function->call('atan2',$l,$r);
  130 }
  131 
  132 ############################################
  133 #
  134 #  Compare two functions for equality
  135 #
  136 sub compare {
  137   my ($l,$r,$flag) = @_; my $self = $l;
  138   if ($l->promotePrecedence($r)) {return $r->compare($l,!$flag)}
  139   $r = Value::Formula->new($r) unless Value::isFormula($r);
  140   Value::Error("Functions from different contexts can't be compared")
  141     unless $l->{context} == $r->{context};
  142 
  143   #
  144   #  Get the test points and evaluate the functions at those points
  145   #
  146   ##  FIXME: Check given points for consistency
  147   my $points = $l->{test_points} || $r->{test_points} || $l->createRandomPoints;
  148   my $lvalues = $l->{test_values} || $l->createPointValues($points,1);
  149   my $rvalues = $r->createPointValues($points);
  150   #
  151   # Note: $l is bigger if $r can't be evaluated at one of the points
  152   return 1 unless $rvalues;
  153 
  154   #
  155   #  Handle parameters
  156   #
  157   $lvalues = $l->{test_values}
  158     if $l->AdaptParameters($r,$self->{context}->variables->parameters);
  159 
  160   #
  161   #  Look through the two lists to see if they are equal.
  162   #  If not, return the comparison of the first unequal value
  163   #    (not good for < and >, but OK for ==).
  164   #
  165   my ($i, $cmp);
  166   foreach $i (0..scalar(@{$lvalues})-1) {
  167     $cmp = $lvalues->[$i] <=> $rvalues->[$i];
  168     return $cmp if $cmp;
  169   }
  170   return 0;
  171 }
  172 
  173 #
  174 #  Create the value list from a given set of test points
  175 #
  176 sub createPointValues {
  177   my $self = shift;
  178   my $points = shift || $self->{test_points} || $self->createRandomPoints;
  179   my $showError = shift;
  180   my @vars   = $self->{context}->variables->variables;
  181   my @params = $self->{context}->variables->parameters;
  182   my @zeros  = @{$self->{parameters} || [split('',"0" x scalar(@params))]};
  183   my $f = $self->{f}; $f = $self->{f} = $self->perlFunction(undef,[@vars,@params]) unless $f;
  184 
  185   my $values = []; my $v;
  186   foreach my $p (@{$points}) {
  187     $v = eval {&$f(@{$p},@zeros)};
  188     if (!defined($v)) {
  189       return unless $showError;
  190       Value::Error("Can't evaluate formula on test point (".join(',',@{$p}).")");
  191     }
  192     push @{$values}, Value::makeValue($v);
  193   }
  194 
  195   $self->{test_points} = $points;
  196   $self->{test_values} = $values;
  197 }
  198 
  199 #
  200 #  Create a list of random points, making sure that the function
  201 #  is defined at the given points.  Error if we can't find enough.
  202 #
  203 sub createRandomPoints {
  204   my $self = shift;
  205   my $num_points = @_[0];
  206   $num_points = int($self->getFlag('num_points',5)) unless defined($num_points);
  207   $num_points = 1 if $num_points < 1;
  208 
  209   my @vars   = $self->{context}->variables->variables;
  210   my @params = $self->{context}->variables->parameters;
  211   my @limits = $self->getVariableLimits(@vars);
  212   my @make   = $self->getVariableTypes(@vars);
  213   my @zeros  = split('',"0" x scalar(@params));
  214   my $f = $self->{f}; $f = $self->{f} = $self->perlFunction(undef,[@vars,@params]) unless $f;
  215   my $seedRandom = $self->{context}->flag('random_seed')? 'PGseedRandom' : 'seedRandom';
  216   my $getRandom  = $self->{context}->flag('random_seed')? 'PGgetRandom'  : 'getRandom';
  217 
  218   $self->$seedRandom;
  219   my $points = []; my $values = [];
  220   my (@P,@p,$v,$i); my $k = 0;
  221   while (scalar(@{$points}) < $num_points && $k < 10) {
  222     @P = (); $i = 0;
  223     foreach my $limit (@limits) {
  224       @p = (); foreach my $I (@{$limit}) {push @p, $self->$getRandom(@{$I})}
  225       push @P, $make[$i++]->make(@p);
  226     }
  227     $v = eval {&$f(@P,@zeros)};
  228     if (!defined($v)) {$k++} else {
  229       push @{$points}, [@P];
  230       push @{$values}, Value::makeValue($v);
  231       $k = 0; # reset count when we find a point
  232     }
  233   }
  234 
  235   Value::Error("Can't generate enough valid points for comparison") if $k;
  236   return ($points,$values) if defined(@_[0]);
  237   $self->{test_values} = $values;
  238   $self->{test_points} = $points;
  239 }
  240 
  241 #
  242 #  Get the array of variable limits
  243 #
  244 sub getVariableLimits {
  245   my $self = shift;
  246   my $userlimits = $self->{limits};
  247   if (defined($userlimits)) {
  248     $userlimits = [[[-$userlimits,$userlimits]]] unless ref($userlimits) eq 'ARRAY';
  249     $userlimits = [$userlimits] unless ref($userlimits->[0]) eq 'ARRAY';
  250     $userlimits = [$userlimits] if scalar(@_) == 1 && ref($userlimits->[0][0]) ne 'ARRAY';
  251     foreach my $I (@{$userlimits}) {$I = [$I] unless ref($I->[0]) eq 'ARRAY'};
  252   }
  253   $userlimits = [] unless $userlimits; my @limits;
  254   my $default;  $default = $userlimits->[0][0] if defined($userlimits->[0]);
  255   my $default = $default || $self->{context}{flags}{limits} || [-2,2];
  256   my $granularity = $self->getFlag('granularity',1000);
  257   my $resolution = $self->getFlag('resolution');
  258   my $i = 0;
  259   foreach my $x (@_) {
  260     my $def = $self->{context}{variables}{$x};
  261     my $limit = $userlimits->[$i++] || $def->{limits} || [];
  262     $limit = [$limit] if defined($limit->[0]) && ref($limit->[0]) ne 'ARRAY';
  263     push(@{$limit},$limit->[0] || $default) while (scalar(@{$limit}) < $def->{type}{length});
  264     pop(@{$limit}) while (scalar(@{$limit}) > $def->{type}{length});
  265     push @limits, $self->addGranularity($limit,$def,$granularity,$resolution);
  266   }
  267   return @limits;
  268 }
  269 
  270 #
  271 #  Add the granularity to the limit intervals
  272 #
  273 sub addGranularity {
  274   my $self = shift; my $limit = shift; my $def = shift;
  275   my $granularity = shift; my $resolution = shift;
  276   $granularity = $def->{granularity} || $granularity;
  277   $resolution = $def->{resolution} || $resolution;
  278   foreach my $I (@{$limit}) {
  279     my ($a,$b,$n) = @{$I}; $b = -$a unless defined $b;
  280     $I = [$a,$b,($n || $resolution || abs($b-$a)/$granularity)];
  281   }
  282   return $limit;
  283 }
  284 
  285 #
  286 #  Get the routines to make the coordinates of the points
  287 #
  288 sub getVariableTypes {
  289   my $self = shift;
  290   my @make;
  291   foreach my $x (@_) {
  292     my $type = $self->{context}{variables}{$x}{type};
  293     if ($type->{name} eq 'Number') {
  294       push @make,($type->{length} == 1)? 'Value::Formula::number': 'Value::Complex';
  295     } else {
  296       push @make, "Value::$type->{name}";
  297     }
  298   }
  299   return @make;
  300 }
  301 
  302 #
  303 #  Fake object for making reals (rather than use overhead of Value::Real)
  304 #
  305 sub Value::Formula::number::make {shift; shift}
  306 
  307 #
  308 #  Find adaptive parameters, if any
  309 #
  310 sub AdaptParameters {
  311   my $l = shift; my $r = shift;
  312   my @params = @_; my $d = scalar(@params);
  313   return 0 if $d == 0; return 0 unless $l->usesOneOf(@params);
  314   $l->Error("Adaptive parameters can only be used for real-valued functions")
  315     unless $l->{tree}->isRealNumber;
  316   #
  317   #  Get coefficient matrix of adaptive parameters
  318   #  and value vector for linear system
  319   #
  320   my ($p,$v) = $l->createRandomPoints($d,1);
  321   my @P = split('',"0" x $d); my ($f,$F) = ($l->{f},$r->{f});
  322   my @A = (); my @b = ();
  323   foreach my $i (0..$d-1) {
  324     my @a = (); my @p = @{$p->[$i]};
  325     foreach my $j (0..$d-1) {
  326       $P[$j] = 1; push(@a,&$f(@p,@P)-$v->[$i]);
  327       $P[$j] = 0;
  328     }
  329     push @A, [@a]; push @b, [&$F(@p,@P)-$v->[$i]];
  330   }
  331   #
  332   #  Use MatrixReal1.pm to solve system of linear equations
  333   #
  334   my $M = MatrixReal1->new($d,$d); $M->[0] = \@A;
  335   my $B = MatrixReal1->new($d,1);  $B->[0] = \@b;
  336   ($M,$B) = $M->normalize($B);
  337   $M = $M->decompose_LR;
  338   if (($d,$B,$M) = $M->solve_LR($B)) {
  339     if ($d == 0) {
  340       #
  341       #  Get parameter values and recompute the points using them
  342       #
  343       my @a; my $i = 0; my $max = $l->getFlag('max_adapt',1E8);
  344       foreach my $row (@{$B->[0]}) {
  345   if (abs($row->[0]) > $max) {
  346     $l->Error("Constant of integration is too large: $row->[0]")
  347       if ($params[$i] eq 'C0');
  348     $l->Error("Adaptive constant is too large: $params[$i] = $row->[0]");
  349   }
  350   push @a, $row->[0]; $i++;
  351       }
  352       $l->{parameters} = [@a];
  353       $l->createPointValues;
  354       return 1;
  355     }
  356   }
  357   $l->Error("Can't solve for adaptive parameters");
  358 }
  359 
  360 sub usesOneOf {
  361   my $self = shift;
  362   foreach my $x (@_) {return 1 if $self->{variables}{$x}}
  363   return 0;
  364 }
  365 
  366 ##
  367 ##  debugging routine
  368 ##
  369 #sub main::Format {
  370 #  my $v = scalar(@_) > 1? [@_]: shift;
  371 #  $v = [%{$v}] if ref($v) eq 'HASH';
  372 #  return $v unless ref($v) eq 'ARRAY';
  373 #  my @V; foreach my $x (@{$v}) {push @V, main::Format($x)}
  374 #  return '['.join(",",@V).']';
  375 #}
  376 
  377 #
  378 #  Random number generator  (replaced by Value::WeBWorK.pm)
  379 #
  380 sub seedRandom {srand}
  381 sub getRandom {
  382   my $self = shift;
  383   my ($m,$M,$n) = @_; $n = 1 unless $n;
  384   return $m + $n*int(rand()*(int(($M-$m)/$n)+1));
  385 }
  386 
  387 #
  388 #  Get the value of a flag from the object itself,
  389 #  or from the context, or from the default context
  390 #  or from the given default, whichever is found first.
  391 #
  392 sub getFlag {
  393   my $self = shift; my $name = shift;
  394   return $self->{$name} if defined($self->{$name});
  395   return $self->{context}{flags}{$name} if defined($self->{context}{flags}{$name});
  396   return $$Value::context->{flags}{$name} if defined($$Value::context->{flags}{$name});
  397   return shift;
  398 }
  399 
  400 ############################################
  401 #
  402 #  Check if the value of a formula is constant
  403 #    (could use shift->{tree}{isConstant}, but I don't trust it)
  404 #
  405 sub isConstant {scalar(%{shift->{variables}}) == 0}
  406 
  407 ###########################################################################
  408 
  409 1;

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9