pca-code/pcatest.py at main · xstarkit/pca-code · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
#!/bin/python

# Author: Michael L. Parker
# Last Modification: 23 Feb 2015

import pylab as p
import numpy as n
from numpy.linalg import svd
from scipy.signal import sawtooth

xs=n.linspace(0,6.*n.pi,200)

def line(xs):
    '''This function just defines a straight line'''
    ys=[]
    for x in xs:
        ys.append(x/8+0.3)
    return ys

# Specify number of 'spectra' to generate
nspec=15

speclist=[]

# This is all plotting, ignore it.
figure1=p.figure(figsize=(20,6))
ax=p.subplot2grid(shape=(3,3),loc=(0,0))
p.plot(line(xs))
ax=p.subplot2grid(shape=(3,3),loc=(1,0))
p.plot(n.sin(xs))
ax=p.subplot2grid(shape=(3,3),loc=(2,0))
p.plot(sawtooth([x*2 for x in xs]))
ax=p.subplot2grid(shape=(3,3),loc=(0,1),rowspan=3)

for i in range(0,nspec):
    # For each spectrum, generate three random numbers.
    # Distributions are flat, centred on zero, but have differing amplitudes.
    randval1=0.2*(n.random.rand()-0.5)
    randval2=0.6*(n.random.rand()-0.5)
    randval3=0.1*(n.random.rand()-0.5)

    # Multiply three different input components (sin wave, sawtooth and stright line) by random values
    sinespec=[randval1 * s for s in n.sin(xs)]
    sinespec2=[randval3 * s for s in sawtooth([x*2 for x in xs])]
    linespec=[randval2 * l for l in line(xs)]

    # Generate some noise
    noisespec=[(n.random.rand()-0.5)*0.1 for x in xs]

    # Add all components together, with noise, and append to spectrum list
    speclist.append([a+b+c+d for a,b,c,d in zip(sinespec,linespec,noisespec,sinespec2)])
    p.plot(speclist[-1])

specarray=n.array(speclist)
specarray=n.transpose(specarray)

specfile='specfile.dat'
n.savetxt(specfile,specarray)

# All the analysis is in this line. This decomposes the spectrum array into constituent components.
a,b,c=svd(specarray)

# Matrix a now contains all the output components, vector b their strengths. Ignore c.
a=n.transpose(a)

# Plot the components
ax=p.subplot2grid(shape=(3,3),loc=(0,2))
p.plot(a[0])
ax=p.subplot2grid(shape=(3,3),loc=(1,2))
p.plot(a[1])
ax=p.subplot2grid(shape=(3,3),loc=(2,2))
p.plot(a[2])
p.show()

#print 'OUTPUT FUNCTIONS'
#for a0,a1,a2 in zip(a[0],a[1],a[2]):
	#print a0,a1,a2


total=0
for val in b:
    total+=val**2

vals=[val**2/total for val in b]
#print '\nEIGENVALUES'
#for val in vals:
	#print val

# p.plot(range(1,len(vals)+1),n.log10(vals),ls='None',marker='x')
# p.show()