From fa53201e99e8c3bca62188f34d8f746217db74d6 Mon Sep 17 00:00:00 2001 From: Alan Coreas Date: Thu, 12 Dec 2024 15:51:23 -0700 Subject: [PATCH 1/7] MLX implementation of Rayleigh Integral --- BabelViscoFDTD/tools/RayleighAndBHTE.py | 177 +++++++++++++++++++++++ BabelViscoFDTD/tools/rayleigh.cpp | 84 +++++++++++ BabelViscoFDTD/tools/rayleighAndBHTE.hpp | 12 ++ 3 files changed, 273 insertions(+) create mode 100644 BabelViscoFDTD/tools/rayleigh.cpp create mode 100644 BabelViscoFDTD/tools/rayleighAndBHTE.hpp diff --git a/BabelViscoFDTD/tools/RayleighAndBHTE.py b/BabelViscoFDTD/tools/RayleighAndBHTE.py index 0070495..e1db321 100644 --- a/BabelViscoFDTD/tools/RayleighAndBHTE.py +++ b/BabelViscoFDTD/tools/RayleighAndBHTE.py @@ -7,14 +7,30 @@ CUDA is automatically selected if running Windows or Linux, while OpenCL is selected in MacOs ''' +import gc import numpy as np import os +from pathlib import Path import sys from sysconfig import get_paths import ctypes import sys import platform +_IS_MAC = platform.system() == 'Darwin' + +def resource_path(): # needed for bundling + """Get absolute path to resource, works for dev and for PyInstaller""" + if not _IS_MAC: + return os.path.split(Path(__file__))[0] + + if getattr(sys, 'frozen', False) and hasattr(sys, '_MEIPASS'): + bundle_dir = os.path.abspath(os.path.join(os.path.dirname(__file__))) + else: + bundle_dir = Path(__file__).parent + + return bundle_dir + KernelCoreSourceBHTE=""" #define Tref 43.0 unsigned int DzDy=outerDimz*outerDimy; @@ -550,6 +566,66 @@ def InitMetal(DeviceName='AMD'): ctx.set_external_gpu(1) prgcl = ctx.kernel('#define _METAL\n'+RayleighOpenCLMetalSource+OpenCLKernelBHTE) +def InitMLX(DeviceName='AMD'): + global ctx + global prgcl_mlx + global clp_mlx + global sel_device + + import mlx.core as mx + clp_mlx = mx + + kernel_files = [ + os.path.join(resource_path(), 'rayleigh.cpp'), + os.path.join(resource_path(), 'BHTE.cpp'), + ] + + with open(os.path.join(resource_path(), 'rayleighAndBHTE.hpp'), 'r') as f: + header = f.read() + + kernel_functions = [{'name': 'ForwardPropagationKernel', + 'file': kernel_files[0], + 'input_names': ["mr2", "c_wvnb_real", "c_wvnb_imag", "MaxDistance", "mr1", + "r2pr","r1pr","a1pr","u1_real","u1_imag","mr1step"], + 'output_names': ["py_data_u2_real","py_data_u2_imag"], + 'atomic_outputs': False}, + {'name': 'BHTEFDTDKernel', + 'file': kernel_files[1], + 'input_names': ["d_input","d_input2","d_bhArr","d_perfArr","d_labels","d_Qarr", + "d_pointsMonitoring","CoreTemp","sonication","outerDimx","outerDimy", + "outerDimz","dt","TotalStepsMonitoring","nFactorMonitoring","n_Step", + "SelJ","StartIndexQ","TotalSteps"], + 'output_names': ["d_output","d_output2","d_MonitorSlice","d_Temppoints"], + 'atomic_outputs': False}] + + preamble = '#define _MLX\n' + build_later = True + + sel_device = mx.default_device() + print('Selecting device: ', sel_device) + + kernels = {} + for kf in kernel_functions: + with open(kf['file'], 'r') as f: + lines = f.readlines() + kernel_code = ''.join(lines[:-1]) # Remove last bracket + + if build_later: + kf['header'] = preamble+header + kf['source'] = kernel_code + kernels[kf['name']] = kf + else: + kernel = mx.fast.metal_kernel(name = kf['name'], + input_names = kf['input_names'], + output_names = kf['output_names'], + atomic_outputs = kf['atomic_outputs'], + header = preamble + header, + source = kernel_code) + + kernels[kf['name']] = kernel + + prgcl_mlx = kernels + def ForwardSimpleCUDA(cwvnb,center,ds,u0,rf,MaxDistance=-1.0,u0step=0): if u0step!=0: mr1=u0step @@ -656,6 +732,105 @@ def ForwardSimpleOpenCL(cwvnb,center,ds,u0,rf,MaxDistance=-1.0,u0step=0): return u2 +def ForwardSimpleMLX(cwvnb,center,ds,u0,rf,MaxDistance=-1.0,u0step=0): + global queue + global prgcl_mlx + global ctx + global clp_mlx + global sel_device + + mr2=rf.shape[0] + + if u0step!=0: + mr1=u0step + assert(mr1*rf.shape[0]==u0.shape[0]) + assert(mr1==center.shape[0]) + assert(mr1==ds.shape[0]) + else: + mr1=center.shape[0] + + # Change to mlx arrays + d_r1pr = clp_mlx.array(center) + d_u1realpr=clp_mlx.array(np.real(u0)) + d_u1imagpr=clp_mlx.array(np.imag(u0)) + d_a1pr = clp_mlx.array(ds) + + u2 = np.zeros((rf.shape[0]),dtype=np.complex64) + + # Build program from source code + knl = clp_mlx.fast.metal_kernel(name = f"{prgcl_mlx['ForwardPropagationKernel']['name']}", + input_names = prgcl_mlx['ForwardPropagationKernel']['input_names'], + output_names = prgcl_mlx['ForwardPropagationKernel']['output_names'], + source = prgcl_mlx['ForwardPropagationKernel']['source'], + header = prgcl_mlx['ForwardPropagationKernel']['header'], + atomic_outputs = prgcl_mlx['ForwardPropagationKernel']['atomic_outputs'], + ) + + # We need to split in small chunks to be sure the kernel does not take too much time + # otherwise the OS will kill it + + NonBlockingstep = int(24000e6) + step = int(NonBlockingstep/mr1) + + if step > mr2: + step = mr2 + if step < 5: + step = 5 + + slice_start = 0 + slice_end = 0 + + while slice_start < mr2: + + slice_end = min(slice_start + step, mr2) + chunk_size = slice_end-slice_start + + print(f"Working on slices {slice_start} to {slice_end} out of {u2.shape[0]}") + + # Grab section of data + d_r2pr = clp_mlx.array(rf[slice_start:slice_end,:]) + d_u2realpr = clp_mlx.zeros_like(d_r2pr) + d_u2imagpr = clp_mlx.zeros_like(d_r2pr) + + # Deploy kernel + d_u2realpr,d_u2imagpr = knl(inputs = [np.int32(slice_end), + np.float32(np.real(cwvnb)), + np.float32(np.imag(cwvnb)), + np.float32(MaxDistance), + np.int32(mr1), + d_r2pr, + d_r1pr, + d_a1pr, + d_u1realpr, + d_u1imagpr, + np.int32(u0step)], + output_shapes = [(chunk_size,),(chunk_size,)], + output_dtypes = [clp_mlx.float32,clp_mlx.float32], + grid=(chunk_size,1,1), + threadgroup=(256, 1, 1), + verbose=False, + stream=sel_device) + clp_mlx.synchronize() + + # Change back to numpy array + u2_real = np.array(d_u2realpr) + u2_imag = np.array(d_u2imagpr) + + # Combine real & imag parts + u2_section = u2_real+1j*u2_imag + + # Update final array + u2[slice_start:slice_end] = u2_section + + # Clean up mlx arrays + del d_u2realpr,d_u2imagpr + gc.collect() + + # Update starting location + slice_start += step + + return u2 + def ForwardSimpleMetal(cwvnb,center,ds,u0,rf,deviceName,MaxDistance=-1.0,u0step=0): os.environ['__BabelMetalDevice'] =deviceName bUseMappedMemory=0 @@ -740,6 +915,8 @@ def ForwardSimple(cwvnb,center,ds,u0,rf,MaxDistance=-1.0,u0step=0,MacOsPlatform= if sys.platform == "darwin": if MacOsPlatform=='Metal': return ForwardSimpleMetal(cwvnb,center,ds,u0,rf,deviceMetal,MaxDistance=MaxDistance,u0step=u0step) + elif MacOsPlatform == 'MLX': + return ForwardSimpleMLX(cwvnb,center,ds,u0,rf,MaxDistance=MaxDistance,u0step=u0step) else: return ForwardSimpleOpenCL(cwvnb,center,ds,u0,rf,MaxDistance=MaxDistance,u0step=u0step) else: diff --git a/BabelViscoFDTD/tools/rayleigh.cpp b/BabelViscoFDTD/tools/rayleigh.cpp new file mode 100644 index 0000000..0b6a0cd --- /dev/null +++ b/BabelViscoFDTD/tools/rayleigh.cpp @@ -0,0 +1,84 @@ +#ifdef _OPENCL +__kernel void ForwardPropagationKernel(const int mr2, + const FloatingType c_wvnb_real, + const FloatingType c_wvnb_imag, + const FloatingType MaxDistance, + const int mr1, + __global const FloatingType *r2pr, + __global const FloatingType *r1pr, + __global const FloatingType *a1pr, + __global const FloatingType *u1_real, + __global const FloatingType *u1_imag, + __global FloatingType *py_data_u2_real, + __global FloatingType *py_data_u2_imag, + const int mr1step + ) +{ + int si2 = get_global_id(0); // Grid is a "flatten" 1D, thread blocks are 1D +#endif +#ifdef _MLX + int si2 = thread_position_in_grid.x; +#endif + FloatingType dx,dy,dz,R,r2x,r2y,r2z; + FloatingType temp_r,tr ; + FloatingType temp_i,ti,pCos,pSin ; + + int offset = mr1step*si2; + + // Ensure index is less than number of detection points + if (si2 < mr2) + { + // Temp variables for real and imag values + temp_r = 0; + temp_i = 0; + + // Detection point x, y, and z coordinates for specific index + r2x=r2pr[si2*3]; + r2y=r2pr[si2*3+1]; + r2z=r2pr[si2*3+2]; + + // loop through each tx element/source point + for (int si1=0; si1 0.0 && R > MaxDistance) continue; + + // Start of Rayleigh Integral calculation + ti=(exp(R*c_wvnb_imag)*a1pr[si1]/R); + tr=ti; + + // Calculate sin and cosine values of distance * real sound speed + #if defined(_METAL) || defined(_MLX) + pSin=sincos(R*c_wvnb_real,pCos); + #else + pSin=sincos(R*c_wvnb_real,ppCos); + #endif + + // Real and imaginary terms of rayleigh integral + tr*=(u1_real[si1+offset]*pCos+u1_imag[si1+offset]*pSin); + ti*=(u1_imag[si1+offset]*pCos-u1_real[si1+offset]*pSin); + + // Summate real and imaginary terms + temp_r += tr; + temp_i += ti; + } + + // Final cumulative real and imaginary pressure at detection point + R = temp_r; + + temp_r = -temp_r*c_wvnb_imag-temp_i*c_wvnb_real; + temp_i = R*c_wvnb_real-temp_i*c_wvnb_imag; + + py_data_u2_real[si2]=temp_r/(2*pi); + py_data_u2_imag[si2]=temp_i/(2*pi); + } +} \ No newline at end of file diff --git a/BabelViscoFDTD/tools/rayleighAndBHTE.hpp b/BabelViscoFDTD/tools/rayleighAndBHTE.hpp new file mode 100644 index 0000000..8508ca1 --- /dev/null +++ b/BabelViscoFDTD/tools/rayleighAndBHTE.hpp @@ -0,0 +1,12 @@ +#if defined(_METAL) || defined(_MLX) +#include +using namespace metal; +#define pi M_PI_F +#endif + +#if !defined(_METAL) && !defined(_MLX) +#define pi 3.141592653589793 +#define ppCos &pCos +#endif + +typedef float FloatingType; From d32c27fa28aab00bf875ef32b7813a75e0a5370d Mon Sep 17 00:00:00 2001 From: Alan Coreas Date: Fri, 3 Jan 2025 09:59:38 -0700 Subject: [PATCH 2/7] Update Rayleigh to loop for all backends + added mlx for BHTE function --- BabelViscoFDTD/tools/BHTE.cpp | 157 +++ BabelViscoFDTD/tools/RayleighAndBHTE.py | 1169 +++++++---------- BabelViscoFDTD/tools/rayleigh.cpp | 46 +- BabelViscoFDTD/tools/rayleighAndBHTE.hpp | 12 + .../Tools -1 - Rayleigh Integral.ipynb | 231 ++-- 5 files changed, 814 insertions(+), 801 deletions(-) create mode 100644 BabelViscoFDTD/tools/BHTE.cpp diff --git a/BabelViscoFDTD/tools/BHTE.cpp b/BabelViscoFDTD/tools/BHTE.cpp new file mode 100644 index 0000000..2a08b57 --- /dev/null +++ b/BabelViscoFDTD/tools/BHTE.cpp @@ -0,0 +1,157 @@ +// Bioheat Transfer Equation + +#ifdef _CUDA +extern "C" __global__ void BHTEFDTDKernel(float *d_output, + float *d_output2, + const float *d_input, + const float *d_input2, + const float *d_bhArr, + const float *d_perfArr, + const unsigned int *d_labels, + const float *d_Qarr, + const unsigned int *d_pointsMonitoring, + const float CoreTemp, + const int sonication, + const int outerDimx, + const int outerDimy, + const int outerDimz, + const float dt, + float *d_MonitorSlice, + float *d_Temppoints, + const int TotalStepsMonitoring, + const int nFactorMonitoring, + const int n_Step, + const int SelJ, + const unsigned int StartIndexQ, + const unsigned TotalSteps) +{ + const int gtidx = (blockIdx.x * blockDim.x + threadIdx.x); + const int gtidy = (blockIdx.y * blockDim.y + threadIdx.y); + const int gtidz = (blockIdx.z * blockDim.z + threadIdx.z); +#endif +#ifdef _OPENCL +__kernel void BHTEFDTDKernel(__global float *d_output, + __global float *d_output2, + __global const float *d_input, + __global const float *d_input2, + __global const float *d_bhArr, + __global const float *d_perfArr, + __global const unsigned int *d_labels, + __global const float *d_Qarr, + __global const unsigned int *d_pointsMonitoring, + const float CoreTemp, + const unsigned int sonication, + const unsigned int outerDimx, + const unsigned int outerDimy, + const unsigned int outerDimz, + const float dt, + __global float *d_MonitorSlice, + __global float *d_Temppoints, + const unsigned int TotalStepsMonitoring, + const unsigned int nFactorMonitoring, + const unsigned int n_Step, + const unsigned int SelJ, + const unsigned int StartIndexQ, + const unsigned TotalSteps) +{ + const int gtidx = get_global_id(0); + const int gtidy = get_global_id(1); + const int gtidz = get_global_id(2); +#endif +#ifdef _METAL +kernel void BHTEFDTDKernel(device float *d_output [[ buffer(0) ]], + device float *d_output2 [[ buffer(1) ]], + device const float *d_input [[ buffer(2) ]], + device const float *d_input2 [[ buffer(3) ]], + device const float *d_bhArr [[ buffer(4) ]], + device const float *d_perfArr [[ buffer(5) ]], + device const unsigned int *d_labels [[ buffer(6) ]], + device const float *d_Qarr [[ buffer(7) ]], + device const unsigned int *d_pointsMonitoring [[ buffer(8) ]], + device float *d_MonitorSlice [[ buffer(9) ]], + device float *d_Temppoints [[ buffer(10) ]], + constant float * floatParams [[ buffer(11) ]], + constant unsigned int * intparams [[ buffer(12) ]], + uint gid[[thread_position_in_grid]]) +{ +#endif +#if defined(_METAL) || defined(_MLX) + #ifdef _MLX + uint gid = thread_position_in_grid.x; + #endif + + #define CoreTemp floatParams[0] + #define dt floatParams[1] + #define sonication intparams[0] + #define outerDimx intparams[1] + #define outerDimy intparams[2] + #define outerDimz intparams[3] + #define TotalStepsMonitoring intparams[4] + #define nFactorMonitoring intparams[5] + #define n_Step intparams[6] + #define SelJ intparams[7] + #define StartIndexQ intparams[8] + #define TotalSteps intparams[9] + const int gtidx = gid/(outerDimy*outerDimz); + const int gtidy = (gid - gtidx*outerDimy*outerDimz)/outerDimz; + const int gtidz = gid - gtidx*outerDimy*outerDimz - gtidy*outerDimz; +#endif + + #define Tref 43.0 + unsigned int DzDy = outerDimz*outerDimy; + unsigned int coord = gtidx * DzDy + gtidy * outerDimz + gtidz; + + float R1,R2,dtp; + if(gtidx > 0 && gtidx < outerDimx-1 && gtidy > 0 && gtidy < outerDimy-1 && gtidz > 0 && gtidz < outerDimz-1) + { + + const unsigned int label = d_labels[coord]; + + d_output[coord] = d_input[coord] + d_bhArr[label] * + (d_input[coord + 1] + d_input[coord - 1] + d_input[coord + outerDimz] + d_input[coord - outerDimz] + + d_input[coord + DzDy] + d_input[coord - DzDy] - 6.0 * d_input[coord]) + + + d_perfArr[label] * (CoreTemp - d_input[coord]); + if (sonication) + { + d_output[coord] += d_Qarr[coord+StartIndexQ]; + } + + R2 = (d_output[coord] >= Tref)?0.5:0.25; + R1 = (d_input[coord] >= Tref)?0.5:0.25; + + if(fabs(d_output[coord]-d_input[coord])<0.0001) + { + d_output2[coord] = d_input2[coord] + dt * pow((float)R1,(float)(Tref-d_input[coord])); + } + else + { + if(R1 == R2) + { + d_output2[coord] = d_input2[coord] + (pow((float)R2,(float)(Tref-d_output[coord])) - pow((float)R1,(float)(Tref-d_input[coord]))) / + ( -(d_output[coord]-d_input[coord])/ dt * log(R1)); + } + else + { + dtp = dt * (Tref - d_input[coord])/(d_output[coord] - d_input[coord]); + + d_output2[coord] = d_input2[coord] + (1 - pow((float)R1,(float)(Tref-d_input[coord]))) / (- (Tref - d_input[coord])/ dtp * log(R1)) + + (pow((float)R2,(float)(Tref-d_output[coord])) - 1) / (-(d_output[coord] - Tref)/(dt - dtp) * log(R2)); + } + } + + if (gtidy==SelJ && (n_Step % nFactorMonitoring ==0)) + { + d_MonitorSlice[gtidx*outerDimz*TotalStepsMonitoring+gtidz*TotalStepsMonitoring+ n_Step/nFactorMonitoring] =d_output[coord]; + } + + if (d_pointsMonitoring[coord]>0) + { + d_Temppoints[TotalSteps*(d_pointsMonitoring[coord]-1)+n_Step]=d_output[coord]; + } + } + else if(gtidx < outerDimx && gtidy < outerDimy && gtidz < outerDimz) + { + d_output[coord] = d_input[coord]; + d_output2[coord] = d_input2[coord]; + } +} \ No newline at end of file diff --git a/BabelViscoFDTD/tools/RayleighAndBHTE.py b/BabelViscoFDTD/tools/RayleighAndBHTE.py index e1db321..e6db5b5 100644 --- a/BabelViscoFDTD/tools/RayleighAndBHTE.py +++ b/BabelViscoFDTD/tools/RayleighAndBHTE.py @@ -31,369 +31,12 @@ def resource_path(): # needed for bundling return bundle_dir -KernelCoreSourceBHTE=""" - #define Tref 43.0 - unsigned int DzDy=outerDimz*outerDimy; - unsigned int coord = gtidx*DzDy + gtidy*outerDimz + gtidz; - - float R1,R2,dtp; - if(gtidx > 0 && gtidx < outerDimx-1 && gtidy > 0 && gtidy < outerDimy-1 && gtidz > 0 && gtidz < outerDimz-1) - { - - const unsigned int label = d_labels[coord]; - - d_output[coord] = d_input[coord] + d_bhArr[label] * ( - d_input[coord + 1] + d_input[coord - 1] + d_input[coord + outerDimz] + d_input[coord - outerDimz] + - d_input[coord + DzDy] + d_input[coord - DzDy] - 6.0 * d_input[coord]) + - + d_perfArr[label] * (CoreTemp - d_input[coord]) ; - if (sonication) - { - d_output[coord]+=d_Qarr[coord+StartIndexQ]; - } - - R2 = (d_output[coord] >= Tref)?0.5:0.25; - R1 = (d_input[coord] >= Tref)?0.5:0.25; - - if(fabs(d_output[coord]-d_input[coord])<0.0001) - { - d_output2[coord] = d_input2[coord] + dt * pow((float)R1,(float)(Tref-d_input[coord])); - } - else - { - if(R1 == R2) - { - d_output2[coord] = d_input2[coord] + (pow((float)R2,(float)(Tref-d_output[coord])) - pow((float)R1,(float)(Tref-d_input[coord]))) / - ( -(d_output[coord]-d_input[coord])/ dt * log(R1)); - } - else - { - dtp = dt * (Tref - d_input[coord])/(d_output[coord] - d_input[coord]); - - d_output2[coord] = d_input2[coord] + (1 - pow((float)R1,(float)(Tref-d_input[coord]))) / (- (Tref - d_input[coord])/ dtp * log(R1)) + - (pow((float)R2,(float)(Tref-d_output[coord])) - 1) / (-(d_output[coord] - Tref)/(dt - dtp) * log(R2)); - } - } - - if (gtidy==SelJ && (n_Step % nFactorMonitoring ==0)) - { - d_MonitorSlice[gtidx*outerDimz*TotalStepsMonitoring+gtidz*TotalStepsMonitoring+ n_Step/nFactorMonitoring] =d_output[coord]; - } - - if (d_pointsMonitoring[coord]>0) - { - d_Temppoints[TotalSteps*(d_pointsMonitoring[coord]-1)+n_Step]=d_output[coord]; - } - } - else if(gtidx < outerDimx && gtidy < outerDimy && gtidz < outerDimz){ - d_output[coord] = d_input[coord]; - d_output2[coord] = d_input2[coord]; - - } - -} -""" - -import pyopencl as cl - -RayleighOpenCLMetalSource=""" -#ifdef _METAL -#include -using namespace metal; -#endif - -#define pi 3.141592653589793 -#define ppCos &pCos - -typedef float FloatingType; - -#ifdef _OPENCL -__kernel void ForwardPropagationKernel( const int mr2, - const FloatingType c_wvnb_real, - const FloatingType c_wvnb_imag, - const FloatingType MaxDistance, - const int mr1, - __global const FloatingType *r2pr, - __global const FloatingType *r1pr, - __global const FloatingType *a1pr, - __global const FloatingType *u1_real, - __global const FloatingType *u1_imag, - __global FloatingType *py_data_u2_real, - __global FloatingType *py_data_u2_imag, - const int mr1step - ) - { - int si2 = get_global_id(0); // Grid is a "flatten" 1D, thread blocks are 1D - - FloatingType dx,dy,dz,R,r2x,r2y,r2z; - FloatingType temp_r,tr ; - FloatingType temp_i,ti,pCos,pSin ; - - if ( si2 < mr2) - { - temp_r = 0; - temp_i = 0; - r2x=r2pr[si2*3]; - r2y=r2pr[si2*3+1]; - r2z=r2pr[si2*3+2]; - - for (int si1=0; si10.0) - if (R>MaxDistance) - continue; - ti=(exp(R*c_wvnb_imag)*a1pr[si1]/R); - - tr=ti; - pSin=sincos(R*c_wvnb_real,ppCos); - - tr*=(u1_real[si1+mr1step*si2]*pCos+u1_imag[si1+mr1step*si2]*pSin); - ti*=(u1_imag[si1+mr1step*si2]*pCos-u1_real[si1+mr1step*si2]*pSin); - - temp_r +=tr; - temp_i +=ti; - } - - R=temp_r; - - temp_r = -temp_r*c_wvnb_imag-temp_i*c_wvnb_real; - temp_i = R*c_wvnb_real-temp_i*c_wvnb_imag; - - py_data_u2_real[si2]=temp_r/(2*pi); - py_data_u2_imag[si2]=temp_i/(2*pi); - } - } -#endif - """ - -OpenCLMetalHeaderBHTE=""" -#ifdef _OPENCL - __kernel void BHTEFDTDKernel( __global float *d_output, - __global float *d_output2, - __global const float *d_input, - __global const float *d_input2, - __global const float *d_bhArr, - __global const float *d_perfArr, - __global const unsigned int *d_labels, - __global const float *d_Qarr, - __global const unsigned int *d_pointsMonitoring, - const float CoreTemp, - const unsigned int sonication, - const unsigned int outerDimx, - const unsigned int outerDimy, - const unsigned int outerDimz, - const float dt, - __global float *d_MonitorSlice, - __global float *d_Temppoints, - const unsigned int TotalStepsMonitoring, - const unsigned int nFactorMonitoring, - const unsigned int n_Step, - const unsigned int SelJ, - const unsigned int StartIndexQ, - const unsigned TotalSteps) - { - const int gtidx = get_global_id(0); - const int gtidy = get_global_id(1); - const int gtidz = get_global_id(2); -#endif -#ifdef _METAL - kernel void BHTEFDTDKernel( device float *d_output [[ buffer(0) ]], - device float *d_output2 [[ buffer(1) ]], - device const float *d_input [[ buffer(2) ]], - device const float *d_input2 [[ buffer(3) ]], - device const float *d_bhArr [[ buffer(4) ]], - device const float *d_perfArr [[ buffer(5) ]], - device const unsigned int *d_labels [[ buffer(6) ]], - device const float *d_Qarr [[ buffer(7) ]], - device const unsigned int *d_pointsMonitoring [[ buffer(8) ]], - device float *d_MonitorSlice [[ buffer(9) ]], - device float *d_Temppoints [[ buffer(10) ]], - constant float * floatParams [[ buffer(11) ]], - constant unsigned int * intparams [[ buffer(12) ]], - uint gid[[thread_position_in_grid]]) - { - - #define CoreTemp floatParams[0] - #define dt floatParams[1] - #define sonication intparams[0] - #define outerDimx intparams[1] - #define outerDimy intparams[2] - #define outerDimz intparams[3] - #define TotalStepsMonitoring intparams[4] - #define nFactorMonitoring intparams[5] - #define n_Step intparams[6] - #define SelJ intparams[7] - #define StartIndexQ intparams[8] - #define TotalSteps intparams[9] - const int gtidx = gid/(outerDimy*outerDimz); - const int gtidy = (gid - gtidx*outerDimy*outerDimz)/outerDimz; - const int gtidz = gid - gtidx*outerDimy*outerDimz - gtidy*outerDimz; - -#endif - """ - -OpenCLKernelBHTE =OpenCLMetalHeaderBHTE + KernelCoreSourceBHTE - Platforms=None queue = None prgcl = None ctx = None - -if sys.platform == "darwin": - - import metalcomputebabel as mc - - # Loads METAL interface - os.environ['__BabelMetal'] =os.path.dirname(os.path.abspath(__file__)) - print('loading',os.path.dirname(os.path.abspath(__file__))+"/libBabelMetal.dylib") - swift_fun = ctypes.CDLL(os.path.dirname(os.path.abspath(__file__))+"/libBabelMetal.dylib") - - swift_fun.ForwardSimpleMetal.argtypes = [ - ctypes.POINTER(ctypes.c_int), - ctypes.POINTER(ctypes.c_float), - ctypes.POINTER(ctypes.c_float), - ctypes.POINTER(ctypes.c_float), - ctypes.POINTER(ctypes.c_int), - ctypes.POINTER(ctypes.c_float), - ctypes.POINTER(ctypes.c_float), - ctypes.POINTER(ctypes.c_float), - ctypes.POINTER(ctypes.c_float), - ctypes.POINTER(ctypes.c_float), - ctypes.POINTER(ctypes.c_char_p), - ctypes.POINTER(ctypes.c_float), - ctypes.POINTER(ctypes.c_float), - ctypes.POINTER(ctypes.c_int), - ctypes.POINTER(ctypes.c_int)] - - - swift_fun.PrintMetalDevices() - print("loaded Metal",str(swift_fun)) - - def StartMetaCapture(deviceName='M1'): - os.environ['__BabelMetalDevice'] =deviceName - swift_fun.StartCapture() - - def Stopcapture(): - swift_fun.Stopcapture() - -else: - - prgcuda = None - - import cupy as cp - - RayleighCUDASource=""" - #include - - #define pi 3.141592653589793 - - typedef float FloatingType; - - #define MAX_ELEMS_IN_CONSTANT 2730 // the total constant memory can't be greater than 64k bytes - - - __device__ __forceinline__ complex cuexpf (complex z) - - { - float res_i,res_r; - sincosf(z.imag(), &res_i, &res_r); - return expf (z.real())*complex (res_r,res_i);; - } - - extern "C" __global__ void ForwardPropagationKernel(int mr2, - complex c_wvnb, - FloatingType MaxDistance, - FloatingType *r2pr, - FloatingType *r1pr, - FloatingType *a1pr, - complex * u1complex, - complex *py_data_u2, - int mr1, - int mr1step) - { - const int si2 = (blockIdx.y*gridDim.x + blockIdx.x)*blockDim.x + threadIdx.x ; // Grid is a "flatten" 1D, thread blocks are 1D - - complex cj=complex(0.0,1); - complex temp,temp2; - - FloatingType dx,dy,dz,R,r2x,r2y,r2z; - if ( si2 < mr2) - { - temp*=0; - - - r2x=r2pr[si2*3]; - r2y=r2pr[si2*3+1]; - r2z=r2pr[si2*3+2]; - - for (int si1=0; si10.0) - if (R>MaxDistance) - continue; - - temp2=cj*c_wvnb; - temp2=temp2*(-R); - temp2=cuexpf(temp2); - temp2=temp2*u1complex[si1+mr1step*si2]; - temp2=temp2*a1pr[si1]/R; - temp=temp+temp2; - } - - temp2=cj*c_wvnb; - temp=temp*temp2; - - py_data_u2[si2]=temp/((float)(2*pi)); - - } - } - """ - - CUDAHeaderBHTE=""" - - extern "C" __global__ void BHTEFDTDKernel( float *d_output, - float *d_output2, - const float *d_input, - const float *d_input2, - const float *d_bhArr, - const float *d_perfArr, - const unsigned int *d_labels, - const float *d_Qarr, - const unsigned int *d_pointsMonitoring, - const float CoreTemp, - const int sonication, - const int outerDimx, - const int outerDimy, - const int outerDimz, - const float dt, - float *d_MonitorSlice, - float *d_Temppoints, - const int TotalStepsMonitoring, - const int nFactorMonitoring, - const int n_Step, - const int SelJ, - const unsigned int StartIndexQ, - const unsigned TotalSteps) - { - const int gtidx = (blockIdx.x * blockDim.x + threadIdx.x); - const int gtidy = (blockIdx.y * blockDim.y + threadIdx.y); - const int gtidz = (blockIdx.z * blockDim.z + threadIdx.z); - """ - +clp = None def SpeedofSoundWater(Temperature): Xcoeff = [0.00000000314643 ,-0.000001478,0.000334199,-0.0580852,5.03711,1402.32] @@ -506,74 +149,15 @@ def GenerateFocusTx(f,Foc,Diam,c,PPWSurface=4): Tx = GenerateSurface(lstep,Diam,Foc) return Tx -def InitCuda(DeviceName=None): +def InitRayleighAndBHTE(DeviceName=None,gpu_backend=None): global prgcuda - devCount = cp.cuda.runtime.getDeviceCount() - if devCount == 0: - raise SystemError("There are no CUDA devices.") - - if DeviceName is not None: - selDevice = None - for deviceID in range(0, devCount): - d=cp.cuda.runtime.getDeviceProperties(deviceID) - if DeviceName in d['name'].decode('UTF-8'): - selDevice=cp.cuda.Device(deviceID) - break - selDevice.use() - AllCudaCode=RayleighCUDASource + CUDAHeaderBHTE + KernelCoreSourceBHTE - prgcuda = cp.RawModule(code= AllCudaCode) - -def InitOpenCL(DeviceName='AMD'): global Platforms global queue global prgcl global ctx - - Platforms=cl.get_platforms() - if len(Platforms)==0: - raise SystemError("No OpenCL platforms") - SelDevice=None - for device in Platforms[0].get_devices(): - print(device.name) - if DeviceName in device.name: - SelDevice=device - if SelDevice is None: - raise SystemError("No OpenCL device containing name [%s]" %(DeviceName)) - else: - print('Selecting device: ', SelDevice.name) - ctx = cl.Context([SelDevice]) - queue = cl.CommandQueue(ctx) - prgcl = cl.Program(ctx, '#define _OPENCL\n'+RayleighOpenCLMetalSource+OpenCLKernelBHTE).build() - -def InitMetal(DeviceName='AMD'): - global ctx - global prgcl - - devices = mc.get_devices() - SelDevice=None - for n,dev in enumerate(devices): - if DeviceName in dev.deviceName: - SelDevice=dev - break - if SelDevice is None: - raise SystemError("No Metal device containing name [%s]" %(DeviceName)) - else: - print('Selecting device: ', dev.deviceName) - - ctx = mc.Device(n) - print(ctx) - if 'arm64' not in platform.platform(): - ctx.set_external_gpu(1) - prgcl = ctx.kernel('#define _METAL\n'+RayleighOpenCLMetalSource+OpenCLKernelBHTE) - -def InitMLX(DeviceName='AMD'): - global ctx - global prgcl_mlx - global clp_mlx + global clp global sel_device - - import mlx.core as mx - clp_mlx = mx + global swift_fun kernel_files = [ os.path.join(resource_path(), 'rayleigh.cpp'), @@ -582,8 +166,63 @@ def InitMLX(DeviceName='AMD'): with open(os.path.join(resource_path(), 'rayleighAndBHTE.hpp'), 'r') as f: header = f.read() + + if gpu_backend == 'CUDA': + import cupy as cp + clp = cp + + preamble = '#define _CUDA\n' + ctx,prgcuda,sel_device = InitCUDA(DeviceName,preamble+header,kernel_files) + + elif gpu_backend == 'OpenCL': + import pyopencl as pocl + clp = pocl + + preamble = '#define _OPENCL\n' + queue,prgcl,ctx = InitOpenCL(DeviceName,preamble+header,kernel_files) - kernel_functions = [{'name': 'ForwardPropagationKernel', + elif gpu_backend == 'Metal': + import metalcomputebabel as mc + + # Loads METAL interface + os.environ['__BabelMetal'] =os.path.dirname(os.path.abspath(__file__)) + print('loading',os.path.dirname(os.path.abspath(__file__))+"/libBabelMetal.dylib") + swift_fun = ctypes.CDLL(os.path.dirname(os.path.abspath(__file__))+"/libBabelMetal.dylib") + + swift_fun.ForwardSimpleMetal.argtypes = [ + ctypes.POINTER(ctypes.c_int), + ctypes.POINTER(ctypes.c_float), + ctypes.POINTER(ctypes.c_float), + ctypes.POINTER(ctypes.c_float), + ctypes.POINTER(ctypes.c_int), + ctypes.POINTER(ctypes.c_float), + ctypes.POINTER(ctypes.c_float), + ctypes.POINTER(ctypes.c_float), + ctypes.POINTER(ctypes.c_float), + ctypes.POINTER(ctypes.c_float), + ctypes.POINTER(ctypes.c_char_p), + ctypes.POINTER(ctypes.c_float), + ctypes.POINTER(ctypes.c_float), + ctypes.POINTER(ctypes.c_int), + ctypes.POINTER(ctypes.c_int)] + + + swift_fun.PrintMetalDevices() + print("loaded Metal",str(swift_fun)) + + def StartMetaCapture(deviceName='M1'): + os.environ['__BabelMetalDevice'] =deviceName + swift_fun.StartCapture() + + def Stopcapture(): + swift_fun.Stopcapture() + else: + assert(gpu_backend=='MLX') + import mlx.core as mx + clp = mx + + preamble = '#define _MLX\n' + kernel_functions = [{'name': 'ForwardPropagationKernel', 'file': kernel_files[0], 'input_names': ["mr2", "c_wvnb_real", "c_wvnb_imag", "MaxDistance", "mr1", "r2pr","r1pr","a1pr","u1_real","u1_imag","mr1step"], @@ -597,21 +236,115 @@ def InitMLX(DeviceName='AMD'): "SelJ","StartIndexQ","TotalSteps"], 'output_names': ["d_output","d_output2","d_MonitorSlice","d_Temppoints"], 'atomic_outputs': False}] + + prgcl,sel_device = InitMLX(DeviceName,preamble+header,kernel_functions,build_later=True) + + +def InitCUDA(DeviceName=None,header='',kernel_files=None,build_later=False): + import cupy as cp + # global prgcuda + + # Obtain list of gpu devices + devCount = cp.cuda.runtime.getDeviceCount() + if devCount == 0: + raise SystemError("There are no CUDA devices.") + + # Select device that matches specified name + if DeviceName is not None: + selDevice = None + for deviceID in range(0, devCount): + d=cp.cuda.runtime.getDeviceProperties(deviceID) + if DeviceName in d['name'].decode('UTF-8'): + selDevice=cp.cuda.Device(deviceID) + break + selDevice.use() + + # Combine kernel codes with header + kernel_codes = [header] + for k_file in kernel_files: + with open(k_file, 'r') as f: + kernel_code = f.read() + kernel_codes.append(kernel_code) + complete_kernel = '\n'.join(kernel_codes) + + # Build program later, return complete kernel for now + if build_later: + prgcuda = complete_kernel + # Build program from source code + else: + + # Windows sometimes has issues finding CUDA + if platform.system()=='Windows': + sys.executable.split('\\')[:-1] + options=('-I',os.path.join(os.getenv('CUDA_PATH'),'Library','Include'), + '-I',str(resource_path()), + '--ptxas-options=-v') + else: + options=('-I',str(resource_path())) + + prgcuda = cp.RawModule(code=complete_kernel,options=options) + + return ctx,prgcuda,selDevice + +def InitOpenCL(DeviceName='AMD',header='',kernel_files=None,build_later=False): + import pyopencl as pocl + + # Obtain list of openCL platforms + Platforms=pocl.get_platforms() + if len(Platforms)==0: + raise SystemError("No OpenCL platforms") + + # Obtain list of available devices and select one + SelDevice=None + for device in Platforms[0].get_devices(): + print(device.name) + if DeviceName in device.name: + SelDevice=device + if SelDevice is None: + raise SystemError("No OpenCL device containing name [%s]" %(DeviceName)) + else: + print('Selecting device: ', SelDevice.name) + + # Create context for selected device + ctx = pocl.Context([SelDevice]) + + # Combine kernel codes with header + kernel_codes = [header] + for k_file in kernel_files: + with open(k_file, 'r') as f: + kernel_code = f.read() + kernel_codes.append(kernel_code) + complete_kernel = '\n'.join(kernel_codes) + + # Build program later, return complete kernel for now + if build_later: + prgcl = complete_kernel + # Build program from source code + else: + prgcl = pocl.Program(ctx,complete_kernel).build() - preamble = '#define _MLX\n' - build_later = True + # Create command queue for selected device + queue = pocl.CommandQueue(ctx) + + # Allocate device memory + mf = pocl.mem_flags + + return queue, prgcl, ctx + +def InitMLX(DeviceName='AMD',header='',kernel_files=None,build_later=False): + import mlx.core as mx sel_device = mx.default_device() print('Selecting device: ', sel_device) kernels = {} - for kf in kernel_functions: + for kf in kernel_files: with open(kf['file'], 'r') as f: lines = f.readlines() kernel_code = ''.join(lines[:-1]) # Remove last bracket if build_later: - kf['header'] = preamble+header + kf['header'] = header kf['source'] = kernel_code kernels[kf['name']] = kf else: @@ -619,125 +352,25 @@ def InitMLX(DeviceName='AMD'): input_names = kf['input_names'], output_names = kf['output_names'], atomic_outputs = kf['atomic_outputs'], - header = preamble + header, + header = header, source = kernel_code) kernels[kf['name']] = kernel - prgcl_mlx = kernels - -def ForwardSimpleCUDA(cwvnb,center,ds,u0,rf,MaxDistance=-1.0,u0step=0): - if u0step!=0: - mr1=u0step - assert(mr1*rf.shape[0]==u0.shape[0]) - assert(mr1==center.shape[0]) - assert(mr1==ds.shape[0]) - else: - mr1=center.shape[0] - mr2=rf.shape[0] - - d_r2pr= cp.asarray(rf) - d_centerpr= cp.asarray(center) - d_dspr= cp.asarray(ds) - d_u0complex= cp.asarray(u0) - d_u2complex= cp.zeros(rf.shape[0],cp.complex64) - - CUDA_THREADBLOCKLENGTH = 512 - MAX_ELEMS_IN_CONSTANT = 2730 - dimThreadBlock=(CUDA_THREADBLOCKLENGTH, 1,1) - - nBlockSizeX = int((mr2 - 1) / dimThreadBlock[0]) + 1 - nBlockSizeY = 1 - - if (nBlockSizeX > 65534 ): - nBlockSizeY = int(nBlockSizeX / 65534) - if (nBlockSizeX %65534 !=0): - nBlockSizeY+=1 - nBlockSizeX = int(nBlockSizeX/nBlockSizeY)+1 - - dimBlockGrid=(nBlockSizeX, nBlockSizeY,1) - - ForwardPropagationKernel = prgcuda.get_function("ForwardPropagationKernel") - - - ForwardPropagationKernel(dimBlockGrid, - dimThreadBlock, - (mr2, - cwvnb, - MaxDistance, - d_r2pr, - d_centerpr, - d_dspr, - d_u0complex, - d_u2complex, - mr1, - u0step)) - - u2=d_u2complex.get() - return u2 - -def ForwardSimpleOpenCL(cwvnb,center,ds,u0,rf,MaxDistance=-1.0,u0step=0): - global queue - global prg - global ctx - mf = cl.mem_flags - - if u0step!=0: - mr1=u0step - assert(mr1*rf.shape[0]==u0.shape[0]) - assert(mr1==center.shape[0]) - assert(mr1==ds.shape[0]) - else: - mr1=center.shape[0] - - - d_r2pr = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=rf) - d_r1pr = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=center) - d_u1realpr=cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=np.real(u0).copy()) - d_u1imagpr=cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=np.imag(u0).copy()) - d_a1pr = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=ds) - - - u2_real = np.zeros((rf.shape[0]),dtype=np.float32) - u2_imag = np.zeros((rf.shape[0]),dtype=np.float32) - - d_u2realpr = cl.Buffer(ctx, mf.WRITE_ONLY, u2_real.nbytes) - d_u2imagpr = cl.Buffer(ctx, mf.WRITE_ONLY, u2_real.nbytes) - + return kernels,sel_device +def ForwardSimple(cwvnb,center,ds,u0,rf,MaxDistance=-1.0,u0step=0,gpu_backend='OpenCL',gpu_device='M1'): #MacOsPlatform='Metal',deviceMetal='6800'): + ''' + MAIN function to call for ForwardRayleigh , returns the complex values of particle speed + cwvnb is the complex speed of sound + center is an [Mx3] array with the position of the decomposed transducer elements + ds is an [M] array with the transducer element surface area + u0 is [M] complex array with the particle speed at eact transducer element + rf is [Nx3] array with the positons where Rayleigh integral will be calculated - knl = prgcl.ForwardPropagationKernel # Use this Kernel object for repeated calls - if u2_real.shape[0] % 64 ==0: - ks=[u2_real.shape[0]] - else: - ks=[int(u2_real.shape[0]/64)*64+64] - knl(queue, ks, [64], - np.int32(rf.shape[0]), - np.float32(np.real(cwvnb)), - np.float32(np.imag(cwvnb)), - np.float32(MaxDistance), - np.int32(mr1), - d_r2pr, - d_r1pr, - d_a1pr, - d_u1realpr, - d_u1imagpr, - d_u2realpr, - d_u2imagpr, - np.int32(u0step)) + Function returns a [N] complex array of particle speed at locations rf - cl.enqueue_copy(queue, u2_real,d_u2realpr) - cl.enqueue_copy(queue, u2_imag,d_u2imagpr) - u2=u2_real+1j*u2_imag - - return u2 - -def ForwardSimpleMLX(cwvnb,center,ds,u0,rf,MaxDistance=-1.0,u0step=0): - global queue - global prgcl_mlx - global ctx - global clp_mlx - global sel_device + ''' mr2=rf.shape[0] @@ -749,137 +382,89 @@ def ForwardSimpleMLX(cwvnb,center,ds,u0,rf,MaxDistance=-1.0,u0step=0): else: mr1=center.shape[0] - # Change to mlx arrays - d_r1pr = clp_mlx.array(center) - d_u1realpr=clp_mlx.array(np.real(u0)) - d_u1imagpr=clp_mlx.array(np.imag(u0)) - d_a1pr = clp_mlx.array(ds) - u2 = np.zeros((rf.shape[0]),dtype=np.complex64) - # Build program from source code - knl = clp_mlx.fast.metal_kernel(name = f"{prgcl_mlx['ForwardPropagationKernel']['name']}", - input_names = prgcl_mlx['ForwardPropagationKernel']['input_names'], - output_names = prgcl_mlx['ForwardPropagationKernel']['output_names'], - source = prgcl_mlx['ForwardPropagationKernel']['source'], - header = prgcl_mlx['ForwardPropagationKernel']['header'], - atomic_outputs = prgcl_mlx['ForwardPropagationKernel']['atomic_outputs'], - ) - - # We need to split in small chunks to be sure the kernel does not take too much time - # otherwise the OS will kill it - - NonBlockingstep = int(24000e6) - step = int(NonBlockingstep/mr1) - - if step > mr2: - step = mr2 - if step < 5: - step = 5 - - slice_start = 0 - slice_end = 0 - - while slice_start < mr2: + # Setup for kernel call + if gpu_backend == 'CUDA': + # Transfer input data to gpu + d_r2pr= clp.asarray(rf) + d_centerpr= clp.asarray(center) + d_dspr= clp.asarray(ds) + d_u0complex= clp.asarray(u0) - slice_end = min(slice_start + step, mr2) - chunk_size = slice_end-slice_start + # Get kernel function + ForwardPropagationKernel = prgcuda.get_function("ForwardPropagationKernel") - print(f"Working on slices {slice_start} to {slice_end} out of {u2.shape[0]}") - - # Grab section of data - d_r2pr = clp_mlx.array(rf[slice_start:slice_end,:]) - d_u2realpr = clp_mlx.zeros_like(d_r2pr) - d_u2imagpr = clp_mlx.zeros_like(d_r2pr) + # Determine kernel call dimensions + CUDA_THREADBLOCKLENGTH = 512 + MAX_ELEMS_IN_CONSTANT = 2730 + dimThreadBlock=(CUDA_THREADBLOCKLENGTH, 1,1) - # Deploy kernel - d_u2realpr,d_u2imagpr = knl(inputs = [np.int32(slice_end), - np.float32(np.real(cwvnb)), - np.float32(np.imag(cwvnb)), - np.float32(MaxDistance), - np.int32(mr1), - d_r2pr, - d_r1pr, - d_a1pr, - d_u1realpr, - d_u1imagpr, - np.int32(u0step)], - output_shapes = [(chunk_size,),(chunk_size,)], - output_dtypes = [clp_mlx.float32,clp_mlx.float32], - grid=(chunk_size,1,1), - threadgroup=(256, 1, 1), - verbose=False, - stream=sel_device) - clp_mlx.synchronize() - - # Change back to numpy array - u2_real = np.array(d_u2realpr) - u2_imag = np.array(d_u2imagpr) + nBlockSizeX = int((mr2 - 1) / dimThreadBlock[0]) + 1 + nBlockSizeY = 1 + + if (nBlockSizeX > 65534 ): + nBlockSizeY = int(nBlockSizeX / 65534) + if (nBlockSizeX %65534 !=0): + nBlockSizeY+=1 + nBlockSizeX = int(nBlockSizeX/nBlockSizeY)+1 - # Combine real & imag parts - u2_section = u2_real+1j*u2_imag + dimBlockGrid=(nBlockSizeX, nBlockSizeY,1) - # Update final array - u2[slice_start:slice_end] = u2_section + elif gpu_backend == 'OpenCL': + mf = clp.mem_flags - # Clean up mlx arrays - del d_u2realpr,d_u2imagpr - gc.collect() + # Transfer input data to gpu + d_r2pr = clp.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=rf) + d_r1pr = clp.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=center) + d_u1realpr=clp.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=np.real(u0).copy()) + d_u1imagpr=clp.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=np.imag(u0).copy()) + d_a1pr = clp.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=ds) - # Update starting location - slice_start += step - - return u2 - -def ForwardSimpleMetal(cwvnb,center,ds,u0,rf,deviceName,MaxDistance=-1.0,u0step=0): - os.environ['__BabelMetalDevice'] =deviceName - bUseMappedMemory=0 - if np.__version__ >="1.22.0": - if 'arm64' in platform.platform() and\ - np.core.multiarray.get_handler_name(center)=="page_data_allocator": - bUseMappedMemory=1 - #We assume arrays were allocated with page_data_allocator to have aligned date + # Get kernel Function + knl = prgcl.ForwardPropagationKernel - - mr2=np.array([rf.shape[0]]) - - if u0step!=0: - mr1=u0step - assert(mr1*rf.shape[0]==u0.shape[0]) - assert(mr1==center.shape[0]) - assert(mr1==ds.shape[0]) - else: - mr1=center.shape[0] - - mr1=np.array([mr1]) - u0step_a=np.array([u0step]) - MaxDistance_a=np.array([MaxDistance]).astype(np.float32) - - ibUseMappedMemory =np.array([bUseMappedMemory]) - cwvnb_real=np.array([np.real(cwvnb)]) - cwvnb_imag=np.array([np.imag(cwvnb)]) - - mr1_ptr = mr1.ctypes.data_as(ctypes.POINTER(ctypes.c_int)) - mr2_ptr = mr2.ctypes.data_as(ctypes.POINTER(ctypes.c_int)) - u0step_ptr = u0step_a.ctypes.data_as(ctypes.POINTER(ctypes.c_int)) - bUseMappedMemory_ptr =ibUseMappedMemory.ctypes.data_as(ctypes.POINTER(ctypes.c_int)) - cwvnb_real_ptr = cwvnb_real.ctypes.data_as(ctypes.POINTER(ctypes.c_float)) - cwvnb_imag_ptr = cwvnb_imag.ctypes.data_as(ctypes.POINTER(ctypes.c_float)) - cwvnb_imag_ptr = cwvnb_imag.ctypes.data_as(ctypes.POINTER(ctypes.c_float)) - MaxDistance_ptr= MaxDistance_a.ctypes.data_as(ctypes.POINTER(ctypes.c_float)) - r1_ptr=center.ctypes.data_as(ctypes.POINTER(ctypes.c_float)) - r2_ptr=rf.ctypes.data_as(ctypes.POINTER(ctypes.c_float)) - a1_ptr=ds.ctypes.data_as(ctypes.POINTER(ctypes.c_float)) - u1_real_ptr=np.real(u0).copy().ctypes.data_as(ctypes.POINTER(ctypes.c_float)) - u1_imag_ptr=np.imag(u0).copy().ctypes.data_as(ctypes.POINTER(ctypes.c_float)) - deviceName_ptr=ctypes.c_char_p(deviceName.encode()) - u2_real = np.zeros(rf.shape[0],np.float32) - u2_imag = np.zeros(rf.shape[0],np.float32) - u2_real_ptr = u2_real.ctypes.data_as(ctypes.POINTER(ctypes.c_float)) - u2_imag_ptr = u2_imag.ctypes.data_as(ctypes.POINTER(ctypes.c_float)) - - - ret = swift_fun.ForwardSimpleMetal(mr2_ptr, + elif gpu_backend == 'Metal': + # Determine if data is aligned + os.environ['__BabelMetalDevice'] = gpu_device + bUseMappedMemory=0 + if np.__version__ >="1.22.0": + if 'arm64' in platform.platform() and\ + np.core.multiarray.get_handler_name(center)=="page_data_allocator": + bUseMappedMemory=1 + #We assume arrays were allocated with page_data_allocator to have aligned date + + # Convert inputs to numpy arrrays + mr1=np.array([mr1]) + mr2=np.array([mr2]) + u0step_a=np.array([u0step]) + MaxDistance_a=np.array([MaxDistance]).astype(np.float32) + ibUseMappedMemory =np.array([bUseMappedMemory]) + cwvnb_real=np.array([np.real(cwvnb)]) + cwvnb_imag=np.array([np.imag(cwvnb)]) + + # Create pointers to use with swift call + mr1_ptr = mr1.ctypes.data_as(ctypes.POINTER(ctypes.c_int)) + mr2_ptr = mr2.ctypes.data_as(ctypes.POINTER(ctypes.c_int)) + u0step_ptr = u0step_a.ctypes.data_as(ctypes.POINTER(ctypes.c_int)) + bUseMappedMemory_ptr =ibUseMappedMemory.ctypes.data_as(ctypes.POINTER(ctypes.c_int)) + cwvnb_real_ptr = cwvnb_real.ctypes.data_as(ctypes.POINTER(ctypes.c_float)) + cwvnb_imag_ptr = cwvnb_imag.ctypes.data_as(ctypes.POINTER(ctypes.c_float)) + cwvnb_imag_ptr = cwvnb_imag.ctypes.data_as(ctypes.POINTER(ctypes.c_float)) + MaxDistance_ptr= MaxDistance_a.ctypes.data_as(ctypes.POINTER(ctypes.c_float)) + r1_ptr=center.ctypes.data_as(ctypes.POINTER(ctypes.c_float)) + r2_ptr=rf.ctypes.data_as(ctypes.POINTER(ctypes.c_float)) + a1_ptr=ds.ctypes.data_as(ctypes.POINTER(ctypes.c_float)) + u1_real_ptr=np.real(u0).copy().ctypes.data_as(ctypes.POINTER(ctypes.c_float)) + u1_imag_ptr=np.imag(u0).copy().ctypes.data_as(ctypes.POINTER(ctypes.c_float)) + deviceName_ptr=ctypes.c_char_p(gpu_device.encode()) + u2_real = np.zeros(rf.shape[0],np.float32) + u2_imag = np.zeros(rf.shape[0],np.float32) + u2_real_ptr = u2_real.ctypes.data_as(ctypes.POINTER(ctypes.c_float)) + u2_imag_ptr = u2_imag.ctypes.data_as(ctypes.POINTER(ctypes.c_float)) + + # Call swift code that already handles looping + ret = swift_fun.ForwardSimpleMetal(mr2_ptr, cwvnb_real_ptr, cwvnb_imag_ptr, MaxDistance_ptr, @@ -894,34 +479,170 @@ def ForwardSimpleMetal(cwvnb,center,ds,u0,rf,deviceName,MaxDistance=-1.0,u0step= u2_imag_ptr, bUseMappedMemory_ptr, u0step_ptr) - if ret ==1: - raise ValueError("Unable to run simulation (mostly likely name of GPU is incorrect)") - - return u2_real+1j*u2_imag + if ret ==1: + raise ValueError("Unable to run simulation (mostly likely name of GPU is incorrect)") -def ForwardSimple(cwvnb,center,ds,u0,rf,MaxDistance=-1.0,u0step=0,MacOsPlatform='Metal',deviceMetal='6800'): - ''' - MAIN function to call for ForwardRayleigh , returns the complex values of particle speed - cwvnb is the complex speed of sound - center is an [Mx3] array with the position of the decomposed transducer elements - ds is an [M] array with the transducer element surface area - u0 is [M] complex array with the particle speed at eact transducer element - rf is [Nx3] array with the positons where Rayleigh integral will be calculated + # Return result + return u2_real+1j*u2_imag + else: + assert(gpu_backend == 'MLX') + + # Convert input data to MLX arrays + d_r1pr = clp.array(center) + d_u1realpr=clp.array(np.real(u0)) + d_u1imagpr=clp.array(np.imag(u0)) + d_a1pr = clp.array(ds) + + # Get kernel function + knl = clp.fast.metal_kernel(name = f"{prgcl['ForwardPropagationKernel']['name']}", + input_names = prgcl['ForwardPropagationKernel']['input_names'], + output_names = prgcl['ForwardPropagationKernel']['output_names'], + source = prgcl['ForwardPropagationKernel']['source'], + header = prgcl['ForwardPropagationKernel']['header'], + atomic_outputs = prgcl['ForwardPropagationKernel']['atomic_outputs'], + ) - Function returns a [N] complex array of particle speed at locations rf + # Determine step size for looping + NonBlockingstep = int(24000e6) + step = int(NonBlockingstep/mr1) + + if step > mr2: + step = mr2 + if step < 5: + step = 5 - ''' - global prgcuda - if sys.platform == "darwin": - if MacOsPlatform=='Metal': - return ForwardSimpleMetal(cwvnb,center,ds,u0,rf,deviceMetal,MaxDistance=MaxDistance,u0step=u0step) - elif MacOsPlatform == 'MLX': - return ForwardSimpleMLX(cwvnb,center,ds,u0,rf,MaxDistance=MaxDistance,u0step=u0step) - else: - return ForwardSimpleOpenCL(cwvnb,center,ds,u0,rf,MaxDistance=MaxDistance,u0step=u0step) - else: - return ForwardSimpleCUDA(cwvnb,center,ds,u0,rf,MaxDistance=MaxDistance,u0step=u0step) + # Loop gpu kernel calls until all outputs have been calculated + detection_point_start_index = 0 + while detection_point_start_index < mr2: + detection_point_end_index = min(detection_point_start_index + step, mr2) + chunk_size = detection_point_end_index - detection_point_start_index + + print(f"Working on detection points {detection_point_start_index} to {detection_point_end_index} out of {u2.shape[0]}") + + # Grab section of data + data_section = rf[detection_point_start_index:detection_point_end_index,:] + + if gpu_backend == 'CUDA': + + # Output section arrays + d_u2complex= clp.zeros(data_section[:,0],clp.complex64) + + # Deploy kernel + ForwardPropagationKernel(dimBlockGrid, + dimThreadBlock, + (mr2, + cwvnb, + MaxDistance, + d_r2pr, + d_centerpr, + d_dspr, + d_u0complex, + d_u2complex, + mr1, + u0step)) + + # Change back to numpy array + u2_section = d_u2complex.get() + + # Update final array + u2[detection_point_start_index:detection_point_end_index] = u2_section + + elif gpu_backend == 'OpenCL': + + # Determine kernel call dimensions + if data_section.shape[0] % 64 ==0: + ks = [data_section.shape[0]] + else: + ks = [int(data_section.shape[0]/64)*64+64] + + # Output section arrays + u2_real = np.zeros_like(data_section[:,0]) + u2_imag = np.zeros_like(data_section[:,0]) + + # Move to gpu + d_r2pr = clp.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=data_section) + d_u2realpr = clp.Buffer(ctx, mf.WRITE_ONLY,data_section.nbytes) + d_u2imagpr = clp.Buffer(ctx, mf.WRITE_ONLY,data_section.nbytes) + + # Deploy kernel + knl(queue, ks, [64], + np.int32(rf.shape[0]), + np.float32(np.real(cwvnb)), + np.float32(np.imag(cwvnb)), + np.float32(MaxDistance), + np.int32(mr1), + d_r2pr, + d_r1pr, + d_a1pr, + d_u1realpr, + d_u1imagpr, + d_u2realpr, + d_u2imagpr, + np.int32(u0step)) + queue.finish() + # Change back to numpy array + clp.enqueue_copy(queue, u2_real,d_u2realpr) + clp.enqueue_copy(queue, u2_imag,d_u2imagpr) + + # Combine real & imag parts + u2_section = u2_real+1j*u2_imag + + # Update final array + u2[detection_point_start_index:detection_point_end_index] = u2_section + + # Clean up pocl arrays + del d_u2realpr,d_u2imagpr + gc.collect() + + elif gpu_backend == 'Metal': + pass + else: + assert(gpu_backend == 'MLX') + + # Grab section of data + d_r2pr = clp.array(data_section) + d_u2realpr = clp.zeros_like(d_r2pr[:,0]) + d_u2imagpr = clp.zeros_like(d_r2pr[:,0]) + + # Deploy kernel + d_u2realpr,d_u2imagpr = knl(inputs = [np.int32(detection_point_end_index), + np.float32(np.real(cwvnb)), + np.float32(np.imag(cwvnb)), + np.float32(MaxDistance), + np.int32(mr1), + d_r2pr, + d_r1pr, + d_a1pr, + d_u1realpr, + d_u1imagpr, + np.int32(u0step)], + output_shapes = [(chunk_size,),(chunk_size,)], + output_dtypes = [clp.float32,clp.float32], + grid=(chunk_size,1,1), + threadgroup=(256, 1, 1), + verbose=False, + stream=sel_device) + clp.synchronize() + + # Change back to numpy array + u2_real = np.array(d_u2realpr) + u2_imag = np.array(d_u2imagpr) + + # Combine real & imag parts + u2_section = u2_real+1j*u2_imag + + # Update final array + u2[detection_point_start_index:detection_point_end_index] = u2_section + + # Clean up mlx arrays + del d_u2realpr,d_u2imagpr + gc.collect() + + # Update starting location + detection_point_start_index += step + + return u2 def getBHTECoefficient( kappa,rho,c_t,h,t_int,dt=0.1): @@ -1248,6 +969,90 @@ def BHTE(Pressure,MaterialMap,MaterialList,dx, MonitorSlice=d_MonitorSlice.get() TemperaturePoints=d_TemperaturePoints.get() + elif Backend=='MLX': + + d_perfArr=clp.array(perfArr) + d_bhArr=clp.array(bhArr) + d_Qarr=clp.array(Qarr) + d_MaterialMap=clp.array(MaterialMap) + d_T0 = clp.array(initTemp) + d_T1 = clp.array(T1) + d_Dose0 = clp.array(Dose0) + d_Dose1 = clp.array(Dose1) + d_MonitorSlice = clp.array(MonitorSlice.nbytes) + d_MonitoringPoints=clp.array(MonitoringPoints) + d_TemperaturePoints=clp.array(TemperaturePoints.nbytes) + + # Build program from source code + knl = clp.fast.metal_kernel(name = f"{prgcl['BHTEFDTDKernel']['name']}", + input_names = prgcl['BHTEFDTDKernel']['input_names'], + output_names = prgcl['BHTEFDTDKernel']['output_names'], + source = prgcl['BHTEFDTDKernel']['source'], + header = prgcl['BHTEFDTDKernel']['header'], + atomic_outputs = prgcl['BHTEFDTDKernel']['atomic_outputs'], + ) + + floatparams = np.array([stableTemp,dt],dtype=np.float32) + d_floatparams= clp.array(floatparams) + + for n in range(TotalDurationSteps): + if n c_wvnb, + FloatingType MaxDistance, + FloatingType *r2pr, + FloatingType *r1pr, + FloatingType *a1pr, + complex * u1complex, + complex *py_data_u2, + int mr1, + int mr1step) +{ + const int si2 = (blockIdx.y*gridDim.x + blockIdx.x)*blockDim.x + threadIdx.x ; // Grid is a "flatten" 1D, thread blocks are 1D +#endif #ifdef _OPENCL __kernel void ForwardPropagationKernel(const int mr2, const FloatingType c_wvnb_real, @@ -16,12 +30,18 @@ __kernel void ForwardPropagationKernel(const int mr2, { int si2 = get_global_id(0); // Grid is a "flatten" 1D, thread blocks are 1D #endif -#ifdef _MLX + #ifdef _MLX int si2 = thread_position_in_grid.x; -#endif + #endif + FloatingType dx,dy,dz,R,r2x,r2y,r2z; + #ifdef _CUDA + complex cj=complex(0.0,1); + complex temp,temp2; + #else FloatingType temp_r,tr ; FloatingType temp_i,ti,pCos,pSin ; + #endif int offset = mr1step*si2; @@ -29,8 +49,12 @@ __kernel void ForwardPropagationKernel(const int mr2, if (si2 < mr2) { // Temp variables for real and imag values + #ifdef _CUDA + temp*=0; + #else temp_r = 0; temp_i = 0; + #endif // Detection point x, y, and z coordinates for specific index r2x=r2pr[si2*3]; @@ -52,6 +76,14 @@ __kernel void ForwardPropagationKernel(const int mr2, // If distance is greater than supplied max distance, ignore calculation if (MaxDistance > 0.0 && R > MaxDistance) continue; + #ifdef _CUDA + temp2=cj*c_wvnb; + temp2=temp2*(-R); + temp2=cuexpf(temp2); + temp2=temp2*u1complex[si1+mr1step*si2]; + temp2=temp2*a1pr[si1]/R; + temp=temp+temp2; + #else // Start of Rayleigh Integral calculation ti=(exp(R*c_wvnb_imag)*a1pr[si1]/R); tr=ti; @@ -69,9 +101,16 @@ __kernel void ForwardPropagationKernel(const int mr2, // Summate real and imaginary terms temp_r += tr; - temp_i += ti; + temp_i += ti; + #endif } + #ifdef _CUDA + temp2=cj*c_wvnb; + temp=temp*temp2; + + py_data_u2[si2]=temp/((float)(2*pi)); + #else // Final cumulative real and imaginary pressure at detection point R = temp_r; @@ -80,5 +119,6 @@ __kernel void ForwardPropagationKernel(const int mr2, py_data_u2_real[si2]=temp_r/(2*pi); py_data_u2_imag[si2]=temp_i/(2*pi); + #endif } } \ No newline at end of file diff --git a/BabelViscoFDTD/tools/rayleighAndBHTE.hpp b/BabelViscoFDTD/tools/rayleighAndBHTE.hpp index 8508ca1..42e2fad 100644 --- a/BabelViscoFDTD/tools/rayleighAndBHTE.hpp +++ b/BabelViscoFDTD/tools/rayleighAndBHTE.hpp @@ -9,4 +9,16 @@ using namespace metal; #define ppCos &pCos #endif +#ifdef _CUDA +#include +#define MAX_ELEMS_IN_CONSTANT 2730 // the total constant memory can't be greater than 64k bytes + +__device__ __forceinline__ complex cuexpf (complex z) +{ + float res_i,res_r; + sincosf(z.imag(), &res_i, &res_r); + return expf (z.real())*complex (res_r,res_i);; +} +#endif + typedef float FloatingType; diff --git a/Tutorial Notebooks/Tools -1 - Rayleigh Integral.ipynb b/Tutorial Notebooks/Tools -1 - Rayleigh Integral.ipynb index 72ff01c..603de3a 100644 --- a/Tutorial Notebooks/Tools -1 - Rayleigh Integral.ipynb +++ b/Tutorial Notebooks/Tools -1 - Rayleigh Integral.ipynb @@ -4,21 +4,13 @@ "cell_type": "code", "execution_count": 1, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "loading /opt/homebrew/Caskroom/miniforge/base/envs/py310/lib/python3.10/site-packages/BabelViscoFDTD/tools/libBabelMetal.dylib\n", - "loaded Metal \n" - ] - } - ], + "outputs": [], "source": [ "%matplotlib inline\n", + "# %cd ..\n", "from matplotlib import pyplot as plt\n", "import numpy as np\n", - "from BabelViscoFDTD.tools.RayleighAndBHTE import ForwardSimple, InitCuda, InitOpenCL , GenerateFocusTx\n", + "from BabelViscoFDTD.tools.RayleighAndBHTE import ForwardSimple, InitRayleighAndBHTE,GenerateFocusTx\n", "from pprint import pprint\n", "from scipy.io import loadmat\n", "import sys\n", @@ -42,15 +34,12 @@ "\n", "For support to pyvista embedded widget, use only **jupyter notebook** as pyvista has still some limitations with jupyter-lab\n", "\n", - "## Preliminary step, initialization of GPU backend\n", - "CUDA will be activated for Windows and Linux, while OpenCL will be activated for MacOS.\n", - "\n", - "Please note that Metal support is already enabled by default and that is the default backend for MacOS\n" + "## Preliminary step, initialization of GPU backend\n" ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "metadata": { "tags": [] }, @@ -59,17 +48,17 @@ "name": "stdout", "output_type": "stream", "text": [ - "Apple M1 Max\n", - "Selecting device: Apple M1 Max\n" + "Selecting device: Device(gpu, 0)\n" ] } ], "source": [ - "if sys.platform == \"darwin\":\n", - " #if using MacOS, OpenCL interface and Metal are available, but OpenCL needs to be initialized\n", - " InitOpenCL('M1')\n", - "else:\n", - " InitCuda()" + "backend = 'MLX'\n", + "# backend = 'Metal'\n", + "# backend='OpenCL'\n", + "# backend = 'CUDA'\n", + "\n", + "InitRayleighAndBHTE('M1',gpu_backend=backend)" ] }, { @@ -82,60 +71,60 @@ }, { "cell_type": "code", - "execution_count": 44, + "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "center (31846, 3) [[ 0.0001173 0.0001173 -0.06399978]\n", - " [-0.0001173 0.0001173 -0.06399978]\n", - " [-0.0001173 -0.0001173 -0.06399978]\n", + "center (221571, 3) [[ 4.40434422e-05 4.40434422e-05 -6.39999697e-02]\n", + " [-4.40434422e-05 4.40434422e-05 -6.39999697e-02]\n", + " [-4.40434422e-05 -4.40434422e-05 -6.39999697e-02]\n", " ...\n", - " [ 0.03213276 -0.00083309 -0.05534249]\n", - " [ 0.03213967 -0.00049989 -0.05534249]\n", - " [ 0.03214313 -0.00016664 -0.05534249]]\n", - "ds (31846, 1) [[8.64576354e-08]\n", - " [8.64576354e-08]\n", - " [8.64576354e-08]\n", + " [ 3.19445148e-02 -3.12452844e-04 -5.54567430e-02]\n", + " [ 3.19454928e-02 -1.87473619e-04 -5.54567430e-02]\n", + " [ 3.19459817e-02 -6.24915253e-05 -5.54567430e-02]]\n", + "ds (221571, 1) [[1.21882787e-08]\n", + " [1.21882787e-08]\n", + " [1.21882787e-08]\n", " ...\n", - " [1.10575255e-07]\n", - " [1.10575255e-07]\n", - " [1.10575255e-07]]\n", - "normal (31846, 3) [[ 0.00183287 0.00183287 -0.99999664]\n", - " [-0.00183287 0.00183287 -0.99999664]\n", - " [-0.00183287 -0.00183287 -0.99999664]\n", + " [1.55696068e-08]\n", + " [1.55696068e-08]\n", + " [1.55696068e-08]]\n", + "normal (221571, 3) [[ 6.88178785e-04 6.88178785e-04 -9.99999526e-01]\n", + " [-6.88178785e-04 6.88178785e-04 -9.99999526e-01]\n", + " [-6.88178785e-04 -6.88178785e-04 -9.99999526e-01]\n", " ...\n", - " [ 0.5020744 -0.01301705 -0.86472646]\n", - " [ 0.50218238 -0.00781079 -0.86472646]\n", - " [ 0.50223637 -0.00260369 -0.86472646]]\n", - "VertDisplay (127384, 3) [[ 0.00000000e+00 0.00000000e+00 -6.40000000e-02]\n", + " [ 4.99133044e-01 -4.88207568e-03 -8.66511610e-01]\n", + " [ 4.99148325e-01 -2.92927530e-03 -8.66511610e-01]\n", + " [ 4.99155965e-01 -9.76430082e-04 -8.66511610e-01]]\n", + "VertDisplay (886284, 3) [[ 0.00000000e+00 0.00000000e+00 -6.40000000e-02]\n", " [ 0.00000000e+00 0.00000000e+00 -6.40000000e-02]\n", - " [ 3.31783877e-04 0.00000000e+00 -6.39991400e-02]\n", + " [ 1.24573608e-04 0.00000000e+00 -6.39998788e-02]\n", " ...\n", - " [ 3.20000000e-02 -3.62594489e-17 -5.54256258e-02]\n", - " [ 3.22851678e-02 -3.34754062e-04 -5.52589891e-02]\n", - " [ 3.22869033e-02 -3.65845413e-17 -5.52589891e-02]]\n", - "FaceDisplay (31846, 4) [[ 0 1 3 2]\n", + " [ 3.18920555e-02 2.05145347e-17 -5.54878076e-02]\n", + " [ 3.19997551e-02 -1.25193908e-04 -5.54256258e-02]\n", + " [ 3.20000000e-02 2.05839699e-17 -5.54256258e-02]]\n", + "FaceDisplay (221571, 4) [[ 0 1 3 2]\n", " [ 4 5 7 6]\n", " [ 8 9 11 10]\n", " ...\n", - " [127372 127373 127375 127374]\n", - " [127376 127377 127379 127378]\n", - " [127380 127381 127383 127382]]\n", + " [886272 886273 886275 886274]\n", + " [886276 886277 886279 886278]\n", + " [886280 886281 886283 886282]]\n", "Beta1 () 0.0\n", "Beta2 () 0.5235987755982988\n" ] } ], "source": [ - "Frequency = 500e3 # Hz\n", + "Frequency = 2000e3 # Hz\n", "MediumSOS = 1500 # m/s - water\n", "MediumDensity=1000 # kg/m3\n", "\n", "ShortestWavelength =MediumSOS / Frequency\n", - "PPW=9\n", + "PPW=6\n", "SpatialStep =ShortestWavelength / PPW\n", "\n", "Focal=6.4e-2\n", @@ -156,7 +145,7 @@ }, { "cell_type": "code", - "execution_count": 45, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -173,18 +162,26 @@ }, { "cell_type": "code", - "execution_count": 46, + "execution_count": 5, "metadata": {}, "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/acoreas/opt/anaconda3/envs/bbmac_arm64_visco/lib/python3.9/site-packages/pyvista/jupyter/notebook.py:58: UserWarning: Failed to use notebook backend: \n", + "\n", + "No module named 'trame'\n", + "\n", + "Falling back to a static output.\n", + " warnings.warn(\n" + ] + }, { "data": { - "application/vnd.jupyter.widget-view+json": { - "model_id": "7a66771c36c74bb6baabdd358cb107fc", - "version_major": 2, - "version_minor": 0 - }, + "image/png": "iVBORw0KGgoAAAANSUhEUgAABAAAAAMACAIAAAA12IJaAABYr0lEQVR4nO3dW8w131kY9nnBJIbQkGKbs004hIOLaUNKi8UhCSROjQkHt0DbtKqEqvQiuXEvSlKUqGpuioRU0ZCoXKAGIdQ2pikGA+ZgQoKJj5jYGJsUSIIdwCkubqSYkGJ7evH62+xvz8yadZzDXr+frL/fveaZtWfP+357P89aa2Y/jOM4AAAAffiwvQ8AAADYjgIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6Mgz9j4AAHL8xsNDfPBY9lxLu3/yWNgxADt4GL19AxzAux8eNng7blQJRIZ9ik8cgANQAABs4d0pA/bDVQ79UJy1xzxLi10yqoXn+kgCaE8BAFBTaqKfqsVb9u5lwGqMwgCgIgUAQI7WiX68Wm/ief3UqgTy6oTn+QgDSKcAAFhxnFw/rOK7eWpX+5YBN+2f6nMNIEgBAPCUs6T7q8rf3BuVAath2UXC0iYlAcA1BQDQu6YZ/yeM4+4VxTZ3/kkKDod9yji+K3jSMmqA601/2Acf0DcFANCdjIz8CHl8LU0v+U2KL1wOlFEGzLarB4DeKACA+9c6d6/Y+2bvyK3LgH2XA5XUBp/mYxG4dwoA4A61y/h3mQU41PcApO5Y5TaghWt+StrVA8D9UQAAd6JF0n/kRT8l792XFU2PP7QuA1rXAIGtVdovjZ/uExO4CwoA4Kx+40nGXytNP3K6H3YHNwAtDKi1FiipH/UAcFIKAOBk6ub95036l1R5Tz/mDUBXA7aZB5ht/wwfpsB5KACAE/iNyfKe8sT9/lL/G1ve/TM+vnw50AY1QEmwSgA4PgUAcFDTpH8oy9rvPuMPKHmjr14JNK0Baq35md2UWi18pk9Y4JAUAMCBzCb9F7L/QpuVAVWu+g3HHGoeICZSMQAchwIA2F8473+UkcHvkvTfPOkY0bjxu/A2ZcAGNUBga9Nrf0sKg3EY/oiPXWBvCgBgN6t5/9GG/K9T9g2qiw3enfOeom4ZEA7bfS1QfHBk46VFJQDsRQEAbC1mvP9RapLdNCk/yCKiFm/ZGX1uWQNkbz1+DfBIJQBsTAEAbOTXHx6GlDR63+z/IOn+qh2/ASB+rwPWAOXthQuBZls+yycysAkFANDWryfetn+vtf5nyfiX1H0rr37Jb2TYZgt+UtsLg1NbVAJAUwoAoIlff3qdT6Ps33cALyl/Z69+98/IsHPVANmNkSXBZ/uMBhpQAAA1/XrBzfu3zP7vL+OfVfj+rgbIa2wxLaASACpSAAAVzOb9FzHZdlJGXpi+d5L9X9vm1p/xwe1qgMDWRuP98cGrLTH9qASAcgoAIF847x+OlPp3mPRPHaoMuPsaIK9l2rj0jCoBIJsCAMixmvoPDVb+nCv7n37P1/VhXDdu/y6c/Yx1y4CmAbuv+290ffANZQCQQQEAJFi6hf/0fWT3sf8tk/5Gz7XBu3PGUxxqKqD19QDVq4Va0wJTKgEgngIAiBL57V2Pbyiroa1H/dtl//suJWr0ft3oYt/44BarfcJbW9znZ7a97rIfZQBQhQIACIn/1t4k7b4TIHuvbXqrosW7dtOpgBPNA2zQWLKCaHXThUoACFAAAPMapf7DIe8K2rq31mq9j7crA/aqAWrtUvdigMJ1/6tbL5QBwCwFAHBr99Q/KTIvvnU/26v7Vl59kU9kZKO1QFsuBNpmEiAy4JEyALihAAA+pF3e/6jdsh9fC3Cjytt6oxpgNTg76z3yQqAWkwCpYYNKAHhCAQA0T/0fVb8vUEZ83d0PrvzN/XTzABXXCB2hBqh7lcUl7HN87kP3FADQtW1S/+FI3wlQsYcTKXyj3/i2P+GAoy0Emm2vcvOf1N+aMgCIpACATr17LvUfF76+qvBt4jjZf1dJ/417qgHCW2utETrXJEB8/HWMMgD6pACA7sym/kmS3jUqfidAXnz5jvcn733/jmuAFpMDVSYBAo0BygAgTAEAHSlP/afC7yB11/2fJfW/ebpx7QAqTrbEUwPEtNdtzK4cwu15kdMYZQD0QwEAXWiR+s+6fkOpu/Jn+7sD7fss45POm75Hty4DdqwB7mwSYHVTXrAyAPqkAIA7t1nqf2N12HtoPPbf6GXvuI6o3Zt1Rs9b1gDZFUK7XH+2vVZVUKVqKoxRBsB9UwDA3dor9Y90ouz/gOex+ht30xogJrikBjjyQqDqkwCRAanBszHKALhXCgC4QwdP/Ye4lPrjx/GfPyS/krtP/acqvom3uNg0PrJFDZBaGxQWBhWrglpj/OVdKQPg/igA4K4cP/V/dMAvBWva2wZ2LAPid2k0zB/eWqUG2HgSoGJmHx8cDlAGwD1RAMCdOEvqPzS76tdXgz2q8p7eaCqg6TxA00mAuo1VVgGlBhfGPG76XDkD3IUP2/sAgApOlP1Hkv1ne9jjbLS7lVN8D0ubqryW+Mbsp2it8Ekf/67ecXdvNdAnMwBwbudK/VukiVL/gML399PNA9SaBJhtL4mseLVAwGbzAIOpADg5BQCc1blS/6FB9l8+otmJLcuAWjVA9esBtlwI1PQGQQHbn3xlAJyUAgDO53Sp/6PqF/6WnIVTnsEC5W/0W15vuhqwzfxA9iRAdljepcwZkSXnf7pJGQCnowCAk8nO/h+/U7Yk8S15s5D9H0f27/E4OWh4a60LhY8zCRBozw6uW2KpAeBcFABwGnmpf+tkN/IdpG4B0O7iy2eP43vOOcESr+RN/15rgOMXAPEBkTGrYRlblQFwFgoAOIek7H+XBDbwVnKE4f8dk/qbp76eh9nx/fcg8wCrYXVH+pc2tZgZyAurUgDU6irjDKsB4BQUAHB08an/EQauxydrjS6OMPa/5Zkpf65LhbDBu3PeUxxnHiBjwUyjVUCz7RkttS7kjQ8rufBAGQAnpQCA4zpX6j8r5qqDU2f/G5z5aU1Vvf+mu2y/Hn11U/WrfuP3rbIKKDImMrJdAfZ8CQYclS8Cg4OKzP6rfOtTOxWP7VAv82HDM/9w9d92/afustndWvftP/6VTsMivzvsgGVtucfDe/u9X04D52UGAA4nPvWf+vhx/Oen+tBtmkdWPxHHObPV37ibzgOcZRKg+oKfyLC8VUDDMHzmOP5S3L/37a+3vm43FQBHowCAY1nN/o+Tg5Y7UfZ/zNNe9+276fUATWuApgVAfPD2BUB88L4FwKAGgIOxBAgOJGbs36dojLpLj46Z/Q+1j63pKy3vOdBDrfpw6QzUXcxT8aqYod5vrfqtdW8OzHIgOBQzAHAIJct+zqtWwpQdvEE/W6r4bp7aVa3rVrMHqo9zO6CKQ/51F1nFhOXNAySdZFMBcARmAGB/Mdn/kceh87TI/ms56ak++B9J9a+DuNkxfvC+8Lmqazo9sq/pgZkKgCMwAwB76nPg/5Hh/3aqvK23mASICdvspp+pu9S9rWf8uv8t7/hZciSp599UAOzIDADsRvZfKywjONDJSU/4s67SqSqvYtrDc4IZ2y4r4Aufccsbue61+01XJb3VPV2/YCoA9mMGAPZxB1/yVWL1RX3cOP5mYn5wnDTrOMrf3891R5rA1h1X/EeGlbzqpLAj3HH10v5vSUJgDwoA2Frnqf9wyOH/ez3Vj0re5Xe5GrgkYPcCYLaxvACICYgPa3GGA5tW25UBsDFLgGBTsv9Im2X/513zE2/L81PlSt/dfyMlNwMteYrNHGch0CPLgWBjCgDYTnz2f8ea3gomVT+/j6O90nY1QOC29FX6ybNlXXTSmlYNAFtSAMBGkrL/jT8JH66ShoeF/+1yYDEOeEjHtNmg7wa/kbrXxZ7IZi8874kKqy81AGzGNQDQXMbA/yluTjKkLxCvvvr/aCsZjm+z6wFOdDXq7KZjXgYQGVOlqw0uA1gKdkkAtGYGANo6WvY/HdSv0tsuZP8ZTvTCs4eoq6wC2n4eIzLnjfwX1/QfZusv3HibqQBoTAEADR1t0X/ThKBWMRDfg0HCPNtct3qoiz1qKbkOOGaJ//Zr/8JbM46n1jUYagBoSgEArfxG1gdYow+9LXOLzZ7rpDdgOYKznLrz/qZijvwIFexmfwkZL1YNAO0oAKCJvOy/kV0OpfWTniWFXbL7Yex+AI8Ochjxml6GG58lb3NHoIolytKRhBvVANCIAgDqy87+q3/W7XtDwIw7CJ1l/c+zEi9SnL2r0tKtlu7GAYvARpcB1BrvT/ozOOzdlpb2yvs3qwaAFtwFCCq7zv5TP7jqftAd6mNzrF0A1E1N6tr+9kobP1H8Xmf8Vtpajdm3/al7emPCtrnhUtKLvWn/PLkKVGUGAGoqyf4rOuBw8hFqm7Nk/8OGcwJHOJMb380p6SrVgy+fO9o/87ClMxzzKswDQF1mAKCO6bKfHYf/z/tR2TT1OeCKlBit36OPPwmw2RcCZA/tRzYm3fi/5zMcaHyBpAVqMAMAFbjkt4q7zP4/tjhfaT0V0HoSIKP/j2+Z5O2YPx7232aLq48a+fkjvdnCeZkBgFJL2f8uMwCn/mw8aQGw2Tlv92a9+yRAo0XqJYvOA+2tLwyI2ZQatuUZTm3POPPmAaCQGQAocqixf7a35a+/3VRAV3/E5S82u4ek7D9e3cvrqzx16wstzANAIQUA5Dta9n+sozmSe8qbz/hbXj3mI7yo6sdQntwf52LrJamvseK4vRoASigAINPRsv9+HOS8H+QwdnTeM9DiyKd9xrS0e/a9OqklplRQA0A2BQAcgs8xZyBei7VAB7x17LUjfCNVC8c55w/B05K61Kfccc4M3CUFAOQ44PD/4Q7orjnb1W38bQBncZDCrFbJVPhapru/9XhvxXAKCgBI9us+cho4zqDsWRzhr/AIxzA8+eP5xDZ3hjnOF4SFHe1esdtccvCgBoAsCgBI0yj7l/6eyB2nGzu+tPJ/AvHzcvf6G9zlYusjrLNSA0AqBQAkSMr+kz7/fHyR4YA3rlnSNDetdR5m+yn5emD/rjejBoAkCgCIZeXPSR3811b+VcH7yji9H7fVS65yk8qSv5/ZJSvtFJ7Wg/9LmXX9ktUAEM83AUOUm+w/8nMm6ePoCCOg+2pxVrN32aCrKqq/gze6s3v5t9Vu8FW1hcF53xm85df9rsZkfC1x6i51v0d52vj5shqIYAYA1k3H/o/2CXO047ljRzvVRytIluz4/buH0vTvZ8evBN7G9OCn59M8AMRQAAB37vhfp3pqm52QvNS5/PAOUvK1mx+rewCf82QAvsV5qzUNAigAYMXS0n8fM0ydbkH/s4oP+OAveByGT4p+jRsUM7sXkNV/Xzev6BefvGFu80pnn+UtJgFgjWsAICR84W/1CffOLwNoNMbZ4ktzD6Xum3ijawBiIksCqlwDsNSetxg9u2VJlasFsk9yrZMffwYKz+e/Lb2BZWYAYNGvrQ0jHerj5VAHk6fi1ZAl8SQ5+BKgzZ7iPv7MAq+i3f3+G/0JmQeAAAUAzFvN/iPdR1rAYe3+B7b7AfQj5i3pOL+OIxyJGgCWKADgQI7wkcm57J7g7H4Awx7/cKavOu88JB35vu8PS8++9MIL/zC8GUJTCgCYcRn+9yF0TPuuAjrUX8XuL+0IZ+PI3651NOFzVfFMbnxtw5J/aBIA5igA4FbS4p9D3ZbuDjKVpZfwnINdz3eso6nnpLnS9r+OvGfc4IuBj/MbnD2S+O9azlvvNNu/GgCmFACwYuPcovDp7jUx/U0f4Zu417+fDcQko4V/xMf8N5B0C6D4lxBz9o55QuAUFADwlIxrf+VMp1D913SE33vre7pXtMvpOsLv6FyO/GsqKbFMAsANBQD8nqXsv/xD8USX+p3FEc7SEY6hrnavaDX/Krx7fd6TduiYVxI36v/6oRoArikA4N6cPSs9+/Fv5lwnqvBo695jZ5tTNz3mc/3KIh05rT7yscG+FADwIf+sYHzoUJcCs6TFr2CvX+tB/pya3r+lbg/XtskLD/I7ulH9LkDHfJmPbo7NJABcKABgGJ5k/xnfdZ+k+n3xArsf+VN5Vbu0slENcK4rxTfudjjebToL/yXe6z+u1n/JJae91iDLz6kBYBgGBQDEa7FGecsOT+SMH9Hb/L5OWtqd8ZgPqPVMY+svAYi/MWj2U1TfBe6VAgCeWvxT8glR/eO521uCHmQYMqPnkx55Xufn/QObtftXqq3asTDe92qKKi6HahIABgUAJC39P9Gn3aPTHXCqvJz7aJn0jt2WdL7Lt8Ye05bD2O1kHEzS38C+K+WuD1UNAAoAuLXxjeoqBjfqYRcnPexHY73ZgIpdBeRd+nnwK4DLB62PdgHAqf9RxGt0wVUnZw8iKQDo2tLwf/bVwC1SIjXA6Tq/fpa8DH6bvP/66U6q6ZHXHSWOXxZ/Ipv9O8poCTMJQOeesfcBwJ7G3LHPwF55fbb2+Ol4wAOr4pjn/MY0QXm4an/o+I6iux/A0PHJj9H0spnNnmuvJ4JjMgNAv94VvPXnPV0N3IOjXQkQeQDj1c/7HkC7p9ilhyr/rg+1sr/Ku8peA/bxtqzk32wSgI4pACBHo1Wq4eAqudTuWe+sjxvnj2vpgJ8zF59XAxzzhGwj++t1D3LStj+M7KsCDnLGAlqX0PuegeOff9jYw7jwuQv37V1XYz+BNChvU2RAZExJfOt+thF/tHmv61xno5bWc1yRYSWD1hmbktojF6iUrFAvvCI5NWaDk1n9DFdpCTz8AlkQXTIDAJU/FOMDImNK4lv3czR1l4swdZzh/4BG2X91dbP/bTrJ663pGd7xQgI4LwUAPXrXZOnniRLHip+aZ/mY3GClwYnORrltXuyOI9b7OuP1xBnnucW8Wd5LWD2S8GyAKwHokwIAVmxwNfBe8wB5z17dca6ZPsLZOLKNT87Rfhd1LwA42qtLcqh/KeH8fuom3z/OC4EtKQDozjvT7/0fsPpBWHG1dEn8am/7Dlv6DN5MSXV0CmdZ/7Oxuq+o/AuAK15fcWM1v5+2/KxJAPqjAKBHVVKE8qcrfJbqOcq4YSUwfa4q11TkBVfc98i2LPNaXABTa9+KNvvHslf/jd6UZsNicvC8C4Ljt0I/FAD05TL8X7cGqJW8Znzctvg8a1cJFI73b1kD3FmisE0RW+W5Vjs59bU3LXpr0eH2T1flWt7sisUkAL1RANCvuoNbdQewD6JuJbCa/R/nYoDL7uMwPOv8dwncMvvfsc9At9us/9lrNVGjN58q75CRg/3lQ/uzzAbAEgUA3Np36DEj4d5guHp8+n/X7fE7Rj5RlZil4GenZ/PveXg4b6JQ/rfRYhHIQRLWQtnp/vbTKXX7Lz+kmB6mRUJM2VD+vNAJXwRGR2Yv/12a993x28Hiw6rsdUDVv0Otypk53endJlHLCC4pADI2lbdH5vrZOwZUKZYqnrSlTSWNtU5vasD04b8rI6IbZgDoyAYT9KsDrvED4buMzJ3ILqvSj3+GU6dcwl0liSyQDr5Y7lAj7rX+biv23+L81Lr8NxyQ9BDungIA5gVSqMLsqu56mOzOj6zF+al1Zo58hiseWKPUsDD7rzuSHalkGPukqiytqT4nkHEMwBIFAL341YeHYfnjp/o8wL5dHTlJjdSoRqpbBhznJNc9mEZn9Tina9jkYArz2tYF1YmS7PLh/8itb3IvILqhAKA7tYaNw0lPxXUO2R+6h8q3MrSbJ6lox0qg4mqfmz4DntNykfS+w/+zr73iiPWhFv8cR8zJTM3KU39BpzhRUNcz9j4A2MG4fEOJiuM/s8+S94zZx1b9RR1T6stscVquc4imJ3zfZOU3nx4irVikbZmi7TL8f5anSK2mjjCWcRPw8HSL7B+mzADQhV+dTOzWGo07xUfLoRarJGl62O1OS4ux+aHx7zGj84P8UVVJWEuy2OrnYYPh/9a/u+y5lOyuAgEPwa033mgVEH0wAwBPWRq2Dwznh0f6jzAPUGX3vSSdnwNOksxmGw9Pjna8ejhMRi4303Q5eNPh/9Yl4gG7atRhXreF1VTGvqv1QEmFdpCaFjZgBoAuJH1KZcx0lw+2xX/wfELZl3ec9BOu4lKTujuWPN349MPtD2PInVXYMvvPs+N6lZJ9I38dGVly0u5Vdsnos3ywfzUg8PCk742QxxeBcf/+6ZMp3dlR3tQvAsv+FrCjfbnV6aYChpRjzn51Zzwt2fLe/TfO/vMK7/L1P5GNMWHVq46S/L7WSVtqr3jeYlrCAUlbrx9+odSIe2cGgI6cYh5gs9nq8ep/Z5F0clpnt6e2wfnZYCQ7da9tfrlNp6EK/8G2PgN1j636hEAn/7ohhgKAvlRZsRquAbK3ZqjV4YnKgA2O83R1UZKSl1b9r7ckoEqWv3Q2ZlPPjHw0/nhqndvqv9xtFlNtX65nz9LAfXARMN0Zo5d5PH4kpF4TXH4AgeetEh/uZ6jUW/YxrD51ld/gBrsfSnl+c6ix/4pFe3nnO+6Y1E/27GV1eX8hSYt/YuIzIuFumAHgzv2TuXu6RQ7vLQWH28ObIgOG9JHauoPW2y8QujxXzJOmrpUqHCk8dUJQ5fgPlSDWSmQLS5qNx+y3/ztMfVeMn0vJO5ikgMIC8g1uBsq9UwBw/9p9arauATK0yBLGp/9Xt9thodvqywaqlAHnqgRqpf51s+pjnsbnPrnos+6xbfk+UL1qCtz4OFL28qfq5636McPZKQDoQuSQ/1JqEh4Jy9i0dFQZnVTcK6nzcRg+/slNxCKHCa/T/XEhJrB7dkxJ/GwPx0xhr+2V+rd+0tWt8e03je96eIiMXGqsNdQ9qzAhzjtpSfMeFXP9jBcbHv4vfx+G+6MAoBdVUpPqms7+t85T332VM41z/xvmHiZpVANUTFXzXlcjFYuTRql/07p042csf5ZxGD4j4naT25+0Kv20q3+SVuckFRgH+YcMG/A9ANyz6QUA8V8FUPd7AFY/sVrf4f7UC1ojDz7jNbY4Le1O9bPH8T1P/qQPlRTWyv5bDGPH7xLfSeS+VQaey2vgiuet5BRFhmW0ZD9cjfz3JUjcL3cBoi9jyg1kKt7/5/FjJLDjakDrAziyyIPPeI0tTst1ylCx53EYfrPlhYnZY9i1Ijcbly085kaLWLLDDpj9l4Sl9pM6VdWiWzgjBQD3bDZXnjYuZYGpNUB5NhmZ3Gc/0dnLgEYFUrvTEk4jHk6eZ2yT/WfvmDQqX2VQP9IRsv9aSs5S+WD/ap9W+8AS1wBw5wo/xassLbjeVDETKsmZYo7kgJqONW5/Wo7wK8h71Um7tFv5E9il1sBwxr6F7wCrnUTKO6tJ7SWNVbL/krVAAUf4hwmtKQC4f/EfSIVDhtebPmV58ehqEpCUJRTmLqf7qIs85pJs/nTnJE/e+Un942yX/Wdk+YUpYPmQfMmzT2OqD0NkvNfF2OafYXa6n/FccAcUAHSh3TzA0qZ3PTyUf2puWQacqxhoPZh3uhMSb5vSqMo4d3aCm6Tu8H8gsuKf0xH+MkveVJOmTfJ2mV3nOfuwpHKA83INAL14fFuffioUXg8Q3pR6PCWdR3YY08lQo6vW4l9vyZk50QlZtdk4buvsP0OLlS2RKWn1oeiS85Y3nxmp4nlLrQfCPeRl/yoB7psZAO5Z9qfpUmPG1Pnq53HMR90uA427j4LXXXFR/lp2PyF5yg/7cfePj7sl4gZj/6n/3Fpk/5E2zv7r7hs4n9nFT95rrDirU3HuCE7N9wBwt37l4WGYG7Vt91UA4RHikm8JSArLjo/sM2/GY9VlnP7mLenIJ+ewcwJV3tYzOqmSwlbP3goT/ZjGI4z9hwO2P3W1JgRSA+IfxkR+kRyJO2UGgDsX/7E0G1k+oBizNXKANnUc9xJf8RPs0uH0f5FPNF79cPnfJzz5lI3/LSx1G6nKWP7sSdjY7G+kvM9G5/MgKWyh6n+QBzx1hW+JtbL/1U5KHsY/C9wf1wBw/6Yj1pfB5rzdlxpXew6PnUceVdLBX+JT98qQUQNc/MbaV1yFT911txlfBJa6V0yHdbudfYqHlHOe0X/dXQpjWg9gz7ZvkwWWZPYxMY1+oanBedl/UgZfctiBmQG4S2YA6MLSaFZM2GxkIHgpfnWvyICksCp7HUTrM9PCc8ZxOktwPU6/9N8huFeLY86YPYjfJSPmk66WXuS90sLsfyksJkkt6S3ceXZM/I6Fp2g2sspvMLXb+BT/vO+KkM0MAF2LGVdOjYyJX+2t0WzA9V4ZOx5B/Jk5wmn5v4MzG6s1wAYq5tbZkbMxv742KRTYN/VFlXRSN/lOjckeEW93iiLT9MLB/siDidmqMKBPCgDu1jRTnM0dp5niUoo5m1MG8tFwqhqToUZmsXllwHD18Xa6SiDy7D06QiVwTLtn/yUBtRLEplls0pHkhWWoeFbzds/I/pNOdbi3yGk02T/3TQHAnZvN72My/qV0PzXXz64Qws+Y11t43+zdd7HBmbnOAE50ZgK2Se8igxtl/3UPMmnfLbP/7AHv1L0K51gKU/nIXSJz+vDW+E7gDrgGgPuX9zm9FDYmDiXGPMvqJ1bSx21S/NLup/j8SzrUrs7MtSpHXv08lwcsPXV8b/H/wKctnza5O2TMK4qsENpl/9XfvsI9ZL/3Jj1R3uu92fF0/66hkBkAulAyDzBtnN19tX22n4oBs/FJuyz1UNjPBuJnA4ayqZLrHgo7aapiNtNiQL3R2H/FXVZ3/8dxlyjUOob4rupOm5TMCURm/6vZfDggqRjIi4S7pADgbt0khTE1wLBcBsSv/g/kl9uXAXm7BPq51i79HdO/ESw1Ka9yWrY8JzH2yvvj4xul/qmb4hsj+8xOylPDqjxRzF6FJUGjv4cNigHZP51QANCRyFH/yH1X25d6jsnyq9wjqHCXyD5vXL4q+PIqrrP5h6t9pyl+uPPIl5AxGzDE9Zza4fD0Gaji+qy2UDGDTA3bbGC7RWPSAeRFFgYkJfTxndStmiqm+/EpvqkA+qQA4J5NU8bIAf6lsGljoH2255hN4T5vYlbDCndJNU7+O/vDkPvhGpNM51U7TWukHcfm23W7WXbbOvuPtGORsHHtVDjMn5HrT1vi0/SSraYC6JYCgPs3m99n3/MnYyogUAMM9b4bOCZyukvSXgfR7rRk73Vq+6b+qzHHGfuPCUs9mUfL/pNKgtlJvPJcP3AYswHh3fOmAmT/3D0FAF2IrAGG4tX/eV8IENgaE5AROd0rY8d9pZ6WyODCvU6h+ih4SWT2AoyNs//yXD9v943H/uODs89S6nRQUg/xJ0T2T+cUANytm3R8mi9GZvxJ6X4gKz14GTA8/bF3lqw3MMFSGHy917WznJmp8rTmIAP/ga2p7ZGqJOupYauR25yi1tn/arqfHZ8XKfunE74HgHuWN2hXvk5gfLLpU8fx8b+f+uSu4WPKJ1ag51VjSnBg9+N/HCYdZ/mLOuOZqfWS44NjwrJ7qJXll4z0Z08IRJ6c1TeKvK2NTtFsY172/0ef/oKFcBKfWjzMRpa8IcN5PYyTbzOBu/F/PTwMc0O2MS3xYfG7x2yKj0kdiq4ydH19e59jyji2imemVm8Z6r6VZ/RWKwPOS8jKB7Bn2ytm/xucn/DWRmP/FYuERg/zIsdh+GIJEvfLEiDu3zRhnW0ZypYDze4+GxzeFHj22U5W+0l60shObj4Vp7f+rGLps7fFzZEqfk/C9LCr3wx09ll27PmYqW35aHe32X/1AmnamD2FEtNb/HkIZP9w3ywBogvjMPyRycxy5Ax1ZOMwDH94brgoPB2/+nEemVskzX2PibvE9Hn939X/RQaHX0LMUcW/zBan5abzYe31Jv2v0UGmJmG1/kSzf+NLm5LiI98NZhtjUvZa/4q3yf7znj0clvHSppn6H7taSBmOXHoo+4cLS4C4Z//o4WF1GU/2cqCkxpJNkQFJYVX2Oogq66kq7nU6GZ8B8blg+OveVnvLS3mT2ksas488NbJ66h/YVP2EhBP0mJa8nD5+x6WfLQHijpkB4M6lfrQMicOBkbuXbIoMiAn7lIXPs6ZDy63Fn5bUV3dzWs54cqZKJhPydlntMGPf8D+l+PgW2X/SKSqfQMhL8aebnj+OS7vkDerPhtXN/rN3NPYPCgDuX8wHTOTHUtKn4+5lwDT4n83MiNxWBWcsBuIPOPvVXWqA856cpcP+xIhhztSkNiY+O6+ttVd8oh///pB6llZj8gLCbzLPn/uNv/3hobAcWn1fLcz+YzqM6SfQf2TncAcUAHQh44NntmVYzgZSB8k2KAMig2ergut9T5Tvlud/SbvfnJyDnKXVjP/Gbyz/AQy1/9hWw8KHvbo1/rmyE/2lfePFn6VwDzHn4XPn7qr59slvPOkdrPDlx/eTVCJGpvJL/2Bl//RGAUAvIgeoVlvCT7H0vLPXBy/tEvnUSUleeSo/TSv3+oyMPC2RXVWpcK7PyTh52ML49H/rlmoZvcUEr/aZlPCtPnsgqS05S+XZ/2pAxbP0jqtcv7xwyq6RZt9aZ3f5gicrkcLx8fl6IMuP2QXulYuAuXO/OLkOOO9rAZKu90292HevK4AbXeRaclfQ8eqy0fD1o0k9twsu97DwqmMun60r4+kKy+P4flIT4qT4+B6ys/9aYXmnscrZyz4hqS2r8fEPA92uhn2J7Ii75nsAuH83aeg0K51tGSZ7DRFh4fbZ5wrHRwYkhU3jMzL11W6v/5ux++q+l60xR570Aq+fd4Ni4OZERb78us/ebseYsJKB//L27Beyceofjtn+LGWfzPjB+5j4itm/fJ+uKADowk1mHJ/fR35f2LSxYnt8wE1YZNZ7nYAe/4t+p5Je5pD46m6KgdOdnFmtk/7I4C1Hu+seed1OuqqRVrP/7HT/5mHG2H/4wODOKADoyHQqYEhvmfYTaMxrnz7jbEA4Gb1J6yNNh+1Pke+mzgbEB093PN3JeVSS0IyJ65HKM9pAQN4oeHx8TGN5bVOxq4onKvuE5LVUzP7jf0ExRYLsnx64BoD794tPXwWQ8dVg5Y0Z7TFb42OGGgPYJ8p3h/3W/e97llKT9XBXLXYpSXmbJrWz7TEtMZsyglucqEYTAtVz/dWA+Ow/NeN//PlLpUbcOzMA3L+bUfDVYf6lUf/ZxipTAdOeb7YGAlZ7uAlLnROYPZjLIRVWFE0lHVvFcf3ZxKHuJRYPVz8kZaiRnVcZSs8Iq5v6L22qnv23OF3hmL1qpLx6IONPtHr2n1EVwH1TANCLcZL0h7P5yIw/vloItIc3pQbEJL5V8t3riuJGyY2ALp2H09zLE60e4WpYlb0i+6zbW61us3OgKkP+qzG1Uvz4+MKcOP6pU8PytlY5gTGNeUVUdrof6CfyJcv+6ZYCgI6spvgx+X3S/MAQHby6KdznTUDSMH+7fDdQIcR3Eu4hafZjiAsu3OssCicKaoVVH+pO2qXiOHdYlUrpXAP/eS1J2X9kKh8TphKgN64BoAvvCF4GsHpVwGxLUmNqcMzWmIDImNm9sgfvd7TNuv/rqYkTKXyv3zjfrZjsppYKeYlsfP+pkS0KgyNk/6mzBxkpfmBT4OcvkxfRATMAdGE69l/rdkCRjanBMVtjAoanP9vcDmg2ODJ+uuNhz0/d/OU4qX94a3n2Xyv1rxi2b5mUPWdSmP0nnXZj/5BBAUAvbnLl2RQ/Zv1P5FW/S8PD2WVASUBGZGDf7B62FFMaTeOHxL2WerjuKvVIkp7xMhdRvedGu2yZ+i+1FzaG1Rr1Xw3I21pYI8XsvhqTWhuEd488eNk/3LAEiF68/ckqoCrrf9rdFTS8aXVrfExecLifpbKnnchnLFzt04OmyW5M8PGz/7zh9ozg7IB2Z6li2AYPM8IuP1v/QyfMANCdm5Qxciqg5K6g08jh6vMm6WZBq1tXOw8ER8aH+4n58Fy6R9BN4+VhzPO2WAc1JJ7J08nOdOLH+1d/gyXZcGoenJ3FVkn9yydJKo76L7XXLZCShvZXnyJyUD/QibF/uFAA0IvrXHOa1Y1zDwMB0z4DYYH2pX5mjyG+w9Sw6ZMm7ZVqnPx3dtN062qfwyZ3BMrY/SAKE52MX0dJ9r996j/bGL/vkpjgFoVQYNMBs//sdD9jHkAlAJYA0ZFfeHi4ydhqrf/ZYPFPlTsCxYdV2Wtf26+DOs6tky6FX/n7e0YPJblsOKBpSVBSDyQ9dWrkXgP/s+0tWpJqg4y0PjJsHIY/LiOiG2YA6MvsSH+7FUHTxkB79qbrgHDMdT9JqeoZx7+TRvqrTH3Mzmlcd1ulPLh0Up5nV9y3yiB3OKZKprvUXjf7rxWWF5B6DuvOhBSO9N+0ZE8LzIZFzhVADxQA9Ci1DIipE8KNkcElmzJiCrPe6afmYUuC7Gqn4isKlwcZXVU0Zk0X1B0Lzxj1T92rRWNSQJWudpkQyAtLraCy0/3UzF7GD5YA0ZdfCH4j2LRlm/U/x7kjUEZ8qqUrgB+Nk0LrOj6m88hjyHB9afJ9yHj3Ty0Vts90lza1G/hPOo2FVUTe6dqgHKo+FZBXDMSELcVY/0NXzADQl+nQfnh5z2zAMIkZFsKW2lNnA/K23sSEw7LjU42T/wYOIBwZ2LfFTYFmE4izFAPlqU38eomS8f6YgIxh3ficPin732YOpOQwtpkJaZr9l0fGVAtyf3pjBoDuvO3hYUgf+99y4H+DS4GTIi/xpxv83v464Ouutj9dl0qy1tt6aj97pf6BTS0G/rfM+1cD4l9gRvzqMP9sY1KuPwQT8bzCIKMq+BNyITrzYXsfAOxjjPikGScPpy2r3c52vhoc+DwObI0JmD5L0uD6+PT/Di7pUG+CS17guHC6pv8bFpKV2R9W+7zZZeq5V1nOc5cznqTfb/zpKvnTzdgU/48x9S+kbuRqP6lbk05U/NmL2Xf1TXI1fik4KXK1h+k/K+iNJUB07fGtv+SmQNOApMaM9pitw9WnWvzVwFWuA4550hbGuFf6KP6mQEl7ZUiqAYbiZOVdV9fAXH7O6zPyqJKqiIwe4jPgQFqc0RJ/MHmRGa86vFfhiVo9JxnVTvzD8k0x5TF0xRIgevTzc6uAHhZ+LmlJagy0hzetbo2PKYnPc3NN8CWPv2lM7bNRcJUdd1T4dp+0e+H4d90MOC+pjTmS+Jjyriqek4qpf0xLUkAgOGPgPybM+h86ZAaAfgXGjCOnAlZblp5laXg+MBIfHtG/7Lj6ogKdBOJXO882Lvx3+jC1z0fxg/0xwUs7XvfQ4iwluS6iqvTWInibEe6k9tTsv+6ZaTEnsE05lJrr37S0mAdY3SVmX7hvZgDo1M8/Wf9QfjVwTCeBxtTgmK0xAUNxtnquIfCDTIBcJjQK+8+bG2nac5WR8taD3PHB2bMBGcH7lkPtBv6nLVUeRk4RrP78+MOflAXRJTMA9G46sj5GPByeTuCmafTSbMC0MRw8Gx/oKnL365iYyPBTXHYfi4uKdsonQOoeRnnSUTdtSR0TTaoTtsz7l9rbTQXEH0NSzF7npNFUQEl8ZKmT+jP0zAwA/Xrrw8NNVpdx68+Y+YHIrvLaKwYUxsd0WLdCCOSgJkDCSt73Ww9sxwRUGRQvLAmSnj01rOJrb1H5VEzuh7UEPTKVz6sEDP/TLTMAdC1msD9ppH9pNH22q0BkoH26aekpZgMCPRTGR3Z4M+1Qpc9we/iiiFoTIEO95T21JA3Sr3bVKDh199TMuPrId0ZMld5qTQgUThGktqTWBkldhcOWnkjKD4MZADr31qs5gKTZgM1uAZR9X6CYgKSwKnvtKOOAW7zGukukLvVG9ffxFiteVsPqLnpZ2lTYGLk1NTIv9W/02rcpBvKC42cPZne5/tnwPz0zA0DXrlOxcTJCH7jDz3TYeHboN9CYFJxxX6CYgJuw1cilvTJ2317GSP/sGH+Vw6ibd7SeTqkSnz22HQgoHMyuFbykfAIh46Ql7ZJdJGSM6ye90piB/8hdYioH6JAZAHr3lsl3AqQO9pes7697JUDdywDqZvOb1QaXdS9NJ0BKdtxd0yU6FcMqTgi0XgiUEV899V/atPGan4yAyFmC1DmBcEnw5ZIf+qYAgJkaYPpwNSByAU+tdT6FlUBS2LDVvX0CS1lKVrlsXPYcpzAYy85bdifZI/olW+PbyyOTDiApJi+gSklQZeB/dZfqlUBqhSD7h0eWAMGH3Cy2ma69GdceDnO7rLYs7RuID+8yPP2BV+Wy4DExPk/JKpGYbre59DnyUFcvLo8JuPx5ZCwICchYLHGJ+fRx/Me3t9eK6mezEfHIsLoTHauRh5oJiayFksqDpK0VK4GliQLonBkAGIYnkwAXhRcEzzY2uko4fmvkEH5JZr/NREEtqcd5lteVoUqdsBSzOoFQPTMuHOMvOdrU4OyAjJIv++UXDvOnxtfK/pcCDP/DYAYAHoWH/5MeLjVGhi01Dk9/huXNCSz1vPQsq8GBZ5m6lAdNi4T4p7g5n582jv9kYeh6Gn/tXIVBYe6TMeBdfWw7sCkpoS8ZIw9LmjNJ3Zq6qelUQGoBU30eIHtCADpnBgA+5B8+Sf5KLgjObgm0Z1whEN5UErnXAH/5WvbtL4Dedz6kyjt7RiftxrwDW8tz4sjR8bDyvD8cUGU2ILsxoxjYZeB/6efHH75CzgPDMCgA4NpsDTB9mBEw25LaGMgmq9z/Jy9VPdf498WON0RafZbLRM3ND9tolPiuhuVlxi1S/5iDSY1skfcHthaWBOVTAUnxhZviywDZP1xYAgS/55Jh3yyVuXxoxK8CWmpZDQs3znYyPP1RF74sOOaa4D8yjr/88BA5ej37iXr8iwFiTsg0+Fr1VxezyKG67M5rpcjVV8LEty/9iZZPYiRF5qX+5aP+kcFnH/ifPULADAA85efWvhZg2lJlfmCpMdCevald2Gon4ya1weVZLs8br/Aa6AO6KWuzO6kb32JEvDD3zR6DzwhuURRVf/mt5wEqpvvXD5faDf/DNQUA3JrWANOHqQGfN46/sHCBaeEVAoH2pICkpPyYmW6k3Rf/1O3zptqp2G2LvVZT2/Br2TL1r570x4TtkvpPGzNiSqYFrn/+knH86au3ytSMf+ln2T/cUADArZ/LvSVo5Lh+SVjJptWteZHXuxx82c+NvNd4xzI+DDZIfAOJ7FK1UHfYO6x8vD8ckFESVG+smOsPcTl9/KbVgf/HH/6UVAee9mF7HwAczh99uiwe1z6Wxqd/vjx8wThOd59tme573bIUPCuwKdxnduT1Ljc7jlftBzQ91NRdUk/RQZS8iqRdVsPCAUubAn9aq7uEIzNOQkxkdkDgtSS9M0RGrp6l2YD4+MB76eymL5u8hS79PHsMlx9k/zBlBgDmvflhZslO4bUBMXulNsZsrXKPoEvk2GxhzLjW/+ym8el9w0+RdDyFLivvby773mAaYbx69vKuqsev5sQVN8Vn+dlHlRpZEpB08CWNqy3hgKStGZtmf55tlP3DLAUALHpzm28GKL8COHt1UGT2mZqhnm5hTHYZc4RXusFdQZP6jym9YvrMy3pnnzo1xY9PlJfETx3kbc2YJchubF0M1Er3l36+bpT9wxK3AYV1j58hN0Oqsw9vgmdbrj+Rpo3xtwRNvd3nuBYQ00k4/sYR0uVZN6fiUcb9QHd5gXUzmrzeVvO/+P6z9539JTYaIM8OiwzOmJQofKUxYXld1U33rx/G7B6ehQAuzABAyM8+vRBos+H/LW8Juv36ny3z5jH9sKu8wConqpbyd/mkBTDlEwKpQ/7hXUoas8NiImvl94FN7Ub9h7VsOzL1D0TGVAKBMuBPS29gmQIAVvxsxF1By+8ZGrljdnvFgOzgEtXXvRx2mdPllU6nfa5PQuRlDxmqzAm0iCnP+1M7idmaEZz9dOXFTJXUf4NJgKVN8SWB7B/CFACw7mfnLgZYfRjTUn4FcIs5gbzI8r32kn2053qZs7LT/cjaI2YUPNzVUoLb4h6gMZvygktqjKTjzC4Sqi/4CT/MmASI/1n2D6tcAwDrLss5rsdfVx9eWgI3gZnuktQ4+xSruwyTT9+Y6wHCYeH+I59rL9OzMcYd59LLPOBrHIrnCmLyxaQnWhrWjew8L3/NCI7vISMsdQYjdVNeWh/T0jr1X9oUUwYAq8wAQJQ3TRYCtRj+P8jY/5bL5Vf7Hyf/vXE9GJy0JCby4Nu9xutXdPPD9SsNtAyJLznJoRLiwKbCLL/WkcdHVhzvD7S3mxxYLb0qpvtLPwcCXiSrgQgKAIj1pqyFQBmFwWxLUmPM1shx7oz092iXwMYoWQV0uhc7lfcxsEHSH9haPrqfnYhnBJ+i8knK7FcDGqX+Sz8//iD7h0gKAEjwpqt7AhWWAXkxgca8+KSY+LBGu2/jju8CVDhdUDEn3mxtTPnSoLCDjPfHt28zD9Bozc9s++UH2T/EUwBAmjdOviE48DAvxW9XCaxuigy4Ca6V7IaX+sQbIzLdQ13cPL1KZLrgZ3g6Zmi8Cii1w5j41d9LxamAWvMAjSLzTkKLM1B9HmDj9T+Xh7J/SKIAgGTTGmDYZF1QoD073W+xBCh7r31lv9JDDfZnaLr+ZzWy7lRAlfmBpIBaXW28BCivPGg9CZD9s+wfUikAIMcbI64HqLX6v1bGXyXXL8x3z5UoV1kLtLvL72uzxT+Ru1QcBQ/MLVSpE5YiI09sxZmNKu2pLRtMAqyu81n6eRyGPyONgXQKAMj0xqenAUoWAlVfGrS6aXVrTEBGZExXtQbUr9fSBN7mUp+raU7/MDns6XKgYe0VVZHUf2Tw7lMBhUuAknapPt4f2FRrtU9MS1KpUJLuL7Xf/Cz7hzwKAMj3hpTrAVYfRrZktJdXAhlJ+REGv1N1MrNxrfX6n5jg1qtfltrj65aYiquwKqi1eKldMXCo9T+PP8j+IZsCAIq8YfL9ADcPay3y2eYi4A1G/U+0aL7ulc17qfgWH5njrqbL1UfH6+b9SauAysOavqhTrP8ZJpl94OfLD7J/KKEAgFJviLs3aOp1w5F7pe4evzUmIDs4sp8qqfP1splpbpo3s7FZTt/0Vj9LWqz/iYlst/olr/OA8ipi+xe7uphnNabKJEB56j/I/qGYAgAqeEP09QCrD2dblhr/nXF8y+SORO0uAs5IeU8x0n+t5IBP92IfNV0FdNipgBZFToukf2lTSZFQnvoPcTl9ZFhqSSD7h3IKAKjj9cVrgUoqgSqXCjSKmd3rLKuAHh3zKucll/5rvbkffElM0i61Uvn4sOyqoHwe4IvG8bVPDxBkrAgKBFSvCsZh+NPj+OMPD4GA/0DSAjUoAKCa10fcG3T1YWRLUmOgPbwpKWZok90m9Tlexd/sFfNOl7rk6VxlTIxG4+KrBUlGSTC7oCvQVaCT1QOIfIr4mLzXGxkceU5SU/+8BT+BTUk/P/4g+4daFABQ0+uX1wKlzgAUJv3ZGX/Tgf8TKTzag5cHJe/7VYa9V8PaLYnJeO2t1ztVebF5Wf60paQwaJf6D7J/qEoBAJW9ruzeoDEBqY1LmWjdXP/ulwO1uPz35uWPc42zrse/6y74uT6S6vHV18MsbUoqEtpF1i0JChcFFab+VRb8XD+MT/3HYXixXAWqUgBAfa+bWwtUvvKn9er/6vXANpn97O19lu75E99b6l6n1m5QfDUyb2tS6rz6FKl/J9llzGpAYYo/21jesn3qfxMs+4fqFADQxOui1wLVGu+vfh1wi7VAN4PWp1BrLdDNSP82Z+B6SqFKV1Uiq5cE5XlzWHnGH4gJVCCth/xXW+IX/wR2LFz8I/uHFhQA0MprE+8LVHHlT/Xh/8hx8Sr5/fVKmy0T5d5G/cOOs1qmSj2w/ctZDWi32KlKS5VR/4x0/+Zn2T80ogCAhl5b475Ahfl93gxAlXogO/44Ci9sOOar3mw2oCRFzsiDk24KFB+QGpn3opa2xtcJeeVBUkAgOHJTUhkg+4d2FADQ1mtz1wLFBIQbv2Ac37z8LcUx/cTsmBo2u+MxE+Wpiod6UySMTzcOy090HTBMfq71hp7Xz2qCGz7CikPmmy3+iQlrOrkR3xiZysfEbzAD8JWSE2hJAQBbCJQBqw9jApIaSzatbr2JrJ4uD5O8OZwuX1y/08UkyqkvM3Wv42iR68eHVZwEWNoUf6jxFdQRXtdsY+oA/+ou8RMCVcqAwcA/bEIBABv5B2uXBKw+jAlYasyLX92aGjacarz/WsUD3uYMXNckm00LNB0vD2xKHe9PyvVjOmy0Nb6YyUj0h2A2H95auClQBhj4h20oAGA7/2B5QU6VywC2uRo4RpXlQLvUCXlPerp6JqxuZhwZkzqon3TnnBavaDVsgzonskJInRaom/qvLvi5PJT9w2YUALCpnwl+TdjNpi8ax9ffhhdNAsxWCOGUNya1TRr+P6N7uryh5B0/ad8tk+OkmYHZsFqLf8JdNUr6I8M2XvyTGiP7hy0pAGAHP5N4h9CYljMuB7re5WiJ8qozHnO8irl+ZEytRTLhrlJf1/U11uGwWiuXApsapf6tF/9I/eGAFACwj5+JWw6UUQkkhQXaV7fuuO5/ejXwxWzj7Kbrh4H3wbzDnl6mfDStpwLKY6oM9u+y+CcQUGvx0rT9S8fx71+9q1QvDEpS/0DASyQhsAcFAOzmNU8vBwqn7BUrgbz2mK0xAVV2OYK6h31JCmdrhpuh6MDPFd/Q87oqGemvOIheK48f4hYIRcaUt+etCNpgBiBml+kPsn/YiwIA9vSaiFsDrQY0vSB4aQB74+H/pkPp0zmBwic9XUmT8TEQv0teWhzYmpoxV7/hT5WY7EmALx3Hn354mA1OyvtXAzJmAJKmAmT/sCMFAOzvNXFfGBwTkBez2h6zNSagyi77qlWKhIf8M1xPCAzpKW9k/+VhJYtnItszkv4qr271eavUMzGNSal/IDh+U9LMgNQfdqcAgEP46eByoFpXAOdNAgTElAQl2e2+a+iTnvrIy/0jtciYVyPzZgBm8+xa63lSO8x7rsJioEpL3sPV1P/F4/jDC1cjjMPwVbIOOAAFABzITwdvEhrZUnJtQN5yoJjdC4NT1V0NP9t/yb4bLGe6PKzVba3IDWYAYp4l6SmSwlJnAJLa82YAWqz8yfhZ6g/HoQCAY/npheVAMS0tUv/yjD81071eG3M6J50HKPwYKBnjD8fUWk4zDUsqEUuqmkb1TExYo4VASan/5QfZPxyKAgCO6Kefngk4yIqg1WIgdc1MuZs07uEqt7u5rnf63+nu8U9armmRUH0eILWHFkl/UpI9/QOIOZ7rWyqFI/MCqtQzTZf6ZESurhSS+sMBKQDgoP7+3FRA+QKhyJaYTatbh6eTsCSnG0R/dJwZgC8YxzdPvkZ61c3nwR8bx59d66TpspmKRULqU18HrJYQFV9ddmNqYZBXFaQO//9ZOQYckgIADm22DJg+jAmI2Wu1PbxpdWtJ8M2ONwP/j1KT75vcLnJmoDC/n05QNDLtfLOpgM3WzEQO8G+55icckLRyKWOwP6alURlw0yj1hyNTAMDR/b2yGwRFtiw1huPH5a2RPVffa1/HmQEoUTFRjow88gxATExGPRPYlJflD8HMftpSOOQ/LKf+g+wfDk8BAOfw9yZfGRaTxLe7I1DGpuzI1X6uc+68/PtmSD78tlhlBmDYpFq4ntwYsvLj2d5Sw/69cXzD00uJas0AxB9VUuRqcPjvpFHSHxkWn+iHH+YN/0v94RQUAHAaf+8qhdp4RVB40+qOkRMF0w5POpp+xmPO/iRoOoJeMcMOxMdfJZy3tTDpj4zMy+zDW+OH/B/J/uEsFABwMj+19l0BR7gYYLN5gOnAf3nNEJkRlmf5jYqcmw7L3+IrDpzHxFRZPb/UQ9Idgapn/EubsvP+imVA6lTATeNXyyXgVBQAcEo/NVkRtPowsiWvPbxpdets/EmH/4fTHnzeh0HJ+vjVgLxlP/GJflJwxbolqf0IZUCgHpD6wxkpAOCsfmp5RdC0pfA64KbzABWX1D9aTb5nAzIGzqtMAlScu7i4Pviku99E9lkSEw4rHDv/knF8TaX7lq5G1q1bUnP6jF0yyoDVeQDZP5yUAgDO7e+uTQWUzAOkBq9ujQmYjT/XgPpmN/qsbuPVPqsB2TMAqdrNYxSuaJqm5rMVXd7anmEh1w90KPWH+6AAgHvwd+O+LqD6TYFi1rpsMw9w6WS8+nl4+vACh3qdr1/387Cccg1PB9RycxhD8CXcbJr2UEvFwmCbhTSFwavrgiq+iozB/piwvOH/1Z/HYfgaaQOcnwIA7sRPFl8cXGtR0OrWmIDryLFBnr2lw04FZL/7ly/1iQnIW0iTuuSpyoqgbYb/ZxsDe7VYAiT7h/ugAIC78pMRFwfHBETuGLkpMiAyJtLNUPoXjePrJqvDp0n5zdKdpPfHFvn9zZxG02H+a3k9F67MiUyLI3vL2OuPj+PjzXZLSpfyYmA1xV+NyZgKWG2X+sM9UQDAHfrJsuuDY/YKNCYFRMaUxMd74Ti+du0S0ir2mhO4Xll03ZLaQ5WwRqt9liq3S3IfGZ/07PGlS3YZENPSIu8fpP5wjxQAcLdenbgoKKkxewlQ0uKfdlaX9Zd3frQFP/EyXn5hAr0UUHH4P35KJ3seI2lTxiqgPzOOP/rwED/2P0SXBEs/S/3hXikA4M69Ovf64KTGQPvspuntGvPS5b2S7PLLEq6vMx6uuqpVOdwsFrppLO82OyActtT+ZQvj99kHMLQc+49v/9Pj+ONrOf20pSQ+8uevlRvAXVMAQBdevZZtV68HwpsiAzIiD+Umyz+sFuvpY8Lq5taBfuJndaoP/5esDsoe7A8/DP8s9YceKACgIz8xuUR4yxuDxg+ZZ88GXPffYhHO9cD/7Ftn3YuYb0bxb2qJ2a3VJQ2ol4Rl1wlJz560S+AXHd5Ut7FK3h9TD0j9oR8KAOjOTywsCpq2tFsCFL+11i7HcZwrBFokzZGRGVuzpykOtfR/2hhZEVXP+wepP3RMAQCd+vG1S4RvWr58HP/u3CLsvHU+Wy7+md5Jc2m1/Wxefhl3z3ivPEKWn6TRWqB2S30a1TB5ryipEogZ6Z/+1WUs7Fna9PjD18kBoEsKAOjajz89G/CnxnFpfuDIswEHGVAvMb1aoGQ5U/XVQRUXAjVaApS3V+T8QOox5w3qx8RkP/yacfz+J/+0pf6AAgD4vTLgImlyINxYsqkw+FAOdUFwlff9woH/QEB5z9l73cQ83qUn3Enh9QDtyoCldqk/oAAAfs+PTdYFDXHZf964frsLgi/7bjA5MPseeqh0P1vkx0PhopqSqwVuIpNWapWUGUmlS16KH17/k7EWSN4PXCgAgFs/tnyVcFJjavAwd+ubVafOsK/V/R6A4Ulvq5a+Ijevt2nMn3z66pHsdD9vuX/TK4ALy4CYywDyHt78/FIf9MDTFADAvB+d3DN0+nC2ZakxvCnpsuAqWXLgVp7j5GH47p+pz9gi18+Wl1iXR2av/8kOjowPHFjSUp/Z9nZlwGy71B+YpQAAQn70auy2MPsvuQygcCFQa+NynZDx7H9iHH8qbki+her1wJeP46uDL6fuqP9QfGlvTCclZUDhov/ITVJ/IEABAER5VeJtQ8ONq5uG9GHy+7gd0MXSLYCmD4fJq56WH4Vv9BVH5UvS/XargCouAcpb/5M9D3DTLu8HYigAgASviv4SsUBjoD28qTB4qYelBT/TrRXtuAqo1itqkcq3Xui/2U0/q5QBSfMA/6FPcyCaAgDI8aq4LwWrOyeQGnazy2bZ9vRd9UQ3BWqxqj4cFsiVswuwmCw/3Hlquj/b2zRNz763z9JDeT+QQQEAFPmRuWuFHyXdLXT11kBJVwnfsWkG2W6aYqiR38fExB9/9UsUhqzl/rPtGyz+ufws7wdKKACACn5k7avEZlsKN92ENR3gr3IXoJinGJZX/2ercsxVcvTsywNWn7HK7T4DAXXLgPBe4QLgP/KpDRRTAAA1/XD7CYGbsJIs+ebWPdmy30aPM1nRYgw++8ra6wskWixJyq5D8sb4Y3ZcHfWX9wMVKQCAJn44Yk5gqTHQPhSsp98l255eUnycpH9JxZy7+pB/0i4lxUDksv5Ae5VRf3k/0IICAGjrh1LmBFY3pd4j6BQJ98HVugwgdXXNam+R67KyDyywqcrKn8BDeT/QlAIA2MgPPT0nkD32H/OtYXeT989+p9h0AmQac9NerukK++HqReVNC5RfA5C0qcpaoGkB8PU+kYFNKACArb0y8csELpsis/9s4ae4GXKO/2/8U5e/hCWznbdYXr8a2e5J4+PbDfkPwRR/2vL4UN4PbEwBAOzmlXFzAuFNl4CMUf9aFwGHxbzJ7j5fUWudT3xwaleXK4OX/Nlx/MEnf1ElS4OS2pPG+G9a5P3AXhQAwP5+cK0SqDvwf7Q1Qkc7mIrB1S8CvuwY+Z2+Gd/2lXT5b8a1v/J+YHcKAOBYfjDu9kFJAZewQ6XaR1arDIifWFjK1L96HH9g7mun458ifHhLz5u9Fmj24Tf4qAWORAEAHNQPFF80fBMWyP5fMo4/FMwye1N3HiCjw9SrgUumGtrd6kfeDxyTAgA4gVc8zKfn8Tm77D7JvgXAZZfNyoAqd/V5JOkHjk8BAJzMK+ZqgXB+X7L4Z3o/n9l+Cu8ClHdUNy1Dg6ercuefpH6ydy+58He6Kf5KgHEYvtEnKXAqCgDgxL5/+VvGLsqz/9WnWFKxErjp4eY7AQKNtZTcJqjW8ZTfqihpU/i1SPqB81IAAHfi+5cX8Wev/6m1cKhKJZBUAwy1y4Aq63wu1VTescXf/Cd1a2SjpB+4DwoA4A79n+nLhGbjx7IJhKny2YCbpT4blAG1lviPTx98oy/9DQcsJfrTXS4P/2OfksDdUQAA9+/vXNUDu39jwHUenJcET5f7BzL+8msD6l4QXF4GZM8ABE74daOMH7h7CgCgO38n/Y6fjW4iVP1C4XaXHbe71c+Qctgt1v/I+IHeKACA3v0f0fVAlTJgnAzYt8vaC9W6BVBgr6TXHrP+J7z7o//EBx/QNwUAwK3vWysJsi8PuNnrsNl/rfH4wI6pa5OyrwCW7gPcUAAArHt5yqqhpdDpvYCOmf0/2qAGqBv/GPCf+lADWKMAAMj0tyOqgm2+F6xE9uqdwmdMXfwzJdcHyKMAAKjpf0+/wjjGJWmefRjfQ+EBxAdPD7Wk0vhzPqoA6lEAAGzkf8utDSLfpsvv+Fl4DIV5vywfYBsKAICj+F8fHoadFgttsBDoP/NxA3AMCgCAs/reyZTCnxvHaeNF3VVA/7mPD4BzUgAAAEBHPmzvAwAAALajAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAAAICOKAAAAKAjCgAAAOiIAgAAADqiAAAAgI4oAAAAoCMKAAAA6IgCAAAAOqIAAACAjigAAACgIwoAAADoiAIAAAA6ogAo8p3f+Z0PT3z6p3/63ocDAAArFABFvuRLvuT7vu/7/vpf/+sPDw8vfelL9z4cAABY8TCO497HcHpf8zVf8453vOMtb3nLR37kR+59LAAAEPKMvQ/g9L77u7/7la985Wte8xrZPwAAx2cGoMiv/dqvfd7nfd6f//N//lu/9Vv3PhYAAFinAEj28DBcztmLX/zid77znW9+85t//+///bseFAAARHERcL7v+q7vetWrXvXSl770Va961Ste8Yq9DwcAANaZAUjz8PChH8Zx+Iqv+Iqf/MmffHz43Oc+953vfOduhwUAAHEUAGmuCwAAADgdS4ASPDz88uzPAABwFmYAElyG/x85cwAAnI4ZgFjTIX+TAAAAnI4ZgFg3w/+PnDwAAM7FDAAAAHTkGXsfwGk8Dva7CxAAAKdmBiCH7B8AgJNSAAAAQEcUAAAA0BEFAAAAdEQBAAAAHVEAAABARxQAAADQEQUAAAB0RAEAAAAdUQAAAEBHFAAAANARBQAAAHREAQAAAB1RAAAAQEcUAAAA0BEFAAAAdEQBAAAAHVEAAABARxQAAADQEQUAAAB0RAEAAAAdUQAAAEBHFAAAANARBQAAAHREAQAAAB1RAAAAQEcUAAAA0BEFQI4XvfWv7X0IAACQ42Ecx72P4Tz+xl98/P8XfenHD8PwY5//V3Y9GgAASKYAiPMk9R+GYfgL33E9A6AMAADgRBQAa55O/a+3KAMAADid7q4B+J7v+Z6Hh4dXv/rVwzC8733ve8ELXvCSl7zkgx/84Ezo3/iLv5f9/4XvuMn+h2H4sc//K5e8/0Vv/WvlFwZ827d92zOf+cx/9a/+1TAMb33rWz/iIz7ila98ZWGfAABwrbsZgPe///2f8Rmf8YIXvOCVr3zl13/917/1rW99wxve8DEf8zG3ccsD/1O1pgJ+6Zd+6bM+67Ne9apXvehFL/riL/7iT/u0T/ve7/3e7N4AAGCquwJgGIZv//Zvf9nLXvZN3/RNL3/5y1//+td/zud8zlObU1L/a1XKgOc///lf+ZVf+dmf/dnf8i3f8va3v/3Zz352Xj8AADCrxwLgfe973/Oe97z3vve9P/ADP/BVX/VVv7chN/W/VlgG/OW//Jdf/vKXv/e97/2bf/NvfuM3fmPeMQAAwJIeC4C3ve1tL3zhC3/nd37nV37lV573vOd9qPV6uX+xSxmQWgO87nWve+ELX/jVX/3Vr3jFK8oPAwAAbnR3EfBv/dZvfe3Xfu2Xf/mXP/OZz/y2b/u2YZhc7FtD9sXBn/iJnzgMw8te9rIqhwEAADf6mgH4wAc+8OIXv/jd7373a1/72r/6V//q//wd/9Ov/vf/xbM/+pnDUC31v5G6IuiHf/iHX/KSl7znPe951rOe1eJ4AADo3DP2PoBNffM3f/Ob3vSmN77xjX/gb33zf/2x7/uOD4zf/lNv+Ws/8sZ2z/iY9D+WAY//DZcBb3vb2z7pkz5J9g8AQCN9zQB8SI2LfVP51jAAAI7gzguAX/6GL/vMT37W8JyPGf7bvzUM+6T+16ZlwB//if/mt9/1W//yH737Hf+D7/wCAKC5ey4APpT9f8K/OXzgg8NHf+Tw4R/+oQ17pP7XLmXA+LsfHD/4wff/y3/9vn/6HjUAAAAbuOdrAD7zk581POcPDcPD8AeeOTyWOXun/o8uFwZ88P97/8MzPnwYho967sfufVAAAHThfm8D+rKvG577nOHf+Mjh8SY/H/jA8C/+5d7H9JTffe9vj8MwPAzP+IPP/Ig/9FF/8HM+8Qv/l/9y74MCAODO3WkB8LKvG57zMcMwDB/4wPD+Dwy/+/7hdz+w9zHd+uDvfuCD//r9v/v//vYwDsM4jOPwUc/9WDUAAABN3eESoF/+hi/7zC9+/ocevO93hmEYfvNf7Hg8S377Xb/1+MNHPWn5fc/66N/3rI/+3L/0VS4GAACgkXu8CPh//K+GcRze9ZvDMPzyr/0/wzB85t/++3sfU8jn/qWvGobhoz/7E4Zh+KjnfuxHfMxH/cQX/nc7HxMAAHfqHgsAAABgwZ1eAwAAAMxRAAAAQEcUAAAA0BEFAAAAdEQBAAAAHVEAAABARxQAAADQEQUAAAB0RAEAAAAdUQAAAEBHFAAAANARBQAAAHREAQAAAB1RAAAAQEcUAAAA0BEFAAAAdEQBAAAAHVEAAABARxQAAADQEQUAAAB0RAEAAAAdUQAAAEBHFAAAANARBQAAAHREAQAAAB1RAAAAQEcUAAAA0BEFAAAAdEQBAAAAHVEAAABARxQAAADQEQUAAAB0RAEAAAAdUQAAAEBHFAAAANARBQAAAHREAQAAAB1RAAAAQEcUAAAA0BEFAAAAdEQBAAAAHVEAAABARxQAAADQEQUAAAB0RAEAAAAdUQAAAEBHFAAAANARBQAAAHREAQAAAB1RAAAAQEcUAAAA0BEFAAAAdEQBAAAAHVEAAABARxQAAADQEQUAAAB0RAEAAAAd+f8BYD9VIVqmVCEAAAAASUVORK5CYII=", "text/plain": [ - "Widget(value='