利用机器学习识别嵌段共聚物自组装的微观相结构

利用机器学习识别嵌段共聚物自组装的微观相结构

  • 对浓度做傅里叶变换是识别微观相结构的极佳方法,由于快速傅里叶变换隐含的周期性,只要在做傅里叶变换时浓度矩阵是晶胞的整数倍,在做傅里叶变换之后可以消除相位的影响,即使是对浓度做了平移变换(即随机初始化跑出来的相),也不会改变傅里叶空间的形态。这一点可以在快速傅里叶变换的公式(以二维的傅里叶变换为例)中看出来,即使做了平移变换,由于周期性,对整个晶胞中所有的点求和后傅里叶空间中每个位置的值并不会发生改变。
  • 在之前的识别中,为了简化运算而切割了矩阵的过中心的三个面做快速傅里叶变换,但现在看来是得不偿失的。三维快速傅里叶变换由于增加了一个轴,因此可以适当降低精度也可以得到较多的判断点;同时切面上做傅里叶变换由于只做了两次求和,只能消除两个方向的相位,如果相发生了平移变换,这个面的浓度和傅里叶矩阵的浓度都会发生变化,因此并不能很好地用于识别。因此仍然要靠三维快速傅里叶变换来获得相的特征。在这里傅里叶空间的大小选为128*128*128,由于傅里叶变换后主要值集中在八个角上,因此截取[0:8,0:8,0:8](最初为了加快计算选择了4*4*4的正方体,但是后来发现这样精度太低,对很多相的识别并不好)这个值比较集中的小正方体来作为相的特征来进行识别。通过验证,单纯地通过峰的位置来识别并不能行得通,因为随着嵌段比例的变化,峰的位置也会发生变化,导致判断出现问题。因此考虑把这个矩阵看作一个高维向量,用支持向量机的方法,来对这个三维矩阵进行分类和识别。
  • 支持向量机(Support Vector Machine)是通过寻找一个超平面来对样本进行分割,并努力使分开的两个类别有最大的间隔,这样才能是的分隔具有更高的可信度。支持向量机通过让离分割面最近的数据点具有最大的距离。

  • 利用支持向量机可以对已知标签的向量进行分类。利用Python的scikit-learn库可以很方便地进行训练和识别。同时为了处理偏离正常位置,即跑出来的相与初始化相不同的点,加入一个惩罚因子,它代表对离群点带来的损失的重视程度,它的值越大,对目标函数的损失越大,通过设置一个值来实现参数调优。经过调试,使用线性核函数有着最好的分类效果。经过调试,识别的精确率很高,并且在csHelix只有16个样本的情况下也能实现很好的识别,同时可以获取一个浓度是什么相的概率(根据点到超平面的距离来计算)。
  • 另外,如果进行判断的相并不在训练过的相的范围内,支持向量机也会在被超平面分割出的区域中,因此也会判断出一个相,这样就发生了错误。为了判断跑出来的相是否在我们已知的范围内,需要进行异常检测,这有四种方法来实现,一类支持向量机(OneClassSVM),孤立森林(IsolationForest), 局部异常因子(LocalOutlierFactor),EllipticEnvelop,相对而言孤立森林算法效率比较高,可以指定样本的异常率来进行计算,因此选择了孤立森林算法。经过实验,在异常很少的样本中进行计算得到的包容性过强,很多其他相的相区都被包在了超平面以内从而达不到预期的效果,因此通过人工地在样本中加入一些异常值后进行计算可以得到较好的异常检测效果。

  • 训练程序,训练后生成一个estimator.pkl和classifier.pkl两个文件,第一个是孤立森林算法计算的异常值检测数据,第二个是支持向量机计算的超平面数据。训练后会自动判断训练文件的相与识别结果是否相符,并打印出来,可以根据这个来调整训练的参数。
def geteig(filename,nx=64,ny=64,nz=64):
    '''
    Get a array with three fft arrays of each component
    '''

    # Read pha file
    phafile = open(filename, 'r')
    phalines = phafile.readlines()
    pha = [[], [], []]
    sf = []
    for item in phalines:
        if not item == '\n':
            list = item.split(' ')
            pha[0].append(float(list[0]))
            pha[1].append(float(list[1]))
            pha[2].append(float(list[2]))
    phafile.close()

    # Convert the pha to a nx*ny*nz matrix
    pha = [np.array(item).reshape(nx, ny, nz) for item in pha]

    # The size of Fourier Space divided by size of pha
    d=2

    # FFT the concentration to a 64*64*64 matrix
    sf = abs(np.fft.fftn(pha[0],[64*d,64*d,64*d]))[0:int(4*d),0:int(4*d),0:int(4*d)]
    # sf = [abs(np.fft.fftn(item,[64*d,64*d,64*d]))[0:int(4*d),0:int(4*d),0:int(4*d)] for item in pha]
    # sf[1] = [abs(np.fft.fftn(item,[1,64*d,64*d]))[:,0:int(4*d),0:int(4*d)].reshape(int(4*d),int(4*d)) for item in pha]
    # sf[2] = [abs(np.fft.fftn(item,[64*d,1,64*d]))[0:int(4*d),:,0:int(4*d)].reshape(int(4*d),int(4*d)) for item in pha]
    
    # Convert all the values to [0,1]
    # for item in sf:
            # item/=item[0][0][0]
    sf /= sf[0][0][0]
    return sf

import numpy as np
import matplotlib.pyplot as plt
import os

# Import datasets, classifiers and performance metrics
from sklearn import datasets, svm, metrics

# Import joblib to save the result
from sklearn.externals import joblib

