/home/mark/model/software/ScrumPy/ScrumPy/Util/Seq.py

"""

ScrumPy -- Metabolic Modelling with Python

Copyright Mark Poolman 1995 - 2002

 This file is part of ScrumPy.

    ScrumPy is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    ScrumPy is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with ScrumPy; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

"""
"""
M.G. Poolman 10/01/01 onwards
Seq.py
some useful things to do to sequences
"""

import math,random

import Sorter

#   First some Vector math type stuff
#   These will be (*much*) slower than the NumPy equivs, but are more flexible
#   with regard to element type


def Mod(seq):
    """ pre: all elements in seq are numeric
       post: return |seq| as a float """

    return math.sqrt(
        sum(map(lambda x:x*x, map(float,seq)))
    )



def DotProd(a,b):
    """ pre: all els in a and b are numeric, len(a)==len(b)
       post: returns the dot (scalar) product of a and b """

    return sum(map(lambda x,y:x*y, a,b))




def CosTheta(a,b):
    """ pre: Mod(a) !=0, Mod(b) != 0
       post: the cosine of the angle between a and b """

    rv = DotProd(a,b)/(Mod(a) * Mod(b))
    if rv > 1.0:
        rv = 1.0
    elif rv < -1.0:
        rv = -1.0
    return rv

def AreParallel(a,b):
    """ pre: len(a) == len(b), all items are of ArbRat type
      post: AreParallel(a,b) =>(a and b are parallel, if parallel relative magintitude) """

    lenab = len(a)
    rat1 = 0

    for i in range(lenab):
        if b[i] !=0:
            if a[i] ==0:
                return False,0
            if rat1 == 0:
                rat1 = a[i]/b[i]
            elif rat1 != a[i]/b[i]:
                    return False,0
        elif a[i] != 0:
            return False,0
    return True,rat1





def Angle(a,b):
    """ pre: Mod(a) != 0, Mod(b) != 0
      post: Angle in radians between a and b"""
    # This function should be in Sci ?

    ct = float(CosTheta(a,b))
    # trap rounding errors resulting in -1.0 < ct < 1.0 (and marginal time saving for ct == +/-1.0)
    if ct >=1.0:
        return 0
    if ct <= -1.0:
        return math.pi
    try:
        return math.acos(ct)
    except:
        # probably won't get here, but we may as well give a diagnostic if we do
        print a, "\n", b,"\n", ct
        raise "Can't calculate Angle"

def SumSq(seq):

    return sum(map(lambda x:x*x,seq))

    #rv = seq[0]**2
    #for i in range(1, len(seq)):
    #    rv += seq[i]**2
    #return rv


def EucDist(a,b):
    " the Euclidean distance between a,b "

    moda = Mod(a)
    modb = Mod(b)
    cosab = DotProd(a,b)/(Mod(a) * Mod(b))

    return math.sqrt(abs(moda*moda + modb*modb - 2*moda*modb*cosab))


def Sum(s):

    sys.stderr.write("Util.Sum ?? WTF - Deprecated")
    raise  "obsolete function"
    rv = 0*s[0]

    for i in s:
        rv+= i
    return rv

def Deriv(x,y):
    """ pre: len(x) == len(y)==3, all x and y numeric
       post: Deriv(x,y) == slope at x[1],y[1], by three point linear derivation """
    return 0.5 * (
        ((y[1] - y[0]) / (x[1] - x[0])) +
        ((y[2] - y[1]) / (x[2] - x[1]))
    )



def Fulcrum(Seq):

    rv = 0
    halflen = len(Seq)/2.0


    for idx in range(len(Seq)):
        rv = rv + (Seq[idx] * (idx-halflen))

    return rv


def AbsFulcrum(Seq):

    rv = 0
    halflen = len(Seq)/2.0


    for idx in range(len(Seq)):
        rv = rv + (abs(Seq[idx]) * (idx-halflen))

    return rv


def Diff(a,b):
    """ pre: len(a) == len(b)
                a[0..i] - b[0..1] valid for i in range len(a)
        post: Diff(a,b)[0..1] = a[0..i]-b[0..i] """

    return map(lambda x,y: x-y, a,b)


