Math Programming

class Math
{
 private $bcmath_enabled = false;
 private $gmp_enabled = false;
 // class methods here
}

Analogy to database abstraction libraries What are the arguments for using a number abstraction library? Are they the same arguments as are used to justify using a database abstraction library in developing web applications? Mostly yes so I don’t have to bother repeating them 🙂 One of the main counterarguments to using a database abstraction library is that you add some bloat and indirection to your code (i.e., the abstraction code). This means there is a market for “lite” classes especially when it comes to interposing a php-based numerical abstraction library between or amidst all your math inputs and outputs. Other work in this areaYour’s truly did a bit of work in this area by porting the Weka comparison methods from Java to PHP. If you study some of the massive Weka codebase, you will notice that these comparison functions are included and used in most of the classes you will encounter. All 6 comparison operations in this class use a single global error tolerance setting to define the level of precision to use:

/**
* The small deviation allowed in floating point comparisons.
*/
const SMALL = 1e-6;

One of the main advantages of this approch is that you can more easily figure out what your practical error bounds are if you know the order of precision your machine is using to represent the quantities involved. Eliminating sources of error Using an arbitrary precision is one way to try to address a common source of error in your numerical computing, namely, truncation-related errors. It is important to recognize, however, that you can also achieve practically useful error bounds on your outputs without using arbitrary precision if your code implements numerically stable approaches to performing the computations involved and error tolerances are specified and used in your comparisons. My definition of numerical stability is very practical: a property of numerical algorithms that allows them to be ported to a new language and/or hardware platform (without any special accomodations) and continue to perform practically useful computations on that new platform. The JAMA package, for example, appears to implement numerically stable matrix manipulation and decomposition algorithms because I and a small team of open-source programmers (see docs) were able to port them from Java to PHP without any significant modifications to the algorithms and achieved practically identical results between the PHP and Java versions. Complex numbers Jesus M. Castagnetto has written a nice Complex numbers package. Would we eventually need to extend a number abstraction library to support arbitrary precision operations with complex numbers as well? Suggestions I think Alix has made some good progress towards a number abstraction library for PHP. There are three main components missing in my opinion: 1) simple examples of class usage, 2) a unit testing class that asserts the expected outcome for each class method and checks the computed outcome against the expected outcome, 3) a non-trivial example of a math solution implemented using the library (e.g., I have a nice guassian elimination port that would be appropriate as an example). In other words, goodies in the demos, tests, and examples subfolders. I’d like to thank Alix for sharing some of his thought-provoking work with me.

quadratic method: 3.13333333333
cubic method: 3.13846153846
adaptive method: 3.14159265371
composite method: 3.14159265297

The method producing the output closest to an exact π value of 3.1415926535 is the adaptive method, beating out the composite method by a tiny margin. Other classes that might be found in a Math_Integration package include:

• Math_Integration_Trapezoid.php
• Math_Integration_Gaussian.php (i.e., gaussian quadrature)
• Math_Integration_Romberg.php

All would have basic, adaptive, and composite variants. Exercises 1. Play around with the $eps parameter in the adaptive Simpson’s rule and the $n parameter in the composite Simpson’s rule. Can the composite method be made to out perform the adaptive method in computing π? 2. Apply these methods to a variety of other functions and note any differences. 3. Implement a Math_Integration_Trapezoid.php class that contains a basic, adaptive, and composite implementation of the trapezoid method of integration. 4. Compare the relative performance of Simpson’s rule to the Trapezoid method of integration on functions that are smooth or contain abrupt changes. 5. How do Lagrange polynomials relate to Simpson’s rule? 6. How do the Newton-Cotes formulas relate to Simpson’s rule and other methods of integration?

Given [1, 3] as input, the function x-3 outputs [1, 0.037037037037]

What I like about this approach is that you can just define and instantiate My_Function as an $f object and pass the $f object as an argument to your numerical analysis function. Once inside your numerical analysis function, you can evaluate your function by invoking the valueAt($x) method which makes your numerical analysis code quite readable. You often need to output the function you are evaluating so why not implement toString() and toHTML() methods at the same time you are implementing your valueAt($x) method. Exercise 1. What happens if you simply call the class Function instead of My_Function? 2. What might the output of a toXML() method look like? See the MathML website for ideas.

