/home/mark/model/software/ScrumPy/ScrumPy/Data/Histo.py

import types,math

from Util import Seq
import Stats,  Plotter

class Histo(dict):

    def __init__(self,vals, nbins=10,minv=None,maxv=None,emptybins=True,name="Histogram"):

        self.vals  = vals[:]
        self.minv  = minv
        self.maxv  = maxv
        self.nbins = nbins
        self.emptybins = emptybins
        self.name = name
        self.style = "boxes"

        self.Plotter = Plotter.GPPlotter()
        self.Plotter.Options['yrange'] = 'set yra [0:]'
        if len(vals)>1:
            self.Update() # might be starting empty


    def __str__(self):
        rv = ""
        ks = self.keys()
        ks.sort()
        for k in ks:
            rv += " ".join([str(k), str(self[k])])+"\n"
        return rv

    def __iadd__(self,items):
        self.vals += items
        return self

    def append(self,items):
        self.vals.append(items)

    def extend(self,items):
        self.vals.extend(items)

    def Add(self, vals):
        self.vals.extend(vals)
        if min(vals) < self.minv:
            self.minv=min(vals)
        if max(vals)> self.maxv:
            self.maxv=max(vals)
        self.Update()

    def AsTuples(self):

        ks = self.keys()
        ks.sort()
        rv = []
        for k in ks:
            rv.append( (k,self[k]))
        return rv

    def AsLists(self):
        keys = self.keys()
        keys.sort()
        vals = map(lambda k: self[k], keys)
        return keys, vals

    def AsList(self):
        return self.AsLists()[1]


    def Plot(self):

        self.Plotter.AddData(self.name, self.AsLists())
        minv, maxv = map(str, (self.minv, self.maxv))
        self.Plotter.Options['xrange'] = 'set xra ['+minv +':'+maxv+']'
        self.Plotter[self.name].Options["Style"] = self.style
        self.Plotter.Plot()


    def Update(self):

        self.vals.sort()

        if self.minv == None:
            self.minv = self.vals[0]

        if self.maxv == None:
            self.maxv = self.vals[-1]

        minv = self.minv = float(self.minv)
        maxv = self.maxv = float(self.maxv)

        binsize  = (maxv - minv)/(self.nbins)
        binsize2 = binsize/2
        curbin = minv+ binsize2

        for n in self.keys():
            del self[n]

        for n in range(self.nbins):
            curbin += binsize
            self[curbin] = 0

        bins = sorted(self.keys())
        curbin = bins.pop(0)
        maxbin = curbin+binsize2

        for v in self.vals:
            while v>=maxbin:
                curbin=bins.pop(0)
                maxbin+=binsize

            self[curbin]+=1


        if self.Plotter.has_key(self.name):
            self.Plotter.RemoveData(self.name)






    def AutoParam(self):

        self.nbins = 3*int(math.ceil(len(self.vals)**0.33))
        self.minv = min(self.vals)
        self.maxv = max(self.vals)
        self.Clear()
        self.Update()


    def Clear(self):

        for k in self.keys()[:]:
            del self[k]

        self.minv=None
        self.maxv=None


    def Save(self, file):

        if type(file)==types.StringType:
            file = open(file,"w")

        file.write(str(self))


    def Normalise(self):

        ks = self.keys()
        range = max(ks) - min(ks)
        Area=sum(self.values())  *range/len(self)
        for k in ks:
            self[k] /=Area




    def ChiSquare(self, sample=None, ref=None,nbins=None,Plot=True,*args,**kwargs):
        """ pre sample = list or None, ref=list or function or None, sample=None or ref=None,nbins>2 """

        if nbins==None:
            nbins=self.nbins
        else:
            self.nbins = nbins

        self.Clear()

        if sample != None: # we are the reference
            other = Histo(sample,nbins)
            minv = min((min(sample), min(self.vals)))
            maxv = max((max(sample), max(self.vals)))
            self.minv = other.minv = minv
            self.maxv = other.maxv = maxv
            self.Update()
            if 0 in self.values():
                print "0 count(s) in reference, need more reference data or fewer bins"
                return

            rv = Stats.ChiSquare(other.AsList(),self.AsList(),True)
        else: # we are the sample

            if type(ref) == types.FunctionType:  # reference is a funtion
                other = Histo([])
                self.Update()
                for k in self.keys():
                    other[k] = ref(k,*args,**kwargs)
                rv = Stats.ChiSquare(self.AsList(),other.AsList(),True)

            else:
                self.Update()
                other = Histo(ref)  # reference is data
                rv =  other.ChiSquare(self.vals[:],nbins=nbins,Plot=False)

        if Plot:
            self.Normalise()
            other.Normalise()
            self.Plotter.AddToPlot(other.AsLists())
            self.Plot()

        return rv



def ChiSquare(obs, ref, nbins,Norm=True):

    maxv = max(obs+ref)
    minv = min(obs +  ref)
    hobs = Histo(obs,nbins,minv=minv,maxv=maxv)
    href = Histo(ref,nbins,minv=minv,maxv=maxv)
    chisq, p = Stats.ChiSquare(hobs.AsList(),href.AsList(),Norm)
    rv = {}
    rv["ChiSq"]=chisq
    rv["p_ChisSq"] = p
    rv["ObsHist"] = hobs
    rv["RefHist"] = href
    return rv




def BoolList2His(List, binsize):

    nbins,rem = divmod(len(List),binsize)

    midbin = binsize/2.
    rv = {}
    for n in range(nbins):
        lo = n*binsize
        hi = lo+binsize
        count = List[lo:hi].count(True)
        rv[lo+midbin]=count

    if rem >0:
        rv[hi+midbin] = List[hi:].count(True)
    return rv