From 1bde91a4a0e683aef79be792831aeacae4699534 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20St=C3=B6lzle?= Date: Mon, 26 May 2025 17:28:36 +0200 Subject: [PATCH 1/5] Rename `pneumatic_planar_pcs` to `pneumatically_actuated_planar_pcs` --- ...r_pcs.py => simulate_pneumatically_actuated_planar_pcs.py} | 4 ++-- src/jsrm/systems/planar_pcs.py | 4 +++- ...tic_planar_pcs.py => pneumatically_actuated_planar_pcs.py} | 4 +++- 3 files changed, 8 insertions(+), 4 deletions(-) rename examples/{simulate_pneumatic_planar_pcs.py => simulate_pneumatically_actuated_planar_pcs.py} (98%) rename src/jsrm/systems/{pneumatic_planar_pcs.py => pneumatically_actuated_planar_pcs.py} (98%) diff --git a/examples/simulate_pneumatic_planar_pcs.py b/examples/simulate_pneumatically_actuated_planar_pcs.py similarity index 98% rename from examples/simulate_pneumatic_planar_pcs.py rename to examples/simulate_pneumatically_actuated_planar_pcs.py index 1aa2a58..cd2634f 100644 --- a/examples/simulate_pneumatic_planar_pcs.py +++ b/examples/simulate_pneumatically_actuated_planar_pcs.py @@ -12,7 +12,7 @@ import jsrm from jsrm import ode_factory -from jsrm.systems import pneumatic_planar_pcs +from jsrm.systems import pneumatically_actuated_planar_pcs num_segments = 1 @@ -48,7 +48,7 @@ strain_selector = jnp.array([True, False, True])[None, :].repeat(num_segments, axis=0).flatten() B_xi, forward_kinematics_fn, dynamical_matrices_fn, auxiliary_fns = ( - pneumatic_planar_pcs.factory( + pneumatically_actuated_planar_pcs.factory( num_segments, sym_exp_filepath, strain_selector, # simplified_actuation_mapping=True ) ) diff --git a/src/jsrm/systems/planar_pcs.py b/src/jsrm/systems/planar_pcs.py index 3465335..fbfebcb 100644 --- a/src/jsrm/systems/planar_pcs.py +++ b/src/jsrm/systems/planar_pcs.py @@ -242,6 +242,7 @@ def actuation_mapping_fn( jacobian_fn: Callable, params: Dict[str, Array], B_xi: Array, + xi_eq: Array, q: Array, ) -> Array: """ @@ -252,6 +253,7 @@ def actuation_mapping_fn( jacobian_fn: function to compute the Jacobian params: dictionary with robot parameters B_xi: strain basis matrix + xi_eq: equilibrium strains as array of shape (n_xi,) q: configuration of the robot Returns: A: actuation matrix of shape (n_xi, n_xi) where n_xi is the number of strains. @@ -359,7 +361,7 @@ def dynamical_matrices_fn( # compute the stiffness matrix K = stiffness_fn(params, B_xi, formulate_in_strain_space=True) # compute the actuation matrix - A = actuation_mapping_fn(forward_kinematics_fn, jacobian_fn, params, B_xi, q) + A = actuation_mapping_fn(forward_kinematics_fn, jacobian_fn, params, B_xi, xi_eq, q) # dissipative matrix from the parameters D = params.get("D", jnp.zeros((n_xi, n_xi))) diff --git a/src/jsrm/systems/pneumatic_planar_pcs.py b/src/jsrm/systems/pneumatically_actuated_planar_pcs.py similarity index 98% rename from src/jsrm/systems/pneumatic_planar_pcs.py rename to src/jsrm/systems/pneumatically_actuated_planar_pcs.py index 8766b5b..0c91d5b 100644 --- a/src/jsrm/systems/pneumatic_planar_pcs.py +++ b/src/jsrm/systems/pneumatically_actuated_planar_pcs.py @@ -15,7 +15,7 @@ def factory( **kwargs ): """ - Factory function for the planar PCS. + Factory function for the pneumatically-actuated planar PCS. Args: num_segments: number of segments segment_actuation_selector: actuation selector for the segments as boolean array of shape (num_segments,) @@ -43,6 +43,7 @@ def actuation_mapping_fn( jacobian_fn: Callable, params: Dict[str, Array], B_xi: Array, + xi_eq: Array, q: Array, ) -> Array: """ @@ -52,6 +53,7 @@ def actuation_mapping_fn( jacobian_fn: function to compute the Jacobian params: dictionary with robot parameters B_xi: strain basis matrix + xi_eq: equilibrium strains as array of shape (n_xi,) q: configuration of the robot Returns: A: actuation matrix of shape (n_xi, n_act) where n_xi is the number of strains and From 1b3d3a9b4f72b66c6779eeecc7c9d36937cf2853 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20St=C3=B6lzle?= Date: Mon, 26 May 2025 17:29:39 +0200 Subject: [PATCH 2/5] Start working on `tendon_acutated_planar_pcs` --- examples/derive_planar_pcs.py | 2 +- .../simulate_tendon_actuated_planar_pcs.py | 247 ++++++++++++++++++ src/jsrm/symbolic_derivation/planar_pcs.py | 31 ++- .../symbolic_expressions/planar_pcs_ns-1.dill | Bin 8068 -> 7890 bytes .../symbolic_expressions/planar_pcs_ns-2.dill | Bin 75819 -> 70480 bytes .../systems/tendon_actuated_planar_pcs.py | 136 ++++++++++ 6 files changed, 410 insertions(+), 6 deletions(-) create mode 100644 examples/simulate_tendon_actuated_planar_pcs.py create mode 100644 src/jsrm/systems/tendon_actuated_planar_pcs.py diff --git a/examples/derive_planar_pcs.py b/examples/derive_planar_pcs.py index 8923f51..23f1664 100644 --- a/examples/derive_planar_pcs.py +++ b/examples/derive_planar_pcs.py @@ -3,7 +3,7 @@ import jsrm from jsrm.symbolic_derivation.planar_pcs import symbolically_derive_planar_pcs_model -NUM_SEGMENTS = 1 +NUM_SEGMENTS = 2 if __name__ == "__main__": sym_exp_filepath = ( diff --git a/examples/simulate_tendon_actuated_planar_pcs.py b/examples/simulate_tendon_actuated_planar_pcs.py new file mode 100644 index 0000000..55025e6 --- /dev/null +++ b/examples/simulate_tendon_actuated_planar_pcs.py @@ -0,0 +1,247 @@ +import cv2 # importing cv2 +from functools import partial +import jax + +jax.config.update("jax_enable_x64", True) # double precision +from diffrax import diffeqsolve, Euler, ODETerm, SaveAt, Tsit5 +from jax import Array, vmap +from jax import numpy as jnp +import matplotlib.pyplot as plt +import numpy as onp +from pathlib import Path +from typing import Callable, Dict + +import jsrm +from jsrm import ode_factory +from jsrm.systems import tendon_actuated_planar_pcs as planar_pcs + +num_segments = 1 + +# filepath to symbolic expressions +sym_exp_filepath = ( + Path(jsrm.__file__).parent + / "symbolic_expressions" + / f"planar_pcs_ns-{num_segments}.dill" +) + +# set parameters +rho = 1070 * jnp.ones((num_segments,)) # Volumetric density of Dragon Skin 20 [kg/m^3] +params = { + "th0": jnp.array(0.0), # initial orientation angle [rad] + "l": 1e-1 * jnp.ones((num_segments,)), + "r": 2e-2 * jnp.ones((num_segments,)), + "rho": rho, + "g": jnp.array([0.0, 9.81]), + "E": 2e3 * jnp.ones((num_segments,)), # Elastic modulus [Pa] + "G": 1e3 * jnp.ones((num_segments,)), # Shear modulus [Pa] + "d": 2e-2 * jnp.array([[1.0, -1.0]]).repeat(num_segments, axis=0), # distance of tendons from the central axis [m] +} +print("params d =\n", params["d"]) +params["D"] = 1e-3 * jnp.diag( + (jnp.repeat( + jnp.array([[1e0, 1e3, 1e3]]), num_segments, axis=0 + ) * params["l"][:, None]).flatten() +) + +# activate all strains (i.e. bending, shear, and axial) +strain_selector = jnp.ones((3 * num_segments,), dtype=bool) + +# define initial configuration +q0 = jnp.repeat(jnp.array([5.0 * jnp.pi, 0.1, 0.2])[None, :], num_segments, axis=0).flatten() +# number of generalized coordinates +n_q = q0.shape[0] + +# set simulation parameters +dt = 1e-4 # time step +ts = jnp.arange(0.0, 2, dt) # time steps +skip_step = 10 # how many time steps to skip in between video frames +video_ts = ts[::skip_step] # time steps for video + +# video settings +video_width, video_height = 700, 700 # img height and width +video_path = Path(__file__).parent / "videos" / f"planar_pcs_ns-{num_segments}.mp4" + + +def draw_robot( + batched_forward_kinematics_fn: Callable, + params: Dict[str, Array], + q: Array, + width: int, + height: int, + num_points: int = 50, +) -> onp.ndarray: + # plotting in OpenCV + h, w = height, width # img height and width + ppm = h / (2.0 * jnp.sum(params["l"])) # pixel per meter + base_color = (0, 0, 0) # black robot_color in BGR + robot_color = (255, 0, 0) # black robot_color in BGR + + # we use for plotting N points along the length of the robot + s_ps = jnp.linspace(0, jnp.sum(params["l"]), num_points) + + # poses along the robot of shape (3, N) + chi_ps = batched_forward_kinematics_fn(params, q, s_ps) + + img = 255 * onp.ones((w, h, 3), dtype=jnp.uint8) # initialize background to white + curve_origin = onp.array( + [w // 2, 0.1 * h], dtype=onp.int32 + ) # in x-y pixel coordinates + # draw base + cv2.rectangle(img, (0, h - curve_origin[1]), (w, h), color=base_color, thickness=-1) + # transform robot poses to pixel coordinates + # should be of shape (N, 2) + curve = onp.array((curve_origin + chi_ps[:2, :].T * ppm), dtype=onp.int32) + # invert the v pixel coordinate + curve[:, 1] = h - curve[:, 1] + cv2.polylines(img, [curve], isClosed=False, color=robot_color, thickness=10) + + return img + + +if __name__ == "__main__": + strain_basis, forward_kinematics_fn, dynamical_matrices_fn, auxiliary_fns = ( + planar_pcs.factory(sym_exp_filepath, strain_selector) + ) + # jit the functions + dynamical_matrices_fn = jax.jit(partial(dynamical_matrices_fn)) + batched_forward_kinematics = vmap( + forward_kinematics_fn, in_axes=(None, None, 0), out_axes=-1 + ) + + # import matplotlib.pyplot as plt + # plt.plot(chi_ps[0, :], chi_ps[1, :]) + # plt.axis("equal") + # plt.grid(True) + # plt.xlabel("x [m]") + # plt.ylabel("y [m]") + # plt.show() + + # Displaying the image + # window_name = f"Planar PCS with {num_segments} segments" + # img = draw_robot(batched_forward_kinematics, params, q0, video_width, video_height) + # cv2.namedWindow(window_name) + # cv2.imshow(window_name, img) + # cv2.waitKey() + # cv2.destroyWindow(window_name) + + x0 = jnp.concatenate([q0, jnp.zeros_like(q0)]) # initial condition + tau = jnp.zeros_like(q0) # torques + + ode_fn = ode_factory(dynamical_matrices_fn, params, tau) + term = ODETerm(ode_fn) + + sol = diffeqsolve( + term, + solver=Tsit5(), + t0=ts[0], + t1=ts[-1], + dt0=dt, + y0=x0, + max_steps=None, + saveat=SaveAt(ts=video_ts), + ) + + print("sol.ys =\n", sol.ys) + # the evolution of the generalized coordinates + q_ts = sol.ys[:, :n_q] + # the evolution of the generalized velocities + q_d_ts = sol.ys[:, n_q:] + + # evaluate the forward kinematics along the trajectory + chi_ee_ts = vmap(forward_kinematics_fn, in_axes=(None, 0, None))( + params, q_ts, jnp.array([jnp.sum(params["l"])]) + ) + # plot the configuration vs time + plt.figure() + for segment_idx in range(num_segments): + plt.plot( + video_ts, q_ts[:, 3 * segment_idx + 0], + label=r"$\kappa_\mathrm{be," + str(segment_idx + 1) + "}$ [rad/m]" + ) + plt.plot( + video_ts, q_ts[:, 3 * segment_idx + 1], + label=r"$\sigma_\mathrm{sh," + str(segment_idx + 1) + "}$ [-]" + ) + plt.plot( + video_ts, q_ts[:, 3 * segment_idx + 2], + label=r"$\sigma_\mathrm{ax," + str(segment_idx + 1) + "}$ [-]" + ) + plt.xlabel("Time [s]") + plt.ylabel("Configuration") + plt.legend() + plt.grid(True) + plt.tight_layout() + plt.show() + # plot end-effector position vs time + plt.figure() + plt.plot(video_ts, chi_ee_ts[:, 0], label="x") + plt.plot(video_ts, chi_ee_ts[:, 1], label="y") + plt.xlabel("Time [s]") + plt.ylabel("End-effector Position [m]") + plt.legend() + plt.grid(True) + plt.box(True) + plt.tight_layout() + plt.show() + # plot the end-effector position in the x-y plane as a scatter plot with the time as the color + plt.figure() + plt.scatter(chi_ee_ts[:, 0], chi_ee_ts[:, 1], c=video_ts, cmap="viridis") + plt.axis("equal") + plt.grid(True) + plt.xlabel("End-effector x [m]") + plt.ylabel("End-effector y [m]") + plt.colorbar(label="Time [s]") + plt.tight_layout() + plt.show() + # plt.figure() + # plt.plot(chi_ee_ts[:, 0], chi_ee_ts[:, 1]) + # plt.axis("equal") + # plt.grid(True) + # plt.xlabel("End-effector x [m]") + # plt.ylabel("End-effector y [m]") + # plt.tight_layout() + # plt.show() + + # plot the energy along the trajectory + kinetic_energy_fn_vmapped = vmap( + partial(auxiliary_fns["kinetic_energy_fn"], params) + ) + potential_energy_fn_vmapped = vmap( + partial(auxiliary_fns["potential_energy_fn"], params) + ) + U_ts = potential_energy_fn_vmapped(q_ts) + T_ts = kinetic_energy_fn_vmapped(q_ts, q_d_ts) + plt.figure() + plt.plot(video_ts, U_ts, label="Potential energy") + plt.plot(video_ts, T_ts, label="Kinetic energy") + plt.xlabel("Time [s]") + plt.ylabel("Energy [J]") + plt.legend() + plt.grid(True) + plt.box(True) + plt.tight_layout() + plt.show() + + # create video + fourcc = cv2.VideoWriter_fourcc(*"mp4v") + video_path.parent.mkdir(parents=True, exist_ok=True) + video = cv2.VideoWriter( + str(video_path), + fourcc, + 1 / (skip_step * dt), # fps + (video_width, video_height), + ) + + for time_idx, t in enumerate(video_ts): + x = sol.ys[time_idx] + img = draw_robot( + batched_forward_kinematics, + params, + x[: (x.shape[0] // 2)], + video_width, + video_height, + ) + video.write(img) + + video.release() + print(f"Video saved at {video_path}") diff --git a/src/jsrm/symbolic_derivation/planar_pcs.py b/src/jsrm/symbolic_derivation/planar_pcs.py index 338a876..ed9977e 100644 --- a/src/jsrm/symbolic_derivation/planar_pcs.py +++ b/src/jsrm/symbolic_derivation/planar_pcs.py @@ -40,6 +40,7 @@ def symbolically_derive_planar_pcs_model( sp.symbols(f"r1:{num_segments + 1}", nonnegative=True) ) # radius of each segment [m] g_syms = list(sp.symbols(f"g1:3")) # gravity vector + d = sp.symbols("d", real=True, nonnegative=True) # # distance of the tendon from the neutral axis # planar strains and their derivatives xi_syms = list(sp.symbols(f"xi1:{num_dof + 1}", nonzero=True)) # strains @@ -59,6 +60,10 @@ def symbolically_derive_planar_pcs_model( chi_sms = [] # Jacobians (positional + orientation) in each segment as a function of the point coordinate s and its time derivative J_sms, J_d_sms = [], [] + # tendon lengths for each segment as a function of the point coordinate s + L_tend_sms = [] + # tendon length jacobians for each segment as a function of the point coordinate s + J_tend_sms = [] # cross-sectional area of each segment A = sp.zeros(num_segments) # second area moment of inertia of each segment @@ -74,13 +79,14 @@ def symbolically_derive_planar_pcs_model( # initialize th_prev = th0 p_prev = sp.Matrix([0, 0]) + L_tend = 0 # tendon length for i in range(num_segments): # bending strain - kappa = xi[3 * i] + kappa_be = xi[3 * i] # shear strain - sigma_x = xi[3 * i + 1] + sigma_sh = xi[3 * i + 1] # axial strain - sigma_y = xi[3 * i + 2] + sigma_ax = xi[3 * i + 2] # compute the cross-sectional area of the rod A[i] = sp.pi * r[i] ** 2 @@ -89,13 +95,13 @@ def symbolically_derive_planar_pcs_model( I[i] = A[i] ** 2 / (4 * sp.pi) # planar orientation of robot as a function of the point s - th = th_prev + s * kappa + th = th_prev + s * kappa_be # absolute rotation of link R = sp.Matrix([[sp.cos(th), -sp.sin(th)], [sp.sin(th), sp.cos(th)]]) # derivative of Cartesian position as function of the point s - dp_ds = R @ sp.Matrix([sigma_x, sigma_y]) + dp_ds = R @ sp.Matrix([sigma_sh, sigma_ax]) # position along the current rod as a function of the point s p = p_prev + sp.integrate(dp_ds, (s, 0.0, s)) @@ -146,12 +152,24 @@ def symbolically_derive_planar_pcs_model( # add potential energy of segment to previous segments U_g = U_g + U_gi + # simplify derived tendon length + L_tend = L_tend + s * (1 + kappa_be * d) * sp.sqrt(sigma_sh**2 + sigma_ax**2) + L_tend_sms.append(L_tend) + print(f"L_tend of segment {i+1}:\n", L_tend) + # take the derivative of the tendon length with respect to the configuration + J_tend = sp.simplify(sp.Matrix([L_tend]).jacobian(xi)) + J_tend_sms.append(J_tend) + print(f"J_tend of segment {i+1}:\n", J_tend) + # update the orientation for the next segment th_prev = th.subs(s, l[i]) # update the position for the next segment p_prev = p.subs(s, l[i]) + # previous tendon length + L_tend = L_tend.subs(s, l[i]) + if simplify_expressions: # simplify mass matrix B = sp.simplify(B) @@ -174,6 +192,7 @@ def symbolically_derive_planar_pcs_model( "r": r_syms, "rho": rho_syms, "g": g_syms, + "d": d, }, "state_syms": { "xi": xi_syms, @@ -197,6 +216,8 @@ def symbolically_derive_planar_pcs_model( "C": C, # coriolis matrix "G": G, # gravity vector "U_g": U_g, # gravitational potential energy + "L_tend_sms": L_tend_sms, # list of tendon lengths for each segment + "J_tend_sms": J_tend_sms, # list of tendon length Jacobians for each segment }, } diff --git a/src/jsrm/symbolic_expressions/planar_pcs_ns-1.dill b/src/jsrm/symbolic_expressions/planar_pcs_ns-1.dill index 9118c75ec1f8f7e596640975ac0089007aa7e4c1..7cf186f98f1d6cc2b75a76463de2516dc3d88393 100644 GIT binary patch literal 7890 zcmcIpO>A99754ja(z>N;UqM7diij+nM}(}zX$&%;fN4tnZd0PPl@OA#pWnNA^M?1Q z_Y=QViR3IuvX}*`Zzb4NC5s5LKw{M$LRO^6REPx<5(`!cL04?p@KeJ1&dj-Y?)$Ob zO1Qc&cV^Ck}t*}<^*KdT; zeRG_eZ-w19I_WMv5>}!@wQ3S6kxO&}6FA zR;Ssxdip*oBtru;*Xx)2VG>?;Sa(`-!@AS|0&J6B>VdY~k@9-Zdp)7a59=FMr(hQE zu_xTJ&X6i@Vd|ZTTyi1iMZQ_?ck7j~w^$8Zy)ZI` zD}#P{s}a78<`oPMquuC>cNpt@x4qqqf}F`$+6{WinEXaJ>_q02Rn}=YZuS-vmF;G^ z-s0vkGS9xu%~&TbUNCbz3TDi6b~clsT&+eX_f7iUrN#9}ngilS?h5@9F_z8q-lR@@ zJM7Zn)%G?9mrTJMYz>-QVYe5VXI^tO`A#deG2JEDOwr0-AG9j{db`zI3>#rHZ1v0C zn~M~Ft=(!jLx4mEdi53ztgAV-sg~3l6yO2=_voSC6;)T-eS+=IE{w{h) zMyx&7pw5NDx!Fm7*H%nl& zV6iSY0w?r~i;T8NQc>DeHozLq__Le6Mpcz{isc}ilKarg zA$EweVi6snJWs*ypDv1o{o=c;W6nv4B_@_*-(x9;*T+t8Ak>dANgN@MEn0u2H;$iG za$HuqB9!zGXd;?;(s1aU8>0CfCPYMgkkz+oc!2#QcQ=w!NyVSCj*cfEmskb;SyC6{ z!|b=*DyRoOw+f2ucLEXRCyCgIWCn|j;iEt20+JwNh)@tD(U5w!9E=m{tE4~(C4L=1 zD1FKqBtr!v3ab_4X-?_Wu_%IoG>0wH(`i-E=V>bBc%7wgNDz)Q&RdA3p_&BlIC5Wh z-ibScoc;n0IJ)XBmvD$zHA7jX|AswdLP_ZK9dJFH28K0c@l{2BQBcMaAP+lQ2#@0n9-i-R^hlONbqQgdE6k4?PTQI#@WJTg9 zD!KEMNL1tK9T)@=-th2roOkq(5EsAwnBLG26int@sYLOXjRS^`#7fo7zS3xy`#5&e zms{o>eP9HcK7V$ar|Bc0QpZ<76kNbpDs7m6;-_Ser7uuo2RSRPA5t%iS<-pfIC=du zNKBf5#PAlo=q$|W2Q*OHM`}Bt5-T>bgW0jHP>%#A0BmhrO+j{o8iJM`27X7L0Eb5= zm`|0xzpbeB0o;CR6dZP&J+v_W>Jw|yzXxgz-2^2JU9y%5ZZ`Q*dN}zR_K{Eynad_o zQzp$dne!k;b}VkenNxVw;D@Ab>MXX0dVyJZKyz=Y%yZ{1_L8Y(%mea~rhVkPmwj6BjF zF*zyi(zWw|dL&zn6it-CiG|DGY~mGQX}Z`Mc(Q)2g|e0g<&c6v^+ ziw>ezT_uKEYpBV?aPYyr_ZbmV+mnU+t7 zQFX-N^d$0}oSqSr*22#|oD$O(5A=@;Oq@zyl{qZ?R-F}Cl5HYw^}xO25bEpq9U-Ao z|B@PF$7H^3&*D{+i~s+EdJNB9Qg~;6I$0O#(`cSH{P?08O}bNczS8Idt}p%ON-S)5 zGfdxc0@6eGn5iBp5uH|Mfvm_ZvcY9;yd!q2Ms8UdQ9}m;=@MJ63zsypYAb~6DIfOMvWMwKcA_ngmy||fh%S}qHqGOd{inBFL zo4Hpv`e7@+i^er$h0bIARf6L_!$nv8m&=Xo_9lmW)-5JD{GRb~EStGfN})v%*Z&mp s8MUV0B`j^wMT&r|>7M{%ToX*y46hM2U195=L;4r=BXB@jc`#V~AD>-on*aa+ literal 8068 zcmb_hU5H#)9pCxb#3Z$(qe4Jvk-&C}IBmAgD5gY{W+R`%huW1*_%Ccm@nqe zxU-UgK^phrD71SELPZd%dB{V3C_*392VZm>D2s^CdGJX@M4_P7`2YRS$31uE&hDnT z%gmmS|L5=jan8N_-Qt^@kIeDU^|-8#^&5jmC#pxAohZgaVJM!6)k8GzZ&sSUL0F;T zYOfuuLoe#Vb}S075v?i?!bUsZir$D-t{uM`i$f~kJ{4n+1b2JgZn)kUw!R(4Tk0@P zzY`96@z&^#xS?_bLqE{;;yA=KRAC@`dZ$RE^j>bAcFIt)j*T0vZ+(%6;Yb~ghK*sE zIFcW?41HnTGLjyG7Lq_K;JY3RZ{kh<6urbsJS`ZjkUE{70%K3HigmO>5M1~3YX}cQeaUX_n7*6kERV%bLHT7sdA!}m;VZ>7<{y$||tTlB^220HPk_!0R6RDMdRyi2J33aF6VzjofS z`@CX+Q1dlgxdwCxcdOIUHA_YPtR1K8SFEVR&vu95I+e*b6)+GGchqORGHLcGKR46b=<)NdjYD&yO z1su&@$h1wN<#(|aV$B?y*@on^-1D!kq1h}(+N}I@NTT1Ey*1;xm{KarD zQP*dgRx|YppYHFP4!(nw1zXY)j*o2YWRC1HjwA(O%{Lv;*a@;k%c3*iQ8B_kd~%(n z%X#_&pEZZ<)|mZ`MvMYWW$(0nL0HMiuwf|;iachvh@YMHr3ZQW8f}ko$`x+L8GTkh zPy0bm{t#r`Y2|w~G|jf4&l>TTUQ7-iDqXsR-EVf5qmUGpDBD=94R1Oz z<>QVNPJo)-w$!t8vn0VM#*#42v2Qc@{R-Mb+)g^R%Kj5feT>aOCX-C%`7ZYUM1;8V%*ei8=P^BJAsM zbn-H7nqcsSBq&;Y<=K5adlols{-R+9QNp7IFQ(b3 zeOP<#Uu%0C9%_SSA7jRhpDuQVQ6YzLXkRVV(BK(kYdo{oOxowWr89!6c@kR8S@2?6`BupTI93{hUTv2MQ|pOsY`4u(89CVbCBgD{B5? zyVn@v>Pt^0xl{BI7UYKfXs(XXgHp4FC#5(zjVF1!l>mx&NF7t3roIpqG_`z}R+&sE zorjH~m%ju>76S^y3wGIPn3KPtjoR-~+Qz0ti!SV7eyS?O9Saix*1EN%AfF+Ipk;%B z-+>um@C1YUWW)Pejcf12?58HdVYAsn4bz>EElJ)5Y7AY55QZ*^WrCY7ezYDgeujO* zm3??wC3;F0ugX1%q{vVCZE@xjo>cfDNt?;V`p_@Hg%8r)3zd5C)M@`GuI&K&(-|h< z7Q0jXr|C))6-Ia5o!}7iiE_MDApf8g$wYH*nQ$6&7=RD&LFQT8Q|RF29${0^{{0x{ z>cTxk7lMc7mq@WGT6vr-BbD4W+06T~b^`c=LT7EM`nJd~k|kE4>9wArWFyW*4G(?W zwjnv#bJHo}g(F?Ix$Lr?NQ0SK3#eSyxv*J#Op~FM$Z9g@=z#{OEIG65rGrMF7t!Z9 z4>%+k;bSwzs?)Zn-ExgKDb93<-J@&SgUGg4G|J5O<{wEe0!&S>$fhd9Kc**oXO^|q zq%3S69H8HtGchTN+(FZjt`Xey?x)`okE`9RhBZ5FJnvlkBu4KMoC+@S87w)iEBh3p z%ZxM6HaV|*(0zbQ>;GLZ>DX#7^u=&0gR+f81<7UCx8y1fZG~9j znGHm@kES>$uNW5BNA|AUv+}B8uodYTO%$^7Id#9;Bf$0onM(GodP-d`_pMd`dV7HC z^Ag*FXbd`VJymDb(|)Nsp>TE)c}C8zh^$!n)rTIN*0`sAsJ{DTLG|IZhZcNn-%e}Y zGw1e6efhpgNF#Eef_{&vamwg}qbz@8=FUoN_|~f2*Yv%AMHSrtd&CE5Sf;eR=L@>U zm+3L=BsD7cJ|9zhhiuK;6$`ZM-Z5ZWb;SBi1Uvud0803jgapf?S0Kd7v@7)b+-g}M zA1PeZ3(l-2nWbs07`j&2C9BIU;eC=5t*Nl&$$~yCd9v&h%d&_ju+e+EE;N!}(e2!8iuFKG1lr1xHaK~vfr!oJ$rZ{v%XFlBcCqfZb@H$vP^hQX`&O6v{R0lF(-%{Q+Aj!Uri~+h?s)!F_U|{(tUI9IeUG7?vH!!PaOK$gSV{Qf`4DwJGOYo<hPhx#a;CI@{{+RxpZaw9(w!8rSp4>TOW`&&+ZKepQ5pgLszy>pWoZv z`SRXk<^0}f_Kpp=E)Jak@E&dw=8KmuUfh1{^wo1ezP-1*xQ#x3e*4O$z1?eH-s6GC zW#I7c-Qg}=++M6)Q7l&k3*95^7=1-X%Z2U5!7Ib1IQ9^|9gW7Na{S6;iuo~&z0Nov z!gh+$Z~p*3 zei}EQ*n4*G&tGRhICSOG<2!pB2N#FVTslv0Ru+dIy0U$FZ}9;^cKOozCwK12C@)<& zeeNQEd}43$&L7~92`4$$2KnvAEsO6GX!g+4XV2~}4t|3EK1)p6I~T5z;~xAp{UxtG zxp=n$x_s&J?JIQiM=w2&nj%a&K$z+9Si~#3P2Q8wZDL8!N+)4^Jk0X0E@KG_~5$)N{kX&Kepv zha1ECu#2&zCTRI3l7C(F45`i7t>IPD+BfnQ)dwqJxUQOgA**MTKwi!tf#H>`t)%gP z9Ig0r)H%RjiSGgPTB(0Vk?9%o0DQfq$K-KoXG5&>4Ykhji=x`?#bFB5?d`p#I{ua6 zSBX;)vAQ|}?YX-8gSSOjzcc(EaY1?sC`B_=6Asb%CECd@ki;76YjJoz<+Se>2Tuv; z+=UC0&Gr`W{gl2hGu!ag;TuEcaf9H z;lFb2@_D2YMCgs-1jDTV!JBWs2{88H7`)H(2PB3U;ydJl7{na@Uy05|35A><0soaW za!p3!Ge$vQ@!m5$c6Si`&*JLnj?iOVeCMgfp$Ad6IiH6INl5-PTn`_k_ZqfJdo5zV z&Us4{5Ym{^Y1MD9B`4IlxM}No= zBFC$&n01xM*S`ig-8fKMb8~n->Ue8-IlkMa5Y+UrOM$@YK^a9}I6)a~ zH5&Gn1hh6V3;3*qqZ^clWik;;?}%X;SG+@cGu%%~;E(Khn~A}#AyeJd>_tZ18eWeS zDZ>t549^fe*)YifTq7&K8huV=e+?rd1A0eDv6keHYZGOt?PIJQYm=}G_M&J^nRq1; zB7ImP+mu9B&PBx13Tiqsy^7GVfiRlg_2*ef;cJE$@#Yv#ck^koz;>#YBFv$jU0979POv@*HnlQ zelRO2mH_F!!tgAJXfG57ZcSAJl5LIhOwql@Az3(1Vk!Nu{~uFN!>6K#fiT+V!N?a| zM~e`6l9Kwo8iekDCJqgk*#zJi1BGL*ZaSFS+3N#KJ?t`2_>oANt_Nzy8M$(?abP_9 z$I@sJ0RR(G0BV=HRc}n(17zU!q}=eg;9BP+kIEFu zpt=kX&^`Lh+F|r8uMmoA3_nFk+WNYZG$7(*q*B^A5S}OCYnYVW`t+W*y~`XCZIK5 z>jVYem+6g;O(kFXG1Crk%X0>>2LOu=HJmL({q)a8ghq^qb1Xsw;qw^&H|Y5#a1}Y>BA|Yr(Lc^N|1M-=}+c9+@Nl#>%@G=iY~`B;%N{_;DqIc!H8ECpa}YA>l&vA}ftbw* z4Q5-hEVvX$<>dC;gkHy{*wU9>yr}h5s!oo7w(*#Aehs&TQYHA>lRpL-HZB z+4-ktIb=3RZh}Bsi(-Q@ju~YPvKhvTv{o}#r!{7%nt|yABZWUmK(m_L!UxFT^`ho0 z&XAmz^bY5%@ZmTXSL+n@O;Ss=DQ6tKS5;DZ{Uia?vK$x4I`P~Ra{~KKpq%+z9EYx& z#4}no0?dgAa>fY{5u$L2&|7|pOqFTEFBy_C-2f5b!tW4*@OvrhmzUrnWfac-l2Ifl zu?s{^Nk$1@Akw(Pnsrjgosoayo1=RyOG+I&1x0kTL#hP$f8%iY3K1GNIkdJ%Q#9HG zxQT8j{3j^pKeJ+Th&$@ynUJ~@e1DZ^3b=xosN-1pAUTmNZ*Yh3J%kvy!~cqJ?6|b9 z<#*=Pt9&3rL3q4KedCB;+N*@y*e#N9>U~EliH;tBA~~7g+D_y!o5uB*uY^J zzAy3)dZt}SGfG)VQT=VKNPry(8-%45XGg+cAgoanG^q`rKa^s2i$0cXrJHm`S3^NZ z;KHK>BT;oPAsc&QQNL-Bi6ew4uP8eLeNCy#UCuU{iolN1W#ep9^CtDV|`uuOZE`+{A z)8A1Vv(6g>418mFKgS4vz%mwBy$@kgp9Y@f{+923X?O`2MPp$g*P#1r>3%q0I15D_ zw^xGhy8#&jjbK#Q!&Q14E{{IrK5bZ9GHG1+Jl%}e(!HdOL%Q`49tG$PQajcfKtqPc z8qVW*0hf6QH)wPWmat3j-h2Zas-VDgnR0FIul?`bEmu00W zs4T|=R4(4nerb3Eu0A#~4MeJja>zvMrJ)OsPonizL;d@)JgnoWM~?0SwzzA&C$KI? z_6fi0Ay;8)T>w1Ybzo6F+4x|(?>578;n!)9RZ@1X>_IjFLK={Gu(qu4TJiH;gWQV1 z0}AbZ!zd@+O{9OEq((;ZL4EdwiMREP_39!Y7uAQUfdS@4{LUql30rKar6?Too2%)wcpV54%q z!S%!9$~T~U3lY7T8;Jyw-9)S#-~){wfb^=FIUgJ3TxytAOCr7ds>fpMEi7`ZU}`B`|DBpO+m zlmtU5GX?;;J@C^7X*=Qv+QDI8jsXEVQ!kN*XbiD12@4WW8gi129kCECW-@v+l~k1d zBav4w7)3mD9av4Of}t3|Iwko7W=f_~NG&7d#(q%Ir-ni(;x~M<@%^BmSA^xm#r8@r@j~O@KO=M z{@uh=bm3j*83SM>{hyCdF)dOQ4S&YEed@_$(w9miSxTtrm8@pDBKgUV*)pyiHD-K0 zIw6-)fr6Hip+b#PqW(#m9bvs5}WTx;O`r-${)3w{!Er`6NsfMvOMICCI z9AkPVhlm;r6-LygrU8rgifFyOl8vs$LR4#mP_&a$#F7P7eQeIG4T6?hD|CP?dChc_ zJ+(NPz5v7OHq>?*U*tm-4{E7_1D)^!X*&Ed{CN(4{;eoJCk^RE#+q%+t4Kc)5s6>lAdl0v z#zFX$|*dNvu0PPG$O=Y za9tXh!ohP;>jHvgi|9gFi8_M!mQxk%mfMM~Sn@|K@)SS`F^fwH#pM8KOM^h8$rxU0 zr!ksXo!0!u`o!n(eegea$KSxaBgAo08K_y|V+rbQlpH`TdpO?V#n?EFYy;6hNq9d9 zCGd@p&^ru+x)v*+Z0=-;%CLS6#ju^BJ(+vaQqbYG^OsIvrL`bBVY~RQ58rS9?XACx zY2E!-Y2lN;mU#z#H|xwfd^u}x<9>Ye@I-9xDR3X;b)vLxpZzDzjic{!o~TxnSpBMcW8FT;HvH6iv|PSR>@I5(=G`~$hlmOPRMP`ALp$ih;2Nm zyNk>89Go`_3RsSg%dF!Yb+5;3*NV(a%L}=B!-ya;%55&)W%3FJAyv?m!%?xPyHB*J-~jsNi?*ZGXN4ku+jVvTMn`w*h8-hzN9 zMT12VM9{S;q{@}7szMjOg#b(=8=j$L3VWuDrB@_s#w!e`g)gIIrcg*t39%=sNChd# za*oVU$~s!zzJHB7-olNHw-B?1BN?p;Tij>&9BC~&Qn-L@^*3UoQj4nxK+Nx_d&lUD zyDJ~2@0V_@NS+BN>8d1{Q!&9%9n!&cZh>yEe;tXb2EdN1eEV(VFqT3vewSqh2S{c^|y4$Tw8QKIvI~oX;@0v z^1#_SZLa6_5jQ$~@GL9;Z0dt3(BKQ(Rps z9Sd|ZIMF$a2RB2Qu5q~(E{56^9Vy%l;dL(PuE?5oYgM%kjwvGp)4-b=@WIEPSf{rR z^*Rn_$zFrkQP`zEX1?h^p6pfI<&?Ncz1_r&@@SKRHijIAoJ4C3okE>5@rVsB3^#{q z#7-ZRuU1TfsGrqRQ+Ib!t6puvW&xn%Ik!e#=^((?1l(d00*`-iz<#xqf6t~OH}46J zY|@O3C1^$lsnK-&LuDB&!%5ZMnT=!7`6w!-I+?D+g&MZ-D~F4YZ?v4tVWfivjCd|O zj~5lEGEE(}kr&Tdm6uW}iUJLZ+mv{bx6X4=K3$5a8?e20hSxK(gUa)JJhf(>`!rc~ z#6Ft|NOA z`otbG+i1TQkXXSZ>~?BEJlJ);s6)OAjThUYvU&rTuuDMzi8^~!lt_R^i4Zs{5gLKC zMyKdON+p`IOl*vobVZ4!qM=r+XoV^4k#trXAgV-t=da#L{*V^tO3R4uEFmmpdHSSn zWJ;+poZ{u>X}c%^H5ve`je0c1<)zVg@`69b1H^s!d6>MAXXfK}7+?3)jB>-KP zCuJsql&5$mf*31KneTuxCTC;y4YjN_M0zh^}OF%CNoq&yW> z$l_5!2$V;eCW6#<>4*<1pd zSiC8pZ#m>Ws^P8$JPt>eo9CDY80KpOw|PacR*TqiHECMVrqmYff^imgsdWqiYSQxA zj;2Vzo>mb-#NymKmjbUXn@^38tZI&^)j3%}hi(@ZaiLN=W?TZBHI(SWZ|%!0oA=LG z*tH(XnXauy&jdl=`yGj>}CoTCveSq&X!^xw7U>EG1t(tM{uCdq{y)N{{gRSG`!% z?YPj3-c*}5vn3p_E#;=6V*vD8+Paj%*zsZLUxA7$BbMkuS{~ar_q@_FIuEj_OX>lT zHLg&~lI3!*BRfvc(du!oN?ZA@PFEejLPD$qv;%;I*vuXtO8`)HP4Xp80j1EL91R)j zG}r+EOLZ>m_*=@hwICfR6=Ys!Pf)c;x3t%rQi`cL?1XKStG4JAs9J{|s3ceIZXGC% z0LMuql;?9Y^>gR#tsI>Qahv+n?;Nk6bnI;NSn0H#IEpOJl}oBX)5xAOY~pg+^Yf(1 z&9%k4+Q_9~Uf^3s+j{Z1E$V?P2eL|T#va-!Dm8dmM0t33e}_ZQH^Z^yt17OKc#cJB zbbp7VnFHGzYbM(^i28K)-%j>}UH-6kI09Cv-7s&q(Q(a!*EnmrI+?}nOw=lz23K3F z3%8}KfX$xD%SL%m&S^_J^vt^Wtc+{cw`b|az%}*axuldX>NqB63bVAmzoU^CS839e zc3$dkuhSJ7%heIo$qX{%rB=taoaf6U%QR*2c9G)C8aTylHM z0Gtp>+MKO4wkp>{Z!~va<>(G>jP?#mECs36zXZfQWS%}by$ z*YRowd)TZc4GYa~)`a#N#2h~Dr7sW8v_s8C)`wf(sC#i__-tC#hfDP8v=Ytf!5R6Q zw)Ei=y_(NX0WST$X0;1yJ+F4!4ruA`vluS@y~(iW17^}y`%Wj98Xss%>Dk24t8Lk# zXMCV3pGPm(m(&c|#eim-l$O7mtaUQJOz-AQ5^T0exyE6$lj$sqlyKT8JP#<_+_Lm^ zI*Ozqms&@!+Kbm_>{8=)#BJfD*Z8u*Z?p|BW$QSbnSilfYF4KKvqe1Mj(WJQB{^%b zjvGW9nxfYP_kslASxs%yEm$*4s*^z{v6u5Xhi+7}SEwV;M&z1w7o0uAfY%yZ0!~sJ z@FsO%Y1S`O1i?XdidH&n+xGS_OfXMXqT??z!7C3ccu_*3HF;~5+=P!Mx$Jq!$)D$K zd;*}hMO1)Cf!KEz3_Hft?YmNj2;flwOU-KK0f41??XyDx0JR6Ttxf?O4`bJ4ra;y0 zMF%QvLDn6@6sV@=ys%Dk)!OJ-G;vrvP6tYT+0lH*R8sgxOM#uQG$2jzLB5(d(Ozq7 zHL=&p$|}dBo2_p>?P{{O$XyqbjlVwzlI;Es&j%r~E&;Kx(+V?0FG{ z;?hjwEr2>Zb^;s4rJ1@t3II@_wX>zpTFv|~)3s}Mi>11dw5^^CWGzSsDixBZMZ5%) zpuDomy}tKdH!p0;b8sH5R3K}Mb`n{Fs&&|bN^;fi)`8LpXeUh1z%qRvl~u7+t3W0~ z{Q~J}(H4`*ks4W?D;EX0TU9#=$5s!Qy=5{vJa)XkD06kZbZl)Ib!c`W$$Zqr7OfFv zbqXZqqqu8@WC_QYbWB11|NHmZxj-gIHM&g&(!>^Bt6DHRgJsR5-q8F$yQ(k^ByGu; zbmVvKXv>eDuXq+Gs__*Gdg^Pj3E0e5N99?N&8*|c^ilGBJn~7o(Zk}+XB6dUsFxto zdwnW+1WV(IHXRR?hg&rsNY#m!yR+XiQB{d*o*txXPhz*Smpn9;%qyCur?9e3uglN0 zqW(mm=igW&r2NyHAeOJ=9bO0ryi{k5X~ZX)J3mNWc?FOnm{jS}wpx#g`t=q$iie0Q zXI3*COQG70(Lwk5R^}^nPRG;n1$to2=;3=V!(Q5*^mMY0J2~(sO%$%=KCiE4y{$Em ze+S;AB&9#~Yk`!KkII+exG2wmq9wd0fA5XqZ;;QI=IB_`QRl3l$5BV!)pficodvR9 z65r+A^L#)elQASRO9|O=A%WwEEhgah6J;mSZu^iAPG+-0&x=0zjAQe8NE(jg2-5K$ z^fTYa1HIPpY#427k||0%&hZ_G`HmZay8Wp~u1*LZN0A*gz7E=~cr_UwSu3K=hlj&g zro%_ap3TUZhXyN1wOb@{VBpfo9rb$857I2#{uZhJHbPN9J+0XvVOfs0_??lZxaFX` z^xhG7)V{sql8as$dAVDZ?@u1Af-q-kX+ z-iuE)D}EV>f5%7cPaUb~kzvswr^PEEz)OYDAmtw+ImN$$S#+~;2ijY$*`NRDS6)nJ zujW6hV4F%B6QKBrjA#UkvKdn${6+_7U-pof@E69%OF5W41xf(qIF#g{8EvBKIoXsp zntS*9t6eI4bD$h06z83QmcyYK-8qF){92fB(jAS=Fco+4D`QJ>%agwZx2EWr+q`PQ z8B}edDr_z|Z%UWU2NLV#)=8mT2GE)2vmu2sb> zK&y&oO=+6Ssa%msUoMo?CKVq9s!-#?GYUWL*jOI^fvS-g-BWTqdXidt=p7?hOd$*fCNne(ygsU7W1XsVvz)ZW;Jn%s=j z80$E3*RMJ*+@hf2pB*@Tr*^4IhyAYq3-SP zCcF@TQ;T!la4&x@6aD_w18Mu2q=QaOtSHt}Ugw#yvFGq@2<%L6I?bfB&!ij4xn9qMuj~Di--9(@t%ecu2Gaa(yW0Ma zIXVvNz-ztB>gd2b>TFJiJx{_Ji(&_>bcBu%0`hJu&yDS#jtLts=qRp(q`B!+3J@_k zRg`xfS4n`p1J~SqY-z5L-{HSwxYo8^9Sj{L{{BpvYdfIN%?pk90NcZ*AOMnsiuSdhX!1OV8zG*~9HD z$M@7>mg6ThUh|f3QV|1qPD&iW?z zmBPk*WG4h+r6J`t;)8dN4gld3gaJ(Q?qa=@H2OMTdwZyuMRdLt>u&Wo(Yjmxozdtf zWqy%Ds5x35TQPk!9BJz4&Z3_e&Fn*Jzap=AE3|g=r1tJk4B5>Sb=d5no!lgmR9C|; z2;#Mw6}*;2s5BY#K|rOsMYBwMvjDC3>RC(75I&RGN16!&xzxDuj3U0e=_s9FC>+Ui zyZx0dJsv=-RLMYD+uZ1kPDZ=CZHSD*K<&Wz4&6nHy`wytHN5u-iq7(-*ZuCyIWf*@ za&obdq}t5GUWOEqIyZJIi!?#a6xXJrDK3|xZf-+Xtjo|+5UnZD8cMjn;FVJ~nL8wW zCnf)ChwKn7_vPj8{++sMFHO}k+G2jIHW9+Vd9QS}MG17gFd8)k`buzhhU(O4pMfWl zUV+E=pGH5E*=Yyx)w(U$>pa{T?LPys$&P^3!5HJWMa6#xV3&%_StrZFdUCq84r2iE zp8?opOh7#Gc`D$XlaxhdernDQY(33~<%i}pV52K|4Oo64RjI=cNPaey3(_*K|0KYs z4%tT?OLg2@j9V9)r6A1`;K?G!)qI|-M!VsnrDC`cFBHQ4CkqRMX&^HPm+o zR5xCc_Z*OI<-wV!jMb&F!>56prO)hv?dhG*XOgehOUL7LkfaGwwtE*Y(W~_`4X)Np z2TlT@;b&dA4qSRdvY1azay}(pQzEfVC<;0 zb5wgZgVd5Tc;@logbLi z$)~KU)U^`8(`8;oPhya>r{d)sW0N;@Sqz1s&tHDeI@yvbBh2G15lULn z_(qJ!dcMlDw*lLH#qW@IOkS4sj&j^XD6M5Pnz}&qW#|YOjbje94#VpkvwmTKqhhk= zC5C31S3E+UNHxx1{wTgi*3p}rQ@u#qWXp=mQZvz1uzhxcp4;if91BuvN_cLW!<-~2 z!8*Pq3RZ39THggh`cRvI5)_~=jrUWj$8$1tZX3U&9?VV*eIp?0bHBC4SFR?ggB!y# zemhR?1Y#WRXY*>Yr|Tm7t?hnm%kx;9C*I1u&`F$b`u3dKJ-yaWH*C=&*>{+*$$ok4 zx~JEgoVaDDIl6Ot&8un6vg4hVq2-AEme_S`&);M+su7Wk?XoOAukHM!*AlA&DrJ$C zEU^A{;?_Cfqta{~E=9RPP1Ji1$VoF*KykpqTN1_sJl5W40f6lhz~Wt}&A&(X)_)Dm zf22~;S(*`uo?iJA-yRayor6RkcjA~)<`%rsF#8!R0?Q3lR7z3OG;@S|=qGz?fkDQ14Lj1?9%FR3_vmH6#lcS`uI?`;D z=aIMrR`q`fD{&a!kKfwE7qJF{OcR#LLOe_}%Ddh__L@$v%O{4G$GBVL9WNptV31LS zf9Fln-YyewIVSE*3@Q_M#_lh3Y;=O9Ugl_+czhVH__94Wjb7O3x3yy=b{bvU;vA*9 z_7kjVBi5M~b&vy~xj5mqSj3I9Khl}lfTLla=BUZ3qy0Dx8zlX!b-jVFz2 zG9BHfpH^rtQF{vBl+z|SjNw?OG>Ndv{Y(`-xJ)Apn-_b`5!+~ zgdQy~>7#hM2wgRGo7B{Mye8qsZrSq3&j4#)<7TYz{g}Rp0j$`s%7#$NB4~_kb5!)` zUg~_9s6b(S-gtT^bDUbG{cc5E%XKf$ryu1A5^ZCNpUn$v(FbxGYrEQQ+qq!0F#R>G zPD>+9)sv{eYa*3`Bso|~oN^qyRk0B!o}Ow#DnF&vH^<#B^M^AI!ibVLMKqE(S-16wyhjbScUk?P*2wdDjm4q>c?@|wJCkSa@{9Mw($y(XwPpj-X-zs- zH5R#4 z;u0CAmG~$KW+lGi4TP7U^*Y^SobGc?s&XiV0M?N(H#c7$^_g#DFUTcNPG`Nvy{u+2 zRY=jCd+|1Xl1|eGw^VsO6{Z+yAWAl35PATZIMq|=F2KwwG+3(VwjIXGu+wY|-&_Nba@_fL$7<4V6!YI~?Wj>tL};|P4$(0N787U^%7uTj1*3N*WdhIz_L%i% ztQ{M;osoz4+dY5wvmd9&+2G%?*b=tNndu)&p}Ur?j&^rsUeaiI()%m4*?T&?onc4y zrnOv$)wrlAAM?Auf*!CCkcO_m^Tck)P6Ff%Kf(pmi@pfwZ3e-=L{-`6%J;T27{ zoQF;Y&g+GWr3EEcme)>ZpqXI!;Xfti6+p$jr2{v0-r~+5y_rb+go@WT+kLdO7uejq zrG>24(a^~=kuCm-(D0}KPZj-j>Ew1;c&x@mPHL z0{?(0^LsBq*t_5B*YQ&c+|h zT#p#*4I@e_HI9))Deg0;V{wek8m=HiEi(PSgBvDHgKobnw^dYuj$W2;s}y>h>F3)D ztqe+(8tV0c5e6lo7{nL`Lp@+5Ad`?lU@kdC(Ku0lJPf|8&LYE*5R(HQr?!~- zi|T@eHp3RmI31jkgaS+ zF0oD|5d`FXX(Tx%tZSVDq&s92@Nvrk(PS2kt0)VZtonb9w=q%m$i=}Qd+4#)7ja9J z`XaC&UG1m@#*N|q;c3e9EKm&VlS%q}DZ^?ABaLX#A@eGG94OeQM5fZg0MtBy!SOn_ z2Z4txkFrZFvGL<&{kT?|mvBeLG7dtIkNQet;UTmmLt~#kN=9BE70j-k-jh?yd*LJo zR%TF4kE2W^CBAl@JP#jDNI(bcnNz4?T{s2oaUPstp^0@HJs51H6l7^Q+!7;F{6Uvs zrnmr*8l~EY^Wd~Gbhx>16;tly|{C|U`*x>92y^@0EqL3>%&Kfz|iUtC&5Md81UUfGlwSAPzKbtM(__{ z^eOs-i{U|_Pk+=5hCWYTWZnY$t4dFmlvzo)CbF^`w!zA|$O@?x;(d}nD$Y$7?-LkR Z<20I8@wN#5a~ORs;#7IIuU)(6{{na4c@_Ww literal 75819 zcmeHwd#q;Hb(inF>?BSU-zK9JD$-DCm`;P;OyY)4)S+C$K<*{CQKm$sax=`>cYNk( zUYR?CAJtT2F?MYDszCgi$Vf&gDlYs1hybNg>q;oJR$_szCjKMh4fm z`1WU(iBD_~KUj`ta73OUuWveP+oG zk4nStt&ewi;NsR|<%;%lMS7uo)H_;Vk=AlyYq9rAcPX^pM;`~RaVc-Va=-R`KiWRR zecq2Yk1V&YE#A0&_4L)Pq=Y>Wo(t{wK6oy4en()1u7tkLmG`nNB!{p|ZQfU0xx2XX zL*~lXV&8-3PTd!V-H*?aQ~oB*k}mig24pGt8wN|?Fi`r2K}g8O$~Fnz-OI-5-tahS zl5~D;v2W|a%VPOk@9oZ=+9nu~B^PfK09-hI^~$+>x3=%RZ|mar)^f4>;cHhfOJsVC^tL!GZO>l2 zc<{J8Cy9UmrHhv?Y(Ysh!1lR|7~tXVjotO`Sa%j3;134) z6#dcFweGR*Q;WU#UfM>>pGkKPcQ?D^bn}VsjktJ~o8uxVo=g{UV|!x-n&}=sOg)0= z8Pofx<+JwpRNm9m_$qzu;XZz5``X6djg`lii$5ijik0qTYml@LUcR!my-ku|IsJfK zTrR%rFT|@KJ^g?fPE>}m;rx5q?(eWYM!Rp{cY$ z=;4KAmT20a4RCor7#CzOgwLRPSw8D{+(jDC3=hJu<>T@2xU{_?F!+kXp!+YMGrSAs>h7W*lTY?1ZGO!KAg-;sde&8qkylnxG( zLOWJ*e}Ix!x~NhcG2tjX(eiI}ze$YHApLDL#0;hOw^REOy(aKP$r_sO_fU@;EB$xl zpFRB{`bSjW{{Vd)YF#b{#V_`+$KW0OE_&R#3l}78FBfn7kbcfXc=yQS)(>60x^+KA z$|os?L(~2H!Nvi6T)f?c%zH0U_&axzqvrk}zjpaN;w3`SjqVuJtpA@^UU>y%>_H() zJ%6w-mgK+fnObaE%bby1GGk3#AW3@@ynwV=WD__Myqfuy)$QM z)!(CB7=F&C&x218A0791gY1X#zYpMl$LW7$2$ADOR?G>N$Jf6G7s$iG#b|al8gye% zZq4=XYeC1G-SgqoP4Z<8PdCY(I6Nt%h)c&Pg02R`z7UDlI%WaiwQ+z!Zde{BQ0W~p zEcX@eP~7f+l$5|<*^ujz!Obp>x~|cajJnx#hHA=L%~D`TDh3x+{uYXm0LpXmPUMBD#Y`XBgF)`MuOTIn-#OtgAOWZ} zXHpqy8nI59N}Pye-7xGG4%L4e!d>rHtr<8EfXV)|a(CbgF|X*r{m)T5NQb}vpVOEA zUo)lZuKz_`dycOW5Z7Z5G&JFo2NZ*n-VPFG@7$x9jZGUGWIb_1G#aRv6&pp8n-RSI ze-@r{`H?q6FTz{+(gb7wyF?Y0utEZ`CMbFZ0I^X0n^dwaxs!(9y2>wjc2Z~x;#!{{)WP|ZK4j+rA({>h59?Njb>n4q&Xg?iLF@w1LgLA7tnej@Tg3Y464ii6LgQh zvv#<9mRAM}YV1Epoz(SpJ!(M2M@t=v3c-I##MdxSJMh@jDbV{Hbf3z`=g^Ex^!_*K zW2hGo$ZV4aEW*n9Bcu z96B^&+>=}x8t5NI^H0(LNJNycDbE(>wgz_~h`fP_)(aS?Vz*R@ob7*w2t_MW zYD_Z+iLnWbS%?OUK_cl;L-2`yNd~~R{!4*Pdc!xPz$_5D{Ht! z#!c`L=?r({)ILz33gh-a23yJ4Qt&c(L#s2|>brp-k zFmOu7%CN=^6;nNpZPEQP9Du8_t^WYwT`y|9;t0uMN$+sH>Ni4LTrES?E2I|Tp7c)< z2cK1yl%0Q)i0gspjt~Y#rmdi{6`?qWpY=r0jPViH3jXi6eV|5+LuS6H)7==kG-f8v|H zdnii`9Xi*AcSj3o>y7>Wzo3CyHPY|UKvX0BDw;`+bpJj=yDTxM8i_Ba#Y~zY(x%K< z(eH$m{`a8QKVZG42;pd(r!tC04huX_!4>?8HV?~M952Igf$=`-7`OX7NSI}zOnKR_8{}CFa{}2xps*e7n zfp?5&s*O0Gl~oyxfIl)h`{$m1i+ZWK>A`-DdKKrwl)C>cVMA0;!Av9ymf%R9B*GX% z*Z&Qo5vM#sohdX!r@2$oj8)aOTBvKaP*?wisG4#Hx3;XL*3x6M2KYHZ9v&+Ka6fJt zaD^n0{?HZ^dO3d#QjUsg9030PT?^M)qj=m`?)>x1x*Ex`=5gPkwrg2V&F#q z>(qs9!5I1!n!VVc35|?oLnj^0+Wr5|WndAE@k545CRP<0rZNv5GS$Bgs6WWKQV#_P zKp#Cc`Yf)JM;@Z)BZ>21-%~Hapupgt+;1(vZc!P9&I&$ zCHhgcp(}_g~a(X=?GU?Wsav^4dMh()84?UN*6s1H^aqR z`urUlu9ZGN36rh#dq_VaeSQ+E8m3Ras9RDiYhiKq-Nb9p{|h@5hiSj*wvW=;kV_Zo z2nlUCt!&)R3_h~h^O1A1Gn10nAJ|1LR6QC!4^f*Fy#K+(H@fd(7vnEjgL^;8z43$C zyP$^0xw__iPjp|zMNxSFI0U)Tt;PHOg09F-I5#;g&nl(6cYw10o7954y52A8WB<1Y z-*KPzh3zo0UH|v!W-zvHnQiQgJ1md_iyT3u)%YcF)c*qvh;0DoqJNgY25&RBz)-w} z+xSBd+VlvJc&s}@l7DX;g5?aw#%-*(gEGiXGsrK{XL3It!Cr^9*S{HTMGc^X9t~Jv z9ru?9)WyyIb7*mcR=rgVEJyc5zmEn0Q?@H@y_vV}f0a53?T8h0DXjcFK3SD9KgTI0 zZqk^?BABu^j3c%X^_4Wwve~}A3eH1|sGh7R)*VJeB0*0o3*T77jsB$}i}4=ve?f%c z$nyPLmoDD^^6)lG4scsj0s6|5X^^BT8dzxxW714A$=W^XnvmP5XnY#ri`)(lFf=!{ zMd*i$$Xq9Kqi1{{=NWb|#y5siv8IDI+1)eR?5;sl&l(w6ssBx{bCsvm2*~4IhZglF zLjX*-!9xoDZE9qdlmM1+WSyap4I&TLqUBvHe!6Qqw<2(ZOnd*sC@0>HLmxPh(3b9+ z&Wq94;5nor_Cm+3gS(mf1N*}!j{L-+9i!}JZHTT2A^-`0KyLsQCA)op>Gb)t;l#mD zlS(+&aKQXSo7}J{5|dlc0C)_ov8XB4<66i%O#Ay9Y>%6il^?;Nv}V#xOl4n1ugA47 zvI4b(9wf)eP{(shT5Zz&q5rDWydND0q(jtr?p|(W6M$lQ!n&K*V7zaUn;VFax()&6UQ=3dDn!#KzDgljEh{%os zaAL=RmS>A1u8k9PIWn&v0*3HsfYbi~DYge*`n&L{2@g?YFzI_InKaZj3(n78MT>AlAwC zU)i_~ml+`Be0+KdN%2pU@%ct{M~V8hM?`Ui$~%CA5^0tf!o^JsQ)x*dw^w!|IR(Is zTwp^jjXWWYP)kc5%)6cq6fSOB@)U!)vaw2y(fP~p4Rc30Mg)K{jLP)}=MRs9BYH-` zJfSO?P{&9xk?Fu&xk?tzqAB!FQ3Q+g>j6&CUy8wS9_T(tBT>f4BuIdQZJQvIIsn2M z@ul!k2;QW$g|s2Dm&+9zOzvv+vKbo7ckp8M6lk)+wPj`=?JkEAlY@8LrzT!2q81`K zwnEDsTcI+@bSC6btp!|Yndwqn<|t%7p+v$-&}q+cFhUEU^X~w3dWT00APLtJZea<| ziQ3wZT53pPt5}1x+5SBWfvr+toNO#fJ+vCeE(FGHAmPsaCUGp=_Vn@xqjT7U0yU`u z!K@DQGDm@&(vt#(&yb&rWUAt#&pCCGq!#JB=eo#(F{C)=IA~(GNRIYPuqHU5E_po` z$uO`K$vCNn=~H*8e+;1tA8)eX1Q#IF;uHQ9>VWYP!e$6%H(7Qxfb9A&9|{EUBm#Di zS+|XIvyl7+Qj)D?6upw%O|D3uvTY=9hAZMZCgz1;#O0K;B$x?LNYKZ7xyj7tuvsw)W{a5G) zkfkAjz)y&7`v1Ckn0kAv`SnW5zUu!dUUP3xzNp@w^1fb4aU}SE|HHSDH-3Qr$8i9? zf|veBQI|d_i7tgMZblqQsfQKeE^xpOa)~RafUP=r4?>Zp$G6vW@qR zn5_~}chiH(FH-keZS@>3>`sEyY=y%fTpZOEK%iERH&C~licSz0<~NZ7Q~|aBkyK13 z;_82Q`m1x|-RUo?ga0?DzyJRhf3M&ChpZFV-m;ZiSsx0?C@tF(!K+D%X$6*#EofPQ z#aG^#$xbZjA_*TfkcF6(f92wrDHO-Dn=6fw*!2*DoI1cFj)>BK^bJG-GW(quc@Vur zaE^^)h9|kfB#sb(P?F;fn*QisyHFnHwJ>pD?(`EhoaPEFlgTY+TW@J%M-SRG>1huS z7Z#8Pd70-qI(gz5x~wS%DL6N!B1KCAOuKl+E-Smz>_vCXWLhMBkydGYRXNmumD`KZ zsm~(wpaU_BqXUKgFJ@7N3o|-nev#eSN9H=T$vxDkav2!JYQLR@vRZZ*;l$2i)oH>a z97Xf2u-$R->DI~p!%~2!sW6%8CR1D zDoEB3F2?gh_Cys8jMHi1yw$VW7<@=85bDT$lubLc(%Qy$Xj$gr-4r5&hB?@yKmCKC z$7#VDH>oyrPOz?jBOLP?tYp)!9?sNKox?|K(+Z3T3vmECSRa21UHpqbg3iYUI$W~q zljB$pPZ&ZqV(hMpNG1qZM9{S$q&AD@v$~2GEi%j)88RY_Af0v_FXpI-9UD{t^9n6p zeG+IHSFmW71rHU(+EQxA!hvW6v~J;Rxb>C6$O;BH61&=+(VDKq&<~Q2%j5A@h!FELPr&jJc->d0xwO< zeo$3{2c?p{ssW#Th{Q6y1@S%d(hp|IUXqty*c^Sx{4CWKWtWrU!szuPW)MfK4Agr_ zZpcBj?4gksVk920!HHpKJq_4t_vB*&Qy{8mHPw{WUBjVp6b01bvA0TIX>@=SE`Y_P z4m|$B0pn_vf6t~OH%}M`HmQ2XNH+t6l-)G&p|T8>;e_hu%*L_k{4NTmx=E`VDAlkT zt{g5J*rk?E(dfTa7*z7SSwcoohJICx%eooY7r&AfNuB0Dlez1elWphLb2jXS4k00|fYy{p=T&d!Ju ziO?t!tOq4RBjCsA6g5b$M3a_@4b?zm43;o*vpRA`D^1ocu7%~Qgz!n>{ngtE4{@o@ z)19@bo-(xU!g+F15GqP?H$TPp~NoPZ#OcEyaX0JROj| zwZUU^yaV$hhq8la9izTvGewSD&j&F(3R^QKfKH)PB-&^>i>{nTHsnjO!sf>D@i2u6 z@|93)ju+{wdhvWmLOW$THwtahWuTEU&6}4I^Vk13A5$BG^ zsmdL#VF-*SC7-QmG?rEi6?l;B0#?aeYYj{>~1a9r0SvGH5bw6;tnrLhMU~47>BoQ>oYG8Gvsi604h? zdOR?`pCUR|TtuJ;Na#3fp7~QAXSXQ?R8I!5WV|cLFhXEFFBjUggSWIWgUa~rI92rAsvSKJBb5mqxgtEO zC|9YCMgogtmB49GDVA+?%sKRID>cucs9Wi)GqD&eOWvfil2^}^XH+74Na&P755~hb zlbtW{>Im2@>QMO1tOSQkN~u<87yy--RxGYe9zr0K5@lItv^29kIx3mS=sZfhNEx|> ztl3KCIO-^jPL9!vI9H_}q!FXbaMgea2ucTS5Fj8{qld>52ug=(m`FoFE_5L9Mt}@u z7;F%Lr81TUp*^I=SXK(spi)8V7V*Z4RHaC_)MHZ&F(ro$RY$o>i%y{`b=aVaa+Ppv zP^kwvAPrO=BEl+qjv%pAEk?(Ic+&Nc-cng^@hoQJm^Y-mdLCGuD;Hti4d@MJ*udqo zHzt$AV*~Xz70h#d%cwJE10`8LP!E8t5}R>W`y{0rJOri(e1ilYdZ6Rojg7pT^z?+{ zVs~T1Y{_K{Sapq6gKeEeF`e++aCq3|+qJRL^YW(0UbT!4YZkn;SrRrPi`kGURXB~V z%vepjxU5~Rsl3>kw&I+!q%odZ4li?LlP>Z(OD|?Qys8%ul9IZp;&8KxM`2d9cQ-Q9 znkr42;)Y9E?KSe9v0UjvjmRLgNGf$)%XwHHSk6hy4PYAL-HnW;u{%77tPqy?b7M5k z$0?Tzz*+GDR!^hNtI0r1p01k_<+u%$LP|qhMXHCpTtiV_s0x3WL5F$wP5eMpS=f08 z14HJihM)#bxsMtqW8NWiw$jk5oDaR$C^w6fG0oQb`(AamWQ3$te?tQ1$(^@0NAb=^ zG)9pEvW$EqB>s_|Dq~U(mEt!XEmayrD>hD8Y7%xy8qlh25C`@egq4^YRxoQ&^D6cc zlKch1TIDtsE-!W-s$zu1i!MzP?_5zMACXz*w@xs?{h;_W0!3dkRq+zKUW>sph+)M7`hLDC2#d4@x$a>?Ma;RFUdwyj2tXkBk zQ(t$IX)vcpXUe^*Eq%H;Udd-80q1zHMnDffR52Ljm7r~mn&W*I!#Uoo40{+bldcjw zjaX_hP?gfNi5ai7Wn(;pfvSAocn0GhpVcraC4VtlYeamR-p!dL*k}=BO|03{nnfeh z*)C$hsdwQaATM+C)YGUaqJrGf8hX`QyfkBTG%iQn5=*-FWCP_q|l1J6ce z&%dkHMVMq6@|mL}u&f&L>gYAWjlz`{0M#j44%WKu?O_;x9Z-1J3LSru3Erdy>B(EG z#C3e$q=-F-?tEdEK zbZhk_K=d8|s3U$}8l^*laDJANB($NNwvz!&RcP=azKZHDU!D!m=~3Slv+NU&s>Iv@ z>w#!$f^4`vIsD`VsqXHw=lLpM7Ngf;6a9kh(>@q!p>Wv|3FJp4;T-;JmS73NJ0%Fi@l_b=aVaa+Ppv zP^kwvAZ;KbWS<;CVyRk!j03e3q^CtE4rMuc_v{h>IFGJbSEo@~nwy1^(X&S;K~g*l zyH<#O0KSA{RAL;c!dJvCxL^p@lv&WgK%TE;XW!`xs7wG3SBx+zme^*|IgBF`?Vvp5 z1Q~&9Vv#nqUP9N187FC$r-ZAIvBQC&EOjKjgI6BXksrFJEkAm#;!T{W#v>B+ve#@A zuo|rn)!&dp+xO3pdXD90knt&&N@xd>vp**YOE|06hNG=-ktQPcnCU zYq;_XBJE&OriWKr!~|1%iyYazLzOeDnT=7Zx?%P@Z>G#1nTuwbBo@6VX7ENWmtilJ z&}&%+q|VR{I63qtOk_}U-~3pzUe}t(ze8_gpN>EEu|SHIhvacKPRjFNWC_&pbUD#K zOPJ5i(WpKgb(#b)9u*vQlvP~=^@#t6D1=A$UCuq^gA&1yDXS6CGojGEgh(HW?c;=M zM0rSU(%5wj`Q&6aE41SR4VwWoj?L$FQUQ)5NP`p9Gq6E?C6S}XK<>{{Owf&{n1z5` zji%tpjp)H(jaheK;3%>YkbE6DdTiumAW|pp29u^#%Eg(kOoxw#J*%EEuNxGrB`l&^ zVc^`z4fT4j4H{R5{!~nsKRvD4AKx7BR?Y-Ty*#BgGE>&N;SBai%0u7H*2p?j`<4?d z(pn!)XutsW(;srq#q|6L4v_z5bEfq=NttLGIAER}7><%u(@~wVvY;4KR-BfU*lv=N z+o#|eAHv?@&(?%OH@wU&`w>U}{Z_F*RiyF*(fGThzX0vSqPLpe$);GNK|S<(ot8>djRUx^k@nysCHUgbOcb?MQ*?GLPJ zh;CF29CXzds>15J)4CKDd?2w5ZjBJSX$YNZPW?3w60WNSj8_#)9;uGvliK@s@p= zYso~>{c7EDI-lAWz>u~BKK(pCyCRv1BVhRRN?41})Q+Qbx9Y!cpa0pd`td7U3hhQtkO-~bJdp(9R+E*6 zTObyH@V}|aIc&J+KWpq&Qq?G$7P&`VghQ4W9h}6)dDUxM@!i*bB#{U7X~PJqZs)%Y z`dSzlXAMVt;Q-JR{CDU8M^WW_F+ zOX_XgW&4D2HK9maqrA>DV`b05<=V>+&OTFF8#anbXP>D$x%9?&=*xUR@-wjJq$7;* z5>|8KzExlTpj`uCjylTPdLF&jA_WpAs-rp?_K<{=7unWJ4ursYjS_nUmF=De1S$}m zhvVv|PcA@U!>NjR*RkBIamgpi@9^I+d~MsVbcRMJ{{BoZLeqp~ifP-IH4WtCCG;CG zK(XG~TuSh7_Tg?&Rr&C603?+-s0U=v4eR@Dv+P*~bsbxtM3r4g)e7*DW2C8Je`{lB zWttQUl>N{2?|i##d*M>fO#|?aO~gRu>0`J;U7LW>lh-hpR2Sn_{W?#=rnYugksHCu zt|A9c8Ui{i>y6g(eE#GbiYL{%g_?N5il;(tq7(azq|H(6B54|>0UcQzn3Y^>Yv0XL z>_WDBn*KXvv{EGRaEqwl3N8iHs=>UP+C37d1iO2LpIOp@R3x^_Z*?|sX);N5kFc`g zHh5~Wx*IiI8?84iH_#eS8oMVVs(o*kp|82!hIo(#%4Tag>D1PZwNWn>uuCh{$7gp_ zcFnwc!0s`X^>`!hOOz^$h9(I@PBr3sm79!GWy~kRDCJbI5ew}O-7G?@y>fCFM)a9P zCYebBxil>R8O~ycdjK`PYk#~{n6Nbj4S7R5C}v|%Z41$)?O>tO!z~urs^F^;lsIaD z0wbt_lD0cJ>oTfHW-gPw?if5w);6%^X%sSeZ7VY8#Zx08=!)g!kOP8Q&H{oGzG1Vi zKlBcI;9(qnyOpnRXbLAYEVofi8k|J`;9FBQ+1^MSjS&;z&K~BuyB8XeS*Nyr zod??xwkSB9;LWq0P2yq4*Up|(%!QPcnP4~S52C(W7H=*W+kSgDM*C0U>8;@r!jpq- zHePt~3O79dH+YA7qs*hH8_QQ6u&E1nmVLdtZ zL4z@n_^-gMiZ3Ld_&gzW-u5ATDE!i>XUSJ%g=d9o)ABo`Dzx>wjTZ5dbg>7KXO{*^ zdNrC$QZugqBFw4|#kkc<&KpIHTTbK}o**P8ghogjFT!*(ASoFT;wP10pI9;>IP@bfVOGQAm(ASFGMSao-gSl~H<*wZPbBTn65-gG4-_XrsW-+{QP z(Ya&{R9YY|(bgOdWOye+6>OeXBVr=hc1Rx)SEu(4Gx{z`j|@4(8m3fL!%GGUXCA|> zWdn)o$aP4-NF~`!zj3&3^dLDXcb<7XcMJa9!SKXS{D>7ok z%ee1)(N$%qUWfW5s4J=g#CJUs!||F1@zfTLjvYbkA@F* zIgA7aD5UXu>geHY89lbJwn8fycsX8pb2>c`6ZKimXt^IWqSx+GN;mcGmbTNM<5ajO z4{5{?_FU)e$|y;h`2&d;AHMTS+t)HKwA_8-_AF(*G#MPPwB1vQ_26W(>3Ar6eZQw^ z41`nd%G0FThq_P`mX0x+pms~*!RA;a`ov0gvL#>Bl2}I#r;%^a%Yffc1;P67zw*i} zkA~VJ{-jFoOkP!%onBt)zfLcJrN^B2UR^%Av}-}>_n!v^T$LdHmXHwSHQf4_h+8Mm zsx%w>b5U*(tdH!Gvm~m3>=3ZGB#Z@kti8`70NVqC#hXu?|IRGeUjXy(D^!Gc6pNmI z^iRCHPFQzNC-Sx*$BZ(!;EjgauNV?q9(+Kv-Lfk|7^=}oZX+$ACEi|nC9dM`=dDeCt?oQEQ-fxwO zHyslqqg6mk%e_%1ZpCb1(Z$3hM zor8E+f0G++syg$dMokZb>MBKr8)u`WQORP2vdnA}%s%Xu{m4OH$kb%ebYdcKOReF^ zQMLH>hY3{GwQ7%QU5=eRiwdE%Ov7CXCuv6+t ztAOgPFK{`R2*(7xUH*m^0oQcYq-p@qa@b_n4enei4gJeX>`SmVOb$BzE9~dRkz-tw z@jAC`T@T-^LqD+~;Ujy?4qY_~oYd5M%`DY7*ftwKZH!)_vM?Ze1@=SuA_g$MAr(Qw z7)4Ma*(xb3);*E$p;4eP8$TNs-CcFwXgH;;scAH2?|J7s2*?(P` zJf&ZkAFUc)y7G);b|;;Bu^!OuE^~C{A`(_bbx4AL`6NJ)f({5$^n_qXah<~*z5FGA z!YCTIw*AY!dN}k_)~g9|gRq=E2Nk(Q1zQVz&fY@*SN$TKryqW8y}|6-Dx)AW-AYyn0d@PJcRjFIkOHP?~R~8 z>!+5NO1pt0uvAWYM+memNU@|r02xeClQI&xuyq@}h7yZJexM zHOM?AZ#Gv58WBHZ(h1-VQrM|AEpDJr=k6N)6D}G=GA+rb#fdK9$+H7A+FG7yKoNY- zp1uckMP#U(Pz;>!Ok?l1Dt66i?)@{|8!+wiWdHF{>-ZfZysU~md2ge>Q*^zv=+(?tCN$bcAGr4>;FC>Q|*zJ zT^^UL=Dc-1)L16xS<_5&LDtSBFSt@8|R^0Nv;DS>ri> z<^7UA=Wj?auNkh%ca)y1%J}PVC+MgxSj40B2v63eQxVls*p`_PYKM9Yua8`ivJWz3 zHj0XxMj5tHW!V<_C}kf*h1`}DlDZkF#A;0&!F2gpLvD_S_pvAA4`nV}4E2T<`6y){ z1Bs&FCr*dr7@5^w!<)0|7jIlQ;W6mg)QfUkMHT9RvUFRe(AzwIx~u?uxq_X08peLC7bA3oirtItRpqWzhp z-P7SSRnAO*SNC*mc~n~BJ6d5pX?>UqZP+zzc=d1jKGcTPn~sN3o;RcLh@i)x_D7W% zeAwZgJM90M^y>CSqa}`?2^^C&_@w>U`F%I_+~1kHf5EWIY?rLYb8-1|Dw4rS+I8z? zoYv^#W{gI1g`n13^dhy9);2R*WCwH@Jwb|n>Mf(?X27n3MVcWXMgSg_g1We1$0HQn z6`E;J>5?5WT%`&r^(V9_noMM7Tyb?6c}tgKkE&wSavaK!m^_93u zY(^9`%?Th*a3@g0x)u}|x-o3@m7|sc8k1QtP*E1ru!{fSZ^Jkh) Array: + """ + Returns the actuation matrix that maps the actuation space to the configuration space. + Args: + forward_kinematics_fn: function to compute the forward kinematics + jacobian_fn: function to compute the Jacobian + params: dictionary with robot parameters + B_xi: strain basis matrix + xi_eq: equilibrium strains as array of shape (n_xi,) + q: configuration of the robot + Returns: + A: actuation matrix of shape (n_xi, n_act) where n_xi is the number of strains and + n_act is the number of actuators + """ + # extract the parameters + l = params["l"] + # map the configuration to the strains + xi = xi_eq + B_xi @ q + + def compute_actuation_matrix_for_segment( + segment_idx, d_sm: Array, + ) -> Array: + """ + Compute the actuation matrix for a single segment. + We assume that each segment is actuated by num_segment_tendons that are routed at a distance of d from the segment's backbone, + respectively, and attached to the segment's distal end. We assume that the motor is located at the base of the robot and that the + tendons are routed through all proximal segments. + The control inputs u1 and u2 are the tensions (i.e., forces) applied by the two tendons. + + Args: + segment_idx: index of the segment + d_sm: distance of the tendons from the segment's backbone (shape: (num_segment_tendons,)) + Returns: + A_sm: actuation matrix of shape (n_xi, num_segment_tendons) + """ + num_segment_tendons = d_sm.shape[0] + + A_sm = [] + for j in range(num_segment_tendons): + d = d_sm[j] + + """ + A_d = [] + for i in range(0, segment_idx + 1): + # length of the i-th segment + l_i = l[i] + # strain of the i-th segment + xi_i = xi[3 * i:3 * (i + 1)] + A_d.append(jnp.array([ + d * l_i * jnp.sqrt(xi_i[1]**2 + xi_i[2]**2), # the derivative of the tendon length with respect to the bending strain of the i-th segment + l_i * xi_i[1] * (1 + d * xi_i[0]) / jnp.sqrt(xi_i[1]**2 + xi_i[2]**2), # the derivative of the tendon length with respect to the shear strain of the i-th segment + l_i * xi_i[2] * (1 + d * xi_i[0]) / jnp.sqrt(xi_i[1]**2 + xi_i[2]**2), # the derivative of the tendon length with respect to the axial strain of the i-th segment + ])) + """ + def compute_A_d(l_i: Array, xi_i: Array) -> Array: + print("l_i.shape", l_i.shape, "xi_i.shape", xi_i.shape) + sigma_norm = jnp.sqrt(xi_i[1] ** 2 + xi_i[2] ** 2) + return jnp.array([ + d * l_i * sigma_norm, + l_i * xi_i[1] * (1 + d * xi_i[0]) / sigma_norm, + l_i * xi_i[2] * (1 + d * xi_i[0]) / sigma_norm, + ]) + A_d = vmap(compute_A_d)(l[:j+1], xi[: 3 * (j + 1)].reshape(-1, 3)) + + # concatenate the derivatives for all segments + A_d = jnp.concatenate(A_d, axis=0) + A_sm.append(A_d) + + # stack the actuation matrices for all tendons + A_sm = jnp.stack(A_sm, axis=1) + print("A_sm.shape", A_sm.shape) + + return A_sm + + segment_indices = jnp.arange(num_segments) + A_sms = vmap(compute_actuation_matrix_for_segment)( + segment_indices, params["d"], + ) + print("A_sms.shape", A_sms.shape) + # concatenate the actuation matrices for all tendons + A = jnp.concatenate(A_sms, axis=1) + print("A.shape", A.shape) + + # apply the actuation_basis + A = A @ actuation_basis + + return A + + return planar_pcs_factory( + *args, actuation_mapping_fn=actuation_mapping_fn, **kwargs + ) From 2d0a009157f8e37794d93da7f6e64ff00f19814a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20St=C3=B6lzle?= Date: Mon, 26 May 2025 18:32:04 +0200 Subject: [PATCH 3/5] Fix some bugs in the tendon actuated planar pcs implementation --- examples/derive_planar_pcs.py | 2 +- .../simulate_tendon_actuated_planar_pcs.py | 39 +++++------ src/jsrm/symbolic_derivation/planar_pcs.py | 1 - .../symbolic_expressions/planar_pcs_ns-1.dill | Bin 7890 -> 7886 bytes .../symbolic_expressions/planar_pcs_ns-2.dill | Bin 70480 -> 70395 bytes .../systems/tendon_actuated_planar_pcs.py | 65 +++++++++--------- 6 files changed, 53 insertions(+), 54 deletions(-) diff --git a/examples/derive_planar_pcs.py b/examples/derive_planar_pcs.py index 23f1664..8923f51 100644 --- a/examples/derive_planar_pcs.py +++ b/examples/derive_planar_pcs.py @@ -3,7 +3,7 @@ import jsrm from jsrm.symbolic_derivation.planar_pcs import symbolically_derive_planar_pcs_model -NUM_SEGMENTS = 2 +NUM_SEGMENTS = 1 if __name__ == "__main__": sym_exp_filepath = ( diff --git a/examples/simulate_tendon_actuated_planar_pcs.py b/examples/simulate_tendon_actuated_planar_pcs.py index 55025e6..cb04fc0 100644 --- a/examples/simulate_tendon_actuated_planar_pcs.py +++ b/examples/simulate_tendon_actuated_planar_pcs.py @@ -36,7 +36,6 @@ "G": 1e3 * jnp.ones((num_segments,)), # Shear modulus [Pa] "d": 2e-2 * jnp.array([[1.0, -1.0]]).repeat(num_segments, axis=0), # distance of tendons from the central axis [m] } -print("params d =\n", params["d"]) params["D"] = 1e-3 * jnp.diag( (jnp.repeat( jnp.array([[1e0, 1e3, 1e3]]), num_segments, axis=0 @@ -53,7 +52,7 @@ # set simulation parameters dt = 1e-4 # time step -ts = jnp.arange(0.0, 2, dt) # time steps +ts = jnp.arange(0.0, 10.0, dt) # time steps skip_step = 10 # how many time steps to skip in between video frames video_ts = ts[::skip_step] # time steps for video @@ -100,34 +99,34 @@ def draw_robot( if __name__ == "__main__": strain_basis, forward_kinematics_fn, dynamical_matrices_fn, auxiliary_fns = ( - planar_pcs.factory(sym_exp_filepath, strain_selector) + planar_pcs.factory(num_segments, sym_exp_filepath, strain_selector) ) + actuation_mapping_fn = auxiliary_fns["actuation_mapping_fn"] # jit the functions dynamical_matrices_fn = jax.jit(partial(dynamical_matrices_fn)) batched_forward_kinematics = vmap( forward_kinematics_fn, in_axes=(None, None, 0), out_axes=-1 ) - # import matplotlib.pyplot as plt - # plt.plot(chi_ps[0, :], chi_ps[1, :]) - # plt.axis("equal") - # plt.grid(True) - # plt.xlabel("x [m]") - # plt.ylabel("y [m]") - # plt.show() - - # Displaying the image - # window_name = f"Planar PCS with {num_segments} segments" - # img = draw_robot(batched_forward_kinematics, params, q0, video_width, video_height) - # cv2.namedWindow(window_name) - # cv2.imshow(window_name, img) - # cv2.waitKey() - # cv2.destroyWindow(window_name) + # test the actuation mapping function + xi_eq = jnp.array([0.0, 0.0, 1.0])[None].repeat(num_segments, axis=0).flatten() + B_xi = strain_basis + # call the actuation mapping function + A = actuation_mapping_fn( + forward_kinematics_fn, + auxiliary_fns["jacobian_fn"], + params, + B_xi, + xi_eq, + jnp.zeros_like(q0), + ) + print("A =\n", A) x0 = jnp.concatenate([q0, jnp.zeros_like(q0)]) # initial condition - tau = jnp.zeros_like(q0) # torques + u = 1e0 * jnp.array([1.0, 1.0])[None].repeat(num_segments, axis=0).flatten() # tendon tensions + print("u =\n", u) - ode_fn = ode_factory(dynamical_matrices_fn, params, tau) + ode_fn = ode_factory(dynamical_matrices_fn, params, u) term = ODETerm(ode_fn) sol = diffeqsolve( diff --git a/src/jsrm/symbolic_derivation/planar_pcs.py b/src/jsrm/symbolic_derivation/planar_pcs.py index ed9977e..99295f0 100644 --- a/src/jsrm/symbolic_derivation/planar_pcs.py +++ b/src/jsrm/symbolic_derivation/planar_pcs.py @@ -192,7 +192,6 @@ def symbolically_derive_planar_pcs_model( "r": r_syms, "rho": rho_syms, "g": g_syms, - "d": d, }, "state_syms": { "xi": xi_syms, diff --git a/src/jsrm/symbolic_expressions/planar_pcs_ns-1.dill b/src/jsrm/symbolic_expressions/planar_pcs_ns-1.dill index 7cf186f98f1d6cc2b75a76463de2516dc3d88393..ed41fb5840369b691a1a82c2e0509e072daf25ff 100644 GIT binary patch literal 7886 zcmb_hOKeb-iqU+ZtR`w=Dz19c`cPtfzmR;AJHg%ujE zc3Y8|`JNqYN2>4>qBX@{SZ_sl`)@`j*NR?^)Qriu&PG@x!JTfW6Ry_>&DX=|t~o)E zKMZ@_=T!svn?g@@?MjhG>9yQCt(2i;9SgTv-}*cc!=ah&59)(3 zaU?%#I{Lz>=_EY?EhK?vz;`QB-pn~~CN#xSv$p00Ok*7Dk@<|pxLqIgnvJkuSqnS;Ffyep!$Ezu6<)&Q z6)cXT?dVIl8P#I1yV;L|f+;q-EgIxZvDOPWBJ+$Dw$W{E^(%?UZoA&>@Z%^lk6hx% zSSE!jn7S1OQ|3|2&Ec-EtwpAAk-oPnjQ&b{NLf(0Lf@1E^X74nX`{Ot_Gs~HcN2>h zQ}Px&!}e;}>qq9{7hOWU)CsLmcL+99wzM~foyMTq?er^QD{P0ILA|$CA@|q2oo+h> zNTi_O>_EYFHLI3YMJ+)9KEVGC`q1ox3e^o$XmtCaY$YpY)ll=asMIVb2TYDhaP*Q% zEcAmMkWtsmBoRzEJ?R}g7I}lTq!Avij6BO|^@l+r$lZy|7p1N_u49|}p6UaMNTluN@w+~_%GzxiYdN=-py@D7jsin@C5>>V^@)c zzg73?_W{skJKt09(-2^}v|zCA#bA{jSofTBu$n3HJi#&x)($8fm%wPjVp(nlM(D@L zQ1y3bhmVigR8e9FOJ#gwDHE))7RUJLr=g_8prnv7C#tnIm(ZCSEqFCb9hjlGz~3&$ zOX`o#4;Mk(X_i#SFRCTQjR%wt*9D4Isy8_jFZ1?Xwtwhn`=7uMIZo>$g2=c@-oTCl zYZ>9&>AT2uAKLx^5^EnrEo9V@mX-CcLF9jStV`5VS)iy5@+rCxEFByNDGL_TA$s!! z;Qr}cNK#)6cTLP039>}Rd~AD6#n}45=^r4{{{oX3LKa)7{z@bEpGC4?R=6UF^g1!( z5uXek`;Lr##8^ZJdA&t!L+l)dJCQU=YWq}gbTIk2DOJ)>lDHTg&U>3xCAH&cRY`H@ zPVk}1Bob?dp|f&;0fzeWNio5Zd|vPnPX6q2G6WeJNF z<>`ddr(;q$0cv(zs3+%CNq?1CLB~Osjv-FCgmKn_EDhB-a1SDP&RHks2y*%vT5x)8 z4A@7PHAt>MI8IEZlS_ZYEje`J{t6|OUw@t`v?h8x+*y^^Dg8UdU1T7m3{Ep24gEI`~({=A>HLAW3v!Bg^ z!)CLE7N*~PVoCbxqU!Ge{rNE_;Ow|t{fFyHnF^zmD1SeJxwddGPzyPw=SZY0G+55024E-ENm8d+xYtNCyONXCJ4N z5l?L7xK=bdZ9MPBwxfxd^M*c#C1-Vw-(TpM5}POO%j*lY({qVkpbC7=tMdxgU^cP- z-}RD?t@Xkl`3_{zcq4HF)tA^(F2<_sNT?7IKDL3__R&xm^*0?$#@7jEs73vvV~~oB z8C?{{<#X!fa!B>b(J+g!*-goH@_OL`DJCifkC z7GE)i_&*Ek0W5b#;hp*EWL>6LqtmqE$5+&Nq}x7iP26ny-hZ0mL zOAr=210f<)X6PbW#LiJA6VL}IVzU|~)Y6C;IuOW|=yFrIq=`jaDolB@pbt}?ESX|i zu*=!zV58fGb7I;Z$T#~zw^b3g6fP|sOP*~^dACRgusv0iAC zo7)CiA)4hVLAF3T=_gC1Ye|QMJ~>aA!uMOhj^21mk4lDlyu%=p5 literal 7890 zcmcIpO>A99754ja(z>N;UqM7diij+nM}(}zX$&%;fN4tnZd0PPl@OA#pWnNA^M?1Q z_Y=QViR3IuvX}*`Zzb4NC5s5LKw{M$LRO^6REPx<5(`!cL04?p@KeJ1&dj-Y?)$Ob zO1Qc&cV^Ck}t*}<^*KdT; zeRG_eZ-w19I_WMv5>}!@wQ3S6kxO&}6FA zR;Ssxdip*oBtru;*Xx)2VG>?;Sa(`-!@AS|0&J6B>VdY~k@9-Zdp)7a59=FMr(hQE zu_xTJ&X6i@Vd|ZTTyi1iMZQ_?ck7j~w^$8Zy)ZI` zD}#P{s}a78<`oPMquuC>cNpt@x4qqqf}F`$+6{WinEXaJ>_q02Rn}=YZuS-vmF;G^ z-s0vkGS9xu%~&TbUNCbz3TDi6b~clsT&+eX_f7iUrN#9}ngilS?h5@9F_z8q-lR@@ zJM7Zn)%G?9mrTJMYz>-QVYe5VXI^tO`A#deG2JEDOwr0-AG9j{db`zI3>#rHZ1v0C zn~M~Ft=(!jLx4mEdi53ztgAV-sg~3l6yO2=_voSC6;)T-eS+=IE{w{h) zMyx&7pw5NDx!Fm7*H%nl& zV6iSY0w?r~i;T8NQc>DeHozLq__Le6Mpcz{isc}ilKarg zA$EweVi6snJWs*ypDv1o{o=c;W6nv4B_@_*-(x9;*T+t8Ak>dANgN@MEn0u2H;$iG za$HuqB9!zGXd;?;(s1aU8>0CfCPYMgkkz+oc!2#QcQ=w!NyVSCj*cfEmskb;SyC6{ z!|b=*DyRoOw+f2ucLEXRCyCgIWCn|j;iEt20+JwNh)@tD(U5w!9E=m{tE4~(C4L=1 zD1FKqBtr!v3ab_4X-?_Wu_%IoG>0wH(`i-E=V>bBc%7wgNDz)Q&RdA3p_&BlIC5Wh z-ibScoc;n0IJ)XBmvD$zHA7jX|AswdLP_ZK9dJFH28K0c@l{2BQBcMaAP+lQ2#@0n9-i-R^hlONbqQgdE6k4?PTQI#@WJTg9 zD!KEMNL1tK9T)@=-th2roOkq(5EsAwnBLG26int@sYLOXjRS^`#7fo7zS3xy`#5&e zms{o>eP9HcK7V$ar|Bc0QpZ<76kNbpDs7m6;-_Ser7uuo2RSRPA5t%iS<-pfIC=du zNKBf5#PAlo=q$|W2Q*OHM`}Bt5-T>bgW0jHP>%#A0BmhrO+j{o8iJM`27X7L0Eb5= zm`|0xzpbeB0o;CR6dZP&J+v_W>Jw|yzXxgz-2^2JU9y%5ZZ`Q*dN}zR_K{Eynad_o zQzp$dne!k;b}VkenNxVw;D@Ab>MXX0dVyJZKyz=Y%yZ{1_L8Y(%mea~rhVkPmwj6BjF zF*zyi(zWw|dL&zn6it-CiG|DGY~mGQX}Z`Mc(Q)2g|e0g<&c6v^+ ziw>ezT_uKEYpBV?aPYyr_ZbmV+mnU+t7 zQFX-N^d$0}oSqSr*22#|oD$O(5A=@;Oq@zyl{qZ?R-F}Cl5HYw^}xO25bEpq9U-Ao z|B@PF$7H^3&*D{+i~s+EdJNB9Qg~;6I$0O#(`cSH{P?08O}bNczS8Idt}p%ON-S)5 zGfdxc0@6eGn5iBp5uH|Mfvm_ZvcY9;yd!q2Ms8UdQ9}m;=@MJ63zsypYAb~6DIfOMvWMwKcA_ngmy||fh%S}qHqGOd{inBFL zo4Hpv`e7@+i^er$h0bIARf6L_!$nv8m&=Xo_9lmW)-5JD{GRb~EStGfN})v%*Z&mp s8MUV0B`j^wMT&r|>7M{%ToX*y46hM2U195=L;4r=BXB@jc`#V~AD>-on*aa+ diff --git a/src/jsrm/symbolic_expressions/planar_pcs_ns-2.dill b/src/jsrm/symbolic_expressions/planar_pcs_ns-2.dill index 2c1f753e7705f4fd98fd567aff2eb0089bc5ce76..3707e05aac3d6989f46c8b0e9b8a4dba5de7132e 100644 GIT binary patch literal 70395 zcmeHwd#q*Gd6(}z?AS?KdlE)LByIE0ba0Cuk7KG1%Ank|llEM4Q!=fq7*~@S&mHgC z%)^;G6Wf77<7({JXYuc>eY$PF+5AVf*Cv z(-*dPaA5D1;R8GKd+7HQPv7_0rOR9Q(bvD)9F3jtyS#Pk{LXXR zU)-6mo!|MHox{WJ^F8Mu*uhD{eDTu7i(8MMx^nI_TRYFq@1UQb-MV~f=eet2+~I*o zWZ-c3bHhD2xHVt9Y*;P}7CJ}RG5WHMmIGV!y_bheI(8p@os7nza{SukhWT-feTH#9 zjxo>dY+aq-wSDE(m94CVJx`s}@q3>-r^xS7R_KW0c1Ip$Sx63Hm&SasI`U9;V{w{Rn=3 z3@4x3d4A{5Kf@N^clpwj+dCV3=ldSJbe_Jf&G((Wy!FJ+{QZLLiA(37-o7uRymaBz zxr_Yssh#;deuRG}oa90qCof*TaC+rY z%DGDyx9{6Jzja~j;+0dEpT3X8fBe$LOBc4FBqFeV?ji_$cDQS}HXI$!0s%fSz~|^g zN7sjEhtJLTK6YsvBY!@hIXv7Lj?u}_4|n0<6&{X*fOtM1#EI>VHE3q|>|w$PpqCu> zFUWVp_d<#3Mf?;Tdl<*hZ(rTmyRr7{&iwmDQnAt*wgyT2)DxGtwzo<0%cq`{gFEwg z{b)M+iBnIC;Y4Mh+xf)K)ziZT@rdDxjlILCH`a#B!_{QZ%=PD!rVcbT_5AS9vxbJ* zaD6yF+`w2;6SVvllK;5q8B)8ko5O3QwXfwPrVm!YaLqLPQdZB5Kwil|f#KDxt)%h4 zuU33Tbq=ss^*dm`QR<&jWO|M~0KZw%WAZq(y&=~5nptP~HBs%>d_RTh*4EBK9slO= zTf`}dcwHTY_Cj5K`z_Jc9}a&^TtF{3SJ2Cd3bjBr+5SnUsj}w#XDOk*N4z^_pK}*3 zNGjWzzx&7SdzsdTXXbbO_{A$*k5l$|p7I(LG#)0-Hg3n)`JZuVw?I&(4~A0{FB&v0#gkiJ{Y8trw(e2tTq zWgp}!<301s)QR(T;T+T&+zNK)45P-685fS9v-$hv8^j0V@esi7#h)LgMoGrAX;@ z`0V&B!IKS>48S$9;%n-2BKtQmLK(0#LW;E{XF{78Lv0^p<#?NnH^82X#*~Ry5<%(1 z3fZP4vT-gVmR2y+k?Bo@h7E+#?5@AeIvT(1crgwZcm#1PYW`h<9eID880XY5&!}cdXmAFbCbKF`#_Pa)jW9dqa$W{M<8yL$>foH@jq3RK zXb_~pXM}pjf5>2_qw!a9>>u+n@~xTG2>2P35(ANP4w4VPYfth)CezW881W`XOti;( zOJc~$r>-@Bs@v*{|X^X=U`dJKgojoB%y^;rT$AK_4pzB z(-t@LmA&FYD9d}r!R#SspZOiLlV_|zPWp-4UC&01ILG!0{B`xkg`&;o&672J;xuVQ&2BPD*1<%FE)*p9PlIs z^hL7K{->FB?jMuP|d zn1}&fMH}B%?iqqs>2ftSfwX~DwlNG$i^7j3Y0@ax|3tat-+*fk%A+wwGMFyMkJCB& z&Dvq~EUyq!HIAPnByD|NOd1gJG14e)>=|Dq;H#K_-2Uv2b-0Th5lxlzL`sV|a-2PM zorsWu8m=@a-_%sC?j*VYY zHt|ELp0ILEIJpfm4q}%Fl&S?$VU-IZQE;OXBFJzyO>PzXq|zN9fUShKbe-(HI_k8I z+D9S?GB7Ew%EU1cGH>kT9KaeT7iI7K&14mVq3JfpdePx4M4EV2Wx#Zd)x^aA8~ba_ zuO-Z|H;^++A#&>}ft+>Zk>@JV$!v?V2N1Ff4~+jIF;@0&5K?u@){yT&%x1JOL5r`3 z#~~e+liMEGsYNXbBq;ft!AuFYs^qJ-?HgM3V)D*<63ST-$(v#2Q6Q5hUBzl zXEKB*bA!QWK{*qB7C$S5t zrX-_`PZMbzVa*1q<96kr_~z*Lx`iw$b=ZUx(ajEN65#)h{o@NnXxik^G9FFHXu1Cu zx}EXgKr#Q86_Z0eQrAuMB@SwyCg2D@qKm__oW&~{+#dW{!iv-5f7361Tw214+j44D z-ltH`=F@Oq*r<$q*(Z`@=6VN*LY#w198 z`%M}JiK{*1M+rqMybg}vPbk3}`3lWn$It9bp}I-aS+^QGqqDZ&2^j=vJWT)+Nk2fC zrruZ7Y?>XLmJXC!I*>4%md3v>>P40>n8-r~SGQz8OU$WYV4b}rtlB?KybC@|;>Cgs ztM&sbFi|*2US$*@f@(BHAD`5NOV@5#TOFHf0+A5LrufHS5>)|jvEqP%-LgIY6M`In zO~5Bpj3mzZ8+^V|fmv0~rN>90aum9c5ceFTzeDhT1bfU`40~<-s19;6O(St(VI6<3 z6oAZOzM8ei&TMXr!;YgO45u88e?r*D*Hs5x2S*ql@F@&aVB+Y;+W0n-Zb4x?IDVY4 zs7oO7*7kdAlLTdmEC{|6gbZRiMV4g@OH==j!$B&iMj*xm?w0a$u4E4pc?mRe!Q?$+ zfkBJeuM-2N{HPE9q)H(XmAzTjl)@=3cD0JxH;K@+n9bMLb!qV)A~W(IU3EW9dnj-F z1YHAg=^|Y!K$|FQ8}DG6&&>Dy#5vh*$kpV(`$Zl!5ti2u@g(M6Ygn!i@8xh{OjwgT zOlLMefc1D9c#NA>KKJ?ID>x{+8qbg;(E0UrejEzdpy1+3nxK0RAjh*bLXtHlJT`Q*96J6oVbF2J2k;3WU&J@BFJ|T# zmElbyc$D{spzmS~bV^vT*zH8}Sp0~jl>!H=40NbhWlSaopIFC<@mHtJrE^qpYak>r zmh9yaE}j0;^fb#1kMsa;3(Otk6XPx(qtc!_##@~vNn=rL50G`eUuv93; z_+Yy4HpBDd@6aHxr0iPRgM0vtX+Yw^+p<3E#m{FQaxVf8D75!Ar<`;)k^ZCbtV3Q= zUq{T4hB*ELvrf(y>R0xsLjw7UK~JN`vOe{$!0DrcwkqI$#>K zM48l#&M&zUB(#y(q(QnlPz=RUMwp0Tyzr!Hfb>42wr7*!L^H(fEGYqlO^C#fbKuO5 zb6Q}qG~)urX<@Vo0YmsW#~JsNWk&eRcqhJf-9yxIw)->82%hc+<;-|>w0;d;t!M;vHmsK9Ryvvi@=H=QQAEMP1N?n?*TOW6efm2lml}n+3SI4TX!jbFw5A z;sTsUQW`G^>A_78o^$ZahybtDl%2Vqac6`%I13c)RBk7@et2B@1a$8q)Qh>1ND$dg zr2D^{+hWlyQ8jPD()LW^2j3>aL)=ZW6O?40D4swSrWm6#SpXoN5yR#lmY8Ujg3E^F zUT#;2mm5!<7P_>Nr57@vP$Kz>>XfHTrPIY{;BYtYCl_!8p^f%bYYmy(dVlu# zBMA0_!aECEmdZvBAu#R}2_uh(Bw9dcC?6qsO*FFNC<%u0qy_-FJ@C^3X*=`>ZO5=L zYd}EG)Jvox7DFsd!h*z;h8$yKD;A=~OsY3isz%vAsJwE)DB_vxz=5PH7>WVBQ<6Vm zrsOJx)G{)z_Jf8#H55V-hb^o)Zy=d)k7~xm*(LROj0;cEenHq=l0_85sOm}jb1GYB zg*p8p(uGgb*E%ahf3oEAU8-%QDyJBNNBEQ~CY-E5m&!ze6lH+|xlhXq6h5`Uqyjlr z;P`R$Hri!KbA$_2&^8Z;WU&NaZ?L{pbn1KY4IdK`?B7c~sSDp=H^2ZG$^PeSp_mqX z{_w1K`_z+b(w9miSxTtb5s2WANPa4hWXq_=oUf}Bav2pUXc;*k9joMTuuMOS_Eqw= zp7Fgb^h0p)6vao4PVxqGX8d2eEGYmzsO$4@cn@QyNb@j1r~)^eAw^vu-3_)|_FpWN z`gml-MP?eGM_>Ffc>44m>=r~`(p1A(Gf{_`2CRBRh%)^VQKBB<*!GGu)8yc^BiZO? ztY$kYMJ!p+)W_z`+F+A8tT35*&GIIDX1+JwxWejo)NU|7<--)uMg$k;4&zIt>GAFO z^F021RTMwwwv?A{%?#*Bd8S*e+FIq{#a5|KSu{DAdOOj*I)ijJymtDsx*1J}`L~rf z`;GHSLw1m{n%YPKVDE1#BJu0%loh5+tv7&*%Q++c&m^+Baz4eV?UZQ?%-hhcp*>Kt zHST{GJ}e5Uatcr6q{UJf1N05Upr-lSN18hsqB5)>LosY;Xiw%|v=p@e z`1wnxuFzT#UC}xJt_L3S|8{0y!?f<9E41)QH%{J3_f$Q04mVZpY&?WJ29IiUPey;3 z*NM`)efFO;Hx_^4&8EDf1_`WqQvV~f1wEJ`{d|qeL|TeYX&lW{>TrBhkPbr4hIZ3o zfb-dLw>E+#6dW!`C*;26XIC^qie3%}(=wr7NadIK7mfp zzKXM-n|tsK%Mt!8L0^9+{y$Re;XpS zu(+LEuHmev8jDC8>ILfr_sRwvYKc!8%7Ox~1vy!Ho!*&!9iI8Y+I>XdnGv@7 zt0dG!4=Wz=*1YK zSo9ft84Q6Gz8VyP8EKLOb83r~GcTE zlG*7aSt?!ncC<8Z3_i)3ZGA2W-gJxzoSoC>dR`xyWE@ElG5D1EAT;J{529l3sZw5) z^ux%wawk{I7KUs+PI0bRQT2CITwN&x3v@9A(YX>2VTN#B<6Qhgw(_4pn9S1APUW3{Z(plsKi{Zek{Rw8=ml zLk>euqBVw2!NY$PBiP@yE=;|mcr3RU< z!-W}k5f>fbFhMGH)S=SB0!AVi4K%UH zqHe(U+DTr|#11OY@A1Hy<9||=!O3NEGwRuj{aWyZjcQy;<|)dRL(nSZQ5_OzSfs&NUu1bT+piv) z+Er(#MTrFHlnB;EO{`PU)}K%ZdXQ3yt}L?>r;Gv0$*oFsY!aDTJ;`Dd1E@;OcYg9N z@`tq4=JjN{^Mp)Cew=P2m->|DIKL|9H5V)*wWr*VPMw#wSjwd>YNhLSEf{`gQWlf2 zXf#Z}8IF$P3C&GmrDY~H#LlYoAmYheZ#gREr8i4SUS5N<6i9h7j(~JhBGyeJjdk)L z8_+v`+4M(N6cFdI%Yl?vMHR9{RNyPKz%CtG1jI$nV#e4bR})`ABJp{;NcU4!6HdK$ zsJJ^LVskM9ch5$Rpxeh-ECniE4$dXeGDM=-&7K%2w@c4FxE;k2y__qG%2EVb443G+ zTDVv~OntN(H1!gBwWM~cqG1+{AgEXpC9b*OqAGGLa5pQsZVpMFQhQhzE@>o+4Q!5t zYrW2LuKbpf8(S{%x$=U8TCOQyWCtnKRKr6HL>x{mx0+)*V3@BB+~yU%S}kH{f{;X9 zl^UnSkXpxKK#f@J6|6|Bh#=Vobd6={F0o3uI!rpxh4}$YW@>>tmliP4X%`l8Vai#R zuGdha3%~Ur*WdN8ugLqJa?Qqxw3j(6u2o96U6fZ_j9YmxiL@zu`6E_L2fNqwE!e}K zMCpnWj{ofK6mwaiH{%)r#$it_P4foj6ULS#@5)d<=)%}vaFA?^Sy-B0xHOs%x|~c? zn?%hOFFKr1xQvR!d0W5mM4q?(e{nIu|3Ob|PQmONvUGpSKP4Z;98Rw{q16Psb~&vn zXv*uPb4qwR1un-%S=Ah35;pA}7@7pMJIN>k1Hid)$x;MB><2wyCuiFHyp%x!A_`#1 z1%R^D0RT%bx|RcQkyi>G*}-H8$aK4A?Py%wE~b?t*QJ~lpqiTN%2jKlW6{K6?KmAM z^JSkM7n@3&pGkMOHa4eZDOc9KiKXP$vwF!Sv4{SJO6fsf#iD5aa?{W;09rI_UCLnW_%LKLpmw5&e9Q9~sANe-SA+Dcm5E!(swNMB^084c)>v(F) z(X}8QC=+D0%$}fXk#1|}rj%l84m+w&a@7`H2CCL!2P(-`yITj!A|P;5g^D~vm(xoR z5=*smbRvWtJ^t-?G}c|BY~=G;L^51PR~A>5i+nvy=pAMF#1*o)NTz_tj@SEAuv*}I zMgy81FDdGQMF3=#+>AZ66I5#O*oEC44l;LRTW5ENgJQo7kxpD86u;wlI7X~chhg4q zqjlAT*ElU!tGJAE5*812)>GH3PDXNYb;at!rSF6ubSP*Dye2)DGoyQ7g8S839eJi7LjPQ7!Mt0SnBvSb!XtxjmU$d@OU z3(9iGF&*&kj>e9&J0ggz(Ax-!!YBpaYG#1dGqo^TH**qTDbg+P)_JMprC`^=Ha#`+ zyHZ0*Uid9Om_h5j`^G?zWvv^#=n6rtJk1o;fvGoCbut!RZH`r!$t5zntQ3vrK4Kb> zTs>Q(ql+Xu?0`V|EE?Mc>ryG76|*Er?I0b-SWm&!eE~>(?M_n`P>DL*PM*viLur=| zNKVNiNXNCfgR1qWv^s!|oVvi(qHB^CfW+4%brmiTcABbU0Z62TE=bgVcj^(DRi1B( z>HVrJK${MMaF?IezSGI2&Ig)OMmBNuYFl>bIUi`s7tzc0B{f5N(%i}4((+f6 zwJro#%oQdHK3k+*6R_FIbRI=YIBgV>2b67YS$aAhMN&{mtyz9@FJ7CmOHGvX8(zxRaW*#rW4qL>E(gpOiGVxm;kGthtM0gg+OV4hItdowSxv!AX9296CDqBG z3+v}|0bNydRH&6_Cvr`?3(lV5z#}_#=6B4w18-7~veIOVXRvTko1&G@+P1yF7-n=& zRievZWP(>-s1QX7h1SwrtK=qp^;O7Tgq-|&wT({z)V5FsL==eqU^3LPy=~u>GAKYq z0W3AERRjQ*>XxcQ0RZL&xjwb6P64Y+x;2?)pz8Lb1C_QQ>keTVsHWzkuugK-+UQs` zaacP}2g-ao(0s>KQuroIf!Y^Dw=__4tZoD=VzRZhn%L`atw#Q*JO+ZP!D1lt!d~Tx zfczF>7*COrk^bVM-0>&KWI*i%5s~ve1}f=zh-Y*)$YN90yb=nst<@9I^`H8C9`WWU zv<{Vrzje1dni5NIYdg-e1+o+Jl>e75kQVMCdr<_TxHOY^3!u)9oxnzMX{PRo0szcs zb2aM3iRIr?=)OkVHM_-9T}awi&jqp;qyuGwv@H@Pm;@D-Rqpkp!VrZ`c^zCtD;3Dv zq8$Szs9J{|s3ceIZXGC#fOf)+3@j7FNm&(3wF+b+)Gv^c7JV@(owgH4%HpbWQGkb4 zwG$q`dW7sPk}2S^#gKmV2|`@-t)1{-$|)kgDzaVf~<$z2u>( zRzY0%NpXYZh5mJ83CWz(Z_=XQ09v@RPW18?u=FSgNH(mi`2$rh! zXj{Fcx0`JHXq9?GvK3R&t_z-h6XE0wOb@{VBpfo9rZ@f4=QaN zDZ}24P}I+GpzDAaqm^rdq+ecUG%8aux)TikC1r}cI~rAIS={~|Dt2o^#|>aVdy!Kq zW|U450OdQI3az(6>O#};1Li4k!-XKta@1z5ZYUNhD@{vk-|m7^*yrRq*iG2uH|YHr zG~37xDF1-PwV4fH4x{Jvsb;w{vX9%Mb@C&y%3yaEDzOb827ehJA5eg|gJ z&FT(xv|6(}|FK8Xx3eME*DfhHd-Ooi}{4z@IV$XDdAvSb7SS z0LU3B$v->UM78T=Q`%_m-P>EcO!nqLxs*_xcL}t0LhoFKGGnzAxVxi+gx=)4c*~eb z*}?7KC}b(R(=l+;HCvbpoBQuG>5}<^#5%ckQs|}w45m5l$2ds3wz+ybaJj&#Y%u)70sFwwHP8--I^+sv0NyrO)4=6RAI)2XB1xDv9Ub-163nM z7duiUOSx7P-QvewbDpuWDrbO>%h0l8bM>3>8Q1ZZJCYUs2vyW;Tg0HGEKg=#sLGs= zO;7D;XF|*B2~Hi2ZKz4SWJ#6QPTcjYjte)qDf%d3E!|VQP^H6e*MH65|Fhfm%e&UB z_BP9-WN3@#T&y}}MCzLrD@(WVSbPwF%NA$da4%n%iGF|C3(|HoN%dS*BKUUi}TqMl>t3MlB_|$O*K9v}e3idoI}i zqfMyY?Z4F)VS3Zwz;PS1Uwz$Xn)8moe^1`zh~}>2ALszku#_zBV1d`k?$d3%&a9<) z@Z_>NAIs)>-6vQySAx;DU9nHdxxb{^nN}3Mu4cx@o!No&2kMeEjoi(@6rb(kf)t^_XZf$qOnt)P-j&V9Wj~yjWm$jQqX~EKSC+ifi zpfA>YBTEKL&mG);=~*Ww9RY*o_|9m#sekFli)Z4cg}Pxpqbz3rnG4@Ska}gc3*uK+ zGnJ0($b!O36}h!w-lgJKzb!oN|0$!5&iWSjmBx?0OSYe3S~FPGh!4?qbO0DXL>Ry% z?=IFmNt4_0n%gYXMZfmIdS|F}3vIa7-^49X{t!1$A=DflZuNIY!;zN7-C6XDG6lPx z*{{eeQC9bxC$)EXV#sfv)QzSB#{40~>uT5qL0B}iLe!E-yfzt&K|rNWB%OHm+FJ?G zYprJkZ0sWy&v|Y;J4Z%gpw3wMK`~Z6 zCx|<@l9tk0XZc(g+H7S-d5AKIEyXRJasaSWumGTTU)@C-?>>m4v)Dm*yR+bm@!jq$ zFFG!LM3R`7QymEa~iPi`cA8u=nSz4TN<|mlAjF~g0zeqKMAm@L-x_YQp#9F zMT}cF{0l*>O~@{^0HpIIz@UjDCBq;1Rt;FFRGYe93}?n^Ph%Fq)h$z(9!eFj$Twp* zqUWA6*38=r;Og|<;WK;SYI<>SwO%?NpMxaGHeI-+&ss0b!PR={z)1i!{HzO?{H(Tn z7cR-y457&GD8IaGmJ(hW5*?BWH08SxG)fP0i}7swq&X!PLvSXaLa1j(orri!zA392 z?h&xK=Mun?P9dG#8eS*!ZYZI+CzlEQTY+Q!aBo0cK(+b1VC<*_R(sn*Fx(VgQ&S^G zT1_64ob49jlgG8;W5XgQ;yJBSmo& zaNWTmHK+`p)p&@)ugybbGsGa}ms_f#W&Ey67rb}qcEUrREV@u{9cRN%zLVO56M0$Y zRrGuYSre3}CX@?Fksf#V&L-FIqFd29YN@>_mV(upuJx*l*oiX6-EO*?>;#Q`mkRod zW&(*^<|te+h?+xXTpMJiHtr~bxh;M~1rf#>pzMY^ItGM9ZBx%( zLP^)4C`&toTo$=_A!@o+$CJRwsx$T*@%+VT!>zg|?r^KFIXc{`>x_n56?}I>u(MSc znE;la+?Crw&d#*7S4d3M4&HGj0jV9l3xY7JC+Z79>RFqDM9w2Aq3PW*2qJIAl+fK; zYZ`DhDk)l*wq(i(K_clS#Y8CaXYq{~*Lr?{XKw>;Q;znhBALsQ-eDlT2!;2=quY!9 zM{yfjM{i+HjVv)^%T97rQ^7t< zmI{t;`r|4T=BAc8++O&K3T_5=BsAyD!lj`Ab7_2EmbzXiGu%7kJ?f*|iJ?0J%H}M} zFl8Kcve)j$<8m91-O{()+LmtJ$sZCSb>&|OQdj=2N&d98wxt5w9m+C!p%?DbLN9J| zcU$5xFYLC&%e-rNwF#Gl!Z@lqy+FPVcG|IjpYEX)+ zT-V0mrYB%hgSB^hM{9lo?xcZh0;HFz2te7x&Av+92JvLlZ0s*ZxlYrz@7^Pqa+m^& z3ml>)VJyI7?R_Nxuss2ozw?y)cY0^`bujiV50e1sb3D$X`f&=6$P7ycw9d!&qhF&e(E#*YTU(WezrP z!^G2z;fkB>dDS(lp?0_S;wGB5I45bY{RAu8OAS{1JFPqbnu`h;+L09HjZg|3$LN~4-g2P%4c{9<#@N@4!=ja6+;7|&-tG|QG^~XFX^Lrx(FR*W$Gz3SMY88t5 z@0Klo{0y)nn>egkq>-O^@CjDzY3m;nrYr&qzBwv-^eA;cN>rdQzTSBGPUb#pm8fGv z;0r?u^gC}sd1?1N%@HKp#u7i9SABw<#@nuT+jcIPEI@r3tJBg5Q}rY&@EWC3kR%5y ziBpWQTNN8&;_0a-r1DcreQVs^VE(YTCyW?*S45S(Tb?ZZ9&wfunXwej2uSq^YkoSP zTdr?;t22#ycBeLx?KCROu~EDsIas*_&?I`$;5XveL1Z~`pm!0@wVaPB!EI&!7X#=oVF%``dFvUOvVRk7R-Ut8_r+NzA zsqQxOEGcyS2rSie+XVpD6J((~0Fc3TM5XRqLc69)0Xjls`bD#gQs&NhTW4sq6{X-} zepu4#O5a>S=s<0DmjXKJzF;!Q-W^cbsW<=aDAXqDD_Kpp@{9Rzlvp zr(J)(v<20>k}?4p0!Pe7GuDm`!p{qO!FF3tj1R~3+Q|PqvIJ;?$V7Rck|;#wTHFtAHKjZ5M_Svrz&?yp45Oxy$-LMG=2!4%*pcs9`a&6yFE?f?a$X)vna^igC6yY-L?OdX)kelK(T7)=4T9ZDLQ__b z=%a_T<);GIBc@8jh>}W;V=7L~$3&bYumvh;P;*Rgb_fH63DD`+Z|!;l(uah>XNCweVDceU{u}w7&pIyd2PBbx;viUP3Z!ly?zGaw%UDsvO*Q;`Id>m8wHbVfcq zJ>ZFv{qNUdhRGc=KJ^f-j%H*6hy_!ZirJmG(@KoFLyXSU>e~NYmBeKkW3uN>{RV=S zk;n(A<5sWD{!Ly8b9yfT-44wZ97jBL(7fp!!{}$4=It6?Q;Oz)^SwWH@^S5dxb;c> z4|rF{IwWCmeR%KqvsB(#ZVlpNN%Z?EnOeXiX;_FM$r^hc783o((qI5u)MKzB!8Rc9 zaD_2;i6yqV9=4zBrBw;HYL#&idi=Cc5(^KTotPN==_(nyxhk05IejNbmhZxehEZlv z?2b#8$VRy0I{7btx-dX)KAt&+8jcI6fIY&46D%}wJUcn`fYS(oFpeyaN>K!pph`i=tQG@69gm%27lHwne30dST6;GoI2b@l3f{|g_U|{(tUI9IeUG7?vH!!PaOK$gSV{Qf`4DwJGOYo<hPhx#a;CI@{{+RxpZaw9(w!8rSp4>TOW`&&+ZKepQ5pgLszy>pWoZv z`SRXk<^0}f_Kpp=E)Jak@E&dw=8KmuUfh1{^wo1ezP-1*xQ#x3e*4O$z1?eH-s6GC zW#I7c-Qg}=++M6)Q7l&k3*95^7=1-X%Z2U5!7Ib1IQ9^|9gW7Na{S6;iuo~&z0Nov z!gh+$Z~p*3 zei}EQ*n4*G&tGRhICSOG<2!pB2N#FVTslv0Ru+dIy0U$FZ}9;^cKOozCwK12C@)<& zeeNQEd}43$&L7~92`4$$2KnvAEsO6GX!g+4XV2~}4t|3EK1)p6I~T5z;~xAp{UxtG zxp=n$x_s&J?JIQiM=w2&nj%a&K$z+9Si~#3P2Q8wZDL8!N+)4^Jk0X0E@KG_~5$)N{kX&Kepv zha1ECu#2&zCTRI3l7C(F45`i7t>IPD+BfnQ)dwqJxUQOgA**MTKwi!tf#H>`t)%gP z9Ig0r)H%RjiSGgPTB(0Vk?9%o0DQfq$K-KoXG5&>4Ykhji=x`?#bFB5?d`p#I{ua6 zSBX;)vAQ|}?YX-8gSSOjzcc(EaY1?sC`B_=6Asb%CECd@ki;76YjJoz<+Se>2Tuv; z+=UC0&Gr`W{gl2hGu!ag;TuEcaf9H z;lFb2@_D2YMCgs-1jDTV!JBWs2{88H7`)H(2PB3U;ydJl7{na@Uy05|35A><0soaW za!p3!Ge$vQ@!m5$c6Si`&*JLnj?iOVeCMgfp$Ad6IiH6INl5-PTn`_k_ZqfJdo5zV z&Us4{5Ym{^Y1MD9B`4IlxM}No= zBFC$&n01xM*S`ig-8fKMb8~n->Ue8-IlkMa5Y+UrOM$@YK^a9}I6)a~ zH5&Gn1hh6V3;3*qqZ^clWik;;?}%X;SG+@cGu%%~;E(Khn~A}#AyeJd>_tZ18eWeS zDZ>t549^fe*)YifTq7&K8huV=e+?rd1A0eDv6keHYZGOt?PIJQYm=}G_M&J^nRq1; zB7ImP+mu9B&PBx13Tiqsy^7GVfiRlg_2*ef;cJE$@#Yv#ck^koz;>#YBFv$jU0979POv@*HnlQ zelRO2mH_F!!tgAJXfG57ZcSAJl5LIhOwql@Az3(1Vk!Nu{~uFN!>6K#fiT+V!N?a| zM~e`6l9Kwo8iekDCJqgk*#zJi1BGL*ZaSFS+3N#KJ?t`2_>oANt_Nzy8M$(?abP_9 z$I@sJ0RR(G0BV=HRc}n(17zU!q}=eg;9BP+kIEFu zpt=kX&^`Lh+F|r8uMmoA3_nFk+WNYZG$7(*q*B^A5S}OCYnYVW`t+W*y~`XCZIK5 z>jVYem+6g;O(kFXG1Crk%X0>>2LOu=HJmL({q)a8ghq^qb1Xsw;qw^&H|Y5#a1}Y>BA|Yr(Lc^N|1M-=}+c9+@Nl#>%@G=iY~`B;%N{_;DqIc!H8ECpa}YA>l&vA}ftbw* z4Q5-hEVvX$<>dC;gkHy{*wU9>yr}h5s!oo7w(*#Aehs&TQYHA>lRpL-HZB z+4-ktIb=3RZh}Bsi(-Q@ju~YPvKhvTv{o}#r!{7%nt|yABZWUmK(m_L!UxFT^`ho0 z&XAmz^bY5%@ZmTXSL+n@O;Ss=DQ6tKS5;DZ{Uia?vK$x4I`P~Ra{~KKpq%+z9EYx& z#4}no0?dgAa>fY{5u$L2&|7|pOqFTEFBy_C-2f5b!tW4*@OvrhmzUrnWfac-l2Ifl zu?s{^Nk$1@Akw(Pnsrjgosoayo1=RyOG+I&1x0kTL#hP$f8%iY3K1GNIkdJ%Q#9HG zxQT8j{3j^pKeJ+Th&$@ynUJ~@e1DZ^3b=xosN-1pAUTmNZ*Yh3J%kvy!~cqJ?6|b9 z<#*=Pt9&3rL3q4KedCB;+N*@y*e#N9>U~EliH;tBA~~7g+D_y!o5uB*uY^J zzAy3)dZt}SGfG)VQT=VKNPry(8-%45XGg+cAgoanG^q`rKa^s2i$0cXrJHm`S3^NZ z;KHK>BT;oPAsc&QQNL-Bi6ew4uP8eLeNCy#UCuU{iolN1W#ep9^CtDV|`uuOZE`+{A z)8A1Vv(6g>418mFKgS4vz%mwBy$@kgp9Y@f{+923X?O`2MPp$g*P#1r>3%q0I15D_ zw^xGhy8#&jjbK#Q!&Q14E{{IrK5bZ9GHG1+Jl%}e(!HdOL%Q`49tG$PQajcfKtqPc z8qVW*0hf6QH)wPWmat3j-h2Zas-VDgnR0FIul?`bEmu00W zs4T|=R4(4nerb3Eu0A#~4MeJja>zvMrJ)OsPonizL;d@)JgnoWM~?0SwzzA&C$KI? z_6fi0Ay;8)T>w1Ybzo6F+4x|(?>578;n!)9RZ@1X>_IjFLK={Gu(qu4TJiH;gWQV1 z0}AbZ!zd@+O{9OEq((;ZL4EdwiMREP_39!Y7uAQUfdS@4{LUql30rKar6?Too2%)wcpV54%q z!S%!9$~T~U3lY7T8;Jyw-9)S#-~){wfb^=FIUgJ3TxytAOCr7ds>fpMEi7`ZU}`B`|DBpO+m zlmtU5GX?;;J@C^7X*=Qv+QDI8jsXEVQ!kN*XbiD12@4WW8gi129kCECW-@v+l~k1d zBav4w7)3mD9av4Of}t3|Iwko7W=f_~NG&7d#(q%Ir-ni(;x~M<@%^BmSA^xm#r8@r@j~O@KO=M z{@uh=bm3j*83SM>{hyCdF)dOQ4S&YEed@_$(w9miSxTtrm8@pDBKgUV*)pyiHD-K0 zIw6-)fr6Hip+b#PqW(#m9bvs5}WTx;O`r-${)3w{!Er`6NsfMvOMICCI z9AkPVhlm;r6-LygrU8rgifFyOl8vs$LR4#mP_&a$#F7P7eQeIG4T6?hD|CP?dChc_ zJ+(NPz5v7OHq>?*U*tm-4{E7_1D)^!X*&Ed{CN(4{;eoJCk^RE#+q%+t4Kc)5s6>lAdl0v z#zFX$|*dNvu0PPG$O=Y za9tXh!ohP;>jHvgi|9gFi8_M!mQxk%mfMM~Sn@|K@)SS`F^fwH#pM8KOM^h8$rxU0 zr!ksXo!0!u`o!n(eegea$KSxaBgAo08K_y|V+rbQlpH`TdpO?V#n?EFYy;6hNq9d9 zCGd@p&^ru+x)v*+Z0=-;%CLS6#ju^BJ(+vaQqbYG^OsIvrL`bBVY~RQ58rS9?XACx zY2E!-Y2lN;mU#z#H|xwfd^u}x<9>Ye@I-9xDR3X;b)vLxpZzDzjic{!o~TxnSpBMcW8FT;HvH6iv|PSR>@I5(=G`~$hlmOPRMP`ALp$ih;2Nm zyNk>89Go`_3RsSg%dF!Yb+5;3*NV(a%L}=B!-ya;%55&)W%3FJAyv?m!%?xPyHB*J-~jsNi?*ZGXN4ku+jVvTMn`w*h8-hzN9 zMT12VM9{S;q{@}7szMjOg#b(=8=j$L3VWuDrB@_s#w!e`g)gIIrcg*t39%=sNChd# za*oVU$~s!zzJHB7-olNHw-B?1BN?p;Tij>&9BC~&Qn-L@^*3UoQj4nxK+Nx_d&lUD zyDJ~2@0V_@NS+BN>8d1{Q!&9%9n!&cZh>yEe;tXb2EdN1eEV(VFqT3vewSqh2S{c^|y4$Tw8QKIvI~oX;@0v z^1#_SZLa6_5jQ$~@GL9;Z0dt3(BKQ(Rps z9Sd|ZIMF$a2RB2Qu5q~(E{56^9Vy%l;dL(PuE?5oYgM%kjwvGp)4-b=@WIEPSf{rR z^*Rn_$zFrkQP`zEX1?h^p6pfI<&?Ncz1_r&@@SKRHijIAoJ4C3okE>5@rVsB3^#{q z#7-ZRuU1TfsGrqRQ+Ib!t6puvW&xn%Ik!e#=^((?1l(d00*`-iz<#xqf6t~OH}46J zY|@O3C1^$lsnK-&LuDB&!%5ZMnT=!7`6w!-I+?D+g&MZ-D~F4YZ?v4tVWfivjCd|O zj~5lEGEE(}kr&Tdm6uW}iUJLZ+mv{bx6X4=K3$5a8?e20hSxK(gUa)JJhf(>`!rc~ z#6Ft|NOA z`otbG+i1TQkXXSZ>~?BEJlJ);s6)OAjThUYvU&rTuuDMzi8^~!lt_R^i4Zs{5gLKC zMyKdON+p`IOl*vobVZ4!qM=r+XoV^4k#trXAgV-t=da#L{*V^tO3R4uEFmmpdHSSn zWJ;+poZ{u>X}c%^H5ve`je0c1<)zVg@`69b1H^s!d6>MAXXfK}7+?3)jB>-KP zCuJsql&5$mf*31KneTuxCTC;y4YjN_M0zh^}OF%CNoq&yW> z$l_5!2$V;eCW6#<>4*<1pd zSiC8pZ#m>Ws^P8$JPt>eo9CDY80KpOw|PacR*TqiHECMVrqmYff^imgsdWqiYSQxA zj;2Vzo>mb-#NymKmjbUXn@^38tZI&^)j3%}hi(@ZaiLN=W?TZBHI(SWZ|%!0oA=LG z*tH(XnXauy&jdl=`yGj>}CoTCveSq&X!^xw7U>EG1t(tM{uCdq{y)N{{gRSG`!% z?YPj3-c*}5vn3p_E#;=6V*vD8+Paj%*zsZLUxA7$BbMkuS{~ar_q@_FIuEj_OX>lT zHLg&~lI3!*BRfvc(du!oN?ZA@PFEejLPD$qv;%;I*vuXtO8`)HP4Xp80j1EL91R)j zG}r+EOLZ>m_*=@hwICfR6=Ys!Pf)c;x3t%rQi`cL?1XKStG4JAs9J{|s3ceIZXGC% z0LMuql;?9Y^>gR#tsI>Qahv+n?;Nk6bnI;NSn0H#IEpOJl}oBX)5xAOY~pg+^Yf(1 z&9%k4+Q_9~Uf^3s+j{Z1E$V?P2eL|T#va-!Dm8dmM0t33e}_ZQH^Z^yt17OKc#cJB zbbp7VnFHGzYbM(^i28K)-%j>}UH-6kI09Cv-7s&q(Q(a!*EnmrI+?}nOw=lz23K3F z3%8}KfX$xD%SL%m&S^_J^vt^Wtc+{cw`b|az%}*axuldX>NqB63bVAmzoU^CS839e zc3$dkuhSJ7%heIo$qX{%rB=taoaf6U%QR*2c9G)C8aTylHM z0Gtp>+MKO4wkp>{Z!~va<>(G>jP?#mECs36zXZfQWS%}by$ z*YRowd)TZc4GYa~)`a#N#2h~Dr7sW8v_s8C)`wf(sC#i__-tC#hfDP8v=Ytf!5R6Q zw)Ei=y_(NX0WST$X0;1yJ+F4!4ruA`vluS@y~(iW17^}y`%Wj98Xss%>Dk24t8Lk# zXMCV3pGPm(m(&c|#eim-l$O7mtaUQJOz-AQ5^T0exyE6$lj$sqlyKT8JP#<_+_Lm^ zI*Ozqms&@!+Kbm_>{8=)#BJfD*Z8u*Z?p|BW$QSbnSilfYF4KKvqe1Mj(WJQB{^%b zjvGW9nxfYP_kslASxs%yEm$*4s*^z{v6u5Xhi+7}SEwV;M&z1w7o0uAfY%yZ0!~sJ z@FsO%Y1S`O1i?XdidH&n+xGS_OfXMXqT??z!7C3ccu_*3HF;~5+=P!Mx$Jq!$)D$K zd;*}hMO1)Cf!KEz3_Hft?YmNj2;flwOU-KK0f41??XyDx0JR6Ttxf?O4`bJ4ra;y0 zMF%QvLDn6@6sV@=ys%Dk)!OJ-G;vrvP6tYT+0lH*R8sgxOM#uQG$2jzLB5(d(Ozq7 zHL=&p$|}dBo2_p>?P{{O$XyqbjlVwzlI;Es&j%r~E&;Kx(+V?0FG{ z;?hjwEr2>Zb^;s4rJ1@t3II@_wX>zpTFv|~)3s}Mi>11dw5^^CWGzSsDixBZMZ5%) zpuDomy}tKdH!p0;b8sH5R3K}Mb`n{Fs&&|bN^;fi)`8LpXeUh1z%qRvl~u7+t3W0~ z{Q~J}(H4`*ks4W?D;EX0TU9#=$5s!Qy=5{vJa)XkD06kZbZl)Ib!c`W$$Zqr7OfFv zbqXZqqqu8@WC_QYbWB11|NHmZxj-gIHM&g&(!>^Bt6DHRgJsR5-q8F$yQ(k^ByGu; zbmVvKXv>eDuXq+Gs__*Gdg^Pj3E0e5N99?N&8*|c^ilGBJn~7o(Zk}+XB6dUsFxto zdwnW+1WV(IHXRR?hg&rsNY#m!yR+XiQB{d*o*txXPhz*Smpn9;%qyCur?9e3uglN0 zqW(mm=igW&r2NyHAeOJ=9bO0ryi{k5X~ZX)J3mNWc?FOnm{jS}wpx#g`t=q$iie0Q zXI3*COQG70(Lwk5R^}^nPRG;n1$to2=;3=V!(Q5*^mMY0J2~(sO%$%=KCiE4y{$Em ze+S;AB&9#~Yk`!KkII+exG2wmq9wd0fA5XqZ;;QI=IB_`QRl3l$5BV!)pficodvR9 z65r+A^L#)elQASRO9|O=A%WwEEhgah6J;mSZu^iAPG+-0&x=0zjAQe8NE(jg2-5K$ z^fTYa1HIPpY#427k||0%&hZ_G`HmZay8Wp~u1*LZN0A*gz7E=~cr_UwSu3K=hlj&g zro%_ap3TUZhXyN1wOb@{VBpfo9rb$857I2#{uZhJHbPN9J+0XvVOfs0_??lZxaFX` z^xhG7)V{sql8as$dAVDZ?@u1Af-q-kX+ z-iuE)D}EV>f5%7cPaUb~kzvswr^PEEz)OYDAmtw+ImN$$S#+~;2ijY$*`NRDS6)nJ zujW6hV4F%B6QKBrjA#UkvKdn${6+_7U-pof@E69%OF5W41xf(qIF#g{8EvBKIoXsp zntS*9t6eI4bD$h06z83QmcyYK-8qF){92fB(jAS=Fco+4D`QJ>%agwZx2EWr+q`PQ z8B}edDr_z|Z%UWU2NLV#)=8mT2GE)2vmu2sb> zK&y&oO=+6Ssa%msUoMo?CKVq9s!-#?GYUWL*jOI^fvS-g-BWTqdXidt=p7?hOd$*fCNne(ygsU7W1XsVvz)ZW;Jn%s=j z80$E3*RMJ*+@hf2pB*@Tr*^4IhyAYq3-SP zCcF@TQ;T!la4&x@6aD_w18Mu2q=QaOtSHt}Ugw#yvFGq@2<%L6I?bfB&!ij4xn9qMuj~Di--9(@t%ecu2Gaa(yW0Ma zIXVvNz-ztB>gd2b>TFJiJx{_Ji(&_>bcBu%0`hJu&yDS#jtLts=qRp(q`B!+3J@_k zRg`xfS4n`p1J~SqY-z5L-{HSwxYo8^9Sj{L{{BpvYdfIN%?pk90NcZ*AOMnsiuSdhX!1OV8zG*~9HD z$M@7>mg6ThUh|f3QV|1qPD&iW?z zmBPk*WG4h+r6J`t;)8dN4gld3gaJ(Q?qa=@H2OMTdwZyuMRdLt>u&Wo(Yjmxozdtf zWqy%Ds5x35TQPk!9BJz4&Z3_e&Fn*Jzap=AE3|g=r1tJk4B5>Sb=d5no!lgmR9C|; z2;#Mw6}*;2s5BY#K|rOsMYBwMvjDC3>RC(75I&RGN16!&xzxDuj3U0e=_s9FC>+Ui zyZx0dJsv=-RLMYD+uZ1kPDZ=CZHSD*K<&Wz4&6nHy`wytHN5u-iq7(-*ZuCyIWf*@ za&obdq}t5GUWOEqIyZJIi!?#a6xXJrDK3|xZf-+Xtjo|+5UnZD8cMjn;FVJ~nL8wW zCnf)ChwKn7_vPj8{++sMFHO}k+G2jIHW9+Vd9QS}MG17gFd8)k`buzhhU(O4pMfWl zUV+E=pGH5E*=Yyx)w(U$>pa{T?LPys$&P^3!5HJWMa6#xV3&%_StrZFdUCq84r2iE zp8?opOh7#Gc`D$XlaxhdernDQY(33~<%i}pV52K|4Oo64RjI=cNPaey3(_*K|0KYs z4%tT?OLg2@j9V9)r6A1`;K?G!)qI|-M!VsnrDC`cFBHQ4CkqRMX&^HPm+o zR5xCc_Z*OI<-wV!jMb&F!>56prO)hv?dhG*XOgehOUL7LkfaGwwtE*Y(W~_`4X)Np z2TlT@;b&dA4qSRdvY1azay}(pQzEfVC<;0 zb5wgZgVd5Tc;@logbLi z$)~KU)U^`8(`8;oPhya>r{d)sW0N;@Sqz1s&tHDeI@yvbBh2G15lULn z_(qJ!dcMlDw*lLH#qW@IOkS4sj&j^XD6M5Pnz}&qW#|YOjbje94#VpkvwmTKqhhk= zC5C31S3E+UNHxx1{wTgi*3p}rQ@u#qWXp=mQZvz1uzhxcp4;if91BuvN_cLW!<-~2 z!8*Pq3RZ39THggh`cRvI5)_~=jrUWj$8$1tZX3U&9?VV*eIp?0bHBC4SFR?ggB!y# zemhR?1Y#WRXY*>Yr|Tm7t?hnm%kx;9C*I1u&`F$b`u3dKJ-yaWH*C=&*>{+*$$ok4 zx~JEgoVaDDIl6Ot&8un6vg4hVq2-AEme_S`&);M+su7Wk?XoOAukHM!*AlA&DrJ$C zEU^A{;?_Cfqta{~E=9RPP1Ji1$VoF*KykpqTN1_sJl5W40f6lhz~Wt}&A&(X)_)Dm zf22~;S(*`uo?iJA-yRayor6RkcjA~)<`%rsF#8!R0?Q3lR7z3OG;@S|=qGz?fkDQ14Lj1?9%FR3_vmH6#lcS`uI?`;D z=aIMrR`q`fD{&a!kKfwE7qJF{OcR#LLOe_}%Ddh__L@$v%O{4G$GBVL9WNptV31LS zf9Fln-YyewIVSE*3@Q_M#_lh3Y;=O9Ugl_+czhVH__94Wjb7O3x3yy=b{bvU;vA*9 z_7kjVBi5M~b&vy~xj5mqSj3I9Khl}lfTLla=BUZ3qy0Dx8zlX!b-jVFz2 zG9BHfpH^rtQF{vBl+z|SjNw?OG>Ndv{Y(`-xJ)Apn-_b`5!+~ zgdQy~>7#hM2wgRGo7B{Mye8qsZrSq3&j4#)<7TYz{g}Rp0j$`s%7#$NB4~_kb5!)` zUg~_9s6b(S-gtT^bDUbG{cc5E%XKf$ryu1A5^ZCNpUn$v(FbxGYrEQQ+qq!0F#R>G zPD>+9)sv{eYa*3`Bso|~oN^qyRk0B!o}Ow#DnF&vH^<#B^M^AI!ibVLMKqE(S-16wyhjbScUk?P*2wdDjm4q>c?@|wJCkSa@{9Mw($y(XwPpj-X-zs- zH5R#4 z;u0CAmG~$KW+lGi4TP7U^*Y^SobGc?s&XiV0M?N(H#c7$^_g#DFUTcNPG`Nvy{u+2 zRY=jCd+|1Xl1|eGw^VsO6{Z+yAWAl35PATZIMq|=F2KwwG+3(VwjIXGu+wY|-&_Nba@_fL$7<4V6!YI~?Wj>tL};|P4$(0N787U^%7uTj1*3N*WdhIz_L%i% ztQ{M;osoz4+dY5wvmd9&+2G%?*b=tNndu)&p}Ur?j&^rsUeaiI()%m4*?T&?onc4y zrnOv$)wrlAAM?Auf*!CCkcO_m^Tck)P6Ff%Kf(pmi@pfwZ3e-=L{-`6%J;T27{ zoQF;Y&g+GWr3EEcme)>ZpqXI!;Xfti6+p$jr2{v0-r~+5y_rb+go@WT+kLdO7uejq zrG>24(a^~=kuCm-(D0}KPZj-j>Ew1;c&x@mPHL z0{?(0^LsBq*t_5B*YQ&c+|h zT#p#*4I@e_HI9))Deg0;V{wek8m=HiEi(PSgBvDHgKobnw^dYuj$W2;s}y>h>F3)D ztqe+(8tV0c5e6lo7{nL`Lp@+5Ad`?lU@kdC(Ku0lJPf|8&LYE*5R(HQr?!~- zi|T@eHp3RmI31jkgaS+ zF0oD|5d`FXX(Tx%tZSVDq&s92@Nvrk(PS2kt0)VZtonb9w=q%m$i=}Qd+4#)7ja9J z`XaC&UG1m@#*N|q;c3e9EKm&VlS%q}DZ^?ABaLX#A@eGG94OeQM5fZg0MtBy!SOn_ z2Z4txkFrZFvGL<&{kT?|mvBeLG7dtIkNQet;UTmmLt~#kN=9BE70j-k-jh?yd*LJo zR%TF4kE2W^CBAl@JP#jDNI(bcnNz4?T{s2oaUPstp^0@HJs51H6l7^Q+!7;F{6Uvs zrnmr*8l~EY^Wd~Gbhx>16;tly|{C|U`*x>92y^@0EqL3>%&Kfz|iUtC&5Md81UUfGlwSAPzKbtM(__{ z^eOs-i{U|_Pk+=5hCWYTWZnY$t4dFmlvzo)CbF^`w!zA|$O@?x;(d}nD$Y$7?-LkR Z<20I8@wN#5a~ORs;#7IIuU)(6{{na4c@_Ww diff --git a/src/jsrm/systems/tendon_actuated_planar_pcs.py b/src/jsrm/systems/tendon_actuated_planar_pcs.py index 5f6d6fd..b6acae8 100644 --- a/src/jsrm/systems/tendon_actuated_planar_pcs.py +++ b/src/jsrm/systems/tendon_actuated_planar_pcs.py @@ -5,7 +5,7 @@ import numpy as onp from typing import Callable, Dict, Optional, Tuple, Union -from .planar_pcs import factory as planar_pcs_factory, stifffness_fn +from .planar_pcs import factory as planar_pcs_factory def factory( num_segments: int, @@ -61,6 +61,8 @@ def actuation_mapping_fn( l = params["l"] # map the configuration to the strains xi = xi_eq + B_xi @ q + # segment indices + segment_indices = jnp.arange(num_segments) def compute_actuation_matrix_for_segment( segment_idx, d_sm: Array, @@ -70,7 +72,8 @@ def compute_actuation_matrix_for_segment( We assume that each segment is actuated by num_segment_tendons that are routed at a distance of d from the segment's backbone, respectively, and attached to the segment's distal end. We assume that the motor is located at the base of the robot and that the tendons are routed through all proximal segments. - The control inputs u1 and u2 are the tensions (i.e., forces) applied by the two tendons. + The positive control inputs u1 and u2 are the tensions (i.e., forces) applied by the two tendons. + At a straight configuration with a positive d1, a positive u1 and zero u2 should cause the bend negatively (to the right) and contract its length. Args: segment_idx: index of the segment @@ -78,53 +81,51 @@ def compute_actuation_matrix_for_segment( Returns: A_sm: actuation matrix of shape (n_xi, num_segment_tendons) """ - num_segment_tendons = d_sm.shape[0] - - A_sm = [] - for j in range(num_segment_tendons): - d = d_sm[j] - + def compute_A_d(d: Array) -> Array: """ - A_d = [] - for i in range(0, segment_idx + 1): - # length of the i-th segment - l_i = l[i] - # strain of the i-th segment - xi_i = xi[3 * i:3 * (i + 1)] - A_d.append(jnp.array([ - d * l_i * jnp.sqrt(xi_i[1]**2 + xi_i[2]**2), # the derivative of the tendon length with respect to the bending strain of the i-th segment - l_i * xi_i[1] * (1 + d * xi_i[0]) / jnp.sqrt(xi_i[1]**2 + xi_i[2]**2), # the derivative of the tendon length with respect to the shear strain of the i-th segment - l_i * xi_i[2] * (1 + d * xi_i[0]) / jnp.sqrt(xi_i[1]**2 + xi_i[2]**2), # the derivative of the tendon length with respect to the axial strain of the i-th segment - ])) + Compute the actuation matrix for a single actuator/tendon with respect to the soft robot's strains. + Args: + d: distance of the tendon from the centerline + Returns: + A_d: actuation matrix of shape (n_xi, ) where n_xi is the number of strains """ - def compute_A_d(l_i: Array, xi_i: Array) -> Array: - print("l_i.shape", l_i.shape, "xi_i.shape", xi_i.shape) + def compute_A_d_wrt_xi_i(i: Array, l_i: Array, xi_i: Array) -> Array: + """ + Compute the actuation matrix for a single actuator with respect to the strains of a single segment. + Args: + i: index of the segment + l_i: length of the segment + xi_i: strains for the segment + Returns: + A_d_segment: actuation matrix for the segment of shape (3, 3) + """ sigma_norm = jnp.sqrt(xi_i[1] ** 2 + xi_i[2] ** 2) - return jnp.array([ + A_d_wrt_xi_i = - jnp.array([ d * l_i * sigma_norm, l_i * xi_i[1] * (1 + d * xi_i[0]) / sigma_norm, l_i * xi_i[2] * (1 + d * xi_i[0]) / sigma_norm, ]) - A_d = vmap(compute_A_d)(l[:j+1], xi[: 3 * (j + 1)].reshape(-1, 3)) - + return jnp.where( + i <= segment_idx, + A_d_wrt_xi_i, + jnp.zeros_like(A_d_wrt_xi_i) + ) + + A_d = vmap(compute_A_d_wrt_xi_i)(segment_indices, l, xi.reshape(-1, 3)) + # concatenate the derivatives for all segments A_d = jnp.concatenate(A_d, axis=0) - A_sm.append(A_d) - - # stack the actuation matrices for all tendons - A_sm = jnp.stack(A_sm, axis=1) - print("A_sm.shape", A_sm.shape) + return A_d + + A_sm = vmap(compute_A_d, in_axes=0, out_axes=1)(d_sm) return A_sm - segment_indices = jnp.arange(num_segments) A_sms = vmap(compute_actuation_matrix_for_segment)( segment_indices, params["d"], ) - print("A_sms.shape", A_sms.shape) # concatenate the actuation matrices for all tendons A = jnp.concatenate(A_sms, axis=1) - print("A.shape", A.shape) # apply the actuation_basis A = A @ actuation_basis From a66de0be6eba51e593c987e18f53d25ca7be6501 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20St=C3=B6lzle?= Date: Mon, 26 May 2025 21:02:51 +0200 Subject: [PATCH 4/5] Fix bugs for multi-segment case --- .../simulate_tendon_actuated_planar_pcs.py | 3 +- pyproject.toml | 2 +- .../systems/tendon_actuated_planar_pcs.py | 31 ++++++------------- 3 files changed, 13 insertions(+), 23 deletions(-) diff --git a/examples/simulate_tendon_actuated_planar_pcs.py b/examples/simulate_tendon_actuated_planar_pcs.py index cb04fc0..08c5a57 100644 --- a/examples/simulate_tendon_actuated_planar_pcs.py +++ b/examples/simulate_tendon_actuated_planar_pcs.py @@ -123,7 +123,8 @@ def draw_robot( print("A =\n", A) x0 = jnp.concatenate([q0, jnp.zeros_like(q0)]) # initial condition - u = 1e0 * jnp.array([1.0, 1.0])[None].repeat(num_segments, axis=0).flatten() # tendon tensions + u = jnp.array([1.0, 1.0])[None].repeat(num_segments, axis=0).flatten() # tendon tensions + # u = 2e-1 * jnp.array([2.0, 0.0, 0.0, 1.0]) print("u =\n", u) ode_fn = ode_factory(dynamical_matrices_fn, params, u) diff --git a/pyproject.toml b/pyproject.toml index 6379774..59937fc 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,7 +17,7 @@ name = "jsrm" # Required # # For a discussion on single-sourcing the version, see # https://packaging.python.org/guides/single-sourcing-package-version/ -version = "0.0.14" # Required +version = "0.0.15" # Required # This is a one-line description or tagline of what your project does. This # corresponds to the "Summary" metadata field: diff --git a/src/jsrm/systems/tendon_actuated_planar_pcs.py b/src/jsrm/systems/tendon_actuated_planar_pcs.py index b6acae8..536824b 100644 --- a/src/jsrm/systems/tendon_actuated_planar_pcs.py +++ b/src/jsrm/systems/tendon_actuated_planar_pcs.py @@ -1,5 +1,5 @@ __all__ = ["factory", "stiffness_fn"] -from jax import Array, lax, vmap +from jax import Array, debug, lax, vmap import jax.numpy as jnp from jsrm.math_utils import blk_diag import numpy as onp @@ -24,17 +24,10 @@ def factory( if segment_actuation_selector is None: segment_actuation_selector = jnp.ones(num_segments, dtype=bool) - # number of input pressures - actuation_dim = segment_actuation_selector.sum() * 2 - - # matrix that maps the (possibly) underactuated actuation space to a full actuation space - actuation_basis = jnp.zeros((2 * num_segments, actuation_dim)) - actuation_basis_cumsum = jnp.cumsum(segment_actuation_selector) - for i in range(num_segments): - j = int(actuation_basis_cumsum[i].item()) - 1 - if segment_actuation_selector[i].item() is True: - actuation_basis = actuation_basis.at[2 * i, j].set(1.0) - actuation_basis = actuation_basis.at[2 * i + 1, j + 1].set(1.0) + # number of system inputs + actuation_dim = segment_actuation_selector.sum() + actuation_dim = 2 * num_segments + actuation_basis = jnp.eye(actuation_dim) def actuation_mapping_fn( forward_kinematics_fn: Callable, @@ -106,26 +99,22 @@ def compute_A_d_wrt_xi_i(i: Array, l_i: Array, xi_i: Array) -> Array: l_i * xi_i[2] * (1 + d * xi_i[0]) / sigma_norm, ]) return jnp.where( - i <= segment_idx, + i * jnp.ones((3, )) <= segment_idx * jnp.ones((3, )), A_d_wrt_xi_i, jnp.zeros_like(A_d_wrt_xi_i) ) - A_d = vmap(compute_A_d_wrt_xi_i)(segment_indices, l, xi.reshape(-1, 3)) + A_d = vmap(compute_A_d_wrt_xi_i)(segment_indices, l, xi.reshape(-1, 3)).reshape(-1) - # concatenate the derivatives for all segments - A_d = jnp.concatenate(A_d, axis=0) return A_d - A_sm = vmap(compute_A_d, in_axes=0, out_axes=1)(d_sm) + A_sm = vmap(compute_A_d, in_axes=0, out_axes=-1)(d_sm) return A_sm - A_sms = vmap(compute_actuation_matrix_for_segment)( + A = vmap(compute_actuation_matrix_for_segment, in_axes=(0, 0), out_axes=0)( segment_indices, params["d"], - ) - # concatenate the actuation matrices for all tendons - A = jnp.concatenate(A_sms, axis=1) + ).transpose((1, 0, 2)).reshape(xi.shape[0], -1) # apply the actuation_basis A = A @ actuation_basis From a2823da75d986d47b3514298b3d40cd1bf926e2b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20St=C3=B6lzle?= Date: Mon, 26 May 2025 21:18:48 +0200 Subject: [PATCH 5/5] Add support for `segment_actuation_selector` --- examples/simulate_tendon_actuated_planar_pcs.py | 5 ++++- src/jsrm/systems/tendon_actuated_planar_pcs.py | 16 ++++++++-------- 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/examples/simulate_tendon_actuated_planar_pcs.py b/examples/simulate_tendon_actuated_planar_pcs.py index 08c5a57..0bf3555 100644 --- a/examples/simulate_tendon_actuated_planar_pcs.py +++ b/examples/simulate_tendon_actuated_planar_pcs.py @@ -44,6 +44,9 @@ # activate all strains (i.e. bending, shear, and axial) strain_selector = jnp.ones((3 * num_segments,), dtype=bool) +# actuation selector for the segments +segment_actuation_selector = jnp.ones((num_segments,), dtype=bool) +# segment_actuation_selector = jnp.array([False, True]) # only the last segment is actuated # define initial configuration q0 = jnp.repeat(jnp.array([5.0 * jnp.pi, 0.1, 0.2])[None, :], num_segments, axis=0).flatten() @@ -99,7 +102,7 @@ def draw_robot( if __name__ == "__main__": strain_basis, forward_kinematics_fn, dynamical_matrices_fn, auxiliary_fns = ( - planar_pcs.factory(num_segments, sym_exp_filepath, strain_selector) + planar_pcs.factory(num_segments, sym_exp_filepath, strain_selector, segment_actuation_selector=segment_actuation_selector) ) actuation_mapping_fn = auxiliary_fns["actuation_mapping_fn"] # jit the functions diff --git a/src/jsrm/systems/tendon_actuated_planar_pcs.py b/src/jsrm/systems/tendon_actuated_planar_pcs.py index 536824b..194a4f5 100644 --- a/src/jsrm/systems/tendon_actuated_planar_pcs.py +++ b/src/jsrm/systems/tendon_actuated_planar_pcs.py @@ -24,11 +24,6 @@ def factory( if segment_actuation_selector is None: segment_actuation_selector = jnp.ones(num_segments, dtype=bool) - # number of system inputs - actuation_dim = segment_actuation_selector.sum() - actuation_dim = 2 * num_segments - actuation_basis = jnp.eye(actuation_dim) - def actuation_mapping_fn( forward_kinematics_fn: Callable, jacobian_fn: Callable, @@ -112,12 +107,17 @@ def compute_A_d_wrt_xi_i(i: Array, l_i: Array, xi_i: Array) -> Array: return A_sm + # compute the actuation matrix for all segments + # will have shape (num_segments, n_xi, num_segment_tendons) A = vmap(compute_actuation_matrix_for_segment, in_axes=(0, 0), out_axes=0)( segment_indices, params["d"], - ).transpose((1, 0, 2)).reshape(xi.shape[0], -1) + ) + + # deactivate the actuation for some segments + A = A[segment_actuation_selector] - # apply the actuation_basis - A = A @ actuation_basis + # reshape the actuation matrix to have shape (n_xi, n_act) + A = A.transpose((1, 0, 2)).reshape(xi.shape[0], -1) return A