How to use intersect method in fMBT

Best Python code snippet using fMBT_python

test_qp_subproblem.py

Source:test_qp_subproblem.py Github

copy

Full Screen

1import numpy as np2from scipy.sparse import csc_matrix3from scipy.optimize._trustregion_constr.qp_subproblem \4 import (eqp_kktfact,5 projected_cg,6 box_intersections,7 sphere_intersections,8 box_sphere_intersections,9 modified_dogleg)10from scipy.optimize._trustregion_constr.projections \11 import projections12from numpy.testing import (TestCase, assert_array_almost_equal,13 assert_array_equal, assert_array_less,14 assert_equal, assert_,15 run_module_suite, assert_allclose, assert_warns,16 dec)17import pytest18class TestEQPDirectFactorization(TestCase):19 # From Example 16.2 Nocedal/Wright "Numerical20 # Optimization" p.452.21 def test_nocedal_example(self):22 H = csc_matrix([[6, 2, 1],23 [2, 5, 2],24 [1, 2, 4]])25 A = csc_matrix([[1, 0, 1],26 [0, 1, 1]])27 c = np.array([-8, -3, -3])28 b = -np.array([3, 0])29 x, lagrange_multipliers = eqp_kktfact(H, c, A, b)30 assert_array_almost_equal(x, [2, -1, 1])31 assert_array_almost_equal(lagrange_multipliers, [3, -2])32class TestSphericalBoundariesIntersections(TestCase):33 def test_2d_sphere_constraints(self):34 # Interior inicial point35 ta, tb, intersect = sphere_intersections([0, 0],36 [1, 0], 0.5)37 assert_array_almost_equal([ta, tb], [0, 0.5])38 assert_equal(intersect, True)39 # No intersection between line and circle40 ta, tb, intersect = sphere_intersections([2, 0],41 [0, 1], 1)42 assert_equal(intersect, False)43 # Outside inicial point pointing toward outside the circle44 ta, tb, intersect = sphere_intersections([2, 0],45 [1, 0], 1)46 assert_equal(intersect, False)47 # Outside inicial point pointing toward inside the circle48 ta, tb, intersect = sphere_intersections([2, 0],49 [-1, 0], 1.5)50 assert_array_almost_equal([ta, tb], [0.5, 1])51 assert_equal(intersect, True)52 # Inicial point on the boundary53 ta, tb, intersect = sphere_intersections([2, 0],54 [1, 0], 2)55 assert_array_almost_equal([ta, tb], [0, 0])56 assert_equal(intersect, True)57 def test_2d_sphere_constraints_line_intersections(self):58 # Interior inicial point59 ta, tb, intersect = sphere_intersections([0, 0],60 [1, 0], 0.5,61 entire_line=True)62 assert_array_almost_equal([ta, tb], [-0.5, 0.5])63 assert_equal(intersect, True)64 # No intersection between line and circle65 ta, tb, intersect = sphere_intersections([2, 0],66 [0, 1], 1,67 entire_line=True)68 assert_equal(intersect, False)69 # Outside inicial point pointing toward outside the circle70 ta, tb, intersect = sphere_intersections([2, 0],71 [1, 0], 1,72 entire_line=True)73 assert_array_almost_equal([ta, tb], [-3, -1])74 assert_equal(intersect, True)75 # Outside inicial point pointing toward inside the circle76 ta, tb, intersect = sphere_intersections([2, 0],77 [-1, 0], 1.5,78 entire_line=True)79 assert_array_almost_equal([ta, tb], [0.5, 3.5])80 assert_equal(intersect, True)81 # Inicial point on the boundary82 ta, tb, intersect = sphere_intersections([2, 0],83 [1, 0], 2,84 entire_line=True)85 assert_array_almost_equal([ta, tb], [-4, 0])86 assert_equal(intersect, True)87class TestBoxBoundariesIntersections(TestCase):88 def test_2d_box_constraints(self):89 # Box constraint in the direction of vector d90 ta, tb, intersect = box_intersections([2, 0], [0, 2],91 [1, 1], [3, 3])92 assert_array_almost_equal([ta, tb], [0.5, 1])93 assert_equal(intersect, True)94 # Negative direction95 ta, tb, intersect = box_intersections([2, 0], [0, 2],96 [1, -3], [3, -1])97 assert_equal(intersect, False)98 # Some constraints are absent (set to +/- inf)99 ta, tb, intersect = box_intersections([2, 0], [0, 2],100 [-np.inf, 1],101 [np.inf, np.inf])102 assert_array_almost_equal([ta, tb], [0.5, 1])103 assert_equal(intersect, True)104 # Intersect on the face of the box105 ta, tb, intersect = box_intersections([1, 0], [0, 1],106 [1, 1], [3, 3])107 assert_array_almost_equal([ta, tb], [1, 1])108 assert_equal(intersect, True)109 # Interior inicial pointoint110 ta, tb, intersect = box_intersections([0, 0], [4, 4],111 [-2, -3], [3, 2])112 assert_array_almost_equal([ta, tb], [0, 0.5])113 assert_equal(intersect, True)114 # No intersection between line and box constraints115 ta, tb, intersect = box_intersections([2, 0], [0, 2],116 [-3, -3], [-1, -1])117 assert_equal(intersect, False)118 ta, tb, intersect = box_intersections([2, 0], [0, 2],119 [-3, 3], [-1, 1])120 assert_equal(intersect, False)121 ta, tb, intersect = box_intersections([2, 0], [0, 2],122 [-3, -np.inf],123 [-1, np.inf])124 assert_equal(intersect, False)125 ta, tb, intersect = box_intersections([0, 0], [1, 100],126 [1, 1], [3, 3])127 assert_equal(intersect, False)128 ta, tb, intersect = box_intersections([0.99, 0], [0, 2],129 [1, 1], [3, 3])130 assert_equal(intersect, False)131 # Inicial point on the boundary132 ta, tb, intersect = box_intersections([2, 2], [0, 1],133 [-2, -2], [2, 2])134 assert_array_almost_equal([ta, tb], [0, 0])135 assert_equal(intersect, True)136 def test_2d_box_constraints_entire_line(self):137 # Box constraint in the direction of vector d138 ta, tb, intersect = box_intersections([2, 0], [0, 2],139 [1, 1], [3, 3],140 entire_line=True)141 assert_array_almost_equal([ta, tb], [0.5, 1.5])142 assert_equal(intersect, True)143 # Negative direction144 ta, tb, intersect = box_intersections([2, 0], [0, 2],145 [1, -3], [3, -1],146 entire_line=True)147 assert_array_almost_equal([ta, tb], [-1.5, -0.5])148 assert_equal(intersect, True)149 # Some constraints are absent (set to +/- inf)150 ta, tb, intersect = box_intersections([2, 0], [0, 2],151 [-np.inf, 1],152 [np.inf, np.inf],153 entire_line=True)154 assert_array_almost_equal([ta, tb], [0.5, np.inf])155 assert_equal(intersect, True)156 # Intersect on the face of the box157 ta, tb, intersect = box_intersections([1, 0], [0, 1],158 [1, 1], [3, 3],159 entire_line=True)160 assert_array_almost_equal([ta, tb], [1, 3])161 assert_equal(intersect, True)162 # Interior inicial pointoint163 ta, tb, intersect = box_intersections([0, 0], [4, 4],164 [-2, -3], [3, 2],165 entire_line=True)166 assert_array_almost_equal([ta, tb], [-0.5, 0.5])167 assert_equal(intersect, True)168 # No intersection between line and box constraints169 ta, tb, intersect = box_intersections([2, 0], [0, 2],170 [-3, -3], [-1, -1],171 entire_line=True)172 assert_equal(intersect, False)173 ta, tb, intersect = box_intersections([2, 0], [0, 2],174 [-3, 3], [-1, 1],175 entire_line=True)176 assert_equal(intersect, False)177 ta, tb, intersect = box_intersections([2, 0], [0, 2],178 [-3, -np.inf],179 [-1, np.inf],180 entire_line=True)181 assert_equal(intersect, False)182 ta, tb, intersect = box_intersections([0, 0], [1, 100],183 [1, 1], [3, 3],184 entire_line=True)185 assert_equal(intersect, False)186 ta, tb, intersect = box_intersections([0.99, 0], [0, 2],187 [1, 1], [3, 3],188 entire_line=True)189 assert_equal(intersect, False)190 # Inicial point on the boundary191 ta, tb, intersect = box_intersections([2, 2], [0, 1],192 [-2, -2], [2, 2],193 entire_line=True)194 assert_array_almost_equal([ta, tb], [-4, 0])195 assert_equal(intersect, True)196 def test_3d_box_constraints(self):197 # Simple case198 ta, tb, intersect = box_intersections([1, 1, 0], [0, 0, 1],199 [1, 1, 1], [3, 3, 3])200 assert_array_almost_equal([ta, tb], [1, 1])201 assert_equal(intersect, True)202 # Negative direction203 ta, tb, intersect = box_intersections([1, 1, 0], [0, 0, -1],204 [1, 1, 1], [3, 3, 3])205 assert_equal(intersect, False)206 # Interior Point207 ta, tb, intersect = box_intersections([2, 2, 2], [0, -1, 1],208 [1, 1, 1], [3, 3, 3])209 assert_array_almost_equal([ta, tb], [0, 1])210 assert_equal(intersect, True)211 def test_3d_box_constraints_entire_line(self):212 # Simple case213 ta, tb, intersect = box_intersections([1, 1, 0], [0, 0, 1],214 [1, 1, 1], [3, 3, 3],215 entire_line=True)216 assert_array_almost_equal([ta, tb], [1, 3])217 assert_equal(intersect, True)218 # Negative direction219 ta, tb, intersect = box_intersections([1, 1, 0], [0, 0, -1],220 [1, 1, 1], [3, 3, 3],221 entire_line=True)222 assert_array_almost_equal([ta, tb], [-3, -1])223 assert_equal(intersect, True)224 # Interior Point225 ta, tb, intersect = box_intersections([2, 2, 2], [0, -1, 1],226 [1, 1, 1], [3, 3, 3],227 entire_line=True)228 assert_array_almost_equal([ta, tb], [-1, 1])229 assert_equal(intersect, True)230class TestBoxSphereBoundariesIntersections(TestCase):231 def test_2d_box_constraints(self):232 # Both constraints are active233 ta, tb, intersect = box_sphere_intersections([1, 1], [-2, 2],234 [-1, -2], [1, 2], 2,235 entire_line=False)236 assert_array_almost_equal([ta, tb], [0, 0.5])237 assert_equal(intersect, True)238 # None of the constraints are active239 ta, tb, intersect = box_sphere_intersections([1, 1], [-1, 1],240 [-1, -3], [1, 3], 10,241 entire_line=False)242 assert_array_almost_equal([ta, tb], [0, 1])243 assert_equal(intersect, True)244 # Box Constraints are active245 ta, tb, intersect = box_sphere_intersections([1, 1], [-4, 4],246 [-1, -3], [1, 3], 10,247 entire_line=False)248 assert_array_almost_equal([ta, tb], [0, 0.5])249 assert_equal(intersect, True)250 # Spherical Constraints are active251 ta, tb, intersect = box_sphere_intersections([1, 1], [-4, 4],252 [-1, -3], [1, 3], 2,253 entire_line=False)254 assert_array_almost_equal([ta, tb], [0, 0.25])255 assert_equal(intersect, True)256 # Infeasible problems257 ta, tb, intersect = box_sphere_intersections([2, 2], [-4, 4],258 [-1, -3], [1, 3], 2,259 entire_line=False)260 assert_equal(intersect, False)261 ta, tb, intersect = box_sphere_intersections([1, 1], [-4, 4],262 [2, 4], [2, 4], 2,263 entire_line=False)264 assert_equal(intersect, False)265 def test_2d_box_constraints_entire_line(self):266 # Both constraints are active267 ta, tb, intersect = box_sphere_intersections([1, 1], [-2, 2],268 [-1, -2], [1, 2], 2,269 entire_line=True)270 assert_array_almost_equal([ta, tb], [0, 0.5])271 assert_equal(intersect, True)272 # None of the constraints are active273 ta, tb, intersect = box_sphere_intersections([1, 1], [-1, 1],274 [-1, -3], [1, 3], 10,275 entire_line=True)276 assert_array_almost_equal([ta, tb], [0, 2])277 assert_equal(intersect, True)278 # Box Constraints are active279 ta, tb, intersect = box_sphere_intersections([1, 1], [-4, 4],280 [-1, -3], [1, 3], 10,281 entire_line=True)282 assert_array_almost_equal([ta, tb], [0, 0.5])283 assert_equal(intersect, True)284 # Spherical Constraints are active285 ta, tb, intersect = box_sphere_intersections([1, 1], [-4, 4],286 [-1, -3], [1, 3], 2,287 entire_line=True)288 assert_array_almost_equal([ta, tb], [0, 0.25])289 assert_equal(intersect, True)290 # Infeasible problems291 ta, tb, intersect = box_sphere_intersections([2, 2], [-4, 4],292 [-1, -3], [1, 3], 2,293 entire_line=True)294 assert_equal(intersect, False)295 ta, tb, intersect = box_sphere_intersections([1, 1], [-4, 4],296 [2, 4], [2, 4], 2,297 entire_line=True)298 assert_equal(intersect, False)299class TestModifiedDogleg(TestCase):300 def test_cauchypoint_equalsto_newtonpoint(self):301 A = np.array([[1, 8]])302 b = np.array([-16])303 _, _, Y = projections(A)304 newton_point = np.array([0.24615385, 1.96923077])305 # Newton point inside boundaries306 x = modified_dogleg(A, Y, b, 2, [-np.inf, -np.inf], [np.inf, np.inf])307 assert_array_almost_equal(x, newton_point)308 # Spherical constraint active309 x = modified_dogleg(A, Y, b, 1, [-np.inf, -np.inf], [np.inf, np.inf])310 assert_array_almost_equal(x, newton_point/np.linalg.norm(newton_point))311 # Box Constraints active312 x = modified_dogleg(A, Y, b, 2, [-np.inf, -np.inf], [0.1, np.inf])313 assert_array_almost_equal(x, (newton_point/newton_point[0]) * 0.1)314 def test_3d_example(self):315 A = np.array([[1, 8, 1],316 [4, 2, 2]])317 b = np.array([-16, 2])318 Z, LS, Y = projections(A)319 newton_point = np.array([-1.37090909, 2.23272727, -0.49090909])320 cauchy_point = np.array([0.11165723, 1.73068711, 0.16748585])321 origin = np.zeros_like(newton_point)322 # newton_point inside boundaries323 x = modified_dogleg(A, Y, b, 3, [-np.inf, -np.inf, -np.inf],324 [np.inf, np.inf, np.inf])325 assert_array_almost_equal(x, newton_point)326 # line between cauchy_point and newton_point contains best point327 # (spherical constrain is active).328 x = modified_dogleg(A, Y, b, 2, [-np.inf, -np.inf, -np.inf],329 [np.inf, np.inf, np.inf])330 z = cauchy_point331 d = newton_point-cauchy_point332 t = ((x-z)/(d))333 assert_array_almost_equal(t, np.full(3, 0.40807330))334 assert_array_almost_equal(np.linalg.norm(x), 2)335 # line between cauchy_point and newton_point contains best point336 # (box constrain is active).337 x = modified_dogleg(A, Y, b, 5, [-1, -np.inf, -np.inf],338 [np.inf, np.inf, np.inf])339 z = cauchy_point340 d = newton_point-cauchy_point341 t = ((x-z)/(d))342 assert_array_almost_equal(t, np.full(3, 0.7498195))343 assert_array_almost_equal(x[0], -1)344 # line between origin and cauchy_point contains best point345 # (spherical constrain is active).346 x = modified_dogleg(A, Y, b, 1, [-np.inf, -np.inf, -np.inf],347 [np.inf, np.inf, np.inf])348 z = origin349 d = cauchy_point350 t = ((x-z)/(d))351 assert_array_almost_equal(t, np.full(3, 0.573936265))352 assert_array_almost_equal(np.linalg.norm(x), 1)353 # line between origin and newton_point contains best point354 # (box constrain is active).355 x = modified_dogleg(A, Y, b, 2, [-np.inf, -np.inf, -np.inf],356 [np.inf, 1, np.inf])357 z = origin358 d = newton_point359 t = ((x-z)/(d))360 assert_array_almost_equal(t, np.full(3, 0.4478827364))361 assert_array_almost_equal(x[1], 1)362class TestProjectCG(TestCase):363 # From Example 16.2 Nocedal/Wright "Numerical364 # Optimization" p.452.365 def test_nocedal_example(self):366 H = csc_matrix([[6, 2, 1],367 [2, 5, 2],368 [1, 2, 4]])369 A = csc_matrix([[1, 0, 1],370 [0, 1, 1]])371 c = np.array([-8, -3, -3])372 b = -np.array([3, 0])373 Z, _, Y = projections(A)374 x, info = projected_cg(H, c, Z, Y, b)375 assert_equal(info["stop_cond"], 4)376 assert_equal(info["hits_boundary"], False)377 assert_array_almost_equal(x, [2, -1, 1])378 def test_compare_with_direct_fact(self):379 H = csc_matrix([[6, 2, 1, 3],380 [2, 5, 2, 4],381 [1, 2, 4, 5],382 [3, 4, 5, 7]])383 A = csc_matrix([[1, 0, 1, 0],384 [0, 1, 1, 1]])385 c = np.array([-2, -3, -3, 1])386 b = -np.array([3, 0])387 Z, _, Y = projections(A)388 x, info = projected_cg(H, c, Z, Y, b, tol=0)389 x_kkt, _ = eqp_kktfact(H, c, A, b)390 assert_equal(info["stop_cond"], 1)391 assert_equal(info["hits_boundary"], False)392 assert_array_almost_equal(x, x_kkt)393 def test_trust_region_infeasible(self):394 H = csc_matrix([[6, 2, 1, 3],395 [2, 5, 2, 4],396 [1, 2, 4, 5],397 [3, 4, 5, 7]])398 A = csc_matrix([[1, 0, 1, 0],399 [0, 1, 1, 1]])400 c = np.array([-2, -3, -3, 1])401 b = -np.array([3, 0])402 trust_radius = 1403 Z, _, Y = projections(A)404 with pytest.raises(ValueError):405 projected_cg(H, c, Z, Y, b, trust_radius=trust_radius)406 def test_trust_region_barely_feasible(self):407 H = csc_matrix([[6, 2, 1, 3],408 [2, 5, 2, 4],409 [1, 2, 4, 5],410 [3, 4, 5, 7]])411 A = csc_matrix([[1, 0, 1, 0],412 [0, 1, 1, 1]])413 c = np.array([-2, -3, -3, 1])414 b = -np.array([3, 0])415 trust_radius = 2.32379000772445021283416 Z, _, Y = projections(A)417 x, info = projected_cg(H, c, Z, Y, b,418 tol=0,419 trust_radius=trust_radius)420 assert_equal(info["stop_cond"], 2)421 assert_equal(info["hits_boundary"], True)422 assert_array_almost_equal(np.linalg.norm(x), trust_radius)423 assert_array_almost_equal(x, -Y.dot(b))424 def test_hits_boundary(self):425 H = csc_matrix([[6, 2, 1, 3],426 [2, 5, 2, 4],427 [1, 2, 4, 5],428 [3, 4, 5, 7]])429 A = csc_matrix([[1, 0, 1, 0],430 [0, 1, 1, 1]])431 c = np.array([-2, -3, -3, 1])432 b = -np.array([3, 0])433 trust_radius = 3434 Z, _, Y = projections(A)435 x, info = projected_cg(H, c, Z, Y, b,436 tol=0,437 trust_radius=trust_radius)438 assert_equal(info["stop_cond"], 2)439 assert_equal(info["hits_boundary"], True)440 assert_array_almost_equal(np.linalg.norm(x), trust_radius)441 def test_negative_curvature_unconstrained(self):442 H = csc_matrix([[1, 2, 1, 3],443 [2, 0, 2, 4],444 [1, 2, 0, 2],445 [3, 4, 2, 0]])446 A = csc_matrix([[1, 0, 1, 0],447 [0, 1, 0, 1]])448 c = np.array([-2, -3, -3, 1])449 b = -np.array([3, 0])450 Z, _, Y = projections(A)451 with pytest.raises(ValueError):452 projected_cg(H, c, Z, Y, b, tol=0)453 def test_negative_curvature(self):454 H = csc_matrix([[1, 2, 1, 3],455 [2, 0, 2, 4],456 [1, 2, 0, 2],457 [3, 4, 2, 0]])458 A = csc_matrix([[1, 0, 1, 0],459 [0, 1, 0, 1]])460 c = np.array([-2, -3, -3, 1])461 b = -np.array([3, 0])462 Z, _, Y = projections(A)463 trust_radius = 1000464 x, info = projected_cg(H, c, Z, Y, b,465 tol=0,466 trust_radius=trust_radius)467 assert_equal(info["stop_cond"], 3)468 assert_equal(info["hits_boundary"], True)469 assert_array_almost_equal(np.linalg.norm(x), trust_radius)470 # The box constraints are inactive at the solution but471 # are active during the iterations.472 def test_inactive_box_constraints(self):473 H = csc_matrix([[6, 2, 1, 3],474 [2, 5, 2, 4],475 [1, 2, 4, 5],476 [3, 4, 5, 7]])477 A = csc_matrix([[1, 0, 1, 0],478 [0, 1, 1, 1]])479 c = np.array([-2, -3, -3, 1])480 b = -np.array([3, 0])481 Z, _, Y = projections(A)482 x, info = projected_cg(H, c, Z, Y, b,483 tol=0,484 lb=[0.5, -np.inf,485 -np.inf, -np.inf],486 return_all=True)487 x_kkt, _ = eqp_kktfact(H, c, A, b)488 assert_equal(info["stop_cond"], 1)489 assert_equal(info["hits_boundary"], False)490 assert_array_almost_equal(x, x_kkt)491 # The box constraints active and the termination is492 # by maximum iterations (infeasible iteraction).493 def test_active_box_constraints_maximum_iterations_reached(self):494 H = csc_matrix([[6, 2, 1, 3],495 [2, 5, 2, 4],496 [1, 2, 4, 5],497 [3, 4, 5, 7]])498 A = csc_matrix([[1, 0, 1, 0],499 [0, 1, 1, 1]])500 c = np.array([-2, -3, -3, 1])501 b = -np.array([3, 0])502 Z, _, Y = projections(A)503 x, info = projected_cg(H, c, Z, Y, b,504 tol=0,505 lb=[0.8, -np.inf,506 -np.inf, -np.inf],507 return_all=True)508 assert_equal(info["stop_cond"], 1)509 assert_equal(info["hits_boundary"], True)510 assert_array_almost_equal(A.dot(x), -b)511 assert_array_almost_equal(x[0], 0.8)512 # The box constraints are active and the termination is513 # because it hits boundary (without infeasible iteraction).514 def test_active_box_constraints_hits_boundaries(self):515 H = csc_matrix([[6, 2, 1, 3],516 [2, 5, 2, 4],517 [1, 2, 4, 5],518 [3, 4, 5, 7]])519 A = csc_matrix([[1, 0, 1, 0],520 [0, 1, 1, 1]])521 c = np.array([-2, -3, -3, 1])522 b = -np.array([3, 0])523 trust_radius = 3524 Z, _, Y = projections(A)525 x, info = projected_cg(H, c, Z, Y, b,526 tol=0,527 ub=[np.inf, np.inf, 1.6, np.inf],528 trust_radius=trust_radius,529 return_all=True)530 assert_equal(info["stop_cond"], 2)531 assert_equal(info["hits_boundary"], True)532 assert_array_almost_equal(x[2], 1.6)533 # The box constraints are active and the termination is534 # because it hits boundary (infeasible iteraction).535 def test_active_box_constraints_hits_boundaries_infeasible_iter(self):536 H = csc_matrix([[6, 2, 1, 3],537 [2, 5, 2, 4],538 [1, 2, 4, 5],539 [3, 4, 5, 7]])540 A = csc_matrix([[1, 0, 1, 0],541 [0, 1, 1, 1]])542 c = np.array([-2, -3, -3, 1])543 b = -np.array([3, 0])544 trust_radius = 4545 Z, _, Y = projections(A)546 x, info = projected_cg(H, c, Z, Y, b,547 tol=0,548 ub=[np.inf, 0.1, np.inf, np.inf],549 trust_radius=trust_radius,550 return_all=True)551 assert_equal(info["stop_cond"], 2)552 assert_equal(info["hits_boundary"], True)553 assert_array_almost_equal(x[1], 0.1)554 # The box constraints are active and the termination is555 # because it hits boundary (no infeasible iteraction).556 def test_active_box_constraints_negative_curvature(self):557 H = csc_matrix([[1, 2, 1, 3],558 [2, 0, 2, 4],559 [1, 2, 0, 2],560 [3, 4, 2, 0]])561 A = csc_matrix([[1, 0, 1, 0],562 [0, 1, 0, 1]])563 c = np.array([-2, -3, -3, 1])564 b = -np.array([3, 0])565 Z, _, Y = projections(A)566 trust_radius = 1000567 x, info = projected_cg(H, c, Z, Y, b,568 tol=0,569 ub=[np.inf, np.inf, 100, np.inf],570 trust_radius=trust_radius)571 assert_equal(info["stop_cond"], 3)572 assert_equal(info["hits_boundary"], True)...

