io.yields.math.concepts.operator.Smoothness.java Source code

Introduction

Here is the source code for io.yields.math.concepts.operator.Smoothness.java

Source

/*
 * Copyright 2015 by Yields.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *        http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package io.yields.math.concepts.operator;

import org.apache.commons.lang.Validate;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.DecompositionSolver;
import org.apache.commons.math3.linear.QRDecomposition;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.stat.descriptive.DescriptiveStatistics;

import java.util.Collection;
import java.util.DoubleSummaryStatistics;
import java.util.List;
import java.util.function.Function;
import java.util.stream.Collectors;
import java.util.stream.IntStream;

/**
 * Matches a polynomial of given order to the normalized data
 * and calculates the absolute distance to the polynome in the least square sense
 *
 * This is a measure of the smoothness of the data
 */
public class Smoothness implements Function<Collection<Tuple>, DescriptiveStatistics> {

    private static final double EPS = 1.e-16;

    private int order;

    public Smoothness(int order) {
        Validate.isTrue(order >= 1);
        this.order = order;
    }

    @Override
    public DescriptiveStatistics apply(Collection<Tuple> tuples) {
        Validate.isTrue(tuples.size() > order);
        //first we normalize the tuples so data fits in the unit square
        List<Tuple> normalizedData = normalize(tuples);
        //calculate error (i.e. distance between tuple y and fitted polynomial
        RealMatrix error = computeDistance(normalizedData);
        //group in stats object
        DescriptiveStatistics stats = new DescriptiveStatistics();
        for (double value : error.getColumn(0)) {
            stats.addValue(Math.abs(value));
        }
        return stats;
    }

    public static List<Tuple> normalize(Collection<Tuple> tuples) {
        DoubleSummaryStatistics statsForX = tuples.stream()
                .collect(Collectors.summarizingDouble(tuple -> tuple.getX()));
        DoubleSummaryStatistics statsForY = tuples.stream()
                .collect(Collectors.summarizingDouble(tuple -> tuple.getY()));
        double minX = statsForX.getMin();
        double maxX = statsForX.getMax();
        double minY = statsForY.getMin();
        double maxY = statsForY.getMax();
        return tuples.stream().map(tuple -> {
            double x = (tuple.getX() - minX) / (maxX - minX);
            double y = (tuple.getY() - minY) / (maxY - minY);
            return new Tuple(x, y);
        }).collect(Collectors.toList());
    }

    private RealMatrix computeDistance(List<Tuple> normalizedData) {
        RealMatrix a = new Array2DRowRealMatrix(normalizedData.size(), order + 1);
        RealMatrix b = new Array2DRowRealMatrix(normalizedData.size(), 1);
        int row = 0;
        for (Tuple tuple : normalizedData) {
            final int currentRow = row;
            IntStream.rangeClosed(0, order)
                    .forEach(power -> a.setEntry(currentRow, power, Math.pow(tuple.getX(), power)));
            b.setEntry(currentRow, 0, tuple.getY());
            row++;
        }
        DecompositionSolver solver = new QRDecomposition(a, EPS).getSolver();
        if (solver.isNonSingular() && !isConstant(b)) {
            RealMatrix solution = solver.solve(b);
            return a.multiply(solution).subtract(b);
        } else {
            return new Array2DRowRealMatrix(normalizedData.size(), 1);
        }
    }

    private boolean isConstant(RealMatrix matrix) {
        double refElement = matrix.getEntry(0, 0);
        return !IntStream.range(0, matrix.getRowDimension()).asDoubleStream()
                .filter(element -> Math.abs(refElement - element) > EPS).findAny().isPresent();
    }
}