Changeset 608 for branches

Timestamp:

12/10/12 11:53:49 (3 years ago)

Author:

mmckerns

Message:

added solve, as_constraint, and as_penalty;
added tests for solve, issolution, as_penalty, and as_constraint;
issolution now works for penalty functions and constraints solvers;
fix for randomstate and a first better handling of solver init in sumt;

Location:

branches/decorate

Files:

• constraints.py (modified) (1 diff)
• penalty.py (modified) (1 diff)
• restarts.py (modified) (3 diffs)
• test_penalty.py (modified) (2 diffs)
• test_restarts.py (modified) (5 diffs)
• wrapper.py (modified) (1 diff)

branches/decorate/constraints.py

@@ -76 +76 @@
 
 
+# EOF

branches/decorate/penalty.py

@@ -419 +419 @@
 
 
-def issolution(penalty, candidate, tol=1e-3):
-    """check if the candidate is a solution to penalty constraints
-
-penalty: a penalty function
-candidate: a candidate solution
-    """
-    if penalty.error(candidate) <= tol: return True
+def issolution(constraints, candidate, tol=1e-3):
+    """Returns whether the candidate is a solution to the constraints
+
+Input:
+    constraints -- a constraints solver function or a penalty function
+    candidate -- list of parameter values prposed to solve the constraints
+    tol -- residual error magnitude for which constraints are considered solved
+    """
+    if hasattr(constraints, 'error'):
+        error = constraints.error(candidate)
+    else: # is a constraints solver
+        error = 0.0
+        constrained = constraints(candidate)
+        for i in range(len(candidate)):
+            error += (constrained[i] - candidate[i])**2  #XXX: better rnorm ?
+        error = error**0.5
+
+    if error <= tol: return True
     return False
 
+
 #XXX: nice if penalty.error could give error for each condition... or total
 
 
+def solve(constraints, guess=None, nvars=None, solver=None, \
+          lower_bounds=None, upper_bounds=None, termination=None):
+    """Use optimization to find a solution to a set of constraints.
+
+Inputs:
+    constraints -- a constraints solver function or a penalty function
+
+Additional Inputs:
+    guess -- list of parameter values proposed to solve the constraints.
+    lower_bounds -- list of lower bounds on solution values.
+    upper_bounds -- list of upper bounds on solution values.
+    nvars -- number of parameter values.
+    solver -- the mystic solver to use in the optimization
+    termination -- the mystic termination to use in the optimization
+
+NOTE: The resulting constraints will likely be more expensive to evaluate
+    and less accurate than writing the constraints solver from scratch.
+    """
+    ndim = 1 #XXX: better, increase in while loop catching IndexError ?
+    if nvars is not None: ndim = nvars
+    elif guess is not None: ndim = len(guess)
+    elif lower_bounds is not None: ndim = len(lower_bounds)
+    elif upper_bounds is not None: ndim = len(upper_bounds)
+
+    def cost(x): return 1.
+
+    #XXX: don't allow solver string as a short-cut?
+    if solver is None or solver == 'diffev':
+        from mystic.solvers import DifferentialEvolutionSolver as TheSolver
+        solver = TheSolver(ndim, min(40, ndim*5))
+    elif solver == 'diffev2':
+        from mystic.solvers import DifferentialEvolutionSolver2 as TheSolver
+        solver = TheSolver(ndim, min(40, ndim*5))
+    elif solver == 'fmin_powell': #XXX: better as the default? (it's not random)
+        from mystic.solvers import PowellDirectionalSolver as TheSolver
+        solver = TheSolver(ndim)
+    elif solver == 'fmin':
+        from mystic.solvers import NelderMeadSimplexSolver as TheSolver
+        solver = TheSolver(ndim)
+    
+    if termination is None:
+        from mystic.termination import ChangeOverGeneration as COG
+        termination = COG()
+    if guess != None:
+        solver.SetInitialPoints(guess) #XXX: nice if 'diffev' also had methods
+    else:
+        solver.SetRandomInitialPoints(lower_bounds, upper_bounds)
+    if lower_bounds or upper_bounds:
+        solver.SetStrictRanges(lower_bounds, upper_bounds)
+    if hasattr(constraints, 'iter') and hasattr(constraints, 'error'):
+        solver.SetPenalty(constraints) #i.e. is a penalty function
+    else: # is a constraints solver
+        solver.SetConstraints(constraints)
+    solver.Solve(cost, termination)
+    soln = solver.bestSolution
+
+    from numpy import ndarray, array
+    if isinstance(soln, ndarray) and not isinstance(guess, ndarray):
+        soln = soln.tolist()
+    elif isinstance(guess, ndarray) and not isinstance(soln, ndarray):
+        soln = array(soln)  #XXX: or always return a list ?
+
+    return soln #XXX: check with 'issolution' ?
+
+
+def as_constraint(penalty, *args, **kwds):
+    """Convert a penalty function to a constraints solver.
+
+Inputs:
+    penalty -- a penalty function
+
+Additional Inputs:
+    lower_bounds -- list of lower bounds on solution values.
+    upper_bounds -- list of upper bounds on solution values.
+    nvars -- number of parameter values.
+    solver -- the mystic solver to use in the optimization
+    termination -- the mystic termination to use in the optimization
+    """
+    def constraint(x): #XXX: better to enable args kwds for penalty ?
+        return solve(penalty, x, *args, **kwds)
+    return constraint
+
+
+def as_penalty(constraint, ptype=None, *args, **kwds):
+    """Convert a constraints solver to a penalty function.
+
+Inputs:
+    constraint -- a constraints solver
+    ptype -- penalty function type [default: quadratic_equality]
+
+Additional Inputs:
+    args -- arguments for the constraints solver [default: ()]
+    kwds -- keyword arguments for the constraints solver [default: {}]
+    k -- penalty multiplier
+    h -- iterative multiplier
+    """
+    def rnorm(x, *argz, **kwdz):
+        error = 0.0
+        constrained = constraint(x, *argz, **kwdz)
+        for i in range(len(x)):
+            error += (constrained[i] - x[i])**2  #XXX: better rnorm ?
+        error = error**0.5
+        return error
+
+    if ptype is None: ptype = quadratic_equality
+
+    @ptype(rnorm, *args, **kwds) #XXX: yes to k,h... but otherwise?
+    def penalty(x):
+        return 0.0
+
+    return penalty
+
+
 # EOF 

