001: /*************************************************************************
002: * *
003: * 1) This source code file, in unmodified form, and compiled classes *
004: * derived from it can be used and distributed without restriction, *
005: * including for commercial use. (Attribution is not required *
006: * but is appreciated.) *
007: * *
008: * 2) Modified versions of this file can be made and distributed *
009: * provided: the modified versions are put into a Java package *
010: * different from the original package, edu.hws; modified *
011: * versions are distributed under the same terms as the original; *
012: * and the modifications are documented in comments. (Modification *
013: * here does not include simply making subclasses that belong to *
014: * a package other than edu.hws, which can be done without any *
015: * restriction.) *
016: * *
017: * David J. Eck *
018: * Department of Mathematics and Computer Science *
019: * Hobart and William Smith Colleges *
020: * Geneva, New York 14456, USA *
021: * Email: eck@hws.edu WWW: http://math.hws.edu/eck/ *
022: * *
023: *************************************************************************/package edu.hws.jcm.functions;
024:
025: import edu.hws.jcm.data.*;
026:
027: /**
028: * A TableFunction is a function that is specified by a table of (x,y)-points.
029: * Values are interpolated between the specified x-values. This can be done
030: * in several differnt ways; the method that is used is controlled by the
031: * "Style" property. Since a TableFunction extends FunctionParserExtension,
032: * a TableFunction can be added to a Parser (provided it has a name), and
033: * it can then be used in expressions parsed by that parser. Note that this
034: * class is meant to be used for functions that are defined by a fairly
035: * small number of points, since each function evaluation involves a linear
036: * search through the list of x-values of the defining points.
037: */
038: public class TableFunction extends FunctionParserExtension {
039:
040: /**
041: * If the style of the function is set to SMOOTH,
042: * then cubic interpolation is used to find the value
043: * of the functions for x-values between the points that define the function.
044: */
045: public static final int SMOOTH = 0;
046:
047: /**
048: * If the style of the function is set to PIECEWISE_LINEAR,
049: * then linear interpolation is used to find the value
050: * of the functions for x-values between the points that define the function.
051: */
052: public static final int PIECEWISE_LINEAR = 1;
053:
054: /**
055: * If the style of the function is set to STEP, then the function is
056: * piecewise constant, and the value of the function at x is taken
057: * from the nearest point in the list of points that define the function.
058: */
059: public static final int STEP = 2;
060:
061: /**
062: * If the style of the function is set to STEP_LEFT, then the function is
063: * piecewise constant, and the value of the function at x is taken
064: * from the nearest point to the left in the list of points that define the function.
065: */
066:
067: public static final int STEP_LEFT = 3;
068:
069: /**
070: * If the style of the function is set to STEP_RIGHT, then the function is
071: * piecewise constant, and the value of the function at x is taken
072: * from the nearest point to the right in the list of points that define the function.
073: */
074: public static final int STEP_RIGHT = 4;
075:
076: private int style; // Type of function, given by one of the above constants.
077: private double[] xCoords = new double[10]; // x-coordinates, in increasing order
078: private double[] yCoords = new double[10]; // the corresponding y-coordinates
079: private CubicSegment[] segments; // interpolation data for SMOOTH style
080: private int pointCt; // number of points (possibly less than xCoords.length)
081:
082: /**
083: * Create a TableFunction with SMOOTH style and no points.
084: */
085: public TableFunction() {
086: this (SMOOTH);
087: }
088:
089: /**
090: * Create a TableFunction with specified style and no points.
091: *
092: * @param style The style for the function: SMOOTH, PIECEWISE_LINEAR, STEP,
093: * STEP_LEFT, or STEP_RIGHT.
094: */
095: public TableFunction(int sytle) {
096: this .style = style;
097: if (style == SMOOTH)
098: segments = new CubicSegment[9];
099: }
100:
101: /**
102: * Copy data from another TableFunction, except that the name of the funcion is
103: * not duplicated. The new TableFunction is nameless.
104: */
105:
106: public void copyDataFrom(TableFunction source) {
107: xCoords = (double[]) source.xCoords.clone();
108: yCoords = (double[]) source.yCoords.clone();
109: style = -1; // Force setStyle to compute data for SMOOTH style.
110: setStyle(source.style);
111: }
112:
113: /**
114: * Set the style of this TableFunction, to specify how values are interpolated
115: * between points on the curve.
116: *
117: * @param style One of the style constants SMOOTH, PIECEWISE_LINEAR, STEP,
118: * STEP_LEFT, STEP_RIGHT. Other values are ignored.
119: */
120: public void setStyle(int style) {
121: if (style == this .style || style < 0 || style > 4)
122: return;
123: this .style = style;
124: if (style == SMOOTH) {
125: segments = new CubicSegment[xCoords.length - 1];
126: for (int i = 0; i < pointCt - 1; i++)
127: segments[i] = new CubicSegment(xCoords[i],
128: xCoords[i + 1], yCoords[i], yCoords[i + 1], 0,
129: 0);
130: for (int i = 0; i < pointCt - 1; i++)
131: interpolateDerivatives(i);
132: } else {
133: segments = null;
134: }
135: }
136:
137: /**
138: * Get the style of this TableFunction, which specifies how values are
139: * interpolated between points on the curve.
140: *
141: * @return The style of this TableFunction. This is one of the constants
142: * SMOOTH, PIECEWISE_LINEAR, STEP, STEP_LEFT, or STEP_RIGHT.
143: */
144: public int getStyle() {
145: return style;
146: }
147:
148: /**
149: * Add points to the table. The x-coordinates of the points are taken from
150: * the xCoords array. The y-coordinate for the i-th point is yCoords[i], if
151: * an i-th position exists in this array. Otherwise, the y-coordinate is
152: * is zero. (Note that if xCoords[i] duplicates an x-value already in
153: * the table, then no new point is added but the corresponging y-value is changed.)
154: *
155: * @param xCoords A list of x-coordinates to be added to the table. If this is
156: * null, then nothing is done.
157: * @param yCoords The value of yCoords[i], if it exists, is the y-coordinate
158: * corresponding to xCoords[i]. Otherwise, the y-coordinate is undefined.
159: * This can be null, in which case all y-coordinates are zero.
160: */
161: public void addPoints(double[] xCoords, double[] yCoords) {
162: if (xCoords == null)
163: return;
164: int ct = xCoords.length;
165: if (yCoords == null)
166: ct = 0;
167: else if (yCoords.length < ct)
168: ct = yCoords.length;
169: for (int i = 0; i < ct; i++)
170: addPoint(xCoords[i], yCoords[i]);
171: for (int i = ct; i < xCoords.length; i++)
172: addPoint(xCoords[i], 0);
173: }
174:
175: /**
176: * Add points to the table. The number of points added is intervals + 1.
177: * The x-coordinates are evenly spaced between xmin and xmax. The y-coordinates
178: * are zero.
179: *
180: * @param intervals The number of intervals. The number of points added is intervals + 1.
181: * The value should be at least 1. If not, nothing is done.
182: * @param xmin The minimim x-coordinate for added points.
183: * @param xmax The maximum x-coodinate for added points. Should be greater than
184: * xmin, for efficiency, but no error occurs if it is not.
185: */
186: public void addIntervals(int intervals, double xmin, double xmax) {
187: if (intervals < 1)
188: return;
189: double dx = (xmax - xmin) / intervals;
190: for (int i = 0; i < intervals; i++)
191: addPoint(xmin + i * dx, 0);
192: addPoint(xmax, 0);
193: }
194:
195: /**
196: * Add a point with the specified x and y coordinates. If a point with the
197: * given x coordinate already exists in the table, then no new point is added,
198: * but the associated y-value is changed.
199: * (If x is Double.NaN, then no change is made and the return value is -1.)
200: *
201: * @param x The x-coordinate of the point to be added or modified.
202: * @param y The y-coordinate of the point.
203: * @return the position of the point in the list of points, where the first point is at position zero.
204: */
205: public int addPoint(double x, double y) {
206: if (Double.isNaN(x))
207: return -1;
208: int pos = 0;
209: while (pos < pointCt && xCoords[pos] < x)
210: pos++;
211: if (pos < pointCt && xCoords[pos] == x) {
212: yCoords[pos] = y;
213: } else {
214: if (pointCt == xCoords.length) {
215: double[] temp = new double[2 * xCoords.length];
216: System.arraycopy(xCoords, 0, temp, 0, xCoords.length);
217: xCoords = temp;
218: temp = new double[2 * yCoords.length];
219: System.arraycopy(yCoords, 0, temp, 0, yCoords.length);
220: yCoords = temp;
221: if (style == SMOOTH) {
222: CubicSegment[] temps = new CubicSegment[xCoords.length - 1];
223: System
224: .arraycopy(segments, 0, temps, 0,
225: pointCt - 1);
226: segments = temps;
227: }
228: }
229: for (int i = pointCt; i > pos; i--) {
230: xCoords[i] = xCoords[i - 1];
231: yCoords[i] = yCoords[i - 1];
232: }
233: xCoords[pos] = x;
234: yCoords[pos] = y;
235: if (style == SMOOTH && pointCt > 0) {
236: if (pos == pointCt)
237: segments[pointCt - 1] = new CubicSegment();
238: else {
239: for (int i = pointCt - 1; i > pos; i--)
240: segments[i] = segments[i - 1];
241: segments[pos] = new CubicSegment();
242: }
243: }
244: pointCt++;
245: }
246: if (style == SMOOTH && pointCt > 0) { // make sure segment data is OK for segments that depend on (x,y)
247: if (pos > 0)
248: segments[pos - 1].setData(xCoords[pos - 1],
249: xCoords[pos], yCoords[pos - 1], yCoords[pos],
250: 0, 0);
251: if (pos < pointCt - 1)
252: segments[pos].setData(xCoords[pos], xCoords[pos + 1],
253: yCoords[pos], yCoords[pos + 1], 0, 0);
254: for (int i = pos - 2; i <= pos + 1; i++)
255: interpolateDerivatives(i);
256: }
257: return pos;
258: }
259:
260: private void interpolateDerivatives(int pos) {
261: // Compute correct derivatives for segments[pos] from data in coordinate arrays.
262: if (pos < 0 || pos > pointCt - 2)
263: return;
264: // pointCt must be >= 2
265: if (pointCt == 2) // pos must be 0
266: segments[0].setDerivativesFromNeighbors(Double.NaN, 0,
267: Double.NaN, 0);
268: else if (pos == 0)
269: segments[0].setDerivativesFromNeighbors(Double.NaN, 0,
270: xCoords[2], yCoords[2]);
271: else if (pos == pointCt - 2)
272: segments[pointCt - 2].setDerivativesFromNeighbors(
273: xCoords[pointCt - 3], yCoords[pointCt - 3],
274: Double.NaN, 0);
275: else
276: segments[pos].setDerivativesFromNeighbors(xCoords[pos - 1],
277: yCoords[pos - 1], xCoords[pos + 2],
278: yCoords[pos + 2]);
279: }
280:
281: /**
282: * Gets the number of points in the table.
283: */
284: public int getPointCount() {
285: return pointCt;
286: }
287:
288: /**
289: * Get the x-coordinate in the i-th point, where the first point
290: * is number zero. Throws an IllegalArgumentException if i is
291: * less than zero or greater than or equal to the number of points.
292: */
293: public double getX(int i) {
294: if (i >= 0 && i < pointCt)
295: return xCoords[i];
296: else
297: throw new IllegalArgumentException(
298: "Point index out of range: " + i);
299: }
300:
301: /**
302: * Get the y-coordinate in the i-th point, where the first point
303: * is number zero. Throws an IllegalArgumentException if i is
304: * less than zero or greater than or equal to the number of points.
305: */
306: public double getY(int i) {
307: if (i >= 0 && i < pointCt)
308: return yCoords[i];
309: else
310: throw new IllegalArgumentException(
311: "Point index out of range: " + i);
312: }
313:
314: /**
315: * Set the y-coordinate in the i-th point to y, where the first point
316: * is number zero. Throws an IllegalArgumentException if i is
317: * less than zero or greater than or equal to the number of points.
318: */
319: public void setY(int i, double y) {
320: if (i >= 0 && i < pointCt)
321: yCoords[i] = y;
322: else
323: throw new IllegalArgumentException(
324: "Point index out of range: " + i);
325: if (style == SMOOTH) {
326: if (i > 0)
327: segments[i - 1].setData(xCoords[i - 1], xCoords[i],
328: yCoords[i - 1], yCoords[i], 0, 0);
329: if (i < pointCt - 1)
330: segments[i].setData(xCoords[i], xCoords[i + 1],
331: yCoords[i], yCoords[i + 1], 0, 0);
332: for (int j = i - 2; j <= i + 1; j++)
333: interpolateDerivatives(j);
334: }
335: }
336:
337: /**
338: * If there is a point in the list with x-coordinate x, then this function returns
339: * the index of that point in the list (where the index of the first point is zero).
340: * If there is no such point, then -1 is returned.
341: */
342: public int findPoint(double x) {
343: int i = 0;
344: while (i < pointCt) {
345: if (x == xCoords[i])
346: return i;
347: else if (x > xCoords[i])
348: i++;
349: else
350: break;
351: }
352: return -1;
353: }
354:
355: /**
356: * Removes the i-th point from the list of points. Throws an IllegalArgumentException if i is
357: * less than zero or greater than or equal to the number of points.
358: */
359: public void removePointAt(int i) {
360: if (i < 0 || i >= pointCt)
361: throw new IllegalArgumentException(
362: "Point index out of range: " + i);
363: pointCt--;
364: for (int j = i; j < pointCt; j++) {
365: xCoords[j] = xCoords[j + 1];
366: yCoords[j] = yCoords[j + 1];
367: }
368: if (style == SMOOTH) {
369: style = -1; // force recompute of data
370: setStyle(SMOOTH);
371: }
372: }
373:
374: /**
375: * Remove all points. The resulting function is undefined everywhere.
376: */
377: public void removeAllPoints() {
378: pointCt = 0;
379: xCoords = new double[10];
380: yCoords = new double[10];
381: }
382:
383: /**
384: * Get the value of the function at x, using interpolation if x lies between
385: * two x-coordinates in the list of points that define the function. If x is
386: * outside the range of x-coords in the table, the value of the function is Double.NaN.
387: */
388: public double getVal(double x) {
389: return computeValue(x, null, 0);
390: }
391:
392: private double computeValue(double x, Cases cases,
393: int derivativeLevel) {
394: // Find the value of the function or one of its derivatives at x.
395: // If cases is not null, then a value is added to cases to help
396: // with continuity computations.
397: if (Double.isNaN(x))
398: return Double.NaN;
399: if (pointCt == 0 || x < xCoords[0] || x > xCoords[pointCt - 1])
400: return Double.NaN;
401: if (pointCt == 1) {
402: if (derivativeLevel > 0)
403: return Double.NaN;
404: else {
405: if (cases != null)
406: cases.addCase(0);
407: return yCoords[0];
408: }
409: }
410: int casenum;
411: double ans;
412: int seg = 0;
413: switch (style) {
414: case STEP: {
415: while (seg < pointCt - 1
416: && x > (xCoords[seg] + xCoords[seg + 1]) / 2)
417: seg++;
418: casenum = seg;
419: if (derivativeLevel == 0)
420: ans = yCoords[seg];
421: else if (x < (xCoords[seg] + xCoords[seg + 1]) / 2
422: || seg == pointCt - 1
423: || yCoords[seg] == yCoords[seg + 1])
424: ans = 0;
425: else
426: ans = Double.NaN;
427: break;
428: }
429: case STEP_RIGHT: {
430: while (seg < pointCt - 1 && x > xCoords[seg])
431: seg++;
432: casenum = seg;
433: if (derivativeLevel == 0)
434: ans = yCoords[seg];
435: else if (x < xCoords[seg] || seg >= pointCt - 1
436: || yCoords[seg] == yCoords[seg + 1])
437: ans = 0;
438: else
439: ans = Double.NaN;
440: break;
441: }
442: case STEP_LEFT: {
443: while (seg < pointCt - 1 && x >= xCoords[seg + 1])
444: seg++;
445: casenum = seg;
446: if (derivativeLevel == 0)
447: ans = yCoords[seg];
448: else if (x > xCoords[seg] || seg == 0
449: || yCoords[seg] == yCoords[seg - 1])
450: ans = 0;
451: else
452: ans = Double.NaN;
453: break;
454: }
455: case PIECEWISE_LINEAR: {
456: while (seg < pointCt - 1 && x > xCoords[seg])
457: seg++;
458: casenum = seg;
459: if (x == xCoords[seg]) {
460: if (derivativeLevel == 0)
461: ans = yCoords[seg];
462: else if (seg == 0) {
463: if (derivativeLevel == 1)
464: ans = (yCoords[1] - yCoords[0])
465: / (xCoords[1] - xCoords[0]);
466: else
467: ans = 0;
468: } else if (seg == pointCt - 1) {
469: if (derivativeLevel == 1)
470: ans = (yCoords[pointCt - 1] - yCoords[pointCt - 2])
471: / (xCoords[pointCt - 1] - xCoords[pointCt - 2]);
472: else
473: ans = 0;
474: } else {
475: double leftslope = (yCoords[seg] - yCoords[seg - 1])
476: / (xCoords[seg] - xCoords[seg - 1]);
477: double rightslope = (yCoords[seg] - yCoords[seg + 1])
478: / (xCoords[seg] - xCoords[seg + 1]);
479: if (Math.abs(leftslope - rightslope) < 1e-12) {
480: if (derivativeLevel == 1)
481: ans = leftslope;
482: else
483: ans = 0;
484: } else
485: ans = Double.NaN;
486: }
487: } else { // x < xCoords[seg] && x > xCoords[seg-1]
488: if (derivativeLevel == 0) {
489: ans = yCoords[seg - 1]
490: + ((yCoords[seg] - yCoords[seg - 1]) / (xCoords[seg] - xCoords[seg - 1]))
491: * (x - xCoords[seg - 1]);
492: } else if (derivativeLevel == 1) {
493: ans = (yCoords[seg] - yCoords[seg - 1])
494: / (xCoords[seg] - xCoords[seg - 1]);
495: } else
496: ans = 0;
497: }
498: break;
499: }
500: default: { // SMOOTH
501: while (seg < pointCt - 2 && x > xCoords[seg + 1])
502: seg++;
503: casenum = seg;
504: if (x == xCoords[seg + 1] && seg < pointCt - 2 && seg > 0) {
505: if (derivativeLevel == 0)
506: ans = yCoords[seg + 1];
507: else if (derivativeLevel == 1)
508: ans = segments[seg].derivativeValue(x, 1);
509: else {
510: double leftslope = segments[seg - 1]
511: .derivativeValue(x, 2);
512: double rightslope = segments[seg].derivativeValue(
513: x, 2);
514: if (Math.abs(leftslope - rightslope) < 1e-12)
515: ans = segments[seg].derivativeValue(x,
516: derivativeLevel);
517: else
518: ans = Double.NaN;
519: }
520: } else { // x < xCoords[seg] && x > xCoords[seg-1]
521: ans = segments[seg].derivativeValue(x, derivativeLevel);
522: }
523: }
524: }
525: if (cases != null)
526: cases.addCase(casenum);
527: return ans;
528: }
529:
530: //---------------- Methods from the Function class -----------------
531:
532: /**
533: * Get the value of the function at the specified parameter value.
534: *
535: * @params param should be an array of length 1 holding the argument of the function.
536: * However if the length is greater than one, the extra arguments are simply ignored.
537: * @cases if cases is non-null, a case value is stored here, for help in continuity computations.
538: */
539: public double getValueWithCases(double[] params, Cases cases) {
540: return computeValue(params[0], cases, 0);
541: }
542:
543: /**
544: * Get the value of the function at the specified parameter value.
545: *
546: * @partam param should be an array of length 1 holding the argument of the function.
547: * However if the length is greater than one, the extra arguments are simply ignored.
548: */
549: public double getVal(double[] params) {
550: return computeValue(params[0], null, 0);
551: }
552:
553: /**
554: * Compute the derivative of this function. The value of the parameter, wrt, must be 1 or an
555: * IllegalArguemntException will be thrown.
556: */
557: public Function derivative(int wrt) {
558: if (wrt != 1)
559: throw new IllegalArgumentException(
560: "Attempt to take the derivative of a function of one argument with respect to argument number "
561: + wrt);
562: return new Deriv(this );
563: }
564:
565: /**
566: * Returns null.
567: * It really should be the constant function zero, but I don't expect this ever to be
568: * called. Since dependsOn(wrt) returns false, it will never be called within the JCM system.
569: */
570: public Function derivative(Variable wrt) {
571: return null;
572: }
573:
574: /**
575: * Returns false.
576: */
577: public boolean dependsOn(Variable wrt) {
578: return false;
579: }
580:
581: /**
582: * Returns the arity of the function, which is 1.
583: */
584: public int getArity() {
585: return 1;
586: }
587:
588: /**
589: * Override method apply() from interface FunctionParserExtension, to handle cases properly.
590: * Not meant to be called directly.
591: */
592: public void apply(StackOfDouble stack, Cases cases) {
593: double x = stack.pop();
594: stack.push(computeValue(x, cases, 0));
595: }
596:
597: //---------------- For creating derivatives ------------------------
598:
599: private static class Deriv extends FunctionParserExtension {
600: // An object of this class represents a derivative function
601: // of a TableFunction.
602:
603: TableFunction derivativeOf; // The function from which this function is derived.
604:
605: int derivativeLevel; // The order of the derivative.
606:
607: Deriv(Deriv f) {
608: derivativeLevel = f.derivativeLevel + 1;
609: derivativeOf = f.derivativeOf;
610: }
611:
612: Deriv(TableFunction f) {
613: derivativeLevel = 1;
614: derivativeOf = f;
615: }
616:
617: public String getName() { // Name comes from the name of the function.
618: String name = derivativeOf.getName();
619: for (int i = 0; i < derivativeLevel; i++)
620: name += "'";
621: return name;
622: }
623:
624: public void setName(String name) {
625: }
626:
627: public double getValueWithCases(double[] params, Cases cases) {
628: return derivativeOf.computeValue(params[0], cases,
629: derivativeLevel);
630: }
631:
632: public double getVal(double[] params) {
633: return derivativeOf.computeValue(params[0], null,
634: derivativeLevel);
635: }
636:
637: public Function derivative(int wrt) {
638: if (wrt != 1)
639: throw new IllegalArgumentException(
640: "Attempt to take the derivative of a function of one argument with respect to argument number "
641: + wrt);
642: return new Deriv(this );
643: }
644:
645: public Function derivative(Variable wrt) {
646: return null;
647: }
648:
649: public boolean dependsOn(Variable wrt) {
650: return false;
651: }
652:
653: public int getArity() {
654: return 1;
655: }
656:
657: public void apply(StackOfDouble stack, Cases cases) {
658: double x = stack.pop();
659: stack.push(derivativeOf.computeValue(x, cases,
660: derivativeLevel));
661: }
662:
663: } // end class Deriv
664:
665: }
|