##
## more generalised sequence functions
##


def Round(seq,sf):
    """ pre: type(seq[0..n]) == FloatType
       post: seq[0..n] rounded to sf sig figs """

    for i in range(len(seq)):
        m,e = math.frexp(seq[i])
        seq[i] = math.ldexp(round(m,sf),e)

def Flatten(seq, eta):
    for i in range(len(seq)):
        if abs(seq[i]) < eta:
            seq[i] = 0.0



def Fill(seq, fun, args, kwargs):
    """ fill seq with fun(args,kwargs) """
    seq[:] = map(lambda x, a=args, k=kwargs, f=fun: f(*a,**k), seq)


def Split(l,sep):
    """ split l in to list of lists about sep, list equivalent of str.split()"""

    rv = []
    l = l[:]
    l.reverse()
    while sep in l:
        idx = l.index(sep)
        slice = l[:idx]
        slice.reverse()
        rv.append(slice)
        del l[:idx+1]
    if len(l) >0:
        l.reverse()
        rv.append(l)
    rv.reverse()
    return rv

def Diffs(a,b):
    """ pre: a[n]-b[n] valid for (n <= 0 < len(a)), len(a) == len(b)
       post: Diffs(a,b)[n] == a[n] - b[n]
    """
    return map(lambda x,y:x-y, a,b)

def AllEq(seq, val):
    """
    pre : a[n] != val valid for (n <= 0 < len(a))
    post: AllEq(seq, val) => all elements in Seq are equal to val
    """
    for el in seq:
        if el != val:
            return 0

    return 1


def AreEqual(s1, s2):
    lens1 = len(s1)

    if lens1 != len(s2):
        return 0

    for i in range(lens1):
        if s1[i] != s2[i]:
            return 0

    return 1


def LimInSeq(seq, fun):
    """
    pre : fun extracts a value from seq w/o changing it
    post: LimInSeq(seq, fun) == index and value in seq found by fun(seq)
    """
    val = fun(seq)
    for idx in range(len(seq)):
        if seq[idx] == val:
            return [idx, val]


def PosOfFirst(seq, fun):
    """pre: fun is a boolean function which can accept any element of seq as input
      post: returns(val,idx) where val and index are the first value and index where fun is true
               returns (None, -1) if no matches """

    for el in seq:
        if fun(el):
            return el, seq.index(el)
    return None,-1



def PosOfValsInSeq(seq, val):
    """
    pre : all elements in seq can be compared with val
    post: PosOfValsInSeq == list of indices to elements in seq of value val
    """
    rv = []
    for idx in range(len(seq)):
        if seq[idx] == val:
            rv.append(idx)

    return rv


def Abs(seq):
    return map(abs, seq)



def PosOfAbsMin(seq):
    seqabs = Abs(seq)
    seqnz = []
    for el in seqabs:
        if el != 0:
            seqnz.append(abs(el))

    return seqabs.index(min(seqnz))

def PosOfMinBeforeZero(seq):
    """pre: seq[0] != 0, seq[n]==0 in range(1, len(seq))
      post: returns index of  smalles element before the first 0 val element """

    rv = 0
    for i in range(1,len(seq)):
        if seq[i] == 0:
            return rv
        else:
            if seq[i] < seq[rv]:
                rv = i


def RemoveFrom(lis, targs):
    """pre : lis and targs are sequences
      post: no elements of targs in lis """

    for el in targs:
        try:
            del lis[lis.index(el)]
        except:
            pass


def NoDups(lis):
    " returns copy of l with duplicate items removed"

    rv = []
    for el in lis:
        if el not in rv:
            rv.append(el)

    return rv

def HasMatches(seq,fun):
    for el in seq:
        if fun(el):
            return 1
    return 0


def IdxListOfMatches(seq, fun):
    """
    pre: all elements in seq can be passed to fun which returns a boolean
    post: returns a list of indices for which fun(seq[idx])==1 (true)
    """
    rv = []
    for idx in range(len(seq)):
        if fun(seq[idx]):
            rv.append(idx)

    return rv


