Commit 024ff07b by Momme Butenschoen

### Initial commit, ported from https://gitlab.ecosystem-modelling.pml.ac.uk/momm/pml-python-tools.

parents
.gitignore 0 → 100644
 *.pyc *~ *.swp __pycache__ PCA.egg-info/ _build
This diff is collapsed.
MANIFEST.in 0 → 100644
 # Include the license file include LICENSE.txt include README.rst # Include the data files # recursive-include data * # If using Python 2.6 or less, then have to include package data, even though # it's already declared in setup.py
PCA/PCA.py 0 → 100644
 from __future__ import print_function from scipy.linalg import svd from numpy.ma import getmaskarray,masked_where,compress_rows,dot,diag from numpy.ma import zeros as mzeros from numpy.ma import array as marray from numpy import any,nonzero,array,rank,sqrt,eye,ones,arange from Map2D import Map2D def standard(d): """ Standardise 2D data along the last dimension. Args: d(numpy.array float): input data """ d=d.transpose() d=(d-d.mean(0))/d.std(0) return d.transpose() def samplemean(d): """ Take out the mean of 2D data along the last dimension. Args: d(numpy.array float): input data """ d=d.transpose() d=d-d.mean(0) return d.transpose() def Varimax(x,norm=None,tol=1.e-10,nit=1000): """ Computex varimax rotation of x. Args: x (numpy.array float): input data norm (optional): normalising function to use on x tol (float): error tolerance for iteration exit nit (integer): maximum number of iterations Returns: numpy.arry float): rotation matrix xrot=dot(x,r) """ if norm!=None: x=norm(x) p,nc=x.shape print("\tcomputing varimax rotation for ",nc," components...") TT=eye(nc) d=-1 for i in arange(nit): z=dot(x,TT) B=dot(x.transpose(), ( z**3-dot( z,diag(dot(ones([p]),z**2).squeeze()) )/p )) u,s,v=svd(B) TT=dot(u,v.transpose()) d2=d d=diag(s).sum() if d20: r=Varimax(u[:,:rotate]) if expand and reductDim==1: print("\texpanding scores...") scores=mzeros([no,self.m2d.columns]) for n,v in enumerate(scores): v[idx1D]=sv[n,:] scores[n,:]=masked_where(Mask,v.data) else: scores=sv self.scores=self._remapScores(scores) if expand and reductDim==0: print("\texpanding loadings...") EV=mzeros([no,self.m2d.rows]) u=u[:,:no].transpose() for n,pc in enumerate(EV): pc[idx1D]=u[n,:] EV[n,:]=masked_where(Mask,pc.data) else: EV=u[:,:no].transpose() self.EV=self._remapEV(EV) self.Sing=s[:no] self.Var=s[:no]**2 #variances of components (eigenvalues) self.VarFrac=s[:no]**2/(s**2).sum() #fractions of total variance #self.rot*X rotates X on varimax dirs: if rotate: self.rot=r.transpose() def __call__(self,n): """ Computes scores for nth prinipal component. Args: n (integer): number of component to compute. Returns: tuple of float numpy.arrays: eigenvector and score of nth component. """ return self.EV[n,:],self.scores[n,:] def squaredDeviations(self): """score squared relative to total variance (contribution to variance per sample and component)""" sD=self.scores()**2/self.V.sum() return sD.reshape([self.EV.shape[0]]+list(self.sampleShape)) def _map(self,data): """ Maps data into 2D shape. Args: data (numpy.array float): input data. Returns: numpy.array float: data mapped onto 2D.shape. """ data=self.m2d(data) print("Mapped data to 2D shape:",data.shape) return data def _remap(self,data): """ Remaps data into full shape. Args: data (numpy.array float): 2D data. Returns: numpy.array float: data mapped onto full shape. """ data=self.m2d.inv(data) print("Remapped data to full shape:",data.shape) return data def _remapEV(self,EV): """ Remaps Eigenvalues into full shape. Args: EV (numpy.array float): 2D data. Returns: numpy.array float: Eigenvalues mapped onto full shape. """ EV=self.m2d.Xinv(EV) return EV def _remapScores(self,scores): """ Remaps scores into full shape. Args: scores (numpy.array float): 2D data. Returns: numpy.array float: scores mapped onto full shape. """ scores=self.m2d.Yinv(scores) return scores
PCA/__init__.py 0 → 100644
 """ PCA package initialisation. """ from .PCA import *