Parent Directory
| 
Revision Log
Update new() and make() methods to accept a context as the first parameter (making it easier to create objects in a given context without having to resort to a separate call to coerce them to the given context after the fact).
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 # Call Parser to make the new Formula 15 # 16 sub new { 17 my $self = shift; 18 my $f = $self->SUPER::new(@_); 19 foreach my $id ('open','close') {$f->{$id} = $f->{tree}{$id}} 20 return $f; 21 } 22 23 # 24 # Create a new Formula with no string 25 # (we'll fill in its tree by hand) 26 # 27 sub blank {shift->SUPER::new(@_)} 28 29 # 30 # with() changes tree element as well 31 # as the formula itself. 32 # 33 sub with { 34 my $self = shift; my %hash = @_; 35 foreach my $id (keys(%hash)) { 36 $self->{tree}{$id} = $hash{$id}; 37 $self->{$id} = $hash{$id}; 38 } 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 ############################################ 55 # 56 # Create a BOP from two operands 57 # 58 # Get the context and variables from the left and right operands 59 # if they are formulas 60 # Make them into Value objects if they aren't already. 61 # Convert '+' to union for intervals or unions. 62 # Make a new BOP with the two operands. 63 # Record the variables. 64 # Evaluate the formula if it is constant. 65 # 66 sub bop { 67 my $bop = shift; 68 my ($self,$l,$r) = Value::checkOpOrder(@_); 69 my $class = ref($self) || $self; 70 my $call = $self->context->{method}{$bop}; 71 my $formula = $self->blank($self->context); 72 my $parser = $formula->{context}{parser}; 73 if (ref($r) eq $class || ref($r) eq $pkg) { 74 $formula->{context} = $r->{context}; 75 $r = $r->{tree}->copy($formula); 76 } 77 if (ref($l) eq $class || ref($l) eq $pkg) { 78 $formula->{context} = $l->{context}; 79 $l = $l->{tree}->copy($formula); 80 } 81 $l = $self->new($l) if (!ref($l) && Value::getType($formula,$l) eq "unknown"); 82 $r = $self->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|Set|Union/ || $r->type =~ m/Interval|Set|Union/); 87 $formula->{tree} = $parser->{BOP}->new($formula,$bop,$l,$r); 88 $formula->{variables} = $formula->{tree}->getVariables; 89 return $formula; 90 } 91 92 sub add {bop('+',@_)} 93 sub sub {bop('-',@_)} 94 sub mult {bop('*',@_)} 95 sub div {bop('/',@_)} 96 sub power {bop('**',@_)} 97 sub cross {bop('><',@_)} 98 99 # 100 # Make dot work for vector operands 101 # 102 sub _dot { 103 my ($l,$r,$flag) = @_; 104 if ($l->promotePrecedence($r)) {return $r->_dot($l,!$flag)} 105 return bop('.',@_) if ($l->type eq 'Vector' || $l->{isVector}) && 106 Value::isValue($r) && ($r->type eq 'Vector' || $r->{isVector}); 107 $l->SUPER::_dot($r,$flag); 108 } 109 110 sub pdot {'('.(shift->stringify).')'} 111 112 # 113 # Call the Parser::Function call function 114 # 115 sub call { 116 my $self = shift; my $name = shift; 117 Parser::Function->call($name,$self); 118 } 119 120 ############################################ 121 # 122 # Form the negation of a formula 123 # 124 sub neg { 125 my $self = shift; 126 my $formula = $self->blank($self->context); 127 $formula->{variables} = $self->{variables}; 128 $formula->{tree} = $formula->{context}{parser}{UOP}->new($formula,'u-',$self->{tree}->copy($formula)); 129 return $formula; 130 } 131 132 # 133 # Form the function atan2 function call on two operands 134 # 135 sub atan2 { 136 my ($self,$l,$r) = Value::checkOpOrder(@_); 137 Parser::Function->call('atan2',$l,$r); 138 } 139 140 sub sin {shift->call('sin',@_)} 141 sub cos {shift->call('cos',@_)} 142 sub abs {shift->call('abs',@_)} 143 sub exp {shift->call('exp',@_)} 144 sub log {shift->call('log',@_)} 145 sub sqrt {shift->call('sqrt',@_)} 146 147 sub twiddle {shift->call('conj',@_)} 148 149 ############################################ 150 # 151 # Compare two functions for equality 152 # 153 sub compare { 154 my ($l,$r) = @_; my $self = $l; 155 $r = $self->Package("Formula")->new($self->context,$r) unless Value::isFormula($r); 156 Value::Error("Functions from different contexts can't be compared") 157 unless $l->{context} == $r->{context}; 158 159 # 160 # Get the test points and evaluate the functions at those points 161 # 162 ## FIXME: Check given points for consistency 163 my $points = $l->{test_points} || $l->createRandomPoints(undef,$l->{test_at}); 164 my $lvalues = $l->{test_values} || $l->createPointValues($points,1,1); 165 my $rvalues = $r->createPointValues($points,0,1,$l->{checkUndefinedPoints}); 166 # 167 # Note: $l is bigger if $r can't be evaluated at one of the points 168 return 1 unless $rvalues; 169 170 my ($i, $cmp); 171 172 # 173 # Handle adaptive parameters: 174 # Get the tolerances, and check each adapted point relative 175 # to the ORIGINAL correct answer. (This will have to be 176 # fixed if we ever do adaptive parameters for non-real formulas) 177 # 178 if ($l->AdaptParameters($r,$self->{context}->variables->parameters)) { 179 my $avalues = $l->{test_adapt}; 180 my $tolerance = $self->getFlag('tolerance',1E-4); 181 my $isRelative = $self->getFlag('tolType','relative') eq 'relative'; 182 my $zeroLevel = $self->getFlag('zeroLevel',1E-14); 183 my $zeroLevelTol = $self->getFlag('zeroLevelTol',1E-12); 184 foreach $i (0..scalar(@{$lvalues})-1) { 185 my $tol = $tolerance; 186 my ($lv,$rv,$av) = ($lvalues->[$i]->value,$rvalues->[$i]->value,$avalues->[$i]->value); 187 if ($isRelative) { 188 if (abs($lv) <= $zeroLevel) {$tol = $zeroLevelTol} 189 else {$tol *= abs($lv)} 190 } 191 return $rv <=> $av unless abs($rv - $av) < $tol; 192 } 193 return 0; 194 } 195 196 # 197 # Look through the two lists of values to see if they are equal. 198 # If not, return the comparison of the first unequal value 199 # (not good for < and >, but OK for ==). 200 # 201 my $domainError = 0; 202 foreach $i (0..scalar(@{$lvalues})-1) { 203 if (ref($lvalues->[$i]) eq 'UNDEF' ^ ref($rvalues->[$i]) eq 'UNDEF') {$domainError = 1; next} 204 $cmp = $lvalues->[$i] <=> $rvalues->[$i]; 205 return $cmp if $cmp; 206 } 207 $l->{domainMismatch} = $domainError; # return this value 208 } 209 210 # 211 # Create the value list from a given set of test points 212 # 213 sub createPointValues { 214 my $self = shift; 215 my $points = shift || $self->{test_points} || $self->createRandomPoints; 216 my $showError = shift; my $cacheResults = shift; 217 my @vars = $self->{context}->variables->variables; 218 my @params = $self->{context}->variables->parameters; 219 my @zeros = (0) x scalar(@params); 220 my $f = $self->{f}; $f = $self->{f} = $self->perlFunction(undef,[@vars,@params]) unless $f; 221 my $checkUndef = scalar(@params) == 0 && (shift || $self->getFlag('checkUndefinedPoints',0)); 222 223 my $values = []; my $v; 224 foreach my $p (@{$points}) { 225 $v = eval {&$f(@{$p},@zeros)}; 226 if (!defined($v) && !$checkUndef) { 227 return unless $showError; 228 Value::Error("Can't evaluate formula on test point (%s)",join(',',@{$p})); 229 } 230 push @{$values}, (defined($v)? Value::makeValue($v,context=>$self->context): $UNDEF); 231 } 232 if ($cacheResults) { 233 $self->{test_points} = $points; 234 $self->{test_values} = $values; 235 } 236 return $values; 237 } 238 239 # 240 # Create the adapted value list for the test points 241 # 242 sub createAdaptedValues { 243 my $self = shift; 244 my $points = shift || $self->{test_points} || $self->createRandomPoints; 245 my $showError = shift; 246 my @vars = $self->{context}->variables->variables; 247 my @params = $self->{context}->variables->parameters; 248 my $f = $self->{f}; $f = $self->{f} = $self->perlFunction(undef,[@vars,@params]) unless $f; 249 250 my $values = []; my $v; 251 my @adapt = @{$self->{parameters}}; 252 foreach my $p (@{$points}) { 253 $v = eval {&$f(@{$p},@adapt)}; 254 if (!defined($v)) { 255 return unless $showError; 256 Value::Error("Can't evaluate formula on test point (%s) with parameters (%s)", 257 join(',',@{$p}),join(',',@adapt)); 258 } 259 push @{$values}, Value::makeValue($v,context=>$self->context); 260 } 261 $self->{test_adapt} = $values; 262 } 263 264 # 265 # Create a list of random points, making sure that the function 266 # is defined at the given points. Error if we can't find enough. 267 # 268 sub createRandomPoints { 269 my $self = shift; my $context = $self->context; 270 my ($num_points,$include) = @_; my $cacheResults = !defined($num_points); 271 $num_points = int($self->getFlag('num_points',5)) unless defined($num_points); 272 $num_points = 1 if $num_points < 1; 273 274 my @vars = $context->variables->variables; 275 my @params = $context->variables->parameters; 276 my @limits = $self->getVariableLimits(@vars); 277 my @make = $self->getVariableTypes(@vars); 278 my @zeros = (0) x scalar(@params); 279 my $f = $self->{f}; $f = $self->{f} = $self->perlFunction(undef,[@vars,@params]) unless $f; 280 my $seedRandom = $context->flag('random_seed')? 'PGseedRandom' : 'seedRandom'; 281 my $getRandom = $context->flag('random_seed')? 'PGgetRandom' : 'getRandom'; 282 my $checkUndef = scalar(@params) == 0 && $self->getFlag('checkUndefinedPoints',0); 283 my $max_undef = $self->getFlag('max_undefined',$num_points); 284 285 $self->$seedRandom; 286 my $points = []; my $values = []; my $num_undef = 0; 287 if ($include) { 288 push(@{$points},@{$include}); 289 push(@{$values},@{$self->createPointValues($include,1,$cacheResults,$self->{checkundefinedPoints})}); 290 } 291 my (@P,@p,$v,$i); my $k = 0; 292 while (scalar(@{$points}) < $num_points+$num_undef && $k < 10) { 293 @P = (); $i = 0; 294 foreach my $limit (@limits) { 295 @p = (); foreach my $I (@{$limit}) 296 {push @p, $self->Package("Real")->make($context,$self->$getRandom(@{$I}))} 297 push @P, $make[$i++]->make($context,@p); 298 } 299 $v = eval {&$f(@P,@zeros)}; 300 if (!defined($v)) { 301 if ($checkUndef && $num_undef < $max_undef) { 302 push @{$points}, [@P]; 303 push @{$values}, $UNDEF; 304 $num_undef++; 305 } 306 $k++; 307 } else { 308 push @{$points}, [@P]; 309 push @{$values}, Value::makeValue($v,context=>$context); 310 $k = 0; # reset count when we find a point 311 } 312 } 313 314 if ($k) { 315 my $error = "Can't generate enough valid points for comparison"; 316 $error .= ':<div style="margin-left:1em">'.($context->{error}{message} || $@).'</div>' 317 if ($self->getFlag('showTestPointErrors')); 318 $error =~ s/ (in \S+ )?at line \d+.*//s; 319 Value::Error($error); 320 } 321 322 return ($points,$values) unless $cacheResults; 323 $self->{test_values} = $values; 324 $self->{test_points} = $points; 325 } 326 327 # 328 # Get the array of variable limits 329 # 330 sub getVariableLimits { 331 my $self = shift; 332 my $userlimits = $self->{limits}; 333 if (defined($userlimits)) { 334 $userlimits = [[[-$userlimits,$userlimits]]] unless ref($userlimits) eq 'ARRAY'; 335 $userlimits = [$userlimits] unless ref($userlimits->[0]) eq 'ARRAY'; 336 $userlimits = [$userlimits] if scalar(@_) == 1 && ref($userlimits->[0][0]) ne 'ARRAY'; 337 foreach my $I (@{$userlimits}) {$I = [$I] unless ref($I->[0]) eq 'ARRAY'}; 338 } 339 $userlimits = [] unless $userlimits; my @limits; 340 my $default; $default = $userlimits->[0][0] if defined($userlimits->[0]); 341 $default = $default || $self->getFlag('limits',[-2,2]); 342 my $granularity = $self->getFlag('granularity',1000); 343 my $resolution = $self->getFlag('resolution'); 344 my $i = 0; 345 foreach my $x (@_) { 346 my $def = $self->{context}{variables}{$x}; 347 my $limit = $userlimits->[$i++] || $def->{limits} || []; 348 $limit = [$limit] if defined($limit->[0]) && ref($limit->[0]) ne 'ARRAY'; 349 push(@{$limit},$limit->[0] || $default) while (scalar(@{$limit}) < $def->{type}{length}); 350 pop(@{$limit}) while (scalar(@{$limit}) > $def->{type}{length}); 351 push @limits, $self->addGranularity($limit,$def,$granularity,$resolution); 352 } 353 return @limits; 354 } 355 356 # 357 # Add the granularity to the limit intervals 358 # 359 sub addGranularity { 360 my $self = shift; my $limit = shift; my $def = shift; 361 my $granularity = shift; my $resolution = shift; 362 $granularity = $def->{granularity} || $granularity; 363 $resolution = $def->{resolution} || $resolution; 364 foreach my $I (@{$limit}) { 365 my ($a,$b,$n) = @{$I}; $b = -$a unless defined $b; 366 $I = [$a,$b,($n || $resolution || abs($b-$a)/$granularity)]; 367 } 368 return $limit; 369 } 370 371 # 372 # Get the routines to make the coordinates of the points 373 # 374 sub getVariableTypes { 375 my $self = shift; 376 my @make; 377 foreach my $x (@_) { 378 my $type = $self->{context}{variables}{$x}{type}; 379 if ($type->{name} eq 'Number') { 380 push @make,($type->{length} == 1)? 'Value::Formula::number': $self->Package("Complex"); 381 } else { 382 push @make, $self->Package($type->{name}); 383 } 384 } 385 return @make; 386 } 387 388 # 389 # Fake object for making reals (rather than use overhead of Value::Real) 390 # 391 sub Value::Formula::number::make {shift; shift; shift->value} 392 393 # 394 # Find adaptive parameters, if any 395 # 396 sub AdaptParameters { 397 my $l = shift; my $r = shift; 398 my @params = @_; my $d = scalar(@params); 399 return 0 if $d == 0; return 0 unless $l->usesOneOf(@params); 400 $l->Error("Adaptive parameters can only be used for real-valued functions") 401 unless $l->{tree}->isRealNumber; 402 # 403 # Get coefficient matrix of adaptive parameters 404 # and value vector for linear system 405 # 406 my ($p,$v) = $l->createRandomPoints($d); 407 my @P = (0) x $d; my ($f,$F) = ($l->{f},$r->{f}); 408 my @A = (); my @b = (); 409 foreach my $i (0..$d-1) { 410 my @a = (); my @p = @{$p->[$i]}; 411 foreach my $j (0..$d-1) { 412 $P[$j] = 1; push(@a,(&$f(@p,@P)-$v->[$i])->value); 413 $P[$j] = 0; 414 } 415 push @A, [@a]; push @b, [(&$F(@p,@P)-$v->[$i])->value]; 416 } 417 # 418 # Use MatrixReal1.pm to solve system of linear equations 419 # 420 my $M = MatrixReal1->new($d,$d); $M->[0] = \@A; 421 my $B = MatrixReal1->new($d,1); $B->[0] = \@b; 422 ($M,$B) = $M->normalize($B); 423 $M = $M->decompose_LR; 424 if (($d,$B,$M) = $M->solve_LR($B)) { 425 if ($d == 0) { 426 # 427 # Get parameter values and recompute the points using them 428 # 429 my @a; my $i = 0; my $max = $l->getFlag('max_adapt',1E8); 430 foreach my $row (@{$B->[0]}) { 431 if (abs($row->[0]) > $max) { 432 $max = Value::makeValue($max); 433 $l->Error("Constant of integration is too large: %s\n(maximum allowed is %s)", 434 $row->[0]->string,$max->string) if ($params[$i] eq 'C0'); 435 $l->Error("Adaptive constant is too large: %s = %s\n(maximum allowed is %s)", 436 $params[$i],$row->[0]->string,$max->string); 437 } 438 push @a, $row->[0]; $i++; 439 } 440 $l->{parameters} = [@a]; 441 $l->createAdaptedValues; 442 return 1; 443 } 444 } 445 $l->Error("Can't solve for adaptive parameters"); 446 } 447 448 sub usesOneOf { 449 my $self = shift; 450 foreach my $x (@_) {return 1 if $self->{variables}{$x}} 451 return 0; 452 } 453 454 ## 455 ## debugging routine 456 ## 457 #sub main::Format { 458 # my $v = scalar(@_) > 1? [@_]: shift; 459 # $v = [%{$v}] if ref($v) eq 'HASH'; 460 # return $v unless ref($v) eq 'ARRAY'; 461 # my @V; foreach my $x (@{$v}) {push @V, main::Format($x)} 462 # return '['.join(",",@V).']'; 463 #} 464 465 # 466 # Random number generator (replaced by Value::WeBWorK.pm) 467 # 468 sub seedRandom {srand} 469 sub getRandom { 470 my $self = shift; 471 my ($m,$M,$n) = @_; $n = 1 unless $n; 472 return $m + $n*int(rand()*(int(($M-$m)/$n)+1)); 473 } 474 475 ############################################ 476 # 477 # Check if the value of a formula is constant 478 # (could use shift->{tree}{isConstant}, but I don't trust it) 479 # 480 sub isConstant { 481 my @vars = (%{shift->{variables}}); 482 return scalar(@vars) == 0; 483 } 484 485 ############################################ 486 # 487 # Provide output formats 488 # 489 sub stringify { 490 my $self = shift; 491 return $self->TeX if $self->getFlag('StringifyAsTeX'); 492 $self->string; 493 } 494 495 ########################################################################### 496 497 1;
| aubreyja at gmail dot com | ViewVC Help |
| Powered by ViewVC 1.0.9 |