[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 3177 - (download) (as text) (annotate)
Fri Feb 18 00:04:09 2005 UTC (14 years, 9 months ago) by dpvc
File size: 12853 byte(s)
Missed one situation when making the previous update.  Fixed it.

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

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9