def NearestIdxVal(l,v):
    """ pre: l is ascending ordered numeric list, v is a value  comparable to elements of l
       post: returns (i,n) i the index of n in l which is closest in value to v """

    def nvi(l,v,idx):
        lenl = len(l)
        if lenl<=3:
            diffs = map(lambda x: abs(x-v), l)
            idx = diffs.index(min(diffs))
            return idx,l[idx]

        mid = lenl/2

        if v >= l[mid]:
            idx,nearest=nvi(l[mid:],v,idx)
            idx += mid
            return idx,nearest
        return nvi(l[:mid+1],v,idx)

    return nvi(l,v,-1)


def IdxOfUnique(seq, fun):
    """
    pre: all elements in seq can be passed to fun which returns a boolean
    post: returns idx for which only fun(idx) is true, else returns None
    """

    match = IdxListOfMatches(seq, fun)
    if len(match) == 1:
        return match[0]
    return None

def NZIdxs(seq,Abs=True):
    "list of indices to nonzero elements of seq, sorted by ascending (Abs => absolute) element value "

    idxs = []
    for i in range(len(seq)):
        if seq[i] !=0:
            idxs.append(i)

    s = Sorter.Sorter(idxs)
    if Abs:
        def fun(a):
            return abs(seq[a])
    else:
        def fun(a):
            return seq[a]

    return s.Sort(fun)

def NumOfMatches(seq, fun):
     """
    pre: all elements in seq can be passed to fun which returns a boolean
    post: returns the number of elements for which fun(seq[idx])==1 (true)
    """

     return len(IdxListOfMatches(seq, fun))


def AllMatch(seq, fun):
    """
    pre: all elements in seq can be passed to fun which returns a boolean
    post: AllMatch(seq, fun) => fun(el) == True for all el in seq
    """

    return NumOfMatches(seq, fun) == len(seq)




def SwapItem(seq, a,b):
    tmp = seq[a]
    seq[a] = seq[b]
    seq[b] = tmp


def Normalise(seq, nidx=-1, KeepSign=1):

    abseq = Abs(seq)

    if nidx < 0:
        nidx = abseq.index(max(abseq))

    val = seq[nidx] * 1# * 1 : force creation of new instance for val
    if KeepSign:
        val = abs(val)

    if val == 0:
        raise "Division by zero in Normalise"

    if val != 1:
        for i in range(len(seq)):
            seq[i] /= val


def NormaliseRound(seq, nidx=-1, dec= 12, KeepSign=1):

    if nidx < 0:
        nidx = seq.index(max(seq))

    if seq[nidx] == 0:
        raise "Division by zero in Normalise"
    else: val = 1.0 / seq[nidx] # force creation of new instance for val
    if KeepSign:
        val = abs(val)

    i = 0
    while i < len(seq) :
        tmp = seq[i]
        if tmp != 0.0 : seq[i] = round (tmp * val, dec)
        else : seq[i] = 0.0
        i = i+1


def Add(a,b):

    rv = a[:]

    for i in range(len(b)):
        rv[i] += b[i]
    return rv



##
##
##
##sequences of rationals
##


def LargestDen(seq):
    """ pre: all items in Seq have GetDen() method, returning comparable +ve values
       post: LargestDen(seq) = largest denominator in seq """

    rv = seq[0].denom()
    for i in seq[1:]:
        d = i.denom()
        if d > rv:
            rv = d

    return rv *1


def RandOrder(seq):
    """ pre: true
       post: return a list with items in random order wrt seq"""


    return random.sample(seq,len(seq))


def Integise(seq):
    """ pre: LargestDen(seq) is defined
       post: seq' = Integise(seq)
             seq' items are integer, maintaining ratios of seq """

    den = LargestDen(seq)

    while den > 1:
        seq = map(lambda x: x* den,seq)
        den = LargestDen(seq)

    #for id in range(len(seq)):
     #   seq[id] = int(seq[id])

    return seq