Parent Directory
| 
Revision Log
Added isZero and isOne checks for Parser::Value objects (i.e., for constants within formulas). These now correctly handle vector and matrices, in particular. The isOne and isZero checks are used in the reduce() method to simplify formulas.
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 ($l,$r,$flag,$bop) = @_; 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 |