# To apply a classifier on your data, we need to flatten the image
# turn the data in a (samples, feature) matrix:

# The images to learn
images=[]
# The labels of each image to identify which phase it belongs to
target=[]

# Read phas from file and convert them into fft matrixs
path=os.listdir('d:/nonfrustrated/nonfrustrated/633')
for item in path:
    images.append(geteig('d:/nonfrustrated/nonfrustrated/633/'+item+'/pha.dat',64,64,64))
    target.append('FCC')
    print(item)
path=os.listdir('d:/nonfrustrated/nonfrustrated/233')
for item in path:
    images.append(geteig('d:/nonfrustrated/nonfrustrated/233/'+item+'/pha.dat',64,64,64))
    target.append('BCC')
    print(item)
path=os.listdir('d:/nonfrustrated/nonfrustrated/15')
for item in path:
    images.append(geteig('d:/nonfrustrated/nonfrustrated/15/'+item+'/pha.dat',64,64,64))
    target.append('A15')
    print(item)
path=os.listdir('d:/nonfrustrated/nonfrustrated/8')
for item in path:
    images.append(geteig('d:/nonfrustrated/nonfrustrated/8/'+item+'/pha.dat',64,128,48))
    target.append('csHelix')
    print(item)
path=os.listdir('d:/nonfrustrated/nonfrustrated/82')
for item in path:
    images.append(geteig('d:/nonfrustrated/nonfrustrated/82/'+item+'/pha.dat',64,128,48))
    target.append('csHelix2')
    print(item)
path=os.listdir('d:/nonfrustrated/nonfrustrated/0')
for item in path:
    images.append(geteig('d:/nonfrustrated/nonfrustrated/0/'+item+'/pha.dat',64,64,64))
    target.append('Random')
    print(item)

# Number of the samples
n_samples = len(images)

# Reshape the samples into a vector by connecting the lines head by foot
data = np.array(images).reshape((n_samples, -1))

# Create a classifier: a support vector classifier, kernel refers to the kernel function, C is a penalty parameter of the error term, tol is the tolerance for stopping criterion
classifier = svm.SVC(kernel='linear',probability=True)
# estimator=svm.OneClassSVM(kernel='linear',nu=0.001)
from sklearn.ensemble import IsolationForest
estimator=IsolationForest(max_samples=n_samples,contamination=0.04,n_estimators=n_samples)
# from sklearn.neighbors import LocalOutlierFactor
# estimator=LocalOutlierFactor(n_neighbors=15,contamination=0.1)
# from sklearn.covariance import EllipticEnvelope
# estimator=EllipticEnvelope(contamination=0.1)

# Train
classifier.fit(data[:], target[:])
estimator.fit(data[:])

# Save the result to a file
with open('classifier.pkl','wb') as file:
    joblib.dump(classifier, file)

with open('estimator.pkl','wb') as file:
    joblib.dump(estimator, file)

# Now predict the all the values:
expected = target[:]
predicted = estimator.predict(data[:])
for i in range(len(target)):
    print(expected[i],predicted[i])

predicted = classifier.predict(data[:])
for i in range(len(target)):
    print(expected[i],predicted[i])

# Print a classification report
print("Classification report for classifier %s:\n%s\n"
      % (classifier, metrics.classification_report(expected, predicted)))
print("Confusion matrix:\n%s" % metrics.confusion_matrix(expected, predicted))
  • 测试程序,可以判断是否是训练过的相,如果是,输出相的种类和概率
def geteig(filename,nx=64,ny=64,nz=64):
    '''
    Get a array with three fft arrays of each component
    '''

    # Read pha file
    phafile = open(filename, 'r')
    phalines = phafile.readlines()
    pha = [[], [], []]
    sf = []
    for item in phalines:
        if not item == '\n':
            list = item.split(' ')
            pha[0].append(float(list[0]))
            pha[1].append(float(list[1]))
            pha[2].append(float(list[2]))
    phafile.close()

    # Convert the pha to a nx*ny*nz matrix
    pha = [np.array(item).reshape(nx, ny, nz) for item in pha]

    # The size of Fourier Space divided by size of pha
    d=2

    # FFT the concentration to a 64*64*64 matrix
    sf = abs(np.fft.fftn(pha[0],[64*d,64*d,64*d]))[0:int(4*d),0:int(4*d),0:int(4*d)]
    # sf = [abs(np.fft.fftn(item,[64*d,64*d,64*d]))[0:int(4*d),0:int(4*d),0:int(4*d)] for item in pha]
    # sf[1] = [abs(np.fft.fftn(item,[1,64*d,64*d]))[:,0:int(4*d),0:int(4*d)].reshape(int(4*d),int(4*d)) for item in pha]
    # sf[2] = [abs(np.fft.fftn(item,[64*d,1,64*d]))[0:int(4*d),:,0:int(4*d)].reshape(int(4*d),int(4*d)) for item in pha]
    
    # Convert all the values to [0,1]
    # for item in sf:
            # item/=item[0][0][0]
    sf /= sf[0][0][0]
    return sf

from sklearn.externals import joblib
import numpy as np

# Load the saved learning result
with open('classifier.pkl','rb') as file:
    classifier=joblib.load(file)

with open('estimator.pkl','rb') as file:
    estimator=joblib.load(file)

# Load a file and convert to fft matrix
matrix = np.array(geteig('./pha.dat',64,64,64)).reshape((1,-1))

# Print the result
# print(estimator.predict(matrix))
if estimator.predict(matrix)[0] > 0:
    print(classifier.predict(matrix)[0])
    print(str(round(max(classifier.predict_proba(matrix)[0]),3)*100)+'%')
else:
    print('not a known phase')

发表评论