1. A setCoefficients method for finding the polynomial coefficients to use for the interpolating polynomial.
2. A setInterpolants method for evaluating the polynomial formula at intermediate points using the computed polynomial coefficients.

You can peruse the code comments for more details about how the class works. <?php /** * Math_Interpolation_Lagrange.php * * The Lagrange interpolation object accepts three arguments: * * 1. A data vector of x values * 2. A data vector of y values * 3. A vector of ix values, in the same range as your x values, that * you want to find the cooresponding interpolated values iy for. * * It uses the data vector of x and y values to compute the polynomial * coefficients c and then uses these coefficients to find the interpolated * values iy for the supplied ix values. All supplied and computed results * are stored as instance variables obtainable by accessor methods or * direct reference. * * @author Paul Meagher * @license LGPL * @created Feb 25/2007 * @version 1.0 */ class Math_Interpolation_Lagrange { /** * @var array of x values */ var $x = array(); /** * @var array of y values */ var $y = array(); /** * @var array of polynomial coefficients */ var $c = array(); /** * @var array of x values to be interpolated */ var $ix = array(); /** * @var array of interpolated y values */ var $iy = array(); /** * @var number of elements in data arrays x and y */ var $n = 0; /** * @var number of elements in interpolation arrays ix and iy */ var $m = 0; /** * Constructor sets up the internal storage variables, calls a * method to compute the polynomial cooefficients, then finds * the interpolated iy values given the supplied ix values. * * @param $x array of x values * @param $y array of y values * @param $ix array of x values to be interpolated */ function Math_Interpolation_Lagrange($x, $y, $ix) { $this->x = $x; $this->y = $y; $this->ix = $ix; $this->n = count($x); $this->m = count($ix); $this->setCoefficients(); $this->setInterpolants(); } /** * Computes the polynomial coefficients. */ function setCoefficients() { for($i=0; $i<$this->n; $i++) { $d[$i] = 1; for($j=0; $j<$this->n; $j++) { if($i != $j) $d[$i] = $d[$i] * ($this->x[$i] - $this->x[$j]); $this->c[$i] = $this->y[$i] / $d[$i]; } } } /** * Computes the interpolated iy values given the ix values * supplied in the constructor. */ function setInterpolants() { for($i=0; $i<$this->m; $i++) { $this->iy[$i] = 0; for($j=0; $j<$this->n; $j++) { $d[$j] = 1; for($k=0; $k<$this->n; $k++) { if($j != $k) $d[$j] = $d[$j] * ($this->ix[$i] - $this->x[$k]); } $this->iy[$i] = $this->iy[$i] + $this->c[$j] * $d[$j]; } } } /** * Recomputes the interpolated iy values given the new ix values. * * @param $ix array of new x values to be interpolated */ function changeInterpolants($ix) { $this->ix = $ix; $this->iy = array(); $this->m = count($ix); $this->setInterpolants(); } /** * @return the polynomial coefficients */ function getCoefficients() { return $this->c; } /** * @return array of interpolated values */ function getInterpolants() { return $this->iy; } /** * @return interpolated values as comma separated string */ function toString() { return implode(", ", $this->iy); } } ?> The lagrange_test.php script demonstrates how to use the class to obtain two sets of interpolant values: <?php include "Math/Interpolation/Lagrange.php"; $x = array(-2, 0, -1, 1, 2); $y = array( 4, 2, -1, 1, 8); // create x values in data range using step size of .5 $ix1 = range(-2, 2, .5); // create x values in data range using step size of .2 $ix2 = range(-2, 2, .2); $MIL = new Math_Interpolation_Lagrange($x, $y, $ix1); echo "<p>Interpolated values using step size of .5 are: ".$MIL->toString()."</p>"; $MIL->changeInterpolants($ix2); echo "<p>Interpolated values using step size of .2 are: ".$MIL->toString()."</p>"; ?> Here is the output the script generates: Interpolated values using step size of .5 are: 4, -1.1875, -1, 0.8125, 2, 1.8125, 1, 1.8125,8 Interpolated values using step size of .2 are: 4, 0.9776, -0.7264, -1.4384, -1.4464, -1, -0.3104, 0.4496, 1.1456, 1.6816, 2, 2.0816, 1.9456, 1.6496, 1.2896, 1,0.9536, 1.3616, 2.4736, 4.5776, 8 Exercise 1. Use a PHP graphing package like JpGraph or Image::Graph to plot your x and y data points and the interpolated points ix and iy. Draw a line through your original data points and your interpolated data points. 2. Use a formula to generate a few data points and then use the Math_Interpolation_Lagrange class to find interpolating values. Come up with a method to measure how accurate your Lagrange interpolating polynomial is by comparing interpolated values to formula generated values. Try using different formulas and compare the relative accuracy of the Lagrange method for those formulas. 3. Use the Runge Function to generate some data points and see how good the Lagrange method is at interpolating values. Note: The Runge function is known to give the Lagrange interpolation method some problems. 4. Create a Math_Interpolation_Chebyshev class to minimize the weaknesses of the Lagrange method for interpolating the Runge function. See this reference for further details on Chebyshev polynomials.