branches/decorate/restarts.py

@@ -1 @@
-# XXX: solver restart live, from log... restart from saved population as initial
-# XXX: extract solver restart loop for penalties from constraints.py line 1785
-# XXX: extract solver restart loop for collapse from collapse_code.py line 67
-
-# XXX: capture and restart from full population
-# XXX: restart from RandomInputs + bestSolution
-
-#def live_restart(solver):
-#def log_restart(logs):
-#def penalty_restart(solver):
-#def collapse_restart(solver):
-
 
 def sumt(cost, solver, term, eps = 1e-4, max_SUMT_iters = 10, **kwds):
@@ -44 +32 @@
 
         if randomstate:
+            import random
             random.setstate(randomstate)
 
         # 'Clear' some of the attributes of the solver
-        solver.population = [[0.0 for i in range(solver.nDim)] for j in range(solver.nPop)]  #XXX no
+        solver.population = [[0.0 for i in range(solver.nDim)] for j in range(solver.nPop)]
         solver.generations = 0
         solver.bestEnergy = 0.0
@@ -53 +42 @@
         solver.trialSolution = [0.0] * solver.nDim
         solver._init_popEnergy  = inf
-        solver.popEnergy = [solver._init_popEnergy] * solver.nPop  #XXX: no
+        solver.popEnergy = [solver._init_popEnergy] * solver.nPop
         solver.energy_history = []
 
-       #XXX: needs __init__ to correct the above 'XXX'
-       #dim = solver.nDim
-       #simplex = dim + 1
-       #solver.popEnergy = [0.0] * simplex
-       #solver.population = [[0.0 for i in range(dim)] for j in range(simplex)]
+        # individual solver __init__
+        if solver.__class__.__name__ == 'NelderMeadSimplexSolver':
+            simplex = solver.nDim + 1
+            solver.popEnergy = [solver._init_popEnergy] * simplex
+            solver.population = [[0.0 for i in range(solver.nDim)] for j in range(simplex)]
+        elif solver.__class__.__name__ == 'PowellDirectionalSolver':
+            solver._direc = None
+        elif solver.__class__.__name__ in ['DifferentialEvolutionSolver','DifferentialEvolutionSolver2']:
+            solver.genealogy = [ [] for j in range(solver.nPop)]
+        else: pass
 
         # Perform the 'inner' optimization

