#run this and proceed to the next block

def make_orthonormal(X):
    Q,R=X.QR();
    n,p=X.dimensions();
    return Q.matrix_from_rows_and_columns(range(n),range(p));

def geodesic(X,Y):
    n,p=X.dimensions();
    XX=make_orthonormal(X);
    YY=make_orthonormal(Y);
    M=YY*(XX.T*YY).inverse()-XX*(XX.T*XX).inverse();
    svd=M.SVD();
    U,SIGMA,V=svd;
    U=U.matrix_from_rows_and_columns(range(n),range(p));
    SIGMA=SIGMA.matrix_from_rows_and_columns(range(p),range(p));
    f(x)=atan(x);
    THETA=SIGMA.apply_map(f);
    var("t");
    cc(x)=cos(t*x);
    ss(x)=sin(t*x);
    COS_THETA=THETA.apply_map(cc);
    SIN_THETA=THETA.apply_map(ss);
    for i in range(p):
        for j in range(p):
            if (i!=j):
                COS_THETA[i,j]=0;
                SIN_THETA[i,j]=0;
    subsets=(Combinations(range(n),p));
    return XX*V*COS_THETA+U*SIN_THETA;

import random
def get_random_color():
    r = lambda: random.randint(0,255)
    return ('#%02X%02X%02X' % (r(),r(),r()))

def plot_deltas(X,Y):
    n,p=X.dimensions();
    G=geodesic(X,Y);
    P=plot(0*t,(t,0,1),color="black");
    sign=1;
    subsets=(Combinations(range(n),p));
    for subset in subsets:
        if (G({t:0}).matrix_from_rows(subset).determinant()<-0.00001):
            sign=-1;
            break;
        if (G({t:0}).matrix_from_rows(subset).determinant()>0.00001):
            sign=1;
            break;
    colors={};
    for subset in subsets:
        colors[tuple(subset)]=get_random_color();
        z(t)=sign*G.matrix_from_rows(subset).determinant();
        P=P+plot(z,(t, 0, 1),color=colors[tuple(subset)]);
    P.show();
    for subset in subsets:
        html("""<font color='"""+colors[tuple(subset)]+"""'> \Delta_"""+str([i+1 for i in subset])+"""</font><br/>""")

# input two n\times p matrices X and Y (n>=p)

#X = matrix(RDF, [[1],[1]]);
#Y = matrix(RDF, [[0],[1]]);

#X=matrix(RDF,[[3,0],[3,2],[2,4]]);
#Y=matrix(RDF,[[0,1],[-1,1],[-1,0]]);

#X=matrix(RDF,[[3,0],[3,2],[2,4],[2,7]]);
#Y=matrix(RDF,[[0,1],[-1,1],[-1,0],[-1,-1]]);

#X=matrix(RDF,[[3,0],[3,2],[2,4],[2,7],[2,8]]);
#Y=matrix(RDF,[[0,1],[-1,1],[-1,0],[-1,-1],[-1,-3]]);

#X=matrix(RDF,[[3,0],[3,2],[2,4],[2,7],[2,8],[1,10]]);
#Y=matrix(RDF,[[0,1],[-1,1],[-1,-1],[-1,-1],[-1,-3],[-1,-5]]);

#X=matrix(RDF,[[3,0],[3,0],[2,4],[2,7],[4,14],[1,10]]);
#Y=matrix(RDF,[[0,1],[0,1],[-1,-1],[-1,-1],[-1,-3],[-1,-5]]);

#X = matrix(RDF, [[5],[1],[4],[2]]);
#Y = matrix(RDF, [[3],[2],[1],[0]]);

#X = matrix(RDF, [[1,0],[1,1],[0,1]]);
#Y = matrix(RDF, [[0,0],[0,1],[1,0]]);



#The cell that I'm obsessed with.
X=matrix(RDF,[[2,0,1],[1.9,0.1,1],[0,2,1],[-1,1,1],[-2,0,1],[0,0,1]]);
Y=matrix(RDF,[[1,0,0],[1,1,0],[0,1,0],[0,0,1],[0,-1,1],[1,-1,1]]);

#plot the minors along the geodesic
plot_deltas(X,Y);