• • array array_map ( callback callback, array arr1 [, array ...] )
    

    array_map() returns an array containing all the elements of arr1 after applying the callback function to each one. The number of parameters that the callback function accepts should match the number of arrays passed to the array_map().

• array array_filter ( array input [, callback callback] )
  

  array_filter() iterates over each value in the input array passing them to the callback function. If the callback function returns true, the current value from input is returned into the result array. Array keys are preserved.

• mixed array_reduce ( array input, callback function [, int initial] )
  

  array_reduce() applies iteratively the callback function to the elements of the array input, so as to reduce the array to a single value. If the optional initial is available, it will be used at the beginning of the process, or as a final result in case the array is empty. If the array is empty and initial is not passed, array_reduce() returns NULL.

The size of an integer is platform-dependent, although a maximum value of about two billion is the usual value (that’s 32 bits signed). PHP does not support unsigned integers.

Elliot Back has written a program to find the maximum representable integer which also provides some insight into how the 32 bits used to represent integers are allocated. From the PHP Manual regarding floating point datatypes:

The size of a float is platform-dependent, although a maximum of ~1.8e308 with a precision of roughly 14 decimal digits is a common value (that’s 64 bit IEEE format).

Chapell (1995) provides some insight into how the 64 bits are allocated:

A double precision variable is usually 8 bytes (16 bits) long. In a typical implementation, 53 bits of the number are devoted to the mantissa, and 11 bits are devoted to the exponent. The 53 bits devoted to the mantissa are enough to represent 15 or 16 significant decimal digits. Similarly, the 11 bits of the exponent are enough to represent numbers as large as 10³⁰⁸ and as small as 10^-308. (p. 474)

Learn more about floating point representation from a C perspective. Consult Mike Cowlishaw’s link collection for everything you ever wanted to know about decimal arithmetic. Exercises 1. Develop a PHP program that finds the maximum floating point number. 2. Does an integer have mantissa and exponent parts? 3. Is there a way to know exactly (instead of just “roughly”) how many bits are used to represent the mantissa and exponent parts of a PHP floating point number.

1. PHP offers integer and float data types (a double precision data type) to accomodate numerical quantities.
2. Some textbooks on numerical analysis use plain integer and double precision data types (aka float) to show-case the math used to estimate relative and absolute round off and truncation error for a group of math algorithms (i.e., different math algorithms for estimating the first and second derivative of a function when that function is linear, quadratic, cubic, nonlinear, etc…).

Deferring to Yakowitz and Szidarovsky (1989), a Truncation Error “is introduced by approximating some ideal mathematical operation by computations directed by a program. It presumes that computations are done without roundoff error.” (p. 12) Another quote from William H. Press et al. (2003) clarifies the importance of trucation error:

