From 0f74d447c497028c07a694441185edca01078eaf Mon Sep 17 00:00:00 2001 From: javoha Date: Sat, 31 May 2025 16:19:18 +0200 Subject: [PATCH 01/12] Test for computation with topography. --- .../test_compute_at_computation_time.py | 60 +++++++++++++++++++ 1 file changed, 60 insertions(+) create mode 100644 test/test_modules/test_compute_at_computation_time.py diff --git a/test/test_modules/test_compute_at_computation_time.py b/test/test_modules/test_compute_at_computation_time.py new file mode 100644 index 00000000..10909939 --- /dev/null +++ b/test/test_modules/test_compute_at_computation_time.py @@ -0,0 +1,60 @@ +import numpy as np + +import gempy as gp +import time + +PLOT = True + + +def test_compute_at_computation_time(): + + # Define the path to data + data_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/' + path_to_data = data_path + "/data/input_data/jan_models/" + + # Create a GeoModel instance + geo_model = gp.create_geomodel( + project_name='EGU_example', + extent=[0, 2500, 0, 1000, 0, 1000], + resolution=[125, 50, 50], + importer_helper=gp.data.ImporterHelper( + path_to_orientations=path_to_data + "model7_orientations.csv", + path_to_surface_points=path_to_data + "model7_surface_points.csv" + ) + ) + + # Map geological series to surfaces + gp.map_stack_to_surfaces( + gempy_model=geo_model, + mapping_object={ + "Fault_Series": ('fault'), + "Strat_Series1": ('rock3'), + "Strat_Series2": ('rock2', 'rock1'), + } + ) + + # Define youngest structural group as fault + gp.set_is_fault(geo_model, ["Fault_Series"]) + + # Compute a solution for the model + start_time = time.perf_counter() + gp.compute_model(geo_model) + end_time = time.perf_counter() + computation_time_model = end_time - start_time + + # Setting a randomly generated topography + gp.set_topography_from_random( + grid=geo_model.grid, + fractal_dimension=2, + d_z=np.array([700, 950]), + topography_resolution=np.array([125, 50]) + ) + + # Recompute model as a new grid was added + start_time = time.perf_counter() + gp.compute_model(geo_model) + end_time = time.perf_counter() + computation_time_topo = end_time - start_time + + print(f"Computation time without topography: {computation_time_model:.2f} seconds") + print(f"Computation time with topography: {computation_time_topo:.2f} seconds") From 248314170f6efc6af6f07c7bb24386e03567a25c Mon Sep 17 00:00:00 2001 From: javoha Date: Sat, 31 May 2025 16:19:18 +0200 Subject: [PATCH 02/12] Test for computation with topography. --- .../test_compute_at_computation_time.py | 60 +++++++++++++++++++ 1 file changed, 60 insertions(+) create mode 100644 test/test_modules/test_compute_at_computation_time.py diff --git a/test/test_modules/test_compute_at_computation_time.py b/test/test_modules/test_compute_at_computation_time.py new file mode 100644 index 00000000..10909939 --- /dev/null +++ b/test/test_modules/test_compute_at_computation_time.py @@ -0,0 +1,60 @@ +import numpy as np + +import gempy as gp +import time + +PLOT = True + + +def test_compute_at_computation_time(): + + # Define the path to data + data_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/' + path_to_data = data_path + "/data/input_data/jan_models/" + + # Create a GeoModel instance + geo_model = gp.create_geomodel( + project_name='EGU_example', + extent=[0, 2500, 0, 1000, 0, 1000], + resolution=[125, 50, 50], + importer_helper=gp.data.ImporterHelper( + path_to_orientations=path_to_data + "model7_orientations.csv", + path_to_surface_points=path_to_data + "model7_surface_points.csv" + ) + ) + + # Map geological series to surfaces + gp.map_stack_to_surfaces( + gempy_model=geo_model, + mapping_object={ + "Fault_Series": ('fault'), + "Strat_Series1": ('rock3'), + "Strat_Series2": ('rock2', 'rock1'), + } + ) + + # Define youngest structural group as fault + gp.set_is_fault(geo_model, ["Fault_Series"]) + + # Compute a solution for the model + start_time = time.perf_counter() + gp.compute_model(geo_model) + end_time = time.perf_counter() + computation_time_model = end_time - start_time + + # Setting a randomly generated topography + gp.set_topography_from_random( + grid=geo_model.grid, + fractal_dimension=2, + d_z=np.array([700, 950]), + topography_resolution=np.array([125, 50]) + ) + + # Recompute model as a new grid was added + start_time = time.perf_counter() + gp.compute_model(geo_model) + end_time = time.perf_counter() + computation_time_topo = end_time - start_time + + print(f"Computation time without topography: {computation_time_model:.2f} seconds") + print(f"Computation time with topography: {computation_time_topo:.2f} seconds") From 993927c489ee7f9400b1348eac4d775442674cc4 Mon Sep 17 00:00:00 2001 From: Deepprakash222 Date: Mon, 2 Jun 2025 19:51:10 +0200 Subject: [PATCH 03/12] added a test_gradient folder to generate sample for c, m=Gempy(c) and dm/dc --- test/test_gradient/generate_samples.py | 211 +++++++++++++++++++++++++ test/test_gradient/helpers.py | 81 ++++++++++ test/test_gradient/prior_model.png | Bin 0 -> 29187 bytes test/test_gradient/save_gempy_data.py | 70 ++++++++ 4 files changed, 362 insertions(+) create mode 100644 test/test_gradient/generate_samples.py create mode 100644 test/test_gradient/helpers.py create mode 100644 test/test_gradient/prior_model.png create mode 100644 test/test_gradient/save_gempy_data.py diff --git a/test/test_gradient/generate_samples.py b/test/test_gradient/generate_samples.py new file mode 100644 index 00000000..9b3d6372 --- /dev/null +++ b/test/test_gradient/generate_samples.py @@ -0,0 +1,211 @@ + +import torch +import numpy as np +from datetime import datetime +from torch.autograd.functional import jacobian +import pyro +import pyro.distributions as dist +from pyro.infer import Predictive + +from pyro.nn import PyroModule + +import gempy as gp +import gempy_engine +from gempy_engine.core.backend_tensor import BackendTensor + + +from helpers import * + + +class GempyModel(PyroModule): + def __init__(self, interpolation_input_, geo_model_test, num_layers, slope, dtype): + super(GempyModel, self).__init__() + + BackendTensor.change_backend_gempy(engine_backend=gp.data.AvailableBackends.PYTORCH) + self.interpolation_input_ = interpolation_input_ + self.geo_model_test = geo_model_test + self.num_layers = num_layers + self.dtype = dtype + self.geo_model_test.interpolation_options.sigmoid_slope = slope + ############################################################################### + # Seed the randomness + ############################################################################### + seed = 42 + np.random.seed(seed) + torch.manual_seed(seed) + # Ensure deterministic behavior + torch.backends.cudnn.deterministic = True + torch.backends.cudnn.benchmark = False + # Setting the seed for Pyro sampling + + pyro.set_rng_seed(42) + + + def create_sample(self): + + """ + This Pyro model represents the probabilistic aspects of the geological model. + It defines a prior distribution for the top layer's location and + computes the thickness of the geological layer as an observed variable. + + + interpolation_input_: represents the dictionary of random variables for surface parameters + geo_model_test : gempy model + + num_layers: represents the number of layers we want to include in the model + + """ + + Random_variable ={} + + # Create a random variable based on the provided dictionary used to modify input data of gempy + counter=1 + for interpolation_input_data in self.interpolation_input_[:self.num_layers]: + + # Check if user wants to create random variable based on modifying the surface points of gempy + if interpolation_input_data["update"]=="interface_data": + # Check what kind of distribution is needed + if interpolation_input_data["prior_distribution"]=="normal": + mean = interpolation_input_data["normal"]["mean"] + std = interpolation_input_data["normal"]["std"] + Random_variable["mu_"+ str(counter)] = pyro.sample("mu_"+ str(counter), dist.Normal(mean, std)) + #print(Random_variable["mu_"+ str(counter)]) + + elif interpolation_input_data["prior_distribution"]=="uniform": + min = interpolation_input_data["uniform"]["min"] + max = interpolation_input_data["uniform"]["min"] + Random_variable["mu_"+ str(interpolation_input_data['id'])] = pyro.sample("mu_"+ str(interpolation_input_data['id']), dist.Uniform(min, max)) + + + else: + print("We have to include the distribution") + + # # Check which co-ordinates direction we wants to allow and modify the surface point data + counter=counter+1 + + + + def GenerateInputSamples(self, number_samples): + + pyro.clear_param_store() + # We can build a probabilistic model using pyro by calling it + + dot = pyro.render_model(self.create_sample, model_args=()) + # Generate 50 samples + num_samples = number_samples # N + predictive = Predictive(self.create_sample, num_samples=num_samples) + samples = predictive() + + samples_list=[] + for i in range(len(self.interpolation_input_)): + samples_list.append(samples["mu_"+str(i+1)].reshape(-1,1)) + ######store the samples ###### + parameters=torch.hstack(samples_list) # (N, p) = number of sample X number of paramter + + return parameters.detach().numpy() + + def GempyForward(self, *params): + index=0 + interpolation_input = self.geo_model_test.interpolation_input + + for interpolation_input_data in self.interpolation_input_[:self.num_layers]: + # Check which co-ordinates direction we wants to allow and modify the surface point data + if interpolation_input_data["direction"]=="X": + interpolation_input.surface_points.sp_coords = torch.index_put( + interpolation_input.surface_points.sp_coords, + (torch.tensor([interpolation_input_data["id"]]), torch.tensor([0])), + params[index]) + elif interpolation_input_data["direction"]=="Y": + interpolation_input.surface_points.sp_coords = torch.index_put( + interpolation_input.surface_points.sp_coords, + (torch.tensor([interpolation_input_data["id"]]), torch.tensor([1])), + params[index]) + elif interpolation_input_data["direction"]=="Z": + interpolation_input.surface_points.sp_coords = torch.index_put( + interpolation_input.surface_points.sp_coords, + (interpolation_input_data["id"], torch.tensor([2])), + params[index]) + else: + print("Wrong direction") + + index=index+1 + + self.geo_model_test.solutions = gempy_engine.compute_model( + interpolation_input=interpolation_input, + options=self.geo_model_test.interpolation_options, + data_descriptor=self.geo_model_test.input_data_descriptor, + geophysics_input=self.geo_model_test.geophysics_input, + ) + + # Compute and observe the thickness of the geological layer + + m_samples = self.geo_model_test.solutions.octrees_output[0].last_output_center.custom_grid_values + return m_samples + + def GenerateOutputSamples(self, Inputs_samples): + + Inputs_samples = torch.tensor(Inputs_samples, dtype=self.dtype) + m_data =[] + dmdc_data =[] + for i in range(Inputs_samples.shape[0]): + params_tuple = tuple([Inputs_samples[i,j].clone().requires_grad_(True) for j in range(Inputs_samples.shape[1])]) + + m_samples = self.GempyForward(*params_tuple) + m_data.append(m_samples) + J = jacobian(self.GempyForward, params_tuple) + J_matrix = torch.tensor([[J[j][i] for j in range(len(J))] for i in range(J[0].shape[0])]) + dmdc_data.append(J_matrix) + + return torch.stack(m_data).detach().numpy() , torch.stack(dmdc_data).detach().numpy() + + + + +def generate_input_output_gempy_data(mesh, number_samples, slope=200, filename=None): + # ---------------- 1️⃣ Check and create a 3D custom grid ---------------- + mesh_coordinates = mesh + data ={} + geo_model_test = create_initial_gempy_model(refinement=3, save=False) + if mesh_coordinates.shape[1]==2: + xyz_coord = np.insert(mesh_coordinates, 1, 0, axis=1) + elif mesh_coordinates.shape[1]==3: + xyz_coord = mesh_coordinates + + gp.set_custom_grid(geo_model_test.grid, xyz_coord=xyz_coord) + geo_model_test.interpolation_options.mesh_extraction = False + + sp_coords_copy_test = geo_model_test.interpolation_input.surface_points.sp_coords.copy() + + ############################################################################### + # 2️⃣ Make a list of gempy parameter which would be treated as a random variable + ############################################################################### + dtype =torch.float64 + test_list=[] + std = 0.06 + test_list.append({"update":"interface_data","id":torch.tensor([1]), "direction":"Z", "prior_distribution":"normal","normal":{"mean":torch.tensor(sp_coords_copy_test[1,2],dtype=dtype), "std":torch.tensor(std,dtype=dtype)}}) + test_list.append({"update":"interface_data","id":torch.tensor([4]), "direction":"Z", "prior_distribution":"normal","normal":{"mean":torch.tensor(sp_coords_copy_test[4,2],dtype=dtype), "std":torch.tensor(std,dtype=dtype)}}) + num_layers = len(test_list) # length of the list + + Gempy = GempyModel(test_list, geo_model_test, num_layers, slope=slope, dtype=torch.float64) + + # ---------------- 3️⃣ Generate the samples for input parameters. c ∼ N(µ, Σ) ---------------- + + start_sample_input_time = datetime.now() + c = Gempy.GenerateInputSamples(number_samples=number_samples) + end_sample_input_time = datetime.now() + total_time_input = end_sample_input_time - start_sample_input_time + + data["input"] = c.tolist() + # ---------------- 3️⃣ Generate the samples for output of gempy and the Jacobian matrix. m= Gempy(c) and dm/dc ---------------- + start_sample_output_time = datetime.now() + m_data, dmdc_data = Gempy.GenerateOutputSamples(Inputs_samples=c) + end_sample_output_time = datetime.now() + total_time_output = end_sample_output_time - start_sample_output_time + + data["Gempy_output"] = m_data.tolist() + data["Jacobian_Gempy"] = dmdc_data.tolist() + + return data, total_time_input, total_time_output + + + diff --git a/test/test_gradient/helpers.py b/test/test_gradient/helpers.py new file mode 100644 index 00000000..a6cb7f04 --- /dev/null +++ b/test/test_gradient/helpers.py @@ -0,0 +1,81 @@ + +import numpy as np + +import matplotlib.pyplot as plt + +import gempy as gp +import gempy_viewer as gpv + + + +def create_initial_gempy_model(refinement,filename='prior_model.png', save=True): + """ Create an initial gempy model objet + + Args: + refinement (int): Refinement of grid + save (bool, optional): Whether you want to save the image + + """ + geo_model_test = gp.create_geomodel( + project_name='Gempy_abc_Test', + extent=[0, 1, -0.1, 0.1, 0, 1], + resolution=[100,10,100], + refinement=refinement, + structural_frame= gp.data.StructuralFrame.initialize_default_structure() + ) + + brk1 = 0.3 + brk2 = 0.5 + + + gp.add_surface_points( + geo_model=geo_model_test, + x=[0.1,0.5, 0.9], + y=[0.0, 0.0, 0.0], + z=[brk1, brk1 , brk1], + elements_names=['surface1','surface1', 'surface1'] + ) + + gp.add_orientations( + geo_model=geo_model_test, + x=[0.5], + y=[0.0], + z=[0.0], + elements_names=['surface1'], + pole_vector=[[0, 0, 0.5]] + ) + geo_model_test.update_transform(gp.data.GlobalAnisotropy.NONE) + + element2 = gp.data.StructuralElement( + name='surface2', + color=next(geo_model_test.structural_frame.color_generator), + surface_points=gp.data.SurfacePointsTable.from_arrays( + x=np.array([0.1,0.5, 0.9]), + y=np.array([0.0, 0.0, 0.0]), + z=np.array([brk2, brk2, brk2]), + names='surface2' + ), + orientations=gp.data.OrientationsTable.initialize_empty() + ) + + + geo_model_test.update_transform(gp.data.GlobalAnisotropy.NONE) + + geo_model_test.structural_frame.structural_groups[0].append_element(element2) + gp.add_orientations( + geo_model=geo_model_test, + x=[0.5], + y=[0.0], + z=[1.0], + elements_names=['surface2'], + pole_vector=[[0, 0, 0.5]] + ) + geo_model_test.structural_frame.structural_groups[0].elements[0], geo_model_test.structural_frame.structural_groups[0].elements[1] = \ + geo_model_test.structural_frame.structural_groups[0].elements[1], geo_model_test.structural_frame.structural_groups[0].elements[0] + + gp.compute_model(geo_model_test) + if save: + picture_test = gpv.plot_2d(geo_model_test, cell_number=5, legend='force') + plt.savefig(filename) + + return geo_model_test \ No newline at end of file diff --git a/test/test_gradient/prior_model.png b/test/test_gradient/prior_model.png new file mode 100644 index 0000000000000000000000000000000000000000..54c15ee67c182ef3ec387893f330ee5a402f30e4 GIT binary patch literal 29187 zcmdqJ2UL@7w=EhgHb4|XK|nx7q>1zz6lqeWBTe9=H|f3jDTq`R6{$h0^xjKUN))7p zUL!cntTor1b3L!LG?b{%Fr9(HVALwg z542&h7?`?S4akAwdxA3p{Dm19406nxQA4h8?h=r3G@zJ&d+|A<)& z+@~%X2J^Zf3j5-3Txs`6diBdOLEqJ1;^1qhZBJTXjVoq**%rU|d{q*1u;kcCrB5+B z8LF;V;VuRXRqt$lbrGv-zN@`I52w5mbsqa2l_|41!b#r4P}&ic{rP?I#m|{DNsh&@ zZmM@6bV3LLn5a1L%zb$}x$a`wU;T(}jO*n}ajbns#x)|Yp5OH!gj0&jr;bhB3kC2#>&~ z3oE)K| zCX-4+2DeEDg@$jC8jHnUO>e}e5`K?zGlv6*MHK3EJlXSVmfSHF&_^O*naXg{LX`GL>l@&*{^;;NJP z^K$N=X;6&4ni$f)36}Q0zXAUJTYurJ{tKUW0#>JKL&RTN#@|%WCpzZ#E_16jM?>UzoQ|GB;hM?v^6<}EslP|x3@rg|vPir!nHr4Tfm+uV>BI z;_H{z1aSv7J+gbh=cQLF`}2_+>hu1yi2fusgscpe9B#fm2T`)YszJIHE4#bkB})s7 zc@8@g`gjMd%;W@IMs<0)d+J^=mKpP%(FDhld)?H8EnvCKPdZ@n)6?F2R0=O?Sk~Gc z9-G%EH+{;EYSukAo$537s#5?dy@(b=?sQ4h$cFLm;Tv3=KQ_%o@~mxH_NLz>4kO7T zn^OUV?Jrca=R@|!F*EOUH1f)}F%fi}Gh>b=F-6s5x2Xhi90G&ylo^Y7eBw1J&8y$e zdMK}T@Ra;sy(W_ufUgtZTdu%{GrN2)OXaRqVWi;d`sDCx&mMjgzLFKS<>@pQGVWXx18p%s%0?6aL5dd+ zeS;!(Fp*2qoOd-PvA$cs`oU%LYcln1P>gsx08e3?%0`+m2L$>F=c>G8+ zBQEx7y+=y8Y}P0fx=s7-(HEZgl85512<2$hXXMIplG2qK9dX@(rk7^+KmcO}d(VgF z^jqO$56DGKyV3ytQ8wo+|J_QxcaY(mag&`0Nv!`H^O-1cW^|*UUcmSbEjFx>=OZ!uMnH6ft*R82u zIB+h-qx*-5#M=$Yc-NpMNCRdt7NaJHJ`p$=^ADeb80e=b(Q}mh@Y>No(lliTjMEo{z{ryNSxv_uldkfE_V) za~=^AOlHd7QhCe2YT&w+onV#}@Hp$&C;HOsxa|R*PofwK{jcnKjVe`bSY`Rf7&_qY z2@Y5~)lpcLHdCxty|O7A4!=;*Bb!X<$eoGeF{jUtN*M5xGNsRns;z(Re?04Ki?+Zs zzxBDnAHyn)M%k;^c3Qb)A02?)m98tqF-a2W64%m6ji@Z5XT6To@UGWVN4Ht3)N5&C zIz25K*8Eg_D*M}<7eBq`=QT(e(;Ew!oe8$ji4;Db?MZCp=N4Ah-rDl~ni_7BrL^nA zB>LL_LDm!_vo>xo0smoGadG<5TeJB}`~J05%77yWY0R*WOMTJkH6EX9b`p1A(x}nb ztvOUf=z)pXpx)s0J>0zLWS!=?!+6?T`u(hV4SDs4hlBN!mMsmL%U!AbY@kDQ+^Hha z%)rG5Cz5kCmq(RdOEO_)yhgxD^UT0lm1DxO{GIiAgA&bqs**jpz~Tb26R?<WqMGx9^$zviduH(jHDk;T91WkH&{lHGvBafA$;`j}c9<2YnY zm5ccYr43DYr>YzVZ~Js8mk-?=ua7rat6yr$rX z3wVDrJGl*A1cQ2GGo7bzvb=MyNT5#y)IV?f)ZCMYzptIl&FttO zlZ|N6T`tvKExnmdSoxVKw0?I`dgS%P9^-yLGru5vw7nGRuM==lH?D6Bc?nV_}oO1byjn-OcmSHpNgq66~2*o zA;|P$+ZwudTew>~hDDQ6d5I@R$2^%@4{1Sa_$fj<$z8>HLdu@uS%%aW57j9iTI0_= z+3Xi^G~{sYq#G8wsou9IYI24RcO}%C?hTQyW+VF0Wah_xZNa#uV$Fq(;Ox;G+A*&- zn$<$_>$ZCFrrP(gGmhHzS_QbH<;9C!g}l@i%3nJd+2!hm()@+8q*k*D z_KUKJqi=RUM2hB2otRGxaA9>s)}sef@=9mIU*p?z>tWoNm8*oOqR?Bspga%~9uPfl zWZd&nc!1qmRfcycS5a>9)0-(~`wy{6&2amJ4s&F4RvJ~F-Cy}(T1gq`jTl&=8Wg^7 z=}??_sREn*3Z(jQgTgp^v+VP;UyjiR_I!GPn_vCGh4|Sbsh?oe zs~mPm6^_YMQrnf-K~naROpfQWjQF+Py0z)r#x~F2s0t|%)i({otVR4ZF*i&HT=;>1~oLc#HGqqd0Dv&Yh& zemh+TQp#@dWMZW{jty3&gWGpfnhUV{k^Xih#b%Z)HG&pR8ic>2xQn|jBRqVG?V zh)-tux9uM{)$?u`zRiKFVG7<8<`Qg#mHQ<=WZfp??u&;V&o;XIR+p2Vx=)#OLAzJW z*tLO1#!JRp!3%+QvZ{ewO>%&MzwX z$Km(>v-Cj}bZO`ZKI6GB9Fr>oD)6{GD!C*MAsIqP|J~w}jrG^vv%3SlDqq!-j5tmN zzy)rd@zNzZ$WgiXW@uAb4mao6VB<=fj zAKdIa<6K@;aD58}6ds&}^%&Pc&f#011#lu+QcVS#AqiRB^!-b!KBKfU8@&_j9H!B^ z!wLMimg;jKd8jNQ0{ELrO~3H{iw-Hg~X@*@lDPvH6fF%hQ+*JW|FcXuM?&_ER^2ROTWt* zQ`3gMiM*L}cF*%AJctc_B3QlU#*oq(OYPOSHk!&lWBom@XTE6*Jl9?*bMHfK$fz7q zVS2j`%DGvB+hSP!P_`WBfl0!<;NbpxPdeGZzQ3OEu=UEjWwm%BKF3LVsjwo$pv272 zbg}A94oOW$Ra<#0*g-E4pt7DY4xO$zfxzjrRkHl4TLYULX(E9~LI+QuhZzdsW^SwZ zB;Zfw7H{J2-zJGG>oDX%86D1DFh8TLPj1e!(y@DO)(MITL`wwqBTaT0rQPN!_f9^)Q;NG+H0Fg*T*EnD zb^p#@u#EUtwsD~T)rcu6-v3;dclMXLDD(Z9+s_owXrs5-y(pS+u{{DZB<9F{`}}Tz zV|t^7jTIKXUg!K)BZLh|*^szzEve=;;r7a^jYMH~jkois2#s(UO!$<&_@l*VZ~c|& zyQ+`MVKYx2eb5T3a)#VAe}|$?^|lV}gx&!WaT#sK0n!z#Gce}(f#`Ft(ijJ^SG>a-Rp+NWQ5Li(Lk4QOoDn_6UZk`~7#oX0AqHEM+MBA1X4!^+Vkc;Pul znCK1rp*OlboP*xtJHO*Sv}j-=kB&W4%oA=iS(JF4llQPcI~&)e{ixOphe3vsXew+DP7I(LH7|F4?E=Pa;T`sOw`{227%;X}vUU={_phV=Y zPaTWVPAcb7$LJ;hO8;wj1LtSdFKo#o>G^#M3H322sr;l-e@e)zDTy@2q6`A#39dXX{g|5 zLQ^?-H(|>lrzX;ngtYk5psA=9&ylGD+9IoQ7m9I~5!uWbIOb)e#9yvAxQcd)acrJqyW+`tcPzT*ONpsb3sQL&DHgH4a{1ZOF)i7aEadA^k z0ZIw{rfc;UE#YkoeVJCPzrG53F7%r0{Q{K6D^?l*);=XRQ~XVE=YcCS;FIk!%#1BB zX?r@~Zmx`@CVen04fuM;=*xl}co0(;z~s}=p)2o!r=iwUcweWH?zBr|CIqnN<;%b z*DAiG^W2K~woP!q0TYcv#qg;nhmtZuxde;(hTb8xFIkZVat{3N^QP9^N~&?bG=6n} zQr#w9rLcv{#aIesttaOZiLvBkFl^b`npT={vPoe-`KpD5?i2Qia&9b zh(%x+sPN-8PR_?1l)zBPJr^SxU}wSzo z97Q)eO_uA8jsYU_62+IeGvXRI1&SVYM!CW=(g0CQdDABtfKM4IpPaLHk$ z_sPS=HOG?HnH#DZ@_V$csJ*M_TQtN);F>=iUDdk0&XQffd6YTrl#lqK@NRQNM68U@ zicwR~yJ5;jIE$ZcFtus|UsS5k0)u?44EN@^6=CrTeXw^MNu)dLj_uS$qtn>Hg2TuN zw5aTQTNuE)8}pqGKt?*#khPjJGZheX7 zCMEO!F6h>6!iK%u<_~#n%$0sCL|T)3(6hccnYStb@sXe09%HS*(*B5sn+P$}#O~iH<$4u4} zQXF|xymrF+f}{Abm6<0M`fb#)aaNQIj3#?A+q9 znp8=phORM}P{K{E1ZXD+8|KY@U2I_(vf24|gShM3t~7j)B*qZ3nc_2($e7l}xXi^C zY+~j!fTE4bkMOl>uOtOGD@0*pN`lCpec0y>`XvF5n}t`eEsuMYEHv%0`b`EQ{qFeV zD;)#(_B-fXg)-xZYqm$hh51vT7jPdaVp5AqHNdhrD;EnApW*rD{VzT~3ps0C&X%Yv0(?lJ;9;HYh3{{N+uVFK%Z4*3L*X+m8x(qxjD3h*iR*!Om`{0dE|S%f6ZP!Tz@R zZ%(RqF4-}HOp{r9@ZoYQMLV5cvKrm@$&}50e5~bV+Qs1sNe=!da{fLg0KfZ5{Ie!T zd^t3@?&P=Kkw*1#J9Jj<2RVo~2AoSQJT^$$EC5UL4cbkXkW;Q6by#Xl5ehM*xu{P= z`>@_?R8W9RBkf%Ra{r)Z=TSV)wI?Ei!g!6>s9w0)t)azzT}AZ3II6_#yljfHuu7AlL`udfDtrz4gH z$(mqgO>4$K?P8*)g2ecn@y0tMIirll+h_xat5Yzils`)migKERhD$)qq1Z5drag-> zeP?@bVa#t&d?x+K&YK+5#_9Q*w(m-c^VitJj%MdwHV|D(9VOr#HIgMbUsAnfdHU4E zdnmyms9Sv$>djsXz2sP*StwqkRN9K_D0+b?x2NSFEyD z7qsZX5AFq$b_)s+gzA_0nA|jGHW1dun%%meP9JglyV3Co%&rvuE?P!F{5BMZ8Svxl zol1jy^szFswW=Fq*-HBxWGSm%W1-BwO;fA(aUlnmJmX@S4qFCG?RJsgGvXc7=mVGk{g3B*-)@A2YH%GU>y*Ga)#DAb2SWcO) zqiJDF+gB})dmXFU#xUQFV1%=+Yo*b+BUqcL^gZ(sQv9Wgp9B2zVp5gu>daf*7m-*9p=Dl=MDn*v2Bc&%kVTaLP_Eokzk)~vGsGO?LR(9y6fIhRP)M1b5D)e zxLP{+>9Bg;%a&qOpAL%_b)@+0nyxxm!?)*GWYJ~y5 zv+FZ6M1_B~7rd^;YuwCRjj8qsBvs6}eTgIW>}B1tAikw^^{zIJyB3z2BKPmn9y#jt z`*wbQKALbY{T#qnjO*r!thP<{s&&anVazd;#qb;TTa#AL;e6W3ky}Wn3~BZMlt5Jx{^gi=6KDRl56Z(kNMy)aBVI#IU5Vdg*&Vb7_6^V zW>U2sY1-oP@vz*mLdGoNw^hKv7eW_`hc0xIFyGU3Az`0Lv7AQh2|T9kBCfc#VJXoC zdjuA)l}BW$buSBCdA8qZy>)@R*p2R{UG(LxZ*WVbj6X`sW4<+IXT4+c`V*N`5D-|M z{+bZHH}JS7SL;({(^d9w6DQU!TTR7K+pl`s3MDEsp`G#TK$b~9PDbLfu710Ofyq@6 zURY^OkPa(<4-V&2N}Q&oy;tgWCMe$*?D{ek3JuFR=ExLAoO$iq+s6HI8LPFmvfQ)? zs$j~J?@pD0CGqWzzE1~KGuR-S2n(*1Bt-7-DzTh)K4ha5gmLPIzjRA^#|gzzJeTVn zdlrLo8p()l+;^#UN3eN~mK)`3F1s(Qlwnu{O3-umuIh1T!vkewabkCB02oYWINV(R zBlj2mw3yBPHYhU@ru74RPs2ME1EtJ@N$J1flJ?=>qam$rtY+n z{EO#rsz*~a5Dc~mY2m4lf-;-9Q65${F;PYMrZVfdQmb5kKcx8q&d05bk8|yHIKa`@ z0oTJ22Nx2>v#iTarNwIIXzaP3+Z)^IE3~f=clknBrZC>Me*VRJbCCW1cDbWf4PI)V z>8!13utg=tvQvEJ^S!=GqyBB|2(tCWa4=_c@Mc&ag_Eto1!FqAFP2#)io!R6n#G2&lah+MJSD>|B;0vd%I?=T;y=f+G3n6}+oSE5 z4g#YXkj&NQ#2eN&#| zdrf|o(q99ywZU?9ZI)QUQykYl?%A=VeeKS|zL3YrYLS51h2hDa_>qjsIk)jkoEdtZ zAN&wls;igG8`s{}l8cjPz68OwS(r}IKIy+HcbJh+TpV_!B$$^B!z?O>k`)R*6wQ2s zN568Ga=OSiWkm?2BfsxS*-J|+QLv60DKd83BQ)-YUlo1#?fGAu6FwO$HqMrP4N+^9 z+-X`(@_OR7a!F8r%k#h1t5;^)OH7j!Jm$J$LK?RNufl(QH*EBD3w{^qgo*n-%4U|b z@4j5ti6`JlFlNR-h$_k1qG>b4M`a*a*YU}Aug6T1u~KXHu=#sRNRKI;o{EHHe#z1+ zFf3=tmRw`ctezyyJyBv*VM9c#3ZU3fvaBUdME^$FgRZOBu8%uad70(x9_Lmwl#srK z8?r_*dqNt!B$i>1Ur3y6B6PBksR zY-o=#AhFp0&4+Q;=zP8Rlx(cHom_C1J)`UPj?ZZw zbOYhvg6itkgjCY39>SxW642c>?vnQM!M@LKvZJ(8e#k1&O5M=Ys5y%3Sfzb&S!ubj zP)o8|vtM-N@T&8L)S4 zU7MvW`Sljs_b;q}YI0|JpjkQLT_PA@#Gh z8k8XGcXvY75v+$5gZ;yHD1w^ruIPmAnGL!U`~juX?Jjfek;!yW(A&3IMqvMTM6#kbm_OROD33Qd^>csVFPJ=ck!RekVI>nw zCE1#qNr=0lnxOb{)3_}vo~zB3yk@5>TJp1rcP^+~YNdi0#)G_IYqFcY{THRhgC!g;y^~} z-z&i#NQt6M4W=OvJqZn{lOa=nQ3>rJ3jWXmsGC9iYpd@1HEQF}gfyx`ly#LWj}FkNfYbEQiIZAUG93`N7510RDf#Q( zfH{2^1$(Aw+npj&N4^)wDl;_lVMzVs9g89q9=N1-w%5_Fwv9PwW&S}N=U z_nN&4i1|d>Ab+zV-I8?qmoyURl;e`J{7(aWfO<*@c?Zgrvq`7T$!iIjzyMJv3YZKz z`RtQNW4ix7adz$QrTZqPQ4D-*cLycx#;aGD#-+R#=^=L&))=#^e>G^M&(9Q$iB)R^ zUBNj=a7R14JQ{>p6)5SDC*PWW2>%Tl&VVma0PLfupf$y)1lee*0 zU>=D&MkD^DNJkDluv5nXcukd-^KTP{Y>UgB+oP}a6`M9JzP7L~Y2J06L>nu;_T3k(gPHdRW_j2lXkcH0f#o`JJXOl4*sVJ3)kUs3){`RKBwpnf2K^e# zwn)7X;Kit8z;D@q2|Ypt%cK4GkNAS&7#1nf9JRy?$Eauui62wF$G9~Iw|nJji|+2T z9|Qnj9W$Gh`~o#eRKGh3fqz?{VsBCtuA2FqiqZ1|KE89HN0Q z&3_BeVp{L{qlKY)v#O*GGVz(LUz5>}drrmbdyUq#c5}>e{V5o;&~Rni<4oEAnm zjFrrfE&<6k-9=A#&_+WlUE|-z^C$7%*x2|C^?y&i|8E;k|I-1dS6K#VMs8>ot~RrP zl}m@tJ)rwW1GSo-Rm^K~Kzu&MBSRA*y92%RM-fN;rpm-Md;RtjYOtiQv^@VL{}R}u zWM+UnG6KlULU^!UfUApf2H4VoMbx1gBOS20EPk-jZ}sEt4Y#>O+xE5DPVA7eYm7?# z&9?7vI5W%3%csoGlH^0dl=hwqnwP@7a;}*a#W(Na?MBP3%@W9zwN8OJt~6f`8ep4c zifaH$q7{HyL|{r(+QLt+*8^Ky+iP}t$aN`b*PHGi#4^Ap4^%4{Q({{bBMRsxi2|Fd zz|RxTc|A|;>ljdNZW0X6(7tUB+EW}ZHeDBG;+KKyW~M|2zQ?~Us0soOeXKgU8X~M> zAUv2*et{)uYGurHhPaK58@}f_Y6xLkSjR#?g6vrgUp0Fl&fdrdjGZCayxg=%prhS{ zY&qg!kvL4aLTcG4@OX(-doA>`&vhjyusM%c2br%}U=PdoTbcB%jXUiFBh?3J65LyV zed~aP0XA9j?O)H=k;YXHmlDkU?eU+q;?d`*?U94;V~uK`eu8~^cDPjD8J(eJZJHS* z$rfX>rJ4yPrjH)cK|Z21z3Jixd`H{)e#)@0+9yCw==UhKWJPd*$>XfvnM+1}d3Fr@ z(JJ;wzgW&L74Q~$#q3qPOxD+*!Y)2s0%NR$;>$6CGB@lmd&ThQG?KSy3&1iMVj5R% zy95TNuDu#qQXe#UujI?0WMYLX3^`6AzqQ{cTh~b9pYyu!F@hor_)=LKIX+irE5;-kb9E^OW@BnUDA4SBpk&YIRKCjLY4!jx#!l{7ob5Kha+3MsQ!y#^i$f+aV>`3WpAw^8ea z*@OS0q&8=uXo?2f7c8CbA>*EXbr>)P&t<&wnewy6``vN>{e1+?%(^V((Pfu^{Grkn zBMU9o5wBrWC}g{#SbQK?vqywhe<>9c&7XypC;a<~X$ws(a8vIan%1~X&OjDe>`+-* zA9O3&t^aNM=0A(${=d1G2+!0BD84Sd;!2BPc5M^{E8#q-Enh)+#~7wx+g0$W2WE zc3he$gDnhO(c`L)rC5o% zkC(-tJ~8o`^l}4U)``?Dz9KMYRRI`GjEhwjc!QaN(2Z>*m=ye z2-Up5)>&};RG!#DpyUPdjOV_uPPKc$S-Gcs46u`uhxs2(NK2+pkJo?J*r$W&Gnx-s3D>c?7Sp+NN|c< z`xW9$ltBBhug|d{2MC;_&Cx`pz2wvLM!%fNPT4s`TSG{))@+!VIE1%6}Rh{Fgdk{-5fK z>>&Q84fMizjVsiSq9H>X@C9SIlzj-K(S=UtS32w<@?;m6i+X!Qb>^o~PxNeO0;;HH ziU~OXBj0c;DOy`w$Dq0paJU9=A$=6IQar~hEe;H~jL`507kf^0OA$;6B_O_}%bo&| z&;|59d=zzl4EePC!00v>1Je{#-26*b;nKV__wlav!wC?y%K&(u^8$D?8h|_$U(LjY zVY2||re+$jz;`Wy?F9|)xmFsc?fW{5{WNfUbd5V`10fyw%e1a~_TGbpbby1T_bm`X zn|TS<5c+-oAC(oK+OcLK2c(anbf95+`0N{03?TbU!>dS-p^bkE0{Ps>l_kBGYul>_ zlAla9`t6j}gNNZ0fUiEIZ50_di4VA%#S!`Ps=KoE=7<%k*Bw0L8Zg}~yz{57bQpAq zaRyHC@10?183AO+Pt^JO9gv1m*f=<=`N5)r$Lrp5>rSX!z-oD<%xlUUeA@ldgExv? z)H1iUA3Z8-uNStG=T_dHBm^@iyEglgDPYHlZ{6RX+G;kCX@qv>r#@g8$^{l0zdB#y zj7atSp(^`+n<5nT{55be`A48i7m*J zB+y30qUQ6M8a#1thfM(*CT{d=qQEw6REotdc~cPSv&f({{>9Q9-d zTKN7`Pr{T;iK@)y@N-u#0NZioz2_%EC7wT#$!u0f6yOp;H0LgL#2Rq=l3yP{EgSIBc$eT6DBGQPURr_y?v!GExr9aVUWlY5N;x%EDKbilg zx-t1XXI$Y0)by_L zGCQ3|V|e?sVLYxscYbyebpsgIB0y0-OGO=jn4PUftzGu)*07(Jf+k`;loWlvI!8p4f#2jZ0Pd4_Vk|Q=Aajo#=fBFJ^sHG-R3a)0 ztx-^we(&t)$D<#j-F7PxjCTU0dEM z-edN4_^$vbT{Zv&jAX!h&&2_5fqxISMGh;LPoU=jb;UM?PM~?)5(l+!Gl$9rd;b92daxrrd|&%t8@MTb>&l%ca?+p<6a=_^(F=6n#Bl56 zw1U!L>2sCSxqqM%bQ=Y@&^+S){zrm@mbzN1qzJ9+!eoPwC+PVpVH0&SoPv`cgYB>V z&wJuc+^2t?8@-OSRgq+sExy#cVb?YbN}B5>~rAp^UvG>N+uY!;P zkN216;B#kV-&Za7^81_lsx|J6tP?RxlfSkJt`s`Wpy{EA{QN~{ zoXebxZ3xH9zXdRBsPgKIhW%+6altvHBxhRPYPVsT=jiQc!CRVams59xZlxW}9auYa zGUo!Eaaa)q?QU>S)a2*i2KnUb!QNIJC_W3zok376Xwm=k>`=!o*x;?8C+(v7-q>&k zAUeZPjT>lH3E%VAoKiOHHV1s?G8(Li#va2XoA{q~#Gk8D>r|-<{Zp@udNjCInV*TQ zc(thz(X=7V9ybqsVw^kYkjy>i>srggmFAuj4w#R-1#Z^%oTWLFNKaq2fWZ4|d;MoY ztdoJ6lSebloRKn?H3H|N8AG;{%PZf+^gyHv0JzZ*gHt6t2bichfB_Vd!3Q(I)mhxX zS>^^z!7&Jt{2dufUDS>B6qopo2QRZ0E_W#*qC}-!Sk#S$OKPgih>t-tRp44YA{wG* zB~C4jlvx$U9(G(LfY?xiWhgJByH@(>3NY|PfQbkwrdBAvZA}HZ)+M(Fyu+)Ja~XVY z>TGqiRm>@L7KxQ4cpq6BI@oyb9`ttX+7 z;o2x#7cH|terkQbCk`w-{^Gj1B)g$PZ5XB7oN{n0v_RQO~g89y&rL(M3 z1AfV&mdS)Ei7BQnWWzUwDgRwenRO4p*VIkXVBpPug`tqQ1fTx;LXzY^aSG=Vz@*p+ z(A%t1{>wWyKR-%3@-S10oFBPGedX?T#-|dXeS2M?>h#5%t)MsMU}3PpKv-C4yfWew zijYV#t#eY)dPG z4}HAtmP=atA=cbKM{GN}F=dLbRPo&1Dkx(`epU8~PZs_WM|vNA2l8cOem8)M(hA^1 zH&8~;_25+#_`5{9>^>ep4X`jYkMJZdn#;CPuMYOi7hMfa_Cxe;5@?7bi{KB?7PHo$3||c)e%biR({6~ zIbfiGse>;uZ94w>hYD9} zptaNtZ&OC&x<8F5JLua}oc*wbSKIy+75nOkkIdnwa-QPp&m62aRvgEV7kVIFaszp% zct3B_3?Qt}GNLkY9`2QAKc}*;8J;EdN$wX#9Gpap7s=1ne!fLL0<2j-zV3^b|bvW(*Q?QYEhe}$6ls&Amn0?P7>zLHYK@Kj`en`)*?E+6*2ewXEl&PD93eFf*^ zZN1~n@jG;G%LgHm@wB6S8sh8NveALjvOsc%-I!aHuSb+V@qI|+RZTcCWsBeGMfdg~ z#_~xI3pP5S`Mj|jg1K`ljj>C%SJcLJl3%fLIIuqaDua5DbGwJJM(^$P3zABF7Qc}<7z*hD@HUltQ+ zOk(5XNmfYy@LrHG+g9dlzmlULKBx|0R(l*98z^OC`Ax}U2lHU%kq(=_br8F{DWIFo zPceKzrxv(0{Ftc|V9(Ttpo0v|cRCD!^g7fcjpeh5uY5{G@54m-2>eDL*j9t1nZ@`a zsAw|MNPlgEfPrJaH$aCYG++)~_|pxO{swo$&$`_8Iy?nKE=K0>xIf(G89May_U!3% z)qntkEk;9h&HvBY%WM1%%jJ-y=~zYhkJMAIDaDQgqeP@E7zt}pKoHmdY54i;e?V+> z^ymPsRKPihd}!EX%JU?KfNu15sMY5&n+G7NP@s8vX8I*tEObVMSCc;}66#U^{3sJT zIKXSfvg5}e3eFCwd&D1_!FPkVTiDtGTgd{JZ{eVuReSUgnM)zZS!v*r?~;Q9Uc({Z)pb<9Ix+kEu2ul)Sy+5#C(eM33dN-Y_ifo z8>lktf%I4fpf{`!h=G^s%??cmq@hT6^&PjF*8S0s;i{o>nZo~&sj#3!GX}IJ2|wFVC!< znY~o$6qvIAUBSErl040M#Tceb&red{4vbB9spmEWnbzNq50FMB!}x6?7E<5)faK8V z0y3^3$i|leY=TbZ7?!OC(q{%ZT*VDILK7N(46w@p{sYj!@h|y+S1b!KsC`wA!;MAX zfgC>r!9hTNtb0vpgAIy+sK`&v_C{zBc58u-ZUqFQ7&=bHRu-QLIskNg*zW_L&1hV8 z5)B9*ND`M&0S%P$!y_NG`%!Emjr#J8o0t&F;Q=t<1eKDOp`$Rl17&%C1E1wG)EsK@ zA7#IXOjUGD7{q`5$EsHmaeE_bt!mgT3ThveP=i+FxpNljeMaCA7qjaB8;LdK>y+>J zw1rx%Re{lz9@5uYYUp9GCs=kc*LeV?PvmA~SjKpDCJ4ku2{emCM`=91WCWIuoicPu z`^kFXPje_M_$Tq(1vGn^NE0E31q91kFR&20Iu)WAX1yR^fUrGB0wAIw9qZq+R+Yd_ zAo-W!=tWIwI)Ij9E^3H>YJbWu163^x6LlUeRaXHy7I;YET~lV<0Ug(dR^#9}9q4o) z5!Qr7)_49_-!+nkjNO1dI=ABz2UG`>7j)o5>=!`O9H7sxLEe`*hd&QxKF z%i}im1x((0Ik-1;R*;L>8Qo*5 zht?(-(M5HM%0neum310G3mip`bx@865C8USh-I=2WSlR?9u z{%AE|Oa^t9xT*d~E}s9h*e=4l1@;#T#7pShBxtV>dlYbU{e$iSEE&W}0{`bf?{{Zs zh#8%q?LhU3|6!aGB|PSS?^4P5Zx~t`KNhq2_M(b)4=5-l7p`SZof4+`!>ZZ*2OHG! z?0M*zi$fL>@Sz=Xy0jLiO$>II2cYF+_ZN_vUMB$A?6H+Gu(0F@FjN9QteyQFZkzs& zNx({#pbvhc7FcGm`|CZ2Cn|v|S%(yC2+nFU2cl| zLsESUEnFIMlWceUG+$I=_4^Bw$evdx#5{)@1%ZP0Z)_~;F{m_8v`WDXDiv^e%X;)V z*Cvv_*Bga@QlR_6aa4^dU=WzbLE%LQbe#g4D-H!rx-9q~Tw8`n1067-TDkFK4mh~#goLyD4E_JgXK)zKf zwtWo|Y5#)o>VJDd$KGIsCPNah5k-NDW{&-7rgK9>n{)**xUr=U721l{)uUCe(zV>y=rweVy z`Nmt@P|IvWbG#9Ih?L8&UB;JL#+(`C;mx*AgT-Ef0RtGzU8%DjlWI>c0{pc2bGpeV zlGYjN6-b`L6qg4tXAQ4g^{~{uVmrrretwD^z<$)}0Xo2A6}1&1FuB=aSRJT9+BI5d zY=fJpragbbMkB-9yjr33;#BDyl4R|E>oIKOcnI4%17DsQJb>KbHv2pz#=FS4)hm}a zLY1%=KXvL$0MG!nxT;4At$+qUH3~jNO{-?(eX_F;%Cltm~X2y~&A!E(@;6_Fu#nX-f2+$q=xfK9e1^?Pm-}_0Ox_?;Kpn3L z_5pU~8oLYUs!kUaZWM7H%vrPZQT$xQ!anU@PT`x8tUzbp<2D}>#k4zUmbQmy*K$nE zinf+vT7oO}m``OaEQU|J7{egXzS!s?KGnl^vH?5s0cW=DoKy&({GySuvsz=$&|@pp zsfkGUY8gk_^mgj9!QAS;I%ZGfP#Dc3P5r{+SoV(q(;(T5Nmy>cuagdKyK1tn-P6`n z(W${SODAtp-hdflccIDvwt$L{btl8(M-6{n(1W9x30imIuOrh#lT{+J3~}DUaQezc z9l=Xp3*xr|N_cvel@wEMq`l-&{9@Vr-PK1ZhE-=`^wfa@i zLbgL8#AYJrT5;y5^(Gj0r09B8GJF3lmxEJSTvfHU@=j`HPgY@y*ax|1+Z3b&uccXFC8Su=@Kfa zDC#Gz%`{DU-d}jyoaPq#q(D@4g~b}EEdO^@{NR?q>{4$DN9X?Jt1%Py|D-kjU32eW ztrSA|iD&WC{}pn_<1QSn{1drLKsw$hF3N@-EKvQOkdW8;t~!f_To_OQfCO0MydK(` zEYN;QW%t$!jUN$Q{W~}Y)AvGnv@!)d@lIbijD4xE3*~VZ=D90~I2ubi+H}(!FvW3L zZ{#?i2BlHL#6LEeFQ6C^`uxY=+lt2aoN0y}XJiU`%Q7z@P2?C3xb;ir&p|qLL4)(x zfp^2p3doukR^>L)A@b5Tllf=0gLh^>i9B?*!PuhnkEDmoZ~cyh`2B>5{+7r+-#}Jq zKXSe7_e%foJ`4`%BUg}?Ek-}rv<_BY6U>39BG)z7zZkpx^A_uII;ch&G)~*roZ+Y)#5NJ)E;Lwi>n}p4~)?gibDtzL%E;s6KpILm4-;yj0L6o zzX78oAZNXhoD`_4G?d(F4os;7r)aW@rZUee?1D0fJ)Z+tR|}Nqe}SV@wy!$ww?5$P zi&K^H$L`uF4yB}~2QKlcC5sE!T2l;6U4A=Ijl@zWeK0);63o&}FM*h@>U99ZU#x$7 zdccjLwomcQ%#ZIvKhcqc&|AB*D3K_KGnJR`F_OW4Kw1azV+&_x?cjI!!X)LI6LG^M zwg^hXyly)PRiYw@X@|GBgZ}RN`Pb8sKYDKaFYjSQ1iCi_VpORDw@gMFwnHG_JzJz)nyCKI{^V{)!sA@(*OBZ6+UyBXgWa&o76B zUX$Tl5!#4aVPNV$DMQtMe^%%9&5!IxbtM1ixSA_e^TZfFpf8f_WVLf2qBigR+$ zTRKXZ^yA$osG~HYFYMCg;HGgJb&FasYWeHpd=vsmmmrD3?Jy z;2}5Cs!d?-a#Zn*(uR9zxhz1=|C8&;7dVz;`N&433lfLms}DT_dvQz|wpo4oraxHn z%g6D1Cd?RcDU6mF8sax2GzO*L0LQk3gZ{x`ya%EKkr-$kB_S#KlN&8s-=w@x)fJ$p zRPMd$X_L3n6!lt?M2gKJ+O~grQ&@*PZ?gjy=l@MEOS9?!5~l{=+W%Ot^S_6r`@i=N zOYutbNBd&<`MGD4YH5h%ce0+sA+E?FxPa4^J@( zywiMdE#J0Zuoqo|n8il?8QPLA+{=7tsNmUY`W?Rl2V+D9q@ry_!$!~p55t^$)6UZ& z4{K{{q5LVOUYRQqYnURPtl;%H;d!HkJ%Pm%%qzX`n!mM>ulhBzA(-Z^k4?K6k&}#3 zRYEhrZWs|90NvOXN-eLiba9BZ)B_tYu4A47c;V_Y2%b7cSw1`(luEOZ4UF>sIa;-Ub-wrCnZ znf7%wf&Dx%bYIFEOz<61xE+Jo*H{;_O(8ui?H_c2qf762ixO|hwKlO7sPd7ol;BKik^wAT*v>Z0d;*enC98C|TlG7qu&6IVMRH7LE)S zl+Tu_=hE63ZO=~z6&k;}@k#JlsfQBD1Wk16A`gc`q4Xv;W4T$m(fCl))NHo>i4z7D z_f9WP`ukJsqbf%&>lX5wAJ*`551a^!s0aT^|JfY#lcwYa1#z(-Q(8d`mr%R9>4Uzf zhU$&#Usz{v?I*29{g@zmsaQNGEBM1ZW3a#j-ut@>aV)ZwS5a@iG_j_YT!Gc@|2nJ&6?eA zaZoMgwpGKk3U>n4(2zZ`RGI&azT)4<1X1*3B}G0$u>f6Pnx5Azhjhj{@g5j<-ITwJ-@~}iL+e=WDf3ZD zP@ZJ$Ia$yf4E79hK+dN!p z=J4@07fg*LpQ~>?Iz%O*o(OzD)eEy65kS=Spm%AV&AsgkuBARy_R06S?nEPBD5`@| zbqW=W%D|gtbEINKQ-BULZS38*G0q>f0Z}mjTph^0y+y^_bHOtQc>b!P0?u}1LK`Gs za`jocM%N03?T%GeBjXiPb`j@$=>~-4WjoaDjywTn+Fh@?+yajZgNd-!#-h&Zi+T%VB0OmvHdNoWp}w{-t4#qK87QJKRcF2_mo_; z%Q5i$dR!zBrudNs0fna_j@-TvWySN}QHFrtFTP))LkNUax&d8LvW1)p!XE+6!y_gh z1^f-jlbdC|SenzT2~eNUeZM~m{f{MahIttFA_H^*H0^D4kqr|fy%7Z)L0l!c48}_4 zf9)Lq4Lk^!Y`=lG5IM*^xptAO`MihkIzRzj1?;r|v?Vu4LZu!@?ypMDdj}OB6G@)n z=}-W^nNLJl{w!te2iP_7d`D->1d1PVO!p)Y=z^i*9#?cGx+AaPQ+y;VrKkGDTU=ZO zV$>0)!H=!cmRlq`4~jVt!Pb!T@_0kfKlXgj@<(vF)@tZ+ve8q+%xDszDiuy(N;4A@W5S z(1A?WcJ^cU3d`*Z zMU|T`0xGKxEg`s!jEtTw4pT&VSlZ5Pd>z#ixmmetz`|l`GPuJ_*A8+?Y{3IFtymD$ zvELRB6>xLapevC5QW8no)?IdeyzSI|^9^>;%M7VMNt9`1rGEkTWfPjr7dw&{R2)L& zrCBjZIcw#-cOr#>HvGBAH+l8)XzF-g?s%-0&!uLMBkDdLd2n6t|eHB$K=eh zvRV;+J+miwRvtQNJ7?o4f3V0$?`t~Wxkf?u{+~(7+TB*&-N3M1q5bSgz}>1v8LVmG+pCk6A3P&($g%3~!qqv7G zH8HO>1jhB%! zM^~)`POvdYnm8;EhoZ0@W7=rj*CYpxBfz)Hu)<@yJ4RxZHuMb{K4L9*k_X4V36&4^ zCKg};X5Hk%^}os$xfJ=x^CnkmNvIjqJ6_GW?IZW%u4^D{Dz9YTZu_Y>qnyF1WhM z&V{vnKiXYs3@5%Q?=kspoZVV?WSTHMM-|U)JZSi?a@=*JAeMGjs|crA^1#rAarr_2o!n%B^Xq1sI3I}#>=tg(Ri{{v2CVFe z4a;mA3s%WJ8n2f2dZZ5KSCRB_QK}P7qJC#=8)fJ*>+V;!$Q7dDq<*Wyyu&P@Geo7K zgY(BUUCt1*Xb?!CoY5S9vl*Yan^g$a$6~EXLUYTF8sdL($wlBhGI2SjI z`z!?nyuqyV?H-(uJk2)=hL+Qc)DJicX@pf2v(A!gy{J(Xv^UOah^RXW(SN)7gf!v; zhG5`SXlEy|VEt4h#r{{=!O@3U_Vz?Ne^^<7J9j zc#f+#kwZ7Fr_6&_UiOQEadOrDf|lquSnPfcW?BZFo$>~P6~<<9CXaTpc#YVr$$m$s zeCgKY-X>0mQxDN)K>vm;_6nF?#zt3J&&+#e-cAJxm7LiRrw zH=XX!KI1P$`iEt)Zx445Rl#TpZWC{9@x96IHGp+I%}cX$9azRCydiG3Yv5)hg!|*Z zM_T31z~V$5r|C<+y()3u$8Ue>>K5Q;6hoaf3Y%LlvDC7h#3h}Lbea#s21aHrI0C~W zqG`BEBon!sT%rjl^NCL%(Lv<}io5tr-_4CmdFrgW2__5-JPSDK~(b>ERQN?JsL$xWLmz0It-Z0p$0O&D}4OTOR9 zsKdJ_@MjJ`Xv!xHOao0TJ(jvt{kE=P;G+iTb)HSr>W5{CPa}Ko9-bePUNy;@qvH1@ zSC@fZE7nJu_%zODM`ELKsfX?=TV1p{TY*C(9){kK?F_1ST^-m%C~R<1z668;QoKo( zf!NkXQ@WFND|okzCcsM9qL-Rn=U4y-9~FAS)NphJh=3$Ia9`}K{+9k}ibJ7K7&zn& z>z^KY!G;R^s!~oAT``a16(MCroB*6zJ`SOV854|~Hj?B{2+}T@S=52rIE4?M69$L~ z6y!8NZlaGD!*up)YI^YuI|WsW5O2uRf!vLE!(Sb_gw9L&J!NAPKXF+LTN3Gd&*(&| z5tb-lNW^~H68a5@+_Zr39%qWT+;_)|QIU#RPmLU2bmEDp;XOYG?`FmYi;qG- zX7k1ij!vkZ>F$|WbA8KzTx){*fE>A9UqzJMj=OR{SSRZFp@|*))O|FY*9#@fsasQn zTd$u8_m-u_lC9-p5F^BAxdo*7qC647Um&p0ELqahdL|xcWe@oT`9j?gRZk^_QSt)> zpnxmLrkSfyoU5AUUOGmzg%f4Mk`?SEzJ*tPY0B&=pag9e{K2teQiO)`Et_@Mt%Zsq zEa(a~d5a!TQ0t`n1V7O(L0{O>(klDTWTU~lmfN^i!>FB#UaDJFyyqkmR8tsa>YT=X z&&T7U;$Blnku`{+DA<#ad=q1?nsVMf1w7GXoi~o=1s}(|P_`dS7Entu`t{Bo)QH$i z623nRyu5${@ryS4@_+%>c?aCY^u{QkXJNK&4Nyw*%w4T0Z3*W9X#!8U)%hcE@XW><{XUwzNi0#6+j?~fYQ)^uaA#gN zNK{)RGZ~EPXT%X82~^zHW#X*D2w{$$nMX_heBbLfK)vIa zn#huC-v0`4H}&s(lem_zOtJ`7p1BsQhbMfBB#!7Ov7&qDtiT~!&_!T+-6RRg4dh*R z^cIpC$aUCDsF{Y78a|r`TGJcA^%GuT^zP%@4e{F(?+~oNSw>llVw+xkf+ph6-`N~b zLKuuMfk=OH;P|>05o8PVB|IJmo$$Ynt0u_-GKumSEFw*A1)vkPu6LtauWq)>XK{Yw zgYpvng;C7`38?~NP`1Lnd<30MVAU8~RHe+PIIIe1ij@Q)EC`f@+9qvGJYqk64Lm>Z z6baDNKDxUH%DbtyyyBeNTx7t)8XEIM_jTSv2=guzz4EZ}jkl|^%6Szkz>wMprU5V0 zJTP8OjYgIT1B?@B4+lsgX#A%9MoiLaXwc1gfYb?0O{Ysl;jDi_XPPn5p;V|#L-zZn zyo4L2%lb)T7~Mt*n0lnJXM>hwqWw^O1?uDaOFu$|<>55gCf%-@^zx;YOEbpTeUB;z z(`EzMsRNQe;2!??_+1os1eD~nJP!q7dRQjJ;JDk1V{{kSU&g@Cy5ED>m&}0w=HdQ^lM5=86jlmv Su7oT4-|8o|Po%4v-T4oXNjz`> literal 0 HcmV?d00001 diff --git a/test/test_gradient/save_gempy_data.py b/test/test_gradient/save_gempy_data.py new file mode 100644 index 00000000..78f5dbb5 --- /dev/null +++ b/test/test_gradient/save_gempy_data.py @@ -0,0 +1,70 @@ +import os,sys +from datetime import datetime +import numpy as np + +import gempy as gp +import gempy_engine +import gempy_viewer as gpv + + + + +from helpers import * +from generate_samples import * + + +import warnings +warnings.filterwarnings("ignore") + + +def main(): + + + # ---------------- 1️⃣ Create the Mesh ---------------- + + + n_points = 32 # number of points along each axis + + x = np.linspace(0, 1, n_points) + z = np.linspace(0, 1, n_points) + + X, Z = np.meshgrid(x, z) # creates grid of shape (n_points, n_points) + + mesh = np.stack([X.ravel(), Z.ravel()], axis=1) # shape: (n_points^2, 2) + + # ---------------- 2️⃣ Generate samples for the input paramter of the gempy, output and it's jacobian ---------------- + data, total_time_input, total_time_output = generate_input_output_gempy_data(mesh=mesh, number_samples=50) + + + filename = "./data_9_parameter.json" + c ,m_data, dmdc_data = np.array(data["input"]),np.array(data["Gempy_output"]), np.array(data["Jacobian_Gempy"]) + print("Shapes-" , "Gempy Input: ", c.shape, "Gempy Output:", m_data.shape, "Jacobian shape:", dmdc_data.shape, "\n") + + print("Time required to generate samples for the input:\n", total_time_input) + print("Time required to generate samples for the output of gempy and it's Jacobian matrix:\n", total_time_output) + + # ---------------- 3️⃣ Save the samples in a file ---------------- + # with open(filename, 'w') as file: + # json.dump(data, file) + + + +if __name__ == "__main__": + + # Your main script code starts here + print("Script started...") + + # Record the start time + start_time = datetime.now() + + main() + # Record the end time + end_time = datetime.now() + + # Your main script code ends here + print("Script ended...") + + # Calculate the elapsed time + elapsed_time = end_time - start_time + + print(f"Elapsed time: {elapsed_time}") From 18a4132ef3b8320b6cf7e636565e52467b98f3b4 Mon Sep 17 00:00:00 2001 From: javoha Date: Thu, 12 Jun 2025 10:10:33 +0200 Subject: [PATCH 04/12] clean up unnecessary test. --- test/test_gradient/generate_samples.py | 211 ------------------------- test/test_gradient/helpers.py | 81 ---------- test/test_gradient/prior_model.png | Bin 29187 -> 0 bytes test/test_gradient/save_gempy_data.py | 70 -------- 4 files changed, 362 deletions(-) delete mode 100644 test/test_gradient/generate_samples.py delete mode 100644 test/test_gradient/helpers.py delete mode 100644 test/test_gradient/prior_model.png delete mode 100644 test/test_gradient/save_gempy_data.py diff --git a/test/test_gradient/generate_samples.py b/test/test_gradient/generate_samples.py deleted file mode 100644 index 9b3d6372..00000000 --- a/test/test_gradient/generate_samples.py +++ /dev/null @@ -1,211 +0,0 @@ - -import torch -import numpy as np -from datetime import datetime -from torch.autograd.functional import jacobian -import pyro -import pyro.distributions as dist -from pyro.infer import Predictive - -from pyro.nn import PyroModule - -import gempy as gp -import gempy_engine -from gempy_engine.core.backend_tensor import BackendTensor - - -from helpers import * - - -class GempyModel(PyroModule): - def __init__(self, interpolation_input_, geo_model_test, num_layers, slope, dtype): - super(GempyModel, self).__init__() - - BackendTensor.change_backend_gempy(engine_backend=gp.data.AvailableBackends.PYTORCH) - self.interpolation_input_ = interpolation_input_ - self.geo_model_test = geo_model_test - self.num_layers = num_layers - self.dtype = dtype - self.geo_model_test.interpolation_options.sigmoid_slope = slope - ############################################################################### - # Seed the randomness - ############################################################################### - seed = 42 - np.random.seed(seed) - torch.manual_seed(seed) - # Ensure deterministic behavior - torch.backends.cudnn.deterministic = True - torch.backends.cudnn.benchmark = False - # Setting the seed for Pyro sampling - - pyro.set_rng_seed(42) - - - def create_sample(self): - - """ - This Pyro model represents the probabilistic aspects of the geological model. - It defines a prior distribution for the top layer's location and - computes the thickness of the geological layer as an observed variable. - - - interpolation_input_: represents the dictionary of random variables for surface parameters - geo_model_test : gempy model - - num_layers: represents the number of layers we want to include in the model - - """ - - Random_variable ={} - - # Create a random variable based on the provided dictionary used to modify input data of gempy - counter=1 - for interpolation_input_data in self.interpolation_input_[:self.num_layers]: - - # Check if user wants to create random variable based on modifying the surface points of gempy - if interpolation_input_data["update"]=="interface_data": - # Check what kind of distribution is needed - if interpolation_input_data["prior_distribution"]=="normal": - mean = interpolation_input_data["normal"]["mean"] - std = interpolation_input_data["normal"]["std"] - Random_variable["mu_"+ str(counter)] = pyro.sample("mu_"+ str(counter), dist.Normal(mean, std)) - #print(Random_variable["mu_"+ str(counter)]) - - elif interpolation_input_data["prior_distribution"]=="uniform": - min = interpolation_input_data["uniform"]["min"] - max = interpolation_input_data["uniform"]["min"] - Random_variable["mu_"+ str(interpolation_input_data['id'])] = pyro.sample("mu_"+ str(interpolation_input_data['id']), dist.Uniform(min, max)) - - - else: - print("We have to include the distribution") - - # # Check which co-ordinates direction we wants to allow and modify the surface point data - counter=counter+1 - - - - def GenerateInputSamples(self, number_samples): - - pyro.clear_param_store() - # We can build a probabilistic model using pyro by calling it - - dot = pyro.render_model(self.create_sample, model_args=()) - # Generate 50 samples - num_samples = number_samples # N - predictive = Predictive(self.create_sample, num_samples=num_samples) - samples = predictive() - - samples_list=[] - for i in range(len(self.interpolation_input_)): - samples_list.append(samples["mu_"+str(i+1)].reshape(-1,1)) - ######store the samples ###### - parameters=torch.hstack(samples_list) # (N, p) = number of sample X number of paramter - - return parameters.detach().numpy() - - def GempyForward(self, *params): - index=0 - interpolation_input = self.geo_model_test.interpolation_input - - for interpolation_input_data in self.interpolation_input_[:self.num_layers]: - # Check which co-ordinates direction we wants to allow and modify the surface point data - if interpolation_input_data["direction"]=="X": - interpolation_input.surface_points.sp_coords = torch.index_put( - interpolation_input.surface_points.sp_coords, - (torch.tensor([interpolation_input_data["id"]]), torch.tensor([0])), - params[index]) - elif interpolation_input_data["direction"]=="Y": - interpolation_input.surface_points.sp_coords = torch.index_put( - interpolation_input.surface_points.sp_coords, - (torch.tensor([interpolation_input_data["id"]]), torch.tensor([1])), - params[index]) - elif interpolation_input_data["direction"]=="Z": - interpolation_input.surface_points.sp_coords = torch.index_put( - interpolation_input.surface_points.sp_coords, - (interpolation_input_data["id"], torch.tensor([2])), - params[index]) - else: - print("Wrong direction") - - index=index+1 - - self.geo_model_test.solutions = gempy_engine.compute_model( - interpolation_input=interpolation_input, - options=self.geo_model_test.interpolation_options, - data_descriptor=self.geo_model_test.input_data_descriptor, - geophysics_input=self.geo_model_test.geophysics_input, - ) - - # Compute and observe the thickness of the geological layer - - m_samples = self.geo_model_test.solutions.octrees_output[0].last_output_center.custom_grid_values - return m_samples - - def GenerateOutputSamples(self, Inputs_samples): - - Inputs_samples = torch.tensor(Inputs_samples, dtype=self.dtype) - m_data =[] - dmdc_data =[] - for i in range(Inputs_samples.shape[0]): - params_tuple = tuple([Inputs_samples[i,j].clone().requires_grad_(True) for j in range(Inputs_samples.shape[1])]) - - m_samples = self.GempyForward(*params_tuple) - m_data.append(m_samples) - J = jacobian(self.GempyForward, params_tuple) - J_matrix = torch.tensor([[J[j][i] for j in range(len(J))] for i in range(J[0].shape[0])]) - dmdc_data.append(J_matrix) - - return torch.stack(m_data).detach().numpy() , torch.stack(dmdc_data).detach().numpy() - - - - -def generate_input_output_gempy_data(mesh, number_samples, slope=200, filename=None): - # ---------------- 1️⃣ Check and create a 3D custom grid ---------------- - mesh_coordinates = mesh - data ={} - geo_model_test = create_initial_gempy_model(refinement=3, save=False) - if mesh_coordinates.shape[1]==2: - xyz_coord = np.insert(mesh_coordinates, 1, 0, axis=1) - elif mesh_coordinates.shape[1]==3: - xyz_coord = mesh_coordinates - - gp.set_custom_grid(geo_model_test.grid, xyz_coord=xyz_coord) - geo_model_test.interpolation_options.mesh_extraction = False - - sp_coords_copy_test = geo_model_test.interpolation_input.surface_points.sp_coords.copy() - - ############################################################################### - # 2️⃣ Make a list of gempy parameter which would be treated as a random variable - ############################################################################### - dtype =torch.float64 - test_list=[] - std = 0.06 - test_list.append({"update":"interface_data","id":torch.tensor([1]), "direction":"Z", "prior_distribution":"normal","normal":{"mean":torch.tensor(sp_coords_copy_test[1,2],dtype=dtype), "std":torch.tensor(std,dtype=dtype)}}) - test_list.append({"update":"interface_data","id":torch.tensor([4]), "direction":"Z", "prior_distribution":"normal","normal":{"mean":torch.tensor(sp_coords_copy_test[4,2],dtype=dtype), "std":torch.tensor(std,dtype=dtype)}}) - num_layers = len(test_list) # length of the list - - Gempy = GempyModel(test_list, geo_model_test, num_layers, slope=slope, dtype=torch.float64) - - # ---------------- 3️⃣ Generate the samples for input parameters. c ∼ N(µ, Σ) ---------------- - - start_sample_input_time = datetime.now() - c = Gempy.GenerateInputSamples(number_samples=number_samples) - end_sample_input_time = datetime.now() - total_time_input = end_sample_input_time - start_sample_input_time - - data["input"] = c.tolist() - # ---------------- 3️⃣ Generate the samples for output of gempy and the Jacobian matrix. m= Gempy(c) and dm/dc ---------------- - start_sample_output_time = datetime.now() - m_data, dmdc_data = Gempy.GenerateOutputSamples(Inputs_samples=c) - end_sample_output_time = datetime.now() - total_time_output = end_sample_output_time - start_sample_output_time - - data["Gempy_output"] = m_data.tolist() - data["Jacobian_Gempy"] = dmdc_data.tolist() - - return data, total_time_input, total_time_output - - - diff --git a/test/test_gradient/helpers.py b/test/test_gradient/helpers.py deleted file mode 100644 index a6cb7f04..00000000 --- a/test/test_gradient/helpers.py +++ /dev/null @@ -1,81 +0,0 @@ - -import numpy as np - -import matplotlib.pyplot as plt - -import gempy as gp -import gempy_viewer as gpv - - - -def create_initial_gempy_model(refinement,filename='prior_model.png', save=True): - """ Create an initial gempy model objet - - Args: - refinement (int): Refinement of grid - save (bool, optional): Whether you want to save the image - - """ - geo_model_test = gp.create_geomodel( - project_name='Gempy_abc_Test', - extent=[0, 1, -0.1, 0.1, 0, 1], - resolution=[100,10,100], - refinement=refinement, - structural_frame= gp.data.StructuralFrame.initialize_default_structure() - ) - - brk1 = 0.3 - brk2 = 0.5 - - - gp.add_surface_points( - geo_model=geo_model_test, - x=[0.1,0.5, 0.9], - y=[0.0, 0.0, 0.0], - z=[brk1, brk1 , brk1], - elements_names=['surface1','surface1', 'surface1'] - ) - - gp.add_orientations( - geo_model=geo_model_test, - x=[0.5], - y=[0.0], - z=[0.0], - elements_names=['surface1'], - pole_vector=[[0, 0, 0.5]] - ) - geo_model_test.update_transform(gp.data.GlobalAnisotropy.NONE) - - element2 = gp.data.StructuralElement( - name='surface2', - color=next(geo_model_test.structural_frame.color_generator), - surface_points=gp.data.SurfacePointsTable.from_arrays( - x=np.array([0.1,0.5, 0.9]), - y=np.array([0.0, 0.0, 0.0]), - z=np.array([brk2, brk2, brk2]), - names='surface2' - ), - orientations=gp.data.OrientationsTable.initialize_empty() - ) - - - geo_model_test.update_transform(gp.data.GlobalAnisotropy.NONE) - - geo_model_test.structural_frame.structural_groups[0].append_element(element2) - gp.add_orientations( - geo_model=geo_model_test, - x=[0.5], - y=[0.0], - z=[1.0], - elements_names=['surface2'], - pole_vector=[[0, 0, 0.5]] - ) - geo_model_test.structural_frame.structural_groups[0].elements[0], geo_model_test.structural_frame.structural_groups[0].elements[1] = \ - geo_model_test.structural_frame.structural_groups[0].elements[1], geo_model_test.structural_frame.structural_groups[0].elements[0] - - gp.compute_model(geo_model_test) - if save: - picture_test = gpv.plot_2d(geo_model_test, cell_number=5, legend='force') - plt.savefig(filename) - - return geo_model_test \ No newline at end of file diff --git a/test/test_gradient/prior_model.png b/test/test_gradient/prior_model.png deleted file mode 100644 index 54c15ee67c182ef3ec387893f330ee5a402f30e4..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 29187 zcmdqJ2UL@7w=EhgHb4|XK|nx7q>1zz6lqeWBTe9=H|f3jDTq`R6{$h0^xjKUN))7p zUL!cntTor1b3L!LG?b{%Fr9(HVALwg z542&h7?`?S4akAwdxA3p{Dm19406nxQA4h8?h=r3G@zJ&d+|A<)& z+@~%X2J^Zf3j5-3Txs`6diBdOLEqJ1;^1qhZBJTXjVoq**%rU|d{q*1u;kcCrB5+B z8LF;V;VuRXRqt$lbrGv-zN@`I52w5mbsqa2l_|41!b#r4P}&ic{rP?I#m|{DNsh&@ zZmM@6bV3LLn5a1L%zb$}x$a`wU;T(}jO*n}ajbns#x)|Yp5OH!gj0&jr;bhB3kC2#>&~ z3oE)K| zCX-4+2DeEDg@$jC8jHnUO>e}e5`K?zGlv6*MHK3EJlXSVmfSHF&_^O*naXg{LX`GL>l@&*{^;;NJP z^K$N=X;6&4ni$f)36}Q0zXAUJTYurJ{tKUW0#>JKL&RTN#@|%WCpzZ#E_16jM?>UzoQ|GB;hM?v^6<}EslP|x3@rg|vPir!nHr4Tfm+uV>BI z;_H{z1aSv7J+gbh=cQLF`}2_+>hu1yi2fusgscpe9B#fm2T`)YszJIHE4#bkB})s7 zc@8@g`gjMd%;W@IMs<0)d+J^=mKpP%(FDhld)?H8EnvCKPdZ@n)6?F2R0=O?Sk~Gc z9-G%EH+{;EYSukAo$537s#5?dy@(b=?sQ4h$cFLm;Tv3=KQ_%o@~mxH_NLz>4kO7T zn^OUV?Jrca=R@|!F*EOUH1f)}F%fi}Gh>b=F-6s5x2Xhi90G&ylo^Y7eBw1J&8y$e zdMK}T@Ra;sy(W_ufUgtZTdu%{GrN2)OXaRqVWi;d`sDCx&mMjgzLFKS<>@pQGVWXx18p%s%0?6aL5dd+ zeS;!(Fp*2qoOd-PvA$cs`oU%LYcln1P>gsx08e3?%0`+m2L$>F=c>G8+ zBQEx7y+=y8Y}P0fx=s7-(HEZgl85512<2$hXXMIplG2qK9dX@(rk7^+KmcO}d(VgF z^jqO$56DGKyV3ytQ8wo+|J_QxcaY(mag&`0Nv!`H^O-1cW^|*UUcmSbEjFx>=OZ!uMnH6ft*R82u zIB+h-qx*-5#M=$Yc-NpMNCRdt7NaJHJ`p$=^ADeb80e=b(Q}mh@Y>No(lliTjMEo{z{ryNSxv_uldkfE_V) za~=^AOlHd7QhCe2YT&w+onV#}@Hp$&C;HOsxa|R*PofwK{jcnKjVe`bSY`Rf7&_qY z2@Y5~)lpcLHdCxty|O7A4!=;*Bb!X<$eoGeF{jUtN*M5xGNsRns;z(Re?04Ki?+Zs zzxBDnAHyn)M%k;^c3Qb)A02?)m98tqF-a2W64%m6ji@Z5XT6To@UGWVN4Ht3)N5&C zIz25K*8Eg_D*M}<7eBq`=QT(e(;Ew!oe8$ji4;Db?MZCp=N4Ah-rDl~ni_7BrL^nA zB>LL_LDm!_vo>xo0smoGadG<5TeJB}`~J05%77yWY0R*WOMTJkH6EX9b`p1A(x}nb ztvOUf=z)pXpx)s0J>0zLWS!=?!+6?T`u(hV4SDs4hlBN!mMsmL%U!AbY@kDQ+^Hha z%)rG5Cz5kCmq(RdOEO_)yhgxD^UT0lm1DxO{GIiAgA&bqs**jpz~Tb26R?<WqMGx9^$zviduH(jHDk;T91WkH&{lHGvBafA$;`j}c9<2YnY zm5ccYr43DYr>YzVZ~Js8mk-?=ua7rat6yr$rX z3wVDrJGl*A1cQ2GGo7bzvb=MyNT5#y)IV?f)ZCMYzptIl&FttO zlZ|N6T`tvKExnmdSoxVKw0?I`dgS%P9^-yLGru5vw7nGRuM==lH?D6Bc?nV_}oO1byjn-OcmSHpNgq66~2*o zA;|P$+ZwudTew>~hDDQ6d5I@R$2^%@4{1Sa_$fj<$z8>HLdu@uS%%aW57j9iTI0_= z+3Xi^G~{sYq#G8wsou9IYI24RcO}%C?hTQyW+VF0Wah_xZNa#uV$Fq(;Ox;G+A*&- zn$<$_>$ZCFrrP(gGmhHzS_QbH<;9C!g}l@i%3nJd+2!hm()@+8q*k*D z_KUKJqi=RUM2hB2otRGxaA9>s)}sef@=9mIU*p?z>tWoNm8*oOqR?Bspga%~9uPfl zWZd&nc!1qmRfcycS5a>9)0-(~`wy{6&2amJ4s&F4RvJ~F-Cy}(T1gq`jTl&=8Wg^7 z=}??_sREn*3Z(jQgTgp^v+VP;UyjiR_I!GPn_vCGh4|Sbsh?oe zs~mPm6^_YMQrnf-K~naROpfQWjQF+Py0z)r#x~F2s0t|%)i({otVR4ZF*i&HT=;>1~oLc#HGqqd0Dv&Yh& zemh+TQp#@dWMZW{jty3&gWGpfnhUV{k^Xih#b%Z)HG&pR8ic>2xQn|jBRqVG?V zh)-tux9uM{)$?u`zRiKFVG7<8<`Qg#mHQ<=WZfp??u&;V&o;XIR+p2Vx=)#OLAzJW z*tLO1#!JRp!3%+QvZ{ewO>%&MzwX z$Km(>v-Cj}bZO`ZKI6GB9Fr>oD)6{GD!C*MAsIqP|J~w}jrG^vv%3SlDqq!-j5tmN zzy)rd@zNzZ$WgiXW@uAb4mao6VB<=fj zAKdIa<6K@;aD58}6ds&}^%&Pc&f#011#lu+QcVS#AqiRB^!-b!KBKfU8@&_j9H!B^ z!wLMimg;jKd8jNQ0{ELrO~3H{iw-Hg~X@*@lDPvH6fF%hQ+*JW|FcXuM?&_ER^2ROTWt* zQ`3gMiM*L}cF*%AJctc_B3QlU#*oq(OYPOSHk!&lWBom@XTE6*Jl9?*bMHfK$fz7q zVS2j`%DGvB+hSP!P_`WBfl0!<;NbpxPdeGZzQ3OEu=UEjWwm%BKF3LVsjwo$pv272 zbg}A94oOW$Ra<#0*g-E4pt7DY4xO$zfxzjrRkHl4TLYULX(E9~LI+QuhZzdsW^SwZ zB;Zfw7H{J2-zJGG>oDX%86D1DFh8TLPj1e!(y@DO)(MITL`wwqBTaT0rQPN!_f9^)Q;NG+H0Fg*T*EnD zb^p#@u#EUtwsD~T)rcu6-v3;dclMXLDD(Z9+s_owXrs5-y(pS+u{{DZB<9F{`}}Tz zV|t^7jTIKXUg!K)BZLh|*^szzEve=;;r7a^jYMH~jkois2#s(UO!$<&_@l*VZ~c|& zyQ+`MVKYx2eb5T3a)#VAe}|$?^|lV}gx&!WaT#sK0n!z#Gce}(f#`Ft(ijJ^SG>a-Rp+NWQ5Li(Lk4QOoDn_6UZk`~7#oX0AqHEM+MBA1X4!^+Vkc;Pul znCK1rp*OlboP*xtJHO*Sv}j-=kB&W4%oA=iS(JF4llQPcI~&)e{ixOphe3vsXew+DP7I(LH7|F4?E=Pa;T`sOw`{227%;X}vUU={_phV=Y zPaTWVPAcb7$LJ;hO8;wj1LtSdFKo#o>G^#M3H322sr;l-e@e)zDTy@2q6`A#39dXX{g|5 zLQ^?-H(|>lrzX;ngtYk5psA=9&ylGD+9IoQ7m9I~5!uWbIOb)e#9yvAxQcd)acrJqyW+`tcPzT*ONpsb3sQL&DHgH4a{1ZOF)i7aEadA^k z0ZIw{rfc;UE#YkoeVJCPzrG53F7%r0{Q{K6D^?l*);=XRQ~XVE=YcCS;FIk!%#1BB zX?r@~Zmx`@CVen04fuM;=*xl}co0(;z~s}=p)2o!r=iwUcweWH?zBr|CIqnN<;%b z*DAiG^W2K~woP!q0TYcv#qg;nhmtZuxde;(hTb8xFIkZVat{3N^QP9^N~&?bG=6n} zQr#w9rLcv{#aIesttaOZiLvBkFl^b`npT={vPoe-`KpD5?i2Qia&9b zh(%x+sPN-8PR_?1l)zBPJr^SxU}wSzo z97Q)eO_uA8jsYU_62+IeGvXRI1&SVYM!CW=(g0CQdDABtfKM4IpPaLHk$ z_sPS=HOG?HnH#DZ@_V$csJ*M_TQtN);F>=iUDdk0&XQffd6YTrl#lqK@NRQNM68U@ zicwR~yJ5;jIE$ZcFtus|UsS5k0)u?44EN@^6=CrTeXw^MNu)dLj_uS$qtn>Hg2TuN zw5aTQTNuE)8}pqGKt?*#khPjJGZheX7 zCMEO!F6h>6!iK%u<_~#n%$0sCL|T)3(6hccnYStb@sXe09%HS*(*B5sn+P$}#O~iH<$4u4} zQXF|xymrF+f}{Abm6<0M`fb#)aaNQIj3#?A+q9 znp8=phORM}P{K{E1ZXD+8|KY@U2I_(vf24|gShM3t~7j)B*qZ3nc_2($e7l}xXi^C zY+~j!fTE4bkMOl>uOtOGD@0*pN`lCpec0y>`XvF5n}t`eEsuMYEHv%0`b`EQ{qFeV zD;)#(_B-fXg)-xZYqm$hh51vT7jPdaVp5AqHNdhrD;EnApW*rD{VzT~3ps0C&X%Yv0(?lJ;9;HYh3{{N+uVFK%Z4*3L*X+m8x(qxjD3h*iR*!Om`{0dE|S%f6ZP!Tz@R zZ%(RqF4-}HOp{r9@ZoYQMLV5cvKrm@$&}50e5~bV+Qs1sNe=!da{fLg0KfZ5{Ie!T zd^t3@?&P=Kkw*1#J9Jj<2RVo~2AoSQJT^$$EC5UL4cbkXkW;Q6by#Xl5ehM*xu{P= z`>@_?R8W9RBkf%Ra{r)Z=TSV)wI?Ei!g!6>s9w0)t)azzT}AZ3II6_#yljfHuu7AlL`udfDtrz4gH z$(mqgO>4$K?P8*)g2ecn@y0tMIirll+h_xat5Yzils`)migKERhD$)qq1Z5drag-> zeP?@bVa#t&d?x+K&YK+5#_9Q*w(m-c^VitJj%MdwHV|D(9VOr#HIgMbUsAnfdHU4E zdnmyms9Sv$>djsXz2sP*StwqkRN9K_D0+b?x2NSFEyD z7qsZX5AFq$b_)s+gzA_0nA|jGHW1dun%%meP9JglyV3Co%&rvuE?P!F{5BMZ8Svxl zol1jy^szFswW=Fq*-HBxWGSm%W1-BwO;fA(aUlnmJmX@S4qFCG?RJsgGvXc7=mVGk{g3B*-)@A2YH%GU>y*Ga)#DAb2SWcO) zqiJDF+gB})dmXFU#xUQFV1%=+Yo*b+BUqcL^gZ(sQv9Wgp9B2zVp5gu>daf*7m-*9p=Dl=MDn*v2Bc&%kVTaLP_Eokzk)~vGsGO?LR(9y6fIhRP)M1b5D)e zxLP{+>9Bg;%a&qOpAL%_b)@+0nyxxm!?)*GWYJ~y5 zv+FZ6M1_B~7rd^;YuwCRjj8qsBvs6}eTgIW>}B1tAikw^^{zIJyB3z2BKPmn9y#jt z`*wbQKALbY{T#qnjO*r!thP<{s&&anVazd;#qb;TTa#AL;e6W3ky}Wn3~BZMlt5Jx{^gi=6KDRl56Z(kNMy)aBVI#IU5Vdg*&Vb7_6^V zW>U2sY1-oP@vz*mLdGoNw^hKv7eW_`hc0xIFyGU3Az`0Lv7AQh2|T9kBCfc#VJXoC zdjuA)l}BW$buSBCdA8qZy>)@R*p2R{UG(LxZ*WVbj6X`sW4<+IXT4+c`V*N`5D-|M z{+bZHH}JS7SL;({(^d9w6DQU!TTR7K+pl`s3MDEsp`G#TK$b~9PDbLfu710Ofyq@6 zURY^OkPa(<4-V&2N}Q&oy;tgWCMe$*?D{ek3JuFR=ExLAoO$iq+s6HI8LPFmvfQ)? zs$j~J?@pD0CGqWzzE1~KGuR-S2n(*1Bt-7-DzTh)K4ha5gmLPIzjRA^#|gzzJeTVn zdlrLo8p()l+;^#UN3eN~mK)`3F1s(Qlwnu{O3-umuIh1T!vkewabkCB02oYWINV(R zBlj2mw3yBPHYhU@ru74RPs2ME1EtJ@N$J1flJ?=>qam$rtY+n z{EO#rsz*~a5Dc~mY2m4lf-;-9Q65${F;PYMrZVfdQmb5kKcx8q&d05bk8|yHIKa`@ z0oTJ22Nx2>v#iTarNwIIXzaP3+Z)^IE3~f=clknBrZC>Me*VRJbCCW1cDbWf4PI)V z>8!13utg=tvQvEJ^S!=GqyBB|2(tCWa4=_c@Mc&ag_Eto1!FqAFP2#)io!R6n#G2&lah+MJSD>|B;0vd%I?=T;y=f+G3n6}+oSE5 z4g#YXkj&NQ#2eN&#| zdrf|o(q99ywZU?9ZI)QUQykYl?%A=VeeKS|zL3YrYLS51h2hDa_>qjsIk)jkoEdtZ zAN&wls;igG8`s{}l8cjPz68OwS(r}IKIy+HcbJh+TpV_!B$$^B!z?O>k`)R*6wQ2s zN568Ga=OSiWkm?2BfsxS*-J|+QLv60DKd83BQ)-YUlo1#?fGAu6FwO$HqMrP4N+^9 z+-X`(@_OR7a!F8r%k#h1t5;^)OH7j!Jm$J$LK?RNufl(QH*EBD3w{^qgo*n-%4U|b z@4j5ti6`JlFlNR-h$_k1qG>b4M`a*a*YU}Aug6T1u~KXHu=#sRNRKI;o{EHHe#z1+ zFf3=tmRw`ctezyyJyBv*VM9c#3ZU3fvaBUdME^$FgRZOBu8%uad70(x9_Lmwl#srK z8?r_*dqNt!B$i>1Ur3y6B6PBksR zY-o=#AhFp0&4+Q;=zP8Rlx(cHom_C1J)`UPj?ZZw zbOYhvg6itkgjCY39>SxW642c>?vnQM!M@LKvZJ(8e#k1&O5M=Ys5y%3Sfzb&S!ubj zP)o8|vtM-N@T&8L)S4 zU7MvW`Sljs_b;q}YI0|JpjkQLT_PA@#Gh z8k8XGcXvY75v+$5gZ;yHD1w^ruIPmAnGL!U`~juX?Jjfek;!yW(A&3IMqvMTM6#kbm_OROD33Qd^>csVFPJ=ck!RekVI>nw zCE1#qNr=0lnxOb{)3_}vo~zB3yk@5>TJp1rcP^+~YNdi0#)G_IYqFcY{THRhgC!g;y^~} z-z&i#NQt6M4W=OvJqZn{lOa=nQ3>rJ3jWXmsGC9iYpd@1HEQF}gfyx`ly#LWj}FkNfYbEQiIZAUG93`N7510RDf#Q( zfH{2^1$(Aw+npj&N4^)wDl;_lVMzVs9g89q9=N1-w%5_Fwv9PwW&S}N=U z_nN&4i1|d>Ab+zV-I8?qmoyURl;e`J{7(aWfO<*@c?Zgrvq`7T$!iIjzyMJv3YZKz z`RtQNW4ix7adz$QrTZqPQ4D-*cLycx#;aGD#-+R#=^=L&))=#^e>G^M&(9Q$iB)R^ zUBNj=a7R14JQ{>p6)5SDC*PWW2>%Tl&VVma0PLfupf$y)1lee*0 zU>=D&MkD^DNJkDluv5nXcukd-^KTP{Y>UgB+oP}a6`M9JzP7L~Y2J06L>nu;_T3k(gPHdRW_j2lXkcH0f#o`JJXOl4*sVJ3)kUs3){`RKBwpnf2K^e# zwn)7X;Kit8z;D@q2|Ypt%cK4GkNAS&7#1nf9JRy?$Eauui62wF$G9~Iw|nJji|+2T z9|Qnj9W$Gh`~o#eRKGh3fqz?{VsBCtuA2FqiqZ1|KE89HN0Q z&3_BeVp{L{qlKY)v#O*GGVz(LUz5>}drrmbdyUq#c5}>e{V5o;&~Rni<4oEAnm zjFrrfE&<6k-9=A#&_+WlUE|-z^C$7%*x2|C^?y&i|8E;k|I-1dS6K#VMs8>ot~RrP zl}m@tJ)rwW1GSo-Rm^K~Kzu&MBSRA*y92%RM-fN;rpm-Md;RtjYOtiQv^@VL{}R}u zWM+UnG6KlULU^!UfUApf2H4VoMbx1gBOS20EPk-jZ}sEt4Y#>O+xE5DPVA7eYm7?# z&9?7vI5W%3%csoGlH^0dl=hwqnwP@7a;}*a#W(Na?MBP3%@W9zwN8OJt~6f`8ep4c zifaH$q7{HyL|{r(+QLt+*8^Ky+iP}t$aN`b*PHGi#4^Ap4^%4{Q({{bBMRsxi2|Fd zz|RxTc|A|;>ljdNZW0X6(7tUB+EW}ZHeDBG;+KKyW~M|2zQ?~Us0soOeXKgU8X~M> zAUv2*et{)uYGurHhPaK58@}f_Y6xLkSjR#?g6vrgUp0Fl&fdrdjGZCayxg=%prhS{ zY&qg!kvL4aLTcG4@OX(-doA>`&vhjyusM%c2br%}U=PdoTbcB%jXUiFBh?3J65LyV zed~aP0XA9j?O)H=k;YXHmlDkU?eU+q;?d`*?U94;V~uK`eu8~^cDPjD8J(eJZJHS* z$rfX>rJ4yPrjH)cK|Z21z3Jixd`H{)e#)@0+9yCw==UhKWJPd*$>XfvnM+1}d3Fr@ z(JJ;wzgW&L74Q~$#q3qPOxD+*!Y)2s0%NR$;>$6CGB@lmd&ThQG?KSy3&1iMVj5R% zy95TNuDu#qQXe#UujI?0WMYLX3^`6AzqQ{cTh~b9pYyu!F@hor_)=LKIX+irE5;-kb9E^OW@BnUDA4SBpk&YIRKCjLY4!jx#!l{7ob5Kha+3MsQ!y#^i$f+aV>`3WpAw^8ea z*@OS0q&8=uXo?2f7c8CbA>*EXbr>)P&t<&wnewy6``vN>{e1+?%(^V((Pfu^{Grkn zBMU9o5wBrWC}g{#SbQK?vqywhe<>9c&7XypC;a<~X$ws(a8vIan%1~X&OjDe>`+-* zA9O3&t^aNM=0A(${=d1G2+!0BD84Sd;!2BPc5M^{E8#q-Enh)+#~7wx+g0$W2WE zc3he$gDnhO(c`L)rC5o% zkC(-tJ~8o`^l}4U)``?Dz9KMYRRI`GjEhwjc!QaN(2Z>*m=ye z2-Up5)>&};RG!#DpyUPdjOV_uPPKc$S-Gcs46u`uhxs2(NK2+pkJo?J*r$W&Gnx-s3D>c?7Sp+NN|c< z`xW9$ltBBhug|d{2MC;_&Cx`pz2wvLM!%fNPT4s`TSG{))@+!VIE1%6}Rh{Fgdk{-5fK z>>&Q84fMizjVsiSq9H>X@C9SIlzj-K(S=UtS32w<@?;m6i+X!Qb>^o~PxNeO0;;HH ziU~OXBj0c;DOy`w$Dq0paJU9=A$=6IQar~hEe;H~jL`507kf^0OA$;6B_O_}%bo&| z&;|59d=zzl4EePC!00v>1Je{#-26*b;nKV__wlav!wC?y%K&(u^8$D?8h|_$U(LjY zVY2||re+$jz;`Wy?F9|)xmFsc?fW{5{WNfUbd5V`10fyw%e1a~_TGbpbby1T_bm`X zn|TS<5c+-oAC(oK+OcLK2c(anbf95+`0N{03?TbU!>dS-p^bkE0{Ps>l_kBGYul>_ zlAla9`t6j}gNNZ0fUiEIZ50_di4VA%#S!`Ps=KoE=7<%k*Bw0L8Zg}~yz{57bQpAq zaRyHC@10?183AO+Pt^JO9gv1m*f=<=`N5)r$Lrp5>rSX!z-oD<%xlUUeA@ldgExv? z)H1iUA3Z8-uNStG=T_dHBm^@iyEglgDPYHlZ{6RX+G;kCX@qv>r#@g8$^{l0zdB#y zj7atSp(^`+n<5nT{55be`A48i7m*J zB+y30qUQ6M8a#1thfM(*CT{d=qQEw6REotdc~cPSv&f({{>9Q9-d zTKN7`Pr{T;iK@)y@N-u#0NZioz2_%EC7wT#$!u0f6yOp;H0LgL#2Rq=l3yP{EgSIBc$eT6DBGQPURr_y?v!GExr9aVUWlY5N;x%EDKbilg zx-t1XXI$Y0)by_L zGCQ3|V|e?sVLYxscYbyebpsgIB0y0-OGO=jn4PUftzGu)*07(Jf+k`;loWlvI!8p4f#2jZ0Pd4_Vk|Q=Aajo#=fBFJ^sHG-R3a)0 ztx-^we(&t)$D<#j-F7PxjCTU0dEM z-edN4_^$vbT{Zv&jAX!h&&2_5fqxISMGh;LPoU=jb;UM?PM~?)5(l+!Gl$9rd;b92daxrrd|&%t8@MTb>&l%ca?+p<6a=_^(F=6n#Bl56 zw1U!L>2sCSxqqM%bQ=Y@&^+S){zrm@mbzN1qzJ9+!eoPwC+PVpVH0&SoPv`cgYB>V z&wJuc+^2t?8@-OSRgq+sExy#cVb?YbN}B5>~rAp^UvG>N+uY!;P zkN216;B#kV-&Za7^81_lsx|J6tP?RxlfSkJt`s`Wpy{EA{QN~{ zoXebxZ3xH9zXdRBsPgKIhW%+6altvHBxhRPYPVsT=jiQc!CRVams59xZlxW}9auYa zGUo!Eaaa)q?QU>S)a2*i2KnUb!QNIJC_W3zok376Xwm=k>`=!o*x;?8C+(v7-q>&k zAUeZPjT>lH3E%VAoKiOHHV1s?G8(Li#va2XoA{q~#Gk8D>r|-<{Zp@udNjCInV*TQ zc(thz(X=7V9ybqsVw^kYkjy>i>srggmFAuj4w#R-1#Z^%oTWLFNKaq2fWZ4|d;MoY ztdoJ6lSebloRKn?H3H|N8AG;{%PZf+^gyHv0JzZ*gHt6t2bichfB_Vd!3Q(I)mhxX zS>^^z!7&Jt{2dufUDS>B6qopo2QRZ0E_W#*qC}-!Sk#S$OKPgih>t-tRp44YA{wG* zB~C4jlvx$U9(G(LfY?xiWhgJByH@(>3NY|PfQbkwrdBAvZA}HZ)+M(Fyu+)Ja~XVY z>TGqiRm>@L7KxQ4cpq6BI@oyb9`ttX+7 z;o2x#7cH|terkQbCk`w-{^Gj1B)g$PZ5XB7oN{n0v_RQO~g89y&rL(M3 z1AfV&mdS)Ei7BQnWWzUwDgRwenRO4p*VIkXVBpPug`tqQ1fTx;LXzY^aSG=Vz@*p+ z(A%t1{>wWyKR-%3@-S10oFBPGedX?T#-|dXeS2M?>h#5%t)MsMU}3PpKv-C4yfWew zijYV#t#eY)dPG z4}HAtmP=atA=cbKM{GN}F=dLbRPo&1Dkx(`epU8~PZs_WM|vNA2l8cOem8)M(hA^1 zH&8~;_25+#_`5{9>^>ep4X`jYkMJZdn#;CPuMYOi7hMfa_Cxe;5@?7bi{KB?7PHo$3||c)e%biR({6~ zIbfiGse>;uZ94w>hYD9} zptaNtZ&OC&x<8F5JLua}oc*wbSKIy+75nOkkIdnwa-QPp&m62aRvgEV7kVIFaszp% zct3B_3?Qt}GNLkY9`2QAKc}*;8J;EdN$wX#9Gpap7s=1ne!fLL0<2j-zV3^b|bvW(*Q?QYEhe}$6ls&Amn0?P7>zLHYK@Kj`en`)*?E+6*2ewXEl&PD93eFf*^ zZN1~n@jG;G%LgHm@wB6S8sh8NveALjvOsc%-I!aHuSb+V@qI|+RZTcCWsBeGMfdg~ z#_~xI3pP5S`Mj|jg1K`ljj>C%SJcLJl3%fLIIuqaDua5DbGwJJM(^$P3zABF7Qc}<7z*hD@HUltQ+ zOk(5XNmfYy@LrHG+g9dlzmlULKBx|0R(l*98z^OC`Ax}U2lHU%kq(=_br8F{DWIFo zPceKzrxv(0{Ftc|V9(Ttpo0v|cRCD!^g7fcjpeh5uY5{G@54m-2>eDL*j9t1nZ@`a zsAw|MNPlgEfPrJaH$aCYG++)~_|pxO{swo$&$`_8Iy?nKE=K0>xIf(G89May_U!3% z)qntkEk;9h&HvBY%WM1%%jJ-y=~zYhkJMAIDaDQgqeP@E7zt}pKoHmdY54i;e?V+> z^ymPsRKPihd}!EX%JU?KfNu15sMY5&n+G7NP@s8vX8I*tEObVMSCc;}66#U^{3sJT zIKXSfvg5}e3eFCwd&D1_!FPkVTiDtGTgd{JZ{eVuReSUgnM)zZS!v*r?~;Q9Uc({Z)pb<9Ix+kEu2ul)Sy+5#C(eM33dN-Y_ifo z8>lktf%I4fpf{`!h=G^s%??cmq@hT6^&PjF*8S0s;i{o>nZo~&sj#3!GX}IJ2|wFVC!< znY~o$6qvIAUBSErl040M#Tceb&red{4vbB9spmEWnbzNq50FMB!}x6?7E<5)faK8V z0y3^3$i|leY=TbZ7?!OC(q{%ZT*VDILK7N(46w@p{sYj!@h|y+S1b!KsC`wA!;MAX zfgC>r!9hTNtb0vpgAIy+sK`&v_C{zBc58u-ZUqFQ7&=bHRu-QLIskNg*zW_L&1hV8 z5)B9*ND`M&0S%P$!y_NG`%!Emjr#J8o0t&F;Q=t<1eKDOp`$Rl17&%C1E1wG)EsK@ zA7#IXOjUGD7{q`5$EsHmaeE_bt!mgT3ThveP=i+FxpNljeMaCA7qjaB8;LdK>y+>J zw1rx%Re{lz9@5uYYUp9GCs=kc*LeV?PvmA~SjKpDCJ4ku2{emCM`=91WCWIuoicPu z`^kFXPje_M_$Tq(1vGn^NE0E31q91kFR&20Iu)WAX1yR^fUrGB0wAIw9qZq+R+Yd_ zAo-W!=tWIwI)Ij9E^3H>YJbWu163^x6LlUeRaXHy7I;YET~lV<0Ug(dR^#9}9q4o) z5!Qr7)_49_-!+nkjNO1dI=ABz2UG`>7j)o5>=!`O9H7sxLEe`*hd&QxKF z%i}im1x((0Ik-1;R*;L>8Qo*5 zht?(-(M5HM%0neum310G3mip`bx@865C8USh-I=2WSlR?9u z{%AE|Oa^t9xT*d~E}s9h*e=4l1@;#T#7pShBxtV>dlYbU{e$iSEE&W}0{`bf?{{Zs zh#8%q?LhU3|6!aGB|PSS?^4P5Zx~t`KNhq2_M(b)4=5-l7p`SZof4+`!>ZZ*2OHG! z?0M*zi$fL>@Sz=Xy0jLiO$>II2cYF+_ZN_vUMB$A?6H+Gu(0F@FjN9QteyQFZkzs& zNx({#pbvhc7FcGm`|CZ2Cn|v|S%(yC2+nFU2cl| zLsESUEnFIMlWceUG+$I=_4^Bw$evdx#5{)@1%ZP0Z)_~;F{m_8v`WDXDiv^e%X;)V z*Cvv_*Bga@QlR_6aa4^dU=WzbLE%LQbe#g4D-H!rx-9q~Tw8`n1067-TDkFK4mh~#goLyD4E_JgXK)zKf zwtWo|Y5#)o>VJDd$KGIsCPNah5k-NDW{&-7rgK9>n{)**xUr=U721l{)uUCe(zV>y=rweVy z`Nmt@P|IvWbG#9Ih?L8&UB;JL#+(`C;mx*AgT-Ef0RtGzU8%DjlWI>c0{pc2bGpeV zlGYjN6-b`L6qg4tXAQ4g^{~{uVmrrretwD^z<$)}0Xo2A6}1&1FuB=aSRJT9+BI5d zY=fJpragbbMkB-9yjr33;#BDyl4R|E>oIKOcnI4%17DsQJb>KbHv2pz#=FS4)hm}a zLY1%=KXvL$0MG!nxT;4At$+qUH3~jNO{-?(eX_F;%Cltm~X2y~&A!E(@;6_Fu#nX-f2+$q=xfK9e1^?Pm-}_0Ox_?;Kpn3L z_5pU~8oLYUs!kUaZWM7H%vrPZQT$xQ!anU@PT`x8tUzbp<2D}>#k4zUmbQmy*K$nE zinf+vT7oO}m``OaEQU|J7{egXzS!s?KGnl^vH?5s0cW=DoKy&({GySuvsz=$&|@pp zsfkGUY8gk_^mgj9!QAS;I%ZGfP#Dc3P5r{+SoV(q(;(T5Nmy>cuagdKyK1tn-P6`n z(W${SODAtp-hdflccIDvwt$L{btl8(M-6{n(1W9x30imIuOrh#lT{+J3~}DUaQezc z9l=Xp3*xr|N_cvel@wEMq`l-&{9@Vr-PK1ZhE-=`^wfa@i zLbgL8#AYJrT5;y5^(Gj0r09B8GJF3lmxEJSTvfHU@=j`HPgY@y*ax|1+Z3b&uccXFC8Su=@Kfa zDC#Gz%`{DU-d}jyoaPq#q(D@4g~b}EEdO^@{NR?q>{4$DN9X?Jt1%Py|D-kjU32eW ztrSA|iD&WC{}pn_<1QSn{1drLKsw$hF3N@-EKvQOkdW8;t~!f_To_OQfCO0MydK(` zEYN;QW%t$!jUN$Q{W~}Y)AvGnv@!)d@lIbijD4xE3*~VZ=D90~I2ubi+H}(!FvW3L zZ{#?i2BlHL#6LEeFQ6C^`uxY=+lt2aoN0y}XJiU`%Q7z@P2?C3xb;ir&p|qLL4)(x zfp^2p3doukR^>L)A@b5Tllf=0gLh^>i9B?*!PuhnkEDmoZ~cyh`2B>5{+7r+-#}Jq zKXSe7_e%foJ`4`%BUg}?Ek-}rv<_BY6U>39BG)z7zZkpx^A_uII;ch&G)~*roZ+Y)#5NJ)E;Lwi>n}p4~)?gibDtzL%E;s6KpILm4-;yj0L6o zzX78oAZNXhoD`_4G?d(F4os;7r)aW@rZUee?1D0fJ)Z+tR|}Nqe}SV@wy!$ww?5$P zi&K^H$L`uF4yB}~2QKlcC5sE!T2l;6U4A=Ijl@zWeK0);63o&}FM*h@>U99ZU#x$7 zdccjLwomcQ%#ZIvKhcqc&|AB*D3K_KGnJR`F_OW4Kw1azV+&_x?cjI!!X)LI6LG^M zwg^hXyly)PRiYw@X@|GBgZ}RN`Pb8sKYDKaFYjSQ1iCi_VpORDw@gMFwnHG_JzJz)nyCKI{^V{)!sA@(*OBZ6+UyBXgWa&o76B zUX$Tl5!#4aVPNV$DMQtMe^%%9&5!IxbtM1ixSA_e^TZfFpf8f_WVLf2qBigR+$ zTRKXZ^yA$osG~HYFYMCg;HGgJb&FasYWeHpd=vsmmmrD3?Jy z;2}5Cs!d?-a#Zn*(uR9zxhz1=|C8&;7dVz;`N&433lfLms}DT_dvQz|wpo4oraxHn z%g6D1Cd?RcDU6mF8sax2GzO*L0LQk3gZ{x`ya%EKkr-$kB_S#KlN&8s-=w@x)fJ$p zRPMd$X_L3n6!lt?M2gKJ+O~grQ&@*PZ?gjy=l@MEOS9?!5~l{=+W%Ot^S_6r`@i=N zOYutbNBd&<`MGD4YH5h%ce0+sA+E?FxPa4^J@( zywiMdE#J0Zuoqo|n8il?8QPLA+{=7tsNmUY`W?Rl2V+D9q@ry_!$!~p55t^$)6UZ& z4{K{{q5LVOUYRQqYnURPtl;%H;d!HkJ%Pm%%qzX`n!mM>ulhBzA(-Z^k4?K6k&}#3 zRYEhrZWs|90NvOXN-eLiba9BZ)B_tYu4A47c;V_Y2%b7cSw1`(luEOZ4UF>sIa;-Ub-wrCnZ znf7%wf&Dx%bYIFEOz<61xE+Jo*H{;_O(8ui?H_c2qf762ixO|hwKlO7sPd7ol;BKik^wAT*v>Z0d;*enC98C|TlG7qu&6IVMRH7LE)S zl+Tu_=hE63ZO=~z6&k;}@k#JlsfQBD1Wk16A`gc`q4Xv;W4T$m(fCl))NHo>i4z7D z_f9WP`ukJsqbf%&>lX5wAJ*`551a^!s0aT^|JfY#lcwYa1#z(-Q(8d`mr%R9>4Uzf zhU$&#Usz{v?I*29{g@zmsaQNGEBM1ZW3a#j-ut@>aV)ZwS5a@iG_j_YT!Gc@|2nJ&6?eA zaZoMgwpGKk3U>n4(2zZ`RGI&azT)4<1X1*3B}G0$u>f6Pnx5Azhjhj{@g5j<-ITwJ-@~}iL+e=WDf3ZD zP@ZJ$Ia$yf4E79hK+dN!p z=J4@07fg*LpQ~>?Iz%O*o(OzD)eEy65kS=Spm%AV&AsgkuBARy_R06S?nEPBD5`@| zbqW=W%D|gtbEINKQ-BULZS38*G0q>f0Z}mjTph^0y+y^_bHOtQc>b!P0?u}1LK`Gs za`jocM%N03?T%GeBjXiPb`j@$=>~-4WjoaDjywTn+Fh@?+yajZgNd-!#-h&Zi+T%VB0OmvHdNoWp}w{-t4#qK87QJKRcF2_mo_; z%Q5i$dR!zBrudNs0fna_j@-TvWySN}QHFrtFTP))LkNUax&d8LvW1)p!XE+6!y_gh z1^f-jlbdC|SenzT2~eNUeZM~m{f{MahIttFA_H^*H0^D4kqr|fy%7Z)L0l!c48}_4 zf9)Lq4Lk^!Y`=lG5IM*^xptAO`MihkIzRzj1?;r|v?Vu4LZu!@?ypMDdj}OB6G@)n z=}-W^nNLJl{w!te2iP_7d`D->1d1PVO!p)Y=z^i*9#?cGx+AaPQ+y;VrKkGDTU=ZO zV$>0)!H=!cmRlq`4~jVt!Pb!T@_0kfKlXgj@<(vF)@tZ+ve8q+%xDszDiuy(N;4A@W5S z(1A?WcJ^cU3d`*Z zMU|T`0xGKxEg`s!jEtTw4pT&VSlZ5Pd>z#ixmmetz`|l`GPuJ_*A8+?Y{3IFtymD$ zvELRB6>xLapevC5QW8no)?IdeyzSI|^9^>;%M7VMNt9`1rGEkTWfPjr7dw&{R2)L& zrCBjZIcw#-cOr#>HvGBAH+l8)XzF-g?s%-0&!uLMBkDdLd2n6t|eHB$K=eh zvRV;+J+miwRvtQNJ7?o4f3V0$?`t~Wxkf?u{+~(7+TB*&-N3M1q5bSgz}>1v8LVmG+pCk6A3P&($g%3~!qqv7G zH8HO>1jhB%! zM^~)`POvdYnm8;EhoZ0@W7=rj*CYpxBfz)Hu)<@yJ4RxZHuMb{K4L9*k_X4V36&4^ zCKg};X5Hk%^}os$xfJ=x^CnkmNvIjqJ6_GW?IZW%u4^D{Dz9YTZu_Y>qnyF1WhM z&V{vnKiXYs3@5%Q?=kspoZVV?WSTHMM-|U)JZSi?a@=*JAeMGjs|crA^1#rAarr_2o!n%B^Xq1sI3I}#>=tg(Ri{{v2CVFe z4a;mA3s%WJ8n2f2dZZ5KSCRB_QK}P7qJC#=8)fJ*>+V;!$Q7dDq<*Wyyu&P@Geo7K zgY(BUUCt1*Xb?!CoY5S9vl*Yan^g$a$6~EXLUYTF8sdL($wlBhGI2SjI z`z!?nyuqyV?H-(uJk2)=hL+Qc)DJicX@pf2v(A!gy{J(Xv^UOah^RXW(SN)7gf!v; zhG5`SXlEy|VEt4h#r{{=!O@3U_Vz?Ne^^<7J9j zc#f+#kwZ7Fr_6&_UiOQEadOrDf|lquSnPfcW?BZFo$>~P6~<<9CXaTpc#YVr$$m$s zeCgKY-X>0mQxDN)K>vm;_6nF?#zt3J&&+#e-cAJxm7LiRrw zH=XX!KI1P$`iEt)Zx445Rl#TpZWC{9@x96IHGp+I%}cX$9azRCydiG3Yv5)hg!|*Z zM_T31z~V$5r|C<+y()3u$8Ue>>K5Q;6hoaf3Y%LlvDC7h#3h}Lbea#s21aHrI0C~W zqG`BEBon!sT%rjl^NCL%(Lv<}io5tr-_4CmdFrgW2__5-JPSDK~(b>ERQN?JsL$xWLmz0It-Z0p$0O&D}4OTOR9 zsKdJ_@MjJ`Xv!xHOao0TJ(jvt{kE=P;G+iTb)HSr>W5{CPa}Ko9-bePUNy;@qvH1@ zSC@fZE7nJu_%zODM`ELKsfX?=TV1p{TY*C(9){kK?F_1ST^-m%C~R<1z668;QoKo( zf!NkXQ@WFND|okzCcsM9qL-Rn=U4y-9~FAS)NphJh=3$Ia9`}K{+9k}ibJ7K7&zn& z>z^KY!G;R^s!~oAT``a16(MCroB*6zJ`SOV854|~Hj?B{2+}T@S=52rIE4?M69$L~ z6y!8NZlaGD!*up)YI^YuI|WsW5O2uRf!vLE!(Sb_gw9L&J!NAPKXF+LTN3Gd&*(&| z5tb-lNW^~H68a5@+_Zr39%qWT+;_)|QIU#RPmLU2bmEDp;XOYG?`FmYi;qG- zX7k1ij!vkZ>F$|WbA8KzTx){*fE>A9UqzJMj=OR{SSRZFp@|*))O|FY*9#@fsasQn zTd$u8_m-u_lC9-p5F^BAxdo*7qC647Um&p0ELqahdL|xcWe@oT`9j?gRZk^_QSt)> zpnxmLrkSfyoU5AUUOGmzg%f4Mk`?SEzJ*tPY0B&=pag9e{K2teQiO)`Et_@Mt%Zsq zEa(a~d5a!TQ0t`n1V7O(L0{O>(klDTWTU~lmfN^i!>FB#UaDJFyyqkmR8tsa>YT=X z&&T7U;$Blnku`{+DA<#ad=q1?nsVMf1w7GXoi~o=1s}(|P_`dS7Entu`t{Bo)QH$i z623nRyu5${@ryS4@_+%>c?aCY^u{QkXJNK&4Nyw*%w4T0Z3*W9X#!8U)%hcE@XW><{XUwzNi0#6+j?~fYQ)^uaA#gN zNK{)RGZ~EPXT%X82~^zHW#X*D2w{$$nMX_heBbLfK)vIa zn#huC-v0`4H}&s(lem_zOtJ`7p1BsQhbMfBB#!7Ov7&qDtiT~!&_!T+-6RRg4dh*R z^cIpC$aUCDsF{Y78a|r`TGJcA^%GuT^zP%@4e{F(?+~oNSw>llVw+xkf+ph6-`N~b zLKuuMfk=OH;P|>05o8PVB|IJmo$$Ynt0u_-GKumSEFw*A1)vkPu6LtauWq)>XK{Yw zgYpvng;C7`38?~NP`1Lnd<30MVAU8~RHe+PIIIe1ij@Q)EC`f@+9qvGJYqk64Lm>Z z6baDNKDxUH%DbtyyyBeNTx7t)8XEIM_jTSv2=guzz4EZ}jkl|^%6Szkz>wMprU5V0 zJTP8OjYgIT1B?@B4+lsgX#A%9MoiLaXwc1gfYb?0O{Ysl;jDi_XPPn5p;V|#L-zZn zyo4L2%lb)T7~Mt*n0lnJXM>hwqWw^O1?uDaOFu$|<>55gCf%-@^zx;YOEbpTeUB;z z(`EzMsRNQe;2!??_+1os1eD~nJP!q7dRQjJ;JDk1V{{kSU&g@Cy5ED>m&}0w=HdQ^lM5=86jlmv Su7oT4-|8o|Po%4v-T4oXNjz`> diff --git a/test/test_gradient/save_gempy_data.py b/test/test_gradient/save_gempy_data.py deleted file mode 100644 index 78f5dbb5..00000000 --- a/test/test_gradient/save_gempy_data.py +++ /dev/null @@ -1,70 +0,0 @@ -import os,sys -from datetime import datetime -import numpy as np - -import gempy as gp -import gempy_engine -import gempy_viewer as gpv - - - - -from helpers import * -from generate_samples import * - - -import warnings -warnings.filterwarnings("ignore") - - -def main(): - - - # ---------------- 1️⃣ Create the Mesh ---------------- - - - n_points = 32 # number of points along each axis - - x = np.linspace(0, 1, n_points) - z = np.linspace(0, 1, n_points) - - X, Z = np.meshgrid(x, z) # creates grid of shape (n_points, n_points) - - mesh = np.stack([X.ravel(), Z.ravel()], axis=1) # shape: (n_points^2, 2) - - # ---------------- 2️⃣ Generate samples for the input paramter of the gempy, output and it's jacobian ---------------- - data, total_time_input, total_time_output = generate_input_output_gempy_data(mesh=mesh, number_samples=50) - - - filename = "./data_9_parameter.json" - c ,m_data, dmdc_data = np.array(data["input"]),np.array(data["Gempy_output"]), np.array(data["Jacobian_Gempy"]) - print("Shapes-" , "Gempy Input: ", c.shape, "Gempy Output:", m_data.shape, "Jacobian shape:", dmdc_data.shape, "\n") - - print("Time required to generate samples for the input:\n", total_time_input) - print("Time required to generate samples for the output of gempy and it's Jacobian matrix:\n", total_time_output) - - # ---------------- 3️⃣ Save the samples in a file ---------------- - # with open(filename, 'w') as file: - # json.dump(data, file) - - - -if __name__ == "__main__": - - # Your main script code starts here - print("Script started...") - - # Record the start time - start_time = datetime.now() - - main() - # Record the end time - end_time = datetime.now() - - # Your main script code ends here - print("Script ended...") - - # Calculate the elapsed time - elapsed_time = end_time - start_time - - print(f"Elapsed time: {elapsed_time}") From 4fc3bf1004cd842a915503a75ca5cb36ca3a7f97 Mon Sep 17 00:00:00 2001 From: javoha Date: Thu, 12 Jun 2025 10:11:38 +0200 Subject: [PATCH 05/12] Update test including timing for compute_at. --- .../test_compute_at_computation_time.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/test/test_modules/test_compute_at_computation_time.py b/test/test_modules/test_compute_at_computation_time.py index 10909939..790a2c88 100644 --- a/test/test_modules/test_compute_at_computation_time.py +++ b/test/test_modules/test_compute_at_computation_time.py @@ -56,5 +56,19 @@ def test_compute_at_computation_time(): end_time = time.perf_counter() computation_time_topo = end_time - start_time - print(f"Computation time without topography: {computation_time_model:.2f} seconds") - print(f"Computation time with topography: {computation_time_topo:.2f} seconds") + # numpy array with random coordinates within the extent of the model + custom_coordinates = np.random.uniform( + low=geo_model.grid.extent[:3], + high=geo_model.grid.extent[3:], + size=(1000, 3) + ) + + start_time = time.perf_counter() + gp.compute_model_at(geo_model, custom_coordinates) + end_time = time.perf_counter() + computation_time_at = end_time - start_time + + print(f"Computation only model dense grid 125*50*50: {computation_time_model:.2f} seconds") + print(f"Computation time with topography 125*50: {computation_time_topo:.2f} seconds") + print(f"Computation compute_at with 1000 custom points: {computation_time_at:.2f} seconds") + From 341a9ae24e2d00d91044aadcb050fdb6ec7af9f3 Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Thu, 12 Jun 2025 10:43:37 +0200 Subject: [PATCH 06/12] [ENH] Making sure that compute_model_at only computes the custom grid --- gempy/API/compute_API.py | 7 +++++-- gempy/API/grid_API.py | 4 ++-- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/gempy/API/compute_API.py b/gempy/API/compute_API.py index acedbf60..25c1fcdd 100644 --- a/gempy/API/compute_API.py +++ b/gempy/API/compute_API.py @@ -82,11 +82,14 @@ def compute_model_at(gempy_model: GeoModel, at: np.ndarray, Returns: np.ndarray: The computed geological model at the specified coordinates. """ + + print("WARNING: This function sets a custom grid and computes the model so be wary of side effects.") set_custom_grid( grid=gempy_model.grid, - xyz_coord=at + xyz_coord=at, + reset=True ) - + sol = compute_model(gempy_model, engine_config, validate_serialization=True) return sol.raw_arrays.custom diff --git a/gempy/API/grid_API.py b/gempy/API/grid_API.py index 53ec5680..22ded97d 100644 --- a/gempy/API/grid_API.py +++ b/gempy/API/grid_API.py @@ -88,11 +88,11 @@ def set_topography_from_file(grid: Grid, filepath: str, crop_to_extent: Union[Se return set_topography_from_subsurface_structured_grid(grid, struct) -def set_custom_grid(grid: Grid, xyz_coord: np.ndarray): +def set_custom_grid(grid: Grid, xyz_coord: np.ndarray, reset: bool = False): custom_grid = CustomGrid(values=xyz_coord) grid.custom_grid = custom_grid - set_active_grid(grid, [Grid.GridTypes.CUSTOM]) + set_active_grid(grid, grid_type=[Grid.GridTypes.CUSTOM], reset=reset) return grid.custom_grid From 43773e744d654fc2303d5a556d77e6e42e01799d Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Thu, 12 Jun 2025 10:53:58 +0200 Subject: [PATCH 07/12] [TEST] Breaking down tests --- ...ime.py => test_compute_times_for_grids.py} | 92 +++++++++++-------- 1 file changed, 55 insertions(+), 37 deletions(-) rename test/test_modules/{test_compute_at_computation_time.py => test_compute_times_for_grids.py} (64%) diff --git a/test/test_modules/test_compute_at_computation_time.py b/test/test_modules/test_compute_times_for_grids.py similarity index 64% rename from test/test_modules/test_compute_at_computation_time.py rename to test/test_modules/test_compute_times_for_grids.py index 790a2c88..1c32f6cd 100644 --- a/test/test_modules/test_compute_at_computation_time.py +++ b/test/test_modules/test_compute_times_for_grids.py @@ -6,12 +6,61 @@ PLOT = True -def test_compute_at_computation_time(): +def test_compute_time_topo_dense_grid(): + geo_model: gp.data.GeoModel = _setup_model() + # Compute a solution for the model + geo_model.grid.active_grids = gp.data.Grid.GridTypes.TOPOGRAPHY + start_time = time.perf_counter() + gp.compute_model(geo_model) + print(f"Computing model on grids: {geo_model.grid.active_grids}") + end_time = time.perf_counter() + computation_time_topo = end_time - start_time + + + geo_model.grid.active_grids = gp.data.Grid.GridTypes.DENSE + start_time = time.perf_counter() + gp.compute_model(geo_model) + print(f"Computing model on grids: {geo_model.grid.active_grids}") + end_time = time.perf_counter() + computation_time_dense = end_time - start_time + + # Recompute model as a new grid was added + geo_model.grid.active_grids = gp.data.Grid.GridTypes.TOPOGRAPHY | gp.data.Grid.GridTypes.DENSE + start_time = time.perf_counter() + gp.compute_model(geo_model) + print(f"Computing model on grids: {geo_model.grid.active_grids}") + end_time = time.perf_counter() + computation_time_topo_dense = end_time - start_time + + + print(f"Computation only model dense grid 125*50*50: {computation_time_dense:.2f} seconds") + print(f"Computation time with topography 125*50: {computation_time_topo:.2f} seconds") + print(f"Computation time with topography and dense grid 125*50*50: {computation_time_topo_dense:.2f} seconds") + + +def test_compute_time_custom_dense_grid(): + geo_model: gp.data.GeoModel = _setup_model() + + # numpy array with random coordinates within the extent of the model + custom_coordinates = np.random.uniform( + low=geo_model.grid.extent[:3], + high=geo_model.grid.extent[3:], + size=(1000, 3) + ) + + start_time = time.perf_counter() + gp.compute_model_at(geo_model, custom_coordinates) + end_time = time.perf_counter() + computation_time_at = end_time - start_time + + print(f"Computation compute_at with 1000 custom points: {computation_time_at:.2f} seconds") + + +def _setup_model(): # Define the path to data data_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/' path_to_data = data_path + "/data/input_data/jan_models/" - # Create a GeoModel instance geo_model = gp.create_geomodel( project_name='EGU_example', @@ -22,26 +71,17 @@ def test_compute_at_computation_time(): path_to_surface_points=path_to_data + "model7_surface_points.csv" ) ) - # Map geological series to surfaces gp.map_stack_to_surfaces( gempy_model=geo_model, mapping_object={ - "Fault_Series": ('fault'), - "Strat_Series1": ('rock3'), - "Strat_Series2": ('rock2', 'rock1'), + "Fault_Series" : ('fault'), + "Strat_Series1": ('rock3'), + "Strat_Series2": ('rock2', 'rock1'), } ) - # Define youngest structural group as fault gp.set_is_fault(geo_model, ["Fault_Series"]) - - # Compute a solution for the model - start_time = time.perf_counter() - gp.compute_model(geo_model) - end_time = time.perf_counter() - computation_time_model = end_time - start_time - # Setting a randomly generated topography gp.set_topography_from_random( grid=geo_model.grid, @@ -49,26 +89,4 @@ def test_compute_at_computation_time(): d_z=np.array([700, 950]), topography_resolution=np.array([125, 50]) ) - - # Recompute model as a new grid was added - start_time = time.perf_counter() - gp.compute_model(geo_model) - end_time = time.perf_counter() - computation_time_topo = end_time - start_time - - # numpy array with random coordinates within the extent of the model - custom_coordinates = np.random.uniform( - low=geo_model.grid.extent[:3], - high=geo_model.grid.extent[3:], - size=(1000, 3) - ) - - start_time = time.perf_counter() - gp.compute_model_at(geo_model, custom_coordinates) - end_time = time.perf_counter() - computation_time_at = end_time - start_time - - print(f"Computation only model dense grid 125*50*50: {computation_time_model:.2f} seconds") - print(f"Computation time with topography 125*50: {computation_time_topo:.2f} seconds") - print(f"Computation compute_at with 1000 custom points: {computation_time_at:.2f} seconds") - + return geo_model From e56a388cc8ac01951725e6adfaffd97bef947441 Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Thu, 12 Jun 2025 12:01:24 +0200 Subject: [PATCH 08/12] [BUG/CLN] Making sure that computing only custom grids works fine --- gempy/core/data/geo_model.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/gempy/core/data/geo_model.py b/gempy/core/data/geo_model.py index b04740db..1e913f86 100644 --- a/gempy/core/data/geo_model.py +++ b/gempy/core/data/geo_model.py @@ -101,14 +101,6 @@ def solutions(self, value: Solutions): self._solutions = value - # * Set solutions per group - if self._solutions.raw_arrays is not None: - for e, group in enumerate(self.structural_frame.structural_groups): - group.kriging_solution = RawArraysSolution( # ? Maybe I need to add more fields, but I am not sure yet - scalar_field_matrix=self._solutions.raw_arrays.scalar_field_matrix[e], - block_matrix=self._solutions.raw_arrays.block_matrix[e], - ) - # * Set solutions per element for e, element in enumerate(self.structural_frame.structural_elements[:-1]): # * Ignore basement element.scalar_field_at_interface = value.scalar_field_at_surface_points[e] From f40ee3f1376efdd08d84c177158fb21c9c600c62 Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Thu, 12 Jun 2025 13:07:08 +0200 Subject: [PATCH 09/12] [WIP] Investigating time explosion when multiple recomputations --- .../test_compute_times_for_grids.py | 121 +++++++++++++++++- 1 file changed, 117 insertions(+), 4 deletions(-) diff --git a/test/test_modules/test_compute_times_for_grids.py b/test/test_modules/test_compute_times_for_grids.py index 1c32f6cd..a4e60f84 100644 --- a/test/test_modules/test_compute_times_for_grids.py +++ b/test/test_modules/test_compute_times_for_grids.py @@ -5,25 +5,54 @@ PLOT = True +def test_compute_time_dense_dense(): -def test_compute_time_topo_dense_grid(): geo_model: gp.data.GeoModel = _setup_model() - # Compute a solution for the model - geo_model.grid.active_grids = gp.data.Grid.GridTypes.TOPOGRAPHY + geo_model.grid.active_grids = gp.data.Grid.GridTypes.DENSE start_time = time.perf_counter() gp.compute_model(geo_model) print(f"Computing model on grids: {geo_model.grid.active_grids}") end_time = time.perf_counter() - computation_time_topo = end_time - start_time + computation_time_dense_I = end_time - start_time + + + start_time = time.perf_counter() + gp.compute_model(geo_model) + print(f"Computing model on grids: {geo_model.grid.active_grids}") + end_time = time.perf_counter() + computation_time_dense_II = end_time - start_time + + + start_time = time.perf_counter() + gp.compute_model(geo_model) + print(f"Computing model on grids: {geo_model.grid.active_grids}") + end_time = time.perf_counter() + computation_time_dense_III = end_time - start_time + + print(f"Computation only model dense grid 125*50*50: {computation_time_dense_I:.2f} seconds") + print(f"Computation only model dense grid 125*50*50: {computation_time_dense_II:.2f} seconds") + print(f"Computation only model dense grid 125*50*50: {computation_time_dense_III:.2f} seconds") +def test_compute_time_topo_dense_grid(): + geo_model: gp.data.GeoModel = _setup_model() + geo_model.grid.active_grids = gp.data.Grid.GridTypes.DENSE start_time = time.perf_counter() gp.compute_model(geo_model) print(f"Computing model on grids: {geo_model.grid.active_grids}") end_time = time.perf_counter() computation_time_dense = end_time - start_time + + # Compute a solution for the model + geo_model.grid.active_grids = gp.data.Grid.GridTypes.TOPOGRAPHY + start_time = time.perf_counter() + gp.compute_model(geo_model) + print(f"Computing model on grids: {geo_model.grid.active_grids}") + end_time = time.perf_counter() + computation_time_topo = end_time - start_time + # Recompute model as a new grid was added geo_model.grid.active_grids = gp.data.Grid.GridTypes.TOPOGRAPHY | gp.data.Grid.GridTypes.DENSE @@ -37,6 +66,11 @@ def test_compute_time_topo_dense_grid(): print(f"Computation only model dense grid 125*50*50: {computation_time_dense:.2f} seconds") print(f"Computation time with topography 125*50: {computation_time_topo:.2f} seconds") print(f"Computation time with topography and dense grid 125*50*50: {computation_time_topo_dense:.2f} seconds") + """ + Computation only model dense grid 125*50*50: 79.56 seconds + Computation time with topography 125*50: 5.34 seconds + Computation time with topography and dense grid 125*50*50: 59.24 seconds + """ def test_compute_time_custom_dense_grid(): @@ -57,6 +91,85 @@ def test_compute_time_custom_dense_grid(): print(f"Computation compute_at with 1000 custom points: {computation_time_at:.2f} seconds") +def test_compute_at_computation_time(): + + # Define the path to data + data_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/' + path_to_data = data_path + "/data/input_data/jan_models/" + + # Create a GeoModel instance + geo_model = gp.create_geomodel( + project_name='EGU_example', + extent=[0, 2500, 0, 1000, 0, 1000], + resolution=[125, 50, 50], + importer_helper=gp.data.ImporterHelper( + path_to_orientations=path_to_data + "model7_orientations.csv", + path_to_surface_points=path_to_data + "model7_surface_points.csv" + ) + ) + + # Map geological series to surfaces + gp.map_stack_to_surfaces( + gempy_model=geo_model, + mapping_object={ + "Fault_Series": ('fault'), + "Strat_Series1": ('rock3'), + "Strat_Series2": ('rock2', 'rock1'), + } + ) + + # Define youngest structural group as fault + gp.set_is_fault(geo_model, ["Fault_Series"]) + + # Compute a solution for the model + start_time = time.perf_counter() + gp.compute_model(geo_model) + print(f"Computing model on grids: {geo_model.grid.active_grids}") + end_time = time.perf_counter() + computation_time_model = end_time - start_time + + return + + # Setting a randomly generated topography + gp.set_topography_from_random( + grid=geo_model.grid, + fractal_dimension=2, + d_z=np.array([700, 950]), + topography_resolution=np.array([125, 50]) + ) + + # Recompute model as a new grid was added + start_time = time.perf_counter() + gp.compute_model(geo_model) + print(f"Computing model on grids: {geo_model.grid.active_grids}") + end_time = time.perf_counter() + computation_time_topo = end_time - start_time + + # numpy array with random coordinates within the extent of the model + custom_coordinates = np.random.uniform( + low=geo_model.grid.extent[:3], + high=geo_model.grid.extent[3:], + size=(1000, 3) + ) + + start_time = time.perf_counter() + gp.compute_model_at(geo_model, custom_coordinates) + print(f"Computing model on grids: {geo_model.grid.active_grids}") + end_time = time.perf_counter() + computation_time_at = end_time - start_time + + print(f"Computation only model dense grid 125*50*50: {computation_time_model:.2f} seconds") + print(f"Computation time with topography 125*50: {computation_time_topo:.2f} seconds") + print(f"Computation compute_at with 1000 custom points: {computation_time_at:.2f} seconds") + + """ + Computation only model dense grid 125*50*50: 19.96 seconds + Computation time with topography 125*50: 60.75 seconds + Computation compute_at with 1000 custom points: 7.88 seconds + """ + + + def _setup_model(): # Define the path to data data_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/' From 25e78fc002bc17c61797db16911d48b990a92d47 Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Thu, 12 Jun 2025 13:12:22 +0200 Subject: [PATCH 10/12] [WIP] Investigating time explosion when multiple recomputations --- test/test_modules/test_compute_times_for_grids.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/test/test_modules/test_compute_times_for_grids.py b/test/test_modules/test_compute_times_for_grids.py index a4e60f84..359dd42c 100644 --- a/test/test_modules/test_compute_times_for_grids.py +++ b/test/test_modules/test_compute_times_for_grids.py @@ -4,6 +4,9 @@ import time PLOT = True +# DENSE_RESOLUTION = [125, 50, 50] +DENSE_RESOLUTION = [20,20,20] + def test_compute_time_dense_dense(): @@ -178,7 +181,7 @@ def _setup_model(): geo_model = gp.create_geomodel( project_name='EGU_example', extent=[0, 2500, 0, 1000, 0, 1000], - resolution=[125, 50, 50], + resolution=DENSE_RESOLUTION, importer_helper=gp.data.ImporterHelper( path_to_orientations=path_to_data + "model7_orientations.csv", path_to_surface_points=path_to_data + "model7_surface_points.csv" From c093e2dfe546e475c81d4781c7a0e6c6c9c6af59 Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Thu, 12 Jun 2025 14:19:57 +0200 Subject: [PATCH 11/12] [TEST] Adding asserts to tests --- .../test_compute_times_for_grids.py | 111 ++---------------- 1 file changed, 11 insertions(+), 100 deletions(-) diff --git a/test/test_modules/test_compute_times_for_grids.py b/test/test_modules/test_compute_times_for_grids.py index 359dd42c..aa84c3f5 100644 --- a/test/test_modules/test_compute_times_for_grids.py +++ b/test/test_modules/test_compute_times_for_grids.py @@ -4,13 +4,12 @@ import time PLOT = True -# DENSE_RESOLUTION = [125, 50, 50] -DENSE_RESOLUTION = [20,20,20] +DENSE_RESOLUTION = [20, 20, 20] def test_compute_time_dense_dense(): - geo_model: gp.data.GeoModel = _setup_model() + geo_model.interpolation_options.evaluation_options.mesh_extraction = True geo_model.grid.active_grids = gp.data.Grid.GridTypes.DENSE start_time = time.perf_counter() @@ -19,23 +18,17 @@ def test_compute_time_dense_dense(): end_time = time.perf_counter() computation_time_dense_I = end_time - start_time - start_time = time.perf_counter() gp.compute_model(geo_model) print(f"Computing model on grids: {geo_model.grid.active_grids}") end_time = time.perf_counter() computation_time_dense_II = end_time - start_time - - start_time = time.perf_counter() - gp.compute_model(geo_model) - print(f"Computing model on grids: {geo_model.grid.active_grids}") - end_time = time.perf_counter() - computation_time_dense_III = end_time - start_time - print(f"Computation only model dense grid 125*50*50: {computation_time_dense_I:.2f} seconds") print(f"Computation only model dense grid 125*50*50: {computation_time_dense_II:.2f} seconds") - print(f"Computation only model dense grid 125*50*50: {computation_time_dense_III:.2f} seconds") + + # Assert that it is not too different + assert abs(computation_time_dense_I - computation_time_dense_II) < 0.2 def test_compute_time_topo_dense_grid(): @@ -47,7 +40,7 @@ def test_compute_time_topo_dense_grid(): print(f"Computing model on grids: {geo_model.grid.active_grids}") end_time = time.perf_counter() computation_time_dense = end_time - start_time - + # Compute a solution for the model geo_model.grid.active_grids = gp.data.Grid.GridTypes.TOPOGRAPHY start_time = time.perf_counter() @@ -56,7 +49,6 @@ def test_compute_time_topo_dense_grid(): end_time = time.perf_counter() computation_time_topo = end_time - start_time - # Recompute model as a new grid was added geo_model.grid.active_grids = gp.data.Grid.GridTypes.TOPOGRAPHY | gp.data.Grid.GridTypes.DENSE start_time = time.perf_counter() @@ -65,16 +57,14 @@ def test_compute_time_topo_dense_grid(): end_time = time.perf_counter() computation_time_topo_dense = end_time - start_time - print(f"Computation only model dense grid 125*50*50: {computation_time_dense:.2f} seconds") print(f"Computation time with topography 125*50: {computation_time_topo:.2f} seconds") print(f"Computation time with topography and dense grid 125*50*50: {computation_time_topo_dense:.2f} seconds") - """ - Computation only model dense grid 125*50*50: 79.56 seconds - Computation time with topography 125*50: 5.34 seconds - Computation time with topography and dense grid 125*50*50: 59.24 seconds - """ - + + # Assert that dense takes longer than topo and that the sum of both is close + assert computation_time_dense > computation_time_topo + assert computation_time_topo_dense > computation_time_dense + def test_compute_time_custom_dense_grid(): geo_model: gp.data.GeoModel = _setup_model() @@ -94,85 +84,6 @@ def test_compute_time_custom_dense_grid(): print(f"Computation compute_at with 1000 custom points: {computation_time_at:.2f} seconds") -def test_compute_at_computation_time(): - - # Define the path to data - data_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/' - path_to_data = data_path + "/data/input_data/jan_models/" - - # Create a GeoModel instance - geo_model = gp.create_geomodel( - project_name='EGU_example', - extent=[0, 2500, 0, 1000, 0, 1000], - resolution=[125, 50, 50], - importer_helper=gp.data.ImporterHelper( - path_to_orientations=path_to_data + "model7_orientations.csv", - path_to_surface_points=path_to_data + "model7_surface_points.csv" - ) - ) - - # Map geological series to surfaces - gp.map_stack_to_surfaces( - gempy_model=geo_model, - mapping_object={ - "Fault_Series": ('fault'), - "Strat_Series1": ('rock3'), - "Strat_Series2": ('rock2', 'rock1'), - } - ) - - # Define youngest structural group as fault - gp.set_is_fault(geo_model, ["Fault_Series"]) - - # Compute a solution for the model - start_time = time.perf_counter() - gp.compute_model(geo_model) - print(f"Computing model on grids: {geo_model.grid.active_grids}") - end_time = time.perf_counter() - computation_time_model = end_time - start_time - - return - - # Setting a randomly generated topography - gp.set_topography_from_random( - grid=geo_model.grid, - fractal_dimension=2, - d_z=np.array([700, 950]), - topography_resolution=np.array([125, 50]) - ) - - # Recompute model as a new grid was added - start_time = time.perf_counter() - gp.compute_model(geo_model) - print(f"Computing model on grids: {geo_model.grid.active_grids}") - end_time = time.perf_counter() - computation_time_topo = end_time - start_time - - # numpy array with random coordinates within the extent of the model - custom_coordinates = np.random.uniform( - low=geo_model.grid.extent[:3], - high=geo_model.grid.extent[3:], - size=(1000, 3) - ) - - start_time = time.perf_counter() - gp.compute_model_at(geo_model, custom_coordinates) - print(f"Computing model on grids: {geo_model.grid.active_grids}") - end_time = time.perf_counter() - computation_time_at = end_time - start_time - - print(f"Computation only model dense grid 125*50*50: {computation_time_model:.2f} seconds") - print(f"Computation time with topography 125*50: {computation_time_topo:.2f} seconds") - print(f"Computation compute_at with 1000 custom points: {computation_time_at:.2f} seconds") - - """ - Computation only model dense grid 125*50*50: 19.96 seconds - Computation time with topography 125*50: 60.75 seconds - Computation compute_at with 1000 custom points: 7.88 seconds - """ - - - def _setup_model(): # Define the path to data data_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/' From 65743a9e3209156382608d7ae460974e2068efe4 Mon Sep 17 00:00:00 2001 From: MigueldelaVarga Date: Fri, 13 Jun 2025 10:36:54 +0200 Subject: [PATCH 12/12] [TEST] Update gravity test to disable scalar gradient computation Modified the approved test file to set `compute_scalar_gradient` to `false`, ensuring consistency with test expectations and configuration. --- .../test_gravity.test_gravity.verify/2-layers.approved.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_modules/_geophysics_TO_UPDATE/test_gravity.test_gravity.verify/2-layers.approved.txt b/test/test_modules/_geophysics_TO_UPDATE/test_gravity.test_gravity.verify/2-layers.approved.txt index 2dc04b80..acbe62f8 100644 --- a/test/test_modules/_geophysics_TO_UPDATE/test_gravity.test_gravity.verify/2-layers.approved.txt +++ b/test/test_modules/_geophysics_TO_UPDATE/test_gravity.test_gravity.verify/2-layers.approved.txt @@ -162,7 +162,7 @@ "mesh_extraction_masking_options": 3, "mesh_extraction_fancy": true, "evaluation_chunk_size": 500000, - "compute_scalar_gradient": true, + "compute_scalar_gradient": false, "verbose": false }, "debug": true,