branches/decorate/test_penalty.py

@@ -1 +1 @@
 from penalty import *
-#from mystic.math import almostEqual
+from mystic.math import almostEqual
 
 def test_penalize():
@@ -29 +29 @@
 
 
+def test_solve():
+
+  from mystic.math.measures import mean
+  def mean_constraint(x, target):
+    return mean(x) - target
+
+  def parameter_constraint(x):
+    return x[-1] - x[0]
+
+  @quadratic_equality(condition=mean_constraint, kwds={'target':5.0})
+  @quadratic_equality(condition=parameter_constraint)
+  def penalty(x):
+    return 0.0
+
+  x = solve(penalty, guess=[2,3,1])
+
+  assert round(mean_constraint(x, 5.0)) == 0.0
+  assert round(parameter_constraint(x)) == 0.0
+  assert issolution(penalty, x)
+
+
+def test_solve_constraint():
+
+  from mystic.math.measures import mean
+  from constraints import with_mean
+
+  @with_mean(1.0)
+  def constraint(x):
+    x[-1] = x[0]
+    return x
+
+  x = solve(constraint, guess=[2,3,1])
+
+  assert almostEqual(mean(x), 1.0, tol=1e-15)
+  assert x[-1] == x[0]
+  assert issolution(constraint, x)
+
+
+def test_as_constraint():
+
+  from mystic.math.measures import mean, spread
+  def mean_constraint(x, target):
+    return mean(x) - target
+
+  def range_constraint(x, target):
+    return spread(x) - target
+
+  @quadratic_equality(condition=range_constraint, kwds={'target':5.0})
+  @quadratic_equality(condition=mean_constraint, kwds={'target':5.0})
+  def penalty(x):
+    return 0.0
+
+  ndim = 3
+  constraints = as_constraint(penalty, solver='fmin_powell')
+  #XXX: this is expensive to evaluate, as there are nested optimizations
+
+  from numpy import arange
+  x = arange(ndim)
+  _x = constraints(x)
+  
+  assert round(mean(_x)) == 5.0
+  assert round(spread(_x)) == 5.0
+
+  def cost(x):
+    return abs(sum(x) - 5.0)
+
+  npop = ndim*3
+  from mystic.solvers import fmin_powell, diffev
+  y = diffev(cost, x, npop, constraints=constraints, disp=False, gtol=10)
+
+  assert round(mean(y)) == 5.0
+  assert round(spread(y)) == 5.0
+  assert round(cost(y)) == 5.0*(ndim-1)
+
+
+def test_as_penalty():
+
+  from mystic.math.measures import mean, spread
+  from constraints import with_mean, with_spread
+
+  @with_spread(5.0)
+  @with_mean(5.0)
+  def constraint(x):
+    return x
+
+  penalty = as_penalty(constraint)
+
+  from numpy import array
+  x = array([1,2,3,4,5])
+  
+  def cost(x):
+    return abs(sum(x) - 5.0)
+
+  from mystic.solvers import fmin_powell
+  y = fmin_powell(cost, x, penalty=penalty, disp=False)
+
+  assert round(mean(y)) == 5.0
+  assert round(spread(y)) == 5.0
+  assert round(cost(y)) == 4*(5.0)
+
+
 
 if __name__ == '__main__':
   test_penalize()
+  test_solve()
+  test_solve_constraint()
+  test_as_constraint()
+  test_as_penalty()
 
 

branches/decorate/test_restarts.py

@@ -5 +5 @@
 random_seed(123)
 
-def test_sumt1(verbose=False):
-
+def test_sumt1(solver, verbose=False):
     def constraint1(x):
         x1,x2 = x[0],x[1]
@@ -27 +26 @@
         print "constraints equations:%s" % (constraints_string.rstrip(),)
 
