Changeset 295


Ignore:
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:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/alta/mystic-0.2a1/constraint_tools.py

    r294 r295  
    230230        return True 
    231231 
    232 def verify_constraints_satisfied(x, constraints, tol=1e-5): 
     232def verify_constraints_consistent(x, constraints, tol=1e-5): 
    233233    """Returns True if x satisfies the constraints function and False if it 
    234234 does not.""" 
     
    330330            print "Target variable is invalid. Returning None." 
    331331            return None 
    332         return str(soln[0]) 
     332        return target + ' = ' + str(soln[0]) 
    333333 
    334334def form_constraints_nonlinear(equations_string, ndim, x0=None, \ 
     
    376376                x[i] = random.uniform(min[i], max[i]) 
    377377            try: 
    378                 consistent = verify_constraints_satisfied(cf(x), cf) 
     378                consistent = verify_constraints_consistent(cf(x), cf) 
    379379                return consistent 
    380380            except ZeroDivisionError: 
     
    408408                    consistent = verify_consistency_unknownx0(cf) 
    409409                else: 
    410                     verify_constraints_satisfied(cf(x0), cf) 
     410                    verify_constraints_consistent(cf(x0), cf) 
    411411                compatible = True 
    412412                if strict: 
     
    473473                varname=varname) 
    474474 
     475def form_constraints_nonlinear_sub(equations_string, ndim, varname='x'): 
     476    """Simplify using substitution, algebraically with help from Sympy. 
     477Algorithm: 
     478Same as using linear algebra to simplify a linear system.... 
     479Treat the list of equations as a matrix.  
     480""" 
     481    eqns = equations_string.splitlines() 
     482    # Remove empty strings: 
     483    actual_eqns = [] 
     484    for j in range(len(eqns)): 
     485        if eqns[j]: 
     486           actual_eqns.append(eqns[j]) 
     487 
     488    neqns = len(actual_eqns) 
     489 
     490    # Sort the list actual_eqns so that any equation containing x1 is first, etc. 
     491    sorted_eqns = [] 
     492    actual_eqns_copy = actual_eqns[:] 
     493    for i in range(ndim): 
     494        variable = varname + str(i+1) 
     495        for eqn in actual_eqns_copy: 
     496            if eqn.find(variable) == -1: 
     497                pass 
     498            else: 
     499                sorted_eqns.append(eqn) 
     500                actual_eqns_copy.remove(eqn) 
     501                break 
     502    actual_eqns = sorted_eqns 
     503 
     504    for i in range(neqns): 
     505        # Trying to use xi as a pivot. Loop through the equations looking for one 
     506        # containing xi. 
     507        target = varname + str(i+1) 
     508        for eqn in actual_eqns[i:]: 
     509            invertedstring = invert_equation(eqn, ndim, varname=varname, target=target) 
     510            if invertedstring: 
     511                break 
     512        # substitute into the remaining equations. the equations' order in the list 
     513        # newsystem is like in a linear coefficient matrix. 
     514        newsystem = ['']*neqns 
     515        j = actual_eqns.index(eqn) 
     516        newsystem[j] = eqn 
     517        othereqns = actual_eqns[:j] + actual_eqns[j+1:] 
     518        for othereqn in othereqns: 
     519            expression = invertedstring.split("=")[1] 
     520            fixed = othereqn.replace(target, '(' + expression + ')') 
     521            k = actual_eqns.index(othereqn) 
     522            newsystem[k] = fixed 
     523        actual_eqns = newsystem 
     524 
     525    # Invert so that it can be fed properly to form_constraints_directly 
     526    simplified = [] 
     527    for eqn in actual_eqns: 
     528        target = varname + str(actual_eqns.index(eqn) + 1) 
     529        simplified.append(invert_equation(eqn, ndim, varname=varname, target=target)) 
     530    #print 'simplified:' 
     531    #print '\n'.join(simplified) 
     532 
     533    # Turn this into a constraints function 
     534    return form_constraints_directly('\n'.join(simplified), ndim, varname=varname) 
     535 
     536 
     537def verify_constraints_satisfied(constraints_string, x, ndim, \ 
     538                                        varname='x', disp=True): 
     539    """For debugging. Makes sure that constraints are *really* satisified. 
     540constraints_string is the constraints string 
     541x is cf(x), i.e. some array that is supposed to satisfy the constraints. 
     542ndim is number of dimensions. 
     543varname is variable name. 
     544disp is False if you don't want it to print each equation.""" 
     545    from mystic.math import approx_equal 
     546    for i in range(ndim): 
     547        variable = varname + str(i+1) 
     548        constraints_string = constraints_string.replace(variable, 'x[' + str(i) + ']') 
     549    constraints_list = constraints_string.splitlines() 
     550 
     551    # Remove empty strings: 
     552    actual_eqns = [] 
     553    for j in range(len(constraints_list)): 
     554        if constraints_list[j]: 
     555           actual_eqns.append(constraints_list[j]) 
     556 
     557    for item in actual_eqns: 
     558        split = item.split('=') 
     559        if disp: 
     560            print eval(split[0]), ' ?= ', eval(split[1]) 
     561        if not approx_equal(eval(split[0]), eval(split[1])): 
     562            return False 
     563    return True 
     564 
    475565#------------------------------------------------------------------- 
    476566 
     
    542632    lb = [-1., -1., -1.]   
    543633    ub = [200., 100., 200.] 
    544     cf = form_constraints_function(string2, len(x0), x0=x0, \ 
    545                                     lower_bounds=lb, upper_bounds=ub) 
     634    #cf = form_constraints_function(string2, len(x0), x0=x0, \ 
     635    #                                lower_bounds=lb, upper_bounds=ub) 
     636    cf = form_constraints_nonlinear_sub(string, len(x0)) 
    546637    c = cf(x0) 
    547     print 'constraints satisfied?', verify_constraints_satisfied(c, cf, tol=1e-3) 
     638    print 'constraints satisfied?', verify_constraints_satisfied(string, c, len(x0)) 
    548639    print 'Strict: constraints(x0) = ', c 
    549640 
     
    558649    cf = form_constraints_function(string, 3) 
    559650    c = cf(x0) 
    560     print 'constraints satisfied?', verify_constraints_satisfied(c, cf, tol=1e-3) 
     651    print 'constraints satisfied?', verify_constraints_satisfied(string, c, len(x0)) 
    561652    print 'Strict: constraints(x0) = ', c 
    562653 
     
    571662    #print 'cf([1.]*4)', cf([1.]*4) 
    572663    #print 'cf(cf([1.]*4))', cf(cf([1.]*4)) 
    573     cf = form_constraints_function(eqnstring, 4) 
     664    cf = form_contsraints_nonlinear_sub(eqnstring, 4)#form_constraints_function(eqnstring, 4) 
    574665    print 'cf([1.]*4)', cf([1.]*4) 
    575666    print 'cf(cf([1.]*4))', cf(cf([1.]*4)) 
    576      
     667 
     668def test8(): 
     669    # Test a nonlinear constraints example with no x0 or bounds and not too many  
     670    # variables in the same equation. 
     671    # Notes: string1 works poorly due to a **(1/2) being evaluated as **(0) by python. 
     672    # stringlinear works fine, and can be compared to form_constraints_linear. 
     673    # string2 works fine. 
     674    string = """ 
     675x2*x3 = 1. + x1*0.9 
     676x3 = x1/x2 - 3. + x4 
     677x4 = 2.*x3 + x1 
     678""" 
     679    stringlinear = """ 
     680x2 + 0.5*x1 = x3*2. + x4 + 2. 
     6815.*x3 - x2 = x1 + x4*3. 
     682""" 
     683    string2 = """ 
     684x1 = x3*x2 + 3. 
     685x2/x4 = 5.  
     686x3 + 4. = 0.5*x1 
     687""" 
     688    x0 = [0.8,1.2,-0.7, 1.3] 
     689    #cf=form_constraints_function(stringlinear, len(x0)) 
     690    cf = form_constraints_nonlinear_sub(string2, len(x0)) 
     691    c = cf(x0) 
     692    print 'consistent?', verify_constraints_consistent(c, cf) 
     693    print 'satisfied?', verify_constraints_satisfied(string2, c, len(x0)) 
    577694 
    578695if __name__ == '__main__': 
     
    581698    #test3() 
    582699    #test4() 
    583     #test5() 
     700    test5() 
    584701    #test6() 
    585     test7() 
     702    #test7() 
     703    #test8() 
Note: See TracChangeset for help on using the changeset viewer.