Unlock 30% off on Manual Testing Annual Plans this Holiday Season.

Copied to Clipboard!

How to use intersect method in fMBT

Best Python code snippet using fMBT_python

test_qp_subproblem.py

Source:test_qp_subproblem.py

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)...

ResourceAreaCharacterization.py

Source:ResourceAreaCharacterization.py

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__":
...

Fractal Images & Chaotic Scattering.py

Source:Fractal Images & Chaotic Scattering.py

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')...

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.