The discrepancy between the true answer and the answer obtained in a practical calculation is called the truncation error. Truncation error would persist even on a hypothetical , “perfect” computer that had an infinitely accurate representation and no roundouff error. As a general rule there is not much that a programmer can do about roundoff error, other than to choose algortithms that do not magnify it unnecesessarily . Truncation error, on the other hand, is entirley under the programmer’s control. In fact, it is only a slight exaggeration to say that clever minimisation of truncation error is practically the entire content of the field of numerical analysis! (p. 33)

Exercises 1. Collect some other definitions of a truncation error and note any differences in the way in which numerical errors are classified or subclassified. 2. Provide an example of a trunction error.

2.5
-8.5
10

Which are the factors in the interpolating polynomial equation: y = 2.5x² – 8.5x + 10 Exercises 1. Add a method to this class that will find the y value given a particular x value (within the range of existing x values). In other words, it will use the computed polynomial equation to interpolate the value of y given a value of x in the domain of the function. 2. Add a check to the evaluation routine above that warns the user when they try to use to extrapolate (x value outside of the domain) rather than interpolate (x value within the domain) a y value.

10
100
1000
10000
100000
1000000
10000000
100000000
1000000000
10000000000
100000000000
1E+012

If you want to render scientific notation number larger than 1E+012 in decimal format, then you should use the sprintf or printf functions to do so as the following program illustrates: <?php for ($i=1; $i <= 20; $i++) { $val = "1.0e$i"; printf("%$i.0f<br />", $val); } ?>

10
100
1000
10000
100000
1000000
10000000
100000000
1000000000
10000000000
100000000000
1000000000000
10000000000000
100000000000000
1000000000000000
10000000000000000
100000000000000000
1000000000000000000
10000000000000000000
100000000000000000000

Exercises 1. Develop a php script to cast numbers in scientific notation to decimal notation when the numbers are extremely small. In other words, your scientfic notation values should have increasing negative exponents. Use float and printf casting to explore the limits of these two casting methods. 2. Develop a php script to determine when the printf or sprintf functions encounter underflow and overflow conditions. In other words, what are the smallest and largest numbers that they can be used to format. 3. What, if any, is the relationship between printf or sprintf underflow and machine epsilon?

If you have a 10000 objects in a set, and you want to check if another element is in there, what do you do? You don’t actually go through the whole list, do you?

I suggested that PHP’s in_array() function might do the trick however Daniel persisted in asking me with whether the lookup was O(n) or O(1) like Python (a language Daniel was more intimately familiar with). In the PHP manual under Arrays there is a comment by S Boisvert that clarified the answer. Essentially, in_array() is (O)n but if you can turn your array lookup problem into a test of whether an array key exists (i.e., isset($array[$key])), you can get O(1) speed because it is a hashed lookup then. To further test this assertion I developed a program that tests array lookups using the in_array() approach versus the isset() approach and was able to conclusively demonstrate that the isset() approach is the O(1) way to do array lookups. Thanks to Daniel for helping to educate me on some of the finer points of PHP arrays and a practical use for the big O notation.

chmod -Rf 777 downloads

All of this points towards the need for a PHP/Math installation system that will create the relevant “downloads” folders, change their permissions, and replace into the “Math/PDL”, “Math/JAMA”, and “Math/REGRESS” system folders the currently downloaded PDL, JAMA, and REGRESS packages. I will be putting together a bash shell script to implement these installation steps. Using a standard set of subfolder namesSome JAMA scripts are currently broken because I standardized on a set of package subfolder names and changed JAMA subfolder names to use these standard names (e.g., “doc” changed to “docs”, “test” changed to “tests”, “util” changed to “utils”). My standard set of subfolder names currently includes:

/examples
/tests
/lib
/utils
/functions
/docs
/docs/includes
/docs/sections
/docs/images

