Changeset 297

Timestamp:

07/07/10 12:59:19 (6 years ago)

Author:

altafang

Message:

Added ability for nonlinear constraints function generator to try simplifying in different ways until bounds are satisfied.

Location:

branches/alta/mystic-0.2a1

Files:

• constraint_tools.py (modified) (11 diffs)
• differential_evolution.py (modified) (8 diffs)

branches/alta/mystic-0.2a1/constraint_tools.py

@@ -332 +332 @@
         return target + ' = ' + str(soln[0])
 
-def form_constraints_nonlinear(equations_string, ndim, x0=None, \
+def form_constraints_nonlinear_old(equations_string, ndim, x0=None, \
                               lower_bounds=None, upper_bounds=None, varname='x'):
     """Use Sympy to simplify/invert equations to find a formulation of 
@@ -473 +473 @@
                 varname=varname)
 
-def form_constraints_nonlinear_sub(equations_string, ndim, varname='x'):
-    """Simplify using substitution, algebraically with help from Sympy.
-Algorithm:
-Same as using linear algebra to simplify a linear system....
-Treat the list of equations as a matrix. 
-"""
+def form_constraints_nonlinear(equations_string, ndim, varname='x', x0=None,\
+                                    lower_bounds=None, upper_bounds=None):
+    """Simplify a linear or nonlinear system using substitution with help from Sympy.
+
+Input parameters:
+equations_string -- A string with constraints equations on each line
+ndim -- Number of parameters.
+
+Optional inputs:
+x0 -- Initial x value.
+lower_bounds -- Lower bounds on x.
+upper_bounds -- Upper bounds on x.
+varname -- A string, default is 'x'.
+
+Returns a constraints function.
+"""
+    # Figure out if trying various permutations is necessary
+    strict = False
+    if (x0 and lower_bounds) or (x0 and upper_bounds):
+        strict = True
+
     eqns = equations_string.splitlines()
     # Remove empty strings:
@@ -485 +500 @@
         if eqns[j]:
            actual_eqns.append(eqns[j])
+    orig_eqns = actual_eqns[:]
 
     neqns = len(actual_eqns)
 
-    # Sort the list actual_eqns so that any equation containing x1 is first, etc.
-    sorted_eqns = []
-    actual_eqns_copy = actual_eqns[:]
-    for i in range(ndim):
-        variable = varname + str(i+1)
-        for eqn in actual_eqns_copy:
-            if eqn.find(variable) == -1:
-                pass
-            else:
-                sorted_eqns.append(eqn)
-                actual_eqns_copy.remove(eqn)
-                break
-    actual_eqns = sorted_eqns
-
-    for i in range(neqns):
-        # Trying to use xi as a pivot. Loop through the equations looking for one
-        # containing xi.
-        target = varname + str(i+1)
-        for eqn in actual_eqns[i:]:
-            invertedstring = invert_equation(eqn, ndim, varname=varname, target=target)
-            if invertedstring:
-                break
-        # substitute into the remaining equations. the equations' order in the list
-        # newsystem is like in a linear coefficient matrix.
-        newsystem = ['']*neqns
-        j = actual_eqns.index(eqn)
-        newsystem[j] = eqn
-        othereqns = actual_eqns[:j] + actual_eqns[j+1:]
-        for othereqn in othereqns:
-            expression = invertedstring.split("=")[1]
-            fixed = othereqn.replace(target, '(' + expression + ')')
-            k = actual_eqns.index(othereqn)
-            newsystem[k] = fixed
-        actual_eqns = newsystem
-
-    # Invert so that it can be fed properly to form_constraints_directly
-    simplified = []
-    for eqn in actual_eqns:
-        target = varname + str(actual_eqns.index(eqn) + 1)
-        simplified.append(invert_equation(eqn, ndim, varname=varname, target=target))
-    #print 'simplified:'
-    #print '\n'.join(simplified)
-
-    # Turn this into a constraints function
-    return form_constraints_directly('\n'.join(simplified), ndim, varname=varname)
+    # Some of the permutations will give the same answer; look into reducing the number
+    # of repeats?
+    for p in permutations(range(ndim)):
+        thisorder = p
+        # Sort the list actual_eqns so that any equation containing x1 is first, etc.
+        sorted_eqns = []
+        actual_eqns_copy = orig_eqns[:]
+        for i in thisorder:#range(ndim):
+            variable = varname + str(i+1)
+            for eqn in actual_eqns_copy:
+                if eqn.find(variable) == -1:
+                    pass
+                else:
+                    sorted_eqns.append(eqn)
+                    actual_eqns_copy.remove(eqn)
+                    break
+        if actual_eqns_copy:
+            for item in actual_eqns_copy:
+                sorted_eqns.append(item)
+        actual_eqns = sorted_eqns
+
+        for i in range(neqns):
+            # Trying to use xi as a pivot. Loop through the equations looking for one
+            # containing xi.
+            target = varname + str(thisorder[i] + 1)#str(i+1)
+            for eqn in actual_eqns[i:]:
+                invertedstring = invert_equation(eqn, ndim, varname=varname, target=target)
+                if invertedstring:
+                    break
+            # substitute into the remaining equations. the equations' order in the list
+            # newsystem is like in a linear coefficient matrix.
+            newsystem = ['']*neqns
+            j = actual_eqns.index(eqn)
+            newsystem[j] = eqn
+            othereqns = actual_eqns[:j] + actual_eqns[j+1:]
+            for othereqn in othereqns:
+                expression = invertedstring.split("=")[1]
+                fixed = othereqn.replace(target, '(' + expression + ')')
+                k = actual_eqns.index(othereqn)
+                newsystem[k] = fixed
+            actual_eqns = newsystem
+            
+        # Invert so that it can be fed properly to form_constraints_directly
+        simplified = []
+        for eqn in actual_eqns:
+            target = varname + str(thisorder[actual_eqns.index(eqn)]+1)#actual_eqns.index(eqn) + 1)
+            simplified.append(invert_equation(eqn, ndim, varname=varname, target=target))
+
+        cf = form_constraints_directly('\n'.join(simplified), ndim, varname=varname)  
+
+        if not strict:
+            return cf
+        else:
+            compatible = constraints_inside_bounds(cf, x0, lower_bounds=\
+                 lower_bounds, upper_bounds=upper_bounds)
+            satisfied = verify_constraints_satisfied(equations_string, cf(x0), ndim, disp=False)
+            if compatible and satisfied:
+                return cf
+    print 'Warning! Initial x value cannot satisfy constraints and bounds, \
+or constraints cannot be enforced correctly. Returning a constraints function \
+anyways, but it may not work correctly.'
+    return cf
 
 
@@ -538 +571 @@
                                         varname='x', disp=True):
     """For debugging. Makes sure that constraints are *really* satisified.
-constraints_string is the constraints string
+constraints_string is the constraints string -- Note: not the function!
 x is cf(x), i.e. some array that is supposed to satisfy the constraints.
 ndim is number of dimensions.
@@ -612 +645 @@
 def test5():
     # Test a nonlinear constraints example with x0 and bounds.
-    #XXX This is a case where it does not work properly because the second line
-    # erases the effects of the first line, and since it's nonlinear, Sympy can't
-    # simplify it. If it were linear, form_constraints_linear would pass both
-    # equations simultaneously to Sympy to simplify, and this would not be
-    # a problem. In order to solve something like this, you would have to solve for 
-    # one variable and then plug that into the other equation. Can get bad if there
-    # are many equations and many dimensions though. Maybe I could try to hack 
-    # something together that tries to eliminate variables?? For example, solves for x1,
-    # then plugs that into x1 everywhere....
     string = """
 x2 - 2.*x1*x3 = 1.
-x3 = x1*x2 + 3.
+x3 = x1/x2 + 3.
 """
     # string2 is linear
@@ -630 +654 @@
 2.*x1 - 3.*x3 = 0.5*x2"""
     x0 = [0.8,1.2,-0.7] 
-    lb = [-1., -1., -1.]  
+    lb = [-1., -1., -1.]
     ub = [200., 100., 200.]
-    #cf = form_constraints_function(string2, len(x0), x0=x0, \
-    #                                lower_bounds=lb, upper_bounds=ub)
-    cf = form_constraints_nonlinear_sub(string, len(x0))
+    cf = form_constraints_function(string, len(x0), x0=x0, \
+                                    lower_bounds=lb, upper_bounds=ub)
     c = cf(x0)
-    print 'constraints satisfied?', verify_constraints_satisfied(string, c, len(x0))
+    print 'constraints satisfied?\n', verify_constraints_satisfied(string, c, len(x0))
     print 'Strict: constraints(x0) = ', c
 
@@ -662 +685 @@
     #print 'cf([1.]*4)', cf([1.]*4)
     #print 'cf(cf([1.]*4))', cf(cf([1.]*4))
-    cf = form_contsraints_nonlinear_sub(eqnstring, 4)#form_constraints_function(eqnstring, 4)
+    cf = form_contsraints_nonlinear(eqnstring, 4)#form_constraints_function(eqnstring, 4)
     print 'cf([1.]*4)', cf([1.]*4)
     print 'cf(cf([1.]*4))', cf(cf([1.]*4))
@@ -669 +692 @@
     # Test a nonlinear constraints example with no x0 or bounds and not too many 
     # variables in the same equation.
-    # Notes: string1 works poorly due to a **(1/2) being evaluated as **(0) by python.
+    # string1 now works fine.
     # stringlinear works fine, and can be compared to form_constraints_linear.
     # string2 works fine.
-    string = """
+    string1 = """
 x2*x3 = 1. + x1*0.9
 x3 = x1/x2 - 3. + x4
@@ -680 +703 @@
 x2 + 0.5*x1 = x3*2. + x4 + 2.
 5.*x3 - x2 = x1 + x4*3.
+x1 = x2
 """
     string2 = """
@@ -687 +711 @@
 """
     x0 = [0.8,1.2,-0.7, 1.3]
+    lb = [-5., -5., -5., -5.]  
+    ub = [10., 10., 10., 10.]
     #cf=form_constraints_function(stringlinear, len(x0))
-    cf = form_constraints_nonlinear_sub(string2, len(x0))
+    cf = form_constraints_nonlinear(string1, len(x0), x0=x0, lower_bounds=lb,\
+            upper_bounds=ub)
     c = cf(x0)
     print 'consistent?', verify_constraints_consistent(c, cf)
-    print 'satisfied?', verify_constraints_satisfied(string2, c, len(x0))
+    print 'satisfied?\n', verify_constraints_satisfied(string1, c, len(x0))
 
 if __name__ == '__main__':
@@ -698 +725 @@
     #test3()
     #test4()
-    test5()
+    #test5()
     #test6()
     #test7()
-    #test8()
+    test8()

branches/alta/mystic-0.2a1/differential_evolution.py

@@ -229 +229 @@
         the current parameter vector.  [default = None]
     disp -- non-zero to print convergence messages.
-
+    usebounceback -- True to use 'bounce back'.
         """
         #allow for inputs that don't conform to AbstractSolver interface
@@ -238 +238 @@
         callback=None        #user-supplied function, called after each step
         disp=0               #non-zero to print convergence messages
+        usebounceback = True # True to use 'bounce back'. 
         if kwds.has_key('strategy'): strategy = kwds['strategy']
         if kwds.has_key('CrossProbability'): \
@@ -244 +245 @@
         if kwds.has_key('callback'): callback = kwds['callback']
         if kwds.has_key('disp'): disp = kwds['disp']
+        if kwds.has_key('usebounceback'): usebounceback = kwds['usebounceback']
         #-------------------------------------------------------------
 
@@ -315 +317 @@
 
                 # Use 'bounce back'. Helps a little when soln is on the boundary.
-                self._keepSolutionWithinRangeBoundary(base)
+                if usebounceback:
+                    self._keepSolutionWithinRangeBoundary(base)
 
                 if costfunction(self.constraints(self.trialSolution[:])) != inf:
@@ -461 +464 @@
         encounters a parameter set that satisfies user-supplied constraints.
     disp -- non-zero to print convergence messages.
+    usebounceback -- True to use 'bounce back'.
         """
         #allow for inputs that don't conform to AbstractSolver interface
@@ -469 +473 @@
         callback=None        #user-supplied function, called after each step
         disp=0               #non-zero to print convergence messages
+        usebounceback = True # True to use 'bounce back'. 
         if kwds.has_key('strategy'): strategy = kwds['strategy']
         if kwds.has_key('CrossProbability'): \
@@ -475 +480 @@
         if kwds.has_key('callback'): callback = kwds['callback']
         if kwds.has_key('disp'): disp = kwds['disp']
+        if kwds.has_key('usebounceback'): usebounceback = kwds['usebounceback']
         #-------------------------------------------------------------
 
@@ -561 +567 @@
                 strategy(self, candidate)
 
-                # Use 'bounce back'. Helps a little when soln is on the boundary.
-                self._keepSolutionWithinRangeBoundary(base)
+                if usebounceback:
+                    # Use 'bounce back'. Helps a little when soln is on the boundary.
+                    self._keepSolutionWithinRangeBoundary(base)
 
                 # Only apply constraints if they don't violate bounds.