var('x1,x2,y1,y2')
X = matrix(SR,[[x1,x2]]).transpose()
A = matrix(QQ,[[1,3],[3,4]])
F = X.transpose()*A*X
F = (F[0][0]).expand()
A,X,F

G1=implicit_plot(F,(x1,-2,2),(x2,-2,2),contours=[j/2 for j in range(20)],color='black')
G1

evr = A.eigenvectors_right()
evr

f1 = vector(evr[0][1][0])
f2 = vector(evr[1][1][0])
f1,f2

f1[1].minpoly()

f2[1].minpoly()

origo = vector([0,0])
er = arrow(origo,f1,color='red') + arrow(origo,f2,color='red')
G1 += er
G1 += circle(origo,1,color='blue')
G1

T = matrix([f1,f2]).transpose()
T

B = T.inverse()*A*T
B

Y = matrix(SR,[[y1,y2]]).transpose()

G = Y.transpose()*B*Y
G = (G[0][0]).expand()
print(B)
print(Y)
print(G)

G2=implicit_plot(G,(y1,-2,2),(y2,-2,2),contours=[j/2 for j in range(20)],color='black')
origo = vector([0,0])
er2 = arrow(origo,vector([1,0]),color='green') + arrow(origo,vector([0,1]),color='green')
G2 += er2
G2 += circle(origo,1,color='blue')
G2

var('a,b')
M=matrix(SR,[[a,1],[1,b]])
M

(X.transpose()*M.subs([a==1,b==2])*X)[0][0]

implicit_plot((X.transpose()*M.subs([a==1,b==2])*X)[0][0], 
              (x1,-2,2),(x2,-2,2),
              contours=[j for j in range(10)],axes=True)

kvf = [implicit_plot((X.transpose()*M.subs([a==1,b==u])*X)[0][0], 
              (x1,-2,2),(x2,-2,2),
              contours=[j for j in range(10)],axes=True) for u in sxrange(-3,3,0.2)]

animate(kvf)

var('x,y,z')

t = z.function(z,y,z)
cm=colormaps.PiYG

implicit_plot3d(x^2+y^2+2*z^2,(x,-4,4),(y,-4,4),(z,-4,4),contour=[16],
                color='seagreen', frame=True, mesh=True, transparent=True, opacity=0.2)

implicit_plot3d(x^2+y^2-2*z^2,(x,-4,4),(y,-4,4),(z,-4,4),contour=[4],
                color='seagreen', frame=True, mesh=True,opacity=0.2)

implicit_plot3d(x^2-y^2-2*z^2,(x,-4,4),(y,-4,4),(z,-4,4),contour=[4],
                color='seagreen', frame=True, mesh=True,opacity=0.2)

implicit_plot3d(x^2+y^2+0*z^2,(x,-4,4),(y,-4,4),(z,-4,4),contour=[4],
                color='seagreen', frame=True, mesh=True,opacity=0.2)

implicit_plot3d(x^2+y^2-1*z,(x,-4,4),(y,-4,4),(z,-4,4),contour=[4],
                color='seagreen', frame=True, mesh=True,opacity=0.2)

implicit_plot3d(x^2+y^2-1*z^2,(x,-4,4),(y,-4,4),(z,-4,4),contour=[0],
                color='seagreen', frame=True, mesh=True,opacity=0.2)

from numpy import mean,cov,double,cumsum,dot,linalg,array
from pylab import plot,subplot,axis,stem,show,figure

def princomp(A):
 """ performs principal components analysis 
     (PCA) on the n-by-p data matrix A
     Rows of A correspond to observations, columns to variables. 

 Returns :  
  coeff :
    is a p-by-p matrix, each column containing coefficients 
    for one principal component.
  score : 
    the principal component scores; that is, the representation 
    of A in the principal component space. Rows of SCORE 
    correspond to observations, columns to components.
  latent : 
    a vector containing the eigenvalues 
    of the covariance matrix of A.
 """
 # computing eigenvalues and eigenvectors of covariance matrix
 M = (A-mean(A.T,axis=1)).T # subtract the mean (along columns)
 [latent,coeff] = linalg.eig(cov(M)) # attention:not always sorted
 score = dot(coeff.T,M) # projection of the data in the new space
 return coeff,score,latent

def doit(A):
    coeff, score, latent = princomp(A.T)

    figure()
    subplot(121)
# every eigenvector describe the direction
# of a principal component.
    m = mean(A,axis=1)
    plot([0, -coeff[0,0]*2]+m[0], [0, -coeff[0,1]*2]+m[1],'--k')
    plot([0, coeff[1,0]*2]+m[0], [0, coeff[1,1]*2]+m[1],'--k')
    plot(A[0,:],A[1,:],'ob') # the data
    axis('equal')
    subplot(122)
    # new data
    plot(score[0,:],score[1,:],'*g')
    axis('equal')
    show()

C = array([ [2.4,0.7,2.9,2.2,3.0,2.7,1.6,1.1,1.6,0.9],
            [2.5,0.5,2.2,1.9,3.1,2.3,2,1,1.5,1.1] ])
doit(C)

# Vi rullar vår egen...

N = 40
xvar = [RR.random_element(-5,5) for _ in range(N)]
yvar = [xvar[j] + RR.random_element(-2,2) for j in range(N)]
B = array([xvar,yvar])

doit(B)

N=100
xvar = [RR.random_element(-5,5) for _ in range(N)]
yvar = [xvar[j]/2 + RR.random_element(-2,2) for j in range(N)]
B = array([xvar,yvar])
doit(B)