[system] / trunk / pg / lib / Value / Formula.pm Repository:
ViewVC logotype

Annotation of /trunk/pg/lib/Value/Formula.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2687 - (view) (download) (as text)

1 : sh002i 2558 ###########################################################################
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 : dpvc 2625 '.' => \&dot,
19 : sh002i 2558 'x' => \&cross,
20 : dpvc 2622 '<=>' => \&compare,
21 : dpvc 2625 'cmp' => \&Value::cmp,
22 : sh002i 2558 '~' => 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 : dpvc 2606 '""' => \&Value::stringify;
33 : sh002i 2558
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 : dpvc 2622 ############################################
51 : sh002i 2558 #
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 : dpvc 2678 my $formula = $pkg->blank; my $parser = $formula->{context}{parser};
67 : sh002i 2558 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 : dpvc 2678 $l = $parser->{Value}->new($formula,$l) unless ref($l) =~ m/^Parser::/;
81 :     $r = $parser->{Value}->new($formula,$r) unless ref($r) =~ m/^Parser::/;
82 : sh002i 2558 $bop = 'U' if $bop eq '+' &&
83 :     ($l->type =~ m/Interval|Union/ || $r->type =~ m/Interval|Union/);
84 : dpvc 2678 $formula->{tree} = $parser->{BOP}->new($formula,$bop,$l,$r);
85 : sh002i 2558 $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 : dpvc 2625 sub cross {bop(@_,'><')}
96 : sh002i 2558
97 : dpvc 2625 #
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 : dpvc 2622 ############################################
109 : sh002i 2558 #
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 : dpvc 2678 $formula->{tree} = $formula->{context}{parser}{UOP}->new($formula,'u-',$self->{tree}->copy($formula));
118 : sh002i 2558 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 : dpvc 2625 if ($l->promotePrecedence($r)) {return $r->compare($l,!$flag)}
128 : sh002i 2558 if ($flag) {my $tmp = $l; $l = $r; $r = $tmp}
129 :     Parser::Function->call('atan2',$l,$r);
130 :     }
131 :    
132 : dpvc 2622 ############################################
133 : dpvc 2621 #
134 : dpvc 2622 # 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 : dpvc 2687 # Handle parameters
156 :     #
157 :     $lvalues = $l->{test_values}
158 :     if $l->AdaptParameters($r,$self->{context}->variables->parameters);
159 :    
160 :     #
161 : dpvc 2622 # 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 : dpvc 2629 foreach $i (0..scalar(@{$lvalues})-1) {
167 : dpvc 2622 $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 : dpvc 2687 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 : dpvc 2622
185 :     my $values = []; my $v;
186 :     foreach my $p (@{$points}) {
187 : dpvc 2687 $v = eval {&$f(@{$p},@zeros)};
188 : dpvc 2622 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 : dpvc 2624 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 : dpvc 2687 my @vars = $self->{context}->variables->variables;
210 :     my @params = $self->{context}->variables->parameters;
211 : dpvc 2669 my @limits = $self->getVariableLimits(@vars);
212 : dpvc 2687 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 : dpvc 2622 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 : dpvc 2666 my (@P,@p,$v,$i); my $k = 0;
221 : dpvc 2622 while (scalar(@{$points}) < $num_points && $k < 10) {
222 : dpvc 2666 @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 : dpvc 2687 $v = eval {&$f(@P,@zeros)};
228 : dpvc 2622 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 : dpvc 2624 return ($points,$values) if defined(@_[0]);
237 : dpvc 2622 $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 : dpvc 2669 my $userlimits = $self->{limits};
247 : dpvc 2666 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 : dpvc 2669 $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 : dpvc 2622 foreach my $x (@_) {
260 : dpvc 2666 my $def = $self->{context}{variables}{$x};
261 :     my $limit = $userlimits->[$i++] || $def->{limits} || [];
262 : dpvc 2669 $limit = [$limit] if defined($limit->[0]) && ref($limit->[0]) ne 'ARRAY';
263 : dpvc 2666 push(@{$limit},$limit->[0] || $default) while (scalar(@{$limit}) < $def->{type}{length});
264 :     pop(@{$limit}) while (scalar(@{$limit}) > $def->{type}{length});
265 : dpvc 2669 push @limits, $self->addGranularity($limit,$def,$granularity,$resolution);
266 : dpvc 2622 }
267 :     return @limits;
268 :     }
269 :    
270 : dpvc 2666 #
271 :     # Add the granularity to the limit intervals
272 :     #
273 : dpvc 2669 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 : dpvc 2666 }
282 : dpvc 2669 return $limit;
283 : dpvc 2666 }
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 : dpvc 2687 #
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;
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 values and recompute the points using them
342 :     #
343 :     my @a; foreach my $r (@{$B->[0]}) {push @a, $r->[0]}
344 :     $l->{parameters} = [@a];
345 :     $l->createPointValues;
346 :     return 1;
347 :     }
348 :     }
349 :     $l->Error("Can't solve for adaptive parameters");
350 :     }
351 :    
352 : dpvc 2666 ##
353 :     ## debugging routine
354 :     ##
355 :     #sub main::Format {
356 :     # my $v = scalar(@_) > 1? [@_]: shift;
357 : dpvc 2671 # $v = [%{$v}] if ref($v) eq 'HASH';
358 : dpvc 2666 # return $v unless ref($v) eq 'ARRAY';
359 :     # my @V; foreach my $x (@{$v}) {push @V, main::Format($x)}
360 :     # return '['.join(",",@V).']';
361 :     #}
362 :    
363 :     #
364 :     # Random number generator (replaced by Value::WeBWorK.pm)
365 :     #
366 : dpvc 2622 sub seedRandom {srand}
367 :     sub getRandom {
368 :     my $self = shift;
369 :     my ($m,$M,$n) = @_; $n = 1 unless $n;
370 :     return $m + $n*int(rand()*(int(($M-$m)/$n)+1));
371 :     }
372 :    
373 :     #
374 :     # Get the value of a flag from the object itself,
375 :     # or from the context, or from the default context
376 :     # or from the given default, whichever is found first.
377 :     #
378 :     sub getFlag {
379 :     my $self = shift; my $name = shift;
380 :     return $self->{$name} if defined($self->{$name});
381 :     return $self->{context}{flags}{$name} if defined($self->{context}{flags}{$name});
382 :     return $$Value::context->{flags}{$name} if defined($$Value::context->{flags}{$name});
383 :     return shift;
384 :     }
385 :    
386 :     ############################################
387 :     #
388 : dpvc 2621 # Check if the value of a formula is constant
389 :     # (could use shift->{tree}{isConstant}, but I don't trust it)
390 :     #
391 :     sub isConstant {scalar(%{shift->{variables}}) == 0}
392 :    
393 : sh002i 2558 ###########################################################################
394 :    
395 :     1;

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9