#3x3 w/ 3 eigenvalues example with eigenspaces highlighted
example1();

#w/ 2 eigenvalues but diagonalizable
example2()

#w/ 2 eigenvalues but not diagonalizable
example3()

#w/ 1 real and 2 complex eigenvalues
example4()

def example1() :
    print 'An example with 3 distinct eigenvalues.'
    a_1=vector([-5, -3, 3]);
    a_2=vector([-1, 2, 1]);
    a_3=vector([-7, -2, 5]);
    A = matrix([a_1, a_2, a_3]).transpose();
    print 'A: \n', A
    print 'Its eigenvalues: \n',A.eigenvalues()
    v_3=A.eigenvectors_right()[0][1][0];
    v_1=A.eigenvectors_right()[1][1][0];
    v_n2=A.eigenvectors_right()[2][1][0];
    print 'A basis of eigenvectors B=(', v_3, v_1, v_n2, ')'
    B=matrix([v_3, v_1, v_n2]).transpose();
    print 'A can be diagonalized:'
    print B.inverse()*A*B
    E_n2 = line([v_n2, -v_n2], color='magenta');
    P_n2=plot(.3*v_n2, color='purple', thickness=10)
    E_1 = line([v_1, -v_1], color='palegreen');
    P_1=plot(.3*v_1, color='forestgreen', thickness=8)
    E_3 = line([v_3, -v_3], color='pink');
    P_3=plot(.3*v_3, color='red', thickness=8)
    print 'A plot of these eigenvectors and their eigenspaces:'
    show(E_n2+E_1+E_3+P_1+P_n2+P_3)
    print 'The previous plot after multiplication by A:'
    P_n2=plot(A*.3*v_n2, color='purple', thickness=5)
    P_1=plot(A*.3*v_1, color='forestgreen', thickness=8)
    P_3=plot(A*.3*v_3, color='red', thickness=3)
    show(E_n2+E_1+E_3+P_1+P_n2+P_3)
    # e_1=vector([1,0,0]);
#     e_2=vector([0,1,0]);
#     e_3=vector([0,0,1]);
#     print 'Compare to the unnatural standard basis picture:'
#     p1=plot(.5*e_1);
#     p2=plot(.5*e_2);
#     p3=plot(.5*e_3);
#     show(E_n2+E_1+E_3+p1+p2+p3)
#     print 'Where multiplication by A is not as pretty:'
#     p1=plot(A*.5*e_1);
#     p2=plot(A*.5*e_2);
#     p3=plot(A*.5*e_3);
#     show(E_n2+E_1+E_3+p1+p2+p3)
    
    
def example2() :
    print 'A diagonalizable example with only 2 distinct eigenvalues.'
    a_1=vector([0, -1, 1]);
    a_2=vector([0.5, 1.5, -0.5]);
    a_3=vector([-1.5, -1.5, 2.5]);
    A = matrix(QQ,[a_1, a_2, a_3]).transpose();
    print 'A: \n', A
    print 'Its eigenvalues: \n',A.eigenvalues()
    v_2=A.eigenvectors_right()[0][1][0];
    v_1a=A.eigenvectors_right()[1][1][0];
    v_1b=A.eigenvectors_right()[1][1][1];
    print 'Note that 2 eigenvalues in the list are the same. (Eigenvalue 1 has algebraic multiplicity 2)'
    print 'A basis of eigenvectors B=(', v_2, v_1a, v_1b, ')'
    B=matrix([v_2, v_1a, v_1b]).transpose();
    print 'A can be diagonalized:'
    print B.inverse()*A*B
    E_2 = line([1.5*v_2, -1.5*v_2], color='magenta');
    P_2=plot(.5*v_2, color='purple', thickness=5)
    var('x y z');
    eq1=-2*x+1*y-3*z==0;
    E_1=implicit_plot3d(eq1  ,(x, -1.5, 1.5), (y, -1.5, 1.5) ,(z, -1.5, 1.5), color='palegreen', opacity=0.4);
    P_1a=plot(v_1a, color='forestgreen', thickness=5)
    P_1b=plot(v_1b, color='forestgreen', thickness=5)
    print 'A plot of the 2-eigenspace (purple) and a plane that the LT rotates by a quarter and scales by the square root of 2 (green)'
    show(E_2+E_1+P_2+P_1a+P_1b)
    print 'The previous plot after multiplication by A:'
    P_1a=plot(A*v_1a, color='forestgreen', thickness=5)
    P_1b=plot(A*v_1b, color='forestgreen', thickness=5)
    P_2=plot(A*.5*v_2, color='purple', thickness=3)
    show(E_2+E_1+P_2+P_1a+P_1b)
    # e_1=vector([1,0,0]);
