edu.mines.jtk.opt.CoordinateTransform Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of edu-mines-jtk Show documentation
Java packages for science and engineering
The newest version!
/****************************************************************************
Copyright 2003, Landmark Graphics and others.
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 edu.mines.jtk.opt;

import java.util.Arrays;
import java.util.logging.Logger;

/**
 * Find a best linear combination of input
 * coordinates to fit output coordinates.
 * Finds a_io (a with subscripts i and o)
 * to best approximate the linear transform.
 *  * out_o = sum_i ( a_io * in_i )
 * 
 * where in_i are input coordinates,

 * out_o are output coordinates,

 * i is the index of each input dimension,

 * and o is the index of each output dimension.

 * 
 * The optimum coefficients minimize this least
 * squares error:
 *  * sum_oj [ sum_i ( a_io * in_ij ) - out_oj ]^2
 * 
 * where in_ij is an input array of different coordinates,

 * out_oj is an output array of corresponding coordinates,

 * and j is the index of pairs of coordinates to be fit
 * in a least-squares sense.
 * 
 * Normal equations (indexed by k) are solved independently for each o:
 *  * sum_ij ( in_kj * in_ij * a_io ) = sum_j ( in_kj * out_oj )
 * 
 * The Hessian is  H = sum_j ( in_kj * in_ij ) 
 * and the gradient  b = - sum_j ( in_kj * out_oj ) 
 * 
 * The solution is undamped and may not behave as you
 * want for degenerate solutions.
 * 
 * If the linear transform needs a translation shift,
 * then include a constant as one of the input coordinates.
 *
 * @author W.S. Harlan
 */
public class CoordinateTransform {
    private static final Logger LOG = Logger.getLogger("edu.mines.jtk.opt");

    private int _nout = 0;
    private int _nin = 0;
    private final double[][] _hessian;
    private final double[][] _b;
    private double[][] _a;
    private double[] _in0 = null;
    private double[] _out0 = null;
    private double[] _inScr = null;
    private double[] _outScr = null;

    /**
     * Constructor sets number of input and output coordinates.
     *
     * @param dimensionOut Number of output coordinates.
     * @param dimensionIn  Number of input coordinates.
     */
    public CoordinateTransform(final int dimensionOut, final int dimensionIn) {
        _nout = dimensionOut;
        _nin = dimensionIn;
        _hessian = new double[_nin][_nin];
        _b = new double[_nout][_nin];
        _inScr = new double[_nin];
        _outScr = new double[_nout];
    }

    /**
     * Add an observation of a set of input and output coordinates
     * You should add enough of these to determine (or overdetermine)
     * a unique linear mapping.
     * To allow translation, include a constant 1 as an input coordinate.
     *
     * @param out A set of observed output coordinates
     *            with an unknown linear relationship to input coordinates.
     * @param in  A set of observed input coordinates
     *            that should be linearly combined to calculate each of the
     *            output coordinates.
     *            To allow translation, include a constant 1.
     */
    public void add(final double[] out, final double[] in) {
        _a = null; // must redo any previous solution

        if (in.length != _nin) {
            throw new IllegalArgumentException("in must have dimension " + _nin);
        }

        if (out.length != _nout) {
            throw new IllegalArgumentException("out must have dimension " + _nout);
        }

        if (_in0 == null) {
            _in0 = in.clone();
        }
        if (_out0 == null) {
            _out0 = out.clone();
        }

        for (int i = 0; i < _nin; ++i) {
            _inScr[i] = in[i] - _in0[i];
        }
        for (int i = 0; i < _nout; ++i) {
            _outScr[i] = out[i] - _out0[i];
        }

        for (int k = 0; k < _nin; ++k) {
            for (int i = 0; i < _nin; ++i) {
                _hessian[k][i] += _inScr[k] * _inScr[i];
            }
            for (int o = 0; o < _nout; ++o) {
                _b[o][k] -= _outScr[o] * _inScr[k];
            }
        }
    }

    /**
     * For a given set of input coordinates,
     * return the linearly predicted output coordinates.
     *
     * @param in A set of input coordinates
     * @return A computed set of output coordinates.
     */
    public double[] get(double[] in) {
        for (int i = 0; i < in.length; ++i) {
            _inScr[i] = in[i] - _in0[i];
        }
        in = null;

        if (_a == null) {
            _a = new double[_nout][_nin];
            for (int o = 0; o < _nout; ++o) {
                final LinearQuadratic lq = new LinearQuadratic(o);
                final QuadraticSolver qs = new QuadraticSolver(lq);
                final ArrayVect1 solution = (ArrayVect1) qs.solve(_nin + 4, null);
                final double[] data = solution.getData();
                System.arraycopy(data, 0, _a[o], 0, _nin);
                solution.dispose();
            }
        }
        final double[] result = new double[_nout];
        for (int o = 0; o < _nout; ++o) {
            for (int i = 0; i < _nin; ++i) {
                result[o] += _a[o][i] * _inScr[i];
            }
        }
        for (int i = 0; i < result.length; ++i) {
            result[i] = result[i] + _out0[i];
        }
        return result;
    }

    // describes normal equations.
    private class LinearQuadratic implements Quadratic {
        private int _o = -1;

        /**
         * Constructor for normal equations
         *
         * @param o Index of output dimension
         */
        private LinearQuadratic(final int o) {
            _o = o;
        }

        @Override
        public void multiplyHessian(final Vect x) {
            final ArrayVect1 m = (ArrayVect1) x;
            final double[] data = m.getData();
            final double[] oldData = data.clone();
            Arrays.fill(data, 0.0);
            for (int i = 0; i < data.length; ++i) {
                for (int j = 0; j < data.length; ++j) {
                    data[i] += _hessian[i][j] * oldData[j];
                }
            }
        }

        @Override
        public Vect getB() {
            return new ArrayVect1(_b[_o].clone(), 1.0);
        }

        @Override
        public void inverseHessian(final Vect x) {
        } // not necessary
    }

}