Commit d8ff8388 authored by Francisco Javier Dominguez Gutierrez's avatar Francisco Javier Dominguez Gutierrez
Browse files

Update pca.py

parent eadb9c4c
......@@ -31,6 +31,13 @@
##
######################################################################################################
##########
# Preamble
##########
import sys, argparse
import os, subprocess
import gc
import numpy as np
from pylab import *
from numpy import *
......@@ -38,22 +45,41 @@ import linecache
import matplotlib.pyplot as plt
from sklearn.decomposition import TruncatedSVD
os.system('clear')
print()
print(" ------------------------------------------------- ")
print("| Fingerprint and visualization of defects |")
print("| in damaged crystal structures |")
print(" ------------------------------------------------- ")
print()
# Check Python version
if sys.version_info[0] < 3:
print("\nRequires Python version 3 or greater\n")
sys.exit()
#####################################################
## The following input parameters are set for the Fe case,
## which is explaining in the paper
##########################################################
num_atom = 314928
i = 0
arr = []
dist_diff = []
# Cleaning quip output file
# Start reading DV at line 13
# Start reading atom DV at line 13
for i in range(13,13+num_atom):
dv_1 = np.array(linecache.getline('../vectors.dat',i)[7::].split())
dv_1 = np.array(linecache.getline('vectors.dat',i)[7::].split())
dv = dv_1.astype(np.float)
arr.append(dv)
## Reading the id of the unclassified atoms in the damaged Fe sample
ref_v = np.loadtxt('id_atoms.dat')
lent1 = len(ref_v)
## Extracting the DV of the unclassified atom from the set of atoms DVs.
Ar = []
for d in range(0,lent1):
i = int(ref_v[d])
......@@ -69,8 +95,8 @@ for i in range(0,50):
A_m = np.mean(Ar.T,axis=1)
Ap = Ar-A_m
# svd
## Appplying PCA by using
## Singular Value Decomposition
svd = TruncatedSVD(n_components=2)
svd.fit(Ap)
result = svd.transform(Ap)
......@@ -80,6 +106,7 @@ print("lengt of the truncated vector: {0}".format(lengt))
for i in range(0,lengt):
thefile.write("{0} {1} {2} \n".format(-result[i,0],' ', -result[i,1]))
## Plotting the results
fig, ax = plt.subplots()
rc('axes', linewidth=2)
plt.scatter(-result[:,0],-result[:,1],alpha=1.00,color="blue",marker="x",s=50,linewidth=2)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment