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

"""
M.G.Poolman 20/06/01 onwards
module Simplex 
do simplex optimisations a la Nelder and Mead
loosely based on the description in Numerical Recipies in C
"""
import random

class Vertex:
    def __init__(self, Params, func, eval=1,**fargs):
        """ pre: func is a function that returns a numerical value
                 Params is a list of numbers that can be used as parameters to func """
                 
        self.Params = Params[:]
        self.nParams = len(Params)
  
        self.func = func
        if eval:
            self.Eval(**fargs)


    def __len__(self):
        return self.nParams

    def __repr__(self):
        rv = ""
        for p in self.Params:
            rv = rv + str(p) + " "
        return rv + "  " + str(self.Value)
   
    def Eval(self, **fargs):
        self.Value = self.func(self.Params, **fargs)
        
    def Copy(self):
        rv = Vertex(self.Params, self.func, eval=0)
        rv.Value = self.Value
        return rv

    def __getitem__(self,n):
       return self.Params[n]
       
    def __setitem__(self,n,value):
       self.Params[n] = value
        
    def __cmp__(self, other):
        
        diff = self.Value - other.Value
        
        if diff >0:
            return 1
        if diff < 0:
            return -1
            
        return 0

    def Randomise(self,fac):
        for pidx in range(len(self.Params)):
            p = self.Params[pidx]
            lo = p*(1-fac)
            hi = p*(1+fac)
            self.Params[pidx]  = random.uniform(lo, hi) 


class Simplex:
    
    def __init__(self, StartParas, Func, **fargs):
        
        self.Vertices = []
        self.fargs = fargs  # optional arguments passed to the user defined objective function
        for n in range(len(StartParas)):
            self.Vertices.append(Vertex(StartParas[n], Func, **fargs))
    
        self.nVertices = len(self.Vertices)    
        self.nParams = self.Vertices[0].nParams # assume all s.v[n] are same (implicit in algorithm)

        self.Vertices.sort()
        self.RecalcPsum()        
        self.WorstIdx = self.nVertices-1        # because indexing starts at 0
        self.Worst2Idx = self.WorstIdx-1
        self.BestIdx = 0

    def __repr__(self):
        rv = ""
        for v in self.Vertices:
            rv = rv + str(v) + "\n"
            
        return rv


    def Eval(self):
        for v in self.Vertices:
            v.Eval(**self.fargs)
        self.Vertices.sort()


    def StopFunc(self):
        """ users can overload this to provide custom stopping criteria """
        self.Iters = self.Iters + 1
        return self.Iters == 100


    def RecalcPsum(self):
        self.Psum = []
        for p in range (self.nParams):
            sum = 0.0 
            for v in range(self.nVertices):
                sum = sum + self.Vertices[v][p] 
            self.Psum.append(sum)
     
     
    def ReplaceWorst(self, Better):

        for p in range(self.nParams):  # keep the Psum list up to date
            self.Psum[p] = self.Psum[p] - self.Vertices[self.WorstIdx][p] + Better[p]
         
        del self.Vertices[self.WorstIdx]
        
        try:
            v = 0                                             # find the point to insert Better 
            while self.Vertices[v].Value < Better.Value: # keeping self.Vertices sorted 
                v = v + 1
            self.Vertices.insert(v, Better)
        except:
            self.Vertices.append(Better)


    def AwayFromWorst(self, Distance):
        """ reflect worst point through face of simplex by Distance and replace if better,
        Equivalent to the "amotry" routine in NRiC. """
     
        fac1 = (1.0 - Distance)/self.nParams
        fac2 = fac1 - Distance

        TryVert = self.Vertices[self.WorstIdx].Copy()
        WorstVal = TryVert.Value

        for p in range(TryVert.nParams):
            TryVert[p]= self.Psum[p] * fac1 - self.Vertices[self.WorstIdx][p]*fac2


        TryVert.Eval(**self.fargs)
        if TryVert.Value < WorstVal:
            self.ReplaceWorst(TryVert)

        return TryVert.Value
     

    def ContractToBest(self):
     
         fac = 0.5
         BestV = self.Vertices[0]
         for v in self.Vertices[1:]:
             for p in range(self.nParams):
                 v[p] = fac*(v[p]+BestV[p])
             v.Eval(**self.fargs)
         self.Vertices.sort()
         self.RecalcPsum()
     

    def Run(self):
        """ do the downhil simplex until self.StopFun() == 1 """
         
        Alpha = -1
        Beta = 0.5
        Gamma = 2.0
        self.Iters = 0
         
        while not self.StopFunc():
            
            BestVal = self.Vertices[self.BestIdx].Value
            WorstVal = self.Vertices[self.WorstIdx].Value
            Worst2Val = self.Vertices[self.Worst2Idx].Value
            
            TryVal = self.AwayFromWorst(Alpha)
            if TryVal <= BestVal:               # better than prev. best
                TryVal = self.AwayFromWorst(Gamma)      # so go some more in same direction
            elif TryVal >= Worst2Val:
                TryVal = self.AwayFromWorst(Beta) # try 1 dim contraction
                if TryVal >= WorstVal:          # still no improvement
                    self.ContractToBest()    
                    
    def Randomise(self,fac):
        for v in self.Vertices:
            v.Randomise(fac)