/*::.
==================================================================================================================================
=================================================¦ Copyright © 2008 Allen Baker ¦=================================================
+------------------------------+
File: PolynomialFit.java
Originator: Allen Baker (2008.05.08 16:18)
LayoutRev: 5
================================================================================================================================== */
/*.
==========================================================================================
Package
------------------------------------------------------------------------------------------ */
package cosmicabyss.com.lib;
/*.
==========================================================================================
Imports
------------------------------------------------------------------------------------------ */
import java.io.*;
import java.util.*;
import Jama.*;
import Jama.util.*;
/*::
======================================================================================================================== *//**
This class determines the polynomial that is a best fit to a set of data points. When the user instantiates an object of
this class, the user specifies the set of data points to fit a polynomial to and also specifies the order of the
polynomial.<P>
This current version implements the least-squares algorithm to determine the polynomial.
<P>
<DL>
<DT>
<B>
Example usage:
</B>
<DD>
<BLOCKQUOTE>
<PRE id="unindent">
// ===============================================
// pick the order of the polynomial you want fitted to the data points
// pick the number of data points you will have
// -----------------------------------------------
int order = 5;
int numDataPoints = 20;
// ===============================================
// these are the data points that a polynomial will be fit to. the x
// axis values and the y axis values for each data point go in these
// arrays.
// -----------------------------------------------
double[] xDataPointComponentValues = new double[numDataPoints];
double[] yDataPointComponentValues = new double[numDataPoints];
// ===============================================
// generate some data points using a known set of polynomial coefficients.
// we will test the polynomial fitting class by seeing if it can figure
// out which coefficients we used to generate the data points
// -----------------------------------------------
cOut.println();
cOut.println("Original Data Points");
for (int i=0; i<xDataPointComponentValues.length; i++)
{
xDataPointComponentValues[i] = i;
yDataPointComponentValues[i] = 0;
for (int j=0; j<order; j++)
{
int coefficientOfThisTerm = j;
yDataPointComponentValues[i] += coefficientOfThisTerm * UMath.pow(xDataPointComponentValues[i],(double)j);
}
cOut.println(xDataPointComponentValues[i] + " " + yDataPointComponentValues[i]);
}
// ===============================================
// this is where we will store the coefficients that we get back from
// the PolynomialFit class object
// -----------------------------------------------
double[] coefficients = new double[order];
// ===============================================
// create an object for the selected order of polynomial and data points
// and tell the object to return the set of coefficients it computes.
// -----------------------------------------------
PolynomialFit pf = new PolynomialFit(order,xDataPointComponentValues,yDataPointComponentValues);
coefficients = pf.coefficients();
// ===============================================
// print out the coefficients
// -----------------------------------------------
cOut.println();
cOut.println("Coefficients");
for (int j=0; j<order; j++)
{
cOut.println(coefficients[j]);
}
// ===============================================
// now regenerate the data points using the coefficients that the object
// returned and see if they are the same as the original ones.
// -----------------------------------------------
cOut.println();
cOut.println("Computed Data Points");
for (int i=0; i<xDataPointComponentValues.length; i++)
{
yDataPointComponentValues[i] = 0;
for (int j=0; j<order; j++)
{
yDataPointComponentValues[i] += coefficients[j] * UMath.pow(xDataPointComponentValues[i],(double)j);
}
cOut.println(xDataPointComponentValues[i] + " " + yDataPointComponentValues[i]);
}
</PRE>
</BLOCKQUOTE>
</DD>
</DT>
<DT>
<B>
View Source:
</B>
<DD>
<A href="PolynomialFit.java.html">
PolynomialFit.java
</A>
</DD>
</DT>
<DT>
<B>
Author:
</B>
<DD>
<A href="mailto:sourcecode.v01@cosmicabyss.com">
Allen Baker
</A>
</DD>
</DT>
</DL>
*//*
======================================================================================================================== */
public class PolynomialFit
{
/*: :METHOD:000:BOOKMARK:
====================================================================================================
[][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][]
==================================================================================================== *//**
This method creates a PolynomialFit class object. The object generates a polynomial of the given
order that fits the given data points.
<P><B>Implementation: </B><A HREF="PolynomialFit.java.html#000">View source</A>
@param
pOrder is the order of the polynomial to generate.
@param
PX is an ArrayList containing the x axis components values for each data point
@param
PY is an ArrayList containing the y axis components values for each data point
*//*
---------------------------------------------------------------------------------------------------- */
public PolynomialFit(int pOrder, ArrayList<Double> pX, ArrayList<Double> pY)
{
setOrder(pOrder);
setDataPoints(pX,pY);
computeBestFitPolynomial();
}
/*: :METHOD:001:BOOKMARK:
====================================================================================================
[][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][]
==================================================================================================== *//**
This method creates a PolynomialFit class object. The object generates a polynomial of the given
order that fits the given data points.
<P><B>Implementation: </B><A HREF="PolynomialFit.java.html#001">View source</A>
@param
pOrder is the order of the polynomial to generate.
@param
PX is an array containing the x axis components values for each data point
@param
PY is an array containing the y axis components values for each data point
*//*
---------------------------------------------------------------------------------------------------- */
public PolynomialFit(int pOrder, double[] pX, double[] pY)
{
setOrder(pOrder);
setDataPoints(pX,pY);
computeBestFitPolynomial();
}
/*: :METHOD:002:BOOKMARK:
====================================================================================================
[][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][]
==================================================================================================== *//**
This method returns the coefficients of the polynomial that this object computed.
<P><B>Implementation: </B><A HREF="PolynomialFit.java.html#002">View source</A>
@return
An array of coefficients
*//*
---------------------------------------------------------------------------------------------------- */
public double[] coefficients()
{
return iCoefficients;
}
/*: :METHOD:003:BOOKMARK:
====================================================================================================
[][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][]
==================================================================================================== *//**
This method returns the coefficients of the polynomial that this object computed.
<P><B>Implementation: </B><A HREF="PolynomialFit.java.html#003">View source</A>
@return
An ArrayList of coefficients
*//*
---------------------------------------------------------------------------------------------------- */
public ArrayList<Double> coefficientsAsList()
{
ArrayList<Double> coefficients = new ArrayList<Double>();
for (int i=0; i<iOrder; i++)
{
coefficients.add(new Double(iCoefficients[i]));
}
return coefficients;
}
/*: :METHOD:004:BOOKMARK:
====================================================================================================
[][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][]
==================================================================================================== *//**
This method returns the polynomial values that are computed using the coefficients
<P><B>Implementation: </B><A HREF="PolynomialFit.java.html#004">View source</A>
@return
An array of polynomial values
*//*
---------------------------------------------------------------------------------------------------- */
public double[] polynomialValues()
{
double[] y = new double[iNumDataPoints];
for (int i=0; i<iNumDataPoints; i++)
{
y[i] = 0;
for (int j=0; j<iOrder; j++)
{
y[i] += iCoefficients[j] * UMath.pow(iX[i],(double)j);
}
}
return y;
}
/*: :METHOD:005:BOOKMARK:
====================================================================================================
[][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][]
==================================================================================================== *//**
This method returns the polynomial values that are computed using the coefficients
<P><B>Implementation: </B><A HREF="PolynomialFit.java.html#005">View source</A>
@return
An ArrayList of polynomial values
*//*
---------------------------------------------------------------------------------------------------- */
public ArrayList<Double> polynomialValuesAsList()
{
/*.
==========================================================================================
Get the y values as an array of doubles
------------------------------------------------------------------------------------------ */
double[] y = polynomialValues();
/*.
==========================================================================================
Create an ArrayList of Double objects from the y values and return the list
------------------------------------------------------------------------------------------ */
ArrayList<Double> yList = new ArrayList<Double>();
for (int i=0; i<iNumDataPoints; i++)
{
yList.add(new Double(y[i]));
}
return yList;
}
/*: :METHOD:006:BOOKMARK:
====================================================================================================
[][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][]
==================================================================================================== *//**
This method sets the ConsoleStream for this object.
<P><B>Implementation: </B><A HREF="PolynomialFit.java.html#006">View source</A>
@return
A reference to this object
@param
pConsole is a ConsoleStream through which this objects sends its output.
*//*
---------------------------------------------------------------------------------------------------- */
public PolynomialFit setOut(ConsoleStream pConsole)
{
cOut = pConsole;
return this;
}
/*: :METHOD:007:BOOKMARK:
====================================================================================================
[][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][]
==================================================================================================== *//**
This boilerplate method is used for testing an instantiated object of this class and may include any
code the developer chooses.
<P><B>Implementation: </B><A HREF="PolynomialFit.java.html#007">View source</A>
*//*
---------------------------------------------------------------------------------------------------- */
public PolynomialFit test() throws Exception
{
cOut.titledPrintf
(
"\"HELLO WORLD!\"",
"%s %s %s",
"I'm an object of the", CLASS_NAME, "class, and I approved this message."
);
return this;
}
/*:.
==============================================================================================================
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@[ Protected ]@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
============================================================================================================== */
/*:.
==============================================================================================================
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@[ Private ]@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
============================================================================================================== */
/*.
==========================================================================================
Class Constants
<BLOCKQUOTE>
<PRE id="unindent">
CLASS_NAME : the name of this class
DFLT_LINE_LEN : the default line length for word wrapping
</PRE>
</BLOCKQUOTE>
------------------------------------------------------------------------------------------ */
private static final XString CLASS_NAME = new XString(PolynomialFit.class.getName());
private static final int DFLT_LINE_LEN = ConsoleMessage.defaultLineLength();
/*.
==========================================================================================
Class variables
cOut : console output.
------------------------------------------------------------------------------------------ */
private static ConsoleStream cOut = ConsoleStream.getSingleton();
/*.
==========================================================================================
Instance variables
<BLOCKQUOTE>
<PRE id="unindent">
iNumDataPoints : the number of data points
iOrder : the order of the polynomial
iX : the data points' x component values
iY : the data points' y component values
iCoefficients : the coefficients of the computed polynomial
iErrorOnCoefficients : the errors of the coefficients
iSigmaX : the x value error for each data point
iSigmaY : the y value error for each data point
</PRE>
</BLOCKQUOTE>
------------------------------------------------------------------------------------------ */
int iNumDataPoints = 0;
int iOrder = 0;
double[] iX = null;
double[] iY = null;
double[] iCoefficients = null;
double[] iErrorOnCoefficients = null;
double[] iSigmaX = null;
double[] iSigmaY = null;
/*: :METHOD:008:BOOKMARK:
====================================================================================================
[][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][]
==================================================================================================== *//**
This method sets the order of the polynomial.
<P><B>Implementation: </B><A HREF="PolynomialFit.java.html#008">View source</A>
@param
pOrder is the order of the polynomial
*//*
---------------------------------------------------------------------------------------------------- */
private void setOrder(int pOrder)
{
iOrder = pOrder;
iCoefficients = new double[iOrder];
iErrorOnCoefficients = new double[iOrder];
}
/*: :METHOD:009:BOOKMARK:
====================================================================================================
[][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][]
==================================================================================================== *//**
This method sets the data point (x,y) values
<P><B>Implementation: </B><A HREF="PolynomialFit.java.html#009">View source</A>
@param
PX is an array of the x values
@param
PY is an array of the y values
*//*
---------------------------------------------------------------------------------------------------- */
private void setDataPoints(double[] pX, double[] pY)
{
iNumDataPoints = pX.length;
iSigmaX = new double[iNumDataPoints];
iSigmaY = new double[iNumDataPoints];
iX = Arrays.copyOf(pX,iNumDataPoints);
iY = Arrays.copyOf(pY,iNumDataPoints);
}
/*: :METHOD:010:BOOKMARK:
====================================================================================================
[][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][]
==================================================================================================== *//**
This method sets the data point (x,y) values
<P><B>Implementation: </B><A HREF="PolynomialFit.java.html#010">View source</A>
@param
PX is an ArrayList of the x values
@param
PY is an ArrayList of the y values
*//*
---------------------------------------------------------------------------------------------------- */
private void setDataPoints(ArrayList<Double> pX, ArrayList<Double> pY)
{
iNumDataPoints = pX.size();
iSigmaX = new double[iNumDataPoints];
iSigmaY = new double[iNumDataPoints];
iX = new double[iNumDataPoints];
iY = new double[iNumDataPoints];
/*.
==========================================================================================
Convert from Double objects in a list to double primitives in an array
------------------------------------------------------------------------------------------ */
for (int i=0; i<iNumDataPoints; i++)
{
iX[i] = pX.get(i).doubleValue();
iY[i] = pY.get(i).doubleValue();
}
}
/*: :METHOD:011:BOOKMARK:
====================================================================================================
[][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][]
==================================================================================================== *//**
This method computes the polynomial that fits the data points.
<P><B>Implementation: </B><A HREF="PolynomialFit.java.html#011">View source</A>
*//*
---------------------------------------------------------------------------------------------------- */
private void computeBestFitPolynomial()
{
/*.
==========================================================================================
NumParams = num coeff + error on each coeff.
------------------------------------------------------------------------------------------ */
double[][] alpha = new double[iOrder][iOrder];
double[] beta = new double[iOrder];
double term = 0;
/*.
==========================================================================================
------------------------------------------------------------------------------------------ */
for (int k=0; k<iOrder; k++)
{
/*.
==========================================================================================
Only need to calculate diagonal and upper half of symmetric matrix.
------------------------------------------------------------------------------------------ */
for (int j=k; j<iOrder; j++)
{
/*.
==========================================================================================
Calc terms over the data points
------------------------------------------------------------------------------------------ */
term = 0.0;
alpha[k][j] = 0.0;
for (int i=0; i<iNumDataPoints; i++)
{
/*.
==========================================================================================
Calculate x^k
------------------------------------------------------------------------------------------ */
double prod1 = UMath.pow(iX[i],k);
/*.
==========================================================================================
Calculate x^j
------------------------------------------------------------------------------------------ */
double prod2 = UMath.pow(iX[i],j);
/*.
==========================================================================================
Calculate x^k * x^j
------------------------------------------------------------------------------------------ */
term = (prod1*prod2);
/*.
==========================================================================================
------------------------------------------------------------------------------------------ */
if (iSigmaY != null && iSigmaY[i] != 0.0)
{
term /= (iSigmaY[i]*iSigmaY[i]);
}
alpha[k][j] += term;
}
alpha[j][k] = alpha[k][j]; // C will need to be inverted.
}
/*.
==========================================================================================
------------------------------------------------------------------------------------------ */
for (int i=0; i < iNumDataPoints; i++)
{
double prod1 = UMath.pow(iX[i],k);
term = (iY[i] * prod1);
if (iSigmaY != null && iSigmaY[i] != 0.0)
{
term /= (iSigmaY[i]*iSigmaY[i]);
}
beta[k] += term;
}
}
/*.
==========================================================================================
Use the Jama QR Decomposition classes to solve for the parameters.
------------------------------------------------------------------------------------------ */
Matrix alpha_matrix = new Matrix(alpha);
QRDecomposition alpha_QRD = new QRDecomposition(alpha_matrix);
Matrix beta_matrix = new Matrix(beta,iOrder);
Matrix param_matrix;
try
{
param_matrix = alpha_QRD.solve(beta_matrix);
}
catch (Exception e)
{
cOut.println(Const.WARNING,"QRD solve failed: " + e);
return;
}
/*.
==========================================================================================
The inverse provides the covariance matrix.
------------------------------------------------------------------------------------------ */
Matrix c = alpha_matrix.inverse();
for (int k=0; k < iOrder; k++)
{
iCoefficients[k] = param_matrix.get(k,0);
/*.
==========================================================================================
Diagonal elements of the covariance matrix provide the square of the parameter errors. Put
in top half of the parametes array.
------------------------------------------------------------------------------------------ */
iErrorOnCoefficients[k] = UMath.sqrt(c.get(k,k));
}
}
/*:.
==============================================================================================================
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@[ Inner Classes ]@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
============================================================================================================== */
/*:.
==============================================================================================================
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@[ Public Static Methods ]@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
============================================================================================================== */
/*: :METHOD:012:BOOKMARK:
====================================================================================================
[][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][]
==================================================================================================== *//**
This method allows this class file to be unit tested as a standalone application. It's the method
that's called when the class is invoked from the command line by using the java application launcher
- "java". Main() is not a required method, but the practice of putting one in each class and
wrapping class test code within it allows easy unit testing of the class; and main does not need to
be removed when testing is complete.
<P>
<DL>
<DT>
<B>
Command line usage:
</B>
<DD>
Java cosmicabyss.com.lib.PolynomialFit
</DD>
</DT>
</DL>
<P><B>Implementation: </B><A HREF="PolynomialFit.java.html#012">View source</A>
@param
pArgs contains the command line arguments with which this class was invoked as an application.
*//*
---------------------------------------------------------------------------------------------------- */
public static void main(String[] pArgs) throws Exception
{
/*.
==========================================================================================
Greetings !
------------------------------------------------------------------------------------------ */
cOut.banner(CLASS_NAME);
/*.
==========================================================================================
Pick the order of the polynomial you want fitted to the data points pick the number of
data points you will have
------------------------------------------------------------------------------------------ */
int order = 5;
int numDataPoints = 367;
/*.
==========================================================================================
These are the data points that a polynomial will be fit to. The x axis values and the y
axis values for each data point go in these arrays.
------------------------------------------------------------------------------------------ */
double[] xDataPointComponentValues = new double[numDataPoints];
double[] yDataPointComponentValues = new double[numDataPoints];
/*.
==========================================================================================
Generate some data points using a known set of polynomial coefficients. We will test the
polynomial fitting class by seeing if it can figure out which coefficients we used to
generate the data points
------------------------------------------------------------------------------------------ */
cOut.println();
cOut.println("Original Data Points");
for (int i=0; i<xDataPointComponentValues.length; i++)
{
xDataPointComponentValues[i] = i;
yDataPointComponentValues[i] = 0;
for (int j=0; j<order; j++)
{
int coefficientOfThisTerm = j;
yDataPointComponentValues[i] += coefficientOfThisTerm * UMath.pow(xDataPointComponentValues[i],(double)j);
}
// cOut.println(xDataPointComponentValues[i] + " " + yDataPointComponentValues[i]);
}
/*.
==========================================================================================
This is where we will store the coefficients that we get back from the PolynomialFit class
object
------------------------------------------------------------------------------------------ */
double[] coefficients = new double[order];
/*.
==========================================================================================
Create an object for the selected order of polynomial and data points and tell the object
to return the set of coefficients it computes.
------------------------------------------------------------------------------------------ */
PolynomialFit pf = new PolynomialFit(order,xDataPointComponentValues,yDataPointComponentValues);
coefficients = pf.coefficients();
/*.
==========================================================================================
Print out the coefficients
------------------------------------------------------------------------------------------ */
cOut.println();
cOut.println("Computed Coefficients");
for (int j=0; j<order; j++)
{
cOut.println(j + " " + UMath.round(coefficients[j],0));
}
/*.
==========================================================================================
Now regenerate the data points using the coefficients that the object returned and see if
they are the same as the original ones.
------------------------------------------------------------------------------------------ */
cOut.println();
cOut.println("Computed Data Points");
for (int i=0; i<xDataPointComponentValues.length; i++)
{
yDataPointComponentValues[i] = 0;
for (int j=0; j<order; j++)
{
yDataPointComponentValues[i] += coefficients[j] * UMath.pow(xDataPointComponentValues[i],(double)j);
}
// cOut.println(xDataPointComponentValues[i] + " " + yDataPointComponentValues[i]);
}
}
} // class PolynomialFit
/*:.
==============================================================================================================
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@[ Notes ]@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
============================================================================================================== */
Copyright © 2012 Allen Baker
| 
Facebook