#import optimizer
import Numeric as N
import LinearAlgebra as LA

#N_MAX = 100
reflect = -1.
halfway = .5
extrapolate = 2.

"""
default test function
"""
def test_function(x):
    y = 42
    for i in range(len(x)):
        y += (x[i]-(i+1))**2
    return y

"""
Performs the 'amoeba'-like downhill simplex method in dimcount dimensions.  used Numerical Recipes in C.
"""
class Amoeba:

    def __init__(self, simplex, values=None, function=test_function,
                 tolerance=0.001, savedata=False):
        y = values
        self.num_vertices = len(simplex)
        self.dimcount = self.num_vertices - 1

        self.simplex = simplex

        self.relativedeviation=0.

        if y is None:
            y = map(function, simplex)
        elif len(simplex) != len(y):
            raise Exception('Vertex count differs from value count')
        self.y = y

        self.evaluationcount = 0
        self.tolerance = tolerance
        self.function = function
        #it probably has to do with the geometric center of the simplex
        self.coord_sums = [None]*self.dimcount
        self.calc_coord_sums()
        self.analyzepoints()

    def optimize(self):
        while self.step() > self.tolerance:
            pass

    def analyzepoints(self):
        simplex = self.simplex
        y = self.y
        num_vertices = self.num_vertices
        dimcount = self.dimcount    
        
        ilow = 0 #index of lowest value
        ihigh = None #index of highest value
        i2ndhigh = None #index of second highest value
        if y[0] > y[1]:
            (ihigh, i2ndhigh) = (0, 1)
        else:
            (ihigh, i2ndhigh) = (1, 0)

        #Loop through vertices to find index values for highest/lowest entries
        for i in range(num_vertices):
            if y[i] < y[ilow]:
                ilow = i
            if y[i] > y[ihigh]:
                i2ndhigh = ihigh
                ihigh = i
            elif y[i] > y[i2ndhigh]:
                if i != ihigh:
                    i2ndhigh = i

        #Things should be floats already, but it's good to be safe
        self.ilow = ilow
        self.ihigh = ihigh
        self.i2ndhigh = i2ndhigh
        self.relativedeviation = float(abs(y[ihigh] - y[ilow]))/abs(y[ihigh]
                                                                    +y[ilow])

        return self.relativedeviation
    
    def calc_coord_sums(self):
        """
        Given a list of (dimcount+1) vectors each with dimcount
        coordinates, returns the list of coordinate sums across
        vectors, i.e. the n'th element is the sum of the n'th
        coordinates of all vectors in p
        """
        for i in range(self.dimcount):
            self.coord_sums[i] = sum([q[i] for q in self.simplex])

    def step(self):
        simplex = self.simplex
        y = self.y
        dimcount = self.dimcount
        num_vertices = self.num_vertices
        coord_sums = self.coord_sums
        function = self.function
        ihigh = self.ihigh
        i2ndhigh = self.i2ndhigh
        ilow = self.ilow
        
        yTry = self.grope(reflect)
        #self.evaluationcount += 1

        if yTry <= y[ilow]:
            yTry = self.grope(extrapolate)
            #self.evaluationcount += 1
        elif yTry >= y[i2ndhigh]:
            ySave = y[ihigh]
            yTry = self.grope(halfway)
            #self.evaluationcount += 1
            if yTry >= ySave:
                for i in range(num_vertices):
                    if i != ilow:
                        for j in range(dimcount):
                            coord_sums[j] = .5 * (simplex[i][j] +
                                                 simplex[ilow][j])
                            simplex[i][j] = coord_sums[j]
                        y[i] = function(coord_sums)
                self.evaluationcount += dimcount
                coord_sums = self.calc_coord_sums()

        return self.analyzepoints()
        
                            
    """
    Extrapolates through or partway to simplex face
    """
    def grope(self, factor):
        y = self.y
        ihigh = self.ihigh
        dimcount = self.dimcount
        simplex = self.simplex
        coord_sums = self.coord_sums
        factor1 = (1. - factor)/dimcount
        factor2 = factor1 - factor

        ptrial = [coord_sums[j]*factor1 - simplex[ihigh][j]*factor2
                  for j in range(dimcount)]

        ytrial = self.function(ptrial)

        if ytrial < y[ihigh]:
            y[ihigh] = ytrial
            for j in range(dimcount):
                coord_sums[j] += ptrial[j] - simplex[ihigh][j]
                simplex[ihigh][j] = ptrial[j]

        self.evaluationcount += 1
        
        return ytrial

def center(simplex):
    vertices = [N.array(point) for point in simplex]
    return sum(vertices)/len(simplex)

def volume(simplex):
    vertices = [N.array(point) for point in simplex]
    differences = [vertices[i]-vertices[i+1] for i in range(len(vertices)-1)]
    return LA.determinant(N.array(differences))

def main():
    simplex = optimizer.get_random_simplex([4,2,1,5,2])
    
    amoeba = Amoeba(simplex, tolerance=0.000001, savedata=True)
    
    amoeba.optimize()
    #print amoeba.simplex
    #print amoeba.y
    return amoeba

if __name__ == '__main__':
     main()