From 656dfb786dba6d7273920185575b7bc7ebefb3d1 Mon Sep 17 00:00:00 2001 From: Lalit Rajendran Date: Fri, 23 Oct 2020 11:27:13 -0400 Subject: [PATCH] initial commit --- ._README.md | Bin 0 -> 4096 bytes ...ityIntegrationUncertaintyQuantification.py | 507 +++++++++++++ LICENSE.txt | 674 ++++++++++++++++++ README.md | 8 + ...onUncertaintyQuantification.cpython-38.pyc | Bin 0 -> 10075 bytes loadmat_functions.py | 35 + sample-data.mat | Bin 0 -> 15560 bytes sample-result.mat | Bin 0 -> 12816 bytes sample-result.png | Bin 0 -> 74536 bytes sample_script.py | 113 +++ 10 files changed, 1337 insertions(+) create mode 100755 ._README.md create mode 100755 DensityIntegrationUncertaintyQuantification.py create mode 100755 LICENSE.txt create mode 100755 README.md create mode 100755 __pycache__/DensityIntegrationUncertaintyQuantification.cpython-38.pyc create mode 100755 loadmat_functions.py create mode 100755 sample-data.mat create mode 100755 sample-result.mat create mode 100755 sample-result.png create mode 100755 sample_script.py diff --git a/._README.md b/._README.md new file mode 100755 index 0000000000000000000000000000000000000000..3ac70716f307be43d2708e79dd9a2a6e27bbc341 GIT binary patch literal 4096 zcmZQz6=P>$Vqox1Ojhs@R)|o50+1L3ClDJkFz{^v(m+1nBL)UWIUt(=a103v;xPZC z1<}E<0H|C5O$#HC4;7b6&d=3LEGWoH)yqjNE-5WeO-V^CNmULA2I=wmIw>AR-#jEP z4WdWEXb6mkz-S1JhQMeDjE2By2#kinXb6mkz-S1JhQMeD;0ggyXA^|MKrSRBvsj@h zwK%`DC^=OjEx#yRAv3QeHLoNyKQA#Sr&1v&HLXM;DJL;68`u|y>Kf7%s{i3$kztVg G{~rMG1S?+v literal 0 HcmV?d00001 diff --git a/DensityIntegrationUncertaintyQuantification.py b/DensityIntegrationUncertaintyQuantification.py new file mode 100755 index 0000000..c2aea57 --- /dev/null +++ b/DensityIntegrationUncertaintyQuantification.py @@ -0,0 +1,507 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Apr 3 15:51:35 2017 + +@author: jiachengzhang +""" + +import numpy as np +import sys +import scipy.sparse as scysparse +import scipy.sparse.linalg as splinalg +import scipy.linalg as linalg + + +def Density_integration_Poisson_uncertainty(Xn,Yn,fluid_mask,grad_x,grad_y,dirichlet_label,dirichlet_value, + uncertainty_quantification=True, sigma_grad_x=None, sigma_grad_y=None, sigma_dirichlet=None): + # Evaluate the density field from the gradient fields by solving the Poisson equation. + # The uncertainty of the density field can be quantified. + """ + Inputs: + Xn,Yn: 2d array of mesh grid. + fluid_mask: 2d array of binary mask of flow field. Boundary points are considered in the flow (mask should be True) + grad_x, grad_y: 2d array of gradient field. + dirichlet_label: 2d array of binary mask indicating the Dirichlet BC locations. Ohterwise Neumann. At least one point should be dirichlet. + dirichlet_value: 2d array of diriclet BC values. + uncertainty_quantification: True to perform the uncertainty quantification, False only perform integration. + sigma_grad_x, sigma_grad_y, sigma_dirichlet: the uncertainty given as std of the input fields. 2d array of fields. + Returns: + Pn: the integrated density field. + sigma_Pn: the uncertainty (std) of the integrated density field. + """ + + Ny, Nx = np.shape(Xn) + dx = Xn[1,1] - Xn[0,0] + dy = Yn[1,1] - Yn[0,0] + invdx = 1.0/dx + invdy = 1.0/dy + invdx2 = invdx**2 + invdy2 = invdy**2 + + fluid_mask_ex = np.zeros((Ny+2,Nx+2)).astype('bool') + fluid_mask_ex[1:-1,1:-1] = fluid_mask + + grad_x_ex = np.zeros((Ny+2,Nx+2)) + grad_x_ex[1:-1,1:-1] = grad_x + grad_y_ex = np.zeros((Ny+2,Nx+2)) + grad_y_ex[1:-1,1:-1] = grad_y + dirichlet_label_ex = np.zeros((Ny+2,Nx+2)).astype('bool') + dirichlet_label_ex[1:-1,1:-1] = dirichlet_label + dirichlet_value_ex = np.zeros((Ny+2,Nx+2)) + dirichlet_value_ex[1:-1,1:-1] = dirichlet_value + + sigma_grad_x_ex = np.zeros((Ny+2,Nx+2)) + sigma_grad_y_ex = np.zeros((Ny+2,Nx+2)) + sigma_dirichlet_ex = np.zeros((Ny+2,Nx+2)) + if uncertainty_quantification==True: + sigma_grad_x_ex[1:-1,1:-1] = sigma_grad_x + sigma_grad_y_ex[1:-1,1:-1] = sigma_grad_y + sigma_dirichlet_ex[1:-1,1:-1] = sigma_dirichlet + + # Generate the linear operator. + j,i = np.where(fluid_mask_ex==True) + Npts = len(j) + fluid_index = -np.ones(fluid_mask_ex.shape).astype('int') + fluid_index[j,i] = range(Npts) + iC = fluid_index[j,i] + iC_label = dirichlet_label_ex[j,i] + iE = fluid_index[j,i+1] + iW = fluid_index[j,i-1] + iN = fluid_index[j+1,i] + iS = fluid_index[j-1,i] + + LaplacianOperator = scysparse.csr_matrix((Npts,Npts),dtype=np.float) + # RHS = np.zeros(Npts) + # Var_RHS = np.zeros(Npts) # variance + + # Also generate the linear operator which maps from the grad_x, grad_y, dirichlet val to the rhs. + Map_grad_x = scysparse.csr_matrix((Npts,Npts),dtype=np.float) + Map_grad_y = scysparse.csr_matrix((Npts,Npts),dtype=np.float) + Map_dirichlet_val = scysparse.csr_matrix((Npts,Npts),dtype=np.float) + + # First, construct the linear operator and RHS as if all they are all Nuemanan Bc + + # if the east and west nodes are inside domain + loc = (iE!=-1)*(iW!=-1) + LaplacianOperator[iC[loc],iC[loc]] += -2.0*invdx2 + LaplacianOperator[iC[loc],iE[loc]] += +1.0*invdx2 + LaplacianOperator[iC[loc],iW[loc]] += +1.0*invdx2 + # RHS[iC[loc]] += (grad_x_ex[j[loc],i[loc]+1] - grad_x_ex[j[loc],i[loc]-1])*(invdx*0.5) + # Var_RHS[iC[loc]] += (invdx*0.5)**2 * (sigma_grad_x_ex[j[loc],i[loc]+1]**2 + sigma_grad_x_ex[j[loc],i[loc]-1]**2) + Map_grad_x[iC[loc],iE[loc]] += invdx*0.5 + Map_grad_x[iC[loc],iW[loc]] += -invdx*0.5 + + # if the east node is ouside domian + loc = (iE==-1) + LaplacianOperator[iC[loc],iC[loc]] += -2.0*invdx2 + LaplacianOperator[iC[loc],iW[loc]] += +2.0*invdx2 + # RHS[iC[loc]] += -(grad_x_ex[j[loc],i[loc]] + grad_x_ex[j[loc],i[loc]-1])*invdx + # Var_RHS[iC[loc]] += invdx**2 * (sigma_grad_x_ex[j[loc],i[loc]]**2 + sigma_grad_x_ex[j[loc],i[loc]-1]**2) + Map_grad_x[iC[loc],iC[loc]] += -invdx + Map_grad_x[iC[loc],iW[loc]] += -invdx + + # if the west node is ouside domian + loc = (iW==-1) + LaplacianOperator[iC[loc],iC[loc]] += -2.0*invdx2 + LaplacianOperator[iC[loc],iE[loc]] += +2.0*invdx2 + # RHS[iC[loc]] += (grad_x_ex[j[loc],i[loc]] + grad_x_ex[j[loc],i[loc]+1])*invdx + # Var_RHS[iC[loc]] += invdx**2 * (sigma_grad_x_ex[j[loc],i[loc]]**2 + sigma_grad_x_ex[j[loc],i[loc]+1]**2) + Map_grad_x[iC[loc],iC[loc]] += invdx + Map_grad_x[iC[loc],iE[loc]] += invdx + + # if the north and south nodes are inside domain + loc = (iN!=-1)*(iS!=-1) + LaplacianOperator[iC[loc],iC[loc]] += -2.0*invdy2 + LaplacianOperator[iC[loc],iN[loc]] += +1.0*invdy2 + LaplacianOperator[iC[loc],iS[loc]] += +1.0*invdy2 + # RHS[iC[loc]] += (grad_y_ex[j[loc]+1,i[loc]] - grad_y_ex[j[loc]-1,i[loc]])*(invdy*0.5) + # Var_RHS[iC[loc]] += (invdy*0.5)**2 * (sigma_grad_y_ex[j[loc]+1,i[loc]]**2 + sigma_grad_y_ex[j[loc]-1,i[loc]]**2) + Map_grad_y[iC[loc],iN[loc]] += invdy*0.5 + Map_grad_y[iC[loc],iS[loc]] += -invdy*0.5 + + # if the north node is ouside domian + loc = (iN==-1) + LaplacianOperator[iC[loc],iC[loc]] += -2.0*invdy2 + LaplacianOperator[iC[loc],iS[loc]] += +2.0*invdy2 + # RHS[iC[loc]] += -(grad_y_ex[j[loc],i[loc]] + grad_y_ex[j[loc]-1,i[loc]])*invdy + # Var_RHS[iC[loc]] += invdy**2 * (sigma_grad_y_ex[j[loc],i[loc]]**2 + sigma_grad_y_ex[j[loc]-1,i[loc]]**2) + Map_grad_y[iC[loc],iC[loc]] += -invdy + Map_grad_y[iC[loc],iS[loc]] += -invdy + + # if the south node is ouside domian + loc = (iS==-1) + LaplacianOperator[iC[loc],iC[loc]] += -2.0*invdy2 + LaplacianOperator[iC[loc],iN[loc]] += +2.0*invdy2 + # RHS[iC[loc]] += (grad_y_ex[j[loc],i[loc]] + grad_y_ex[j[loc]+1,i[loc]])*invdy + # Var_RHS[iC[loc]] += invdy**2 * (sigma_grad_y_ex[j[loc],i[loc]]**2 + sigma_grad_y_ex[j[loc]+1,i[loc]]**2) + Map_grad_y[iC[loc],iC[loc]] += invdy + Map_grad_y[iC[loc],iN[loc]] += invdy + + # Then change the boundary conidtion at locatiosn of Dirichlet. + loc = (iC_label==True) + LaplacianOperator[iC[loc],:] = 0.0 + LaplacianOperator[iC[loc],iC[loc]] = 1.0*invdx2 + # RHS[iC[loc]] = dirichlet_value_ex[j[loc],i[loc]] * invdx2 + # Var_RHS[iC[loc]] = sigma_dirichlet_ex[j[loc],i[loc]]**2 * invdx2**2 + Map_grad_x[iC[loc],:] = 0.0 + Map_grad_y[iC[loc],:] = 0.0 + Map_dirichlet_val[iC[loc],iC[loc]] = 1.0*invdx2 + + LaplacianOperator.eliminate_zeros() + Map_grad_x.eliminate_zeros() + Map_grad_y.eliminate_zeros() + Map_dirichlet_val.eliminate_zeros() + + # Solve for the field. + grad_x_vect = grad_x_ex[j,i] + grad_y_vect = grad_y_ex[j,i] + dirichlet_val_vect = dirichlet_value_ex[j,i] + RHS = Map_grad_x.dot(grad_x_vect) + Map_grad_y.dot(grad_y_vect) + Map_dirichlet_val.dot(dirichlet_val_vect) + # Solve the linear system + Pn_ex = np.zeros((Ny+2,Nx+2)) + p_vect = splinalg.spsolve(LaplacianOperator, RHS) + Pn_ex[j,i] = p_vect + + # Uncertainty propagation + sigma_Pn_ex = np.zeros((Ny+2,Nx+2)) + if uncertainty_quantification==True: + # Propagate to get the covariance matrix for RHS + Cov_grad_x = scysparse.diags(sigma_grad_x_ex[j,i]**2,shape=(Npts,Npts),format='csr') + Cov_grad_y = scysparse.diags(sigma_grad_y_ex[j,i]**2,shape=(Npts,Npts),format='csr') + Cov_dirichlet_val = scysparse.diags(sigma_dirichlet_ex[j,i]**2,shape=(Npts,Npts),format='csr') + Cov_RHS = Map_grad_x*Cov_grad_x*Map_grad_x.transpose() + Map_grad_y*Cov_grad_y*Map_grad_y.transpose() + \ + Map_dirichlet_val*Cov_dirichlet_val*Map_dirichlet_val.transpose() + + Laplacian_inv = linalg.inv(LaplacianOperator.A) + Cov_p = np.matmul(np.matmul(Laplacian_inv, Cov_RHS.A), Laplacian_inv.T) + Var_p_vect = np.diag(Cov_p) + + sigma_Pn_ex[j,i] = Var_p_vect**0.5 + + return Pn_ex[1:-1,1:-1], sigma_Pn_ex[1:-1,1:-1] + + + +def Density_integration_WLS_uncertainty(Xn,Yn,fluid_mask,grad_x,grad_y,dirichlet_label,dirichlet_value, + uncertainty_quantification=True, sigma_grad_x=None, sigma_grad_y=None, sigma_dirichlet=None): + # Evaluate the density field from the gradient fields by solving the WLS system. + # The uncertainty of the density field is also quantified and returned. + """ + Inputs: + Xn,Yn: 2d array of mesh grid. + fluid_mask: 2d array of binary mask of flow field. Boundary points are considered in the flow (mask should be True) + grad_x, grad_y: 2d array of gradient field. + dirichlet_label: 2d array of binary mask indicating the Dirichlet BC locations. Ohterwise Neumann. At least one point should be dirichlet. + dirichlet_value: 2d array of diriclet BC values. + uncertainty_quantification: True to perform the uncertainty quantification, False only perform integration. + sigma_grad_x, sigma_grad_y, sigma_dirichlet: the uncertainty given as std of the input fields. 2d array of fields. + Returns: + Pn: the integrated density field. + sigma_Pn: the uncertainty (std) of the integrated density field. + """ + + Ny, Nx = np.shape(Xn) + dx = Xn[1,1] - Xn[0,0] + dy = Yn[1,1] - Yn[0,0] + invdx = 1.0/dx + invdy = 1.0/dy + invdx2 = invdx**2 + invdy2 = invdy**2 + + fluid_mask_ex = np.zeros((Ny+2,Nx+2)).astype('bool') + fluid_mask_ex[1:-1,1:-1] = fluid_mask + + dirichlet_label_ex = np.zeros((Ny+2,Nx+2)).astype('bool') + dirichlet_label_ex[1:-1,1:-1] = dirichlet_label + dirichlet_value_ex = np.zeros((Ny+2,Nx+2)) + dirichlet_value_ex[1:-1,1:-1] = dirichlet_value + + grad_x_ex = np.zeros((Ny+2,Nx+2)) + grad_x_ex[1:-1,1:-1] = grad_x + grad_y_ex = np.zeros((Ny+2,Nx+2)) + grad_y_ex[1:-1,1:-1] = grad_y + + sigma_grad_x_ex = np.zeros((Ny+2,Nx+2)) + sigma_grad_y_ex = np.zeros((Ny+2,Nx+2)) + sigma_dirichlet_ex = np.zeros((Ny+2,Nx+2)) + sigma_grad_x_ex[1:-1,1:-1] = sigma_grad_x + sigma_grad_y_ex[1:-1,1:-1] = sigma_grad_y + sigma_dirichlet_ex[1:-1,1:-1] = sigma_dirichlet + + # Generate the index for mapping the pressure and pressure gradients. + j,i = np.where(fluid_mask_ex==True) + Npts = len(j) + fluid_index = -np.ones(fluid_mask_ex.shape).astype('int') + fluid_index[j,i] = range(Npts) + + # Generate the mask for the gradients + fluid_mask_Gx_ex = np.logical_and(fluid_mask_ex[:,1:],fluid_mask_ex[:,:-1]) + fluid_mask_Gy_ex = np.logical_and(fluid_mask_ex[1:,:],fluid_mask_ex[:-1,:]) + + # Generate the linear operator and the mapping matrix for generating the rhs. + # For Gx + jx,ix = np.where(fluid_mask_Gx_ex==True) + Npts_x = len(jx) + iC = fluid_index[jx,ix] + iE = fluid_index[jx,ix+1] + Operator_Gx = scysparse.csr_matrix((Npts_x,Npts),dtype=np.float) + Map_Gx = scysparse.csr_matrix((Npts_x,Npts),dtype=np.float) + Operator_Gx[range(Npts_x),iC] += -invdx + Operator_Gx[range(Npts_x),iE] += invdx + Map_Gx[range(Npts_x),iC] += 0.5 + Map_Gx[range(Npts_x),iE] += 0.5 + + # For Gy + jy,iy = np.where(fluid_mask_Gy_ex==True) + Npts_y = len(jy) + iC = fluid_index[jy,iy] + iN = fluid_index[jy+1,iy] + Operator_Gy = scysparse.csr_matrix((Npts_y,Npts),dtype=np.float) + Map_Gy = scysparse.csr_matrix((Npts_y,Npts),dtype=np.float) + Operator_Gy[range(Npts_y),iC] += -invdy + Operator_Gy[range(Npts_y),iN] += invdy + Map_Gy[range(Npts_y),iC] += 0.5 + Map_Gy[range(Npts_y),iN] += 0.5 + + # For Dirichlet BC + j_d, i_d = np.where(dirichlet_label_ex==True) + Npts_d = len(j_d) + iC = fluid_index[j_d,i_d] + Operator_d = scysparse.csr_matrix((Npts_d,Npts),dtype=np.float) + Map_d = scysparse.eye(Npts_d,Npts_d,format='csr') * invdx + Operator_d[range(Npts_d),iC] += 1.0*invdx + # dirichlet value vector and cov + dirichlet_vect = dirichlet_value_ex[j_d,i_d] + dirichlet_sigma_vect = sigma_dirichlet_ex[j_d,i_d] + cov_dirichlet = scysparse.diags(dirichlet_sigma_vect**2, format='csr') + # Generate the vector and cov for pgrad. + pgrad_x_vect = grad_x_ex[j,i] + pgrad_y_vect = grad_y_ex[j,i] + pgrad_vect = np.concatenate((pgrad_x_vect,pgrad_y_vect)) + cov_pgrad_x = sigma_grad_x_ex[j,i]**2 + cov_pgrad_y = sigma_grad_y_ex[j,i]**2 + cov_pgrad_vect = np.concatenate((cov_pgrad_x, cov_pgrad_y)) + cov_pgrad = scysparse.diags(cov_pgrad_vect, format='csr') + + # Construct the full operator. + Operator_GLS = scysparse.bmat([[Operator_Gx],[Operator_Gy],[Operator_d]]) + + # Construct the full mapping matrics and get the rhs. + Map_pgrad = scysparse.bmat([[Map_Gx, None],[None, Map_Gy],[scysparse.csr_matrix((Npts_d,Npts),dtype=np.float), None]]) + Map_dirichlet = scysparse.bmat([[scysparse.csr_matrix((Npts_x+Npts_y,Npts_d),dtype=np.float)],[Map_d]]) + rhs = Map_pgrad.dot(pgrad_vect) + Map_dirichlet.dot(dirichlet_vect) + + # Evaluate the covriance matrix for the rhs + cov_rhs = Map_pgrad * cov_pgrad * Map_pgrad.transpose() + Map_dirichlet * cov_dirichlet * Map_dirichlet.transpose() + + # Solve for the WLS solution + weights_vect = cov_rhs.diagonal()**(-1) + weights_matrix = scysparse.diags(weights_vect,format='csr') + # Operator_WLS = weights_matrix * Operator_GLS + # rhs_WLS = weights_matrix.dot(rhs) + sys_LHS = Operator_GLS.transpose() * weights_matrix * Operator_GLS + sys_rhs = (Operator_GLS.transpose() * weights_matrix).dot(rhs) + + # Get the solution from lsqr + # p_vect_wls = splinalg.lsqr(Operator_WLS,rhs_WLS)[0] + # Pn_WLS = np.zeros(fluid_mask_ex.shape) + # Pn_WLS[j,i] = p_vect_wls + + # Solve for the WLS solution + p_vect_wls = splinalg.spsolve(sys_LHS, sys_rhs) + Pn_WLS_ex = np.zeros(fluid_mask_ex.shape) + Pn_WLS_ex[j,i] = p_vect_wls + + # Perform the uncertainty propagation + sigma_Pn_ex = np.zeros((Ny+2,Nx+2)) + if uncertainty_quantification == True: + cov_sys_rhs = (Operator_GLS.transpose() * weights_matrix) * cov_rhs * (Operator_GLS.transpose() * weights_matrix).transpose() + sys_LHS_inv = linalg.inv(sys_LHS.A) + Cov_p = np.matmul(np.matmul(sys_LHS_inv, cov_sys_rhs.A), sys_LHS_inv.T) + Var_p_vect = np.diag(Cov_p) + sigma_Pn_ex[j,i] = Var_p_vect**0.5 + + return Pn_WLS_ex[1:-1,1:-1], sigma_Pn_ex[1:-1,1:-1] + + + +def Density_integration_WLS_uncertainty_weighted_average(Xn,Yn,fluid_mask,grad_x,grad_y,dirichlet_label,dirichlet_value, + uncertainty_quantification=True, sigma_grad_x=None, sigma_grad_y=None, sigma_dirichlet=None): + # Evaluate the density field from the gradient fields by solving the WLS system. + # The uncertainty of the density field is also quantified and returned. + # The gradient interpolation (from grid points to staggered location) is done by a weighted average approach + # which minimizes the sum of squared bias error and random error. + """ + Inputs: + Xn,Yn: 2d array of mesh grid. + fluid_mask: 2d array of binary mask of flow field. Boundary points are considered in the flow (mask should be True) + grad_x, grad_y: 2d array of gradient field. + dirichlet_label: 2d array of binary mask indicating the Dirichlet BC locations. Ohterwise Neumann. At least one point should be dirichlet. + dirichlet_value: 2d array of diriclet BC values. + uncertainty_quantification: True to perform the uncertainty quantification, False only perform integration. + sigma_grad_x, sigma_grad_y, sigma_dirichlet: the uncertainty given as std of the input fields. 2d array of fields. + Returns: + Pn: the integrated density field. + sigma_Pn: the uncertainty (std) of the integrated density field. + """ + + Ny, Nx = np.shape(Xn) + dx = Xn[1,1] - Xn[0,0] + dy = Yn[1,1] - Yn[0,0] + invdx = 1.0/dx + invdy = 1.0/dy + invdx2 = invdx**2 + invdy2 = invdy**2 + + fluid_mask_ex = np.zeros((Ny+2,Nx+2)).astype('bool') + fluid_mask_ex[1:-1,1:-1] = fluid_mask + + dirichlet_label_ex = np.zeros((Ny+2,Nx+2)).astype('bool') + dirichlet_label_ex[1:-1,1:-1] = dirichlet_label + dirichlet_value_ex = np.zeros((Ny+2,Nx+2)) + dirichlet_value_ex[1:-1,1:-1] = dirichlet_value + + grad_x_ex = np.zeros((Ny+2,Nx+2)) + grad_x_ex[1:-1,1:-1] = grad_x + grad_y_ex = np.zeros((Ny+2,Nx+2)) + grad_y_ex[1:-1,1:-1] = grad_y + + sigma_grad_x_ex = np.zeros((Ny+2,Nx+2)) + sigma_grad_y_ex = np.zeros((Ny+2,Nx+2)) + sigma_dirichlet_ex = np.zeros((Ny+2,Nx+2)) + sigma_grad_x_ex[1:-1,1:-1] = sigma_grad_x + sigma_grad_y_ex[1:-1,1:-1] = sigma_grad_y + sigma_dirichlet_ex[1:-1,1:-1] = sigma_dirichlet + + # Generate the index for mapping the pressure and pressure gradients. + j,i = np.where(fluid_mask_ex==True) + Npts = len(j) + fluid_index = -np.ones(fluid_mask_ex.shape).astype('int') + fluid_index[j,i] = range(Npts) + + # Generate the mask for the gradients + fluid_mask_Gx_ex = np.logical_and(fluid_mask_ex[:,1:],fluid_mask_ex[:,:-1]) + fluid_mask_Gy_ex = np.logical_and(fluid_mask_ex[1:,:],fluid_mask_ex[:-1,:]) + + # Generate the linear operator and the mapping matrix for generating the rhs. + # For Gx + jx,ix = np.where(fluid_mask_Gx_ex==True) + Npts_x = len(jx) + iC = fluid_index[jx,ix] + iE = fluid_index[jx,ix+1] + Operator_Gx = scysparse.csr_matrix((Npts_x,Npts),dtype=np.float) + Operator_Gx[range(Npts_x),iC] += -invdx + Operator_Gx[range(Npts_x),iE] += invdx + # The Mapping Gx that maps grided gradients to staggered location is by weighted average. + Map_Gx = scysparse.csr_matrix((Npts_x,Npts),dtype=np.float) + V_C = grad_x_ex[jx,ix] + V_E = grad_x_ex[jx,ix+1] + sigma_C = sigma_grad_x_ex[jx,ix] + sigma_E = sigma_grad_x_ex[jx,ix+1] + weight_C = ((V_C-V_E)**2 + 2*sigma_E**2) / (2*(V_C-V_E)**2 + 2*sigma_C**2 + 2*sigma_E**2) + weight_E = 1.0 - weight_C + Map_Gx[range(Npts_x),iC] += weight_C + Map_Gx[range(Npts_x),iE] += weight_E + + # For Gy + jy,iy = np.where(fluid_mask_Gy_ex==True) + Npts_y = len(jy) + iC = fluid_index[jy,iy] + iN = fluid_index[jy+1,iy] + Operator_Gy = scysparse.csr_matrix((Npts_y,Npts),dtype=np.float) + Operator_Gy[range(Npts_y),iC] += -invdy + Operator_Gy[range(Npts_y),iN] += invdy + # The Mapping Gy that maps grided gradients to staggered location is by weighted average. + Map_Gy = scysparse.csr_matrix((Npts_y,Npts),dtype=np.float) + V_C = grad_y_ex[jy,iy] + V_N = grad_y_ex[jy+1,iy] + sigma_C = sigma_grad_y_ex[jy,iy] + sigma_N = sigma_grad_y_ex[jy+1,iy] + weight_C = ((V_C-V_N)**2 + 2*sigma_N**2) / (2*(V_C-V_N)**2 + 2*sigma_C**2 + 2*sigma_N**2) + weight_N = 1.0 - weight_C + Map_Gy[range(Npts_y),iC] += weight_C + Map_Gy[range(Npts_y),iN] += weight_N + + # For Dirichlet BC + j_d, i_d = np.where(dirichlet_label_ex==True) + Npts_d = len(j_d) + iC = fluid_index[j_d,i_d] + Operator_d = scysparse.csr_matrix((Npts_d,Npts),dtype=np.float) + Map_d = scysparse.eye(Npts_d,Npts_d,format='csr') * invdx + Operator_d[range(Npts_d),iC] += 1.0*invdx + # dirichlet value vector and cov + dirichlet_vect = dirichlet_value_ex[j_d,i_d] + dirichlet_sigma_vect = sigma_dirichlet_ex[j_d,i_d] + cov_dirichlet = scysparse.diags(dirichlet_sigma_vect**2, format='csr') + # Generate the vector and cov for pgrad. + pgrad_x_vect = grad_x_ex[j,i] + pgrad_y_vect = grad_y_ex[j,i] + pgrad_vect = np.concatenate((pgrad_x_vect,pgrad_y_vect)) + cov_pgrad_x = sigma_grad_x_ex[j,i]**2 + cov_pgrad_y = sigma_grad_y_ex[j,i]**2 + cov_pgrad_vect = np.concatenate((cov_pgrad_x, cov_pgrad_y)) + cov_pgrad = scysparse.diags(cov_pgrad_vect, format='csr') + + # Construct the full operator. + Operator_GLS = scysparse.bmat([[Operator_Gx],[Operator_Gy],[Operator_d]]) + + # Construct the full mapping matrics and get the rhs. + Map_pgrad = scysparse.bmat([[Map_Gx, None],[None, Map_Gy],[scysparse.csr_matrix((Npts_d,Npts),dtype=np.float), None]]) + Map_dirichlet = scysparse.bmat([[scysparse.csr_matrix((Npts_x+Npts_y,Npts_d),dtype=np.float)],[Map_d]]) + rhs = Map_pgrad.dot(pgrad_vect) + Map_dirichlet.dot(dirichlet_vect) + + # Evaluate the covriance matrix for the rhs + cov_rhs = Map_pgrad * cov_pgrad * Map_pgrad.transpose() + Map_dirichlet * cov_dirichlet * Map_dirichlet.transpose() + + # Solve for the WLS solution + weights_vect = cov_rhs.diagonal()**(-1) + weights_matrix = scysparse.diags(weights_vect,format='csr') + # Operator_WLS = weights_matrix * Operator_GLS + # rhs_WLS = weights_matrix.dot(rhs) + sys_LHS = Operator_GLS.transpose() * weights_matrix * Operator_GLS + sys_rhs = (Operator_GLS.transpose() * weights_matrix).dot(rhs) + + # Get the solution from lsqr + # p_vect_wls = splinalg.lsqr(Operator_WLS,rhs_WLS)[0] + # Pn_WLS = np.zeros(fluid_mask_ex.shape) + # Pn_WLS[j,i] = p_vect_wls + + # Solve for the WLS solution + p_vect_wls = splinalg.spsolve(sys_LHS, sys_rhs) + Pn_WLS_ex = np.zeros(fluid_mask_ex.shape) + Pn_WLS_ex[j,i] = p_vect_wls + + # Perform the uncertainty propagation + sigma_Pn_ex = np.zeros((Ny+2,Nx+2)) + if uncertainty_quantification == True: + cov_sys_rhs = (Operator_GLS.transpose() * weights_matrix) * cov_rhs * (Operator_GLS.transpose() * weights_matrix).transpose() + sys_LHS_inv = linalg.inv(sys_LHS.A) + Cov_p = np.matmul(np.matmul(sys_LHS_inv, cov_sys_rhs.A), sys_LHS_inv.T) + Var_p_vect = np.diag(Cov_p) + sigma_Pn_ex[j,i] = Var_p_vect**0.5 + + return Pn_WLS_ex[1:-1,1:-1], sigma_Pn_ex[1:-1,1:-1] + + + + + + + + + + + + + + + + + + + + + + + diff --git a/LICENSE.txt b/LICENSE.txt new file mode 100755 index 0000000..e62ec04 --- /dev/null +++ b/LICENSE.txt @@ -0,0 +1,674 @@ +GNU GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The GNU General Public License is a free, copyleft license for +software and other kinds of works. + + The licenses for most software and other practical works are designed +to take away your freedom to share and change the works. By contrast, +the GNU General Public License is intended to guarantee your freedom to +share and change all versions of a program--to make sure it remains free +software for all its users. We, the Free Software Foundation, use the +GNU General Public License for most of our software; it applies also to +any other work released this way by its authors. You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +them if you wish), that you receive source code or can get it if you +want it, that you can change the software or use pieces of it in new +free programs, and that you know you can do these things. + + To protect your rights, we need to prevent others from denying you +these rights or asking you to surrender the rights. Therefore, you have +certain responsibilities if you distribute copies of the software, or if +you modify it: responsibilities to respect the freedom of others. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must pass on to the recipients the same +freedoms that you received. You must make sure that they, too, receive +or can get the source code. And you must show them these terms so they +know their rights. + + Developers that use the GNU GPL protect your rights with two steps: +(1) assert copyright on the software, and (2) offer you this License +giving you legal permission to copy, distribute and/or modify it. + + For the developers' and authors' protection, the GPL clearly explains +that there is no warranty for this free software. For both users' and +authors' sake, the GPL requires that modified versions be marked as +changed, so that their problems will not be attributed erroneously to +authors of previous versions. + + Some devices are designed to deny users access to install or run +modified versions of the software inside them, although the manufacturer +can do so. This is fundamentally incompatible with the aim of +protecting users' freedom to change the software. The systematic +pattern of such abuse occurs in the area of products for individuals to +use, which is precisely where it is most unacceptable. Therefore, we +have designed this version of the GPL to prohibit the practice for those +products. If such problems arise substantially in other domains, we +stand ready to extend this provision to those domains in future versions +of the GPL, as needed to protect the freedom of users. + + Finally, every program is threatened constantly by software patents. +States should not allow patents to restrict development and use of +software on general-purpose computers, but in those that do, we wish to +avoid the special danger that patents applied to a free program could +make it effectively proprietary. To prevent this, the GPL assures that +patents cannot be used to render the program non-free. + + The precise terms and conditions for copying, distribution and +modification follow. + + TERMS AND CONDITIONS + + 0. Definitions. + + "This License" refers to version 3 of the GNU General Public License. + + "Copyright" also means copyright-like laws that apply to other kinds of +works, such as semiconductor masks. + + "The Program" refers to any copyrightable work licensed under this +License. Each licensee is addressed as "you". "Licensees" and +"recipients" may be individuals or organizations. + + To "modify" a work means to copy from or adapt all or part of the work +in a fashion requiring copyright permission, other than the making of an +exact copy. The resulting work is called a "modified version" of the +earlier work or a work "based on" the earlier work. + + A "covered work" means either the unmodified Program or a work based +on the Program. + + To "propagate" a work means to do anything with it that, without +permission, would make you directly or secondarily liable for +infringement under applicable copyright law, except executing it on a +computer or modifying a private copy. Propagation includes copying, +distribution (with or without modification), making available to the +public, and in some countries other activities as well. + + To "convey" a work means any kind of propagation that enables other +parties to make or receive copies. Mere interaction with a user through +a computer network, with no transfer of a copy, is not conveying. + + An interactive user interface displays "Appropriate Legal Notices" +to the extent that it includes a convenient and prominently visible +feature that (1) displays an appropriate copyright notice, and (2) +tells the user that there is no warranty for the work (except to the +extent that warranties are provided), that licensees may convey the +work under this License, and how to view a copy of this License. If +the interface presents a list of user commands or options, such as a +menu, a prominent item in the list meets this criterion. + + 1. Source Code. + + The "source code" for a work means the preferred form of the work +for making modifications to it. "Object code" means any non-source +form of a work. + + A "Standard Interface" means an interface that either is an official +standard defined by a recognized standards body, or, in the case of +interfaces specified for a particular programming language, one that +is widely used among developers working in that language. + + The "System Libraries" of an executable work include anything, other +than the work as a whole, that (a) is included in the normal form of +packaging a Major Component, but which is not part of that Major +Component, and (b) serves only to enable use of the work with that +Major Component, or to implement a Standard Interface for which an +implementation is available to the public in source code form. A +"Major Component", in this context, means a major essential component +(kernel, window system, and so on) of the specific operating system +(if any) on which the executable work runs, or a compiler used to +produce the work, or an object code interpreter used to run it. + + The "Corresponding Source" for a work in object code form means all +the source code needed to generate, install, and (for an executable +work) run the object code and to modify the work, including scripts to +control those activities. However, it does not include the work's +System Libraries, or general-purpose tools or generally available free +programs which are used unmodified in performing those activities but +which are not part of the work. For example, Corresponding Source +includes interface definition files associated with source files for +the work, and the source code for shared libraries and dynamically +linked subprograms that the work is specifically designed to require, +such as by intimate data communication or control flow between those +subprograms and other parts of the work. + + The Corresponding Source need not include anything that users +can regenerate automatically from other parts of the Corresponding +Source. + + The Corresponding Source for a work in source code form is that +same work. + + 2. Basic Permissions. + + All rights granted under this License are granted for the term of +copyright on the Program, and are irrevocable provided the stated +conditions are met. This License explicitly affirms your unlimited +permission to run the unmodified Program. The output from running a +covered work is covered by this License only if the output, given its +content, constitutes a covered work. This License acknowledges your +rights of fair use or other equivalent, as provided by copyright law. + + You may make, run and propagate covered works that you do not +convey, without conditions so long as your license otherwise remains +in force. You may convey covered works to others for the sole purpose +of having them make modifications exclusively for you, or provide you +with facilities for running those works, provided that you comply with +the terms of this License in conveying all material for which you do +not control copyright. Those thus making or running the covered works +for you must do so exclusively on your behalf, under your direction +and control, on terms that prohibit them from making any copies of +your copyrighted material outside their relationship with you. + + Conveying under any other circumstances is permitted solely under +the conditions stated below. Sublicensing is not allowed; section 10 +makes it unnecessary. + + 3. Protecting Users' Legal Rights From Anti-Circumvention Law. + + No covered work shall be deemed part of an effective technological +measure under any applicable law fulfilling obligations under article +11 of the WIPO copyright treaty adopted on 20 December 1996, or +similar laws prohibiting or restricting circumvention of such +measures. + + When you convey a covered work, you waive any legal power to forbid +circumvention of technological measures to the extent such circumvention +is effected by exercising rights under this License with respect to +the covered work, and you disclaim any intention to limit operation or +modification of the work as a means of enforcing, against the work's +users, your or third parties' legal rights to forbid circumvention of +technological measures. + + 4. Conveying Verbatim Copies. + + You may convey verbatim copies of the Program's source code as you +receive it, in any medium, provided that you conspicuously and +appropriately publish on each copy an appropriate copyright notice; +keep intact all notices stating that this License and any +non-permissive terms added in accord with section 7 apply to the code; +keep intact all notices of the absence of any warranty; and give all +recipients a copy of this License along with the Program. + + You may charge any price or no price for each copy that you convey, +and you may offer support or warranty protection for a fee. + + 5. Conveying Modified Source Versions. + + You may convey a work based on the Program, or the modifications to +produce it from the Program, in the form of source code under the +terms of section 4, provided that you also meet all of these conditions: + + a) The work must carry prominent notices stating that you modified + it, and giving a relevant date. + + b) The work must carry prominent notices stating that it is + released under this License and any conditions added under section + 7. This requirement modifies the requirement in section 4 to + "keep intact all notices". + + c) You must license the entire work, as a whole, under this + License to anyone who comes into possession of a copy. This + License will therefore apply, along with any applicable section 7 + additional terms, to the whole of the work, and all its parts, + regardless of how they are packaged. This License gives no + permission to license the work in any other way, but it does not + invalidate such permission if you have separately received it. + + d) If the work has interactive user interfaces, each must display + Appropriate Legal Notices; however, if the Program has interactive + interfaces that do not display Appropriate Legal Notices, your + work need not make them do so. + + A compilation of a covered work with other separate and independent +works, which are not by their nature extensions of the covered work, +and which are not combined with it such as to form a larger program, +in or on a volume of a storage or distribution medium, is called an +"aggregate" if the compilation and its resulting copyright are not +used to limit the access or legal rights of the compilation's users +beyond what the individual works permit. Inclusion of a covered work +in an aggregate does not cause this License to apply to the other +parts of the aggregate. + + 6. Conveying Non-Source Forms. + + You may convey a covered work in object code form under the terms +of sections 4 and 5, provided that you also convey the +machine-readable Corresponding Source under the terms of this License, +in one of these ways: + + a) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by the + Corresponding Source fixed on a durable physical medium + customarily used for software interchange. + + b) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by a + written offer, valid for at least three years and valid for as + long as you offer spare parts or customer support for that product + model, to give anyone who possesses the object code either (1) a + copy of the Corresponding Source for all the software in the + product that is covered by this License, on a durable physical + medium customarily used for software interchange, for a price no + more than your reasonable cost of physically performing this + conveying of source, or (2) access to copy the + Corresponding Source from a network server at no charge. + + c) Convey individual copies of the object code with a copy of the + written offer to provide the Corresponding Source. This + alternative is allowed only occasionally and noncommercially, and + only if you received the object code with such an offer, in accord + with subsection 6b. + + d) Convey the object code by offering access from a designated + place (gratis or for a charge), and offer equivalent access to the + Corresponding Source in the same way through the same place at no + further charge. You need not require recipients to copy the + Corresponding Source along with the object code. If the place to + copy the object code is a network server, the Corresponding Source + may be on a different server (operated by you or a third party) + that supports equivalent copying facilities, provided you maintain + clear directions next to the object code saying where to find the + Corresponding Source. Regardless of what server hosts the + Corresponding Source, you remain obligated to ensure that it is + available for as long as needed to satisfy these requirements. + + e) Convey the object code using peer-to-peer transmission, provided + you inform other peers where the object code and Corresponding + Source of the work are being offered to the general public at no + charge under subsection 6d. + + A separable portion of the object code, whose source code is excluded +from the Corresponding Source as a System Library, need not be +included in conveying the object code work. + + A "User Product" is either (1) a "consumer product", which means any +tangible personal property which is normally used for personal, family, +or household purposes, or (2) anything designed or sold for incorporation +into a dwelling. In determining whether a product is a consumer product, +doubtful cases shall be resolved in favor of coverage. For a particular +product received by a particular user, "normally used" refers to a +typical or common use of that class of product, regardless of the status +of the particular user or of the way in which the particular user +actually uses, or expects or is expected to use, the product. A product +is a consumer product regardless of whether the product has substantial +commercial, industrial or non-consumer uses, unless such uses represent +the only significant mode of use of the product. + + "Installation Information" for a User Product means any methods, +procedures, authorization keys, or other information required to install +and execute modified versions of a covered work in that User Product from +a modified version of its Corresponding Source. The information must +suffice to ensure that the continued functioning of the modified object +code is in no case prevented or interfered with solely because +modification has been made. + + If you convey an object code work under this section in, or with, or +specifically for use in, a User Product, and the conveying occurs as +part of a transaction in which the right of possession and use of the +User Product is transferred to the recipient in perpetuity or for a +fixed term (regardless of how the transaction is characterized), the +Corresponding Source conveyed under this section must be accompanied +by the Installation Information. But this requirement does not apply +if neither you nor any third party retains the ability to install +modified object code on the User Product (for example, the work has +been installed in ROM). + + The requirement to provide Installation Information does not include a +requirement to continue to provide support service, warranty, or updates +for a work that has been modified or installed by the recipient, or for +the User Product in which it has been modified or installed. Access to a +network may be denied when the modification itself materially and +adversely affects the operation of the network or violates the rules and +protocols for communication across the network. + + Corresponding Source conveyed, and Installation Information provided, +in accord with this section must be in a format that is publicly +documented (and with an implementation available to the public in +source code form), and must require no special password or key for +unpacking, reading or copying. + + 7. Additional Terms. + + "Additional permissions" are terms that supplement the terms of this +License by making exceptions from one or more of its conditions. +Additional permissions that are applicable to the entire Program shall +be treated as though they were included in this License, to the extent +that they are valid under applicable law. If additional permissions +apply only to part of the Program, that part may be used separately +under those permissions, but the entire Program remains governed by +this License without regard to the additional permissions. + + When you convey a copy of a covered work, you may at your option +remove any additional permissions from that copy, or from any part of +it. (Additional permissions may be written to require their own +removal in certain cases when you modify the work.) You may place +additional permissions on material, added by you to a covered work, +for which you have or can give appropriate copyright permission. + + Notwithstanding any other provision of this License, for material you +add to a covered work, you may (if authorized by the copyright holders of +that material) supplement the terms of this License with terms: + + a) Disclaiming warranty or limiting liability differently from the + terms of sections 15 and 16 of this License; or + + b) Requiring preservation of specified reasonable legal notices or + author attributions in that material or in the Appropriate Legal + Notices displayed by works containing it; or + + c) Prohibiting misrepresentation of the origin of that material, or + requiring that modified versions of such material be marked in + reasonable ways as different from the original version; or + + d) Limiting the use for publicity purposes of names of licensors or + authors of the material; or + + e) Declining to grant rights under trademark law for use of some + trade names, trademarks, or service marks; or + + f) Requiring indemnification of licensors and authors of that + material by anyone who conveys the material (or modified versions of + it) with contractual assumptions of liability to the recipient, for + any liability that these contractual assumptions directly impose on + those licensors and authors. + + All other non-permissive additional terms are considered "further +restrictions" within the meaning of section 10. If the Program as you +received it, or any part of it, contains a notice stating that it is +governed by this License along with a term that is a further +restriction, you may remove that term. If a license document contains +a further restriction but permits relicensing or conveying under this +License, you may add to a covered work material governed by the terms +of that license document, provided that the further restriction does +not survive such relicensing or conveying. + + If you add terms to a covered work in accord with this section, you +must place, in the relevant source files, a statement of the +additional terms that apply to those files, or a notice indicating +where to find the applicable terms. + + Additional terms, permissive or non-permissive, may be stated in the +form of a separately written license, or stated as exceptions; +the above requirements apply either way. + + 8. Termination. + + You may not propagate or modify a covered work except as expressly +provided under this License. Any attempt otherwise to propagate or +modify it is void, and will automatically terminate your rights under +this License (including any patent licenses granted under the third +paragraph of section 11). + + However, if you cease all violation of this License, then your +license from a particular copyright holder is reinstated (a) +provisionally, unless and until the copyright holder explicitly and +finally terminates your license, and (b) permanently, if the copyright +holder fails to notify you of the violation by some reasonable means +prior to 60 days after the cessation. + + Moreover, your license from a particular copyright holder is +reinstated permanently if the copyright holder notifies you of the +violation by some reasonable means, this is the first time you have +received notice of violation of this License (for any work) from that +copyright holder, and you cure the violation prior to 30 days after +your receipt of the notice. + + Termination of your rights under this section does not terminate the +licenses of parties who have received copies or rights from you under +this License. If your rights have been terminated and not permanently +reinstated, you do not qualify to receive new licenses for the same +material under section 10. + + 9. Acceptance Not Required for Having Copies. + + You are not required to accept this License in order to receive or +run a copy of the Program. Ancillary propagation of a covered work +occurring solely as a consequence of using peer-to-peer transmission +to receive a copy likewise does not require acceptance. However, +nothing other than this License grants you permission to propagate or +modify any covered work. These actions infringe copyright if you do +not accept this License. Therefore, by modifying or propagating a +covered work, you indicate your acceptance of this License to do so. + + 10. Automatic Licensing of Downstream Recipients. + + Each time you convey a covered work, the recipient automatically +receives a license from the original licensors, to run, modify and +propagate that work, subject to this License. You are not responsible +for enforcing compliance by third parties with this License. + + An "entity transaction" is a transaction transferring control of an +organization, or substantially all assets of one, or subdividing an +organization, or merging organizations. If propagation of a covered +work results from an entity transaction, each party to that +transaction who receives a copy of the work also receives whatever +licenses to the work the party's predecessor in interest had or could +give under the previous paragraph, plus a right to possession of the +Corresponding Source of the work from the predecessor in interest, if +the predecessor has it or can get it with reasonable efforts. + + You may not impose any further restrictions on the exercise of the +rights granted or affirmed under this License. For example, you may +not impose a license fee, royalty, or other charge for exercise of +rights granted under this License, and you may not initiate litigation +(including a cross-claim or counterclaim in a lawsuit) alleging that +any patent claim is infringed by making, using, selling, offering for +sale, or importing the Program or any portion of it. + + 11. Patents. + + A "contributor" is a copyright holder who authorizes use under this +License of the Program or a work on which the Program is based. The +work thus licensed is called the contributor's "contributor version". + + A contributor's "essential patent claims" are all patent claims +owned or controlled by the contributor, whether already acquired or +hereafter acquired, that would be infringed by some manner, permitted +by this License, of making, using, or selling its contributor version, +but do not include claims that would be infringed only as a +consequence of further modification of the contributor version. For +purposes of this definition, "control" includes the right to grant +patent sublicenses in a manner consistent with the requirements of +this License. + + Each contributor grants you a non-exclusive, worldwide, royalty-free +patent license under the contributor's essential patent claims, to +make, use, sell, offer for sale, import and otherwise run, modify and +propagate the contents of its contributor version. + + In the following three paragraphs, a "patent license" is any express +agreement or commitment, however denominated, not to enforce a patent +(such as an express permission to practice a patent or covenant not to +sue for patent infringement). To "grant" such a patent license to a +party means to make such an agreement or commitment not to enforce a +patent against the party. + + If you convey a covered work, knowingly relying on a patent license, +and the Corresponding Source of the work is not available for anyone +to copy, free of charge and under the terms of this License, through a +publicly available network server or other readily accessible means, +then you must either (1) cause the Corresponding Source to be so +available, or (2) arrange to deprive yourself of the benefit of the +patent license for this particular work, or (3) arrange, in a manner +consistent with the requirements of this License, to extend the patent +license to downstream recipients. "Knowingly relying" means you have +actual knowledge that, but for the patent license, your conveying the +covered work in a country, or your recipient's use of the covered work +in a country, would infringe one or more identifiable patents in that +country that you have reason to believe are valid. + + If, pursuant to or in connection with a single transaction or +arrangement, you convey, or propagate by procuring conveyance of, a +covered work, and grant a patent license to some of the parties +receiving the covered work authorizing them to use, propagate, modify +or convey a specific copy of the covered work, then the patent license +you grant is automatically extended to all recipients of the covered +work and works based on it. + + A patent license is "discriminatory" if it does not include within +the scope of its coverage, prohibits the exercise of, or is +conditioned on the non-exercise of one or more of the rights that are +specifically granted under this License. You may not convey a covered +work if you are a party to an arrangement with a third party that is +in the business of distributing software, under which you make payment +to the third party based on the extent of your activity of conveying +the work, and under which the third party grants, to any of the +parties who would receive the covered work from you, a discriminatory +patent license (a) in connection with copies of the covered work +conveyed by you (or copies made from those copies), or (b) primarily +for and in connection with specific products or compilations that +contain the covered work, unless you entered into that arrangement, +or that patent license was granted, prior to 28 March 2007. + + Nothing in this License shall be construed as excluding or limiting +any implied license or other defenses to infringement that may +otherwise be available to you under applicable patent law. + + 12. No Surrender of Others' Freedom. + + If conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot convey a +covered work so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you may +not convey it at all. For example, if you agree to terms that obligate you +to collect a royalty for further conveying from those to whom you convey +the Program, the only way you could satisfy both those terms and this +License would be to refrain entirely from conveying the Program. + + 13. Use with the GNU Affero General Public License. + + Notwithstanding any other provision of this License, you have +permission to link or combine any covered work with a work licensed +under version 3 of the GNU Affero General Public License into a single +combined work, and to convey the resulting work. The terms of this +License will continue to apply to the part which is the covered work, +but the special requirements of the GNU Affero General Public License, +section 13, concerning interaction through a network will apply to the +combination as such. + + 14. Revised Versions of this License. + + The Free Software Foundation may publish revised and/or new versions of +the GNU General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + + Each version is given a distinguishing version number. If the +Program specifies that a certain numbered version of the GNU General +Public License "or any later version" applies to it, you have the +option of following the terms and conditions either of that numbered +version or of any later version published by the Free Software +Foundation. If the Program does not specify a version number of the +GNU General Public License, you may choose any version ever published +by the Free Software Foundation. + + If the Program specifies that a proxy can decide which future +versions of the GNU General Public License can be used, that proxy's +public statement of acceptance of a version permanently authorizes you +to choose that version for the Program. + + Later license versions may give you additional or different +permissions. However, no additional obligations are imposed on any +author or copyright holder as a result of your choosing to follow a +later version. + + 15. Disclaimer of Warranty. + + THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY +APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT +HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY +OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, +THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM +IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF +ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. Limitation of Liability. + + IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS +THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY +GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE +USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF +DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD +PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), +EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF +SUCH DAMAGES. + + 17. Interpretation of Sections 15 and 16. + + If the disclaimer of warranty and limitation of liability provided +above cannot be given local legal effect according to their terms, +reviewing courts shall apply local law that most closely approximates +an absolute waiver of all civil liability in connection with the +Program, unless a warranty or assumption of liability accompanies a +copy of the Program in return for a fee. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +state the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + +Also add information on how to contact you by electronic and paper mail. + + If the program does terminal interaction, make it output a short +notice like this when it starts in an interactive mode: + + Copyright (C) + This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, your program's commands +might be different; for a GUI interface, you would use an "about box". + + You should also get your employer (if you work as a programmer) or school, +if any, to sign a "copyright disclaimer" for the program, if necessary. +For more information on this, and how to apply and follow the GNU GPL, see +. + + The GNU General Public License does not permit incorporating your program +into proprietary programs. If your program is a subroutine library, you +may consider it more useful to permit linking proprietary applications with +the library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. But first, please read +. diff --git a/README.md b/README.md new file mode 100755 index 0000000..0bcd05c --- /dev/null +++ b/README.md @@ -0,0 +1,8 @@ +This code package implements the density integration methodology for BOS outlined in: + +Rajendran, L. K., Zhang, J., Bane, S., & Vlachos, P. (2020). Uncertainty-based weighted least squares density integration for background-oriented schlieren. Experiments in Fluids. + +Please cite the above paper if you use this code package for your work. + +sample_script.py is a sample python script that loads the sample dataset in 'sample-data.mat' and calls the density integration + uncertainty quantification function to perform the calculations. +It saves the result to 'sample-result.mat' and a figure to 'sample-result.png' diff --git a/__pycache__/DensityIntegrationUncertaintyQuantification.cpython-38.pyc b/__pycache__/DensityIntegrationUncertaintyQuantification.cpython-38.pyc new file mode 100755 index 0000000000000000000000000000000000000000..dd537eb8db25f859fcb145d063d7cb2028217c13 GIT binary patch literal 10075 zcmeHNNpsuC6$ZdX5=ByrX>0KYZ`!O~mSxF{B+IsBTb5j!<5N+ksReQkd8(_gL&!r!P8{)?irfUo&;fGT7KRaB*-n3ec7l8x|dbSIXLi>|I| zCngD^u+yAv@gxb9n^8`vQ8o5kG25!Ps&TamrKFlvo7EPSQmV8e45q%gomet7#nM!{ zT#}#QyI3(~c}kv`o|~SSo0^s<$0ueIi3MfX+A14!^0uzzx3p5>i!G&8NPLSRg|!#? zR=RqBH*k&PyM(VPHw%I#?2Ak|5I>XbRy)ZMb&mj)dW3GlPBD?!F-l%fs}YtuE!b_i zrg3d&ZEEy{V0W-~e(i)eiM6WHE`cQpH>Cq{KcdDCMEeAj*a;fx?<5yDrCZedA$vQ+ zA}sS>mM2-1oqQilFNRYr#!kJDz1_v)tn0ljyN?`0&ymtz)}%InF8y9K{>YkHccXvS z!g>zL60EmT2H#Q4B|P!_cG^B;#n|cJi1t}%;_Sq>$Y;IfKx`H4b5;{O`x}v+W2e~} znsdL@BiQHJDOZZJ?8?7v1MGs@^m&9`B+Ll*AiKiOvq7v*j+;_U=n>Rb){zlrgouD~`T7fF z<>@0w=r10l@ALeJ+GBku1OL80)Sk7q3Ly)MSo(lSL!Q)BC{n+t&`7B#Z=Ak1)HOUpCFIK_ao~?UNYab;wS!>#&$% zGinkMkNcx#qxHWk*jL9JC7Z-TI=}h1BnSyvepssPTIL)|viwPD_|wvyJgLfxVJKC( zyeaQ!=9XMAbahnbCY!}wUCr$%=I8Zxuk?~)RAo|7v&~}pjl8LAMYOwH-Yu!5td#YV zWx}^6=gTEiS2Y71NH59Omgcq?;O)(=@@`R;Uup6SV^r%i2T2m3|=Pq zJHGh)$NaZ=J_INFs$4E|qi?^%-wP*#(B6tpOtq>Nt-x&s;m=;Hzs~!T1-~GyIgz|+ z3`$PS+)^r<6Q}#7Sc6@TRH_(l)Jg0M%_y5r9E%A~o!FZ#JQ+@;sFj>37MJP73^0Gq zX)*Ivv!WQL<|N<~kD+Dgdrk~bn_@X_T2bG@V`ynP*Iz^}TTZiCDMDB*I8A27EEiwH zzp5(*(`m8b)vT1^E$*3MaIYP4(TT(J&Ti2WUpP@R9=zd5pOhTw(~^^@?Ep^PT`S(T z>a+#x!+q+)0=}@F-q6izSkq22wBq*aYVh2!{v2toitnBysr0QnF+`3czf_$#uT8qw z>ZFsZ1)bCOoEA6kTCqyS&W7k4WPAW5^J1y_2n$?K4bf67aoZ7fC%RU#OsCZi2U`ym zlD?E})|b2)ll1#pN>qHC#8c9ji*EsN{8}^c2snVwH_-;)zo6p3+dcn(ZKYZ48JP zOQ@XBwtA7i*7BBDspfbkG-6o0vdHt5jcjbagqXAO3YR&p{sTwC#8UZn(0x!w_f)Ry zKHEeh@|Oxmk)stBBYDM1yi|;w+r?YKc!DS4UwdOO%f(&%=!_{E_5mZ8A6qxd+lbs8 zyZdxwjBmIjCU!SNGe^9wVNAjLu9~_zW*YgiH$`*gRoNW7=k6O1YrEV>wf*#2-F7-! zsXG0hJpv52#oT&XH_dV>7h1LL2-EzUp8L+U*p5%L`;pFw8R?Xm5$TIET~Vh`0yTUh zov_THMG_xsiI0?tab2ukl<5iRtg6&~`)8E*VcTmbz>L-3jCn+-$-p$SZ(v{GVkLypw+hC;#`&Da79| zr(mB#p5lG83DN&bHsODkL%@A*Mg{_>e4fn#QcfU{Bq+;65^&J-o^71l@)m(D<$?kU zY-QL++4&%|alyWb%*F*|HYnS0fnDU8hs$m`$W&ahTDBAHl9zcvJ`I_Nz9VG)tPi=6 zOUQ*>cC!PaQ?swR;%odfQ+*3H*&z4qeL0VxpQU(aX22dq`~!%82r&<^!65zto=I`* zS9z>Mb@2}$Asb=Ch-N4hf9-VXYp#0n`{%lm`WAki%$IxizMRLO3Fd#)9z*=2h<_aM zkFqg0ejeL+UAz;JC2_hKS5s>Gs|agD#%7XDp!e}wyf{fh?#4a(QxC&jMO z>56-j?CDUB>U!wxH-lKy>^hsF9HV`s{$ys)vbvL)J;$wQ*&MgN$$Po!^>E8-*>7Wg z1(6NoJ5dk|5qP}KrTuE6fHlfS3ew?GlaAZoAzQA$#D;KI8|OJNz7pYa5*{hWd3I;N zo!!7`bspC#enwqjci1hQ>+bA#@KfJBn?`Nkt>Fy2$QJiI4+OOdS=PHe>bqXlOYAmT zu_f5w#%$lK)!}^s{oS`t;6!+j-N%#l0OLpgl`@0&GJC+5{c|9fH)L$>74Os=H^ysT z*&_Ai_vuDxKeSHbgn4hf3w|HL7iB}c*(3H45{W(9?m^#=A?sz2VX?|q`6+b8@0~3- zo;le{qlBMA+YBPvMuZ?nKuK()iGY59IPW%EPeGNCV!A2K^krJ#BQ5 zh`6-TNpON7LvWJd6hRljX)TruI71b4N=Y@kN!&xwOS*_w)r|Ap3dd8N{50ayH~|pN zK0{_QDVvE9m%%uUvm`r5aDm_=!6kyr1Xl?90E17Aeku`%Wm?F4FFimisV1|H_{>Ci99KmgZI|MfgZV}8AED$Ub+$ErQ1japr`v9jsM2(mC zh^f{|IG>faiKgy3abmBzJ>vnjUZ&QqJ}b>FLrnw}S7{sETn7v1=2c^*!4ev7S#=`Y zITa@=fNQHdiNK-i#5kdK(jig_BIsmlb(fbCmvvJ4F#U9rl@ROoYcB0|60VLnYo%tM ztrMzLopepZEn9-RlMI4fUfpn7C`|6cN%4sN6oL~mwoIpq+EF3-M$-#hm_Y6>?bp0R zT$5QfbE_*GIHsXM_6e67=Rjzk77+4}S2u8bj7NMn$*FJH>NO>9Z9JxXvr0fT&V@B< za_K07OFl#7_D?W0{{T3O+}22KnB3M;+sLpu93`4d1eY|uAes1Z>U*?|tjX8RrTnH`})l#!Tq34K4~6s%K>Q~1Yz=ZLvA_HT_vj8EpQSR zWmoB%5xP*X$um?N3TXl>ycQk?(Hq9=;}O83GH*Mj+* z4$V>^KQw=VRcQXMv#Ewr`19w>{V`k*&0mea&)7FW-)BJIXEA3p>_&*b&(_V?97dmH zvuqasf{^m)`%N~7-e-@c@4mlVpznSeq#5-6HeWg9NwV*R==*$#zAxZjr`bGP;Dmg! zfxh2upzllE`Yv1I*7tZX_js%lXs%D+LA2SNN8caR(f28jzSoU9pzq5e`aTZ&K3hlM zZ+i57g)M_vFM__W;ChSG_lF?T_d$!7T{_HG*d5eX+#2ZnBlaku?~i%ZkG-f@*#l0c zS7HAE^!-V#4(|^^me)9af5O&4p`T*>K7C(jPuaRppt-z}zCR0IT8-DdvPaaD-=`a) z{hZVHC!D@-z?VzkH`sGd-#0jYe*t+fdx3~PVjrn#yhM8D_s-TE$vS)1DER^OotV3^ zPQZ!$Gn58j{vT6h<70~T6M|0(vIIXN2zapZ8C8D@$VOh~mI&@=n_PCfoCD z5Ai0({0saFnZrre59gs^p;Cr(QsSTs;HUMR^!xNbn4kU{qZo|7SX*1G!?e?s%ck literal 0 HcmV?d00001 diff --git a/loadmat_functions.py b/loadmat_functions.py new file mode 100755 index 0000000..721cb44 --- /dev/null +++ b/loadmat_functions.py @@ -0,0 +1,35 @@ +# The following functions convert an object to a struct so that it can be saved to a mat file +import scipy.io as sio + +def loadmat(filename): + ''' + this function should be called instead of direct spio.loadmat + as it cures the problem of not properly recovering python dictionaries + from mat files. It calls the function check keys to cure all entries + which are still mat-objects + ''' + data = sio.loadmat(filename, struct_as_record=False, squeeze_me=True) + return _check_keys(data) + +def _check_keys(dict): + ''' + checks if entries in dictionary are mat-objects. If yes + todict is called to change them to nested dictionaries + ''' + for key in dict: + if isinstance(dict[key], sio.matlab.mio5_params.mat_struct): + dict[key] = _todict(dict[key]) + return dict + +def _todict(matobj): + ''' + A recursive function which constructs from matobjects nested dictionaries + ''' + dict = {} + for strg in matobj._fieldnames: + elem = matobj.__dict__[strg] + if isinstance(elem, sio.matlab.mio5_params.mat_struct): + dict[strg] = _todict(elem) + else: + dict[strg] = elem + return dict \ No newline at end of file diff --git a/sample-data.mat b/sample-data.mat new file mode 100755 index 0000000000000000000000000000000000000000..bed2163e336018c3eca4df6e12a514963f0eadfc GIT binary patch literal 15560 zcmeHuc|4Wd+qWSq88VZ39+HfiU27XNj~Swnq0ACRlcYErNdu{*5fv4Z1{P9DDup6M zN=l|Zo1gof_w>H!cY6N#y}v)sd7jVbtj}lfYp=c6eP8Pu)>_|n@9%OVc{-9TWc4+* zWSvN!E5gE~LSd9&w z80)Sw*4LHQ(bCak`40-Kof9u{nHDPx3nw8L5wFF+--hG&`cgvv`eJ2Sz!Dv_C6eve z---4=71n<#M4Nvp|NHBIGWb`&p+EZVH^;EB{OmUt7CxeFNO*jB@TREHgn+1^4WUsi z|GpppsF+arFTVfZQ2$dN{LviZ^KX7667NU;{VIPg|D7Lsi0|V!#RmLeUlR4-@|Cz= zGF@HQwvhs}BERI0(zygCknjwjlc13iUQ{KlUrHLsF3HTt?cDV1@SYhH*z^s zfpTK)>sP*1XtYrMB(;GGtc5v+%hyn$aQ_jHNlPk3R*sy?HKIaSwV2lpe&YMv-WB~H zDKN0$NbRRl3Os&uknQCi3S1A>xqS5{1%x9FrE9w=@C$kW8ed*2#B*524FDC!Q?nl7 zbSgZ$rGU5Irh-D|SETxq3XBax&kDP#&?(_MY|u}IV@EE(E9{_xra7ODa1|BEF=|br zn~8qe=CjQj68-Pl*I+(N0h#<`RuB6qpfRml)zD4>owL_c$0-zeS(+CnDny0-jW4@6 z9H@}qEOplTG!-%rC8kAvqe22F*=0(S2K%`<7n7W6aNF(8J&G3%KA(OIZ{le1$oxK?cuDMy3Xy!qwvkEoEMz1tCrsIceOmBSwi|3pbu7wi)GSS-!jr=Z zZAU!_dtB5td5O5Fp6h!Fs)Qci4^Lq$6~gNm?P(#}E39I^p{%7rT1&R-fp!{PD!IJD z!jca46%(J@FVn#;Zq>ZOXFAAc$GtsUO9y{@{{2Z=bSQQ`r4ZynhXXxkJV~`QxXh^T zc5|n}G&y-%nN9_RVBq3^wfU)a?|#Q1Jd;q2DaN^0UXa4&OQ z*;-12>ZDu?g8({IblmXF`%Z@^Of|Rj&J2i^xlrkLoB`aoFw&4qOEiMBTrv|<1c%lG~z4;i2*mU*CNlmUEW z=!(=h11xq$mQ>X+V8<$z)Br~YMDtNNR4&rNxp+@FgHD4_7F)t*eQ02X@99YP5b-wJ z_4bB45s$|A5#j+<7;y=DnlwO#Yg}h;=?2n3LH$F)RVg~S3huF9(n*JdpF;D#?O=dt z58Du9DH9$H8HSt;V1fyIpk75V6YAEss$Jq>LhknD-S;jMez+`>yrzQ=sU2BGiJwChm1OpZoI||T+n9x{g$LHZX3*C$flWv60zaX@mnpD5KZ^nV z>SZo2eRR;s93|B&(!plQV0S0h+GJ4AzrhVkDn+@M2?)VAhZb_~$It?YU4J_GWuyC2>u#suen$(n6~ zOwg|UJUCI!fSv;(1NTfA;P4q2B4$(k;P2VxAS5y$XXD2=t zN`;db-d)rqQDOVIukM#aRFD%3oKIOngEF7@wmp|<@buj0x(`G?v-YkJS>8s69sY{5 zYaJL+Uf@qI%_QOhVp0#{7!YFI3G>>WAyhV-IyAtExjeot%`Y|d5NS=)EbD}{A#ng6n z4-JI(9X}NALkC3DB=J?!LGZ~t|D2C>koJB0HUAAAQe-rwJj>~D{%LYut|yT<>%3%A zzYzHze^J@Pmj>+XC|qnKRB#SF*y5Ey+`s)^YnTcZ_~VB6gc0-n_LhB5oS6UQ-#ROvr-Ilc%OI}xRFL@u>SBc55`7j& zl6a`FWb0tH$x%^hOihbzJ!_zXNd~o@O4#+}RQleXgkJ|HBJ5I#ykfI64)eG`jHi+nr_8q#>(<++ zLn>>jaOs%*S%w7R|0227%oz&YyInAq*Gz$xLS79O)f6ZQ@YTEioB{%0&)&O4tWOy0 zy5tEGap1kbYt3>hR7Bb(_}NjR2JDpT9I3FBEX6yPLWMi58jHJ#as2V|n9^BdJsG$;rDj%`*ra9ec<=|!tZ+EcfCL4Kl=ai|3CTn&-nXi z{`qJA{b&983)kQOAKc*c-{6KH_~EaXKk&-mC4bq^Nd5zUSR84)yswTAt>}5Zkg2~E+2veR$eWc#69t7| z+z+Xu;0D7Q?5u*eI>M>(JVm4!np<7>R1C?J!zSySX258T{&l$v+*rSX=e^k$Q7rua zJf)&v1y?J)$x>@k#G5jTl;yU`V%vpr%-&6MxNY$LkqhC=aCm2gn?XWR-O!+`vN}9aUENusA~pZF^~{|6o85T{5`AERfMfpAM!Uy5XRUD(WKcA9d$W`h#e6QX>!+o^+-szh`RHtOM>_8i(}fqJ-DURAH|v)bV!ROpBozWqiO`WM&ab13ka`;*jP(eZ=#AUSz{FBh;txSiL0>&^LiAUQv9e z=q-O@zfisrI=uOeq`!w2+8r&|Sv(?$j*7o;q9lmm8QGFsTlZ+=*z`h8X0#z*s_4l7 z*$i+{5sCRq)eJ}YyJmkqX@vEKbR)2=9@ctw!8<%?C0_2|zGOsFA5BmWcWA6LMH$lI zZlF#=RM{c6U2;~ap33{VvCtAJjycDwye6Szv#e?9CI)ErIrqJ;S<0wa(6Ms8v<&W| zOBA#ftilfmJV?>a2s^Tmaju@Vz!}Lo&%UQv;-!nW@86I?!rI%~o}`MH;azyUSW%7+ zZZTS+bBlEqYL4qXtW|4)N?mB1{c3E`r4WsPvSoHCV7++erVM*@!tka2t9~0~Fg$rm z=d?MBW~n_;(W!ze*YJM#W5{Ff3k4;r$>!KLLUn4m+#1&{Ie189za8dQK0h6{$_{&O zOe%Wwiw#~`uf*~y!vZT_aWAh;H^9D=czvqPG*9s2qiFLrOv|epQA@&@^C4w-q?>U6raB=<;Y(2T zQdji&$pFh5QET+M`(QP!(ne1f=+O=Rw6I9wMw&*BH4Z^e;hPq_;B%JE_Z~I7VIAWo zl>s^K`1y&=TXr*Dv2}L6OCE;<)tSQTi_@ zt#H7AHN%mx^SJ5KOA31KxKL^MO?SFGHV{c0d|cy-<;G4dTWjHfPxl0wNL(P{jn|CV z8!op%g*Y-~ma-b9>**C|*}0&7zH*0)Bi+yg19_XnkKEAGtS_FGtDO=5#R4+uR{hvGKx#9aVT=eQ?L>#ym4V^P| z!JpctkExc}V%GSbhecaR_>f7E)^MsRipw8w)pxN$J-Xh7L2n(<-m~6i)KC|+|KObb zgC=KmGr0Por>H$D9J(}LrEh`Ad;8o+H!Meoo)2?nC8^>)>?5~*wpid5{>=$I8|?7| zzdcZ2>Wt6ax|MUc-5Eo1)B4X(9Psx+BgfQR)>u?Px^3ABz!H6_D$f?ILNlGbx$p8U z&{DAl(It;;klBc>8qYpEl=ptnv{Kp*rLyh(B1E%BkG5JwO-cg_t)#z-FOf&7U#!l* z5KzLGW4W5Nz0L4-CG*~gN!HlMSB5Jn#188fZx&OVv%_q~am9v-wm8#z{$z@%1-4sq zdC@*eL)>+Mc|A2$3o$bWVnvEfQR~pQtrMQR778e}=l&}h`j9czyppeBJ82>i7X_p0E)bKc_Um4awRkiRe%P3>)6K-R`{ zwiBUZeMY#gr(ygK;ScA$&%yB>=Gbr(cn8>;-~-V)Mf4r|_?7abnImFK*y=7$k?2)I zX>%*%E$fz})G8mH$;;Y^HD^N*Ym^?^CUV-tIZFpQn`SgfXX>EH({mh;ZmS{gbvy6P zyc0k=RnpfMPcFszfty=~22^oUbkh}O6>aRaOerdHsXmsH?65qYp^eqHjKs~x=;CY2 zsl8mKE3uxqk$JtDEG|p@$l=N=jL1jwU)6W>q2TS0l0V5RA>NbEG-mY)xx8lmD>ha1 z{dS5%SFSSJlD=FrtzRBh9pbhtj^sjy>WXF`*Ky$TO&=dCCkWv5Mg5@>wn{kUvH!+u zUS)jW_wk41ce40I$5H3VMnxRbG8pf$M-EFKJelGBiVOcI_#x>J_#t5LXzGJ6yeQ&O zVW;X_3H18*21Z)968g-)LZFsK7hPza99ePG5K$aFqO^n&nidJm9v&s3>!Fpx4gDl^ zPJaHilZPehN)5dRlXhA8@)~fU7=-K5et-8Gy$gH%e z$3$R_x`q^6AnTDadk$(zdoJ>8vO()QU7=~s~zVr*^D-DTW-u2~i+ z|KqY_4V~tQ=cxYeSvvzH;l(&+*Siuy`pUHT>GEhK`()y$Rt2;mB)Cs2NE`XhM(rO` zH$hQ?X=8^Mk&w3Xl-|Z}8?-oXG@Rn&fcUlA4NptEAX)akYNN4k=uFGmT}5AA(N)f` zMzSR?=q-IHaob8q#2d4~=F&?W)crM{ruQ77kF&Q=eC0Mq%J+XslD(mWhDANAir%TB z+a@Oz`%+YpxsTo*fgnB9wTS1M&MH$R^jNVqFUAse*9oks8MZ@>R7Je;sxy)kIVGI^ z*bU{Y?rCFP>yBoHx3->n=#I4XBMXX>+>prfgW8L;oskE&=#7xHLw062mYB>~qVF!N z(mz$3p>jR?Y;KSaGPkAomy4((m9?AC^4Tk+=XaAczm(~rG5Wn74zEp+bE&d_$!;rj zhkNYYrKhXWprdz1-bp9)^61$8+$mS&xT|J?+(vhlKj?I?!NDDQhp~_Lmb#)I);!)( zUKjL?qCXb)el?1F=P6K8Y>BRjS;RLsn4%-31D95;)kSBhf}{nN>ge7Yi*wKK%As$o zrj$jBS0e8>;w|d#hG^hP_2W2SgifCjg5~p8XbjPvthd^uQLZN^x(1xkDar~-<|St& zq+jd1Q_Tgj@@C57VJF16o@4&l$R4TJeb|(2V}mlc`F-5}8KE1+w#sii4bck_n3`O! zjUqy~aK#Yg(V4{2H+xVV#h#suUvga)`8CsDz;jJRf3fSu`7C{O@b(iKCfN)f*Lbu0 z)iDdyay#9kroa*zl77kGCU1=vi?#FYjI>6x857Lk{;KU)obJvuNg8Ibo!l+}0T_CV=b{Gy6VH(IBWj zZ?Djmg}5%HFi4!T~fg!{e?-rvlo zfhXP;s2YlDV%*NUyWd9(YxQh@bg^_fzCIkeNhU=N8@urO8`Y}dU~}O_M>|=}TC}f3 zFJ&2Sv(IsE8eNPH3cpn}$qd&T>lh{`%?sg|AEq9Vn-%bj*4~?M*J|OH8_V+|Obs#S zo~iy+V1lp5*aM##2|o}WvgNU{!W$kPVm_?0#tLys)WkGv%ze>fPMp^Y_hmb7h`UF^ zi+5>mTC&|7%k3IydK_$sQ()2d$1~b^_j~PWlf$am>CCF^wJCCVADh7=c}H2?X{h_K zaaJ1}((_e&yA5%}RZfeUaRT41E$n2Iw82`F7v4`Fwa0t!#kqJoIpMwWnhYHe7cB3; zdXAmxj4u}C=I7No;nvve-IJd7Sp4onZmD`3e3Ut@ZgmlGNdFN(X0jnJY~9<@*P)HG zclU!#hJ&$3nY(ZiBU{H5!2&2jVGhN{2^mUvmQ`Za2{9rk;8D)+m% z3ofypq|vsx;hW2QT^AL(;~E!&DkC6%$FTilM4gVH$H`ORnM-4sV3 zn31k#Sl~pB)pq%xt?|dJSETwb*FGPMrH15gJS6y}tX52)emsTvf4oWJr~eVf>W~M z24pRq&fS$i2{Olyn7mdShYv1HfzuNdILH1}k$*KA(wAwZe|S3%e>OT$D>!ZE0fP5( zqWW_;!3(-z5}_@uF%7yo-`~zU(P3Rlub;%C8Q_f_PL9mFPd|zhO10 zGq9Zn=~bT74(FbNQpU*8`Hz!e5ftqe8BFklOzoMU+a`cT_T3@TvT-=0w9R_;D<EP0PnAs^whBC|J)!t$gkbEp&U#VdNj@ra8ysAL(`8uc%m<0s?PJyxJ3X=u> z+5V2}Wfx+(I8bh@(e`~cbP#8LeX=ug9F(uDYKav9YB|TDVR%MWyOs%}@7D=w-zE5?DFHsEV&h8Fjf5{|MVAU6E54Q(X@O-EU2Mj#vv$} z4mZqOr)~yQ;Yu|(oj-}=C;cm3s#fZ`Ov0=8JV#7b$w21luPGXz1X9-N)g>f?|H-*J zW%>*Ya?arr*q5;YB?akC1*fwhX6&dR$!7xEyawO9S`oa}TL7zgImhxq0Bbi{$j^Ge0wsi zAU!q%8A{{N^RA9UtE>3Jdz&Y~amBYY2RjJd^r!NW`#c-_FW-PwIzarww{ehjm3V!e zHV&oFg~qmsE=1#{YgTU5;6OF2+yu0>7(m+}bt8CS1_Bb|B%Ec*u)?5Z$v%TQfQ@IJ zQZ_RI1&fAnCi+)uU+3oek_`pyCOe0hPDAB(^POtolKgFMGl|l1h z;3A|fd1a+%$~2stvr^efqr>!T2JU@01tyn^4~{-yMWnbJ2aT0xLI2j^%lf8CFhAU* zX%amFQK=KHFO}vYMmj;8wvi57131>%KjT1F;CBWUFGg7zFWiyoG^F_y$5d-C{7FBb zjMYn}Pb|o;JWAwn?i9RB6p7kD^&RNBUU|1}PeMWGH<1~F-)QC(ksX{m0t)x*n=0I< z;BK^vR)+;05;xnH4GXZLy5ON!hqV-lE_>_sHk$%9UYy?o3EpgbwuV6K&KVH=7>qux zWBe3`@_lVgi&%q0&{&bhv-5uy_T734XbT zm1iH5!I9CQFnMDJ3eF5f`q1f+P&^USJ~Red?!iv6^Ta$oI7AC8VZimu88P-Ftf+X! zs)P04W`5HDr{VCIReyaB@t@&{U!O+&&%+U4kFwSsZzp)`1^v|JXz<9Zbw1i)#xTx|X_rd;AqTI4tj4AN+v~mJ?O%sbiDi`B8kO6TwqI zDQze;8^VMQYNg+`qZT1K4aJ^kSLR@X?tJR)Y06J=Oz)XK(V|R&Q1;ZiYqJdSs3b*) zk(h8G#-ygIjSZQ2`JEksX}IO%dQp4>70kbOobKAs0RHFzm8MHnpefaG-y-nG)$;I> zjAS}Uc-++D>7hbX#GAuzA*{%1d(zf73s_KIS)s}UPufrVk5!LVwYiXiaZG>Nm-ln< zvPjC?Wg#8B^-fT(h|Iuj_^?~^3=0Z&>u?kHn}mk0CNAS#8mx$%%1S#!1JiS9{7S^R z2*-D|w1hA^bfh#r&N)Q}dim)sJd!NvCZ*m_aCQM&uN?H^5it&aOM*YVX`7X0%RZlk z_s2MyU79r5^1fj3(!NP(|2|eoDVu_Z+9zsj2Wjwq*%it+`B9KxbVpEm-#qBC?O=9z z&Oo{yk8b=17BsHg5f-k-fuw4V?4ArFc=wzuXU&Ln7=;pcJJP#m;KP)z!oBRD?H8X& z@O!;^8X8c6G?t)2rBI26@{Ms|>$aXG84~9^>U~FJni){HM$tj|8XJldzU*d2V?kkG zBJK?3O~c-`Ijau%&Vi{}gzVwd^N_(S5M(4(au->1TCn%5ZzB#PMg!7g5)0b}%>q%hS*8=rZD&N;x^bZ5+RH<q~9=liVnS_H{`O}_U zvk)MaU8~AE`BVH0?msPqlZwjV3>x4VG=HZkx zyU#}A{0WGhxb~)b8UmvwuAbxCwJ(9_AcgK|1Ibz)n@=aREHAn{U=6NTL2r`h&Y_o5( zO+%{ovYKN7lknoQQ|H3gF)-NCOF6xXIOk9mB{C|@fl6OqN!~-?rRB`z6$GqCfM<<2;n zNhs%(nakiO!&!Fz>Ox{&nX0OYY9j9choSe|$58Q6YmdkYINz+Ol;=7J`UQ-vvLa%i zkZUBwAvyuBN|8Cs1j&&6MTYXxbOu(YCTppbOu*x7)WuJHsBpN&lr(sp4iR}DHg}FQ zp?o$_@5ON@Je8{1hkMx2>(V!4{`pL>zKnOd^8ReUp@77L2Lvt>(^`I3SdI!u2hNdb zM4oZf=1r7xk|D+`{(@oqC={i%xgYkO1P{$oxz}+_NM0wC$Wz3C8=KCiSb34*wcVJN z&1Dv3-?Y2QuZjY1L5t<|!JF^E&$Kbg2i~Z3Y%ZmUVoNp9H&^ zh%@RFbFesV4fC=24Af8T&AIPDfuU|@Wz#wuO!cuO$vICz$10h{cehW&!~FyE2^Fk} zHH6DbTb&j4xkT1;Da=78B|Sb{d=9!Dx2<6C7O~8wFd@ccmiX&Zk+Z8V+mLc*=ruyG#MIbMuw1* zsFb8a9priTZsc3u!+Eby=e@k|cU|A<{o`_+fk5DPKMN+D& zT54)qDk@THDrzch|A$}fCKdwtH5DE{oaA`vfuAr7%%&u^bN$H^bJ4m?H~8&k9+LLee$EO|ERA&)ww_N{zv}( z$b%nz{K1PK=krsX{^%S2#{a#K`8)o<_fvoK@Atm$@BI6{zx=!Y{NCsOU4Q>qzy0_9 z^ZWesNB{fp-oM}ahUNd3-oVDjhwlb^2e@hgyCtnL43;8h&q|(aqE$T*t|xz36Q!uD36=A<C02Qwl6wN_a!enNaGukT!gd35tvz5AHQH!Ee6* zyN4}I_;~H*Q0aLlBtLxVzuSTdXC_5m^m|4DS!vwDh4I>SH23A)*Cc3oXvK?f9E(3e zZu$(tOC)ezO$>V_O@a2uxEARG3VLgH&*U^=f`{+y3&D*{==ObFsJDUz6GTi;)$U+H zSnNkR6|`Bc`xX@%vA{v}#;*O(nNZB|%wC84cpY(O)UuF*k7cHE0}2%6{wB7d4*Bu( z`E!OH5qy4|1incrx>M2b6DIR>#Re#lk7<@s&0qpOM%6kZY!0JHPh+F7M6s16Pvxs3a`$jyEG5XK|je%f_=b|&0-cW#%88wG~t#cE=$ zWAXc}+-f#qeAHX^`_u(fpe5@o(<8xzo!(}O=JiZy%ia>bFOmh`d%_Y{un)kykCJo) z(*d}!BmIJ@*#OMq+4R&K_p#;5VCv9z7L4S0cW=V;$-Ht!dJgi#t=J`L*Cq;*D$*y6 z;`#pU|M0DP$NuLe^vy}HAMi)Nr`)CSh%>>xXy^0k4NS0RF;Wt@vOxBfOjI4)04UfL zG&C3w0GnFMWLezDr!D+nA0qyo&#|hGQ7qsgkKR&$!34*H1uEO-G2v#PbzC6ki9Tg0FeZe4$|qlI3sQnn**x*0d)Z zJ|2r7{*hSJ9Z$m0wD+AyW>OG)s?KQdWeS7>XMgoKWx}X3`w4+2CaB)fm_LB~(4MO? z*#Z6B+sQq71KP5Of)5+<^UM4;54Ylb1(&a$EXRClf7F(qd4&Ssn75o}$m@?6_b@k} z8jEi)s^x6GlLV9L+jBaxZfxQ)+&i>}0X5Ew5Z0oTe|M8wF@1K zA70xRzD*nJq_Tm2IM%DlURT|EK9g`TS0*GB>j?d3eqAE^p%`7_8#6*dl#oB@$S|R5 z>0Y%aJl{JVm6;3i^FTFDnaB8^TT9y3Jk0yUw>7!hhe^n;4r8rMCt*glf!ieKvG{yp zT}9eLBy?&TO=&Q}x*M(Bw-D=Zv#FTH{w@-7JDv*9!}{yT`Bv$TJ=P`jZrc8E3W9yS za-+9U;J>$DEYK12R^XvV<6H_{-7+M#J4i@ht!bEdf&}q(FVsr0?iB@ou~C^c7XN5o zhF96pXZU)g-})XG31NEfH;x#RAg6bKvu^?kYt{%Zxt2r1^f@^*?p2YXy_A)UU8UcF>EqP7{0JDMK>Jpt!FhovUCzoCrck^&Bpkjn`1mMGWP!Y zIUoP3=Wl(4@xR$eOu=_UeAcdabNRWO`0wi{wg_Kd`6_%T{8jen^ZJ3?DDO4Nza|(o zB@#w;ngii#;Y7Bs!9a)!RymzMGYB^JH93?k27zbp>9_OMg5VIRSySY~AebY0evy}K z5Gc>rc5Ipw1PN-Ei&uOK1osufaS5G)5X^fkc-4d8Km7Yw_zfeYMk1|&pzGxu7~td$ zWj*T#`JOO98lqiW4>3SUaG!?B2?mt=CFv;UGhirj*14h_2C&ypjuXgaK-iXJYn)OT zkQt%#M7@Xsy7c9OF${0$CX;(Nz6$&$zU2x1VY8t?c<;GUvALK5`<33;M)-Sx`lNY9 z^^xwNsozKMP;!UGM<#r|%Igk1PP=?dgxz7Im1ONJes>UJf57kR<_^{_nLQFo9F$uZ)d4h~xj;#@)Ol zZZOF9gfnuj3p}q7cpEI>1UHXQDyZgngpZNs6H-haVKLiThyCVGz`Zx?h<1t#R9<{= zpC5Ulej`&gs)7OJ#nbQ_*Zq-a|BC;GjRJm-GC`1(d{;&15(Dx(47le?xWkEqB{QP9 zUBT{E*a6!iC-C}^5xbId1m>--dS7uTh-{XB$$aPpk`H4S6mhx2+@)D(8>SNJvlRB3C5o+~iRy_;9uae+Xq zAo&D;7s$V2e{gFA@*^@>itY~6@99f_Mn0~J&^?hc*&Plz=2ktz z_;4F6N)ShXi`$BK869fFR+&MSv^g{-mU+A?r1MAqY@HLdwOE|g$;f^Dh zhmxN+#MCKi{f7KaP7=I5u#W)|YuEGK!E=eLPk#}g#{jh_O_NL?Gr*KX%)h3F0ZrvE zb7UGBV5z%v=myCE6P0ki;$7ZAvmEls@BKA?iI<^kCDz}K0(q=sF9Kou;gvf-;rDii zO}^-Wb+uL{#AfP(An;DOTzmxcL?y0a@?4!D&_Ay$D5Dw#97Efe?!r75*}ddt2o(rN z)ho(;+XJDSCz$Qh{a@GLAMyGBwm@H6Po{KKf#Ufe(NMQ`i1=Q{~^X}-#_El zZJ)P+>(giG4&U3Niu3wMjgiB=*q7QgDMra)U$^VJ&Uw}IBpm3Ao}TuQ1lqZwtH-gw z-X6Yl=^WG}_FDLDutc4Lc_(Ib=|vLgnxz*QTd^-J*d$?VPQs$Y_Q_%5s9Std?=WJG z?FZu&NtN9*Nf3E@B|8lJ*2XQ&7t%Pt&)$_+%tD;=a}JCuW1kCKSKhpZb4}X2aHVdX zuj-}+KaoK_LS=qoOD6XH&l+lna~Gg4!5XOe@P>qSp$|5394BF^!;Ew}>`U({9+Od& z9@`I#N-IRHHO78fi^E3*^#~&BnC&y12a3HM)+=HE?OoU}{u<|V(W^YYfPFtK^9zZ^ zd4N}`W@r!U5=k|cq4r-XD9Uq-{cxRv@@*X}^_?g%&3=3-3+I6bE5?R{I1gm7R|zLs zkHzoSwqujPIG0i}O&HvzKu& z&O_;fC)bDJ+{^j&M)N4@6nCH7C|pDxV1j`{AKP;ZoU<-$bHh1d#wg!PL!2u*hu_od zwvEN_eI?TDa+U;ATG(triUOtkX|E&mQ2(hccfY%W3BoHTrgY+5+~m8+D;wv&F5A`9 zv~dn3Hh12*jq|t18lwhX{M@4)AXNiMZ+!qR_dCosJ6XOu- zskdbW`TH~e$nm6)EiXt=Xm6G04WyvL`Q2`7)B_ZkWyE%L45U9o{2SI zSg_{Vx#sp213;XWr>KSC$suVbIM{hayKkVt zWSXo>-=E+=96Bm}7~?<8(P}=)26dyjz}Qp@=gyS7=VeikIPEUC{|V|r*X{2s$zz-a z9=Xi2!Z_6I%k4}?E4q}Y_jw5mP6vC^@8Ns8c9U%ja6T`~jh`bkNI{h_w`GhS@`~`& zQ~#6rH%^PVDtQNlit~|A_aDdbq8_1`y}51T zxv}S8QWCLdBkFX~WvQBvP>+bbzinp^=AT&c9j6k^zofP8Bb!$;VTpJ2l~Ig`Y}TeK zf5eHg;&FDzypZ@1$~@P~LVZJMv`b8q z33Gj~>+zsH=CS^|=>XOx^O@@a^H$MBzhMf-o9CjW&P!Y;e}d^NNnCfvvelJJMq~Z| zlDj{_N0fwt&#v79t4OfZez2qS5D8V5jZHeJ-!+z+QmK*@9Joey$09#!A>!7?Ak=9l z=WX5;LV<|-hL$k&^JW{r>uQX%4b{159r|1Sv_QcYapkLK(gbnc_p-VB^F_wux2CJK z@Z$ZWTsfR}ANS!h%BNL^`_!W9XG(1#p+wa7&{xzM_-+(>+aXVtCn?M`#X8v&5pA#) z`7GC}tYP_pgb!jqqB>caM~&~Ddoj*jI^27DF|Kij?1OxGuU)Tj_7v_Nd;k1g5B{p> z|K#5__;FMJ)FZxi5&wNUgssS7^I;7~us^*>aCnLXaJF}xiw~iLtpL0>cc6pn&C$ZH z*90huHhM-LAt3eK$+K+x3E&rB!}8opz}?G}d$h|5=yJb$MB_36!pbL~tQVvMX{r=0 zZ$=0IXv@=ks_0;lEqH%8z+p@sKiP0=dm;Y)hm#VI?WfRTRs0J1%i#p#7j}5t?Z$J`K&JWdJzJKzZ_RF<|hEPtGYxN0-6n% zZQGCg(mz|58YWH$iJZ*0Ef?uy@WnHgOPgHjF!At#$QL`%507a+4hspms4T-(sdtbDNh_O7GTYGaGh(X^F#f0dXD>gC@n|b`KO!}efTN$T)H0S5 zP*G*G(^rCkAZ>jXdR zARktnX|vfAc+dCW4%tV*$oe4qTiU@Pr1-POQHs zCrE&eW%jIiT>@Tm*UwnvNPu{&hWw(QvHtrh$t;LBAz4q*spzIb=M@N=`iKUDa!)!|b7NfP8eDSFx@|iDG9JG#r=VV5 z&4qZfz__kJz*A)#2@`b!;_33$WnaeP3$Kmkv_hUf2vF>1i~fXvQBJ`evd9lHYIet00jwX}AC#?rM}u-7`VlezQDRiqp4xQ@QGiECCjLm+n*| zzYXe{M_wTRj>QnH0E~OAPWDrSJ{ss=(iJ+5d7qeCn3#dQ*uvsVIbuaX&n4j(0jmfw zZP<`F7kRaq{q~4@Edg`<17duxj>S(XPdvRXlny(!M2tF1kvAG#&&B-+sFe5AI5mra zfY!n)?dbjb8T5Yboh;0ysW>BfKRCdmbGw*~ik1rrV^nTI>Ci>YIEJO=3rNiu5W=_R!beQ2* hY?_@zht9|A2~9f(crR{zv#7^$O#S(*`i7rz{s;Q9RJZ^D literal 0 HcmV?d00001 diff --git a/sample-result.png b/sample-result.png new file mode 100755 index 0000000000000000000000000000000000000000..9f42b86fe3debccb3429007db72a94e37a76a3b7 GIT binary patch literal 74536 zcmce;1yogg_cpq*8wHh6K?Fpkq#G+YAmy#mg(jv{K zVbgKv+MILV^M3#D-f_n_?zoKcz6W*hwSMatb3V^ApSipq$x56)L3sj&LYcD_@2~^#pF>u6rJf#xAID=lp77`6FC>+%P$=T($p46vL=z3+FZry+ zl&s~=^{j25Tk4`rpIg5$F}F4`)Vlgo*V4+++>DczgO!8%s)4oj3jsE^|NJ|vxurf^ z5VemN3Uw7Fb^p#&yYQI-C&#B_5qo=U4Rj%NAF6KC(U4qy^!3q<(Y>%E=HFsm-|=>y zxG1AbrmTGX_KLSyABo!A3_E$}cQHyAvo1WhWlTw{V6!ky{^;uI`i@7zo5bzx_s=f^ zbvQUDwj{z0Yw7f_zpxye!#fH!jofJ})DxLF6&PZhZWVzmV6M+T1P=p0){o{^2Z{E7 z|90kz(*AFgEtB7nDM6ubdx#&sxc`>|Tmmp>$WP0cR=2~1-(1i5o?#~ZhO#!W*4)2P z^u5@FC-?vI3g!RdB@=aie{&1H_57~GJS@>|k2@j;iDZwAzlwe^{or%~@rN<5I0>;C=yb!8u~X>Z)T@vA#HeK<)^hBp zyx?sV%2i?d#)iRkNAh%MiegTLzUcaQ;*~UWrxg6rqo^zA;hD882W+~hE)A`_@8ZgL zrR79-uylsak-|mGO@cPHS8ONagVo#e%Xds=A_ShfY_>;L_%Sgtwg2#^&9$9Z7v5bT z*_<_ZNqJuBZus*fBRX40ssYbdvT#>;;oCvXN~;8Zsh&+bV0CUFC4^IddOXa$t0P%H zLRb^VcB5C2x-o=P;@Puj=o~$ns?D__7YN@^~v@rU|n&Z7iS!Me%{4)cdx-=`!c zNv;jpMmaBfYfP4H&7I7k=6ex+hFT;wJDUSPKco@xR8cW8r%6DXxBtaGzN($|TG1}| z)0bAd%Xgpg^YahaYFJcX=q<6&IYT9sgvDY%uU2O>%Y43dM{j@ z;=5F)TVp+5OPq%7ll}gqpST}Y+Vr&1Y;P`8iCK*g1!njynf32)MEIGWcCUL!Twm56 zV@9zT?r!43QU&KrmMexI29S`F;v3DK<)fT84ebVS8RhBh+BNsza7u-8=)RMU6iRMu zYx}(Reav!dV`7?tj!H8H(-x0%cd_mKbzF&hhLN#xr6{#92HR7sr?2mkpD#W? zSZagsR7^=vmmR6~>n@y#?E6}K5W+(|T;DjoQ2gN0GaZ3)B8EuK%y$13Kr`4WBtU0RXF&xQg8D9!c@r?3a4Bv;h zLVMfG4`gM#*2?#GK1U(nDcTqh<2LHJVLDLsOi{xrH$gVyE6fgmX&>Cgw5HCHPAY#8 zlR`p5V%;voo1B{ZFh{S^($bbF`d;+Glefdr1uwoJN1~2Cc=9codP|M(C7b0&Ufd@| z^~dnjK>G2>0U~6F1i=9e5wehZmvwl`ibr;>4#g07c>UYORxOW_M3to-OgLSA%n!5qM7XarwR%+)F!B+(@W=C zN^46uY5qJ zF0dG~g61I#Mgf6LG!6NwX$Z^@q@?`o$BgC&h)~w|;BkIsY5#0$3hJ9RTgKSf*eq>0 zvd?yRmtQJX&nQthuPQE1M{%?H zLtPL7sGa+{yH#U3oA~afC}zC;o%EZFg9H^_8qy${H6GjI?C&4`o&HVQi-q5J#zsf6 zYi1}^^wN?g9@m{AW8Bf%$!L39?%C5LF?q8+If?J7MX9pPhAZAef_bO32BBM~KA5e$ z#D3YjS$J)-k$0$LrYC0%kHH5w=LZ9JGf|bB%bXLk3v+@e`R>lw_6xeueJ&wZNBlrcFmJ}#1F&#hyjlUfcQ^W1PMKg74r z%;%NYx=vpacnS*}OsA`>D-{5X*pY8O^dbxb$94V7trl-M19%@t-oHmtry1p4JF{`W zx;kw3+dFU0tt*%G*QJ}sc*2NS+V#C&xeI0rxyyx8rt4RSGiSyz9XB>n-OK ztRo$4DZNcOfSE0Ihjhi%v~0c~CQC|9W!4(-pmt6oI(S4fz|zX<)2C1I0RaIEFQb|h z8^d^wM}PP?T|SFKeKUMVPYsu-hSZ1BWJ-#jsbTwN?rRNr!2mvk#v1rvmkm$i9deXs zQb=uqV^7B&$sHkRn+}T@M|02J;Kv6VqwPW@5X1EE*Gqii$+|`>CV}9m z%nO;!hfT(nvieqShIZy`bDvipVNuGw4Je@;#+GY6^<)6z|1i(;M6;r)X_{Wc=PzHs z-W;2lNGT}z^`k@sE}7*SlxhyC{*cwnnh&jPqlYtVoRK_Sv| zx=%qtMs0CYH<_6`B845EDk&v-kun%Xd5)GpQVA*bzX7q*3=XSH8we2j661kbpp8p0@*$KW%qs%Xq3S zo>IU%$*dHwuv1^Y?M-Dr>i;4PK02LU&)WX#@JZ_?{t1KcuMSgk>iGn*s!=Yg*7{L* z%kU1w=|?#!OifLhV!~lYFMKDz7%oD1dI{zYPoGid*egmRQF4N|a~%*yX}kQc2r2*! zGhJnmkB^t>^}EC?YhH>EqNJoW>2uuS*KZ0>X38je60DAynxcX|JK)}&9mC2eOn9i^ z{f81~un@Hz^@LoSV?hC&D(Noy24lmd!@78;= zQ?KJ$eYM48_VCGzNkv5hB_D%>quMMEqFO?@ zmC}<=$`Vr?N4Yl92RYssw+w2PK1NDH0v|~X2q4*6@DM$f!3GcxQcnAf_K6cGR*Hvq zShdtFEV2+7y1u^tA`Ct_U{QPJoT>Qx%ZyUU`VKm+4R}thPq^&rM_4uS7j4ZIk1g7T zrp@9pkW0c+Z5@~kM;cV?LqjkPk`W&w7e!+6O_E963iak8S6-~gv(*`t?%c_#wcXY! zX~jeejB9K4V@hJ-kZ697pP|%=C^2l``}X!{ z905!>3Pb0YEyEL=9$t}pCu+9I+1uj1{%M`I*0-zud!=~&MoNlr7jCI+yUl=K?!miO z%wV7sM0;gJaoph&cTVN-@8yqA>Bkvt+~_#c#NO*hF9owH|NQwE8ihD)A-JDJ-8hJ0cG8!m%Q&3O^a8W^e`R?Y6Fn?MJgsW`!G&?6m?QKs`3fX0LrKuzW z=4$siin<+Kd5oF!ha?jebk*Kv>Xtu=Y;tWQJu30kCvwDuLH)toc!)|>Z?E9x6Ha%6QQ+=2cj8sv<$QV&u&wkRT0 zdo;0;Yth5SlFoXoBo(OOvoDPj*l4VELrUMY*ar1TACuF&PM)S1Y{o-dTAV%0`6P8i z@)|9F<2sUmS_*F+0K9gUgf_E*z~_IRP%=AgBA#`GydFqgQHa4tQkL(Rh^gxYtsw7} zC1X<1LG+_IPWWyq#D{PtNP#!<52)TfhC-da$U@DXM?=f+wg2~RSeP{&9fR_~c3Okl_;1Hrgd{emk^Xvwya)jvbfBMhW3<@t` zDXvlL*a^V^^k4urJUVm?t!_3aNJuV?N?~1M1#~rkl8H~hiL2?PNaMXRWu=`Pb9&yo z*KN}}7C*hM=Mo|3=x>6JWHP{plNdEP_s}_0+)~MCo0%^f<8U|{6x}^DU(7)JU5Ulz z^{~bf_5&x=bSAq8MmXgq0{$2^?haN*EL=*P_$u-o8(Ch|YancG=B8ci)8&opsWcF_ zRg;Uw|Cw%Z)<1fB*NcSfvOl?-E*0Uaf;7zf0#VKt6#pK12!%MvI|6AfYA@UT9JrTzUC!HRF2))}GRUaE@vT*v}XH3n9Y9I^Bl;@2uq# z#o3b0%%WEG6lh#VdhwjVX=)Kw0OnmF7ITWYhk>Ie|`3@}}$T z=(x{}$KzQT3Q0o;9beJ_vIy7sL?U{emE>e!EDEXs6fM57Y}GDynPRk$T$v#>>Q~ z6Zy*_?A_KUM5}*7k~gv}sDB-ET?kQyr?F!TEHSXT1&a35?f$+x@gv?9E%MqRAHs4p z(>ny3|3(A{lq-u%3rl{ZQ~u6(*YBn9Jp&g&;o6r7^b?J2gJ~1;(2aygz zR=k6Dm3cekjbX?is;cgvrsR?Y71TPA2Uu7)P_j`xhTG)VEz@6_G+JAL@9+cj&WenT zjMMFjeu&xwB#9rR3JF!LOYz>Q3YyuyR?yw&5bE~CSr&vO>6aWd9w@R21@<#nx~2}v zS(;4E+Ml(2+j4fOyu4=>pnLlTMb$6ftVNQ=FK63?*ZO07a`cUNHu1wefQS)M=)}pB zx%$mi7cO2j`cviATju0Y&v-)cSg4Cw)yHK2 zjk`fk4XS$uyL_|Bwk&q6KG>rhp_|?64rv$J)cfOh25g7Yf;F6#sT~*Y*bN}S2W|?0 z_CdgPnb5eNQEvDfUFB@@@vmO55_5%QS^)I>-qKnj-ew*;Nv3XV)6*3oFokuUaGKxA9E~z9U!#;= z#|~#d6JJec!kwU)3bSPis5Ks67FhK3=Lm?Sy5MY8O=q7%38`wY?|CeK`{)_Uw4+$@ ze&rr<+=f*X7kXYNujip!Utd4HHIYWdp#5AVPctnDWnMm`T(MS~8;&Bg{e|*K;QL<~ z3rxF7RguL=jp5kOPeuj?j~pFKAao<@IhSFM^;Fx%wfLeUbU;UnB2JLv2wv!#dGor0 z=-ks&G6I@$vOnHlZ;y_%w#{zET#qkc^67E_VPS7zb)85Y$7qnJ+J8Lfic)&4o?zzo zr>IM5`k@mAJ^r+Vd{zfTh21dPxp%jNR9I8Sv{(8}?wJ%l5Vt%`Nk;nWh`fnYXb=Cx zev#OT;ITTsN`uk@ft)Dq@7lEuUCSFjU7U>(COqZE4y%Cbe$~#*-@6wY;&zV4?S9Ts z57Fx#U82g@5&ICQ*5GItgyZwa4t)#>p~w;erR>6m#=R_!$B$c09Lk3>XeHh!nGcn* zeMu1*f=FFX(N~dvqyN6Ze+p8rIxaMPU?>7aN%O(MYvSz7XlN7g9 zD7HRYr0T!t6{O3pOqRCO-JL(quiqFNk0_qiUzTYg3o4+|XuG<)x-3&-2m)2xXU6yW z)aQfsLO>`fC@7Zfgk+jFpD8Ht@C%X$eEQUBvIV>eG6b9ut7iMH*?insutpd4?!r;f z&YEAkkTEMoPft%TJW*Gt!jluBCXJ5ioB9#X?%kt7_Tqw>pfT?grvOhp95PCNOPw4) z#v3;b+R6$FlxW3$ki|9)b8(cd&*n1jV%}NJFYiPI|Nj2|q%9mC^9c@37CiO9A?D!i zpyF!!+kAqi`^41%;le_&wWN?KLx}u2TJfal=-Y^90tB4fY13NDSAXTqWxnW%C>Ql7 zBLmh@sYtM29yjd88s21K>HPVLiPk-@8Hy-CN3y6elQe2pjog^2_~4aZ;|xQ9K*KzU z8V@4(3P2GcB=n9v9H4Upsyd>91ZcAHZyc`)PyGF|(m2KZHzeA;rxO0UHD9jbapp?f z!f+*iIwiyCn`@0==oJBLN#IDt`OS@uohcb*Uj^q1?#YR4#AWB?;MYg|aGSUc;(Q8z z%QytEMD5JzZ-VHW4q7J?0)I9)6Xzcv2*9jNww~L{G3vZ26T%Tcdxjk4Dhv9|U$n$p zGXz4WbO55~Kwjo5o&}B9INPiv=`jYD&@^gULUhyLc_o%Qw?LRHnGv$Z6@knK8j$Q; z_*>j`)rp&2xkuMYF2icrdnpfkizcm8Z>K$gzB&on=x_YMq=@5hw%u{z0Cz~z9` zvY8b-#-2V?RFn?D;HT-OdOrOBA57^8Du(`R?E+=QW2hS7^Eaz>n zqz+wT3vdrib<`s*#}a?^1l?txCk%805hhr5#h0><#Ibme5g_H>-X|&3COA}NEov|M zbS=9$cNmnWvufH=xfa~C!mEBqf>`Q}d!zUC$*e3buXVAh7CaIY6JrrE>q=F|K&cZ8 zlY6~SHT?$$i4Pv{6Tf(zf5;_Y&tZ4PZH#AY@MKMh@8lhG;QSU&4daak+Ir>f1G#!Q zWp%$Ej5}vYMq)^L@!e@mNZN!nyTjwB_vCV-^mpm3au2v?uz4fM0OVcyeT>Hp#foyZ zfS>jC$kA0dx7Ax*?KAJLXDd%dV*hk^hK8M!fjG&@{eVL($C4tZEggJb>n8*nuTY z6U&IagG87-)S( z(fcAx+XtF^Tsl9_vY@e&5G_N6?3WCwVYj%FBGEid$=^R_pCxCf-UvP`P9d`A2PGv| zgk)J)Z9Yp;K(wv++Z+wPk;LUYXU6^R&9%$o~K^ zbkYdM_u>+t`@T?HCw&YyTG4?-Hk$j^`%ur9sBw95Y@ew&ZyIIwoq;Pvz!ly-7ev1$+>fpItbG&f0hH_9 zb$8^1{l)A4u)cwz`xl)2VsEV(@#=oQhL;&+95p zG_mR}99x)*=ckVhvE8XGKNXM~t)h@Q#wJp1I;oy&$4G9T5X9C9`qk3fEjYBm$D4me z@zn#uDIPQt_G9~QAl&Xlyl`vVpPFl%EJe2{a=4rGgF#+o-N?z%9jhBg1fgjE0m+d$ zJnMB%X`JA3m|^!4giS+t-GcWosrwy=lwrs5hsZ2`{x>epW)4ohim+!C4I*`)4kx%3 zGHdDgAN5TO9}@U(pHI)YSd6fz0V%9w9C=|Vtd`yWe>mSzp>yWjUel@-5# z2*%|rN{P`mQArpV2sKm!Cg);Z$9T{gU%tci^{vbwuswChxZqUT9ht5Z&q-{nZhScZ z@srN@!H##l>xyoDGa6STKOTw2k8(IVQsro5yKBhMV;U7!wi@*~+7H~xq3E`4YZ$VB z?k+K@M3jS9J%T%Iyi;Q&rA4y#@;0Fyfz4GZ zL{P@>6SBmO|HHz73UC$Td?;Xd#z_StNJ7|QH7PUm5mfX$R%iR<3gM3H4-7m@2HH+I zGVl^04g-laLacv6;%U8O`u`; zZn95L@N32FtwmUf_ZF7dIEe4mC-_g>d z*>`;g(3KXCDdMv(-h@UZG!t|9UM^5da0 zClwP?Op|;vlN3lpUo$d@_f+NONq4RS+|*Ldrj9-*kiUtam&)j7%ILn7HTSeUuvmbw zC8QDcffg8K8*5-eaoNk0FY2>HQ@y)(aX4-Jrmilv61yxgl2k~`%F241!sW8H!iQ)W z7#RLIw3{V47{i55?(Z$8@o%uwm6OL(eSA*Y3q@Wp^#Il?gQ2*Wot+(ON)6=%FnaXm z9T%QR?)Vf3RRx2=68qsOVzEz&*Zr-?oR1-l5cb16~QF}5%QHt0O zpoSjYVaFSRU@#bXD356<5U^ZuK`lgB(HVm9g&-_^%RS<$P~-GgX|~ch;Jh*_c_v=S z$uNjy`UPRi{B&7h7@_8~OehlJCcGJePun8BdGj z#7cdGDup5y&PYw_kd*4}?NM%1#o`~HJk&U7+pr%u@8kqpsr(%{aZg1i$V?+8v>GJ< z4H?!*KzeqSy*`;u{L8(P7m`FV0d)N)7swPSkfX$qr`alWi@3GNJDt0=rJluHsqcNN zoHy-$4z(YvU(Z480isxZkx<~C3+pm<2W!2vtNZZxOOxMT)e$k)HB3l5%0>AL1dQT0 zeLDlNJizVvzNZ7U{3J*Rffx2jSc7jiMP;r$~43KM}1Yn9CkbDj~ZuX3E|`Sy`1 zaojU1UWaFub!e;kV9mRwS?DJ#aPAIH?B_LpAtd3PegAKKtNJiv)Z>u96S`nC*{=>L1yT zkiA`56mR{!tn_%?1@1;bR{xYL{yU-m_IB5y=b|G|#w%(B`3uk-91-O~&oj%t%R=7? zncdwj7>;YjFdS5e{d6a|2czdaR^e!7bZ$k$?>ELn2uIe}XS2g=OfOJOGx~46RH4TEE_R?JJhr$VFWx8{BjOef;<3 z`e#ew^U<$!R@UsPLRETigge2p4;8TaiueUSevCs%FW7Lv8FX%Tb$%!#CB5i@FlObE zxVXDsbL2oAh@h^qgktar;i*h=yN{)oZ;!E^IPw}Z9-_bEYwmatv3-JPC;<}rq`w4} zKi%O2L~8b9p9Gd0ICFLDFJ2MGJk9U7^k%MaBb5_s^VVPXM#dYz{+H-pAmq6%)StXY2!v=mWyoN$TZ@ae)HpAO$QIBWD{(+z5t*r# zRajW<`0?XklagLwMc@QV?vrv>8ll5ck(2R!_P!lIm%7N4Nta&sTg2kad_=WUyRtN2 zWuXIBoOd*1@w=E9z!z`{CgX1j{*M>g;$Ny z#!k(H9W|fevrs3II3=-$3EN9&zIn3^ir)GvmO$$%*|qMkPrN=#Yy{e8(Xjg^Yk&yi zSerss8QZ$-P&L22u~j8YhlcI4qq=i@y!`1077HY7RPBQ>L>_N!-4&!CgjH&uBeJsSo_xey{xKQS?!2$44OdJOR8y&n=6d2IVR^xVz{I z>~YV>D*bOi?A!CeLURE-dLl3U+e}`>dt^R9%b5ll8_Xf+oInX4HdR(u4u}JJ#OPP% zbHt=lNi3SmR+`JT))@`vCTB*LhxCEG}xbnAj|@E(Wmnx zqdb@H&yPv7n``qKLwmbxra7eW*+v6zNULx#0$R5)U^-6D%)I4@`%zJG*9nX_d45^U z(x5tkt54y{lK{l_Tk5bTj!n+U=mwu&I!qvbvh^O~M->6j=fWz;Uiak_WuHOfhr$(? z-GU|y=oTVT88YulQLe0j7{415M5({wVy_!`opfdgOI5;nOxs{8unucx0Bm|i^y&iW zF*8}USg;w){>r)sQ;2vOtpS<72m_a*5ikdI@yoP4lb#r`@fugWIey~Mfq=MwMGGvL z;&zsT%R3PFj>BAGqKqFI$~FDl1)s9tD+Hlbbze`QD1Ni>npmU-xxZ;257?!o10tPP zl`ot>4+gM@;M{o`SgRBikNA8B9p*d1Sk54YF2&@6uLhAG2441=g8beE(jy+Hv4;j% zPBPsTSTKFT^FmM^0YFoVI+xV=QByk)-g$ilG3bJ`XuhM{B==d1_b22ER$$q+`mU{tlV3@0}txbuGyAP@CBOM|nBe^YzS((H+9l`Jc z(yWiw)gu+s(DeC&o4+{;e5g-nr?r=AsUS4n{AASiu!z?&y-_+?b>+LcS!iKl8-NCF z4h&5kX!9ZVocRkbYG5tWOOYbV*g1A2u~RpZJdJFqy+L0|>Ek&X*^}4FOY737F%0A> zW9GM_-oy-uc?iyMWeDog1q+722qmuW8-*b8a?%n3|{ zTfOHo`E@nfLTMZPT)76V*PxPkhmhAFhQxpcTP9LSIp8D6c**jK(@lc&9eca;djv}$ z*x+rz`R6`gK-f8;QG+P}GwLs~0?VEPEVFt80OhwK3-8H1#I_r~ysz=Y-Kh2SSkqbwL^~*c{Mg+=qU&;NyxoBG z*1WTsPE-r~+d-znl(|oWg<~KEGzpCq04G>UQ`eJ-i|AYCjKzyzcY&>=eju1RW~6tP zbBnif&&9Ez6;qtD1=-k82D2i5W-%(My9wvgGj31UwQrKOJ_PW8^3DH1>L<{ahgD?3 zLjZ~8zw6GH%D(?00*J9|A9xlZ^8ZyZL{RQyeS&wEASdl3S2j&4v#>)KuSNh>GP_Q& zglQALnWy@W<3RndNVt*psedbW`v1#~{{Hy)|5oDuDwMO<2h@+BiW%Sm*qMHVaD5n= zz+nyN5axWgR>`z-I}Y$D)J3&p8u#RYKe=XH*`4242&Hy6`3eX0M4G2}#zlCEyA#M$ z=%w_l)_`7URse4}QXP^OJSWQo9QEs|%e?4}{-IifQls*oHGQ>jCgZ85w-nia6T_Dn zG5@Ah1mfgbb&PqrXqC|p8Xq=aDXhf%0^Y=xEH?#0L&w}eF>M}H)-@E>Y(7D$=?i3T z@Xo40W^Pz(dcV<#hFI3v`Plb2baf zgG$ia8Y^K_&yxQmn^}0Viqy1R8<4*W!NAaqw*4u&kJ(zyNC#>vM|jQplm7hCO4BG; zzr<_i^Gg9v4bs^J+6dAyKqtI88H@A*Q9+FzvCB=BAuE0pSrNl~11z)v0}=Cckqma= zniM$m&Ph8t?HG28IId?ROh99Av-^RpnLw-Q{}yp4szFGI<|jO2&i>BZSA`@D4Go8N z7ynczKI?+crw|UY_C&dEP-Q|*7nUahKQh^rd`9FR+} zP+{?$g(L=Lbpy*TD6)KG|7%Qy3kpj=x%G2Z;m403K&!+%ULrn7D!boq3(W{lF0dHE znhgjLAP(Yl{r#GF8iF<0fHUH1q@OPKAz(8~9+S^t%)l>H663)6(4uBL5D%q2o1euvIAl#Uge5&W zz)efg<}=ACubJm#N&RYL71iyWc-Ykgp8D@tbf3>r<);;gC=fH_H1zi>tF!UU5v-NYNgCL*#?3c0(?4- z^G|t03}qZNPYtTq&08zp)QM@uoIZ^RH%{ap^%T2Uqe;WO$UgH`*{i4P;6w)vFP0_X zLLAEvO*eL)@WQjcaWxvzxv5Si+OvHnytjgEu?`aM6MNO7R*ZwtPAsowOE6JR#zf&D zuzJDvb|`b&e=sH7B7+{TyIb+C%4Lho?q}fAB#7_k@|F==m^_67|(4> z^M>;;N`9x#?7c_yFqJ9H+A5*6_-(&@;Ltk$;`huo2*74qAtUd6MsMtIqca=?+V9co zHvVmz^{*)TJ$lwOMr~g^+0Q|LNjQgD>gL4B@3=$X+kZ#a)z66C?J9|tCIulTdtO7| zrVWSF2t}YX%n%~>$9LdPU3N&N+XiSgw{5}xB&9abQtNs?!C`69>MVgO%m}Nmu)Xe_aP7)PLQs?KD5(GfxOrS3M_4J2gETbH!1?$ zLP8gm*wIk{&U2a-H8e7yB~y+ySQ;F2P~nA1$%fE(IPon4)fh~r08M@K;K=Bj(FTd4EnXVwcbG21rME+UITwRBiV(7YV%(MLPX}de zWAM))J)<_@UzpxrZtjL2vP-4lG{?;5c3+D$AF1&LWYR_JC)74xdlJ0N=a3Wi8@dva z6iuKv5)@$vQ?S^lFx!aRAY*HHtvmzZHPTABPi?@eCW1=EJrEm5`kZF=bK|0{oY4M|cKbOMNCeD0z5{8QPFtsGO^laSq z>8wL`o&~3PH)~Kx-cqWc$-u;FJ|eNDF=0}>NlVe8x$t_~v)CgGJv;o(>uD^$jRZ=L z4!kUoDMy@NYn%vpGSW~Msh&S%P)h8l4{c4GYhR}Cn?x1M8ZJ0MadNAn^CV?zr_!wl zA}8|Q-GX&*(#rhYln`x*$*{QP?5+c4#4`HFuf5Ntemp(-|q70d)j1+*WV1Vdkmbik-3 z_mx&nWrZBvUTb%L?b2r#ljSN>xupF(VbJ)+5PC;|N}geZ)R=|=PqW(HtK?31NRh%q zQ6|Bs`#4c}6$37SBIvv^?K^du?-8`vzK6cwew{w>#0-wkJcqUkD8_So{#}eeC3>BS zNgDb~;NYFpCkHmfZ8<6_)rPooz$;gW!;1^;LEb}pM-h`nV;AC(r8)?QTU+^#vO+yk zgPDr0ID)5f`yR#c_PK9A|4uPLoe=sPW?q&&4-69;EjuBgq6j+QRp(phjkRxcFD^X?*ZMhbIQVfy=g>{{vn zJq1ar6m5BQ0dTy=cTW44sSNLR6Ku zfQ+J~Qqv73@hR3HoyvE_?2Mw3u7+e`j(`-8LgJ=3IyuyB z6pToVVPV5A|A(MQ`q;r1Y=TuoZv6i|t-PXSnkxJ+DrVM{I=ysa@s!(ooFzm0>FNtQ z5hA7dp*|tMs{OXZhTUpNQvy*h%fwco4=M&)Il#xoNK1PUj6!NvtAGe_UsKYbK$i)G z!Y^PQt)E!Vof1J{B`ZNe_y{dPP-S8`e@wjaP3EloF9e)IFH3ujnD@$S>h0Um&? zXDj{ifug0KOg11t8-aOD0JD;0RycMMY(Wf>b}@~sFb-=gP{^UJh@6n>f-YcndQ`PK z6mrh3?>Hmf{Lq7(2@Ttn>U(~qVPs@vIA~^3w6`w=v4&m>*ke)x&|XLuABHz!F=Zfc zA|)b1^Tdqu_etqXc3|?dfd>1+0ka5~-3( z)kG6d4;vtd&q#CIV;^W&eKc>F+0=jED}lj6v1Os5NzI5#Dmfm^B?j2~` zy)Fg%bQ|LLLVBcIC8&}7u`=vI{Z%nT9kl945K9l0L+jfENlEZ9RWo*r1~<96mSfWPGa47O;#-2H^{>)8{H>fvvCMy=0xpFFavboforn)lwDo}qHr;Y+lK zPwh3H<>L2pFdf;8NcI8aGhmNeWlt~}6w-UWMSkK3f&s+x>1-_ya2_(^Q+4%0#teyB zAG^lJ5pG_;?*7y~=Tb0(0=&59{yP5?IIAI0vLzIEDg||p7H7qCdPCEGn=kL6!ImSU zX~75Sq(JBNwrc{-&=4SZKpZrrg1jXBC+j~9z|_O11HTX&R%>iVl~JgyMo*i2o%-!6 z!cw4F`}^$66peCcQ*{Pz?oY&zD*rVa$I|w;w%$?uut4bLgwaSual#1#*XZLuNI#Gj zzJg3irf&CzF8RxsJ5i2HUPn%mUzI{C&LEd2K|5A5cw7T08CdX`QgDY^S(mqRa6IX* zCJe;BWL(eYtxv~AWS+zbj+)62;QT^7oQNz8!VOX{rgC0?F4cCLxBsfto?oehpZC|M zy`#I7pPwJO@aR_IjT`DD=Q8f4XbxJ6EoP%bKZ+M6U-CA)ftO~Nv6Sve)onbfp`zH; z|8YPphlPp#BzFzT;cSxu=JWSFZxLnOx)n7NHAKFD<%bPJ4B0A-vviLe*JlV(HJuOXg1x-oob}TvQ<?(nh}QrKO;PrN-)SXg*QFmL>8a?=J( z;4PP3$3sVsRCIP;D8=9_;B_mHpFQimQQPwogzR_UzTGz9qBO`wUlx4%kOK_Vr0BGZ zl)gtGtKNrZ6Xk~(jKHhcufHgyoV$G42jV#|pSf|{`3WIL#^c073n(xXx5P`6J&Y2$ z0!1~eNega889b*t1N<|L1R#XrLe^Ps#Z>>%TFT_bEb#y`9PLJyc+7Q%q}GSwe2?F} zdE*IZ$sbhWyLazaL4AZSUpe0_;)RtJDa5Kj;3~X!^X7x~%}oGsC$t|veCRXTo0|f( z`fQ?H)J5>JUU+@vGzrMOdRq>V8}94q#G$UjvA^u>(#Va>vPQ_TPCqo8U4_{#q{*%Cp&lzM>QlUF-Y$)v!w-gi!s&f^W}T$SmSQl zs_JU59R22ZNZ)wU)0>O%x|4eHKo~S^xN2yFu7-c01!FPO#|QCz7JsUQCro8Rr^Dr- zksGNAI6|xBHGp6bK~Oz(=#VFPqIfsidE*lVhHT2#|C|I7hnOAIhI3?O-3GN%9UV%v zb$;AWMznnC*t#T^w0-Y45*^K?tb_3#Mg?k!M94CrJfVeB^kP0c`KgN9T8UJZ{H9+o zw^mJI8_|HF_|4h7(6@B_(*t^VK4Tw7X6EzY#rXpjFRq>;sQ-#Rdh}r(EYW@UExIl@ zZA2hK@PUH|FJY7#e`Yx2na+`tYPVfHf1VB;m}cxyJ&Pi{bLY-=Q{*i?j zzr%>zquIDiGS&aGFLywR5;Z(53HQ;Gx0(Kv)&h2RcA8M2tH~R#^r)!YhY`sf>)X>cBz6D)!G6V+J3LvnHP+2Pa;&Hbll^Z|7Z5)S&hwy-h2;E%Xlr{Y# zav*d!JV!i6;PN*kI_Qg@Y>W3+b_ixwd!(wWI?sdpJv$4;7yeIGQC0O6fDmhV1=1in z%^@gtczpg${>I&%XI%;XEt0w_^euV1g38!<=+lKvPRyY}gGAl#X-mlYLJ zLTM_1-DayKxFk&sHhb}LnzPi@et==l($LVXVfY%=Vg8|5_yKCvwu%%)KTM)O4|Mri?G>N=BAowcT_{I0xsAjoO3el022;H#CoJ_vj3@T#RO z=+`OU9C%XaRK})O^ak1=W(G?yQTis_abhrfz%slX${oPiCHcOZQc6A555k34XlQn_ zU`rdbdI`Co@akushX5lNpiSZrEHjtN8k7@5d^x%o7Z-^h36jj>aK9RQUohHPoGu)jUA~w_K;Z0f5okqSmaI3w5e>E8_p#XQ=iRuA{Leu@tgqRV4B`FQi4y^;T$HC_K#$bb-$R|cqF@#h<9mQOIm9ah zZf$^<`t=03;A>2=Xch=%QV$gvPjw60Mh)NK$iODA}l znH!jQL7+T#`t;ie@c}-`eG{dlrYsQR4(2MBLx&b(I3ETt*|AfntgC3ajb&`gYJD$x zKxYE^UUO4pBlj+QL|9mTg%8??@YG)1PqhE_(Bj%6G-YXFd=@d}@bavdDeJq5?OyP% zkI$v)1%Pdk>S0YrcbW==kdROnyjiLWUfESO;wNghXj2XC8>Ak@)EAL_0Po3pRGQ&= z@vi6g_BPL=`xTndxp1<@bNL)rQZ?5s>N%DQj z6q0gs`Ap`<-RCDt%AzI&U2M@Q37Zt35wafYr5J7nMK4biSdJKl->`)Z z;7;Qddrg*~&bUei7x%r3jn)BFz zKN)!*0$%v>L-Nq!!_T3S#11?yu(3>f{m^MVJ;di8%l~0{bD!nmKpFqsOG^40%8%rI%S%fO@Dehujis@L)zun}-6gVyC{b#F zD>Wcppin?4?%Dlcw7m&9m+RU$e3wF{q9wB=n$1E|lw^pc$(UKnlqpljLJ<;5q>@HL zrcjwHQ%RF4GKVr`3d!{SF08$uy`N`)$M?S9`yR(W*0Y{g{_g*MU-x;P=Wn`Px4$^^ z<2>N|qPK6CfnO+rA#nhezYQ&^i#EaPR#MMk-M2&U8_X@k((F0LO`=Z7Lj}07_weCW zP(6~Nur?o$;1X?$Dx~=m3^Gmb@inPg2se&RbfI)npHoV=x z+j(Z_G0o*d+`U3Z32#Mmm|*N%7jZU*s2xsQ&|`(+Lr{yc&l$ruJ#tzc7R{Yta9t7F z6c!eC4wt$+rM4YhL_kPLF&Tyt3&aWJYcH_S!7{#3ih@~i5OS7wm=O2k6~$QRkCoaE zF=i1~%^(6;ba$12$(j9FwD4)rd2g4Vvh^y2SCf&M*$+r7m6@%oeq-Abm@Ebc1}G|b zYP2Cr{NxI(+f+1cP-L0_`Ct<_XONx#ZMyS!i?8G@D?Jaq85&!zeMexXujL-A=2|>Y zBghZ2`@hC#qsa{cWaeQPc(XN51$U%0TDS z2l?yoO73Pn6qF9`Wx;N8=gfOvxbzEc;|L+cSXv*3{gX#71gW#PpKx?P3gaU{WlWT&@6gJ`6jTOdI756=R&S@l8n~e5;DanHJ zO}DoPjur)?N}Sk(2V5g}-9LUR7b3cx!5+e@Q?q+%pz8kptEkr)yF+Z8GE1#e^2}Gs zK#{j#H|U=Cbq{NMiXr9PND{2hP8`!SS&#M5WYhZZzFB)0~M&N=a&&b3jr2isdV8AZXGd{%-b35zV zJOd0%_Txh*(mci;7VO>$vR|IQ!K&zm=@+m|q9yfX642UGjG!^p^F(>Pe#cR3F%}wq zx*Zw0feHXd#@S$cMOIC2roxWQtgj3tXp6;@+Un2D*DatU>i%er^9^{10PlE&`ay>qSIr)bhl| z$)1KV_l5wwGX9c&0_(E~BmTQld)z_Od>Nw0AZ@A60drMPRfS|8Ku5Zp3=Iu~`V0&V z%*h~}eY#Q1%M+V%P{YLdID2h#A}IlwQ_J4H+h<|18GPyGmKLLSFHaATiYCQ6ly5+1 z<+;!=*FCf6LngyHY|O51ZpP3?mQ_|tmfU~vz!=@(@3)52=TYTc&z@~7-x#X^MRvJ; zwM(x{mn}9Hypt+1y2&M26n)$Iu=;HxGT)z2eRBNy!Uuk#6(PQO2yi9i@|E#O)THsbt{FcCg^8uC@vOD32 zN{XJm=pQCR>2zCv`2|>h8IxI6*$P#E&>hD6zi?+PTZ8+!Y0DNds{)U4F^a8uovSdm zc)t%HG#EIxl_k{;3=O@6a)X}NPsJURBxGLehqqkY|GKoa9TW9$>_KF#J^tbhIB7Kz zgIk!O#L>Fhx4v9_JlB;8g5Wwhp9i}Q^u4=fMlWzlVLe8tk3P2!yndY(pGN9;%v06n zXB!h6b6ogDL>NI2N6C8bfMESHt|G^v=F^9>4Dkg*_6`p6kJNAdv}nc)6Bifv!7}qZ z$>?Z!jJtsSnOnq&0@M?PJ!`F$6f05#wKL3@8@TsLm`+0p-%iufL(dBkqT)tyfIu6? za)^_N;kl$8UTfd4`jTnI3KXmMREp{vsnc9f3=1wJ-k?;#w(+bUq`piT<6eW186dk2 z@D+1Wyz&r;ih|_}g;s}#LSNB>E)Cxy=Jhm2tdfLZSQOt}BW14(vhO;(&8J>R)}=nj zCFRw`eT_{m*c`;@Vb!%B)}!Zqi~$dr-qt3s37d!@cZigBpFR6x@SR^s(V@&cnf;0K zafL7k3_OyXTEZ>kysf0F`nXH28PGwfXlSj()};s7^732yB6KB=!~B%dnd_gP~ivcQ96Vldy0pii4>8s5KhfA~a}Nq^|+Sp7LnW zKkca)2nJP#2MwF+M)t zM>El1{gJpq?~7rbL(3(=bYrJjnU(O6iU_bx;`*&I%KBi1r%wCqhe@6R)fZvku zNdPj4tv;j}%HcKmY^^(kD2rC0F;{Z2WIq~NoC|;2Z?HQb(4OxEm8(e?3a?|( zuU}{2FY(^9Z{K1M4-XNDR+Drymu|VX6hDdv8XSC^Qr=P(j3G5&{yg`Pg6VJ)Q+Y|ods#BHBZqS-n1Mn(Ls7Wp zu(&@Odrl%JJP4@e{rmS@NSBYu;Xq&j?Rb`jaO8bU-FcYM2#|d90oTT>1q7!Ao!KA? z6(sYJ zZNZdeXGOFxrhp1FqescxT?f`;2t5xiGdhGj2|0Tes`$6l1O;?E5mXBG7B|nj`w1fx z(;-e$vY1&|ObUG#*kSEpx){J>F@b{J#$jM7B|HA@FM{WC<*vIc86rN%R>XfHTnOoW zF~dOZ!buh&?)4ysmRF69c?mgaau4|7MZ}`kf%k_Q^O=8aZ0w-)ro&iwM>a7p&FBEE z5UO(gzflUE#*ByfqQy=uMPGbKQ-wjGzFOE?xMmp>WldBDQV3F^+|J2D zL?;(U*WkCegbjcf1zfr$f9A|K5RQ!81AGn z6ql?Sgj%6ef<#XF#A^9_AIYqUJy=v{1~CYx7^l$5)1gJc&F2o!-aH3fn5|}K4@45t zm)~!k^qhWtS}ZFoTT=Yx%a>i;EfSG3B zfp6cw86k`ofVcvh%HrSX?9|YMBdkV2U3~@e;?AAC?M90NR}zIW`Yx_YJNvYLp6&H) zDcd8h6v00->Gxh>IG2(lj7M=173Q_sE%irx6%_nb8vK_0bp_1CNmlZGe|~|uDh2U!(u&}E@N zcUEqC_!4tqk&#FXfKsqCzVdOWx_>faIq~}yV(VxDQ?>C?Nx=#y3}R4jf2m5zk>jkk zbjtG&-{s-ynFdf5gCz!vw&d)lmnhu>r~+isQ9jrXT;1Jcth;9ChBn%+jXm0YSgA)M zG8EG@Wgquk7nC@B?AU6wN&a{*W1eEIRk`nZlE2IqIjKGdDHr!YNUtw&FK;-y-RHzb zI67W6Hy0%2pe^@Py)G1oDi?C29j`BaZ_|R}^%9z8@LJ=8fkLrZedHXPgMk%(s?({ir>APuW*zj;ptOm^xDq)R}%Xd z8emWdI^PNmI@24#aDpB6b_HTIC?7h}CHui|oLTf}&z|&e2}c@>_SJZtb#bfYqmCwa z*S(2<7t~w9)4%bqyBitT_F_TfMT1knILd4O&_B!OC~sQk5ZNY}srPx6 z&>{zX1*hhn>yK(|G2DvHeNBeWypFOfjHDwoDrk1@ z&#R!(UK$;GaP#oW_+T1sr*8tA>RsaXyN?XNZOzHac{_Sz;^s}Rh{Q{{!yxbqak(Kj zfKKFHU7ra|S@huI99$f;){u7c@^Ojq%s=$0DT>{5D16_GwOY-KccQx#78d?dDY5qN zH$g>gccIT8x*zvzbR1BIZ_0T${dV&5lEu^y&L_~vhvMhwzmuBE2&iMB#>9mpB%a?k@qMmLzAbvxbDp? zFwO#Q-WAwKNFV+s3YP9-8#F)zK-L=W_syZ?nLoz0#%_>>K)_D^AN!X7x;wp4O_;&C zPgq8V6Lm+zi<^gMB^4#(A_O-6(9xr2L-n2DEy$(_Ox_RNH3Uwfm?U`+{V!JHck!Yr zMt|++4(s7W@dr%DD=4@edt*kTE$)~?Rw}Nrq)v+uW){93J2=oI0^oUtRgW3j2h;-ADe6q4@up#)N0Y?rc8T8v_6G4oP_Ug13EL*lLEG8xhkk%ap!%BA82+h&g z;at0DtSX{JOS;6h!*ATM92!Ac4Fx07)a`{J|99P6zq4yojH_4sqZ>8C z*#c;><&~A`;qeK;EvNDMvM(RRGuKX6M!cTI1TpC9x(JNN7&W&f3h%Ee^k5|^DEI-u zMhQ^;8;KP*ZyJuz-%s|n8tsL%aPi&9N?p&Rtfnb-b=t2>N_3CAjO~-9NVfw7N56j> zT2)@R)(?@n@BDKQW=K64ZAZ*2>wj>sPNX0#Xa5>bBq7c%r&F zFSZ6BUM5Xzg$?g+WKraw0EuwZa#ztM@b`EMlqI&ZGP0vVLF7#kf` z!~j}sb15iDumpk#1b8s#=jV@0YAMK2 zAaJ-R!FcY91!}ta>QZm_e&%0kkFI!w;WFNe`lRCKMxo62Sc|+uLMzB*s(3H$c`T+k zU}7RdxRH{!s@QULk1kvMYn5td{EYa`-uVm0TbHKjuv3RFb%)DOhLa~xW^4_F2F_fI zt_OM5+QC{|u3?CrWz68B(l~s0VY22jRIvw#0|7~($5{gG0UHU6<3^OJ+}y9vFirlr zeK`itcYWeH^05Ey)-=g~Z)?M+1c6T{D0Imv{zWZD@?D!G3hQqiR$}wX`TGZ6OfB*c z4zruW9G8kj`iLl?K~y0&UOPe_c}@Z~hYElZZyzl#JLq|}r9_ow<`iA`ZraTKH6GQc zhdL8~rF#72+;}~A#_gzE*^L#|ec#@l)F}7Zn_ZzLeDkr{0lJO@ZlUixZd8nK`8MTt zi0!!9yGM`XBL%O=-F}}`VwOyq+Sk}6h26-`m+=<1lC55@^SOMAY1#DkJu5 zd3gL89qN=B(|aqt&SJSW`<>bKyaMo;hsRE+VwgbH0)?<{&gX`_H6So>frk6Ol=*5- zy?*`L2!sNl@=}bYai4#rLOu>mtn}IY|5(uJsAtcg^YZh*zRv^n;^Ot|*L7sDV57z8 z7c;E#Zrz{E*N}$2ex0RsQniwk=(p{7yI{U;g{zcI)aW@GSFJh^wBat08_EaxDl;=P zTRRO6jXe~y0{SnQVaque8_OeUQ}2hAI5txN`tafzBprXgDp508KE)sp(O4I73`e=A zD++hA`s$m9)(=*l+hU-9hQ3A_eF9H4|9`!^1rZSAJ;b zwEBK49}e0g>}E^%tT8Z~amP|y-ZKK@w!BxRV&6(c50`r^i+flx3bTgY#Q>(e; zW4jetv47`dPqWnQ$6pc$4;3v~1fu{ljl*0|S z?-lNj6-{Ivs&pkRprvI^_q}>HaQkx#uJG6FvrOH??KG-ha}1t}#LFcTcc)M8tQxSt zfzPcf(M}AtG}Y#lUN!%r3zo9+KUs^XLq0HM{?W|;@Yy(pYv06EwQu4%9%?X2mVF#z z9P&+h{=d?0FY@{GsEkv0LZ))ICI0unAR&K_BXDcL;T~LlfG(fa$RT|O=@oOo_we$Y z8hrNtX&1CeI~^aAMmE`0MF;*_;xddsPecCw$ypAQUFV-yd_y`f2d~hzHJ=K6E<3+f z%{TnsV37V`A$8QPy6ncm)+g5`V)q{rw32-;Z1COo?L@QbG5_3xo@%X*DSfH5H9X13on-p*GV)JZ6G!0vL2K0hu=^aXW=?5f^`%5keOL=E!HdNvVR;g(R0v?Kpun#pd^7H<2W_P1ScHP2 z;yN+0Pd~i~njZyOD1;V*mGy@=F+X_wW?t_=~w zdEDK!My8Au@M{-f1|%&6lrh#q=8!muXC}6w0tFF&I71$!Q1c3dmSONK{mmP^0LcMG-9=)Ubs2N=(c$ zL&Ij;v2*z5-%n%9zu@k=to47Vk@AWdfl`mOu!JlGULrjG>0xcG6%wkfymxBn--XL}Rpo_0sE>z zl_~^I6n5jrrMc~;F7rL41d}-0YA-}^(7=cwMv?Hm1%Ikjdh1u>u7wHWwa*kCtyae2 zWjo*)4MkQKMP~a!Xu*vk>tKl$yb&ILo(x1AH_}kBNngY|vi|z)<945GfHi^m@$vF5 zpa^h*)&zga6sPS`mUPddBDw&eo`ym*eGbadwJ=>hnU<1lk|�_#$QzWr3G(s!MZV zh@h3H&sZ-Zv6{Mqt0!u~ratH*Fi=>18YKuWygF#m)hWQ&7ksOBkawV_jLZo%px$Nw0m7MHs zDbig2EWimpCLaHo=k;tWZcQ8q2#QL_#LT=1ol;d*RSO)rMfmMPZXOPh@QOcji5Sy6 zlM~^o^^x=Dd)55o*YeMfhcEd4-Sa?=$F3a)*Cb-I6`R_{>!e@Zi}QNXxTnU0_05j* z8!2k-8?0pMuO_Sg9_@7GIq^;Ic>H9?&t7dXhBD?wFy z9FmY62Jm^2i;IhC&I2h4m!Dq=31*&=j7M2vu2BsgUK(4<(JieV9qZX0yJTmYXC}N! zVI?Q}P{0~#>;~KR@vmPmk+?vha(p8`S%r4Hjy|cU^{e!{@4etCZ)|KdOWWqt@sIvC zg7mhSeD!+wYHLU7OoB0ItD93zUMn>k9(J2Gg&C)P1VFSg`= zSf=`4n}t!LrvO4IzjrSa3(Ezxe)sO&p-}B8%jUH(69hmnx`4t71O#3Lzher20&qpb z{tfn>rx>B;L1?WpCMXR290)%4Aw(N_-=omiw_!mLD2+_)z&LU-jk)4@!QO6f)Am`s zPNi)QNLBkX%pd?W+y9k*lHIh)wS74O&cPOJ+Ps-`p|^3+BjM7WM+(yBL{a?5^xR7J zO3XEGHo>IREScum49Prxilr9D?;Dv&9p-=*C-N&fE&{x4#o5{%bj=tc5!~uLa_scq zy%n$L;*W(=r%`|oahlF=igRCBjD91V!qiI3>$;@uR&ihp>NY8}pI>AQ8iiCKyY{O8 zspkrBYGLq?Lo@h&eGNXCDN_;s;|n$iuKAveCMnnd%75v{im$Y9D_@_0_{79; z5~T?dDxmhR*d;`2hVrD?zWn@ssoNdVzT#^L`eu?rzx_WB`j`c@gWvy;ll~PTO~D`k z6{`LRJS|lM!8Fm$Ks~+)6Pt#z7J^HG9A48+bEGSOzU+vFvSpu!tp|{4Wp07E| z)0pMybqssAW20u^-oZ?yds*m)R)UgKLQs&-Bn-6!5D~Xwh!@Q&e*KzsF*|Opn1yCA zQg|Rvt0W5e|6PFfbr~l^i@hd+vp}=F43RoSXBhJb=TZ_|LU7TxO{=%WWLd<=h+lXc z{-|;}k6Y!2Or-ig9MQ;hNsvA`cPY5YYz<0QqP6&Y&Ba$g}Lum`RQ?q1GAls+-_%i-hN@fD-(|I6P&ys25xW* z(CB%c&G_bE4No%yi?Bl#HZ+{651HT8P9OOfz4%Dq5&CmzGKiXoe#w%A?1FzH9YfPa z1YGB`JL}O7+X4PHeSk%NT%!DUV%w8Lkc8{RWd6=+==(y>6B~pWyGF(-Umd^^_Q#O7 z4^)RL<9RWlt$4g#U6a*;Vg98N9ASdUIediWCHSsaXmolX8zP7EOLfr>uqQ0-sOLDo zh2K=5Z+NBmZyKEFrBAvWD0?G%Q*@Ir=q6{KJ&ub>j?UZ~ko%$zkOJzwE{mSIHozAF zor^U3$B#pv2;r3Rr`u3Z0qO{9@bmYNF6P;Bh8|ETgdE|bE)-D1jUG)X zw>^2SlW$W=x$AClT$+2kFdzbka80i)TX65M)UHzXw6En|2$h?7}Iy5z|Ed z4KT$jmu3TgZ!4al&uek6-yU>Or;~-?fRQ{{^qUASX2EW#a-gvke4={bfP@uPN(5~c ze2%#PlP2ILN1m`IPda7cr$YS2Q@_mvSQu)??YCbcJM@yc)E9N_px1~R4pMChR|gO< zRF~TaD?BZ*F$h~@l2!i=hxN+&LqG?!tr+awt@y*|1jFHGeG4AY)@j{rd}dDrGtf&2 zGW+VEE3UIR z;INFBY&cvJv>P|^-zFo-FgFB3Ozi9yE^=fmt+>65bGGftLT4=sni}K+74O;NW^TR# ze{x`rr8++F8~yx_b^8FGA;#8Vx_=#8ZS3>@m9DG)?g?b$ZdNIPGd3=yv z`*QJuB=aYVg8eH!oMd!9?*Z~4s1F^@NW zrFuhNGx-kJ$K2#Y-W)`$fHN_H+TrC6kD z^d3;F)rQ2zf{DXN8{`)*ECbm zaAtUT7`^l|$V~{!jtL6h*E$GlZr;2Jq|hYY6Vv;hkEhh?+uPbiu%{*&7Oa88`4w!s zcE1OVp*s2cH{L=fl)aZdga%#%5!)K-3Z9tJ6DS%mJ&%78Y|DUuiq z{eykz^Ime>qMu0BOqMsm9nR)>*1rCxrfdpolz+t#}G?A5Wo*+xOOq zmjA{G3ceWT#9zH%{(@HBToi;jx;AVEiXdk&aoe&SLEjGh_d=++h};=CRxmo5##c<+ z*34f5{U#o@zqJ6o%$l?GRMzWRTnm}n%OUeDg8kSBJv)J9uz!r?-;$o6{b6Sbi-_=D zy;ZXXF$Q}eN)lE4PlbfdaisfLFq`N?NLFw}(E?aTf)M1q96Ld7cEo38*V)3NlHt50 z>ukkx4gW(Lug9UR-$~dv}5B#>gpCU z(-q82yIVJ8EeF_EguG2Z5SKfVD2o^B`qajUOoJ>h&DG0xzCNW<9&_+P0+@#E&KuUx zsXIB3k7CrmdOE$+s_&a@tiw#MYtEaA3A@T;`OQw0$g?u{gE2+>@)vOJJ@WO*M&`_* zs|w ze>z)DY-3}95sLUe7f_I46hYz{0zpM9Z`CXTfZ9`K(ljVhBv40H+pWbFuo zB};##=wNSuPBU>I9AGN$ii)(L83a%N#lm9F;{6ZyM}-h%1mD6)vq2PrH{Jv)jhPUp zn|nkx_K~ZxU4(6>t5w1P4+)ausZ*OkW)lw_W*JxyX?OQ=E&?YU^8Vqaz`!DGu1v5P zLX1VE>yWvXWG`E|`}lD##SjjE3?HUA|7W{h3xyz*V2EKjO>uP0CA4En`Z+9w*F^C; zhvC~4tmfrig|tv2;`ghhkKS>XeSHgGuAW2Bz|7P(g!zABk5AoPg3}7k9P|t^=n*y< zPss#gFo4>x08)vVmN1yA!7@OEkPzydR6joUy*iN61kxVOSR_^qH7Mj3ZNcP!- z4&uVNWX34^)^M6yz*TH(%NW#K;Q`S@0A(p5jF&?is zCCaPKM;_LFegF49Q8HV?A@v0E{f%vIShj!*0!7m7ahSQ$p{e=bHn&#Fm`@s&U!H%? z0rmSVcj2ZS(j!C#L7C_7k@zM{)`Ss!f`SdOGO}r2nEwghfBz5jyJ^TudjhkF6<>Vz z7N7o(+q)$&AUJ4d8^bM@V9XBlzkXiurDjhX*?aTGx*z{LFz+;N z7lHI9Q+7w__7dZx+N~X{5h$;*iP;rr=J2$lwiOA=n<#s99GN7%^d=H9GZDkr1tzqh*Q6`FmCgd{H7#fXCfsf%3p{h8kjL2 zxD(eY4!kACZeOyL5CDblKS$+ePL2#XKsA;SMN_Tj8>?pde-^T%RDvnl?004C;SW%I z#&FvZ_B{=GCSF9)J6q4M3+i1!Fdb{N$8boRAl@y6)&XRqeB$FbLEsKdQ&yw|ws@vn zJDgquv-r6=Wj*jpu+W8UYQBf#()$bhKORhLg_6i|*Y27#cHK_;S;tLZ#9kL!%m1Xke_>E>h0_GZLBt{yK@mt*Dz-VZn-pn?u4q#Hjp`*y?BEqi! zs5>}1XZ~fS@I+h4&NewUWh0}37DrRzTzxiwb^*e-))O_`^hD9x`>j}b<&M8Ghg zi@2(QSg)Yq#@h~c`RQP8QIlB{29BAeK^&XIDLIb6=2r8|x#3~O`ksGmSI+xVuwC`O zG5pu=U?XBV*9$x{^83ovvwgrA|1(@9Ish-nd=MaKO zr(3}wFoz8Qf25a8+L%*zc3B*ZlYs7keyP=qi-|F=IX!ONfwtUyko4o&m*8zSQ)`s% zUr*lK>hvAE@|PGe3Hn{}D2ZewsTa_sUQ5#&o5#pZR6VPhB@SC{m=N^)&wJa|Ur5>R zt*9bclJw5&JnchbCRTX}PxQIvxM_oBgjfojTvRKKC3x7iJ_ncicNV`3v}oR>hMRF14RMv zfmC-xFXCa%BVi||h94#;jXW@1BD%(Gk6Xd7>zVlvy+(y^4LZC`Mn@}e zkMF%)X`RMhLp3fVUTzb^uXC5DS0D3CX4>U7)-U;9N$p)I)u|k5er>k9I7DXxY8ZD? zfqTtd`_OkYq$0x7y8Byl%Jf4Ix$~^H&&xDfsma6ryjGqk`#NqJ*VL*vOuinfZ zd9u$ZBI5JJPJ}t*GKy#QcVfYF#KsxYFw8YI7bMtqKJWWEVt^tXkn~R5>sK>d`aF1V zJl6!>4-!u`E1mm*0IOO3#E&)(?V&lcguP_$d%Xbm7R>Ir{D zG|yl^4+(GJ&D~=*o>@WvrxUG8^oZsnpaW!zw}yyT zv8zkUJs2qwK|F}fK-ladLql7U&PvBhcz9vZdyTJm!1igyvj8HmVmB99F!yksBpxzo zOu2X?yG}jraDw0AgsROsa9?q7a)LU1aO7u4$Knvh&%Hk+I1XeCB2E7p{v$jgy1~z< zZsIWSQgnR;H{oM7(?`fk7nb{!*15(5{z)U%`We@gpJZrt+@jTKsWYtk!{ML%*OrML z$HOyaIiqL%Cudi@mWE*edsT}EF-aqVmg8Qv3krWeOsXz=xFO7`GN*J16B0K!Lc+reG{J`lZ)K!V zh46l7Wo4bsNQOq`m$`<8o4dOdr|rp;bQm9fMe7||Rg# z(ncl#10+!g_<7*%##q021NH)`qR2w)dLLW*DAsE!MI?Hdmx=2c$t$~)HBA9EE32wT zR6a-UJn7nDJ6*Bi;O$`Wu*T3i6GsbG2)Q-7h#h^TTSKop@+NsM<-v?`?t6YZQOTkG z@Dh|6=+OsA>8u@;Fd|q7a~l>xpQanw*3w!k>WuU=@gY1?LusICUT0ykX}~ zF0!n`*Q*ka=U0C-4*{B%ED=qy@1tf_U}Yk6aO|gQTJcoM2muaJQnm&e6TlE3YhhA6 z2OyDcXBULyB8bX3AE{srD_^!tI@tttNlq*S|zMe-E+O_|N(Q-7Dk#1PDc zZ|8r&emN2Utkbtiqh(gBLb7a-+E# z7#k~v*Lb82wFb};1t=YY58y*Lnzxe=J;MGbbo<1cw!g24yP8^!A;cB{@gdhU3e%bb5x* zyx3sn>5$`*n3FUz<&b$XY2?_kGI!56Gi75rE2vSi8Q1#EcMA79tKVv0k4%#Zm}!KZ zLi*W@>3-X=!wFAMi4%PRiZT$GmTE1-Vgyyfg*q7e$;)4UOomh##VRKHXm4x(Te+(x z-*~RPHRID;w`08D_iI5vmqyy!BOt~K1TW1M@hb^dM$cco=w>j|xq;n~GJ?1sj{0@L zJ{~L(u(h?#`=W^{^XwdUE^y8YTWsho1Z~`B=Sj!wSIMDA(DBO3nn2i-R|072IL$$B z?&9Ovpy!?q?Pfn8`KLF}vDeq$#FZn$M)A*aUENEt)j@9?2PXm9abQFmgzaY`g+}ip zGI`+B;i2rHn1yP61C{x7mYK8^*kkMAfM4+)K8Zq@6OU7s z0tZJ%pd~R1>RRjrV*?g2yR{=86dc=5lIQ@1lqPL&UF?)>($AHtAXM~@wS41Dj{I*U zY$#P!52RG2eEt0SFQl6ln^L5hJo<`%)Ia)os->mH2M(mCY>tm&y{Fj;SAg*41~ktS>Z*s=9C0G`F!K+{yQ+n@VXV|*ElwaU$& z$z{B8ghLz0BNRe1xxD08Ujr+F4zN>>5K*N{|IfHc=m>n8f3E%W`?qojD#a(!rlzL& z#KgUDT2ykEKuH3LR?p3Th$V>+Q(K!2+IO>G^-!fE5Ah>5p1YN5AL%zSd>uNx?M-du zWxv^V+jT-f?o-K?j;|wBxu3&*WUxW^Y z$Rq(_PL;r%6T>pVJzL<~_FQ*QIz8 zdxZME9V?aZjH;DQ^lY5ym(jVh*_QoVvrJHUUxVd`JA2sRv?m^3uz0BRMmP!-@j>C0 z4Nd@*!k}Dg$+j58d_g31P~#BUK*bG-Wgws7bf%?tKR9fRC=%DfV%9;N+xZuG6rwmI zMB$GR8CId=yUon_$y|V=)LK#8kd-9(WM(%W`pl~zgEm2PmE)BXeaFTgW?yjOT$oVK z|B}23u!gwcs{^Gy|7=5A;{EpsYY{u$lpBci*W(~~L9!r$t&nGUKMED$Rb}N0$mpKA z4R0V607~P_ni>|0n98sT(1K3f0~x7P_g}}4LeW$vVj6U1DCgG`Q6<#!ZX-3F`?s1O zusmR_P*>B`*GE!@zp~neTk^L{J-bU~`Biu6O7&kGOG}9PIdmX8uR}=h80E7=Dw@aD zIczL#ox803DIGC6`S{*%p^Q89BgYco4PNpnSh;*8`$~ZeJFdQ^yFM#)3qTS14O8sS zqnBwL&V1)lJz?=z=E0h>u>2gmqmgH+wRU|o>Iu2!Nr_riHU!7`OZ%MY__=_8g+wbK}DFrfllu@REy zls?aZ8hiK=xq~t*5ARj(xmI%Rd@EVZ{*P(*{;5AU;!k7o7vr~zM)FvgmO0Fi2@-+L zl6xmm3oy}|j)zHJ5KwHY)gZsKl<)t}q~ zFY|$gOKFa@Ebk*tKzDF(hf+P8WA4oK1aS|LL6_Q%z7a7EGS+pe3qTm39YiXT0yYll zsFzY?9LA7JLK(nD9NSTfx9AE?eH7Xfc__>0cyn5+F3&x?(k&a60)#$DC}SA4&Gq-{ zwaqH>#e!!WBc>l>?MyELFWpMUN%Z|^oSd9C!Qy!4%u{tY;zA!Dj^H$U{2IkE`x( zZhP*1+}1#g+1T1+?_L`GM{5y$=2Bkb%5TP5g7gdw+mgl2UMkTyNays7q}M!pdfk2f zUwU6G(=D!y)%Tt*bFra5eCK!cT(KXJPfWUc)1O{5uiUQ&oqxEtc5t_7XiB`%joq8+ z;hkv1@omZg8UIr*SefyECH8Ny^_1Av@lMB2zH>5kSy!0!j`3Zq{3~iVjx}ef1LN_V z>662_i?pQ&CJ)j(yj3^<{2e3U_hG9zi}%Ur(p|?>jTS&pcLBz!EbT{`a}sU#`}*xg zAIvXfhwtQzMv> zaviTm85{OSo6Fz2yB}vrNcnUiM_Er#FE9EK>1oB|pQ71`3CYgL@a-}exBOcR05?TZ zRA=RLyHu5Mg=`TzZe;i!2(9kx=C;K*9@e zJoNZ!>Nnxl#c1XvsidToQBNMG82~6gaR+`nIyz*7@FG{~(N5TY!GynC1B&(`YXQpR zB09QP*hq~k`ns6qjo{CRq=ydF7af9vA{qoY%FfE`X)O#F+ng$pCRr2K9a)T4(G=G=QhTfvMWEF7m;-3}p3_9hO)&2}lm&q}olv z4(}=?#loq$0z_LVi)S6(vKmB6MNWX^oW)B_`QxTLB3I6B=IGEbyoVMN3ZBpI;_9IP zaoV9T7BvH>Bm)WEbE18QXBZ$*tnGC_L{tg7GK2`7Cp;5IcdTdyO-&a32n~g}738u- zGxx(PjKmSb&mr5y7ZC`c^L(Mf-Y<&VOB8xIx$iuP;k({5(+I9qf<1&e!}t}pOvHq# zyQ2v>535GB({h`@`A40$d0?hdOJ6_vNYc?Q%Y_ct>muu74e_2peUMOr8HMDeK!!t9 zx#VO^R~(5%cukmu3fstdcvqjFx6-K zyRqNvBEqcae&0H=&2rsr?eg5TvrtUfo~68ZMC?wg=e)nErI1ZNhaiL`WsM{yvssy% z`obGTj{~|21nOx~K-U$z;(qhKo#U3-Yxt?XR@=W0UizLzyoqSPS#U{F$YcHJ%DV^? z1^S9Bk`}~h;$&Zf*Ns7kVSYb5a?&Uo2W*jzK1}v{WaM2$OCtKac>k90Vq#Cn{{%q zZx|@>-yhRhA)^n;CCDqGJYEq|Q6>PUrJT083wF;emmxC>HV%?@4&SYE0|b`$qcC=n zoWN3~hL9ZvPNI3fU_jqc7T+ao&Asukap_g_fU{Xtpe^sVV^}6d*^OVvFWiDG=w8=x z5{Q{%V52ER7OHM?2&5B&^Fa0eQ_LD{&-~*>vZe{zaA>3HTuW<5CEm(Nd~z*Y-e11y zkC~afYSTvZzLd8O*F#n~)NZtL(%4q67kj`rcSO3XY@!1Y`yQ1OC%`~_fz3ni&fP=# zJQ~tH;o7dWEQgdF{~Y8@YFq@hl%lUm5=$ihrf>*M}-IS2^;1 zwM=HeZ}FThpL|l>zWvg>ZsT)aaG+yh-Wxd=gfFJHjXmGDEGGz0wyK)#ZRk@`GtoIr z<`qW>|EZFP?}={FNfL|8M>^d9OPP)iY(k*k?m&KhIDz;t{#S?STJ!f>dU`=%UbB$pqOh$kRdIvImLZq=?A}q8K{l^Y{T}RcD*3bG zTAn(VeF2BUE+6l=T5(_|jV_xRr@t%wbZQ$E_*GlSe*N-Oyk{}UMgc|mw0R6x`%#N+ z9NV{VzXSD#UZ>D6*Uf!f{ePn*ft}G01 z)L&_|e&fcyaT;u2(JzuH^y$epyb2i|**T5{n@hj087`vZxr6L=fR&xcIsWT(-Q!Ab zZ=4ESv+~N>`>S4YbuvIYd|%@TpWpxg37dvOb8r#k;-+gaR)#I6qcg>Z{Tue%8_9+V zN4vQpR;{Y0|G3DwMdIUyQKu7gUrKmxaOiqMmx>9SC>FK}VUdBtBhA7_|1$__aC8(S{BQ2=?g+}~gLH)=spL?2nL-&u`GE9o zfw3&I(jC!vWFY`U-_dh`kaFTTZcuxdutj8;&&_4$)NH&Ti)oFH7?8z={;xP>^1;4^-i;1ByX77%nF2Lyn8>tt`Y}b5kT65NgsRm)}Ab)%^_DwU{ooe(Bt+P^Rs* zNJJ@J$Ju)y8!x(y$)XdClT`y4*cZl;E4d@3Qzw?jr7q{XtgU%&+`qB6VfJ?BN zFb0;V9;wsIk;cPZED><_n-S~m7~5U{(8;1;Y>m_H6*B(5=On9*n?w52p9a=%Bw&}B zjYEg>!%n4w%s}bz^^GHX^jO;Q5mwe>`O&<9@{ymTzd{dd-&8aAkP>ms7Vd7`8F78JrXN{%E!pNqD(6&OBkyQ8Bc$kj#oqhGKd zp)OSgiLO@?%QpCXz}fKH)(QvnMN36b;K$L^#5eId2$uh+ZkqWX&KiO40Q>^+fY{sH zBlV#OFi?dhC)u)9f;*KuyOgoMl$DhWiF_9EQABx5qPL(~%@~AA42{)7(mCNjk;u(i z!B|)rZ*f#OU$kCI<^YqtxZT;E$Oes(Si>8f9%DY>p0cOP)caMlUiAj#yg+=6EWP5S zb|NVTU7L$T-Y6gZ3$j<+QBkj%5F`>f-daJycS(*B;{*c@%PgOwM~zZ$*x=}JAGv|V zV1tt^e2+@P@=db&W~Wk|7{8R{SRyWE|CS)|VgNdEzRH5A$k}U~pN!0Z{k)rNTW4Wu z`I4k-AsPN!LbavEf-}kF+#)$S!r*KgR#?)3pIqg5Z=$aA>F0Z} zhlGTfcZdsi$o-SG4enGyIEY}4L84IWt>)yT*MEKN`pLnNW%CXVC=<7>^ies-~ zI_aK5gjv;z*VPe*`u57ro7NvC@i2+)5aCVdA3b_hDT@JdScecGHogSH%4D^<#TOmH zvj9t?vs3(CakI|pH_|flYoJaL2h`pMoFCYRMa0Y%utr1-is`ac%leI67c`vhtnPzOq4oRsO?sY!DypilNQfcM zZy*Qf*m?tI1(%x!Lqvi#5H>If*E6E^S#T{u|GX9k>BD;8lvb=*LGpcot&qcmlskZ- zM(uQS$S5UaHC~~b%?&{G5O!p2Yb#vL_9jes_5z8S*6bJv0VrLSV^wkd2;ys`q@;+# z0K4CI@9CW++K(z&Zmh>M7bxEe2{`;;Ec=xbz6_`*JJijJu+$3+;r z-Ayx#37V8CkD(Ji79)-pA|Wpb!@=2~po;0JJg&|E(_GqoZG^Q`ZvH?}w&&szuxSQLXn!g0=()PEcjcd`VTQsWL`MW#Br})$MLlH9mUfRbLO-(q|~mnthG*Rq-Ot>{A2H2pN}h3Q0V)wXGiJ7 zc+QEF!4$y9V*JMca3}5HLhvo6=f%6j0FoH97B48Pjp{c)1STXW7XzTj*FT4IBU;I} zhNOV{-f>GSD-!hpLiy^|tJ=Nw^;@k*#{XZ~CyN9Bd$O8N{o00T?R)QGY8XINa}l_G z#N6Z1ilx|DZ9eu`rX-DsKe-kjz7Pc%=lSI1v)moMBpp5(#h1COQ+de5i(jr&+05yD zeFfCI3PwhEK;y826<8-CqJjLx@)l7LdyGOtLU3&DokY&_I!Q@Q41)mx>NanXl{NSZ z?d>{XF-W~t0DcpQP36*s3$%H|c>^z(nLxa=&}Wm7P}XmWEUB9AVF)rcfw-%Zu=6&= zaj~EVt-=+zkz?pS_grg0H|7r?!*68vzsCA5#>lCFJm%LJyFXb^;&@qG2>k^h0O41^ zy>l}Ou_mwAG|qih6~$!(Cf9^A^ElZD$p5YT@+AZ_XVg7!N zjM#m^GGU-N5M!7jp_5pU>Givnjzks?7SgN?09cdx#c7oerz8KksFI~#McB9ENoB8! zu=rgqSSZp+$D-~GzkQ^cxR@FQuP5QIcH~F|RIV(zX;i^afxU-Fhv!D&+$iZ)nRFh0sUe7hM^}$kqkQ6j ziQ>Cc=_YKDOEf7hD6YPkSRM$c=%+i!Hj^X@^pY;_I(>@z`rOcbpZT+Zg2F(rmo8H` znNHDje7;m4>Q)$f1ZN@{Ed;804!eIcsK$z#6j{I0&a zwj#{y6<4`K6kqvU{D0X!OUIdvresfER?vl-i4SqNSYl^8@DIY+;Otk}LUwfEAD?6) zdIX+`FuZ%kpizm15n5G%h|MFJXP9exhh?T9YDdel5*VNDsMAv|9MGg_YI?6I|3w65 zB+2}^xCQtlNw#~tPFi*M2_#i_xvY{W2O^N8kTJy7)xWhEOZw70ud=!l9&Qq-%4=DR zYYmk?++cE+34L#T+v z6hh8uM~gm9rSHQH`v1mF)^#{-?;JhmIOq}?<4{>MDb-D`NtaAYXmS8VnQP|AY^bge{G07hDp_8VpQpifEZyJIr63Txpgox@clELiLIAl;=9Z*;%V+dfaHjYvN? zZ}02vB_1LXMfb{@ssU*Xa)r~Nx zWByNY%sd@X3%p0FCqPv|Vsq{ClBG+}qj;&SLDx%aOfk^{V`Sb3*aP(jo7`dw=xHfJ zJZ9dvd48(dAsk}Uuz}=&F_Ll?a-i$a$egb1in<&H^2Rx6iu2#`?SH-6)O9w!z%jh5#Cy!POt->0v&rU*D9<=Qr|a66LN9QVCXaAnb4C z0Pn+zrca$V)Yo=lq|?EnABSQ24M4plSwTyA{aisGkV6BRQ#?I?!W3UjwgR}MMhs#n zpxapnE+P%*huA|F@&C~F9&kDLZ~OQ~WMIdeIYD`VeiZ@X5|~{+?EfB@r!;NLf>}`lhjHX-xK*ZN^F0{%z#(|n z8ndGKsI#IlC;_8Zc7TE<4M-DnW7-LVdOomZ+Ll^|gwjCCU?)Inx%ANiVBQ z#(I69oYgAWDuDy3j`5ETdwbQaiUfG!xnEHXT^Cn|5Fn8JIvY8s@uB8Ix;6hrB1F2i zsJ~MFoCsz9=R~OZIk1!ZK6QVNoywXcI{)|Iw;Cm)0?LD$RhV91v#Tu>Q+mM8HrfA1UW*wDOVmNE90Hp(j+=(sIVYp1wW~Cp?WS zDG*&YIlK&keha=8FzP)Md4z>!5mcIvqM@(hAQhG-V;%b(!?WT1b06}8@0fR!Po_W- z(E|QG0*DfYiE&Mh$cey!z?3aVPgFxAM0W_8Ufc{&tRG;9VGFl{%#D3cUK_A1-Z^h( zdJ?US2}vvWb1IVFRInz*Kj+nWTY&61w9&*-pNE&Xarzx#EDwEn`xy1^R;WY=Zmd!+ zU&5!4ACGBkN1YEsCmJbl;HgAtM5s5+Ml50VwSRwUi&f6lFG(VUfb6lNuoYTaKd;pu zXCaxZTnK|Lpp0J;h=`*QQiP{IQ%V~&5rcF&K+;5}GhVrJg@}X{ zKu3{_M#HH(M$CLDDta60E_7dB)p8*@_}nmU`yGf?rdLDdl4t^FAq*S_ATW}dJEU6o zihn_@F<|sg_zp0sWIApg>J|Y4?DF!ihb)bsaGw}y5)vJK>i#%@Rl@(KXJ+m|%RF&9 z0W*Udi};=5Ssd<_dIwhj)WP)_b%Zf%dt0wSmnCZhO;(JzDV`n;g_u0zS|VT1Mw~cN zyyiq@C%ba$Y04_2R|9vNa?nq-7|F{3a8>cEm=&}RY4>GR3$nFXQ@3TGuLCF`04r`` zRcz@8s7?GE6j{6-u73+DFW|irXhKMzxv;gT@Sr4$GWXWRFOoQ3d-8yM3!jtI!c>o_+0S<%$L#rt;(JZj{eppCuppOEE~9FbRkXFuPRDSKcUW^}mz6-5L5VF_;D znG_njlMKZB2IR}*6QLh~X%mH_K*48-(U7xfS0xWGV!T1=>+36jGz`KNPW~EDcT*hD zr*Z0BwQ|zT<7)J2o{^Kg15}l)%1xNrXWrLv4@$F!jxEc9BZ>mX%pb&wqJp2H}ybboV(*onGe%32VlMYo(3tNHj zqK|=e7aTI{dQ_G=|0tOnixQnL%Lvl&w}?u z9_o+D(SD!m+)P21E&RsJaA`TCth_8msWE!XeI?wrv5MdsC!Ct!ac-dYxDCJ;ph#tx z$(pk@@Kt%n@%BkdXkSZT9~o_%{R^cPg@SFWWY&b7fJksink4Cs&q`*gK$G<}B;;$n zXZYb^;E*H%pit@Os6g~d6RwnaZy?SSw4EqwP$*Y@=fZ;^NFa%ffyT!jRpU|&>iHn9 z-VSJ&*v5kP1paweSy}HFe;*$cP3dK6O{gYt=K9*IBEo^+1h}y>-D)^mRodIjiwI?K zyeaZfKQmx~>z(Dh$dIz0S0^gR$pxK00J9g=kuvm&n6-JN^;ry(6CAChND28%D=9`*}p-A({vR@|5Lw2g>BE{iay zmM73XHy`T*Y^Ke&cmCZdU9@#|%UMx>zD3@_i8VqJC!GuR_oDO?lgw%EMrn{~K!lkl zr7vPPE!NdfWD-B%FZUV}b0h>|U*LA%SthIj;$%4#bm%tzyccW-WC&^+!b|KK2qDh{P9H4kX$cqE&8E}SPxdh%4c<&o9 z@JNtWf^^%8ISrue-;-$j9eJN?{!aHJw}?j$0Dqv~(BXMEExj?Rbk;(zJWPVD2xOtFtMfgYCgx7#~WO| z%UALVXD>88rXve!X$ke{=jTU8QDT;qgJub)E|lt9A~z74?BTWn5KUyeB!5R$ItBgv z7|JrzoK3O;Jiz@=H{OkkA$AykR}ek{ZIW4mcTjNuz~vJ4BGABR2jg-fwE={&JNPiN zDErC!X48mhJpk9@AIni~HunsWdx0Sf!_=lUc*DR7l#EX%jzq-26R;k!HU_xKhisJH z{|gevv)zXT1a>epFCc!=aCC!E4{EtBxhVg_B zh(suE00c1Rk*pV|B&4ALBgi+Id4{}7 zDGLkb%Dg!}45>t)nY*%XL>=rAi?{#zMZ^6yTwzeFK~=%CDHzcQG*}vvDj9+6sJGF! z>VQGbM~^xT(Vine?RFXFF~sQ`$2X28GS0+t00#)!4)_+a{zW|Q>JdhUOU}F;+euW5 z+jx1ab_^@g{AO1EXbVenI}9;6B#z;|BQ@zURl1n8u>+$L2wsG`5DDrvSk6cuq-vO@ zNDmDhR)`NVDnK&XjCd8{#4LUfeb{A8725>s0sPTGon8X0JXA4c{LVTPlx&1S#`4l0 zkQ!*K)?=SjPsfVl($JuSD6PY_8|547PVnw;1xFhy^iM-DNapB4dvNM1Hve&)_T&^E zFefO2q7uk?Mw_3%uuCj+uFJOgla|sad!(i0G)1i4WT2+km9k-lx*sfWml66k zC4fv$LhO4Aw9G)M0N3-iCh{=PdU&27yttviD-Z;&H;CEs_DrRZ|j35XqF5uRUE;!#wr<^RAazQNb+q1bhhpdf504$T_#d(=54 z1Mi1W6uZiFO}*PQle@A$j_y=wmV$I{3q$MlRipEtl-16$%ywU@*k)v#n%1J?R@!&x z)o`ua!QN#}vd}kAGT|FK=aSiiBl>o9In)hC+PmutYTGu7|z1^X;%<7#j0D6@!n&JdA7Y79dL?@!RB zIHF)stgEO}+$na;g*r_k_AZTrch$W`iODUWJRm@Eu;lkdN=DLGsD6YYxRJF4Qw+qR zp@)koLvUD`QD5pS=@AYXv2FPMaULOq?USRSFLUsoJ=W;3JLxm3@&Ar2CK(!;8i>3n z6pseC{W!M&XVzyM(J77m=6j$;69Kk}yt$__YvZDoj68$7C=%|vaCG# z%-M?H6A~;i&cx{uVC)KjRb=i6!t_3Ba~E=RAE5-`hh&xr(hH>9r|}uiEPl=ZDIz*P z{?4gz5$4^6ql>@gAHVWkOmc)bG1H=nH8>ca!M+av>`Y+s7HIh8*BRY?Oi`jrbdPNgpe26j5p z(WO0%vCAqRTzzhk=rzO5e%aqbCQ;U?dvxt`x-!4RR1IiUC>|Na3N~OS<{C?~+}+R#&A)H-P@s+MV&tw>+7+M<*9) zdSCK9%zc3WVau@1kAzFj|t*0u1kusyIKfqDorILqcT zz@!kgY4Qh3Lk8FlRcAIM=FOa0;d~r$kB=>CljWw3=URG~ia2gvLUEYc{&4E<*1gFt zE7dXjX&)!23_d2bpF9V;lF)sKv%0i21Td~PApq zfWy#z5?Z9mx0uVW1{bBSo}Ogaj3HmtBf&Wk&Ubq&8!m~;?d|P_>R@mm?j>KK_w230 z<2b^O&aM%+JK?rYNpw*nkZ^XN!f*(#5+?A4_u#?NbQvhk+y_)tR0M1%k`0R96*T19 z=c-zwA#uY7+JR$yHf4v{d`gJ@fDe<$&f87v z72h^q{}?C2Uij_#q`pd-(rVC;Av86=AuE^nG`&z>`SdNGlpdk%cHok!WU zrvaB)x{dleinO~RWqbkr=+BaZk;9K z)vJpTS`FX-_~FAt=HkG^@(DQ;(KTE%-z$&zk{kre+x~ zX|+1HpI3C*mh*KVz9?iI_Lcv1lmCeSy{glj`1E~I zadw_tiLv%vG=VNp;GVF@#aj7#b2>kqw! zGkhneyeegx$HB)a9#|Pz3@IF}BX5DlY6I`DgSZc9aulG_sO04Hs0$)wR#F`16@>54 z(6{sC6K2s@r>?!(pXB2o+r!WHfWPb}*8VvFTW#7NJGN@Z2OF2n-10EbG<$n($)za# zJF7XfMkJ}1Zah!l`ZKp?uI`gK&5GL}Ql&Kv%*+DiTaE0)Ks28e6C2f0Zfv%m9oKz! zhpu1n*x}Ke$|jL~^qp(trR8^s8Hd?6i68b%ciLZ>?zBptyCRK0DdTLxQT5G7a^vF% zCyEP&`vz^8OQ*!Knt6Y~*#cXz{L*dKa!ps~6?*C{%@hj7=NGx~?P5!JK^9AKU{Bj> zTtCI!%a^l*jY7$L+{!7Ox$3yvx>L*5=TRuiY~NQs#R>`Tb)D}&*w}I4%SOC(%+1`X4w2=*|L_3@2vR6b zzBjKj(%^S_-pZG#FMe&c{dn2!O~zqoE>qv3EAh?xOT$7Za~vC%BbpqT`4zg1>y7}o z_PO7a zHJl2yt=&4LLwiq(ZdUay{i$VNAm*W2KbhBKo!xC;oKw@dIs*3WBU4*BOgKy9r3cj- zIc|Tv%6+ow-O!C~$qn+EG3@`)f?5e-fWNV}*URC`Ans ztuM2i15qKZpciJ?;M9KYLp60Q1wX#(?$KYN>L$!zG!Efm%{@9cS>3ca7SzUQ99x>n zgbF;NFa5`g+s@9Thx7`+U4vSRqvHqOZ9$=$T_uBK>9^Vic-G~B!9!3g4bhum%p}TBX$--vGDvHT$jr=qn1bhX?ZwewJMAC5+Ii?m03NxV z&-!@jfby#<>D?ND;!$KC35QCMxi=RC-Nfsab|Gq0-#6aB+hO@@3|$xJWo4hC;awU` z+qv`SAA#I>3cR1287454bYHBk?^~BVeZ9+xET;OUUiF$wZAMcuUvZy+K;5~1QvuE_ z1DrUA`_hktfC~{C^eY*nKWIB)Diti$ICuk=yAy9V=7bb5TAxPBn$=^6m%Ugt_o*+< zE%uMweTsL4HP!r}$CzSsMz;{CgK!IAjo425%!_iThz z0u4+aWU+_H{9Xpj7vY8~H|FlIdz`ZdH72Oe3sCX~c?JZ;l-P|VMJFW4D=8_tPbzwc z1E|5worsbNS}gAbJU|)GoxNXP9EUm+$G$9%9UvOITht&_Mb!&;10AT`Aa(VTDT4hp z6g+H^70;3FFdVs2s%I@P?*l{yoifaH3<3eP8r#B1HEVMSCVrGCg zRsweCz5e|Kjz%5aN_?Ce5DSi7TM&P?QA5P7LJx4dk9~Y10ZnV;)>};t6o6k2CyEI< zq%<;UfOUkIVh9j9a%h*6@k()^9AP6GnS#VwLKepb)Um^%Ow;1zkTe5V#jcWxZkY=~ zO!Lv+`wWCiCz7iSq0x+z=MfR2vp-kPhM?asPBCeBlu3u=Ttu5OT+7hI{!C9-*L1M> zpfzKAd%L7(4)$RTw%^=AvAoL{y>DfJBSPESG?9a7!vYt=%u1<)jq56}W1Y(=&>4&$ywzo|_G66*Vo9Oral?Si z3o~-^!9ib9^p*k2?7@s-z#Rv3$E6Y`jbQ78rKM^{#4ZGjN=ZwL2@gL>tN`HE=O}mq zf!P7;lH;>SU0vO)Q`Xj5B?AHi0_T;LV}Rr)OyCQiQEY+Gf6n^Ge98VC4xSwkp7_Bu zyGKF*ekbT6s^@bW8Ic5 z-4~d)Y>8=WYh!e82IQ<^gM~`85@`QGG#TJRO4A>cl5zq0jVxvfU*zzyoBF}j9nd@# z(L3(F}w{>M4Z*}dBCT-H8LBObK%^c6bsToU#NCul-@RYooPw;1$SoYito)v7zp zM0c>Rsa`9dw@x)?@SJG9iGqh{{j%-@lZqRf40C6W?J+a+*Nah(TpypfD*om|C#S1^ zdJlXKIoy0R)b2N){B9FDWPF#&j3FJ*?yKxE6|inB`1bR&#y|4STn%q;Z+JBXm4l~= zE3xe}AzNBiH3%i_Fh&{%d)MQHs z)p@F19qb<~hFKTC(gEI9Iyyk|6Ordd5MyJ3hsL)Km zt+)l3sd7jUc1ZxA+{eKohrbKE$ET1E_F|ByiN!Gqr6f&PkT~JsjJ1*O37B>~L{BZO zV)n@t0kIm!;!?m2?}>(Xp~%ck0a6lWg5_|;HkiAOFJqAfY}lKk@Rs3sj?m|;4z{}+ zm~sU+S5;+Y80hp4BzyEwl`+GTN&;6HG=@mXAY`LJ>@X0b)5?!b5*4Qe?n z`}l5fANrj{MIEMnnh*R9(Ym#ruS~_xa628cjMBf37iy-)QHaBe)5W|V3yXHgT(#jm z_jLuy7)u{GZK;YtryT-2R04og!7 zq>`w3$uf$=%@0>^KS%jCW(Ee_gF@+G!pYDg{#0wKmW`eL0yax^f?5KecJUMbb4XB- zu=j95;Cc+=(Nf_0Xl02+FR>09vIMSubznh~p@gmnXm-jl@92@2S13BfyKr?uzW3V~ z-jn`I3vhoapnO{^G?1NFRfR1d$5l;DOdvP0wfjt=3~5mde-0-8H5AIU*tYT2{Ic5o zY=sA`Gz!p`Vf$^FS}EE~1kI3ER0ENY$0rOPgQjZIsxPB`2soN`7i0-$V100XlAbuA zs;a9a)X)IBC z_NMT@-(=A5e@@fu_9dydoKq?c?H5LetWr`_kql+Z)OXD^DW5k_+u=;t$s0d-(Rszy znIiFE+5;%_4jKCEhA|+Qxp;1|pmVLgg-g=bi_a<(=!OR991}`NU6J7LJaL3+36HQR zQvKqN`tZ=u<2dUTIc~ z;DVE<9nuC9k}QCiG??r*6=*dTm<~ccMB|scz@??K0x#8_fuUj(K}#iD?ZgTfXRYh99- zdrjpC;1g|#rRty~1-Tyv5nv|>F&y{5-jcQptBKA?GnDCcOhOHVdRI{W-dVw*`*Sn$ zUi(gCX9+&NLy0-Z=kB3BBGnPD#hZW7C}7K4p-fN4Og^AVI;iUfY{n81a9Jdf=*U9c z&LXsTuLLRGk=3T4UjbB(9ZPU04D~L*{>Uzg6s^?yj{K71`2Z_Gcx5ha1E*SZ_^TU36AorX--d z;aZp{#bK-jP?ZjJ9FRB%y}KE*{g66@_{Ar%w7C;Sy6;~US&VFJWufthmlAnBsl&{y ztdfYI$3}k;nx4U;;VD_OieyBg|3RZ5OuJ%MQTX`^5gp%Aec(?+p*YbhB>jccwxZv;TF zZnPlnuC%x~5&J}}sjbhm)GOyooEE~z;>5&TUmAX6j zsBGkP7Dsu{ppH_ZSiWoOs8Vw6N}u@AKxBF+{U&40Z~M*eD)GEx5t`X2#Uc}Q9;QB# zO`Wf%avcj>%gOyAs_5YHA}xIlh86wo?_aV73s7>lHr!itZ2z>?3)+eH^vP2eA?~*~ z+hqIN8PJw82gg5sM?H^Qy^pULo>}QM`NII!wu4qDz6T@xbBGne!^j1etFnDJ{Ee#=_C9RLW>Mz|Fc9Ro9*J$T>Sizb zfu^Xcni@T@@25_mK98ycENM9G1opA9odyj%P-hR(7>@T}(Q5xLH;WD-HR-Cw_TwHb zwA=IVUq(78H;Rm7*Cy;_e?9|$4D8q~@u{Sw1`>*)e5RQ*p3wh>y?_5c#Ldsc<0%mq zqvN6yHJT+p&2N0|7z3XBwpBM<#ESBkvTIUvnWs5U*RYpg>dOhT@7?>N(H>qaq=Jx{ zhT1s!%16(T3EaA$64#bIwVbQ!3=6kAOU-ZF&5zaW=<<6t&pje@$rS;-TN`MmH+;oX z^#N{YROp6eHUd+G^ZukVs(=C?Q?$0C$`^fTQ}9gv42353uSd=A}mfdTL~lPe!fK|M$k>L38dgc-kl7E zHG{BS!`eXzh4oNQ-*<7TOt7EUn;x~FO|$#u72MH3)7O9Nh-1z>Dn23?wedPWX-EIW z<<+BNjFu;0U|#{ppOS$dzb-`*sW);`PiCp0BZ_#m808cLeQFQ? z`#FMPo5VhV=_gO##(GIMY8H|x_hwgq^zdN_(J-6zbjN3LMCAnJQDaF-@58I-TC}Fg z-kw(=W4BJRq&zO5yC+V6uhb? zcgBolPFMsL&;+f%)+{xL&%9)t1th4xe+jd2*PrF=-*63IoJq0yBP%o06vmbHz7O2T zu?!WVzNkcS9e%TnbI-wp@=&5yqjGVyQ=M4jvr4tK;m)i`pLh~rL$l3#tbAQG8U_;O zW82w9>jPE|><#UuGmQv_cfDnewsokpzpe08*;i{7&xYeqH0!y?6m56)Y3A>8O5eBH zhV#&-@RO_93i4)7^Vr%&G~3fm_dAy^X3M%Qk>=n8XQ7xkTC%_>RHbXPATEt*92PXiPCd@ ziwM~p#~SYJ2cQ-+#N)vGH?^cWygql@twm7c^po)vk<-lnMi@bDK8)2K0Jl;+lLni37CCK;=5pI7!P*k zad2<~aG@w<{1}Ke(4cb06NyPEJfC#1j+K(NJ&4^pR1;Cr(Vx~P!S(p+l`FVsjO92n zhqIcer3vo_3DrNUwXp&=R`TJte|9IB5)Mio=u44+lt0>g&J<$slG6y{yKzBQ-@QCN z$#>#ZvgFW4Cyyx|&{RgiP1E`^ESZy`#5I5xh1ev3FFvoN6pg|`rU)b8qoARV1`!8+ zwH$`GVz-#WcyN%$3lW=>XPHL*E7~>YpVRkR3}gD#|nFN(RhU(`tzWLp}38N{IeS!Q$fA}Ij(AH zv9hp;p+)vy1oVxB+$7=ym%g*6$hcL`5y=QNM)$5!vF=5AM%l>Bpk-saSX32aa%kg0AA zAC8Z>V2tiv?gjm5!PS6Y>3>Bt}iP zkyT2qOn%wVtY@Ku(#z23n)AQv;W^Zaej$y4-4WFCT2XSY!!u!qBA7=r58)saE zU;S{jyQ}L3!e~$OfZX2HTvZhg<4rwiB(vX}G7Wkqvq}(_+{kG3h#{zQd$^-7m-Vre z5eSs#5B&U)!5(j81zH^ySh!RGAHkG^3;z#F8V9OAlS z+rio*j!iof|6xa;n`8vLl6UL5cFvjdiKA|t^|K&)!GyPgM@|d!Yij0EuC2eOu!6sV z?_K6k75|41zYb?_?r7DzDPT1eMbK%WI>>k|Ic|D|0jTX(Zl3TaOcTLlh(+G8Ux|$< zY-be2CrfL7BRFkh8(YN}(}i*;{-ZV?qSnUa-TzV>+ZM@`$1|^Aw{DNHaPwe#T%79y zI?DdBv4Hur1BA4(itAtXsaD<-IT-MA0kh~P%XS6>pxMn+ZJ0_yw`^Gj38{nt3j6)3 zhq!NE-0@+T)Ap;Tw11H52@_z}uU;8?rC^Bxk1}X*_k$JadJKHsQbD~6)9nLpZf+w) zs~|8hzW-XMg=W;E#8;VD6jgl^hm*5&a#W$r?o1MLT0@Nq%k|VZYmKDJd!jqzF>oe; zrs0v1dT1yCMs*%XRM@iI4$SvD1pyKOP(uOKawH4Bzn&HGAr7-5Zm94j{#$~3UK%n{ ztF}arby>Zy5B+_{KJ;-RA*|7uCFRc};R(_s91^}RWPhVufBbs}rn_PN391f*N>?1~ z*u<3g37+o&FC$?IiRdDM+`kO`I&|{ zAQ#}{I*2p~sRY7L7s^i{q!;kW^Lpn~98}YVIn3)F&zw6K3ITcTVEUwNGKie|eS7x~ zd|cRIn1O zxHCcQvDcw@@%{WR)!wysz@PT+-Ysg9EWX-Je5XN5}SPx^MPFPydohnaJ^#OP@yPp59*rjGl zk7{Ifb+u5ge$1pbr_P}J`VH%=^qx#kS@KR@%k~iZ)@8cOm|3jr{LbEsOIDPdP;B46 z%s3_%ygyz0tbC2-w7kKk;P#z2+r-xJ%SLSuzKQg7Vsqj3r(X<^0B<@C#Iq4F2yJ4# z@gY93R`!9FV`Yevi>8NWb&-rB3-t-*ual0p4 zf7oDby413%r^dD7rGm^ZBPvdyF2oWkjCgfG0yrs=4&Tmwf2200e&jMJ`fk zK13z5HX^yWW@eH;jqvO6n)%^WhgXX@dfwn0>G|y0Goq4#+a-?h^XK=#O&uUPy|CRY zOz3W1YXSu{2;49@?-!Qp&N(fCMuZG$Cj) za&8O0|7W#T4uNEf&A4%+?Pwz@xL6_}{6bE84z>DjDgQ9Ag+%*9Ae#aK1c|;^gVKVq z56pxamZP1}X>}@NZL;VLo-nN!cw=|8*r)O(ls7f?nWVe9tpF~)j{rA9+}!7I!gaq` zhJ!|E0DW}wptJ=JA4V-8*oJqHgtKR1W~;aIeEvz+ThzPwTf1?~@8)Bad7xBf_Oh`J z@i+YqpEp_h2cKuK1zg@}e>1o+%kTm_6{+-P-!5=Wy{esEx!I zeT(Gh0nvKy-puX&s+}A|&vQ*r^VBs-8fdC<8f~!?{y++=9YQ>HPj)P~lRGRBpQr7S z%e!?`-kZnjzqAxIH_z|b35AROJ*5>%vugC~KhVoF843R`gsP+?lWJ79`W(^nzB*nx z{vbEGsm~WkP$htPs4{S;_XK!l{F+<8)UQ7`yrswiC#s>y$^3P))ylE6vkP%YWg@%} zZ(@RK1KQR)9PLRE%;Ap!6%6QGWLyOHdmbLb&`GRkWbA7{Q4HJ>yCjCwa5a#c(%^-j zObfks{?Aa~#sa9C<4Kz{&>QPhD3!=_$M4m69p%>$p4EYAraFU~%kTfP;N=SKME-vZ z`jXSP1Ll-fTpB^T2W1tgRuGNkaqnF59WGiY4(U<>FO_Wz@*O8zQf zv3_oy&qX6;rse-Vn!BLEnWPZ&3jE(v@t}Vcpld|F%tc13?udl8w(A?olPXynvH-UD zvS?d7k7{r03g;l(5ykp3(38c!T*pyoR$@jDFM99=Cptoft^=bvV{H8LiF72|K_E+b z{Z)VfWl_++Ks0in?4Yk|i`OOr)$IghBo70;qXXzVt%$x5b6!#LIjR8_!a5RA8&$B( z={cnXNj1V`n#*(6z$l*jVLNMy*hATa`8N9x11+MUeq(XK$MEwWZ_m}%ql|W zSUfnE*&Sbj4#H4G&n~olc`zff(}BAofSZS=-}oxGh7TqEOAcN0PAb;ODqa6LulE08)Z1)LRobG;Jvr$(lv9A|I97$ZPF zrK0muMdxzGgP7VA1D*iNSrvlaLC@M(TiwtQOC+FUa}rSJSfb@|KgzN83WYUnQ(oq$ zlXlu8>xElroXbWsXAU+fC~G~*cTL$Pks^n}LP_X=wOy#QoqDohg>>*NfZrw|lQ?`| zxhP1ja@Ifa=N6r@*3FlGC1qCF>|WQQvUu&9SHhzJuIC|>+g(Yp`zXHBX0HDtQ;uexf#Y5J2 zOspW#@vcf5g^fvUk|+^^3Lv2*bkP01e?h!BCzq?DOlG8R`D{u}-qK1KdJSVfmL((3 zt|;WJVA;+Bq!wMqau8Vmltm{MYncxc6rn2m;|LL+4kI<&)rGPKk+;wA?s! zZ37`hSzMOX*lQjGv;$6oP(+mi3mGawf|rRv@L!~jXKVpnEq3}7uwCD=qwl6ax0V}Q zz)4Fqpub{lwx69=O^UNh!rUr$#h1+KU~VZk^8tAWmGLGv+0WV%=_#^7)USjm0jR1! z`hraokeVzuFsb@s2i1{MqdItB$H=G#9byS+wb6l-iWc%EA~%41i5%yKS-pH27da%# zNlyHs*Ss5UZZ!n$!^HISkbRd!Vm>is1Ksjsg^=X_ISw)LZkFYHj_t%Uei%kk-DqL> z6)T3MKxU7iAYtC`2u;?OO!R0ddK0J?o%6Q6dTL(G(^Kuj@?c$p>?U*?k$>wN7*xSf zcqn0ZB7q$I%n~pYc!Yi_)0Pjbi=q35d3%gu_OM*Zbu{fKV7%;c0-OT)pnD<+-rOhO z78z~<3u#+AI{SktL{Mu4+`NG52gyRhyYrDU!SVPVRaj(;Mx{!`0SJv`0xLQzQECfL zLRuIN*q4=xdd}$sGoY}^(9ZQzt>cUf?N-rh0ri! zoE*yTgh_E(LMPec3lToa6hUWg+M80^-?IBJEkMo0+adD?;pt8ptDYAHLtepCZ&OSD zD9V#u7g5whE-$rdo}>MC2p{NvsijWJPzWc@}frxRxL!TX2gzh8@WO=8B}yKB4FI|V z)99Y9LEQL})CDjcb9wpqS{AG$1_lPJ+3C@2yG!~IqNCNO zW?%+;v9m3}NDdV*64Kt$284k}!NL8O4b|16w`eJwKzMX}kx!-m71_cdhc!hHQvvCD zbssgSL&5`urpqFCGSD!ybGJk}&YG&ax-PF1X4rqSxCrFbyf|DWWW~gsad42ZAwjf= z@Qi?RBOc0&Q=aH@9*I!1=Kwkr@35Mo6v6eNPqxCP9L&u6Ow^zX3Q-Am!!dIie#P;Xf??2o z=)wVjz~iRHyfqzehx;Bf6S*F;?TG^WDKTjFWf2xj@yS%A0tKo}3Id7ZmIN0_^HLzk z%Z98dXe87A^(&>5L1vo`N;;Esd$v_JtZu7=OzeG~XISbWmuS5;`viSx?=)95gv1eM zK?g?$ZX64!t)gu$zzY3%>kz@avsQ$CE{pZgH)z+UuS~XK(FwTjWbClRNUdHml1Kws&p&LW# zwHHMwI%HE^6QZKdVELm)<-ssS*pn2pv;qGG*Y91?c^=6iEKxgBve9LPU>Si6dQ<|^G!NAq-xl^ZWSm2RE%Zkd60-z=J2_O3WNaoOMBKQH znvS5#5cfp*f%zi_o9+ODB;z@RP@(sq6h>qG z$C9Xt>67jQbwS*3E+q@@ZBUTBq@r2E{{48;+O1jf*C$pk-CnaQsn@Tq2Z+vFXU+Nt zNAyw!+czie-1=D#Ai0a3oQ2=C>O*GdLeKm*mJ`koRPIl^lkEHr)ej6&+eYJyQ9WBN9y4q`5A4Grg(7c5kaFFEl1dXcRShcFSv+l zz;?5f%`s5t)bJ&&o9pawzFIqze|aJpk{k!tVj^tcDe;&7VWu2dE6-noQZoUf%^Z zoGC6YE{s_ZgM$85sCv5aB@y2X4i<4+!CMJJf^{89>g%YCn(C;EdC#d;Q0)DsshRgq z9&Hz>d17gz!k!59~L1L(C}PH1t>g>s;-;-UDD3f4|km)Ux=r0b54_XCX&%%!QEJ5^WdyFwo+1On9Qq{S|FG?z9M0YLyIYV=?4c z7vly54OsyZKMz2;N&5flH8>Psp1=^#PP42knQHEk*y@*>PY9O+RQJXi6SkbP&ACNwUD}XvF2(E zCo(^nRO0pQ&RL`UKO5`kXPg{r-NVjag#bEaG&^EMGHQJdc%_;Ca6Q_Wr;$yF3{GcNz5rPUT1_93r0O=FXUq92rqZtxILN!k@Hm#MHk#F_Oki(D1G|gK$WZj; zfe{&UCUH;Vs({uQ6<`i%1_&S|H=0olnsWrdR3JA0{M`_AiSWkiB?c0PX;Mo#_(39+ zAtW8KSX}g{@TZ^s8O%SgLj~0df0_zN431nye_s&hLVWH3B2fmo^Rb^FG>FpA0ks@# zG(ERC2zSLna;(_>RtXex)Po^N*)zXnWE;`h3 z$AYFclpG)CI>c^IkYvDNh6w+VW$xQJbE)QWpGdH=sh{<6=QQUnP}B;POl)ym$?j93 zZ{0UTsL}Oh_4zW+X#_Or*nZ4`4i^n;Ki{1>2F1%2280s8i< z`jl?C!q#9yNP_*ikNsKHP!OaHqm1qlp7|j>^nQ!cP`dk7MHD{ekZld8B-lz3-VCU} zi7l>*C}DR=AkihlS!!fLPc+#SLzU-bIx+3J3ed_cw^5XbO2FxW?&|?Rfgt$KF5!>0 zkNebTMM>c-+;WiI176fZP=>Jr&OrBRP$u>%`8^?}al|MF|4sRf{D_Ox@lc&zfbv@Y zG(Zmmx}qIV2Qp31FbsPUL&q&?I*iea&xlxA!jl~9FOTjR(^ElV+i8;7DUJ_5HoSc) zK$kkgbAgF2)Qm%H5%g=IWQ5e+>bA}@PK2?a=Yma^)w5XvE(_%ca|Y<@ zeW6s7&Yhj^AO3u}e@x#R=d!JxAQ5%5biU!d3l4aoGA%jz#>l+~56+?RC3_)3ZIiN+ z(j!bBks#T}aJeMs$4*l#24U#5ahD_QFaU4#iycTC{2$ff63;l0=7hl@108gFIHf$pKejHDHwVWbaGR5TjRJ-~?{ zGvh{euCyO>WW=;xiE@PTA(`1BnV=stXhf=j7c$t|u9T*?;lQh9L>5u+O<&sKlKhRSEATSingEBwc1Of6LKdZIE1yKSdRPoG}ar$~DKZ`8}&EVSrwYJ*_ZqQ9O!z+mxW zD*q2TDQT$Ya=yyNyb10Id<0@}+pzb|YRm8}8p>e`SVJ87*e|A^85kG<=3zBbt7Y&b7~?R&+{=(0XFWPef!vUm!@uDmLMCu#uSx_F zxXpxFqV9jz(sR@}OuasP&D=@D*8PTlMYi`>&%cO^QN=fgi)u`TdPPKgP^H*a|UyZr8{XzR^56XlLKVxl%}`dC{ncosF;pceqJ; zr|F>^Kes_dm@UwU`R94SX@5daIFAB!)Hx`K5tAw2R82}Tz;1THwIABx_*62k8A2^|VvlspVusUqkWJbK24*RWY#WHjkzCCB zyp{u=X@k+XgKxGjY*%dt`9F0|9(Gk~n~n-Q%&K`kDiv;Va54F25_C>R#$DMOld@q3 z6ftdc=gytd`uYdx;*}4m(tU^E&@ch5T+0o|m*-NQE<9(2=D#WwU*SxH zj#y0K%u=?7momg5Vhqi?DVcrhJqUp}U%osrG&*`bp+?5nBjP}UQu@u-+f7ZXlst*8 za^t7(-tDlfsjpw=^|7Ube>p8pktVlv-7NQ|W@0W1D#=0*a|U4bT9zDK1C-k-cXnnA z`trB(oKX7qzYBUV~>$-GqCvbZn9AE5fBZ}qr=IhrGc}vaMDenFI9njSU z&s(E6E3eyD-8qbt{ zzGvgO`QvP2+b50U#|8eZg=z1nF7Wq94CAR}G$9AY^+XuQWnQs4gS_{*l8nc|plE}}bKL_;qRy0!yAA|qWFY>Q-OC-*&$3!m@tqs9AscvuAI@ zb%Rg@A3lBhbwdbtz(~ps4RCV5Lxu84;3Tm*+snOaIVRiu`oN7p)M)d2@k_?o_`s0W zbh_;m4e6_@DL)TVhx43GGirX1s<&+-6AE{}06MryDg3to$B{(LjwyYHs@WJ6LdMu( z-o$q+=lvMeH!Q7ZXD3AFy`dMM!kIm^2F3VwPP(Nx$VAD zx_{iM7fov?A(L8Y(SMbYN1Kd56+iGkE=9eV1$rT!)H zigfCHPb{_*C;Z>ju`y+03OjC!`T;-u^AN`5uqbu4cDW08^~KP`EXE`Y=*-BB#xx;| zzJ}Kz?N>1}y1{Dt6$bjk+_Wp|M`MP`BG|6;m|_ zPMDbOVmRWrYlm!JO7q&{u(|8%db>>fgskiuOr>$YJ_(j_+uprPL_=2-XJ2>$beEy! zj~=~YA_y;F5uG4+ZW??wbHddhCZ0cl0}%^Wd|qNgrUn#6sV4k_Zua{?P-0IjIup^% zX|1}c# z8AnkOY3r&KDp%~jm4651vA@8!d_)9>UPca=+X5NtNT=0fIT&ACWvh<#H*-SW0S_-u zOw;^@=ynd`KfrhXWI|0L#ECz;DsZ*)=FNMFh6)`6gU68jo>Es2 z!AKz8PlFR{Q*?unJ7SuZF{UcMsMYs;{CF1{TpvPNK50{|D;v{i`I1*<<|bp)M})5{dP_m|$X&0z-=JXv{Xl zz@xPjlaq}Ti1&CnuMHX?!zTe^47Y;|y#k%8LA#%;40>c{Wx;{oJfX0BxpV$aXi(Aw zC`YZVCVb5hN*JfP4;;7#1tJbJ_qq5vyivExhNrQyk@8hOq0^%>eiVAfd(c_LBDSge zh=l3rs^Q}ME{pmL6^Vb6owO2cn11L^8P@ykazYvpXLVL1F}`&l!Iv++BV6wD9N>DJ zHi9oTCgxcTw2t$L<}u0r?zV#mH*R>a5Y2+4ue|+4iv8V9;X6!5b8mR?uG>r-wctkD zMT-jJs>IE+Y%6E5oIAR@=+oRyik=sBb7Cip@4(?c{|J(E%un>~GUkyfwgT$;9mFb< zh#i&j0BPuHDXV)vOec|1WO-*iSFC}j5`{`XrP5Vyh`Y)-`(`2?N60i8@ zHKijPK_>eQtc7l*}S}?McBR#LZ7)J4sC~H* zT8THn6)`~9mb0+zNJ>{D+_OJ_md0G9j3B11ThUL%Le46{lZy}tMqdf-Pz(6S$94SU z5g4gk?mcp3EB4b5A*2DQj>{SwT_GCJa+~gq3@&YK?7e=>0^0K1wfdM?6&O?e?HXo& zyf!4}b-4W?eXEegCz;52aZ;R!H3Zvx(VY?v>ot90qw zwX*m-iS90YXO)+foB+;yY61K8lYkqzP65lHSpND)dusE5v!9^S8Q2mBF1U2ji{G~> zA>RSGxJAm-#W4hU3M{Zln-1)zgx5cRTf9B@wi>XQ4+;x=R`~l4=;9*ax$M9K&Ls-i zSy-|pndjWr~xbf@@m`vz)_=9x^b*oSF2U-CAR;OMbBuuz>LNZ zY<0466MeaLBL&i zSNHt$>z&pAwMy>a-pA{&RQ2*7FL3+UC+P{SXZ9WMS)%hy`Kh7u#9Jjw|9`8PLfbD& zO1z)e+dBX+-GzYthj~G#>Vm-~le`3Q!v;u9Tvo^ox^of?OnRjwzzrK9u_W^j8Lc5& cTk$R4UMBL`o=Dj$U@x4()78&qol`;+0J1ZB4FCWD literal 0 HcmV?d00001 diff --git a/sample_script.py b/sample_script.py new file mode 100755 index 0000000..6fce400 --- /dev/null +++ b/sample_script.py @@ -0,0 +1,113 @@ +#!/usr/bin/env python2 +# -*- coding: utf-8 -*- + +import numpy as np +import matplotlib +# matplotlib.use('Agg') +import matplotlib.pyplot as plt +import sys +import os +from scipy.io import savemat +from scipy.io import loadmat +import timeit + +from DensityIntegrationUncertaintyQuantification import Density_integration_Poisson_uncertainty +from DensityIntegrationUncertaintyQuantification import Density_integration_WLS_uncertainty +from DensityIntegrationUncertaintyQuantification import Density_integration_WLS_uncertainty_weighted_average + +# Load the data: +# all loaded variables is in standard physical units: X, Y: [m], rho_x, rho_y: [kg/m^4] +data = loadmat('sample-data.mat', squeeze_me=True) + +# set the density uncertainty at the boundary points [kg/m^3] +sigma_rho_dirichlet = 0.01 + +# calculate density and uncertainty (Poisson) +rho_poisson, sigma_rho_poisson = Density_integration_Poisson_uncertainty(data['X'], data['Y'], data['mask'], + data['rho_x'], data['rho_y'], + data['dirichlet_label'], data['rho_dirichlet'], + uncertainty_quantification=True, + sigma_grad_x=data['sigma_rho_x'], sigma_grad_y=data['sigma_rho_y'], + sigma_dirichlet=sigma_rho_dirichlet) + +# calculate density and uncertainty (WLS) +rho_wls, sigma_rho_wls = Density_integration_WLS_uncertainty(data['X'], data['Y'], data['mask'], + data['rho_x'], data['rho_y'], + data['dirichlet_label'], data['rho_dirichlet'], + uncertainty_quantification=True, + sigma_grad_x=data['sigma_rho_x'], sigma_grad_y=data['sigma_rho_y'], + sigma_dirichlet=sigma_rho_dirichlet) + +# save the results to file +savemat(file_name='sample-result.mat', mdict={'X': data['X'], 'Y': data['Y'], 'rho_poisson': rho_poisson, 'sigma_rho_poisson': sigma_rho_poisson, + 'rho_wls': rho_wls, 'sigma_rho_wls': sigma_rho_wls}, long_field_names=True) + +# Plot the results +fig1 = plt.figure(1, figsize=(12,8)) +plt.figure(1) +ax1 = fig1.add_subplot(3,2,1) +ax2 = fig1.add_subplot(3,2,2) +ax3 = fig1.add_subplot(3,2,3) +ax4 = fig1.add_subplot(3,2,4) +ax5 = fig1.add_subplot(3,2,5) +ax6 = fig1.add_subplot(3,2,6) + +# plot x gradient +plt.axes(ax1) +plt.pcolor(data['X'], data['Y'], data['rho_x'], vmin=-60, vmax=60) +plt.colorbar() +plt.title('rho_x') + +# plot y gradient +plt.axes(ax2) +plt.pcolor(data['X'], data['Y'], data['rho_y'], vmin=-60, vmax=60) +plt.colorbar() +plt.title('rho_y') + +# plot density (poisson) +plt.axes(ax3) +plt.pcolor(data['X'], data['Y'], rho_poisson, vmin=1.2, vmax=1.5) +plt.colorbar() +plt.title('rho, Poisson') + +# plot density uncertainty (poisson) +plt.axes(ax4) +plt.pcolor(data['X'], data['Y'], sigma_rho_poisson, vmin=0.0, vmax=0.01) +plt.colorbar() +plt.title('sigma rho, Poisson') + +# plot density (wls) +plt.axes(ax5) +plt.pcolor(data['X'], data['Y'], rho_wls, vmin=1.2, vmax=1.5) +plt.colorbar() +plt.title('rho, WLS') + +# plot density uncertainty (wls) +plt.axes(ax6) +plt.pcolor(data['X'], data['Y'], sigma_rho_wls, vmin=0.0, vmax=0.01) +plt.colorbar() +plt.title('sigma rho, WLS') + +plt.tight_layout() + +# save plot to file +plt.savefig('sample-result.png') +plt.close() + + + + + + + + + + + + + + + + + +