// 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.

/*
 * MeetJoin.java
 *
 * Created on October 12, 2005, 2:10 PM
 *
 * Copyright 2005-2007 Daniel Fontijne, University of Amsterdam
 * fontijne@science.uva.nl
 *
 */

package subspace.basis;

/**
 * Computes the meet and join.
 * Usage:
 * new MeetJoin(a, b).getMeet();
 * Or:
 * MeetJoin MJ = new MeetJoin(a, b);
 * Multivector M = MJ.getMeet();
 * Multivector J = MJ.getJoin();
 *
 * @author  fontijne
 */
public class MeetJoin implements InnerProductTypes {
    protected Multivector m_meet;
    protected Multivector m_join;


    public static void main(String args[]) {
        // test code for meet / join
        int dim = 7;

        for (int i = 0; i < 1000; i++) {
            Multivector a = Multivector.getRandomBlade(dim, (int)(0.6 * Math.random() * dim), 1.0);
            Multivector b = Multivector.getRandomBlade(dim, (int)(0.6 * Math.random() * dim), 1.0);

            if (Math.random() < 0.25) {
            // make a subspace of b:
            b = a.op(Multivector.getRandomBlade(dim, (int)(0.5 * Math.random() * dim + 0.99), 1.0));
            }

            if (Math.random() < 0.25) {
            // use basis elements for 'a'
            a = new Multivector(new BasisBlade((int)(Math.random() * ((1 << 7) - 1)), 1.0));
            }
            if (Math.random() < 0.25) {
            // use basis elements for 'b'
            b = new Multivector(new BasisBlade((int)(Math.random() * ((1 << 7) - 1)), 1.0));
            }

            if (Math.random() < 0.5) {
            // swap a & b:
            Multivector tmp = b;
            b = a;
            a = tmp;
            }

            if (a.isNull(1e-3)) continue;
            if (b.isNull(1e-3)) continue;

    //      System.out.println("a = " + a + ",");
    //      System.out.println("b = " + b + ",");

            MeetJoin MJ = null;
            try {
            MJ = new MeetJoin(a, b);
            } catch (Exception ex) {
            MJ = new MeetJoin(a, b);
            }

            Multivector T;
            try {
            if ((T = MJ.getMeet().ip(a, LEFT_CONTRACTION).compress(1e-7)).isNull())
                throw new java.lang.Exception("Meet is not subspace of 'a'");
            if ((T = MJ.getMeet().ip(b, LEFT_CONTRACTION).compress(1e-7)).isNull())
                throw new java.lang.Exception("Meet is not subspace of 'b'");

            if ((T = a.ip(MJ.getJoin(), LEFT_CONTRACTION).compress(1e-7)).isNull())
                throw new java.lang.Exception("'a' is not subspace of join");
            if ((T = b.ip(MJ.getJoin(), LEFT_CONTRACTION).compress(1e-7)).isNull())
                throw new java.lang.Exception("'b' is not subspace of join");
            } catch (Exception ex)  {
            MJ = new MeetJoin(a, b);
            }

            System.out.println("Loop " + i + " ");
        }

    }

    /** Creates a new instance of MeetJoin */
    public MeetJoin(Multivector a, Multivector b) {
        computeMeetJoin(a, b);
    }

    public Multivector getMeet() {
        return m_meet;
    }
    public Multivector meet() {
        return m_meet;
    }

    public Multivector getJoin() {
        return m_join;
    }
    public Multivector join() {
        return m_join;
    }

    /** @return the delta product of 'a' and 'b' */
    public static Multivector deltaProduct(Multivector a, Multivector b) {
        Multivector D = a.gp(b);
        D = D.compress();
        return D.extractGrade(D.topGradeIndex());
    }