-    def costfunc(x):
+    def cost(x):
         x1 = x[0]
         x2 = x[1]
         return  x1**4 - 2.*x1**2*x2 + x1**2 + x1*x2**2 - 2.*x1 + 4.
 
-    ndim = 2
-    x0 = [3., 2.]
-    npop = 25
-    from mystic.solvers import DifferentialEvolutionSolver
-    from mystic.solvers import NelderMeadSimplexSolver
+    #XXX: Rand1Bin needs randomstate to work properly when no random_seed ?
+   #import random
+   #randomstate = random.getstate() #XXX: this should be an internal
+    solver.SetPenalty(penalty)
     from mystic.termination import VTR
-    #solver = DifferentialEvolutionSolver(ndim, npop)
-    solver = NelderMeadSimplexSolver(ndim)
-    solver.SetInitialPoints(x0)
-    solver.SetPenalty(penalty)
     term = VTR()
-    end_solver = sumt(costfunc, solver, term, disp=verbose)
-    soln = end_solver.Solution()
+    _solver = sumt(cost, solver, term, disp=verbose)#, randomstate=randomstate)
+    soln = _solver.Solution()
     satisfied = issolution(penalty, soln)
     if verbose: 
@@ -53 +47 @@
 
 
-def test_sumt2(verbose=False):
+def test_sumt2(solver, verbose=False):
     def constraint1(x):
         x1,x2,x3 = x[0],x[1],x[2]
@@ -80 +74 @@
         print "constraints equations:%s" % (constraints_string.rstrip(),)
 
-    def costfunc(x):
+    def cost(x):
         return (x[0] - 1.)**2 + (x[1]-2.)**2 + (x[2]-3.)**4
 
-    ndim = 3
-    x0 = [0.]*ndim
-    from mystic.solvers import DifferentialEvolutionSolver
-    from mystic.solvers import NelderMeadSimplexSolver
+    #XXX: Rand1Bin needs randomstate to work properly when no random_seed ?
+   #import random
+   #randomstate = random.getstate() #XXX: this should be an internal
+    solver.SetPenalty(penalty)
     from mystic.termination import VTR
-    #solver = DifferentialEvolutionSolver(ndim, npop)
-    solver = NelderMeadSimplexSolver(ndim)
-    solver.SetInitialPoints(x0)
-    solver.SetPenalty(penalty)
     term = VTR()
-    end_solver = sumt(costfunc, solver, term, disp=verbose)
-    soln = end_solver.Solution()
+    _solver = sumt(cost, solver, term, disp=verbose)#, randomstate=randomstate)
+    soln = _solver.Solution()
     satisfied = issolution(penalty, soln)
     if verbose:
@@ -104 +94 @@
 
 if __name__ == '__main__':
-    test_sumt1()
-    test_sumt2()
+    from mystic.solvers import DifferentialEvolutionSolver
+    from mystic.solvers import DifferentialEvolutionSolver2
+    from mystic.solvers import PowellDirectionalSolver
+    from mystic.solvers import NelderMeadSimplexSolver
+    verbose = True
+
+    x0 = [3., 2.]; ndim = 2; npop = 25
+   #solver = DifferentialEvolutionSolver(ndim, npop)
+   #solver = DifferentialEvolutionSolver2(ndim, npop)
+   #solver = PowellDirectionalSolver(ndim)
+    solver = NelderMeadSimplexSolver(ndim)
+    solver.SetInitialPoints(x0)
+    test_sumt1(solver, verbose=verbose)
+
+    x0 = [0., 0., 0.]; ndim = 3; npop = 40
+   #solver = DifferentialEvolutionSolver(ndim, npop)
+   #solver = DifferentialEvolutionSolver2(ndim, npop)
+   #solver = PowellDirectionalSolver(ndim)
+    solver = NelderMeadSimplexSolver(ndim)
+    solver.SetInitialPoints(x0)
+    test_sumt2(solver, verbose=verbose)
 
 

branches/decorate/wrapper.py

@@ -1 +1 @@
 #XXX: be mindful if the decorators restrict to functions that expect arrays
 #     compare against the originals for restrictions
+
+#XXX: when_decorated registers methods to 'populate up' when is decorated ?
+
 
 def wrap(outer=lambda x:x, args=None, kwds=None): #XXX: *args, **kwds ?