Surface NMR processing and inversion GUI
您最多选择25个主题 主题必须以字母或数字开头,可以包含连字符 (-),并且长度不得超过35个字符

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768
  1. ##########################################################
  2. #
  3. #
  4. #
  5. #
  6. #
  7. #from scipy import mean
  8. import numpy
  9. def pca(A, remove=[]):
  10. ''' The input matrix A should be a 2D numpy array in column-major
  11. order. Each column is a dataset and PCA will be applied across
  12. columns
  13. '''
  14. #nrow = len(A[:,0])
  15. #ncol = len(A[0,:])
  16. nrow,ncol = numpy.shape(A)
  17. # Cast into matrix type for easier math
  18. A = numpy.matrix(A)
  19. # Allocate Covariance Matrix
  20. covMatrix = numpy.matrix(numpy.zeros((nrow, nrow)))
  21. # Compute Means of each row. Each column must be normalised
  22. meanArray = []
  23. for i in range(nrow):
  24. meanArray.append(numpy.mean(A[i].tolist()[0]) )
  25. A[i] -= meanArray[i]
  26. meanArray = numpy.array(meanArray)
  27. # Generate Covariance Matrix
  28. covMatrix = numpy.cov(A)
  29. # Compute Eigen Values, Eigen Vectors
  30. eigs = numpy.linalg.eig(covMatrix)
  31. K = eigs[1].T
  32. #print K
  33. #print "Eigen Values"
  34. #print eigs[0]
  35. # Zero requested components
  36. for i in remove:
  37. K[i] = numpy.zeros(len(K))
  38. # Make Transform Matrices
  39. transMatrix = K*A
  40. # Return Necssary Stuff
  41. return numpy.array(transMatrix), K, meanArray
  42. #return K, A, meanArray
  43. #def invpca(K, A, means):
  44. def invpca(transMatrix, K, means):
  45. '''Converts a PCA rotated dataset back to normal. Input parameters
  46. are the components to discard in re-creation.
  47. '''
  48. K = numpy.matrix(K)
  49. # Transform and Untransform the data
  50. #transMatrix = K*A
  51. untransMatrix = K.T*transMatrix
  52. # Correct for normalisation
  53. for i in range(len(transMatrix)):
  54. untransMatrix[i] += means[i]
  55. return untransMatrix