Full Screen

Full Screen

ResourceAreaCharacterization.py

Source:ResourceAreaCharacterization.py Github

copy

Full Screen

1import config2import arcpy3import time45# This script does intersections of multiple inputs to determine the areas of treated stormwater in a given resource area.6# Output is provided as a time stamped gdb.7# Currently all intermediate steps are saved out as feature classes8# ==> Program is meant to be run using the ArcPro version of Python <==91011def characterize_resource():1213 #create new time stamped gdb14 gdb_name = time.strftime("ResourceAreas_%m_%d_%Y_%H_%M")1516 arcpy.CreateFileGDB_management(config.output_folder, gdb_name, "10.0")17 gdb_path = config.output_folder + r"/" + gdb_name + r".gdb/"1819 working_resource_areas_path = gdb_path + config.resource_area_characterization20 imp_intersect_path = gdb_path + 'imp_intersect'21 imp_intersect_dissolved_path = gdb_path + 'imp_intersect_dissolved'2223 mem = 'memory/'24 #feature_classes_in_memory = []2526 #this copies the resource area and removes extra fields27 arcpy.Dissolve_management(config.resource_areas, gdb_path + config.resource_area_characterization, 'NewSiteNum')2829 working_resource_areas_layer = arcpy.MakeFeatureLayer_management(gdb_path + config.resource_area_characterization, 'working_resource_areas_layer')30 #add all the new fields31 arcpy.AddField_management(working_resource_areas_layer,'Resource_Area_sqft', 'DOUBLE')32 arcpy.CalculateField_management(working_resource_areas_layer,'Resource_Area_sqft',0)3334 arcpy.AddField_management(working_resource_areas_layer,'ImpA_sqft', 'DOUBLE')35 arcpy.CalculateField_management(working_resource_areas_layer,'ImpA_sqft',0)3637 arcpy.AddField_management(gdb_path + config.resource_area_characterization,'Total_Treated_sqft', 'DOUBLE')38 arcpy.CalculateField_management(gdb_path + config.resource_area_characterization,'Total_Treated_sqft',0)3940 arcpy.AddField_management(gdb_path + config.resource_area_characterization,'Total_Untreated_sqft', 'DOUBLE')41 arcpy.CalculateField_management(gdb_path + config.resource_area_characterization,'Total_Untreated_sqft',0)4243 arcpy.AddField_management(gdb_path + config.resource_area_characterization,'ImpA_BMP_Treated_sqft', 'DOUBLE')44 arcpy.CalculateField_management(gdb_path + config.resource_area_characterization,'ImpA_BMP_Treated_sqft',0)4546 arcpy.AddField_management(gdb_path + config.resource_area_characterization,'ImpA_UIC_Treated_sqft', 'DOUBLE')47 arcpy.CalculateField_management(gdb_path + config.resource_area_characterization,'ImpA_UIC_Treated_sqft',0)4849 arcpy.AddField_management(working_resource_areas_layer,'Combined_ImpA_sqft', 'DOUBLE')50 arcpy.CalculateField_management(working_resource_areas_layer, 'Combined_ImpA_sqft', 0)5152 arcpy.AddField_management(gdb_path + config.resource_area_characterization,'Sanitary_ImpA_sqft', 'DOUBLE')53 arcpy.CalculateField_management(gdb_path + config.resource_area_characterization,'Sanitary_ImpA_sqft', 0)5455 ##Impervious Area565758 #intersect impervious areas and resource areas in memory59 #imp_intersect = arcpy.Intersect_analysis([working_resource_areas_path, config.modeled_impervious_areas], imp_intersect_path)60 imp_intersect = arcpy.PairwiseIntersect_analysis([working_resource_areas_path, config.modeled_impervious_areas], imp_intersect_path)6162 #dissolve on NewSiteNum63 imp_intersect_dissolved = arcpy.Dissolve_management(imp_intersect_path, imp_intersect_dissolved_path, 'NewSiteNum')6465 #make a feature layer66 imp_intersect_dissolved_layer = arcpy.MakeFeatureLayer_management(imp_intersect_dissolved, 'imp_intersect_dissolved_layer')6768 # join the impervious area to the resource area and transfer area to 'ImpA_sqft'69 joined = arcpy.AddJoin_management(working_resource_areas_layer, 'NewSiteNum', imp_intersect_dissolved_layer, 'NewSiteNum')70 arcpy.CalculateField_management(working_resource_areas_layer, "ImpA_sqft", "!imp_intersect_dissolved.Shape_Area!")71 arcpy.RemoveJoin_management(working_resource_areas_layer)7273 ##BMP Areas74 #get intersection of bmp and resource areas, dissolve the result, make a feature layer so it can be joined75 bmp_intersect_resource = arcpy.Intersect_analysis([working_resource_areas_layer, config.bmp_areas], gdb_path + 'bmp_intersect_resource')76 bmp_intersect_resource_dissolved = arcpy.Dissolve_management(bmp_intersect_resource, gdb_path + 'bmp_intersect_resource_dissolved', 'NewSiteNum')77 bmp_intersect_resource_dissolved_layer = arcpy.MakeFeatureLayer_management(bmp_intersect_resource_dissolved, 'bmp_intersect_resource_dissolved_layer')7879 ##BMPs intersected with impervious areas80 bmp_intersect_impervious = arcpy.Intersect_analysis([bmp_intersect_resource_dissolved, imp_intersect_dissolved_layer], gdb_path + 'bmp_intersect_impervious')81 bmp_intersect_impervious_dissolved = arcpy.Dissolve_management(bmp_intersect_impervious, gdb_path + 'bmp_intersect_impervious_dissolved', 'NewSiteNum')82 bmp_intersect_impervious_dissolved_layer = arcpy.MakeFeatureLayer_management(bmp_intersect_impervious_dissolved, 'bmp_intersect_impervious_dissolved_layer')8384 ##UIC Areas85 #get intersection of uic and resource areas, dissolve the result, make a feature layer so it can be joined8687 uic_intersect_resource = arcpy.Intersect_analysis([working_resource_areas_layer, config.uic_areas], gdb_path + 'uic_intersect_resource')88 uic_intersect_dissolved = arcpy.Dissolve_management(uic_intersect_resource, gdb_path + 'uic_intersect_dissolved', 'NewSiteNum')89 uic_intersect_dissolved_layer = arcpy.MakeFeatureLayer_management(uic_intersect_dissolved, 'uic_intersect_dissolved_layer')9091 ##UICs intersected with impervious areas92 uic_intersect_impervious = arcpy.Intersect_analysis([uic_intersect_dissolved, imp_intersect_dissolved_layer], gdb_path + 'uic_intersect_impervious')93 uic_intersect_impervious_dissolved = arcpy.Dissolve_management(uic_intersect_impervious, gdb_path + 'uic_intersect_impervious_dissolved', 'NewSiteNum')94 uic_intersect_impervious_dissolved_layer = arcpy.MakeFeatureLayer_management(uic_intersect_impervious_dissolved, 'uic_intersect_impervious_dissolved_layer')959697 arcpy.AddJoin_management(working_resource_areas_layer, 'NewSiteNum', bmp_intersect_impervious_dissolved_layer, 'NewSiteNum')98 arcpy.CalculateField_management(working_resource_areas_layer, "ImpA_BMP_Treated_sqft", "!bmp_intersect_impervious_dissolved.Shape_Area!")99 #copy out the join to see whats up100 arcpy.CopyFeatures_management(working_resource_areas_layer, gdb_path + 'watershed_resource_combined_intersect_dissolved_bmp_join')101 arcpy.RemoveJoin_management(working_resource_areas_layer)102103104 arcpy.AddJoin_management(working_resource_areas_layer, 'NewSiteNum', uic_intersect_impervious_dissolved_layer, 'NewSiteNum')105 arcpy.CalculateField_management(working_resource_areas_layer, "ImpA_UIC_Treated_sqft", "!uic_intersect_impervious_dissolved.Shape_Area!")106 #copy out the join to see whats up107 arcpy.CopyFeatures_management(working_resource_areas_layer, gdb_path + 'watershed_resource_combined_intersect_dissolved_uic_join')108109 arcpy.RemoveJoin_management(working_resource_areas_layer)110111112 ##UICs intersected with impervious areas113 ##Combined Sewer114 #Make layer to filter by sewer area, dissolve, make layer on the dissolved result so it can be joined115 combined_sewer_area = arcpy.MakeFeatureLayer_management(config.sewer_basin_areas, 'combined_sewer_area', "TYPE = 'C'")116 combined_sewer_area_dissolved = arcpy.Dissolve_management(combined_sewer_area, "memory/combined_dissolved", "TYPE")117 combined_sewer_area_dissolved_layer = arcpy.MakeFeatureLayer_management(combined_sewer_area_dissolved, 'combined_sewer_area_dissolved_layer')118119120 #Sanitary Sewer121 #Make layer to filter by sewer area, dissolve, make layer on the dissolved result so it can be joined122 sanitary_sewer_area = arcpy.MakeFeatureLayer_management(config.sewer_basin_areas, 'sanitary_sewer_area', "TYPE = 'S'")123 sanitary_sewer_area_dissolved = arcpy.Dissolve_management(sanitary_sewer_area, "memory/sanitary_dissolved", "TYPE")124 sanitary_sewer_area_dissolved_layer = arcpy.MakeFeatureLayer_management(sanitary_sewer_area_dissolved, 'sanitary_sewer_area_dissolved_layer')125126127 #intersect watershed areas and combined sewer areas128 watershed_resource_combined_intersect = arcpy.Intersect_analysis([working_resource_areas_layer, combined_sewer_area_dissolved_layer], "memory/watershed_resource_combined_intersect")129 watershed_resource_combined_intersect_dissolved = arcpy.Dissolve_management(watershed_resource_combined_intersect, "memory/watershed_resource_combined_intersect_dissolved",'NewSiteNum' )130 arcpy.CopyFeatures_management(watershed_resource_combined_intersect_dissolved, gdb_path + 'watershed_resource_combined_intersect_dissolved')131132 watershed_resource_combined_intersect_dissolved_layer = arcpy.MakeFeatureLayer_management(watershed_resource_combined_intersect_dissolved, 'watershed_resource_combined_intersect_dissolved_layer')133 ##134 #135 #intersect new combined areas with impervious areas and dissolve136 watershed_resource_combined_intersect_impervious = arcpy.Intersect_analysis([watershed_resource_combined_intersect_dissolved, imp_intersect_dissolved_layer], "memory/watershed_resource_combined_intersect_impervious")137 watershed_resource_combined_intersect_impervious_dissolved = arcpy.Dissolve_management(watershed_resource_combined_intersect_impervious, "memory/watershed_resource_combined_intersect_impervious_dissolved", 'NewSiteNum')138 arcpy.CopyFeatures_management(watershed_resource_combined_intersect_impervious_dissolved, gdb_path + 'watershed_resource_combined_intersect_impervious_dissolved')139 watershed_resource_combined_intersect_impervious_dissolved_layer = arcpy.MakeFeatureLayer_management(watershed_resource_combined_intersect_impervious_dissolved, 'watershed_resource_combined_intersect_impervious_dissolved_layer')140141 #intersect watershed areas and sanitary sewer areas142 watershed_resource_sanitary_intersect = arcpy.Intersect_analysis([working_resource_areas_layer, sanitary_sewer_area_dissolved_layer], "memory/watershed_resource_sanitary_intersect")143 watershed_resource_sanitary_intersect_dissolved = arcpy.Dissolve_management(watershed_resource_sanitary_intersect, "memory/watershed_resource_sanitary_intersect_dissolved",'NewSiteNum' )144145 arcpy.CopyFeatures_management(watershed_resource_sanitary_intersect_dissolved, gdb_path + 'watershed_resource_sanitary_intersect_dissolved')146147 watershed_resource_sanitary_intersect_dissolved_layer = arcpy.MakeFeatureLayer_management(watershed_resource_sanitary_intersect_dissolved, 'watershed_resource_sanitary_intersect_dissolved_layer')148 #149 #150 #intersect new sanitary areas with impervious areas and dissolve151 watershed_resource_sanitary_intersect_impervious = arcpy.Intersect_analysis([watershed_resource_sanitary_intersect_dissolved, imp_intersect_dissolved_layer], "memory/watershed_resource_sanitary_intersect_impervious")152 watershed_resource_sanitary_intersect_impervious_dissolved = arcpy.Dissolve_management(watershed_resource_sanitary_intersect_impervious, "memory/watershed_resource_sanitary_intersect_impervious_dissolved", 'NewSiteNum')153 arcpy.CopyFeatures_management(watershed_resource_sanitary_intersect_impervious_dissolved, gdb_path + 'watershed_resource_sanitary_intersect_impervious_dissolved')154 watershed_resource_sanitary_intersect_impervious_dissolved_layer = arcpy.MakeFeatureLayer_management(watershed_resource_sanitary_intersect_impervious_dissolved, 'watershed_resource_sanitary_intersect_impervious_dissolved_layer')155156 #join combined and sanitary impervious areas and transfer the areas157158 arcpy.AddJoin_management(working_resource_areas_layer, 'NewSiteNum', gdb_path + 'watershed_resource_combined_intersect_impervious_dissolved', 'NewSiteNum')159 arcpy.CalculateField_management(working_resource_areas_layer, 'Combined_ImpA_sqft', '!watershed_resource_combined_intersect_impervious_dissolved.Shape_Area!')160 arcpy.RemoveJoin_management(working_resource_areas_layer)161162 arcpy.AddJoin_management(working_resource_areas_layer, 'NewSiteNum', gdb_path + 'watershed_resource_sanitary_intersect_impervious_dissolved', 'NewSiteNum')163 arcpy.CalculateField_management(working_resource_areas_layer, 'Sanitary_ImpA_sqft', '!watershed_resource_sanitary_intersect_impervious_dissolved.Shape_Area!')164 arcpy.RemoveJoin_management(working_resource_areas_layer)165166167168 #prior to calculating total treated, need to merge UIC and BMP areas and then dissolve to remove duplicates169 merged_uic_and_bmp = arcpy.Merge_management([bmp_intersect_impervious_dissolved_layer, uic_intersect_impervious_dissolved_layer], 'memory/merged_uic_and_bmp')170 merged_uic_and_bmp_dissolved = arcpy.Dissolve_management(merged_uic_and_bmp, 'memory/merged_uic_and_bmp_dissolved', 'NewSiteNum')171 merged_uic_and_bmp_dissolved_layer = arcpy.MakeFeatureLayer_management(merged_uic_and_bmp_dissolved, 'merged_uic_and_bmp_dissolved_layer')172 arcpy.CopyFeatures_management(merged_uic_and_bmp_dissolved_layer, gdb_path + 'merged_uic_and_bmp_dissolved')173174 new_join = arcpy.AddJoin_management(working_resource_areas_layer, 'NewSiteNum', gdb_path + 'merged_uic_and_bmp_dissolved', 'NewSiteNum')175 arcpy.CopyFeatures_management(new_join, gdb_path + 'new_join')176177 arcpy.CalculateField_management(working_resource_areas_layer, 'Total_Treated_sqft', '!merged_uic_and_bmp_dissolved.Shape_Area!')178 arcpy.RemoveJoin_management(working_resource_areas_layer)179180 arcpy.CalculateField_management(working_resource_areas_layer, 'Total_Untreated_sqft', "!ImpA_sqft!- !Total_Treated_sqft!", "PYTHON3")181 arcpy.CalculateField_management(working_resource_areas_layer, 'Resource_Area_sqft', "!Shape_Area!", "PYTHON3")182183 #create copies of feature classes184 arcpy.CopyFeatures_management(combined_sewer_area_dissolved_layer, gdb_path + 'combined_sewer_area_dissolved_layer')185 arcpy.CopyFeatures_management(sanitary_sewer_area_dissolved_layer, gdb_path + 'sanitary_sewer_area_dissolved_layer')186187if __name__ == "__main__": ...