Developers don’t have to use all these subfolder names for any given project, only the subfolders which are required. Other subfolder naming options might be added in the future. One argument for doing things in this way is that the download system can be developed with less ad-hocery (i.e., no need to individually customize which subfolders are to be included for each package; interate through a standard set of subfolder names instead). A standardized set of subfolder names also allows developers and users to develop a consistent set of expectations about how current and future PHP/Math packages might be browsed online, code downloaded, and code installed.

include("BunnyRegex.php");
$pattern = new BunnyRegex();
$pattern->bol() // ^
 ->digit() // d
 ->exactly(4) // {4}
 ->string('/') // /
 ->digit() // d
 ->exactly(2) // {2}
 ->string('/') // /
 ->digit() // d
 ->exactly(2) // {2}
 ->eol(); // $
$pattern->match('2006/02/28'); // returns true
$pattern->numberOfMatches('2006/02/28'); // returns 1

Implementing fluent interfaces in PHP In Making Fluent Interfaces Readable in PHP, Jonnay deconstructs fluent interfaces in these terms:

Fluent Interfaces are a method of programming where setters, and sometimes behavioural methods that normally return void/null actually return an object instance, usually (but not always) the object being operated on. This is so that you can get around the effect where you make a temporary variable to alter an object, but don’t end up using it again.

Making JAMA port more fluent The Java-based JAMA linear algebra library is already quite fluent as illustrated by this line of computation: R = L.times(U).minus(M.getMatrix(p,0,n-1)); Perhaps fluency holds the secret to helping me fully port the JAMA library so that it can be similiarly fluent. Currently, this is how I compute the value of $R:

$S = $L->times($U);
$R = $S->minus($M->getMatrix($p,0,$n-1));

Ideally it would be:

$R = $L->times($U)->minus($M->getMatrix($p,0,$n-1));

Just like the JAVA version. Aspect Oriented Programming AOP Jonnay has written a followup article Aspect Oriented Programming in PHP as a contrast to other languages which embeds this line of development and thinking in the AOP tradition. Now that you know about fluency of the AOP sort, you can use this knowledge to create some very powerful API’s for your php-based mathematical classes.

:
include "Math/PDL/BinomialDistribution.php";
include "Math/JAMA/QRDecomposition.php";
include "Math/CHI/CHI2D.php";
include "Math/IT/JointEntropy.php";
include "Math/BAYES/NiaveBayes.php";
:

This is similiar to the goal of PEAR/Math developers and it is my hope that the PHP-based packages developed and promoted here will complement those packages. If anyone looks at the php-based math packages I have developed to date you will notice that I often adopt a miminimalist coding style that does not agree on important presentational aspects with PEAR coding standards. Some of my preferences and hangups include:

1. I prefer 2 spaces for indentation – humans have excellent vernier acuity so we arguably don’t need more space to resolve differences in code-block indentation.
2. I prefer the KR style of placing the beginning brace “{” at the end of a code line rather that on its own line.
3. I only use braces if I have to.
4. Sometimes it scans better when you have zero space between the variable and operator elements of an algebraic equation.
5. I will sometimes use one letter variable names $A to denote an important class member rather than something longer and more descriptive like $Matrix.
6. Sometimes, a math algorithm is comprehended better if it is made to physically occupy as little space as possible on a sheet of paper rather than to be too liberally loaded with line space, comments, and space-hogging variable names (i.e., flow-breaking elements).

My coding style tries to economize on the amount of physical space the code occupies on a sheet of paper if it were to be printed. I call it the “ecomony of space” coding style. The argument for using this coding style is that it is less arbitrary than other coding styles. There is a consistent economy-of-space principle behind the code aesthetics the style proposes. Given these presentational differences in coding style, I still try to follow many of the PEAR best practices when developing PHP/MATH object-oriented packages. I expect the site to also differ in some significant respects from the PEAR:MATH effort in that there may be an educational component to this site and other varieties of math-related code might be actively promoted and developed (i.e., php math function libraries, php math snippets, datamining tools, web interfaces, AI experiments, game programming components, financial analysis tools, 3D rendering, genomic analysis, natural language processing, survey analysis, etc…).