/home/mark/model/software/ScrumPy/ScrumPy/Structural/Decompose.py

import sys

import ScrumPy
from ScrumPy.Util import Sci
from ScrumPy.Data import DataSets


def Consistify(E, Obs):
    "delete and reorder columns in Modes and Obs so as to leave them consistent"

    for c in E.cnames[:]:
        if  not c  in Obs.cnames:
            E.DelCol(c)

    for c in Obs.cnames[:]:
        if not c in E.cnames:
            Obs.DelCol(c)

    Obs.ColReorder(E.cnames)  # column ordering of E and Obs now identical



def DecomposeVec(E, v,Irrevs):

    sciv = Sci.vector(v)
    E = E.Copy(float)
    sciE = Sci.matrix(E.rows)
    g = Sci.pinv(sciE)

    e = v*g # <<<------  the decomposition : e = v.E#
    evals = e[0].tolist()[0]

    Again = False
    Erownames = E.rnames[:]
    for i in range(len(evals)):
        if evals[i] <0 and Erownames[i] in Irrevs:
            Again = True
            E.DelRow(Erownames[i])
    if Again:
        print "Again"
        return  DecomposeVec(E,v, Irrevs)
    else:
        vrec = e*sciE
        RelErr = Sci.Mod(vrec - v)/Sci.Mod(v)
        return E.rnames, evals, RelErr, vrec



def DecomposeSet(Modes, Obs):

    E = Modes.mo.Copy(float)
    Obs = Obs.Copy()
    

    Consistify(E,Obs)

    rv = []
    Irrevs = Modes.Irrevs()

    for v in Obs:
        rv.append(DecomposeVec(E,v,Irrevs))
        print ".",

    return rv


class Decompose(DataSets.DataSet):

    def __init__(self,Modes, Obs):

        
        self.Recs = []                                # a list of recovered fluxes
        self.Modes = Modes
        ItemNames=["RelErr"] + Modes.mo.rnames
        ScrumPy.Data.DataSets.DataSet.__init__(self, ItemNames=ItemNames)  # create self as  an empty data set for the assignments
        
        Raw_e =  DecomposeSet(Modes,Obs)            # get a list of result tuples

        for e in Raw_e:
            names,vals,relerr = e[0],e[1],e[2]
            self.NewRow()
            self[-1,0] = relerr
            for n in range(len(names)):
                self[-1,names[n]] = vals[n]
      
        #self.Recs.append(e[3])                     # hang on to the recovered fluxes - ? Why we need ??
            
        self.rnames = Obs.rnames[:]
        for c in Obs.cnames:                           # populate self with unused obs (e.g. time)
            if c not in Modes.mo.cnames:               # not an observation used in the calculation
                self.NewCol(Obs.GetCol(c),c)           # so copy it in




    def FromDS(self,ds,Copy=1):

        if Copy:
            ds = ds.Copy()
        self.rows = ds.rows
        self.rnames = ds.rnames
        self.cnames = ds.cnames
        self.Conv = ds.Conv


    def Filter(self, *args, **kwargs):

        rv = Decompose(self.Modes)
        rv.FromDS(
            ScrumPy.Data.DataSets.DataSet.Filter(self, *args, **kwargs), Copy=0
        )
        return rv