Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Created a script for SIMPLISMA implementation in python #10

Open
wants to merge 28 commits into
base: 0.4.X
Choose a base branch
from
Open
Changes from 1 commit
Commits
Show all changes
28 commits
Select commit Hold shift + click to select a range
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Add files via upload
ClarkAH authored Apr 9, 2019
commit ec43fb8c8da0def50238e31b476ecc3da9acc9be
67 changes: 67 additions & 0 deletions pymcr/tests/test_svd_simplisma.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,67 @@
import numpy as np
import random
import matplotlib.pyplot as plt
from pymcr.simplisma import svd, simplisma

#generate some sample data D
#create random normalised gaussian functions

#Number of Spectral Components
nPure = 5
#maximum number of components for SVD
nSVD = 15
#Allowed Noise Percentage
noise = 5
#manual
manual = False

x0 = np.zeros(nPure)
sigma = np.zeros(nPure)
for i in range(nPure):
x0[i] = random.uniform(-100, 100)
sigma[i] = random.uniform(3, 25)

x = np.linspace(start = -120, stop = 120, num = 2000)

gx = np.zeros((len(x),nPure))
plt.subplot(5, 1, 1)
plt.subplots_adjust(left=0.1, bottom=0.075, right=0.95, top=0.9, wspace=0.2, hspace=0.5)

for i in range(nPure):
gx[:,i] = np.exp(-(x-x0[i])**2/(2*sigma[i]**2))/np.sqrt(2*np.pi*sigma[i]**2)
plt.plot(gx[:,i])
plt.title('Real Components')

#create D with random normalised linear combination of gaussian functions
nspec = 200
D = np.zeros((len(x), nspec))
idx = list(range(nPure))

C_r = np.zeros((nspec, nPure))

for i in range(nspec):
randj = np.zeros(nPure)
random.shuffle(idx)
for j in range(nPure):
if j < nPure-1:
randj[j] = random.uniform(0, 1-np.sum(randj))
D[:,i] = D[:,i]+(gx[:,idx[j]]*randj[j])
C_r[i,j] = randj[j]
elif j == nPure-1:
randj[j] = 1-np.sum(randj)
D[:,i] = D[:,i]+(gx[:,idx[j]]*randj[j])
C_r[i,j] = randj[j]

#run SVD
eigens, explained_variance_ratio = svd.svd(D, nSVD)
nPure = np.int(input('Number of Principle Components for SIMPLISMA :'))

#Run Simplisma
S, C_u, _, _ = simplisma.pure(D, nPure, noise, False)

#Run Simplisma with constrained LCF
S, C_u, C_c, LOF = simplisma.pure(D, nPure, noise, True)