image
 
image
PolynomialFit.java


/*::.
==================================================================================================================================
=================================================¦ 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&lt;xDataPointComponentValues.length; i++)
                     {
                     xDataPointComponentValues[i] = i;
                     yDataPointComponentValues[i] = 0;
                     for (int j=0; j&lt;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&lt;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&lt;xDataPointComponentValues.length; i++)
                     {
                     yDataPointComponentValues[i] = 0;
                     for (int j=0; j&lt;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  ]@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
   ============================================================================================================== */