Changeset 295

Timestamp:

07/07/10 10:46:39 (6 years ago)

Author:

altafang

Message:

Added a function for forming nonlinear constraints by simplifying the nonlinear system using substitution (using Sympy). Still needs more work, but it works for simple test cases.

File:

• branches/alta/mystic-0.2a1/constraint_tools.py (modified) (9 diffs)

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

@@ -230 +230 @@
         return True
 
-def verify_constraints_satisfied(x, constraints, tol=1e-5):
+def verify_constraints_consistent(x, constraints, tol=1e-5):
     """Returns True if x satisfies the constraints function and False if it
  does not."""
@@ -330 +330 @@
             print "Target variable is invalid. Returning None."
             return None
-        return str(soln[0])
+        return target + ' = ' + str(soln[0])
 
 def form_constraints_nonlinear(equations_string, ndim, x0=None, \
@@ -376 +376 @@
                 x[i] = random.uniform(min[i], max[i])
             try:
-                consistent = verify_constraints_satisfied(cf(x), cf)
+                consistent = verify_constraints_consistent(cf(x), cf)
                 return consistent
             except ZeroDivisionError:
@@ -408 +408 @@
                     consistent = verify_consistency_unknownx0(cf)
                 else:
-                    verify_constraints_satisfied(cf(x0), cf)
+                    verify_constraints_consistent(cf(x0), cf)
                 compatible = True
                 if strict:
@@ -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. 
+"""
+    eqns = equations_string.splitlines()
+    # Remove empty strings:
+    actual_eqns = []
+    for j in range(len(eqns)):
+        if eqns[j]:
+           actual_eqns.append(eqns[j])
+
+    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)
+
+
+def verify_constraints_satisfied(constraints_string, x, ndim, \
+                                        varname='x', disp=True):
+    """For debugging. Makes sure that constraints are *really* satisified.
+constraints_string is the constraints string
+x is cf(x), i.e. some array that is supposed to satisfy the constraints.
+ndim is number of dimensions.
+varname is variable name.
+disp is False if you don't want it to print each equation."""
+    from mystic.math import approx_equal
+    for i in range(ndim):
+        variable = varname + str(i+1)
+        constraints_string = constraints_string.replace(variable, 'x[' + str(i) + ']')
+    constraints_list = constraints_string.splitlines()
+
+    # Remove empty strings:
+    actual_eqns = []
+    for j in range(len(constraints_list)):
+        if constraints_list[j]:
+           actual_eqns.append(constraints_list[j])
+
+    for item in actual_eqns:
+        split = item.split('=')
+        if disp:
+            print eval(split[0]), ' ?= ', eval(split[1])
+        if not approx_equal(eval(split[0]), eval(split[1])):
+            return False
+    return True
+
 #-------------------------------------------------------------------
 
@@ -542 +632 @@
     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_function(string2, len(x0), x0=x0, \
+    #                                lower_bounds=lb, upper_bounds=ub)
+    cf = form_constraints_nonlinear_sub(string, len(x0))
     c = cf(x0)
-    print 'constraints satisfied?', verify_constraints_satisfied(c, cf, tol=1e-3)
+    print 'constraints satisfied?', verify_constraints_satisfied(string, c, len(x0))
     print 'Strict: constraints(x0) = ', c
 
@@ -558 +649 @@
     cf = form_constraints_function(string, 3)
     c = cf(x0)
-    print 'constraints satisfied?', verify_constraints_satisfied(c, cf, tol=1e-3)
+    print 'constraints satisfied?', verify_constraints_satisfied(string, c, len(x0))
     print 'Strict: constraints(x0) = ', c
 
@@ -571 +662 @@
     #print 'cf([1.]*4)', cf([1.]*4)
     #print 'cf(cf([1.]*4))', cf(cf([1.]*4))
-    cf = form_constraints_function(eqnstring, 4)
+    cf = form_contsraints_nonlinear_sub(eqnstring, 4)#form_constraints_function(eqnstring, 4)
     print 'cf([1.]*4)', cf([1.]*4)
     print 'cf(cf([1.]*4))', cf(cf([1.]*4))
-    
+
+def test8():
+    # 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.
+    # stringlinear works fine, and can be compared to form_constraints_linear.
+    # string2 works fine.
+    string = """
+x2*x3 = 1. + x1*0.9
+x3 = x1/x2 - 3. + x4
+x4 = 2.*x3 + x1
+"""
+    stringlinear = """
+x2 + 0.5*x1 = x3*2. + x4 + 2.
+5.*x3 - x2 = x1 + x4*3.
+"""
+    string2 = """
+x1 = x3*x2 + 3.
+x2/x4 = 5. 
+x3 + 4. = 0.5*x1
+"""
+    x0 = [0.8,1.2,-0.7, 1.3]
+    #cf=form_constraints_function(stringlinear, len(x0))
+    cf = form_constraints_nonlinear_sub(string2, len(x0))
+    c = cf(x0)
+    print 'consistent?', verify_constraints_consistent(c, cf)
+    print 'satisfied?', verify_constraints_satisfied(string2, c, len(x0))
 
 if __name__ == '__main__':
@@ -581 +698 @@
     #test3()
     #test4()
-    #test5()
+    test5()
     #test6()
-    test7()
+    #test7()
+    #test8()