    /**
     * <p>Both meet and join are computed simultaneously by this
     * algorithm by Ian Bell, which turns out to be somewhat
     * more efficient than computing the join and deriving the
     * meet from that.
     *
     * @return the meet and join of this with 'b'
     * The first array element returned is the meet (intersection),
     * the second array element is the join (union).
     */
    protected void computeMeetJoin(Multivector a, Multivector b) {

        double la = a.largestCoordinate();
        double lb = b.largestCoordinate();

        double smallEpsilon = 10e-9;
        double largeEpsilon = 10e-4;

        // step one: check for near-zero input
        if ((la < smallEpsilon) || (lb < smallEpsilon)) {
            m_meet = new Multivector();
            m_join = new Multivector();
        }

        // check grade of input
        int ga = a.grade();
        int gb = b.grade();
        if (ga < 0) { // if we are not handed homogeneous multivectors, take the grade parts with the largest norm
            a = a.largestGradePart();
            ga = a.grade();
        }
        if (gb < 0) {
            b = b.largestGradePart();
            gb = b.grade();
        }

        // normalize (approximately) and swap (optionally)
        Multivector ca, cb;
        if (ga <= gb) {
            // just normalize:
            ca = a.gp(1.0 / la);
            cb = b.gp(1.0 / lb);
        }
        else {
            // also swap:
            ca = b.gp(1.0 / lb);
            cb = a.gp(1.0 / la);

            int tempg = ga;
            ga = gb;
            gb = tempg;
        }

        // compute delta product & 'normalize'
        Multivector d, _d = deltaProduct(ca, cb);
        int gd = _d.grade();
        double ld = _d.largestCoordinate();
        d = _d.gp(1.0 / ld);

    //  System.out.println("Delta = " + d + ",");

        // if delta product is scalar, we're done:
        if (gd == 0) {
            // meet = 1
            m_meet = ca;
            // join = computed from meet
            m_join = ca.op(ca.versorInverse().ip(cb, LEFT_CONTRACTION));

            return;
        }

        // if grade of delta product is equal to ga + gb, we're done, too
        if (gd == ga + gb) {
            // a and b entirely disjoint
            // meet = 1
            m_meet = new Multivector(1.0);
            // join = computed from meet
            m_join = ca.op(cb);

            return;
        }


        // dimension of space we are working in:
        int dim = Math.max(ca.spaceDim(), cb.spaceDim());
        Multivector I = new Multivector(new BasisBlade((1 << dim)-1, 1.0));

        // init join to pseudoscalar
        Multivector j = I;
        int Ej = dim - ((ga + gb + gd) >> 1); // compute excessity of join

        // check join excessity
        if (Ej == 0) {
            // join computed
            m_join = j;
            // meet = computed from join
            m_meet = cb.ip(j.versorInverse(), LEFT_CONTRACTION).ip(ca, LEFT_CONTRACTION);
        }

        // init meet
        Multivector m = new Multivector(1.0);
        int Em = ((ga + gb - gd) >> 1); // compute excessity of meet

        // init s, the dual of the delta product:
        Multivector s = d.ip(I.versorInverse(), LEFT_CONTRACTION);

        // precompute inverse of ca
        Multivector cai = ca.versorInverse();


        // todo: maybe we can improve: search only the largest basis blade of the not-delta product?
        for (int i = 0; i < dim; i++) {
            // compute next factor 'c'
            Multivector c;

            {
            // project 'tmpc' onto 's' (the dual of the delta product)
            // project using MHIP because 's' may just be a scalar some times?
            Multivector tmpc = new Multivector(new BasisBlade(1 << i, 1.0)).ip(s, MHIP);
            c = tmpc.ip(s, MHIP); // no need to invert 's' here
            }

            // todo: then this naughty step could be avoided:

            // check if 'c' is an OK candidate:
            if (c.largestCoordinate() < largeEpsilon)
            continue;

            // compute projection, rejection of 'c' wrt to 'ca'
            Multivector cp, cr; // c projected, c rejected
            Multivector tmpc = c.ip(ca, LEFT_CONTRACTION);
            cp = tmpc.ip(cai, LEFT_CONTRACTION); // use correct inverse because otherwise cr != c - cp
            cr = c.subtract(cp);

            // if 'c' has enough of it in 'ca', then add to meet
            if (cp.largestCoordinate() > largeEpsilon) {
            m = m.op(cp);
            Em--; // reduce excessity of meet
            if (Em == 0) { // has the meet been computed???
                m_meet = m;
                // join = computed from meet
                m_join = ca.op(ca.versorInverse().ip(cb, LEFT_CONTRACTION));
                return;
            }
            }

            if (cr.largestCoordinate() > largeEpsilon) {
            j = cr.ip(j, LEFT_CONTRACTION);
            Ej--; // reduce excessity of join
            if (Ej == 0) { // has the join been computed???
                m_join = j;
                // meet = computed from join
                m_meet = cb.ip(j.versorInverse(), LEFT_CONTRACTION).ip(ca, LEFT_CONTRACTION);
                return;
            }
            }
        }
        throw new java.lang.ArithmeticException("meet & join algorithm failed!");
    }



}