00001
00002
00003 """
00004
00005 ScrumPy -- Metabolic Modelling with Python
00006
00007 Copyright Mark Poolman 2006 -
00008
00009 This file is part of ScrumPy.
00010
00011 ScrumPy is free software; you can redistribute it and/or modify
00012 it under the terms of the GNU General Public License as published by
00013 the Free Software Foundation; either version 2 of the License, or
00014 (at your option) any later version.
00015
00016 ScrumPy is distributed in the hope that it will be useful,
00017 but WITHOUT ANY WARRANTY; without even the implied warranty of
00018 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00019 GNU General Public License for more details.
00020
00021 You should have received a copy of the GNU General Public License
00022 along with ScrumPy; if not, write to the Free Software
00023 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00024
00025 """
00026
00027 """
00028 Common entry point for "scientific" functions, especially on vectors and matrices.
00029 Currently using numpy, but this may change in future.
00030 """
00031
00032
00033 import math,sys
00034
00035 import numpy
00036
00037 from ScrumPy.Util import Seq
00038
00039
00040
00041
00042 matrix = numpy.matrix
00043 vector = numpy.array
00044 svd = numpy.linalg.svd
00045 pinv = numpy.linalg.pinv
00046
00047
00048
00049
00050
00051
00052
00053 TOL=1e-12
00054
00055 def PolyFit(Xdat, Ydat, deg=1):
00056 "polynomial fit of X,Y of degree deg"
00057
00058 return numpy.polyfit(Xdat,Ydat,deg)
00059
00060
00061
00062 def null(A,tol=TOL,AsList=False):
00063 "orthonormal nullspace, elements are floats"
00064
00065 u,w,vt=svd(A,full_matrices=1)
00066 v,i = Seq.PosOfFirst(w.tolist(),lambda x: x <= tol)
00067 if i == -1:
00068 i = len(w)
00069 rv = vt[i:,].transpose()
00070
00071 if AsList:
00072 return rv.tolist()
00073 return rv
00074
00075 def mul(a,b,AsList=False):
00076
00077 return a*b
00078
00079
00080 def Mod(vec):
00081 return math.sqrt(numpy.vdot(vec,vec))
00082
00083
00084 def CosTheta(veca,vecb):
00085 rv = numpy.vdot(veca,vecb)/(Mod(veca)*Mod(vecb))
00086 if rv >1.0:
00087 rv = 1.0
00088 elif rv < -1.0:
00089 rv = -1.0
00090
00091 return rv
00092
00093
00094 AbsCosTheta = lambda a,b: abs(CosTheta(a,b))
00095
00096 def AreParallel(a,b, tol=TOL):
00097 return abs(1-CosTheta(a,b)) <= tol
00098
00099 def AreOrthogonal(a,b,tol = TOL):
00100 return abs(CosTheta(a,b)) <= tol
00101
00102
00103 Theta = lambda a,b : math.acos(CosTheta(a,b))
00104 AbsTheta = lambda a,b : abs(Theta(a,b))
00105 Q1Theta = lambda a,b : math.acos(AbsCosTheta(a,b))
00106
00107
00108
00109
00110
00111
