Changeset 295
- Timestamp:
- 07/07/10 10:46:39 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/alta/mystic-0.2a1/constraint_tools.py
r294 r295 230 230 return True 231 231 232 def verify_constraints_ satisfied(x, constraints, tol=1e-5):232 def verify_constraints_consistent(x, constraints, tol=1e-5): 233 233 """Returns True if x satisfies the constraints function and False if it 234 234 does not.""" … … 330 330 print "Target variable is invalid. Returning None." 331 331 return None 332 return str(soln[0])332 return target + ' = ' + str(soln[0]) 333 333 334 334 def form_constraints_nonlinear(equations_string, ndim, x0=None, \ … … 376 376 x[i] = random.uniform(min[i], max[i]) 377 377 try: 378 consistent = verify_constraints_ satisfied(cf(x), cf)378 consistent = verify_constraints_consistent(cf(x), cf) 379 379 return consistent 380 380 except ZeroDivisionError: … … 408 408 consistent = verify_consistency_unknownx0(cf) 409 409 else: 410 verify_constraints_ satisfied(cf(x0), cf)410 verify_constraints_consistent(cf(x0), cf) 411 411 compatible = True 412 412 if strict: … … 473 473 varname=varname) 474 474 475 def form_constraints_nonlinear_sub(equations_string, ndim, varname='x'): 476 """Simplify using substitution, algebraically with help from Sympy. 477 Algorithm: 478 Same as using linear algebra to simplify a linear system.... 479 Treat 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 537 def verify_constraints_satisfied(constraints_string, x, ndim, \ 538 varname='x', disp=True): 539 """For debugging. Makes sure that constraints are *really* satisified. 540 constraints_string is the constraints string 541 x is cf(x), i.e. some array that is supposed to satisfy the constraints. 542 ndim is number of dimensions. 543 varname is variable name. 544 disp 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 475 565 #------------------------------------------------------------------- 476 566 … … 542 632 lb = [-1., -1., -1.] 543 633 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)) 546 637 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)) 548 639 print 'Strict: constraints(x0) = ', c 549 640 … … 558 649 cf = form_constraints_function(string, 3) 559 650 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)) 561 652 print 'Strict: constraints(x0) = ', c 562 653 … … 571 662 #print 'cf([1.]*4)', cf([1.]*4) 572 663 #print 'cf(cf([1.]*4))', cf(cf([1.]*4)) 573 cf = form_con straints_function(eqnstring, 4)664 cf = form_contsraints_nonlinear_sub(eqnstring, 4)#form_constraints_function(eqnstring, 4) 574 665 print 'cf([1.]*4)', cf([1.]*4) 575 666 print 'cf(cf([1.]*4))', cf(cf([1.]*4)) 576 667 668 def 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 = """ 675 x2*x3 = 1. + x1*0.9 676 x3 = x1/x2 - 3. + x4 677 x4 = 2.*x3 + x1 678 """ 679 stringlinear = """ 680 x2 + 0.5*x1 = x3*2. + x4 + 2. 681 5.*x3 - x2 = x1 + x4*3. 682 """ 683 string2 = """ 684 x1 = x3*x2 + 3. 685 x2/x4 = 5. 686 x3 + 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)) 577 694 578 695 if __name__ == '__main__': … … 581 698 #test3() 582 699 #test4() 583 #test5()700 test5() 584 701 #test6() 585 test7() 702 #test7() 703 #test8()
Note: See TracChangeset
for help on using the changeset viewer.