PoincareHalfSpaceIdealCenter.prejava

/*
 * Given a set of points
 * on the boundary plane of the poincare half-space model,
 * find the point in the poincare half-space
 * such that, when the whole thing is transformed
 * into the poincare ball with that point at the origin,
 * the euclidean centroid of the points
 * is the euclidean center of the ball.
 */
#include "macros.h"
import com.donhatchsw.util.VecMath;
import com.donhatchsw.util.NewtonSolver;
import com.donhatchsw.util.MyMath;
import com.donhatchsw.util.Arrays;
public class PoincareHalfSpaceIdealCenter
{
    // vs are n-dimensional verts on the plane-at-infinity
    // in the n+1-dimensional poincare half-space model.
    // Compute the point in the n+1-dimensional poincare half-space model
    // that, when centered in the n+1-dimensional poincare disk model,
    // centers the euclidean average of the verts.
    public static double[] PoincareHalfSpaceIdealCenter(double vs[][])
    {
        System.out.println("        in pseudoCentroid");

        //
        // Initial guess
        // is a generalization
        // of the in-center of 3 ideal points:
        // Take a weighted average of the vertices,
        // each vertex weighted by the variance
        // of the other n-1 vertices,
        // and offset from the plane
        // by the weighted average (same weights) of the standard deviations
        // (I think this is roughly right).
        // This immediately gives the right answer for an ideal triangle, anyway.
        // XXX can we modify to give the right answer for an ideal regular simplex??
        //

        final int n = vs.length;
        final int nDims = vs[0].length;
        // not, dim of answer will be nDims+1

        // M = mean
        double M[] = VecMath.average(vs);
        // S = n times variance, in each dimension separately
        double S[] = new double[nDims]; // zeros
        FORI (i, n)
        {
            double v[] = vs[i];
            FORI (iDim, nDims)
                S[iDim] += SQR(v[iDim]-M[iDim]);
        }

        // We want the variance of each subset
        // of n-1 points of vs.
        // From http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
        // To compute the new S = n*newVariance
        // given s = (n-1)*oldVariance and a new sample x:
        //     S = s + (x-m)*(x-M)
        // so,
        //     s = S - (x-m)*(x-M)
        // where m and M are the old and new means:
        //   M = m + (x-m)/n
        //     = (1-1/n)*m + (1/n)*x
        // so,
        //   m = (M-x/n)/((n-1)/n)
        //     = (M*n-x)/(n-1)
        // The simplest way to think about it
        // is to do each dimension separately;
        // the full variance is the sum
        // of the variances in each dimension.
        double weights[] = new double[n];
        FORI (i, n)
        {
            double v[] = vs[i];
            double s = 0.;
            FORI (iDim, nDims)
            {
                double miDim = (M[iDim]*n-v[iDim])/(n-1.);
                double siDim = S[iDim] - (v[iDim]-miDim)*(v[iDim]-M[iDim]);
                s += siDim;
            }
            weights[i] = s / (n-1.);
        }
        PRINTARRAY(weights);
        VecMath.vxs(weights, weights, 1./VecMath.sum(weights));
        PRINTARRAY(weights);

        double weightedAvg[] = VecMath.vxm(weights, vs);
        PRINTARRAY(weightedAvg);

        double sum = 0.;
        FORI (i, n)
            sum += weights[i] * VecMath.distsqrd(vs[i], weightedAvg);
        double height = Math.sqrt(sum);

        // initial guess is weightedAvg with height appended
        double initialGuess[] = Arrays.append(weightedAvg, height);
        PRINTARRAY(initialGuess);


        // vs0 will be vs with 0 appended to each vertex
        final double vs0[][] = new double[n][];
        FORI (i, n)
            vs0[i] = Arrays.append(vs[i], 0.);


        // For the solve, solve for log(height) rather than height.
        double logInitialGuess[] = VecMath.copyvec(initialGuess);
        logInitialGuess[nDims] = Math.log(initialGuess[nDims]);
        PRINTARRAY(logInitialGuess);

        final boolean doItInLogSpace = true;
        if (!doItInLogSpace)
            logInitialGuess = initialGuess;

        final double eps = 1e-6; // for finite difference
        NewtonSolver.Fun fun = new NewtonSolver.Fun(nDims+1) {

            public void f(double logGuess[], double answerImagesSum[])
            {
                //System.out.println("            in f");
                //System.out.print("                "); PRINTARRAY(logGuess);


                double guess[] = VecMath.copyvec(logGuess);
                guess[nDims] = Math.exp(logGuess[nDims]);
                if (!doItInLogSpace)
                    guess = logGuess;
                //System.out.print("                "); PRINTARRAY(guess);

                double reflectionRadius = guess[nDims];
                CHECK_GT(reflectionRadius, 0); // guaranteed since it's the exp of something... if we were given it directly, then newton might overshoot, which would be a disaster
                // (image poincare disk radius is half that)

                double reflectionCenter[] = VecMath.copyvec(guess);
                reflectionCenter[nDims] = -reflectionRadius;

                // The following inverts the reflection circle
                // and scales it down to unit size at the origin:
                //       v -> reflectionRadius * (v-reflectionCenter)/(v-reflectionCenter).length2()
                // which means the poincare half-plane
                // will get mapped to a disk of radius 1/2 centered at 0,0,...,1/2.
                // To turn that into a disk centered at 0,
                // we then subtract 0,0,...,1/2.
                double image[] = new double[nDims+1]; // scratch
                FORI (i, n)
                {
                    double v0[] = vs0[i];
                    VecMath.vmv(image, v0, reflectionCenter); // image = v - reflectionCenter
                    VecMath.vxs(image, image, reflectionRadius/VecMath.normsqrd(image));
                    // image is now on circle of radius 1/2 centered at 0,0,...,1/2
                    image[nDims] -= .5;
                    // image is now on circle of radius 1/2 centered at origin
                    VecMath.vpv(answerImagesSum, answerImagesSum, image);
                }
                //System.out.print("                "); PRINTARRAY(answerImagesSum);
                //System.out.print("            out f");
            } // f()
            // used in finite difference computation of jacobian
            public double eps()
            {
                //System.out.println("        in eps");
                //System.out.println("        out eps");
                return eps;
            }
        }; // fun



        int minIterations = 10; // probably not needed, but used to be the default and I haven't tested this since it was made explicit
        int maxIterations = 20;
        double logAnswer[] = VecMath.copyvec(logInitialGuess);
        PRINT(logAnswer.length);
        PRINT(logInitialGuess.length);
        PRINTARRAY(logInitialGuess);
        NewtonSolver.solve(logAnswer,
                           new double[nDims+1], // target is zero
                           fun,
                           minIterations,
                           maxIterations,
                           false); // adaptiveFlag XXX could experiment with this if needed
        double answer[] = VecMath.copyvec(logAnswer);
        answer[nDims] = Math.exp(logAnswer[nDims]);
        if (!doItInLogSpace)
            answer = logAnswer;
        return answer;
    } // PoincareHalfSpaceIdealCenter


    // little test program
    public static void main(String args[])
    {
        System.out.println("in main");
        {
            double vs[][] = {
                {-1.},
                {0.},
                {2.},
            };
            PRINTARRAY(vs);
            double center[] = PoincareHalfSpaceIdealCenter(vs);
            PRINTARRAY(center);
        }
        {
            double vs[][] = {
                {-1.,.3,4},
                {.4,.6,.7},
                {.22,.743,.37},
            };
            PRINTARRAY(vs);
            double center[] = PoincareHalfSpaceIdealCenter(vs);
            PRINTARRAY(center);
        }
        System.out.println("out main");
    } // main
} // public class PoincareHalfSpaceIdealCenter