public
Last active

  • Download Gist
PolynomialFitter.java
Java
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213
/***************************************************************************
* Copyright (C) 2009 by Paul Lutus, Ian Clarke *
* lutusp@arachnoid.com, ian.clarke@gmail.com *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program; if not, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
 
import java.util.ArrayList;
import java.util.List;
 
/**
* A class to fit a polynomial to a (potentially very large) dataset.
*
* @author Paul Lutus <lutusp@arachnoid.com>
* @author Ian Clarke <ian.clarke@gmail.com>
*
*/
 
/*
* Changelog:
* 20100130: Add note about relicensing
* 20091114: Modify so that points can be added after the curve is
* created, also some other minor fixes
* 20091113: Extensively modified by Ian Clarke, main changes:
* - Should now be able to handle extremely large datasets
* - Use generic Java collections classes and interfaces
* where possible
* - Data can be fed to the fitter as it is available, rather
* than all at once
*
* The code that this is based on was obtained from: http://arachnoid.com/polysolve
*
* Note: I (Ian Clarke) am happy to release this code under a more liberal
* license such as the LGPL, however Paul Lutus (the primary author) refuses
* to do this on the grounds that the LGPL is not an open source license.
* If you want to try to explain to him that the LGPL is indeed an open
* source license, good luck - it's like talking to a brick wall.
*/
public class PolynomialFitter {
 
private final int p, rs;
 
private long n = 0;
 
private final double[][] m;
 
private final double[] mpc;
/**
* @param degree
* The degree of the polynomial to be fit to the data
*/
public PolynomialFitter(final int degree) {
assert degree > 0;
p = degree + 1;
rs = 2 * p - 1;
m = new double[p][p + 1];
mpc = new double[rs];
}
 
/**
* Add a point to the set of points that the polynomial must be fit to
*
* @param x
* The x coordinate of the point
* @param y
* The y coordinate of the point
*/
public void addPoint(final double x, final double y) {
assert !Double.isInfinite(x) && !Double.isNaN(x);
assert !Double.isInfinite(y) && !Double.isNaN(y);
n++;
// process precalculation array
for (int r = 1; r < rs; r++) {
mpc[r] += Math.pow(x, r);
}
// process RH column cells
m[0][p] += y;
for (int r = 1; r < p; r++) {
m[r][p] += Math.pow(x, r) * y;
}
}
 
/**
* Returns a polynomial that seeks to minimize the square of the total
* distance between the set of points and the polynomial.
*
* @return A polynomial
*/
public Polynomial getBestFit() {
final double[] mpcClone = mpc.clone();
final double[][] mClone = new double[m.length][];
for (int x = 0; x < mClone.length; x++) {
mClone[x] = m[x].clone();
}
 
mpcClone[0] += n;
// populate square matrix section
for (int r = 0; r < p; r++) {
for (int c = 0; c < p; c++) {
mClone[r][c] = mpcClone[r + c];
}
}
gj_echelonize(mClone);
final Polynomial result = new Polynomial(p);
for (int j = 0; j < p; j++) {
result.add(j, mClone[j][p]);
}
return result;
}
private double fx(final double x, final List<Double> terms) {
double a = 0;
int e = 0;
for (final double i : terms) {
a += i * Math.pow(x, e);
e++;
}
return a;
}
private void gj_divide(final double[][] A, final int i, final int j, final int m) {
for (int q = j + 1; q < m; q++) {
A[i][q] /= A[i][j];
}
A[i][j] = 1;
}
 
private void gj_echelonize(final double[][] A) {
final int n = A.length;
final int m = A[0].length;
int i = 0;
int j = 0;
while (i < n && j < m) {
// look for a non-zero entry in col j at or below row i
int k = i;
while (k < n && A[k][j] == 0) {
k++;
}
// if such an entry is found at row k
if (k < n) {
// if k is not i, then swap row i with row k
if (k != i) {
gj_swap(A, i, j);
}
// if A[i][j] is not 1, then divide row i by A[i][j]
if (A[i][j] != 1) {
gj_divide(A, i, j, m);
}
// eliminate all other non-zero entries from col j by
// subtracting from each
// row (other than i) an appropriate multiple of row i
gj_eliminate(A, i, j, n, m);
i++;
}
j++;
}
}
 
private void gj_eliminate(final double[][] A, final int i, final int j, final int n, final int m) {
for (int k = 0; k < n; k++) {
if (k != i && A[k][j] != 0) {
for (int q = j + 1; q < m; q++) {
A[k][q] -= A[k][j] * A[i][q];
}
A[k][j] = 0;
}
}
}
 
private void gj_swap(final double[][] A, final int i, final int j) {
double temp[];
temp = A[i];
A[i] = A[j];
A[j] = temp;
}
 
 
public static class Polynomial extends ArrayList<Double> {
private static final long serialVersionUID = 1692843494322684190L;
 
public Polynomial(final int p) {
super(p);
}
 
public double getY(final double x) {
double ret = 0;
for (int p=0; p<size(); p++) {
ret += get(p)*(Math.pow(x, p));
}
return ret;
}
 
@Override
public String toString() {
final StringBuilder ret = new StringBuilder();
for (int x = size() - 1; x > -1; x--) {
ret.append(get(x) + (x > 0 ? "x^" + x + " + " : ""));
}
return ret.toString();
}
}
}

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.