[system] / trunk / pg / macros / parserImplicitEquation.pl Repository: Repository Listing bbplugincoursesdistsnplrochestersystemwww # View of /trunk/pg/macros/parserImplicitEquation.pl

Mon Jun 11 18:16:40 2007 UTC (12 years, 8 months ago) by gage
File size: 13141 byte(s)
```Fixing docementation so that it can be read from the web.
```

```    1 loadMacros("Parser.pl");
2
3 sub _parserImplicitEquation_init {}; # don't reload this file
4
6
7 ######################################################################
8 #
9 #  This is a Parser class that implements an answer checker for
10 #  implicitly defined equations.  The checker looks for the zeros of
11 #  the equation and tests that the student and professor equations
12 #  both have the same solutions.
13 #
14 #  This type of check is very subtle, and there are important issues
15 #  that you may have to take into account.  The solutions to the
16 #  equations are found numerically, and so they will not be exact;
17 #  that means that there are tolerances that may need to be adjusted
18 #  for your particular equation.  Also, it is always possible for the
19 #  student to represent the function in a form that will exceed those
20 #  tolerances, and so be marked as incorrect.  The answer checker
21 #  attempts to set the parameters based on the values of the functions
22 #  involved, but this may not work perfectly.
23 #
24 #  The method used to locate the solutions of A=B is by finding zeros
25 #  of A-B, and it requires this function to take on both positive and
26 #  negative values (that is, it can only find transverse intersections
27 #  of the surface A-B=0 and the plane at height 0).  For example, even
28 #  though the solutions of (x^2+y^1-1)^2=0 form a circle, the
29 #  algorithm will fail to find any solutions for this equation.
30 #
31 #  In order to locate the zeros, you may need to change the limits so
32 #  that they include regions where the function is both positive an
33 #  negative (see below).  The algorithm will avoid discontinuities, so
34 #  you can specify things like x-y=1/(x+y) rather than x^2-y^2=1.
35 #
36 #  Because the solutions are found using a random search, it is
37 #  possible the randomly chosen starting points don't locate enough
38 #  zeros for the function, in which case the check will fail.  This
39 #  can happen for both the professor's function and the student's,
40 #  since zeros are found for both.  This means that a correct answer
41 #  can sometimes be marked incorrect depending on the random points
42 #  chosen initially.  These points also affect the values selected for
43 #  the tolerances used to determine when a function's value is zero,
44 #  and so can affect whether the student's function is marked as
45 #  correct or not.
46 #
47 #  If an equation has several components or branches, it is possible
48 #  that the random location of solutions will not find zeros on some
49 #  of the branches, and so might incorrectly mark as correct an
50 #  equation that only is zero on one of the components.  For example,
51 #  x^2-y^2=0 has solutions along the line y=x and y=-x, so it is
52 #  possible that x-y=0 or x+y=0 will be marked as correct if the
53 #  random points are unluckily chosen.  One way to reduce this problem
54 #  is to increase the number of solutions that are required (by
55 #  setting the ImplicitPoints flag in the Context).  Another is to
56 #  specify the solutions yourself, so that you are sure there are
57 #  points on each component.
58 #
59 #  These problems should be rare, and the values for the various
60 #  parameters have been set in an attempt to minimize the possibility
61 #  of these errors, but they can occur, and you should be aware of
62 #  them, and their possible solutions.
63 #
64 #
65 #  Usage examples:
66 #
67 #     Context("ImplicitEquation");
68 #     \$f = ImplicitEquation("x^2 = cos(y)");
69 #     \$f = ImplicitEquation("x^2 - 2y^2 = 5",limits=>[[-3,3],[-2,2]]);
70 #     \$f = ImplicitEquation("x=1/y",tolerance=>.0001);
71 #
72 #  Then use
73 #
74 #     ANS(\$f->cmp);
75 #
76 #  to get the answer checker for \$f.
77 #
78 #  There are a number of Context flags that control the answer checker.
79 #  These include:
80 #
81 #    ImplicitPoints => 7                   (the number of solutions to test)
82 #    ImplicitTolerance => 1E-6             (relative tolerance value for when
83 #                                           the a tested function is zero)
84 #    ImplicitAbsoluteMinTolerance => 1E-3  (the minimum tolerance allowed)
85 #    ImplicitAbsoluteMaxTolerance => 1E-3  (the maximum tolerance allowed)
86 #    ImplicitPointTolerance => 1E-9        (relative tolerance for how close
87 #                                           the solution point must be to an
88 #                                           actual solution)
89 #    BisectionTolerance => .01             (extra factor used for the tolerance
90 #                                           when finding the solutions)
91 #    BisectionCutoff => 40                 (maximum number of bisections to
92 #                                           perform when looking for a solution)
93 #
94 #  You may set any of these using Context()->flags->set(...).
95 #
96 #  In addition to the Context flags, you can set some values within
97 #  the ImplicitEquation itself:
98 #
99 #    tolerance                (the absolute tolerance for zeros of the function)
100 #    bisect_tolerance         (the tolerance used when searching for zeros)
101 #    point_tolerance          (the absolute tolerance for how close to an
102 #                              actual solution the located solution must be)
103 #    limits                   (the domain to use for the function; see the
104 #                              documentation for the Formula object)
105 #    solutions                (a reference to an array of references to arrays
106 #                              that contain the coordinates of the points
107 #                              that are the solutions of the equation)
108 #
109 #  These can be set in the in the ImplicitEquation() call that creates
110 #  the object, as in the examples below:
111 #
112 #  For example:
113 #
114 #    \$f = ImplicitEquation("x^2-y^2=0",
115 #            solutions => [[0,0],[1,1],[-1,1],[-1,-1],[1,-1]],
116 #            tolerance => .001
117 #         );
118 #
119 #
120 #    \$f = ImplicitEquation("xy=5",limits=>[-3,3]);
121 #
122 #  The limits value can be set globally within the Context, if you wish,
123 #  and the others can be controlled by the Context flags discussed
124 #  above.
125 #
126 ######################################################################
127
128 =cut
129
130 #
131 #  Set up the context for ImplicitEquations and activate it
132 #
133 \$context{ImplicitEquation} = Context("Numeric")->copy;
134 \$context{ImplicitEquation}->variables->are(x=>'Real',y=>'Real');
135 \$context{ImplicitEquation}{precedence}{ImplicitEquation} = Context()->{precedence}{special};
136 Context("ImplicitEquation");
137 Parser::BOP::equality::Allow;
138 \$context{ImplicitEquation}->flags->set(
139   ImplicitPoints => 10,
140   ImplicitTolerance => 1E-6,
141   ImplicitAbsoluteMinTolerance => 1E-3,
142   ImplicitAbsoluteMaxTolerance => 1,
143   ImplicitPointTolerance => 1E-9,
144   BisectionTolerance => .01,
145   BisectionCutoff => 40,
146 );
147
148 #
149 #  Syntactic sugar for creating ImplicitEquations
150 #
151 sub ImplicitEquation {ImplicitEquation->new(@_)}
152
153 #
154 #  Create the ImplicitEquation package
155 #
156 package ImplicitEquation;
157 our @ISA = qw(Value::Formula);
158
159 sub new {
160   my \$self = shift; my \$class = ref(\$self) || \$self;
161   my \$f = shift; return \$f if ref(\$f) eq \$class;
162   \$f = main::Formula(\$f);
163   Value::Error("Your formula doesn't look like an implicit equation")
164     unless \$f->type eq 'Equality';
165   my \$F = (Value::Formula->new(\$f->{tree}{lop}) -
166      Value::Formula->new(\$f->{tree}{rop}))->reduce;
167   \$F = main::Formula(\$F) unless Value::isFormula(\$F);
168   Value::Error("Your equation must be real-valued") unless \$F->isRealNumber;
169   Value::Error("Your equation should not be constant") if \$F->isConstant;
171     if (\$F->usesOneOf(\$F->{context}->variables->parameters));
172   \$F = bless \$F, \$class;
173   my %options = (@_);  # user can supply limits, tolerance, etc.
174   foreach my \$id (keys %options) {\$F->{\$id} = \$options{\$id}}
175   \$F->{F} = \$f; \$F->{isValue} = \$F->{isFormula} = 1;
176   \$F->createPoints unless \$F->{solutions};
177   return \$F;
178 }
179
181
182 #
183 #  Override the comparison method.
184 #
185 #  Turn the right-hand equation into an ImplicitEquation.  This creates
186 #  the test points (i.e., finds the solution points).  Then check
187 #  the professor's function on the student's test points and the
188 #  student's function on the professor's test points.
189 #
190
191 =cut
192
193 sub compare {
194   my (\$l,\$r,\$flag) = @_; my \$tolerance;
195   if (\$l->promotePrecedence(\$r)) {return \$r->compare(\$l,!\$flag)}
196   my @params; @params = (limits=>\$l->{limits}) if \$l->{limits};
197   \$r = ImplicitEquation->new(\$r,@params);
198   Value::Error("Functions from different contexts can't be compared")
199     unless \$l->{context} == \$r->{context};
200
201   #
202   #  They are not equal if couldn't get solutions for one of them
203   #
204   return 1 unless \$l->{solutions} && \$r->{solutions};
205
206   #
207   #  Test the right-hand function on the solutions of the left-hand one
208   #  and vice-versa
209   #
210   my \$rzeros = \$r->createPointValues(\$l->{solutions});
211   my \$lzeros = \$l->createPointValues(\$r->{solutions});
212   return 1 unless \$lzeros && \$rzeros;
213
214   #
215   #  Check that the values are, in fact, zeros
216   #
217   \$tolerance = \$r->getFlag('tolerance',1E-3);
218   foreach my \$v (@{\$rzeros}) {return 1 unless abs(\$v) < \$tolerance}
219   \$tolerance = \$l->getFlag('tolerance',1E-3);
220   foreach my \$v (@{\$lzeros}) {return 1 unless abs(\$v) < \$tolerance}
221
222   return 0; # equal
223 }
224
225 #
226 #  Use the original equation for these (but not for perl(), since we
227 #  need that to use perlFunction).
228 #
229 sub string {shift->{F}->string}
230 sub TeX    {shift->{F}->TeX}
231
232 sub cmp_class {'an Implicit Equation'}
233 sub showClass {shift->cmp_class}
234
235 #
236 #  Locate points that satisfy the equation
237 #
238 sub createPoints {
239   my \$self = shift;
240   my \$num_points = int(\$self->getFlag('ImplicitPoints',10));
241   \$num_points = 1 if \$num_points < 1;
242
243   #
244   #  Get some positive and negative test points (try up to 5 times)
245   #
246   my (\$OK,\$p,\$n,\$z,\$f,@zero); my \$k = 5;
247   while (!\$OK && --\$k) {(\$OK,\$p,\$n,\$z,\$f) = \$self->getPositiveNegativeZero(\$num_points)}
248   Value::Error("Can't find any solutions to your equation") unless \$OK;
249   my (\$P,@intervals) = \$self->getIntervals(\$p,\$n);
250
251   #
252   #  Get relative tolerance values and make them absolute
253   #
254   my \$minTolerance = \$self->getFlag('ImplicitAbsoluteMinTolerance');
255   my \$maxTolerance = \$self->getFlag('ImplicitAbsoluteMaxTolerance');
256   my \$tolerance = \$f * \$self->getFlag('ImplicitTolerance');
257   \$tolerance = \$minTolerance if \$tolerance < \$minTolerance;
258   \$tolerance = \$maxTolerance if \$tolerance > \$maxTolerance;
259   \$self->{tolerance} = \$tolerance unless \$self->{tolerance};
260   \$self->{bisect_tolerance} = \$self->{tolerance} * \$self->getFlag('BisectionTolerance')
261     unless \$self->{bisect_tolerance};
262   \$self->{point_tolerance} = \$P * \$self->getFlag('ImplicitPointTolerance')
263     unless \$self->{point_tolerance};
264
265   #
266   #  Locate solutions to be used for comparison test
267   #
268   @zero = @{\$z}; @zero = \$zero[0..\$num_points-1] if (scalar(@zero) > \$num_points);
269   for (\$i = 0; scalar(@zero) < \$num_points && \$i < scalar(@intervals); \$i++) {
270     my \$Q = \$self->Bisect(\$intervals[\$i],\$intervals[\$i]);
271     push(@zero,\$Q) if \$Q;
272   }
273   Value::Error("Can't find enough solutions for an effective test")
274     unless scalar(@zero) == \$num_points;
275
276   #
277   #  Save the solutions to the equation
278   #
279   \$self->{solutions} = [@zero];
280 }
281
282 #
283 #  Get random points and sort them by sign of the function.
284 #  Also, get the average function value, and indicate if
285 #  we actually did find both positive and negative values.
286 #
287 sub getPositiveNegativeZero {
288   my \$self = shift; my \$n = shift;
289   my (\$p,\$v) = \$self->SUPER::createRandomPoints(3*\$n);
290   my (@pos,@neg,@zero);
291   my \$f = 0; my \$k = 0;
292   foreach my \$i (0..scalar(@{\$v})-1) {
293     if (\$v->[\$i] == 0) {push(@zero,\$p->[\$i])} else {
294       \$f += abs(\$v->[\$i]); \$k++;
295       if (\$v->[\$i] > 0) {push(@pos,\$p->[\$i])} else {push(@neg,\$p->[\$i])}
296     }
297   }
298   \$f /= \$k if \$k;
299   return (scalar(@pos) && scalar(@neg),[@pos],[@neg],[@zero],\$f);
300 }
301
302 #
303 #  Make a list of positive and negative points sorted by
304 #  the distance between them.  Also return the average distance
305 #  between points.
306 #
307 sub getIntervals {
308   my \$self = shift; my \$pos = shift; my \$neg = shift;
309   my @intervals = (); my \$D = 0;
310   foreach my \$p (@{\$pos}) {
311     foreach my \$n (@{\$neg}) {
312       my \$d = abs(Value::Point->make(@{\$p}) - Value::Point->make(@{\$n}));
313       push(@intervals,[\$p,\$n,\$d]); \$D += \$d;
314     }
315   }
316   @intervals = main::PGsort(sub {\$_-> < \$_->},@intervals);
317   return(\$D/scalar(@intervals),@intervals);
318 }
319
320 #
321 #  Use bisection algorithm to find a point where the function is zero
322 #  (i.e., where the original equation is an equality)
323 #  If we can't find a point (e.g., we are at a discontinuity),
324 #  return an undefined value.
325 #
326 sub Bisect {
327   my \$self = shift;
328   my \$tolerance  = \$self->getFlag('bisect_tolerance',1E-5);
329   my \$ptolerance = \$self->getFlag('point_tolerance',1E-9);
330   my \$m = \$self->getFlag('BisectionCutoff',30); my (\$P,\$f);
331   my \$P0 = Value::Point->make(@{\$_}); my \$P1 = Value::Point->make(@{\$_});
332   my (\$f0,\$f1) = @{\$self->createPointValues([\$P0->data,\$P1->data],1)};
333   for (my \$i = 0; \$i < \$m; \$i++) {
334     \$P = (\$P0+\$P1)/2; \$f = \$self->createPointValues([\$P->data]);
335     return unless ref(\$f);
336     \$f = \$f->;
337     return [\$P->value] if abs(\$f) < \$tolerance && abs(\$P1-\$P0) < \$ptolerance;
338     if (\$f > 0) {\$P0 = \$P; \$f0 = \$f} else {\$P1 = \$P; \$f1 = \$f}
339   }
340   return;
341 }
342
343 1;
```