#     e_2=vector([0,1,0]);
#     e_3=vector([0,0,1]);
#     print 'Compare to the unnatural standard basis picture:'
#     p1=plot(.5*e_1);
#     p2=plot(.5*e_2);
#     p3=plot(.5*e_3);
#     show(E_2+E_1+p1+p2+p3)
#     print 'Where multiplication by A is not as pretty:'
#     p1=plot(A*.5*e_1);
#     p2=plot(A*.5*e_2);
#     p3=plot(A*.5*e_3);
#     show(E_2+E_1+p1+p2+p3)

def example3() :
    print 'A nondiagonalizable example with real eigenvalues.'
    a_1=vector([0, -1, 1]);
    a_2=vector([-0.5, 1, 0]);
    a_3=vector([-2.5, -2, 3]);
    A = matrix(QQ,[a_1, a_2, a_3]).transpose();
    print 'A: \n', A
    print 'Its eigenvalues: \n',A.eigenvalues()
    v_2=A.eigenvectors_right()[0][1][0];
    v_1=A.eigenvectors_right()[1][1][0];
    v=vector([1/2, -1/2, -1/2])
    print 'Note that 2 eigenvalues in the list are the same. (Eigenvalue 1 has algebraic multiplicity 2)'
    print 'A maximal linearly independent set of eigenvectors B=(', v_2, v_1, ')'
    B=matrix([v_2, v_1, v]).transpose();
    print 'A can be only partially diagonalized:'
    print B.inverse()*A*B
    E_2 = line([1.5*v_2, -1.5*v_2], color='magenta');
    P_2=plot(.5*v_2, color='purple', thickness=5)
    E_1= line([1.5*v_1, -1.5*v_1], color='palegreen');
    P_1=plot(v_1, color='forestgreen', thickness=5)
    P_v=plot(v, thickness=5)
    print 'A plot of these eigenvectors and their eigenspaces:'
    show(E_2+E_1+P_2+P_1+P_v)
    print 'The previous plot after multiplication by A:'
    P_1=plot(A*v_1, color='forestgreen', thickness=5)
    P_v=plot(A*v, thickness=5)
    P_2=plot(A*.5*v_2, color='purple', thickness=3)
    show(E_2+E_1+P_2+P_1+P_v)
    
def example4() :
    print 'An example with 1 real and 1 complex conjugate pair of eigenvalues.'
    a_1=vector([4, 0, -2]);
    a_2=vector([-1, 1, 1]);
    a_3=vector([3, -1, -1]);
    A = matrix(QQ,[a_1, a_2, a_3]).transpose();
    print 'A: \n', A
    print 'Its eigenvalues: \n',A.eigenvalues()
    v_2=A.eigenvectors_right()[0][1][0];
    v_r1=.5*vector([-1, -1, 1]);
    v_r2=.5*vector([1, -1, -1]);
    print 'Note that 2 eigenvalues are complex conjugates.'
    print 'A maximal list of (real) eigenvectors B=(', v_2, ')'
    B=matrix(QQ, [v_2, v_r1, v_r2]).transpose();
    print 'However, the rest of A can be "rotationalized":'
    print B.inverse()*A*B
    E_2 = line([1.5*v_2, -1.5*v_2], color='magenta');
    P_2=plot(.5*v_2, color='purple', thickness=5)
    var('x y z');
    eq1=x+z==0;
    E_r=implicit_plot3d(eq1  ,(x, -1.5, 1.5), (y, -1.5, 1.5) ,(z, -1.5, 1.5), color='palegreen', opacity=0.4);
    P_r1=plot(v_r1, color='forestgreen', thickness=5)
    P_r2=plot(v_r2, color='forestgreen', thickness=5)
    print 'A plot of these eigenvectors and their eigenspaces:'
    show(E_2+E_r+P_2+P_r1+P_r2)
    print 'The previous plot after multiplication by A:'
    P_r1=plot(A*v_r1, color='forestgreen', thickness=5)
    P_r2=plot(A*v_r2, color='forestgreen', thickness=5)
    P_2=plot(A*.5*v_2, color='purple', thickness=3)
    show(E_2+E_r+P_2+P_r1+P_r2)