Full Screen

Full Screen

Fractal Images & Chaotic Scattering.py

Source:Fractal Images & Chaotic Scattering.py Github

copy

Full Screen

1## This code will model chaotic scattering in a 3 disk system2from numpy import *3import matplotlib.pyplot as plt4import time as time5#########################################################################################################6###### This Portion Is For The RunTimeWarning Created By Complex Numbers in IntersectTimes Function######7#########################################################################################################8import warnings9warnings.simplefilter("ignore", RuntimeWarning)10 11######################################12######FUNCTIONS FOR PROGRAM###########13######################################14def Intersect(xcoordinates,ycoordinates,paramsr,theta,x0,y0): # This function returns the minimum intersection time and circle hit and intersection points15 16 #########FINDING LEAST TIME###############17 18 numintersects = numdisks * 2 # we intersect at each disk twice19 timeintersect = zeros(numintersects)20 21 22 # Unpack parameters we will be using23 24 r1 = paramsr[0]25 r2 = paramsr[1]26 r3 = paramsr[2]27 28 x1 = xcoordinates[0]29 x2 = xcoordinates[1]30 x3 = xcoordinates[2]31 y1 = ycoordinates[0]32 y2 = ycoordinates[1]33 y3 = ycoordinates[2]34 35 # create a,b,c for quadratic formula for finding time36 37 a1 = cos(theta)**2 + sin(theta)**238 a2 = a139 a3 = a140 41 b1 = 2*x0*cos(theta) + 2*y0*sin(theta) - 2*x1*cos(theta) - 2*y1*sin(theta)42 b2 = 2*x0*cos(theta) + 2*y0*sin(theta) - 2*x2*cos(theta) - 2*y2*sin(theta)43 b3 = 2*x0*cos(theta) + 2*y0*sin(theta) - 2*x3*cos(theta) - 2*y3*sin(theta)44 45 c1 = x0**2 + y0**2 + x1**2 + y1**2 - 2*x1*x0 - 2*y1*y0 - r1**246 c2 = x0**2 + y0**2 + x2**2 + y2**2 - 2*x2*x0 - 2*y2*y0 - r2**247 c3 = x0**2 + y0**2 + x3**2 + y3**2 - 2*x3*x0 - 2*y3*y0 - r3**248 49 #first disk intersect points50 51 timeintersect[0] = (-b1 + sqrt(b1**2 - 4*a1*c1))/(2*a1)52 timeintersect[1] = (-b1 - sqrt(b1**2 - 4*a1*c1))/(2*a1)53 54 #second disk intersect points55 56 timeintersect[2] = (-b2 + sqrt(b2**2 - 4*a2*c2))/(2*a2)57 timeintersect[3] = (-b2 - sqrt(b2**2 - 4*a2*c2))/(2*a2)58 59 #third disk intersect points60 61 timeintersect[4] = (-b3 + sqrt(b3**2 - 4*a3*c3))/(2*a3)62 timeintersect[5] = (-b3 - sqrt(b3**2 - 4*a3*c3))/(2*a3)63 64 # We are only concerned with finding the minimum t to the point we are crossing65 66 # If t is complex, we will change that number to a really large one67 # if t is approximately 0 or negative, we will change that to a really large number68 69 reallylargenumber = 10**870 reallysmallnumber = 10 ** -871 72 timeintersect[isnan(timeintersect)] = reallylargenumber 73 timeintersect[timeintersect < reallysmallnumber] = reallylargenumber74 75 IntersectTime = min(timeintersect)76 IntersectCircle = 077 78 if IntersectTime == timeintersect[0] or IntersectTime == timeintersect[1]:79 IntersectCircle = 180 if IntersectTime == timeintersect[2] or IntersectTime == timeintersect[3]:81 IntersectCircle = 282 if IntersectTime == timeintersect[4] or IntersectTime == timeintersect[5]:83 IntersectCircle = 384 if IntersectTime == reallylargenumber:85 IntersectCircle = 086 87 #### FIND INTERSECTING POINTS ####88 89 IntersectPointX = LastImpactX + IntersectTime*cos(theta)90 IntersectPointY = LastImpactY + IntersectTime*sin(theta)91 92 IntersectPoints = [IntersectPointX,IntersectPointY]93 return [IntersectTime,IntersectCircle,IntersectPoints]94def flipvectorgettheta(ImpactPoints,LastImpactX,LastImpactY,IntersectCircle): #this function flips the vector and computes our new angle theta95 96 xp = ImpactPoints[0]97 yp = ImpactPoints[1]98 99 x0 = LastImpactX100 y0 = LastImpactY101 102 103 CircleIntersect = IntersectCircle104 105 if CircleIntersect == 1 : # we hit upper disk106 xc = x1107 yc = y1108 109 if CircleIntersect == 2: # we hit right disk110 xc = x2111 yc = y2112 113 if CircleIntersect == 3: # we hit bottom disk114 xc = x3115 yc = y3116 if CircleIntersect == 0: # We Have Not Hit Any Circles, We return 'NaN' because we cant compute it117 return float('NaN')118 119 incomingx = xp-x0120 incomingy = yp-y0121 122 normalx = xp-xc123 normaly = yp-yc124 125 normalmag = sqrt(normalx**2 + normaly**2)126 127 vectorincoming = array([incomingx,incomingy])128 normal_direction = array([(normalx),(normaly)])129 normal_direction = normal_direction/normalmag130 131 vectorreflect = vectorincoming - 2*normal_direction*(dot(vectorincoming,normal_direction))132 133 ycomponent = vectorreflect[1]134 xcomponent = vectorreflect[0]135 136 theta = arctan(ycomponent/xcomponent)137 138 if xcomponent <= 0.0 and ycomponent >= 0.0:139 theta += pi140 elif xcomponent <= 0.0 and ycomponent <= 0.0:141 theta += -pi142 143 return theta144 145#####################146####Create Disks#####147#####################148numdisks = 3149paramsr = zeros(numdisks)150for i in range(numdisks):151 paramsr[i] = 3152d = 5 # distance between disks153# Disk Up #154x1 = -d/6.0 * sqrt(3.0)155y1 = d/2.0156# Disk Right #157x2 = d/3.0 * sqrt(3.0)158y2 = 0159# Disk Left # 160x3 = -d/6.0 * sqrt(3.0)161y3 = -d/2.0162xcoordinates = zeros(numdisks)163ycoordinates = zeros(numdisks)164xcoordinates[0] = x1165xcoordinates[1] = x2166xcoordinates[2] = x3167ycoordinates[0] = y1168ycoordinates[1] = y2169ycoordinates[2] = y3170#####################171####Shoot Particle#####172#####################173# Parameters of Problem #174xp = -4 # starting point of particle175ystart = -0.5 # beginning impact parameter176yend = 0.5 # ending impact parameter177interval = 100000 # points between ystart and yend and how many times we will run the code178impact = linspace(ystart,yend,interval)179theta = 0180# WHAT ARE WE GRAPHING? #181ExitAngle = zeros(interval) # How many exit angles we have is how many times we run the code182BounceCount = zeros(interval) # How Many Times We Bounce For Each Trajectory183TimeSpent = zeros(interval) # Time Spent For Each Trajectory184NumParticles = linspace(0,interval,interval)185##### STARTING ALGORITHM #####186IntersectPointsHeld = zeros([16,2]) # We Use This For Plotting A Single Trajectory187SimulationCounter = 1 # Simulation Counter188Counter = 0 # Initialze Counter For Things We Will Be Graphing189PlotNow = False # This is only used for Plotting A Single Trajectory190for i in impact: 191 exitParticle = False192 193 # Every time we run the code we need to initialize the beginning values #194 LastImpactX = xp195 LastTheta = theta196 LastImpactY = i197 198 199 # if LastImpactY == impact[91]: # scatter number200 #PlotNow = True201 #IntersectPointsHeld[0] = [LastImpactX,LastImpactY]202 #else:203 #PlotNow = False204 #IntersectPointsHeld[15] = IntersectPointsHeld[14]205 206 IntersectArray = Intersect(xcoordinates,ycoordinates,paramsr,LastTheta,LastImpactX,LastImpactY)207 IntersectTime = IntersectArray[0]208 IntersectCircle = IntersectArray[1]209 IntersectPoints = IntersectArray[2]210 211 if IntersectCircle == 0:212 exitParticle = True213 print("We Have Exited")214 timespentthistrajectory = IntersectTime215 216 print("Running Simulation %s" % SimulationCounter)217 218 BounceCounter = 1219 220 timespentthistrajectory = 0221 222 while exitParticle == False:223 224 225 226 print("Projectile has bounced %s times" % BounceCounter) # Sanity Check #227 228 # Use Intersect Function #229 230 IntersectArray = Intersect(xcoordinates,ycoordinates,paramsr,LastTheta,LastImpactX,LastImpactY) 231 232 # Unpack Intersect Function #233 234 IntersectTime = IntersectArray[0]235 IntersectCircle = IntersectArray[1]236 IntersectPoints = IntersectArray[2]237 238 239 240 if PlotNow == True: # IF WE WANT TO PLOT A TRAJECTORY WE USE THIS LINE OF CODE #241 IntersectPointsHeld[BounceCounter] = IntersectPoints242 243 244 if IntersectCircle == 0: # This is our exit condition245 ExitAngle[Counter] = LastTheta246 exitParticle = True247 print("We Have Exited")248 249 250 251 NewAngle = flipvectorgettheta(IntersectPoints,LastImpactX,LastImpactY,IntersectCircle)252 253 LastImpactX = IntersectPoints[0]254 LastImpactY = IntersectPoints[1]255 LastTheta = NewAngle256 257 if exitParticle == False:258 timespentthistrajectory += IntersectTime259 260 261 262 BounceCounter += 1263 264 265 TimeSpent[Counter] = timespentthistrajectory266 267 268 BounceCount[Counter] = BounceCounter269 270 Counter += 1271 SimulationCounter += 1272 273#####################274######Plots##########275#####################276ax=plt.gca()277plt.plot(IntersectPointsHeld[:,0],IntersectPointsHeld[:,1],)278circle1 = plt.Circle((x1,ycoordinates[0]), radius=1, color='g', fill=False)279circle2 = plt.Circle((x2,ycoordinates[1]), radius=1, color='b', fill=False)280circle3 = plt.Circle((x3,ycoordinates[2]), radius=1, color='r', fill=False)281ax.add_patch(circle1)282ax.add_patch(circle2)283ax.add_patch(circle3)284plt.axis('scaled')285plt.title('Trajectory Plot')286plt.show()287plt.figure()288plt.plot(impact,ExitAngle)289plt.xlabel('Impact Parameter (''$\it{b}$)')290plt.ylabel('Exit Angle ('r'$\theta$)')291plt.title('Output ('r'$\theta$) vs Input Ordinate')292plt.show()293plt.figure()294plt.plot(impact,BounceCount)295plt.title("Bounce Count vs Input Ordinate")296plt.xlabel('Impact Parameter (''$\it{b}$)')297plt.ylabel('Number Bounces')298plt.show()299plt.figure()300plt.plot(impact,TimeSpent)301plt.xlabel('Impact Parameter (''$\it{b}$)')302plt.ylabel('Time Spent')303plt.title('Time Spent vs Input Ordinate')...

Full Screen

Full Screen

Automation Testing Tutorials

Learn to execute automation testing from scratch with LambdaTest Learning Hub. Right from setting up the prerequisites to run your first automation test, to following best practices and diving deeper into advanced test scenarios. LambdaTest Learning Hubs compile a list of step-by-step guides to help you be proficient with different test automation frameworks i.e. Selenium, Cypress, TestNG etc.

LambdaTest Learning Hubs:

YouTube

You could also refer to video tutorials over LambdaTest YouTube channel to get step by step demonstration from industry experts.

Run fMBT automation tests on LambdaTest cloud grid

Perform automation testing on 3000+ real desktop and mobile devices online.

Try LambdaTest Now !!

Get 100 minutes of automation test minutes FREE!!

Next-Gen App & Browser Testing Cloud

Was this article helpful?

Helpful

NotHelpful