package polyofdm.matrixgen;
import java.util.Arrays;
import Jama.Matrix;
public class Newton {
public static void main(String[] args) {
int all = 24;
int missed = 12;
generateMatrix(all, missed);
}
public static Matrix generateMatrix(int all, int missed) {
double[] weights = Utils.createWeightings(all);
System.out.println("denominator-weights: " + Arrays.toString(weights));
double[] roots = findOrthogonalRoots(all, missed, weights);
System.out.println("Found orthogonal Roots: " + Arrays.toString(roots));
Matrix matrix = new Matrix(all-missed, all);
for (int poly = 0; poly < matrix.getRowDimension(); poly++) {
for (int x = 0; x < matrix.getColumnDimension(); x++) {
matrix.set(poly, x, calcProduct(poly, x, roots)*weights[x]);
}
}
// Utils.makeOrthoGramSchmidt(matrix);
Utils.normalize(matrix);
System.out.println("Found Matrix");
matrix.print(15, 5);
return matrix;
}
private static double calcProduct(int poly, int x, double[] roots) {
double product = 1;
for (int i = 0; i < roots.length; i++) {
if (i!=poly){
product = product * (x-roots[i]);
}
}
// Only for better chart Design
// All Max values up
// Does not effect orthogonality
if ((roots.length-poly)%2==1)
return -product;
else
return product;
}
private static double[] findOrthogonalRoots(int all, int missed,
double[] weights) {
double additionalSpace = 1;
double step = ((double)all-additionalSpace)/((double)(all-missed));
double offset = ((all-1) - ((all-missed-1) * step))/2d;
double[] roots = new double[all-missed];
for (int i=0; i< all-missed;i++) {
//System.out.println("i=" + i + " pos=" + i*step);
roots[i]=offset + i*step;
}
for (int i = 0; i < 5; i++) {
roots = moveToCenter(roots, weights);
roots = doMinimize(roots, weights);
}
return roots;
}
private static double[] doMinimize(double[] rootsxx2, double[] weights) {
OrthogonalityNorm orthoNorm = new OrthogonalityNorm(rootsxx2, weights);
int counter = 0;
do {
counter++;
// System.out.println("counter=" + counter);
// System.out.println("orthoNorm.getMeasure()=" + orthoNorm.getMeasure());
double[] roots = new double[orthoNorm.getRoots().length];
double gradientNormSquared = 0;
for (int i = 0; i < orthoNorm.getDeviations().length; i++) {
gradientNormSquared += orthoNorm.getDeviations()[i] * orthoNorm.getDeviations()[i];
}
for (int i = 0; i < roots.length; i++) {
roots[i] = orthoNorm.getRoots()[i] - orthoNorm.getMeasure()*orthoNorm.getDeviations()[i] / gradientNormSquared;
}
orthoNorm = new OrthogonalityNorm(roots, weights);
} while (orthoNorm.getMeasure() > 1E-28);
System.out.println("##### Newton.doMinimize() finisched ####### Rounds=" + counter + " Orthogonality=" + orthoNorm.getMeasure());
return orthoNorm.getRoots();
}
private static double[] moveToCenter(double[] roots, double[] weights) {
double sum = 0;
for (double d : roots) {
sum += d;
}
double avg = sum / roots.length;
double shouldBe = (weights.length-1)/2d;
double diff = avg-shouldBe;
System.out.println("moveToCenter avg=" + avg + " shouldBe="+ shouldBe+ " moving all roots with diff="+ diff);
double[] newRoots = new double[roots.length];
for (int i = 0; i < newRoots.length; i++) {
newRoots[i] = roots[i] -diff;
}
return newRoots;
}
}