From 3f81f0f0e1416588a99b291ec57bf1ba7d607124 Mon Sep 17 00:00:00 2001 From: larsbratholm Date: Tue, 7 Aug 2018 17:39:25 +0100 Subject: [PATCH 01/25] Made base representations --- qml/ml/representations/representations.py | 99 +++++++++++++++++++++++ 1 file changed, 99 insertions(+) diff --git a/qml/ml/representations/representations.py b/qml/ml/representations/representations.py index da2556ff5..d404dbb14 100644 --- a/qml/ml/representations/representations.py +++ b/qml/ml/representations/representations.py @@ -596,3 +596,102 @@ def generate_acsf(nuclear_charges, coordinates, elements = [1,6,7,8,16], nRs2 = else: return fgenerate_acsf(coordinates, nuclear_charges, elements, Rs2, Rs3, Ts, eta2, eta3, zeta, rcut, acut, natoms, descr_size) + + +class BaseRepresentation(object): + """ + Base class for representations + """ + + def fit(self, X, y = None): + """ + The fit routine is needed for scikit-learn pipelines + """ + return self + + def transform(self, **params): + """ + The transform routine is needed for scikit-learn pipelines. + Just returns the generate routine, since generate is a more + natural name. + """ + return self.generate(**params) + + def generate(self, **params): + raise NotImplementedError + + +class MolecularRepresentation(BaseRepresentation): + """ + Base class for molecular representations + """ + + representation_type = "molecular" + +class AtomicRepresentations(BaseRepresentation): + """ + Base class for molecular representations + """ + + representation_type = "atomic" + +class AtomCenteredSymmetryFunctions(AtomicRepresentation): + """ + Variant of atom-centered symmetry functions used in https://doi.org/10.1039/C7SC04934J + """ + + def __init__(self): + pass + + def generate(self): + pass + +def generate_acsf(nuclear_charges, coordinates, elements = [1,6,7,8,16], nRs2 = 3, nRs3 = 3, nTs = 3, eta2 = 1, + eta3 = 1, zeta = 1, rcut = 5, acut = 5, gradients = False): + """ + Generate the variant of atom-centered symmetry functions used in https://doi.org/10.1039/C7SC04934J + + :param nuclear_charges: List of nuclear charges. + :type nuclear_charges: numpy array + :param coordinates: Input coordinates + :type coordinates: numpy array + :param elements: list of unique nuclear charges (atom types) + :type elements: numpy array + :param nRs2: Number of gaussian basis functions in the two-body terms + :type nRs2: integer + :param nRs3: Number of gaussian basis functions in the three-body radial part + :type nRs3: integer + :param nTs: Number of basis functions in the three-body angular part + :type nTs: integer + :param eta2: Precision in the gaussian basis functions in the two-body terms + :type eta2: float + :param eta3: Precision in the gaussian basis functions in the three-body radial part + :type eta3: float + :param zeta: Precision parameter of basis functions in the three-body angular part + :type zeta: float + :param rcut: Cut-off radius of the two-body terms + :type rcut: float + :param acut: Cut-off radius of the three-body terms + :type acut: float + :param gradients: To return gradients or not + :type gradients: boolean + :return: Atom-centered symmetry functions representation + :rtype: numpy array + """ + + Rs2 = np.linspace(0, rcut, nRs2) + Rs3 = np.linspace(0, acut, nRs3) + Ts = np.linspace(0, np.pi, nTs) + n_elements = len(elements) + natoms = len(coordinates) + + descr_size = n_elements * nRs2 + (n_elements * (n_elements + 1)) // 2 * nRs3*nTs + + if gradients: + return fgenerate_acsf_and_gradients(coordinates, nuclear_charges, elements, Rs2, Rs3, + Ts, eta2, eta3, zeta, rcut, acut, natoms, descr_size) + else: + return fgenerate_acsf(coordinates, nuclear_charges, elements, Rs2, Rs3, + Ts, eta2, eta3, zeta, rcut, acut, natoms, descr_size) + + From f735085c8a770710cb3844f29addf7e03c6b6f56 Mon Sep 17 00:00:00 2001 From: larsbratholm Date: Thu, 9 Aug 2018 13:35:25 +0100 Subject: [PATCH 02/25] started CM and data class --- qml/data/dataprovider.py | 1 - qml/data/xyzdataprovider.py | 1 - qml/ml/representations/representations.py | 150 +++++++++++----------- 3 files changed, 72 insertions(+), 80 deletions(-) diff --git a/qml/data/dataprovider.py b/qml/data/dataprovider.py index 53fa92c8f..d5265f732 100644 --- a/qml/data/dataprovider.py +++ b/qml/data/dataprovider.py @@ -44,4 +44,3 @@ def get_properties(self, idx=None): def read_database(self, db_filename): self.compounds = connect(db_filename) - diff --git a/qml/data/xyzdataprovider.py b/qml/data/xyzdataprovider.py index f0eca2a9f..4bce0b494 100644 --- a/qml/data/xyzdataprovider.py +++ b/qml/data/xyzdataprovider.py @@ -40,4 +40,3 @@ def add_structures(self, xyz_filenames): print(i, xyz_filename, self.properties[i]) compound = read(xyz_filename) self.compounds.write(compound) - diff --git a/qml/ml/representations/representations.py b/qml/ml/representations/representations.py index d404dbb14..7100103f7 100644 --- a/qml/ml/representations/representations.py +++ b/qml/ml/representations/representations.py @@ -40,34 +40,6 @@ from .facsf import fgenerate_acsf, fgenerate_acsf_and_gradients -def vector_to_matrix(v): - """ Converts a representation from 1D vector to 2D square matrix. - :param v: 1D input representation. - :type v: numpy array - :return: Square matrix representation. - :rtype: numpy array - """ - - if not (np.sqrt(8*v.shape[0]+1) == int(np.sqrt(8*v.shape[0]+1))): - print("ERROR: Can not make a square matrix.") - exit(1) - - n = v.shape[0] - l = (-1 + int(np.sqrt(8*n+1)))//2 - M = np.empty((l,l)) - - index = 0 - for i in range(l): - for j in range(l): - if j > i: - continue - - M[i,j] = v[index] - M[j,i] = M[i,j] - - index += 1 - return M - def generate_coulomb_matrix(nuclear_charges, coordinates, size = 23, sorting = "row-norm"): """ Creates a Coulomb Matrix representation of a molecule. Sorting of the elements can either be done by ``sorting="row-norm"`` or ``sorting="unsorted"``. @@ -635,63 +607,85 @@ class AtomicRepresentations(BaseRepresentation): representation_type = "atomic" -class AtomCenteredSymmetryFunctions(AtomicRepresentation): - """ - Variant of atom-centered symmetry functions used in https://doi.org/10.1039/C7SC04934J - """ +#class AtomCenteredSymmetryFunctions(AtomicRepresentation): +# """ +# Variant of atom-centered symmetry functions used in https://doi.org/10.1039/C7SC04934J +# """ +# +# def __init__(self): +# pass +# +# def generate(self): +# pass - def __init__(self): - pass +class CoulombMatrix(MolecularRepresentation): + """ Coulomb Matrix representation of a molecule. + Sorting of the elements can either be done by ``sorting="row-norm"`` or ``sorting="unsorted"``. + A matrix :math:`M` is constructed with elements - def generate(self): - pass + .. math:: -def generate_acsf(nuclear_charges, coordinates, elements = [1,6,7,8,16], nRs2 = 3, nRs3 = 3, nTs = 3, eta2 = 1, - eta3 = 1, zeta = 1, rcut = 5, acut = 5, gradients = False): - """ - Generate the variant of atom-centered symmetry functions used in https://doi.org/10.1039/C7SC04934J + M_{ij} = + \\begin{cases} + \\tfrac{1}{2} Z_{i}^{2.4} & \\text{if } i = j \\\\ + \\frac{Z_{i}Z_{j}}{\\| {\\bf R}_{i} - {\\bf R}_{j}\\|} & \\text{if } i \\neq j + \\end{cases}, + + where :math:`i` and :math:`j` are atom indices, :math:`Z` is nuclear charge and + :math:`\\bf R` is the coordinate in euclidean space. + If ``sorting = 'row-norm'``, the atom indices are reordered such that + + :math:`\\sum_j M_{1j}^2 \\geq \\sum_j M_{2j}^2 \\geq ... \\geq \\sum_j M_{nj}^2` + + The upper triangular of M, including the diagonal, is concatenated to a 1D + vector representation. + + If ``sorting = 'unsorted``, the elements are sorted in the same order as the input coordinates + and nuclear charges. + + The representation is calculated using an OpenMP parallel Fortran routine. - :param nuclear_charges: List of nuclear charges. - :type nuclear_charges: numpy array - :param coordinates: Input coordinates - :type coordinates: numpy array - :param elements: list of unique nuclear charges (atom types) - :type elements: numpy array - :param nRs2: Number of gaussian basis functions in the two-body terms - :type nRs2: integer - :param nRs3: Number of gaussian basis functions in the three-body radial part - :type nRs3: integer - :param nTs: Number of basis functions in the three-body angular part - :type nTs: integer - :param eta2: Precision in the gaussian basis functions in the two-body terms - :type eta2: float - :param eta3: Precision in the gaussian basis functions in the three-body radial part - :type eta3: float - :param zeta: Precision parameter of basis functions in the three-body angular part - :type zeta: float - :param rcut: Cut-off radius of the two-body terms - :type rcut: float - :param acut: Cut-off radius of the three-body terms - :type acut: float - :param gradients: To return gradients or not - :type gradients: boolean - :return: Atom-centered symmetry functions representation - :rtype: numpy array """ - Rs2 = np.linspace(0, rcut, nRs2) - Rs3 = np.linspace(0, acut, nRs3) - Ts = np.linspace(0, np.pi, nTs) - n_elements = len(elements) - natoms = len(coordinates) + short_name = "cm" - descr_size = n_elements * nRs2 + (n_elements * (n_elements + 1)) // 2 * nRs3*nTs + def __init__(self, size=23, sorting="row-norm"): + """ + :param size: The size of the largest molecule supported by the representation + :type size: integer + :param sorting: How the atom indices are sorted ('row-norm', 'unsorted') + :type sorting: string + """ - if gradients: - return fgenerate_acsf_and_gradients(coordinates, nuclear_charges, elements, Rs2, Rs3, - Ts, eta2, eta3, zeta, rcut, acut, natoms, descr_size) - else: - return fgenerate_acsf(coordinates, nuclear_charges, elements, Rs2, Rs3, - Ts, eta2, eta3, zeta, rcut, acut, natoms, descr_size) + def generate(nuclear_charges, coordinates): + :param nuclear_charges: Nuclear charges of the atoms in the molecule + :type nuclear_charges: numpy array + :param coordinates: 3D Coordinates of the atoms in the molecule + :type coordinates: numpy array + :return: 1D representation - shape (size(size+1)/2,) + :rtype: numpy array + + # Raise warning is size too small and increase it + # Data options: Filenames, Data object + + if (sorting == "row-norm"): + return fgenerate_coulomb_matrix(nuclear_charges, \ + coordinates, size) + + elif (sorting == "unsorted"): + return fgenerate_unsorted_coulomb_matrix(nuclear_charges, \ + coordinates, size) + + else: + print("ERROR: Unknown sorting scheme requested") + raise SystemExit + +class Data(object): + """ + Temporary data class which should be replaced at some point + """ + + def __init__(self): + pass From b3c6168ca10974e3589551f930ee5f6cd85353c2 Mon Sep 17 00:00:00 2001 From: larsbratholm Date: Fri, 10 Aug 2018 11:57:01 +0100 Subject: [PATCH 03/25] Working on generate routine --- qml/ml/representations/representations.py | 1155 +++++++++++---------- 1 file changed, 615 insertions(+), 540 deletions(-) diff --git a/qml/ml/representations/representations.py b/qml/ml/representations/representations.py index 7100103f7..f24d9fd31 100644 --- a/qml/ml/representations/representations.py +++ b/qml/ml/representations/representations.py @@ -40,534 +40,534 @@ from .facsf import fgenerate_acsf, fgenerate_acsf_and_gradients -def generate_coulomb_matrix(nuclear_charges, coordinates, size = 23, sorting = "row-norm"): - """ Creates a Coulomb Matrix representation of a molecule. - Sorting of the elements can either be done by ``sorting="row-norm"`` or ``sorting="unsorted"``. - A matrix :math:`M` is constructed with elements - - .. math:: - - M_{ij} = - \\begin{cases} - \\tfrac{1}{2} Z_{i}^{2.4} & \\text{if } i = j \\\\ - \\frac{Z_{i}Z_{j}}{\\| {\\bf R}_{i} - {\\bf R}_{j}\\|} & \\text{if } i \\neq j - \\end{cases}, - - where :math:`i` and :math:`j` are atom indices, :math:`Z` is nuclear charge and - :math:`\\bf R` is the coordinate in euclidean space. - If ``sorting = 'row-norm'``, the atom indices are reordered such that - - :math:`\\sum_j M_{1j}^2 \\geq \\sum_j M_{2j}^2 \\geq ... \\geq \\sum_j M_{nj}^2` - - The upper triangular of M, including the diagonal, is concatenated to a 1D - vector representation. - - If ``sorting = 'unsorted``, the elements are sorted in the same order as the input coordinates - and nuclear charges. - - The representation is calculated using an OpenMP parallel Fortran routine. - - :param nuclear_charges: Nuclear charges of the atoms in the molecule - :type nuclear_charges: numpy array - :param coordinates: 3D Coordinates of the atoms in the molecule - :type coordinates: numpy array - :param size: The size of the largest molecule supported by the representation - :type size: integer - :param sorting: How the atom indices are sorted ('row-norm', 'unsorted') - :type sorting: string - - :return: 1D representation - shape (size(size+1)/2,) - :rtype: numpy array - """ - - if (sorting == "row-norm"): - return fgenerate_coulomb_matrix(nuclear_charges, \ - coordinates, size) - - elif (sorting == "unsorted"): - return fgenerate_unsorted_coulomb_matrix(nuclear_charges, \ - coordinates, size) - - else: - print("ERROR: Unknown sorting scheme requested") - raise SystemExit - -def generate_atomic_coulomb_matrix(nuclear_charges, coordinates, size = 23, sorting = "distance", - central_cutoff = 1e6, central_decay = -1, interaction_cutoff = 1e6, interaction_decay = -1, - indices = None): - """ Creates a Coulomb Matrix representation of the local environment of a central atom. - For each central atom :math:`k`, a matrix :math:`M` is constructed with elements - - .. math:: - - M_{ij}(k) = - \\begin{cases} - \\tfrac{1}{2} Z_{i}^{2.4} \\cdot f_{ik}^2 & \\text{if } i = j \\\\ - \\frac{Z_{i}Z_{j}}{\\| {\\bf R}_{i} - {\\bf R}_{j}\\|} \\cdot f_{ik}f_{jk}f_{ij} & \\text{if } i \\neq j - \\end{cases}, - - where :math:`i`, :math:`j` and :math:`k` are atom indices, :math:`Z` is nuclear charge and - :math:`\\bf R` is the coordinate in euclidean space. - - :math:`f_{ij}` is a function that masks long range effects: - - .. math:: - - f_{ij} = - \\begin{cases} - 1 & \\text{if } \\|{\\bf R}_{i} - {\\bf R}_{j} \\| \\leq r - \Delta r \\\\ - \\tfrac{1}{2} \\big(1 + \\cos\\big(\\pi \\tfrac{\\|{\\bf R}_{i} - {\\bf R}_{j} \\| - - r + \Delta r}{\Delta r} \\big)\\big) - & \\text{if } r - \Delta r < \\|{\\bf R}_{i} - {\\bf R}_{j} \\| \\leq r - \Delta r \\\\ - 0 & \\text{if } \\|{\\bf R}_{i} - {\\bf R}_{j} \\| > r - \\end{cases}, - - where the parameters ``central_cutoff`` and ``central_decay`` corresponds to the variables - :math:`r` and :math:`\Delta r` respectively for interactions involving the central atom, - and ``interaction_cutoff`` and ``interaction_decay`` corresponds to the variables - :math:`r` and :math:`\Delta r` respectively for interactions not involving the central atom. - - if ``sorting = 'row-norm'``, the atom indices are ordered such that - - :math:`\\sum_j M_{1j}(k)^2 \\geq \\sum_j M_{2j}(k)^2 \\geq ... \\geq \\sum_j M_{nj}(k)^2` - - if ``sorting = 'distance'``, the atom indices are ordered such that - - .. math:: - - \\|{\\bf R}_{1} - {\\bf R}_{k}\\| \\leq \\|{\\bf R}_{2} - {\\bf R}_{k}\\| - \\leq ... \\leq \\|{\\bf R}_{n} - {\\bf R}_{k}\\| - - The upper triangular of M, including the diagonal, is concatenated to a 1D - vector representation. - - The representation can be calculated for a subset by either specifying - ``indices = [0,1,...]``, where :math:`[0,1,...]` are the requested atom indices, - or by specifying ``indices = 'C'`` to only calculate central carbon atoms. - - The representation is calculated using an OpenMP parallel Fortran routine. - - :param nuclear_charges: Nuclear charges of the atoms in the molecule - :type nuclear_charges: numpy array - :param coordinates: 3D Coordinates of the atoms in the molecule - :type coordinates: numpy array - :param size: The size of the largest molecule supported by the representation - :type size: integer - :param sorting: How the atom indices are sorted ('row-norm', 'distance') - :type sorting: string - :param central_cutoff: The distance from the central atom, where the coulomb interaction - element will be zero - :type central_cutoff: float - :param central_decay: The distance over which the the coulomb interaction decays from full to none - :type central_decay: float - :param interaction_cutoff: The distance between two non-central atom, where the coulomb interaction - element will be zero - :type interaction_cutoff: float - :param interaction_decay: The distance over which the the coulomb interaction decays from full to none - :type interaction_decay: float - :param indices: Subset indices or atomtype - :type indices: Nonetype/array/string - - - :return: nD representation - shape (:math:`N_{atoms}`, size(size+1)/2) - :rtype: numpy array - """ - - - if indices == None: - nindices = len(nuclear_charges) - indices = np.arange(1,1+nindices, 1, dtype = int) - elif type("") == type(indices): - if indices in NUCLEAR_CHARGE: - indices = np.where(nuclear_charges == NUCLEAR_CHARGE[indices])[0] + 1 - nindices = indices.size - if nindices == 0: - return np.zeros((0,0)) - - else: - print("ERROR: Unknown value %s given for 'indices' variable" % indices) - raise SystemExit - else: - indices = np.asarray(indices, dtype = int) + 1 - nindices = indices.size - - - if (sorting == "row-norm"): - return fgenerate_local_coulomb_matrix(indices, nindices, nuclear_charges, - coordinates, nuclear_charges.size, size, - central_cutoff, central_decay, interaction_cutoff, interaction_decay) - - elif (sorting == "distance"): - return fgenerate_atomic_coulomb_matrix(indices, nindices, nuclear_charges, - coordinates, nuclear_charges.size, size, - central_cutoff, central_decay, interaction_cutoff, interaction_decay) - - else: - print("ERROR: Unknown sorting scheme requested") - raise SystemExit - -def generate_eigenvalue_coulomb_matrix(nuclear_charges, coordinates, size = 23): - """ Creates an eigenvalue Coulomb Matrix representation of a molecule. - A matrix :math:`M` is constructed with elements - - .. math:: - - M_{ij} = - \\begin{cases} - \\tfrac{1}{2} Z_{i}^{2.4} & \\text{if } i = j \\\\ - \\frac{Z_{i}Z_{j}}{\\| {\\bf R}_{i} - {\\bf R}_{j}\\|} & \\text{if } i \\neq j - \\end{cases}, - - where :math:`i` and :math:`j` are atom indices, :math:`Z` is nuclear charge and - :math:`\\bf R` is the coordinate in euclidean space. - The molecular representation of the molecule is then the sorted eigenvalues of M. - The representation is calculated using an OpenMP parallel Fortran routine. - - :param nuclear_charges: Nuclear charges of the atoms in the molecule - :type nuclear_charges: numpy array - :param coordinates: 3D Coordinates of the atoms in the molecule - :type coordinates: numpy array - :param size: The size of the largest molecule supported by the representation - :type size: integer - - :return: 1D representation - shape (size, ) - :rtype: numpy array - """ - return fgenerate_eigenvalue_coulomb_matrix(nuclear_charges, - coordinates, size) - -def generate_bob(nuclear_charges, coordinates, atomtypes, size=23, asize = {"O":3, "C":7, "N":3, "H":16, "S":1}): - """ Creates a Bag of Bonds (BOB) representation of a molecule. - The representation expands on the coulomb matrix representation. - For each element a bag (vector) is constructed for self interactions - (e.g. ('C', 'H', 'O')). - For each element pair a bag is constructed for interatomic interactions - (e.g. ('CC', 'CH', 'CO', 'HH', 'HO', 'OO')), sorted by value. - The self interaction of element :math:`I` is given by - - :math:`\\tfrac{1}{2} Z_{I}^{2.4}`, - - with :math:`Z_{i}` being the nuclear charge of element :math:`i` - The interaction between atom :math:`i` of element :math:`I` and - atom :math:`j` of element :math:`J` is given by - - :math:`\\frac{Z_{I}Z_{J}}{\\| {\\bf R}_{i} - {\\bf R}_{j}\\|}` - - with :math:`R_{i}` being the euclidean coordinate of atom :math:`i`. - The sorted bags are concatenated to an 1D vector representation. - The representation is calculated using an OpenMP parallel Fortran routine. - - :param nuclear_charges: Nuclear charges of the atoms in the molecule - :type nuclear_charges: numpy array - :param coordinates: 3D Coordinates of the atoms in the molecule - :type coordinates: numpy array - :param size: The maximum number of atoms in the representation - :type size: integer - :param asize: The maximum number of atoms of each element type supported by the representation - :type asize: dictionary - - :return: 1D representation - :rtype: numpy array - """ - - n = 0 - atoms = sorted(asize, key=asize.get) - nmax = [asize[key] for key in atoms] - ids = np.zeros(len(nmax), dtype=int) - for i, (key, value) in enumerate(zip(atoms,nmax)): - n += value * (1+value) - ids[i] = NUCLEAR_CHARGE[key] - for j in range(i): - v = nmax[j] - n += 2 * value * v - n /= 2 - - return fgenerate_bob(nuclear_charges, coordinates, nuclear_charges, ids, nmax, n) - -def get_slatm_mbtypes(nuclear_charges, pbc='000'): - """ - Get the list of minimal types of many-body terms in a dataset. This resulting list - is necessary as input in the ``generate_slatm_representation()`` function. - - :param nuclear_charges: A list of the nuclear charges for each compound in the dataset. - :type nuclear_charges: list of numpy arrays - :param pbc: periodic boundary condition along x,y,z direction, defaulted to '000', i.e., molecule - :type pbc: string - :return: A list containing the types of many-body terms. - :rtype: list - """ - - zs = nuclear_charges - - nm = len(zs) - zsmax = set() - nas = [] - zs_ravel = [] - for zsi in zs: - na = len(zsi); nas.append(na) - zsil = list(zsi); zs_ravel += zsil - zsmax.update( zsil ) - - zsmax = np.array( list(zsmax) ) - nass = [] - for i in range(nm): - zsi = np.array(zs[i],np.int) - nass.append( [ (zi == zsi).sum() for zi in zsmax ] ) - - nzmax = np.max(np.array(nass), axis=0) - nzmax_u = [] - if pbc != '000': - # the PBC will introduce new many-body terms, so set - # nzmax to 3 if it's less than 3 - for nzi in nzmax: - if nzi <= 2: - nzi = 3 - nzmax_u.append(nzi) - nzmax = nzmax_u - - boas = [ [zi,] for zi in zsmax ] - bops = [ [zi,zi] for zi in zsmax ] + list( itl.combinations(zsmax,2) ) - - bots = [] - for i in zsmax: - for bop in bops: - j,k = bop - tas = [ [i,j,k], [i,k,j], [j,i,k] ] - for tasi in tas: - if (tasi not in bots) and (tasi[::-1] not in bots): - nzsi = [ (zj == tasi).sum() for zj in zsmax ] - if np.all(nzsi <= nzmax): - bots.append( tasi ) - mbtypes = boas + bops + bots - - return mbtypes #, np.array(zs_ravel), np.array(nas) - -def generate_slatm(coordinates, nuclear_charges, mbtypes, - unit_cell=None, local=False, sigmas=[0.05,0.05], dgrids=[0.03,0.03], - rcut=4.8, alchemy=False, pbc='000', rpower=6): - """ - Generate Spectrum of London and Axillrod-Teller-Muto potential (SLATM) representation. - Both global (``local=False``) and local (``local=True``) SLATM are available. - - A version that works for periodic boundary conditions will be released soon. - - NOTE: You will need to run the ``get_slatm_mbtypes()`` function to get the ``mbtypes`` input (or generate it manually). - - :param coordinates: Input coordinates - :type coordinates: numpy array - :param nuclear_charges: List of nuclear charges. - :type nuclear_charges: numpy array - :param mbtypes: Many-body types for the whole dataset, including 1-, 2- and 3-body types. Could be obtained by calling ``get_slatm_mbtypes()``. - :type mbtypes: list - :param local: Generate a local representation. Defaulted to False (i.e., global representation); otherwise, atomic version. - :type local: bool - :param sigmas: Controlling the width of Gaussian smearing function for 2- and 3-body parts, defaulted to [0.05,0.05], usually these do not need to be adjusted. - :type sigmas: list - :param dgrids: The interval between two sampled internuclear distances and angles, defaulted to [0.03,0.03], no need for change, compromised for speed and accuracy. - :type dgrids: list - :param rcut: Cut-off radius, defaulted to 4.8 Angstrom. - :type rcut: float - :param alchemy: Swith to use the alchemy version of SLATM. (default=False) - :type alchemy: bool - :param pbc: defaulted to '000', meaning it's a molecule; the three digits in the string corresponds to x,y,z direction - :type pbc: string - :param rpower: The power of R in 2-body potential, defaulted to London potential (=6). - :type rpower: float - :return: 1D SLATM representation - :rtype: numpy array - """ - - c = unit_cell - iprt=False - if c is None: - c = np.array([[1,0,0],[0,1,0],[0,0,1]]) - - if pbc != '000': - # print(' -- handling systems with periodic boundary condition') - assert c != None, 'ERROR: Please specify unit cell for SLATM' - # ======================================================================= - # PBC may introduce new many-body terms, so at the stage of get statistics - # info from db, we've already considered this point by letting maximal number - # of nuclear charges being 3. - # ======================================================================= - - zs = nuclear_charges - na = len(zs) - coords = coordinates - obj = [ zs, coords, c ] - - iloc = local - - if iloc: - mbs = [] - X2Ns = [] - for ia in range(na): - # if iprt: print ' -- ia = ', ia + 1 - n1 = 0; n2 = 0; n3 = 0 - mbs_ia = np.zeros(0) - icount = 0 - for mbtype in mbtypes: - if len(mbtype) == 1: - mbsi = get_boa(mbtype[0], np.array([zs[ia],])) - #print ' -- mbsi = ', mbsi - if alchemy: - n1 = 1 - n1_0 = mbs_ia.shape[0] - if n1_0 == 0: - mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) - elif n1_0 == 1: - mbs_ia += mbsi - else: - raise '#ERROR' - else: - n1 += len(mbsi) - mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) - elif len(mbtype) == 2: - #print ' 001, pbc = ', pbc - mbsi = get_sbop(mbtype, obj, iloc=iloc, ia=ia, \ - sigma=sigmas[0], dgrid=dgrids[0], rcut=rcut, \ - pbc=pbc, rpower=rpower) - mbsi *= 0.5 # only for the two-body parts, local rpst - #print ' 002' - if alchemy: - n2 = len(mbsi) - n2_0 = mbs_ia.shape[0] - if n2_0 == n1: - mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) - elif n2_0 == n1 + n2: - t = mbs_ia[n1:n1+n2] + mbsi - mbs_ia[n1:n1+n2] = t - else: - raise '#ERROR' - else: - n2 += len(mbsi) - mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) - else: # len(mbtype) == 3: - mbsi = get_sbot(mbtype, obj, iloc=iloc, ia=ia, \ - sigma=sigmas[1], dgrid=dgrids[1], rcut=rcut, pbc=pbc) - - if alchemy: - n3 = len(mbsi) - n3_0 = mbs_ia.shape[0] - if n3_0 == n1 + n2: - mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) - elif n3_0 == n1 + n2 + n3: - t = mbs_ia[n1+n2:n1+n2+n3] + mbsi - mbs_ia[n1+n2:n1+n2+n3] = t - else: - raise '#ERROR' - else: - n3 += len(mbsi) - mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) - - mbs.append( mbs_ia ) - X2N = [n1,n2,n3]; - if X2N not in X2Ns: - X2Ns.append(X2N) - assert len(X2Ns) == 1, '#ERROR: multiple `X2N ???' - else: - n1 = 0; n2 = 0; n3 = 0 - mbs = np.zeros(0) - for mbtype in mbtypes: - if len(mbtype) == 1: - mbsi = get_boa(mbtype[0], zs) - if alchemy: - n1 = 1 - n1_0 = mbs.shape[0] - if n1_0 == 0: - mbs = np.concatenate( (mbs, [sum(mbsi)] ), axis=0 ) - elif n1_0 == 1: - mbs += sum(mbsi ) - else: - raise '#ERROR' - else: - n1 += len(mbsi) - mbs = np.concatenate( (mbs, mbsi), axis=0 ) - elif len(mbtype) == 2: - mbsi = get_sbop(mbtype, obj, sigma=sigmas[0], \ - dgrid=dgrids[0], rcut=rcut, rpower=rpower) - - - if alchemy: - n2 = len(mbsi) - n2_0 = mbs.shape[0] - if n2_0 == n1: - mbs = np.concatenate( (mbs, mbsi), axis=0 ) - elif n2_0 == n1 + n2: - t = mbs[n1:n1+n2] + mbsi - mbs[n1:n1+n2] = t - else: - raise '#ERROR' - else: - n2 += len(mbsi) - mbs = np.concatenate( (mbs, mbsi), axis=0 ) - else: # len(mbtype) == 3: - mbsi = get_sbot(mbtype, obj, sigma=sigmas[1], \ - dgrid=dgrids[1], rcut=rcut) - - if alchemy: - n3 = len(mbsi) - n3_0 = mbs.shape[0] - if n3_0 == n1 + n2: - mbs = np.concatenate( (mbs, mbsi), axis=0 ) - elif n3_0 == n1 + n2 + n3: - t = mbs[n1+n2:n1+n2+n3] + mbsi - mbs[n1+n2:n1+n2+n3] = t - else: - raise '#ERROR' - else: - n3 += len(mbsi) - mbs = np.concatenate( (mbs, mbsi), axis=0 ) - - return mbs - -def generate_acsf(nuclear_charges, coordinates, elements = [1,6,7,8,16], nRs2 = 3, nRs3 = 3, nTs = 3, eta2 = 1, - eta3 = 1, zeta = 1, rcut = 5, acut = 5, gradients = False): - """ - Generate the variant of atom-centered symmetry functions used in https://doi.org/10.1039/C7SC04934J - - :param nuclear_charges: List of nuclear charges. - :type nuclear_charges: numpy array - :param coordinates: Input coordinates - :type coordinates: numpy array - :param elements: list of unique nuclear charges (atom types) - :type elements: numpy array - :param nRs2: Number of gaussian basis functions in the two-body terms - :type nRs2: integer - :param nRs3: Number of gaussian basis functions in the three-body radial part - :type nRs3: integer - :param nTs: Number of basis functions in the three-body angular part - :type nTs: integer - :param eta2: Precision in the gaussian basis functions in the two-body terms - :type eta2: float - :param eta3: Precision in the gaussian basis functions in the three-body radial part - :type eta3: float - :param zeta: Precision parameter of basis functions in the three-body angular part - :type zeta: float - :param rcut: Cut-off radius of the two-body terms - :type rcut: float - :param acut: Cut-off radius of the three-body terms - :type acut: float - :param gradients: To return gradients or not - :type gradients: boolean - :return: Atom-centered symmetry functions representation - :rtype: numpy array - """ - - Rs2 = np.linspace(0, rcut, nRs2) - Rs3 = np.linspace(0, acut, nRs3) - Ts = np.linspace(0, np.pi, nTs) - n_elements = len(elements) - natoms = len(coordinates) - - descr_size = n_elements * nRs2 + (n_elements * (n_elements + 1)) // 2 * nRs3*nTs - - if gradients: - return fgenerate_acsf_and_gradients(coordinates, nuclear_charges, elements, Rs2, Rs3, - Ts, eta2, eta3, zeta, rcut, acut, natoms, descr_size) - else: - return fgenerate_acsf(coordinates, nuclear_charges, elements, Rs2, Rs3, - Ts, eta2, eta3, zeta, rcut, acut, natoms, descr_size) +#def generate_coulomb_matrix(nuclear_charges, coordinates, size = 23, sorting = "row-norm"): +# """ Creates a Coulomb Matrix representation of a molecule. +# Sorting of the elements can either be done by ``sorting="row-norm"`` or ``sorting="unsorted"``. +# A matrix :math:`M` is constructed with elements +# +# .. math:: +# +# M_{ij} = +# \\begin{cases} +# \\tfrac{1}{2} Z_{i}^{2.4} & \\text{if } i = j \\\\ +# \\frac{Z_{i}Z_{j}}{\\| {\\bf R}_{i} - {\\bf R}_{j}\\|} & \\text{if } i \\neq j +# \\end{cases}, +# +# where :math:`i` and :math:`j` are atom indices, :math:`Z` is nuclear charge and +# :math:`\\bf R` is the coordinate in euclidean space. +# If ``sorting = 'row-norm'``, the atom indices are reordered such that +# +# :math:`\\sum_j M_{1j}^2 \\geq \\sum_j M_{2j}^2 \\geq ... \\geq \\sum_j M_{nj}^2` +# +# The upper triangular of M, including the diagonal, is concatenated to a 1D +# vector representation. +# +# If ``sorting = 'unsorted``, the elements are sorted in the same order as the input coordinates +# and nuclear charges. +# +# The representation is calculated using an OpenMP parallel Fortran routine. +# +# :param nuclear_charges: Nuclear charges of the atoms in the molecule +# :type nuclear_charges: numpy array +# :param coordinates: 3D Coordinates of the atoms in the molecule +# :type coordinates: numpy array +# :param size: The size of the largest molecule supported by the representation +# :type size: integer +# :param sorting: How the atom indices are sorted ('row-norm', 'unsorted') +# :type sorting: string +# +# :return: 1D representation - shape (size(size+1)/2,) +# :rtype: numpy array +# """ +# +# if (sorting == "row-norm"): +# return fgenerate_coulomb_matrix(nuclear_charges, \ +# coordinates, size) +# +# elif (sorting == "unsorted"): +# return fgenerate_unsorted_coulomb_matrix(nuclear_charges, \ +# coordinates, size) +# +# else: +# print("ERROR: Unknown sorting scheme requested") +# raise SystemExit +# +#def generate_atomic_coulomb_matrix(nuclear_charges, coordinates, size = 23, sorting = "distance", +# central_cutoff = 1e6, central_decay = -1, interaction_cutoff = 1e6, interaction_decay = -1, +# indices = None): +# """ Creates a Coulomb Matrix representation of the local environment of a central atom. +# For each central atom :math:`k`, a matrix :math:`M` is constructed with elements +# +# .. math:: +# +# M_{ij}(k) = +# \\begin{cases} +# \\tfrac{1}{2} Z_{i}^{2.4} \\cdot f_{ik}^2 & \\text{if } i = j \\\\ +# \\frac{Z_{i}Z_{j}}{\\| {\\bf R}_{i} - {\\bf R}_{j}\\|} \\cdot f_{ik}f_{jk}f_{ij} & \\text{if } i \\neq j +# \\end{cases}, +# +# where :math:`i`, :math:`j` and :math:`k` are atom indices, :math:`Z` is nuclear charge and +# :math:`\\bf R` is the coordinate in euclidean space. +# +# :math:`f_{ij}` is a function that masks long range effects: +# +# .. math:: +# +# f_{ij} = +# \\begin{cases} +# 1 & \\text{if } \\|{\\bf R}_{i} - {\\bf R}_{j} \\| \\leq r - \Delta r \\\\ +# \\tfrac{1}{2} \\big(1 + \\cos\\big(\\pi \\tfrac{\\|{\\bf R}_{i} - {\\bf R}_{j} \\| +# - r + \Delta r}{\Delta r} \\big)\\big) +# & \\text{if } r - \Delta r < \\|{\\bf R}_{i} - {\\bf R}_{j} \\| \\leq r - \Delta r \\\\ +# 0 & \\text{if } \\|{\\bf R}_{i} - {\\bf R}_{j} \\| > r +# \\end{cases}, +# +# where the parameters ``central_cutoff`` and ``central_decay`` corresponds to the variables +# :math:`r` and :math:`\Delta r` respectively for interactions involving the central atom, +# and ``interaction_cutoff`` and ``interaction_decay`` corresponds to the variables +# :math:`r` and :math:`\Delta r` respectively for interactions not involving the central atom. +# +# if ``sorting = 'row-norm'``, the atom indices are ordered such that +# +# :math:`\\sum_j M_{1j}(k)^2 \\geq \\sum_j M_{2j}(k)^2 \\geq ... \\geq \\sum_j M_{nj}(k)^2` +# +# if ``sorting = 'distance'``, the atom indices are ordered such that +# +# .. math:: +# +# \\|{\\bf R}_{1} - {\\bf R}_{k}\\| \\leq \\|{\\bf R}_{2} - {\\bf R}_{k}\\| +# \\leq ... \\leq \\|{\\bf R}_{n} - {\\bf R}_{k}\\| +# +# The upper triangular of M, including the diagonal, is concatenated to a 1D +# vector representation. +# +# The representation can be calculated for a subset by either specifying +# ``indices = [0,1,...]``, where :math:`[0,1,...]` are the requested atom indices, +# or by specifying ``indices = 'C'`` to only calculate central carbon atoms. +# +# The representation is calculated using an OpenMP parallel Fortran routine. +# +# :param nuclear_charges: Nuclear charges of the atoms in the molecule +# :type nuclear_charges: numpy array +# :param coordinates: 3D Coordinates of the atoms in the molecule +# :type coordinates: numpy array +# :param size: The size of the largest molecule supported by the representation +# :type size: integer +# :param sorting: How the atom indices are sorted ('row-norm', 'distance') +# :type sorting: string +# :param central_cutoff: The distance from the central atom, where the coulomb interaction +# element will be zero +# :type central_cutoff: float +# :param central_decay: The distance over which the the coulomb interaction decays from full to none +# :type central_decay: float +# :param interaction_cutoff: The distance between two non-central atom, where the coulomb interaction +# element will be zero +# :type interaction_cutoff: float +# :param interaction_decay: The distance over which the the coulomb interaction decays from full to none +# :type interaction_decay: float +# :param indices: Subset indices or atomtype +# :type indices: Nonetype/array/string +# +# +# :return: nD representation - shape (:math:`N_{atoms}`, size(size+1)/2) +# :rtype: numpy array +# """ +# +# +# if indices == None: +# nindices = len(nuclear_charges) +# indices = np.arange(1,1+nindices, 1, dtype = int) +# elif type("") == type(indices): +# if indices in NUCLEAR_CHARGE: +# indices = np.where(nuclear_charges == NUCLEAR_CHARGE[indices])[0] + 1 +# nindices = indices.size +# if nindices == 0: +# return np.zeros((0,0)) +# +# else: +# print("ERROR: Unknown value %s given for 'indices' variable" % indices) +# raise SystemExit +# else: +# indices = np.asarray(indices, dtype = int) + 1 +# nindices = indices.size +# +# +# if (sorting == "row-norm"): +# return fgenerate_local_coulomb_matrix(indices, nindices, nuclear_charges, +# coordinates, nuclear_charges.size, size, +# central_cutoff, central_decay, interaction_cutoff, interaction_decay) +# +# elif (sorting == "distance"): +# return fgenerate_atomic_coulomb_matrix(indices, nindices, nuclear_charges, +# coordinates, nuclear_charges.size, size, +# central_cutoff, central_decay, interaction_cutoff, interaction_decay) +# +# else: +# print("ERROR: Unknown sorting scheme requested") +# raise SystemExit +# +#def generate_eigenvalue_coulomb_matrix(nuclear_charges, coordinates, size = 23): +# """ Creates an eigenvalue Coulomb Matrix representation of a molecule. +# A matrix :math:`M` is constructed with elements +# +# .. math:: +# +# M_{ij} = +# \\begin{cases} +# \\tfrac{1}{2} Z_{i}^{2.4} & \\text{if } i = j \\\\ +# \\frac{Z_{i}Z_{j}}{\\| {\\bf R}_{i} - {\\bf R}_{j}\\|} & \\text{if } i \\neq j +# \\end{cases}, +# +# where :math:`i` and :math:`j` are atom indices, :math:`Z` is nuclear charge and +# :math:`\\bf R` is the coordinate in euclidean space. +# The molecular representation of the molecule is then the sorted eigenvalues of M. +# The representation is calculated using an OpenMP parallel Fortran routine. +# +# :param nuclear_charges: Nuclear charges of the atoms in the molecule +# :type nuclear_charges: numpy array +# :param coordinates: 3D Coordinates of the atoms in the molecule +# :type coordinates: numpy array +# :param size: The size of the largest molecule supported by the representation +# :type size: integer +# +# :return: 1D representation - shape (size, ) +# :rtype: numpy array +# """ +# return fgenerate_eigenvalue_coulomb_matrix(nuclear_charges, +# coordinates, size) +# +#def generate_bob(nuclear_charges, coordinates, atomtypes, size=23, asize = {"O":3, "C":7, "N":3, "H":16, "S":1}): +# """ Creates a Bag of Bonds (BOB) representation of a molecule. +# The representation expands on the coulomb matrix representation. +# For each element a bag (vector) is constructed for self interactions +# (e.g. ('C', 'H', 'O')). +# For each element pair a bag is constructed for interatomic interactions +# (e.g. ('CC', 'CH', 'CO', 'HH', 'HO', 'OO')), sorted by value. +# The self interaction of element :math:`I` is given by +# +# :math:`\\tfrac{1}{2} Z_{I}^{2.4}`, +# +# with :math:`Z_{i}` being the nuclear charge of element :math:`i` +# The interaction between atom :math:`i` of element :math:`I` and +# atom :math:`j` of element :math:`J` is given by +# +# :math:`\\frac{Z_{I}Z_{J}}{\\| {\\bf R}_{i} - {\\bf R}_{j}\\|}` +# +# with :math:`R_{i}` being the euclidean coordinate of atom :math:`i`. +# The sorted bags are concatenated to an 1D vector representation. +# The representation is calculated using an OpenMP parallel Fortran routine. +# +# :param nuclear_charges: Nuclear charges of the atoms in the molecule +# :type nuclear_charges: numpy array +# :param coordinates: 3D Coordinates of the atoms in the molecule +# :type coordinates: numpy array +# :param size: The maximum number of atoms in the representation +# :type size: integer +# :param asize: The maximum number of atoms of each element type supported by the representation +# :type asize: dictionary +# +# :return: 1D representation +# :rtype: numpy array +# """ +# +# n = 0 +# atoms = sorted(asize, key=asize.get) +# nmax = [asize[key] for key in atoms] +# ids = np.zeros(len(nmax), dtype=int) +# for i, (key, value) in enumerate(zip(atoms,nmax)): +# n += value * (1+value) +# ids[i] = NUCLEAR_CHARGE[key] +# for j in range(i): +# v = nmax[j] +# n += 2 * value * v +# n /= 2 +# +# return fgenerate_bob(nuclear_charges, coordinates, nuclear_charges, ids, nmax, n) +# +#def get_slatm_mbtypes(nuclear_charges, pbc='000'): +# """ +# Get the list of minimal types of many-body terms in a dataset. This resulting list +# is necessary as input in the ``generate_slatm_representation()`` function. +# +# :param nuclear_charges: A list of the nuclear charges for each compound in the dataset. +# :type nuclear_charges: list of numpy arrays +# :param pbc: periodic boundary condition along x,y,z direction, defaulted to '000', i.e., molecule +# :type pbc: string +# :return: A list containing the types of many-body terms. +# :rtype: list +# """ +# +# zs = nuclear_charges +# +# nm = len(zs) +# zsmax = set() +# nas = [] +# zs_ravel = [] +# for zsi in zs: +# na = len(zsi); nas.append(na) +# zsil = list(zsi); zs_ravel += zsil +# zsmax.update( zsil ) +# +# zsmax = np.array( list(zsmax) ) +# nass = [] +# for i in range(nm): +# zsi = np.array(zs[i],np.int) +# nass.append( [ (zi == zsi).sum() for zi in zsmax ] ) +# +# nzmax = np.max(np.array(nass), axis=0) +# nzmax_u = [] +# if pbc != '000': +# # the PBC will introduce new many-body terms, so set +# # nzmax to 3 if it's less than 3 +# for nzi in nzmax: +# if nzi <= 2: +# nzi = 3 +# nzmax_u.append(nzi) +# nzmax = nzmax_u +# +# boas = [ [zi,] for zi in zsmax ] +# bops = [ [zi,zi] for zi in zsmax ] + list( itl.combinations(zsmax,2) ) +# +# bots = [] +# for i in zsmax: +# for bop in bops: +# j,k = bop +# tas = [ [i,j,k], [i,k,j], [j,i,k] ] +# for tasi in tas: +# if (tasi not in bots) and (tasi[::-1] not in bots): +# nzsi = [ (zj == tasi).sum() for zj in zsmax ] +# if np.all(nzsi <= nzmax): +# bots.append( tasi ) +# mbtypes = boas + bops + bots +# +# return mbtypes #, np.array(zs_ravel), np.array(nas) +# +#def generate_slatm(coordinates, nuclear_charges, mbtypes, +# unit_cell=None, local=False, sigmas=[0.05,0.05], dgrids=[0.03,0.03], +# rcut=4.8, alchemy=False, pbc='000', rpower=6): +# """ +# Generate Spectrum of London and Axillrod-Teller-Muto potential (SLATM) representation. +# Both global (``local=False``) and local (``local=True``) SLATM are available. +# +# A version that works for periodic boundary conditions will be released soon. +# +# NOTE: You will need to run the ``get_slatm_mbtypes()`` function to get the ``mbtypes`` input (or generate it manually). +# +# :param coordinates: Input coordinates +# :type coordinates: numpy array +# :param nuclear_charges: List of nuclear charges. +# :type nuclear_charges: numpy array +# :param mbtypes: Many-body types for the whole dataset, including 1-, 2- and 3-body types. Could be obtained by calling ``get_slatm_mbtypes()``. +# :type mbtypes: list +# :param local: Generate a local representation. Defaulted to False (i.e., global representation); otherwise, atomic version. +# :type local: bool +# :param sigmas: Controlling the width of Gaussian smearing function for 2- and 3-body parts, defaulted to [0.05,0.05], usually these do not need to be adjusted. +# :type sigmas: list +# :param dgrids: The interval between two sampled internuclear distances and angles, defaulted to [0.03,0.03], no need for change, compromised for speed and accuracy. +# :type dgrids: list +# :param rcut: Cut-off radius, defaulted to 4.8 Angstrom. +# :type rcut: float +# :param alchemy: Swith to use the alchemy version of SLATM. (default=False) +# :type alchemy: bool +# :param pbc: defaulted to '000', meaning it's a molecule; the three digits in the string corresponds to x,y,z direction +# :type pbc: string +# :param rpower: The power of R in 2-body potential, defaulted to London potential (=6). +# :type rpower: float +# :return: 1D SLATM representation +# :rtype: numpy array +# """ +# +# c = unit_cell +# iprt=False +# if c is None: +# c = np.array([[1,0,0],[0,1,0],[0,0,1]]) +# +# if pbc != '000': +# # print(' -- handling systems with periodic boundary condition') +# assert c != None, 'ERROR: Please specify unit cell for SLATM' +# # ======================================================================= +# # PBC may introduce new many-body terms, so at the stage of get statistics +# # info from db, we've already considered this point by letting maximal number +# # of nuclear charges being 3. +# # ======================================================================= +# +# zs = nuclear_charges +# na = len(zs) +# coords = coordinates +# obj = [ zs, coords, c ] +# +# iloc = local +# +# if iloc: +# mbs = [] +# X2Ns = [] +# for ia in range(na): +# # if iprt: print ' -- ia = ', ia + 1 +# n1 = 0; n2 = 0; n3 = 0 +# mbs_ia = np.zeros(0) +# icount = 0 +# for mbtype in mbtypes: +# if len(mbtype) == 1: +# mbsi = get_boa(mbtype[0], np.array([zs[ia],])) +# #print ' -- mbsi = ', mbsi +# if alchemy: +# n1 = 1 +# n1_0 = mbs_ia.shape[0] +# if n1_0 == 0: +# mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) +# elif n1_0 == 1: +# mbs_ia += mbsi +# else: +# raise '#ERROR' +# else: +# n1 += len(mbsi) +# mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) +# elif len(mbtype) == 2: +# #print ' 001, pbc = ', pbc +# mbsi = get_sbop(mbtype, obj, iloc=iloc, ia=ia, \ +# sigma=sigmas[0], dgrid=dgrids[0], rcut=rcut, \ +# pbc=pbc, rpower=rpower) +# mbsi *= 0.5 # only for the two-body parts, local rpst +# #print ' 002' +# if alchemy: +# n2 = len(mbsi) +# n2_0 = mbs_ia.shape[0] +# if n2_0 == n1: +# mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) +# elif n2_0 == n1 + n2: +# t = mbs_ia[n1:n1+n2] + mbsi +# mbs_ia[n1:n1+n2] = t +# else: +# raise '#ERROR' +# else: +# n2 += len(mbsi) +# mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) +# else: # len(mbtype) == 3: +# mbsi = get_sbot(mbtype, obj, iloc=iloc, ia=ia, \ +# sigma=sigmas[1], dgrid=dgrids[1], rcut=rcut, pbc=pbc) +# +# if alchemy: +# n3 = len(mbsi) +# n3_0 = mbs_ia.shape[0] +# if n3_0 == n1 + n2: +# mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) +# elif n3_0 == n1 + n2 + n3: +# t = mbs_ia[n1+n2:n1+n2+n3] + mbsi +# mbs_ia[n1+n2:n1+n2+n3] = t +# else: +# raise '#ERROR' +# else: +# n3 += len(mbsi) +# mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) +# +# mbs.append( mbs_ia ) +# X2N = [n1,n2,n3]; +# if X2N not in X2Ns: +# X2Ns.append(X2N) +# assert len(X2Ns) == 1, '#ERROR: multiple `X2N ???' +# else: +# n1 = 0; n2 = 0; n3 = 0 +# mbs = np.zeros(0) +# for mbtype in mbtypes: +# if len(mbtype) == 1: +# mbsi = get_boa(mbtype[0], zs) +# if alchemy: +# n1 = 1 +# n1_0 = mbs.shape[0] +# if n1_0 == 0: +# mbs = np.concatenate( (mbs, [sum(mbsi)] ), axis=0 ) +# elif n1_0 == 1: +# mbs += sum(mbsi ) +# else: +# raise '#ERROR' +# else: +# n1 += len(mbsi) +# mbs = np.concatenate( (mbs, mbsi), axis=0 ) +# elif len(mbtype) == 2: +# mbsi = get_sbop(mbtype, obj, sigma=sigmas[0], \ +# dgrid=dgrids[0], rcut=rcut, rpower=rpower) +# +# +# if alchemy: +# n2 = len(mbsi) +# n2_0 = mbs.shape[0] +# if n2_0 == n1: +# mbs = np.concatenate( (mbs, mbsi), axis=0 ) +# elif n2_0 == n1 + n2: +# t = mbs[n1:n1+n2] + mbsi +# mbs[n1:n1+n2] = t +# else: +# raise '#ERROR' +# else: +# n2 += len(mbsi) +# mbs = np.concatenate( (mbs, mbsi), axis=0 ) +# else: # len(mbtype) == 3: +# mbsi = get_sbot(mbtype, obj, sigma=sigmas[1], \ +# dgrid=dgrids[1], rcut=rcut) +# +# if alchemy: +# n3 = len(mbsi) +# n3_0 = mbs.shape[0] +# if n3_0 == n1 + n2: +# mbs = np.concatenate( (mbs, mbsi), axis=0 ) +# elif n3_0 == n1 + n2 + n3: +# t = mbs[n1+n2:n1+n2+n3] + mbsi +# mbs[n1+n2:n1+n2+n3] = t +# else: +# raise '#ERROR' +# else: +# n3 += len(mbsi) +# mbs = np.concatenate( (mbs, mbsi), axis=0 ) +# +# return mbs +# +#def generate_acsf(nuclear_charges, coordinates, elements = [1,6,7,8,16], nRs2 = 3, nRs3 = 3, nTs = 3, eta2 = 1, +# eta3 = 1, zeta = 1, rcut = 5, acut = 5, gradients = False): +# """ +# Generate the variant of atom-centered symmetry functions used in https://doi.org/10.1039/C7SC04934J +# +# :param nuclear_charges: List of nuclear charges. +# :type nuclear_charges: numpy array +# :param coordinates: Input coordinates +# :type coordinates: numpy array +# :param elements: list of unique nuclear charges (atom types) +# :type elements: numpy array +# :param nRs2: Number of gaussian basis functions in the two-body terms +# :type nRs2: integer +# :param nRs3: Number of gaussian basis functions in the three-body radial part +# :type nRs3: integer +# :param nTs: Number of basis functions in the three-body angular part +# :type nTs: integer +# :param eta2: Precision in the gaussian basis functions in the two-body terms +# :type eta2: float +# :param eta3: Precision in the gaussian basis functions in the three-body radial part +# :type eta3: float +# :param zeta: Precision parameter of basis functions in the three-body angular part +# :type zeta: float +# :param rcut: Cut-off radius of the two-body terms +# :type rcut: float +# :param acut: Cut-off radius of the three-body terms +# :type acut: float +# :param gradients: To return gradients or not +# :type gradients: boolean +# :return: Atom-centered symmetry functions representation +# :rtype: numpy array +# """ +# +# Rs2 = np.linspace(0, rcut, nRs2) +# Rs3 = np.linspace(0, acut, nRs3) +# Ts = np.linspace(0, np.pi, nTs) +# n_elements = len(elements) +# natoms = len(coordinates) +# +# descr_size = n_elements * nRs2 + (n_elements * (n_elements + 1)) // 2 * nRs3*nTs +# +# if gradients: +# return fgenerate_acsf_and_gradients(coordinates, nuclear_charges, elements, Rs2, Rs3, +# Ts, eta2, eta3, zeta, rcut, acut, natoms, descr_size) +# else: +# return fgenerate_acsf(coordinates, nuclear_charges, elements, Rs2, Rs3, +# Ts, eta2, eta3, zeta, rcut, acut, natoms, descr_size) class BaseRepresentation(object): @@ -592,7 +592,6 @@ def transform(self, **params): def generate(self, **params): raise NotImplementedError - class MolecularRepresentation(BaseRepresentation): """ Base class for molecular representations @@ -649,25 +648,54 @@ class CoulombMatrix(MolecularRepresentation): short_name = "cm" - def __init__(self, size=23, sorting="row-norm"): + def __init__(self, size=23, sorting="row-norm", data=False): """ :param size: The size of the largest molecule supported by the representation :type size: integer :param sorting: How the atom indices are sorted ('row-norm', 'unsorted') :type sorting: string + :param data: Optional Data object containing all molecules used in training \ + and/or prediction + :type data: Data object """ - def generate(nuclear_charges, coordinates): - :param nuclear_charges: Nuclear charges of the atoms in the molecule - :type nuclear_charges: numpy array - :param coordinates: 3D Coordinates of the atoms in the molecule - :type coordinates: numpy array + self.size = size + self.sorting = sorting + self.data = data + - :return: 1D representation - shape (size(size+1)/2,) + # TODO Make it possible to pass data in other ways as well + def generate(self, X, y=None): + """ + :param X: Data object or indices to use from the \ + Data object passed at initialization + :type X: Data object or array of indices + :param y: Dummy argument for scikit-learn + :type y: NoneType + :return: representations of shape (N, size(size+1)/2) :rtype: numpy array + """ + + if isinstance(X, Data): + # TODO Replace with isnone + if isinstance(X.natoms, type(None)): + print("Error: Empty Data object passed to CoulombMatrix representation") + raise SystemExit + if self.size > max(X.natoms): + print("Warning: Maximum size of system increased from %d to %d" + % (self.size, max(X.natoms))) + self.size = max(X.natoms) + + nuclear_charges = X.nuclear_charges + coordinates = X.coordinates + energies = X.energies + elif data: + + else: + error + # Raise warning is size too small and increase it - # Data options: Filenames, Data object if (sorting == "row-norm"): return fgenerate_coulomb_matrix(nuclear_charges, \ @@ -683,9 +711,56 @@ def generate(nuclear_charges, coordinates): class Data(object): """ - Temporary data class which should be replaced at some point + Temporary data class which should be replaced at some point by the ASE-interface """ - def __init__(self): - pass + def __init__(self, property_type = "energy", filenames = None): + """ + :param property_type: What kind of property will be predicted ('energy') + :type property_type: string + """ + + self.property_type = property_type + + self.ncompounds = 0 + self.coordinates = None + self.nuclear_charges = None + self.energies = None + self.properties = None + self.natoms = None + + if isinstance(filenames, list): + self.parse_xyz_files(filenames) + + def parse_xyz_files(self, filenames): + """ + Parse a list of xyz files. + """ + + self.ncompounds = len(filenames) + self.coordinates = np.empty(self.ncompounds, dtype=object) + self.nuclear_charges = np.empty(self.ncompounds, dtype=object) + self.natoms = np.empty(self.ncompounds, dtype = int) + + if self.property_type == "energy": + self.energies = np.empty(self.ncompounds, dtype=float) + + for i, filename in enumerate(filenames:) + with open(filename, "r") as f: + lines = f.readlines() + + self.natoms[i] = int(lines[0]) + self.nuclear_charges[i] = np.empty(natoms, dtype=int) + self.coordinates[i] = np.empty((natoms, 3), dtype=float) + + if self.property_type == "energy": + self.energies[i] = float(lines[1]) + + for j, line in enumerate(lines[2:self.natoms+2]): + tokens = line.split() + + if len(tokens) < 4: + break + nuclear_charges[i,j] = NUCLEAR_CHARGE[tokens[0]] + coordinates[i,j] = np.asarray(tokens[1:4], dtype=float) From ae9c6a9561e0f60e60f50ed7424549e86ff604d5 Mon Sep 17 00:00:00 2001 From: larsbratholm Date: Tue, 14 Aug 2018 18:32:30 +0100 Subject: [PATCH 04/25] Working basic example --- qml/__init__.py | 1 + qml/aglaia/aglaia.py | 6 +- qml/data/__init__.py | 2 +- qml/data/compound.py | 2 +- qml/data/data.py | 70 + qml/ml/arad/arad.py | 2 +- qml/ml/kernels/fkernels.f90 | 49 +- qml/ml/kernels/kernels.py | 105 +- qml/ml/kernels/wrappers.py | 2 +- qml/ml/representations/fchl.py | 2 +- qml/ml/representations/representations.py | 1257 +++++++++--------- qml/models/__init__.py | 3 +- qml/models/{mlmodel.py => basemodel.py} | 31 +- qml/models/kernelridge.py | 184 +-- qml/utils/__init__.py | 23 + qml/{ml/representations => utils}/alchemy.py | 0 qml/{aglaia => utils}/utils.py | 0 setup.py | 1 + test/test_armp.py | 2 +- test/test_mrmp.py | 2 +- test/test_neural_network.py | 2 +- test/test_qmlearn.py | 76 ++ 22 files changed, 1075 insertions(+), 747 deletions(-) create mode 100644 qml/data/data.py rename qml/models/{mlmodel.py => basemodel.py} (64%) create mode 100644 qml/utils/__init__.py rename qml/{ml/representations => utils}/alchemy.py (100%) rename qml/{aglaia => utils}/utils.py (100%) create mode 100644 test/test_qmlearn.py diff --git a/qml/__init__.py b/qml/__init__.py index 6e902c35d..9fb20b8f7 100644 --- a/qml/__init__.py +++ b/qml/__init__.py @@ -36,6 +36,7 @@ from . import ml from . import models from . import aglaia +from . import utils __author__ = "Anders S. Christensen" __copyright__ = "Copyright 2016" diff --git a/qml/aglaia/aglaia.py b/qml/aglaia/aglaia.py index e3b114ad7..73fb06991 100644 --- a/qml/aglaia/aglaia.py +++ b/qml/aglaia/aglaia.py @@ -30,13 +30,13 @@ import tensorflow as tf from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error from sklearn.base import BaseEstimator -from qml.aglaia.symm_funct import generate_parkhill_acsf -from qml.aglaia.utils import InputError, ceil, is_positive_or_zero, is_positive_integer, is_positive, \ +from .symm_funct import generate_parkhill_acsf +from ..utils import InputError, ceil, is_positive_or_zero, is_positive_integer, is_positive, \ is_bool, is_positive_integer_or_zero, is_string, is_positive_integer_array, is_array_like, is_none, \ check_global_representation, check_y, check_sizes, check_dy, check_classes, is_numeric_array, is_non_zero_integer, \ is_positive_integer_or_zero_array, check_local_representation -from qml.aglaia.tf_utils import TensorBoardLogger +from .tf_utils import TensorBoardLogger try: from qml.data import Compound diff --git a/qml/data/__init__.py b/qml/data/__init__.py index fb70a3ae2..ceca518b6 100644 --- a/qml/data/__init__.py +++ b/qml/data/__init__.py @@ -20,5 +20,5 @@ # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE # SOFTWARE. -from .xyzdataprovider import XYZDataProvider +from .data import Data from .compound import Compound diff --git a/qml/data/compound.py b/qml/data/compound.py index fad0937b1..c1d9cad01 100644 --- a/qml/data/compound.py +++ b/qml/data/compound.py @@ -25,7 +25,7 @@ import numpy as np import collections -from ..ml.representations.alchemy import NUCLEAR_CHARGE +from ..utils.alchemy import NUCLEAR_CHARGE from ..ml.representations import generate_coulomb_matrix from ..ml.representations import generate_atomic_coulomb_matrix diff --git a/qml/data/data.py b/qml/data/data.py new file mode 100644 index 000000000..9175ae048 --- /dev/null +++ b/qml/data/data.py @@ -0,0 +1,70 @@ + +from __future__ import print_function + +import glob +import numpy as np +from ..utils.alchemy import NUCLEAR_CHARGE + +class Data(object): + """ + Temporary data class which should be replaced at some point by the ASE-interface. + This could in principle also be replaced by a dictionary + + """ + + def __init__(self, filenames=None):#, property_type = "energy"): + #""" + #:param property_type: What kind of property will be predicted ('energy') + #:type property_type: string + #""" + + #self.property_type = property_type + + self.ncompounds = 0 + self.coordinates = None + self.nuclear_charges = None + #self.energies = None + #self.properties = None + self.natoms = None + + if isinstance(filenames, str): + filenames = sorted(glob.glob(filenames)) + if isinstance(filenames, list): + self._parse_xyz_files(filenames) + + def set_energies(self, energies): + self.energies = energies + + def _parse_xyz_files(self, filenames): + """ + Parse a list of xyz files. + """ + + self.ncompounds = len(filenames) + self.coordinates = np.empty(self.ncompounds, dtype=object) + self.nuclear_charges = np.empty(self.ncompounds, dtype=object) + self.natoms = np.empty(self.ncompounds, dtype = int) + + #if self.property_type == "energy": + # self.energies = np.empty(self.ncompounds, dtype=float) + + for i, filename in enumerate(filenames): + with open(filename, "r") as f: + lines = f.readlines() + + natoms = int(lines[0]) + self.natoms[i] = natoms + self.nuclear_charges[i] = np.empty(natoms, dtype=int) + self.coordinates[i] = np.empty((natoms, 3), dtype=float) + + #if self.property_type == "energy": + # self.energies[i] = float(lines[1]) + + for j, line in enumerate(lines[2:natoms+2]): + tokens = line.split() + + if len(tokens) < 4: + break + + self.nuclear_charges[i][j] = NUCLEAR_CHARGE[tokens[0]] + self.coordinates[i][j] = np.asarray(tokens[1:4], dtype=float) diff --git a/qml/ml/arad/arad.py b/qml/ml/arad/arad.py index 64ebac7d2..df36f9520 100644 --- a/qml/ml/arad/arad.py +++ b/qml/ml/arad/arad.py @@ -33,7 +33,7 @@ from .farad_kernels import fget_atomic_kernels_arad from .farad_kernels import fget_atomic_symmetric_kernels_arad -from ..representations.alchemy import PTP +from qml.utils.alchemy import PTP def getAngle(sp,norms): epsilon = 10.* np.finfo(float).eps diff --git a/qml/ml/kernels/fkernels.f90 b/qml/ml/kernels/fkernels.f90 index cdc4e0b4f..7cca75426 100644 --- a/qml/ml/kernels/fkernels.f90 +++ b/qml/ml/kernels/fkernels.f90 @@ -358,19 +358,58 @@ subroutine fgaussian_kernel(a, na, b, nb, k, sigma) allocate(temp(size(a, dim=1))) -!$OMP PARALLEL DO PRIVATE(temp) + !$OMP PARALLEL DO PRIVATE(temp) COLLAPSE(2) do i = 1, nb do j = 1, na temp(:) = a(:,j) - b(:,i) - k(j,i) = exp(inv_sigma * sum(temp*temp)) + k(j,i) = exp(inv_sigma * dot_product(temp,temp)) enddo enddo -!$OMP END PARALLEL DO + !$OMP END PARALLEL DO deallocate(temp) end subroutine fgaussian_kernel +subroutine fgaussian_kernel_symmetric(x, n, k, sigma) + + implicit none + + double precision, dimension(:,:), intent(in) :: x + + integer, intent(in) :: n + + double precision, dimension(:,:), intent(inout) :: k + double precision, intent(in) :: sigma + + double precision, allocatable, dimension(:) :: temp + double precision :: val + + double precision :: inv_sigma + integer :: i, j + + inv_sigma = -0.5d0 / (sigma*sigma) + + k = 1.0d0 + + allocate(temp(n)) + + !$OMP PARALLEL DO PRIVATE(temp, val) SCHEDULE(dynamic) + do i = 1, n + do j = i+1, n + temp = x(:,j) - x(:,i) + val = exp(inv_sigma * dot_product(temp,temp)) + k(j,i) = val + k(i,j) = val + enddo + enddo + !$OMP END PARALLEL DO + + deallocate(temp) + + +end subroutine fgaussian_kernel_symmetric + subroutine flaplacian_kernel(a, na, b, nb, k, sigma) implicit none @@ -389,13 +428,13 @@ subroutine flaplacian_kernel(a, na, b, nb, k, sigma) inv_sigma = -1.0d0 / sigma -!$OMP PARALLEL DO + !$OMP PARALLEL DO do i = 1, nb do j = 1, na k(j,i) = exp(inv_sigma * sum(abs(a(:,j) - b(:,i)))) enddo enddo -!$OMP END PARALLEL DO + !$OMP END PARALLEL DO end subroutine flaplacian_kernel diff --git a/qml/ml/kernels/kernels.py b/qml/ml/kernels/kernels.py index b92db417e..533197a1a 100644 --- a/qml/ml/kernels/kernels.py +++ b/qml/ml/kernels/kernels.py @@ -24,7 +24,10 @@ import numpy as np -from .fkernels import fgaussian_kernel +from qml.utils import is_none +from qml.data import Data + +from .fkernels import fgaussian_kernel, fgaussian_kernel_symmetric from .fkernels import flaplacian_kernel from .fkernels import flinear_kernel from .fkernels import fsargan_kernel @@ -250,7 +253,7 @@ def get_local_kernels_gaussian(A, B, na, nb, sigmas): nma = len(na) nmb = len(nb) - + sigmas = np.asarray(sigmas) nsigmas = len(sigmas) @@ -296,3 +299,101 @@ def get_local_kernels_laplacian(A, B, na, nb, sigmas): nsigmas = len(sigmas) return fget_local_kernels_laplacian(A.T, B.T, na, nb, sigmas, nma, nmb, nsigmas) + + +class BaseKernel(object): + """ + Base class for kernels + """ + + def fit(self, X): + """ + The fit routine is needed for scikit-learn pipelines. + This is actually never used but has to be here. + """ + raise NotImplementedError + + def transform(self, X, y = None): + """ + The transform routine is needed for scikit-learn pipelines. + """ + raise NotImplementedError + + def generate(self, **params): + raise NotImplementedError + + def _check_data_object(self, X): + # Check that we have a data object + if not isinstance(X, Data): + print("Error: Expected Data object as input in %s" % self.__class__.__name__) + raise SystemExit + + +class GaussianKernel(BaseKernel): + """ + Gaussian kernel + """ + + def __init__(self, sigma=1): + self.sigma = sigma + + self.representations = None + + def _set_representations(self, rep): + self.representations = rep + + def transform(self, X): + + self._check_data_object(X) + + # Kernel between representation stored in fit and representation + # given in data object. The order matters + kernel = self.generate(X.representations, self.representations) + + X.kernel = kernel + + return X + + def fit_transform(self, X, y=None): + + self._check_data_object(X) + + # Store representation for future transform calls + self._set_representations(X.representations) + + kernel = self.generate(X.representations) + + X.kernel = kernel + + return X + + + # TODO atomic + def generate(self, X, Y=None): + """ + Create a gaussian kernel from representations `X`. Optionally + an asymmetric kernel can be constructed between representations + `X` and `Y`. + + :param X: representations of shape (n_samplesX, representationX_size) + :type X: array + :param X: representations of shape (n_samplesY, representationY_size) + + :return: Gaussian kernel matrix of shape (n_samplesX, n_samplesX) if \ + Y=None else (n_samplesX, n_samplesY) + :rtype: array + """ + + # Note: Transposed for Fortran + n = X.shape[0] + + if is_none(Y): + # Do symmetric matrix + K = np.empty((n, n), order='F') + fgaussian_kernel_symmetric(X.T, n, K, self.sigma) + else: + m = Y.shape[0] + K = np.empty((n, m), order='F') + fgaussian_kernel(X.T, n, Y.T, m, K, self.sigma) + + return K diff --git a/qml/ml/kernels/wrappers.py b/qml/ml/kernels/wrappers.py index d129e6d0b..2bade4bd1 100644 --- a/qml/ml/kernels/wrappers.py +++ b/qml/ml/kernels/wrappers.py @@ -1,6 +1,6 @@ # MIT License # -# Copyright (c) 2016-2017 Anders Steen Christensen, Felix A. Faber, Lars A. Bratholm +# Copyright (c) 2016-2016 Anders Steen Christensen, Felix A. Faber, Lars A. Bratholm # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal diff --git a/qml/ml/representations/fchl.py b/qml/ml/representations/fchl.py index 7c1c3f18c..a8cb8d8a4 100644 --- a/qml/ml/representations/fchl.py +++ b/qml/ml/representations/fchl.py @@ -30,7 +30,7 @@ from .ffchl_module import fget_atomic_kernels_fchl from .ffchl_module import fget_atomic_symmetric_kernels_fchl -from .alchemy import get_alchemy +from qml.utils.alchemy import get_alchemy def generate_representation(coordinates, nuclear_charges, max_size=23, neighbors=23, cut_distance = 5.0, cell=None): diff --git a/qml/ml/representations/representations.py b/qml/ml/representations/representations.py index f24d9fd31..8223a3749 100644 --- a/qml/ml/representations/representations.py +++ b/qml/ml/representations/representations.py @@ -22,8 +22,15 @@ from __future__ import print_function + import numpy as np import itertools as itl +import glob +import copy + +from qml.data import Data + +from qml.utils import is_none, is_positive_integer_or_zero_array from .frepresentations import fgenerate_coulomb_matrix from .frepresentations import fgenerate_unsorted_coulomb_matrix @@ -32,7 +39,7 @@ from .frepresentations import fgenerate_eigenvalue_coulomb_matrix from .frepresentations import fgenerate_bob -from .alchemy import NUCLEAR_CHARGE +from qml.utils.alchemy import NUCLEAR_CHARGE from .slatm import get_boa from .slatm import get_sbop @@ -40,534 +47,534 @@ from .facsf import fgenerate_acsf, fgenerate_acsf_and_gradients -#def generate_coulomb_matrix(nuclear_charges, coordinates, size = 23, sorting = "row-norm"): -# """ Creates a Coulomb Matrix representation of a molecule. -# Sorting of the elements can either be done by ``sorting="row-norm"`` or ``sorting="unsorted"``. -# A matrix :math:`M` is constructed with elements -# -# .. math:: -# -# M_{ij} = -# \\begin{cases} -# \\tfrac{1}{2} Z_{i}^{2.4} & \\text{if } i = j \\\\ -# \\frac{Z_{i}Z_{j}}{\\| {\\bf R}_{i} - {\\bf R}_{j}\\|} & \\text{if } i \\neq j -# \\end{cases}, -# -# where :math:`i` and :math:`j` are atom indices, :math:`Z` is nuclear charge and -# :math:`\\bf R` is the coordinate in euclidean space. -# If ``sorting = 'row-norm'``, the atom indices are reordered such that -# -# :math:`\\sum_j M_{1j}^2 \\geq \\sum_j M_{2j}^2 \\geq ... \\geq \\sum_j M_{nj}^2` -# -# The upper triangular of M, including the diagonal, is concatenated to a 1D -# vector representation. -# -# If ``sorting = 'unsorted``, the elements are sorted in the same order as the input coordinates -# and nuclear charges. -# -# The representation is calculated using an OpenMP parallel Fortran routine. -# -# :param nuclear_charges: Nuclear charges of the atoms in the molecule -# :type nuclear_charges: numpy array -# :param coordinates: 3D Coordinates of the atoms in the molecule -# :type coordinates: numpy array -# :param size: The size of the largest molecule supported by the representation -# :type size: integer -# :param sorting: How the atom indices are sorted ('row-norm', 'unsorted') -# :type sorting: string -# -# :return: 1D representation - shape (size(size+1)/2,) -# :rtype: numpy array -# """ -# -# if (sorting == "row-norm"): -# return fgenerate_coulomb_matrix(nuclear_charges, \ -# coordinates, size) -# -# elif (sorting == "unsorted"): -# return fgenerate_unsorted_coulomb_matrix(nuclear_charges, \ -# coordinates, size) -# -# else: -# print("ERROR: Unknown sorting scheme requested") -# raise SystemExit -# -#def generate_atomic_coulomb_matrix(nuclear_charges, coordinates, size = 23, sorting = "distance", -# central_cutoff = 1e6, central_decay = -1, interaction_cutoff = 1e6, interaction_decay = -1, -# indices = None): -# """ Creates a Coulomb Matrix representation of the local environment of a central atom. -# For each central atom :math:`k`, a matrix :math:`M` is constructed with elements -# -# .. math:: -# -# M_{ij}(k) = -# \\begin{cases} -# \\tfrac{1}{2} Z_{i}^{2.4} \\cdot f_{ik}^2 & \\text{if } i = j \\\\ -# \\frac{Z_{i}Z_{j}}{\\| {\\bf R}_{i} - {\\bf R}_{j}\\|} \\cdot f_{ik}f_{jk}f_{ij} & \\text{if } i \\neq j -# \\end{cases}, -# -# where :math:`i`, :math:`j` and :math:`k` are atom indices, :math:`Z` is nuclear charge and -# :math:`\\bf R` is the coordinate in euclidean space. -# -# :math:`f_{ij}` is a function that masks long range effects: -# -# .. math:: -# -# f_{ij} = -# \\begin{cases} -# 1 & \\text{if } \\|{\\bf R}_{i} - {\\bf R}_{j} \\| \\leq r - \Delta r \\\\ -# \\tfrac{1}{2} \\big(1 + \\cos\\big(\\pi \\tfrac{\\|{\\bf R}_{i} - {\\bf R}_{j} \\| -# - r + \Delta r}{\Delta r} \\big)\\big) -# & \\text{if } r - \Delta r < \\|{\\bf R}_{i} - {\\bf R}_{j} \\| \\leq r - \Delta r \\\\ -# 0 & \\text{if } \\|{\\bf R}_{i} - {\\bf R}_{j} \\| > r -# \\end{cases}, -# -# where the parameters ``central_cutoff`` and ``central_decay`` corresponds to the variables -# :math:`r` and :math:`\Delta r` respectively for interactions involving the central atom, -# and ``interaction_cutoff`` and ``interaction_decay`` corresponds to the variables -# :math:`r` and :math:`\Delta r` respectively for interactions not involving the central atom. -# -# if ``sorting = 'row-norm'``, the atom indices are ordered such that -# -# :math:`\\sum_j M_{1j}(k)^2 \\geq \\sum_j M_{2j}(k)^2 \\geq ... \\geq \\sum_j M_{nj}(k)^2` -# -# if ``sorting = 'distance'``, the atom indices are ordered such that -# -# .. math:: -# -# \\|{\\bf R}_{1} - {\\bf R}_{k}\\| \\leq \\|{\\bf R}_{2} - {\\bf R}_{k}\\| -# \\leq ... \\leq \\|{\\bf R}_{n} - {\\bf R}_{k}\\| -# -# The upper triangular of M, including the diagonal, is concatenated to a 1D -# vector representation. -# -# The representation can be calculated for a subset by either specifying -# ``indices = [0,1,...]``, where :math:`[0,1,...]` are the requested atom indices, -# or by specifying ``indices = 'C'`` to only calculate central carbon atoms. -# -# The representation is calculated using an OpenMP parallel Fortran routine. -# -# :param nuclear_charges: Nuclear charges of the atoms in the molecule -# :type nuclear_charges: numpy array -# :param coordinates: 3D Coordinates of the atoms in the molecule -# :type coordinates: numpy array -# :param size: The size of the largest molecule supported by the representation -# :type size: integer -# :param sorting: How the atom indices are sorted ('row-norm', 'distance') -# :type sorting: string -# :param central_cutoff: The distance from the central atom, where the coulomb interaction -# element will be zero -# :type central_cutoff: float -# :param central_decay: The distance over which the the coulomb interaction decays from full to none -# :type central_decay: float -# :param interaction_cutoff: The distance between two non-central atom, where the coulomb interaction -# element will be zero -# :type interaction_cutoff: float -# :param interaction_decay: The distance over which the the coulomb interaction decays from full to none -# :type interaction_decay: float -# :param indices: Subset indices or atomtype -# :type indices: Nonetype/array/string -# -# -# :return: nD representation - shape (:math:`N_{atoms}`, size(size+1)/2) -# :rtype: numpy array -# """ -# -# -# if indices == None: -# nindices = len(nuclear_charges) -# indices = np.arange(1,1+nindices, 1, dtype = int) -# elif type("") == type(indices): -# if indices in NUCLEAR_CHARGE: -# indices = np.where(nuclear_charges == NUCLEAR_CHARGE[indices])[0] + 1 -# nindices = indices.size -# if nindices == 0: -# return np.zeros((0,0)) -# -# else: -# print("ERROR: Unknown value %s given for 'indices' variable" % indices) -# raise SystemExit -# else: -# indices = np.asarray(indices, dtype = int) + 1 -# nindices = indices.size -# -# -# if (sorting == "row-norm"): -# return fgenerate_local_coulomb_matrix(indices, nindices, nuclear_charges, -# coordinates, nuclear_charges.size, size, -# central_cutoff, central_decay, interaction_cutoff, interaction_decay) -# -# elif (sorting == "distance"): -# return fgenerate_atomic_coulomb_matrix(indices, nindices, nuclear_charges, -# coordinates, nuclear_charges.size, size, -# central_cutoff, central_decay, interaction_cutoff, interaction_decay) -# -# else: -# print("ERROR: Unknown sorting scheme requested") -# raise SystemExit -# -#def generate_eigenvalue_coulomb_matrix(nuclear_charges, coordinates, size = 23): -# """ Creates an eigenvalue Coulomb Matrix representation of a molecule. -# A matrix :math:`M` is constructed with elements -# -# .. math:: -# -# M_{ij} = -# \\begin{cases} -# \\tfrac{1}{2} Z_{i}^{2.4} & \\text{if } i = j \\\\ -# \\frac{Z_{i}Z_{j}}{\\| {\\bf R}_{i} - {\\bf R}_{j}\\|} & \\text{if } i \\neq j -# \\end{cases}, -# -# where :math:`i` and :math:`j` are atom indices, :math:`Z` is nuclear charge and -# :math:`\\bf R` is the coordinate in euclidean space. -# The molecular representation of the molecule is then the sorted eigenvalues of M. -# The representation is calculated using an OpenMP parallel Fortran routine. -# -# :param nuclear_charges: Nuclear charges of the atoms in the molecule -# :type nuclear_charges: numpy array -# :param coordinates: 3D Coordinates of the atoms in the molecule -# :type coordinates: numpy array -# :param size: The size of the largest molecule supported by the representation -# :type size: integer -# -# :return: 1D representation - shape (size, ) -# :rtype: numpy array -# """ -# return fgenerate_eigenvalue_coulomb_matrix(nuclear_charges, -# coordinates, size) -# -#def generate_bob(nuclear_charges, coordinates, atomtypes, size=23, asize = {"O":3, "C":7, "N":3, "H":16, "S":1}): -# """ Creates a Bag of Bonds (BOB) representation of a molecule. -# The representation expands on the coulomb matrix representation. -# For each element a bag (vector) is constructed for self interactions -# (e.g. ('C', 'H', 'O')). -# For each element pair a bag is constructed for interatomic interactions -# (e.g. ('CC', 'CH', 'CO', 'HH', 'HO', 'OO')), sorted by value. -# The self interaction of element :math:`I` is given by -# -# :math:`\\tfrac{1}{2} Z_{I}^{2.4}`, -# -# with :math:`Z_{i}` being the nuclear charge of element :math:`i` -# The interaction between atom :math:`i` of element :math:`I` and -# atom :math:`j` of element :math:`J` is given by -# -# :math:`\\frac{Z_{I}Z_{J}}{\\| {\\bf R}_{i} - {\\bf R}_{j}\\|}` -# -# with :math:`R_{i}` being the euclidean coordinate of atom :math:`i`. -# The sorted bags are concatenated to an 1D vector representation. -# The representation is calculated using an OpenMP parallel Fortran routine. -# -# :param nuclear_charges: Nuclear charges of the atoms in the molecule -# :type nuclear_charges: numpy array -# :param coordinates: 3D Coordinates of the atoms in the molecule -# :type coordinates: numpy array -# :param size: The maximum number of atoms in the representation -# :type size: integer -# :param asize: The maximum number of atoms of each element type supported by the representation -# :type asize: dictionary -# -# :return: 1D representation -# :rtype: numpy array -# """ -# -# n = 0 -# atoms = sorted(asize, key=asize.get) -# nmax = [asize[key] for key in atoms] -# ids = np.zeros(len(nmax), dtype=int) -# for i, (key, value) in enumerate(zip(atoms,nmax)): -# n += value * (1+value) -# ids[i] = NUCLEAR_CHARGE[key] -# for j in range(i): -# v = nmax[j] -# n += 2 * value * v -# n /= 2 -# -# return fgenerate_bob(nuclear_charges, coordinates, nuclear_charges, ids, nmax, n) -# -#def get_slatm_mbtypes(nuclear_charges, pbc='000'): -# """ -# Get the list of minimal types of many-body terms in a dataset. This resulting list -# is necessary as input in the ``generate_slatm_representation()`` function. -# -# :param nuclear_charges: A list of the nuclear charges for each compound in the dataset. -# :type nuclear_charges: list of numpy arrays -# :param pbc: periodic boundary condition along x,y,z direction, defaulted to '000', i.e., molecule -# :type pbc: string -# :return: A list containing the types of many-body terms. -# :rtype: list -# """ -# -# zs = nuclear_charges -# -# nm = len(zs) -# zsmax = set() -# nas = [] -# zs_ravel = [] -# for zsi in zs: -# na = len(zsi); nas.append(na) -# zsil = list(zsi); zs_ravel += zsil -# zsmax.update( zsil ) -# -# zsmax = np.array( list(zsmax) ) -# nass = [] -# for i in range(nm): -# zsi = np.array(zs[i],np.int) -# nass.append( [ (zi == zsi).sum() for zi in zsmax ] ) -# -# nzmax = np.max(np.array(nass), axis=0) -# nzmax_u = [] -# if pbc != '000': -# # the PBC will introduce new many-body terms, so set -# # nzmax to 3 if it's less than 3 -# for nzi in nzmax: -# if nzi <= 2: -# nzi = 3 -# nzmax_u.append(nzi) -# nzmax = nzmax_u -# -# boas = [ [zi,] for zi in zsmax ] -# bops = [ [zi,zi] for zi in zsmax ] + list( itl.combinations(zsmax,2) ) -# -# bots = [] -# for i in zsmax: -# for bop in bops: -# j,k = bop -# tas = [ [i,j,k], [i,k,j], [j,i,k] ] -# for tasi in tas: -# if (tasi not in bots) and (tasi[::-1] not in bots): -# nzsi = [ (zj == tasi).sum() for zj in zsmax ] -# if np.all(nzsi <= nzmax): -# bots.append( tasi ) -# mbtypes = boas + bops + bots -# -# return mbtypes #, np.array(zs_ravel), np.array(nas) -# -#def generate_slatm(coordinates, nuclear_charges, mbtypes, -# unit_cell=None, local=False, sigmas=[0.05,0.05], dgrids=[0.03,0.03], -# rcut=4.8, alchemy=False, pbc='000', rpower=6): -# """ -# Generate Spectrum of London and Axillrod-Teller-Muto potential (SLATM) representation. -# Both global (``local=False``) and local (``local=True``) SLATM are available. -# -# A version that works for periodic boundary conditions will be released soon. -# -# NOTE: You will need to run the ``get_slatm_mbtypes()`` function to get the ``mbtypes`` input (or generate it manually). -# -# :param coordinates: Input coordinates -# :type coordinates: numpy array -# :param nuclear_charges: List of nuclear charges. -# :type nuclear_charges: numpy array -# :param mbtypes: Many-body types for the whole dataset, including 1-, 2- and 3-body types. Could be obtained by calling ``get_slatm_mbtypes()``. -# :type mbtypes: list -# :param local: Generate a local representation. Defaulted to False (i.e., global representation); otherwise, atomic version. -# :type local: bool -# :param sigmas: Controlling the width of Gaussian smearing function for 2- and 3-body parts, defaulted to [0.05,0.05], usually these do not need to be adjusted. -# :type sigmas: list -# :param dgrids: The interval between two sampled internuclear distances and angles, defaulted to [0.03,0.03], no need for change, compromised for speed and accuracy. -# :type dgrids: list -# :param rcut: Cut-off radius, defaulted to 4.8 Angstrom. -# :type rcut: float -# :param alchemy: Swith to use the alchemy version of SLATM. (default=False) -# :type alchemy: bool -# :param pbc: defaulted to '000', meaning it's a molecule; the three digits in the string corresponds to x,y,z direction -# :type pbc: string -# :param rpower: The power of R in 2-body potential, defaulted to London potential (=6). -# :type rpower: float -# :return: 1D SLATM representation -# :rtype: numpy array -# """ -# -# c = unit_cell -# iprt=False -# if c is None: -# c = np.array([[1,0,0],[0,1,0],[0,0,1]]) -# -# if pbc != '000': -# # print(' -- handling systems with periodic boundary condition') -# assert c != None, 'ERROR: Please specify unit cell for SLATM' -# # ======================================================================= -# # PBC may introduce new many-body terms, so at the stage of get statistics -# # info from db, we've already considered this point by letting maximal number -# # of nuclear charges being 3. -# # ======================================================================= -# -# zs = nuclear_charges -# na = len(zs) -# coords = coordinates -# obj = [ zs, coords, c ] -# -# iloc = local -# -# if iloc: -# mbs = [] -# X2Ns = [] -# for ia in range(na): -# # if iprt: print ' -- ia = ', ia + 1 -# n1 = 0; n2 = 0; n3 = 0 -# mbs_ia = np.zeros(0) -# icount = 0 -# for mbtype in mbtypes: -# if len(mbtype) == 1: -# mbsi = get_boa(mbtype[0], np.array([zs[ia],])) -# #print ' -- mbsi = ', mbsi -# if alchemy: -# n1 = 1 -# n1_0 = mbs_ia.shape[0] -# if n1_0 == 0: -# mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) -# elif n1_0 == 1: -# mbs_ia += mbsi -# else: -# raise '#ERROR' -# else: -# n1 += len(mbsi) -# mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) -# elif len(mbtype) == 2: -# #print ' 001, pbc = ', pbc -# mbsi = get_sbop(mbtype, obj, iloc=iloc, ia=ia, \ -# sigma=sigmas[0], dgrid=dgrids[0], rcut=rcut, \ -# pbc=pbc, rpower=rpower) -# mbsi *= 0.5 # only for the two-body parts, local rpst -# #print ' 002' -# if alchemy: -# n2 = len(mbsi) -# n2_0 = mbs_ia.shape[0] -# if n2_0 == n1: -# mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) -# elif n2_0 == n1 + n2: -# t = mbs_ia[n1:n1+n2] + mbsi -# mbs_ia[n1:n1+n2] = t -# else: -# raise '#ERROR' -# else: -# n2 += len(mbsi) -# mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) -# else: # len(mbtype) == 3: -# mbsi = get_sbot(mbtype, obj, iloc=iloc, ia=ia, \ -# sigma=sigmas[1], dgrid=dgrids[1], rcut=rcut, pbc=pbc) -# -# if alchemy: -# n3 = len(mbsi) -# n3_0 = mbs_ia.shape[0] -# if n3_0 == n1 + n2: -# mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) -# elif n3_0 == n1 + n2 + n3: -# t = mbs_ia[n1+n2:n1+n2+n3] + mbsi -# mbs_ia[n1+n2:n1+n2+n3] = t -# else: -# raise '#ERROR' -# else: -# n3 += len(mbsi) -# mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) -# -# mbs.append( mbs_ia ) -# X2N = [n1,n2,n3]; -# if X2N not in X2Ns: -# X2Ns.append(X2N) -# assert len(X2Ns) == 1, '#ERROR: multiple `X2N ???' -# else: -# n1 = 0; n2 = 0; n3 = 0 -# mbs = np.zeros(0) -# for mbtype in mbtypes: -# if len(mbtype) == 1: -# mbsi = get_boa(mbtype[0], zs) -# if alchemy: -# n1 = 1 -# n1_0 = mbs.shape[0] -# if n1_0 == 0: -# mbs = np.concatenate( (mbs, [sum(mbsi)] ), axis=0 ) -# elif n1_0 == 1: -# mbs += sum(mbsi ) -# else: -# raise '#ERROR' -# else: -# n1 += len(mbsi) -# mbs = np.concatenate( (mbs, mbsi), axis=0 ) -# elif len(mbtype) == 2: -# mbsi = get_sbop(mbtype, obj, sigma=sigmas[0], \ -# dgrid=dgrids[0], rcut=rcut, rpower=rpower) -# -# -# if alchemy: -# n2 = len(mbsi) -# n2_0 = mbs.shape[0] -# if n2_0 == n1: -# mbs = np.concatenate( (mbs, mbsi), axis=0 ) -# elif n2_0 == n1 + n2: -# t = mbs[n1:n1+n2] + mbsi -# mbs[n1:n1+n2] = t -# else: -# raise '#ERROR' -# else: -# n2 += len(mbsi) -# mbs = np.concatenate( (mbs, mbsi), axis=0 ) -# else: # len(mbtype) == 3: -# mbsi = get_sbot(mbtype, obj, sigma=sigmas[1], \ -# dgrid=dgrids[1], rcut=rcut) -# -# if alchemy: -# n3 = len(mbsi) -# n3_0 = mbs.shape[0] -# if n3_0 == n1 + n2: -# mbs = np.concatenate( (mbs, mbsi), axis=0 ) -# elif n3_0 == n1 + n2 + n3: -# t = mbs[n1+n2:n1+n2+n3] + mbsi -# mbs[n1+n2:n1+n2+n3] = t -# else: -# raise '#ERROR' -# else: -# n3 += len(mbsi) -# mbs = np.concatenate( (mbs, mbsi), axis=0 ) -# -# return mbs -# -#def generate_acsf(nuclear_charges, coordinates, elements = [1,6,7,8,16], nRs2 = 3, nRs3 = 3, nTs = 3, eta2 = 1, -# eta3 = 1, zeta = 1, rcut = 5, acut = 5, gradients = False): -# """ -# Generate the variant of atom-centered symmetry functions used in https://doi.org/10.1039/C7SC04934J -# -# :param nuclear_charges: List of nuclear charges. -# :type nuclear_charges: numpy array -# :param coordinates: Input coordinates -# :type coordinates: numpy array -# :param elements: list of unique nuclear charges (atom types) -# :type elements: numpy array -# :param nRs2: Number of gaussian basis functions in the two-body terms -# :type nRs2: integer -# :param nRs3: Number of gaussian basis functions in the three-body radial part -# :type nRs3: integer -# :param nTs: Number of basis functions in the three-body angular part -# :type nTs: integer -# :param eta2: Precision in the gaussian basis functions in the two-body terms -# :type eta2: float -# :param eta3: Precision in the gaussian basis functions in the three-body radial part -# :type eta3: float -# :param zeta: Precision parameter of basis functions in the three-body angular part -# :type zeta: float -# :param rcut: Cut-off radius of the two-body terms -# :type rcut: float -# :param acut: Cut-off radius of the three-body terms -# :type acut: float -# :param gradients: To return gradients or not -# :type gradients: boolean -# :return: Atom-centered symmetry functions representation -# :rtype: numpy array -# """ -# -# Rs2 = np.linspace(0, rcut, nRs2) -# Rs3 = np.linspace(0, acut, nRs3) -# Ts = np.linspace(0, np.pi, nTs) -# n_elements = len(elements) -# natoms = len(coordinates) -# -# descr_size = n_elements * nRs2 + (n_elements * (n_elements + 1)) // 2 * nRs3*nTs -# -# if gradients: -# return fgenerate_acsf_and_gradients(coordinates, nuclear_charges, elements, Rs2, Rs3, -# Ts, eta2, eta3, zeta, rcut, acut, natoms, descr_size) -# else: -# return fgenerate_acsf(coordinates, nuclear_charges, elements, Rs2, Rs3, -# Ts, eta2, eta3, zeta, rcut, acut, natoms, descr_size) +def generate_coulomb_matrix(nuclear_charges, coordinates, size = 23, sorting = "row-norm"): + """ Creates a Coulomb Matrix representation of a molecule. + Sorting of the elements can either be done by ``sorting="row-norm"`` or ``sorting="unsorted"``. + A matrix :math:`M` is constructed with elements + + .. math:: + + M_{ij} = + \\begin{cases} + \\tfrac{1}{2} Z_{i}^{2.4} & \\text{if } i = j \\\\ + \\frac{Z_{i}Z_{j}}{\\| {\\bf R}_{i} - {\\bf R}_{j}\\|} & \\text{if } i \\neq j + \\end{cases}, + + where :math:`i` and :math:`j` are atom indices, :math:`Z` is nuclear charge and + :math:`\\bf R` is the coordinate in euclidean space. + If ``sorting = 'row-norm'``, the atom indices are reordered such that + + :math:`\\sum_j M_{1j}^2 \\geq \\sum_j M_{2j}^2 \\geq ... \\geq \\sum_j M_{nj}^2` + + The upper triangular of M, including the diagonal, is concatenated to a 1D + vector representation. + + If ``sorting = 'unsorted``, the elements are sorted in the same order as the input coordinates + and nuclear charges. + + The representation is calculated using an OpenMP parallel Fortran routine. + + :param nuclear_charges: Nuclear charges of the atoms in the molecule + :type nuclear_charges: numpy array + :param coordinates: 3D Coordinates of the atoms in the molecule + :type coordinates: numpy array + :param size: The size of the largest molecule supported by the representation + :type size: integer + :param sorting: How the atom indices are sorted ('row-norm', 'unsorted') + :type sorting: string + + :return: 1D representation - shape (size(size+1)/2,) + :rtype: numpy array + """ + + if (sorting == "row-norm"): + return fgenerate_coulomb_matrix(nuclear_charges, \ + coordinates, size) + + elif (sorting == "unsorted"): + return fgenerate_unsorted_coulomb_matrix(nuclear_charges, \ + coordinates, size) + + else: + print("ERROR: Unknown sorting scheme requested") + raise SystemExit + +def generate_atomic_coulomb_matrix(nuclear_charges, coordinates, size = 23, sorting = "distance", + central_cutoff = 1e6, central_decay = -1, interaction_cutoff = 1e6, interaction_decay = -1, + indices = None): + """ Creates a Coulomb Matrix representation of the local environment of a central atom. + For each central atom :math:`k`, a matrix :math:`M` is constructed with elements + + .. math:: + + M_{ij}(k) = + \\begin{cases} + \\tfrac{1}{2} Z_{i}^{2.4} \\cdot f_{ik}^2 & \\text{if } i = j \\\\ + \\frac{Z_{i}Z_{j}}{\\| {\\bf R}_{i} - {\\bf R}_{j}\\|} \\cdot f_{ik}f_{jk}f_{ij} & \\text{if } i \\neq j + \\end{cases}, + + where :math:`i`, :math:`j` and :math:`k` are atom indices, :math:`Z` is nuclear charge and + :math:`\\bf R` is the coordinate in euclidean space. + + :math:`f_{ij}` is a function that masks long range effects: + + .. math:: + + f_{ij} = + \\begin{cases} + 1 & \\text{if } \\|{\\bf R}_{i} - {\\bf R}_{j} \\| \\leq r - \Delta r \\\\ + \\tfrac{1}{2} \\big(1 + \\cos\\big(\\pi \\tfrac{\\|{\\bf R}_{i} - {\\bf R}_{j} \\| + - r + \Delta r}{\Delta r} \\big)\\big) + & \\text{if } r - \Delta r < \\|{\\bf R}_{i} - {\\bf R}_{j} \\| \\leq r - \Delta r \\\\ + 0 & \\text{if } \\|{\\bf R}_{i} - {\\bf R}_{j} \\| > r + \\end{cases}, + + where the parameters ``central_cutoff`` and ``central_decay`` corresponds to the variables + :math:`r` and :math:`\Delta r` respectively for interactions involving the central atom, + and ``interaction_cutoff`` and ``interaction_decay`` corresponds to the variables + :math:`r` and :math:`\Delta r` respectively for interactions not involving the central atom. + + if ``sorting = 'row-norm'``, the atom indices are ordered such that + + :math:`\\sum_j M_{1j}(k)^2 \\geq \\sum_j M_{2j}(k)^2 \\geq ... \\geq \\sum_j M_{nj}(k)^2` + + if ``sorting = 'distance'``, the atom indices are ordered such that + + .. math:: + + \\|{\\bf R}_{1} - {\\bf R}_{k}\\| \\leq \\|{\\bf R}_{2} - {\\bf R}_{k}\\| + \\leq ... \\leq \\|{\\bf R}_{n} - {\\bf R}_{k}\\| + + The upper triangular of M, including the diagonal, is concatenated to a 1D + vector representation. + + The representation can be calculated for a subset by either specifying + ``indices = [0,1,...]``, where :math:`[0,1,...]` are the requested atom indices, + or by specifying ``indices = 'C'`` to only calculate central carbon atoms. + + The representation is calculated using an OpenMP parallel Fortran routine. + + :param nuclear_charges: Nuclear charges of the atoms in the molecule + :type nuclear_charges: numpy array + :param coordinates: 3D Coordinates of the atoms in the molecule + :type coordinates: numpy array + :param size: The size of the largest molecule supported by the representation + :type size: integer + :param sorting: How the atom indices are sorted ('row-norm', 'distance') + :type sorting: string + :param central_cutoff: The distance from the central atom, where the coulomb interaction + element will be zero + :type central_cutoff: float + :param central_decay: The distance over which the the coulomb interaction decays from full to none + :type central_decay: float + :param interaction_cutoff: The distance between two non-central atom, where the coulomb interaction + element will be zero + :type interaction_cutoff: float + :param interaction_decay: The distance over which the the coulomb interaction decays from full to none + :type interaction_decay: float + :param indices: Subset indices or atomtype + :type indices: Nonetype/array/string + + + :return: nD representation - shape (:math:`N_{atoms}`, size(size+1)/2) + :rtype: numpy array + """ + + + if indices == None: + nindices = len(nuclear_charges) + indices = np.arange(1,1+nindices, 1, dtype = int) + elif type("") == type(indices): + if indices in NUCLEAR_CHARGE: + indices = np.where(nuclear_charges == NUCLEAR_CHARGE[indices])[0] + 1 + nindices = indices.size + if nindices == 0: + return np.zeros((0,0)) + + else: + print("ERROR: Unknown value %s given for 'indices' variable" % indices) + raise SystemExit + else: + indices = np.asarray(indices, dtype = int) + 1 + nindices = indices.size + + + if (sorting == "row-norm"): + return fgenerate_local_coulomb_matrix(indices, nindices, nuclear_charges, + coordinates, nuclear_charges.size, size, + central_cutoff, central_decay, interaction_cutoff, interaction_decay) + + elif (sorting == "distance"): + return fgenerate_atomic_coulomb_matrix(indices, nindices, nuclear_charges, + coordinates, nuclear_charges.size, size, + central_cutoff, central_decay, interaction_cutoff, interaction_decay) + + else: + print("ERROR: Unknown sorting scheme requested") + raise SystemExit + +def generate_eigenvalue_coulomb_matrix(nuclear_charges, coordinates, size = 23): + """ Creates an eigenvalue Coulomb Matrix representation of a molecule. + A matrix :math:`M` is constructed with elements + + .. math:: + + M_{ij} = + \\begin{cases} + \\tfrac{1}{2} Z_{i}^{2.4} & \\text{if } i = j \\\\ + \\frac{Z_{i}Z_{j}}{\\| {\\bf R}_{i} - {\\bf R}_{j}\\|} & \\text{if } i \\neq j + \\end{cases}, + + where :math:`i` and :math:`j` are atom indices, :math:`Z` is nuclear charge and + :math:`\\bf R` is the coordinate in euclidean space. + The molecular representation of the molecule is then the sorted eigenvalues of M. + The representation is calculated using an OpenMP parallel Fortran routine. + + :param nuclear_charges: Nuclear charges of the atoms in the molecule + :type nuclear_charges: numpy array + :param coordinates: 3D Coordinates of the atoms in the molecule + :type coordinates: numpy array + :param size: The size of the largest molecule supported by the representation + :type size: integer + + :return: 1D representation - shape (size, ) + :rtype: numpy array + """ + return fgenerate_eigenvalue_coulomb_matrix(nuclear_charges, + coordinates, size) + +def generate_bob(nuclear_charges, coordinates, atomtypes, size=23, asize = {"O":3, "C":7, "N":3, "H":16, "S":1}): + """ Creates a Bag of Bonds (BOB) representation of a molecule. + The representation expands on the coulomb matrix representation. + For each element a bag (vector) is constructed for self interactions + (e.g. ('C', 'H', 'O')). + For each element pair a bag is constructed for interatomic interactions + (e.g. ('CC', 'CH', 'CO', 'HH', 'HO', 'OO')), sorted by value. + The self interaction of element :math:`I` is given by + + :math:`\\tfrac{1}{2} Z_{I}^{2.4}`, + + with :math:`Z_{i}` being the nuclear charge of element :math:`i` + The interaction between atom :math:`i` of element :math:`I` and + atom :math:`j` of element :math:`J` is given by + + :math:`\\frac{Z_{I}Z_{J}}{\\| {\\bf R}_{i} - {\\bf R}_{j}\\|}` + + with :math:`R_{i}` being the euclidean coordinate of atom :math:`i`. + The sorted bags are concatenated to an 1D vector representation. + The representation is calculated using an OpenMP parallel Fortran routine. + + :param nuclear_charges: Nuclear charges of the atoms in the molecule + :type nuclear_charges: numpy array + :param coordinates: 3D Coordinates of the atoms in the molecule + :type coordinates: numpy array + :param size: The maximum number of atoms in the representation + :type size: integer + :param asize: The maximum number of atoms of each element type supported by the representation + :type asize: dictionary + + :return: 1D representation + :rtype: numpy array + """ + + n = 0 + atoms = sorted(asize, key=asize.get) + nmax = [asize[key] for key in atoms] + ids = np.zeros(len(nmax), dtype=int) + for i, (key, value) in enumerate(zip(atoms,nmax)): + n += value * (1+value) + ids[i] = NUCLEAR_CHARGE[key] + for j in range(i): + v = nmax[j] + n += 2 * value * v + n /= 2 + + return fgenerate_bob(nuclear_charges, coordinates, nuclear_charges, ids, nmax, n) + +def get_slatm_mbtypes(nuclear_charges, pbc='000'): + """ + Get the list of minimal types of many-body terms in a dataset. This resulting list + is necessary as input in the ``generate_slatm_representation()`` function. + + :param nuclear_charges: A list of the nuclear charges for each compound in the dataset. + :type nuclear_charges: list of numpy arrays + :param pbc: periodic boundary condition along x,y,z direction, defaulted to '000', i.e., molecule + :type pbc: string + :return: A list containing the types of many-body terms. + :rtype: list + """ + + zs = nuclear_charges + + nm = len(zs) + zsmax = set() + nas = [] + zs_ravel = [] + for zsi in zs: + na = len(zsi); nas.append(na) + zsil = list(zsi); zs_ravel += zsil + zsmax.update( zsil ) + + zsmax = np.array( list(zsmax) ) + nass = [] + for i in range(nm): + zsi = np.array(zs[i],np.int) + nass.append( [ (zi == zsi).sum() for zi in zsmax ] ) + + nzmax = np.max(np.array(nass), axis=0) + nzmax_u = [] + if pbc != '000': + # the PBC will introduce new many-body terms, so set + # nzmax to 3 if it's less than 3 + for nzi in nzmax: + if nzi <= 2: + nzi = 3 + nzmax_u.append(nzi) + nzmax = nzmax_u + + boas = [ [zi,] for zi in zsmax ] + bops = [ [zi,zi] for zi in zsmax ] + list( itl.combinations(zsmax,2) ) + + bots = [] + for i in zsmax: + for bop in bops: + j,k = bop + tas = [ [i,j,k], [i,k,j], [j,i,k] ] + for tasi in tas: + if (tasi not in bots) and (tasi[::-1] not in bots): + nzsi = [ (zj == tasi).sum() for zj in zsmax ] + if np.all(nzsi <= nzmax): + bots.append( tasi ) + mbtypes = boas + bops + bots + + return mbtypes #, np.array(zs_ravel), np.array(nas) + +def generate_slatm(coordinates, nuclear_charges, mbtypes, + unit_cell=None, local=False, sigmas=[0.05,0.05], dgrids=[0.03,0.03], + rcut=4.8, alchemy=False, pbc='000', rpower=6): + """ + Generate Spectrum of London and Axillrod-Teller-Muto potential (SLATM) representation. + Both global (``local=False``) and local (``local=True``) SLATM are available. + + A version that works for periodic boundary conditions will be released soon. + + NOTE: You will need to run the ``get_slatm_mbtypes()`` function to get the ``mbtypes`` input (or generate it manually). + + :param coordinates: Input coordinates + :type coordinates: numpy array + :param nuclear_charges: List of nuclear charges. + :type nuclear_charges: numpy array + :param mbtypes: Many-body types for the whole dataset, including 1-, 2- and 3-body types. Could be obtained by calling ``get_slatm_mbtypes()``. + :type mbtypes: list + :param local: Generate a local representation. Defaulted to False (i.e., global representation); otherwise, atomic version. + :type local: bool + :param sigmas: Controlling the width of Gaussian smearing function for 2- and 3-body parts, defaulted to [0.05,0.05], usually these do not need to be adjusted. + :type sigmas: list + :param dgrids: The interval between two sampled internuclear distances and angles, defaulted to [0.03,0.03], no need for change, compromised for speed and accuracy. + :type dgrids: list + :param rcut: Cut-off radius, defaulted to 4.8 Angstrom. + :type rcut: float + :param alchemy: Swith to use the alchemy version of SLATM. (default=False) + :type alchemy: bool + :param pbc: defaulted to '000', meaning it's a molecule; the three digits in the string corresponds to x,y,z direction + :type pbc: string + :param rpower: The power of R in 2-body potential, defaulted to London potential (=6). + :type rpower: float + :return: 1D SLATM representation + :rtype: numpy array + """ + + c = unit_cell + iprt=False + if c is None: + c = np.array([[1,0,0],[0,1,0],[0,0,1]]) + + if pbc != '000': + # print(' -- handling systems with periodic boundary condition') + assert c != None, 'ERROR: Please specify unit cell for SLATM' + # ======================================================================= + # PBC may introduce new many-body terms, so at the stage of get statistics + # info from db, we've already considered this point by letting maximal number + # of nuclear charges being 3. + # ======================================================================= + + zs = nuclear_charges + na = len(zs) + coords = coordinates + obj = [ zs, coords, c ] + + iloc = local + + if iloc: + mbs = [] + X2Ns = [] + for ia in range(na): + # if iprt: print ' -- ia = ', ia + 1 + n1 = 0; n2 = 0; n3 = 0 + mbs_ia = np.zeros(0) + icount = 0 + for mbtype in mbtypes: + if len(mbtype) == 1: + mbsi = get_boa(mbtype[0], np.array([zs[ia],])) + #print ' -- mbsi = ', mbsi + if alchemy: + n1 = 1 + n1_0 = mbs_ia.shape[0] + if n1_0 == 0: + mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) + elif n1_0 == 1: + mbs_ia += mbsi + else: + raise '#ERROR' + else: + n1 += len(mbsi) + mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) + elif len(mbtype) == 2: + #print ' 001, pbc = ', pbc + mbsi = get_sbop(mbtype, obj, iloc=iloc, ia=ia, \ + sigma=sigmas[0], dgrid=dgrids[0], rcut=rcut, \ + pbc=pbc, rpower=rpower) + mbsi *= 0.5 # only for the two-body parts, local rpst + #print ' 002' + if alchemy: + n2 = len(mbsi) + n2_0 = mbs_ia.shape[0] + if n2_0 == n1: + mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) + elif n2_0 == n1 + n2: + t = mbs_ia[n1:n1+n2] + mbsi + mbs_ia[n1:n1+n2] = t + else: + raise '#ERROR' + else: + n2 += len(mbsi) + mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) + else: # len(mbtype) == 3: + mbsi = get_sbot(mbtype, obj, iloc=iloc, ia=ia, \ + sigma=sigmas[1], dgrid=dgrids[1], rcut=rcut, pbc=pbc) + + if alchemy: + n3 = len(mbsi) + n3_0 = mbs_ia.shape[0] + if n3_0 == n1 + n2: + mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) + elif n3_0 == n1 + n2 + n3: + t = mbs_ia[n1+n2:n1+n2+n3] + mbsi + mbs_ia[n1+n2:n1+n2+n3] = t + else: + raise '#ERROR' + else: + n3 += len(mbsi) + mbs_ia = np.concatenate( (mbs_ia, mbsi), axis=0 ) + + mbs.append( mbs_ia ) + X2N = [n1,n2,n3]; + if X2N not in X2Ns: + X2Ns.append(X2N) + assert len(X2Ns) == 1, '#ERROR: multiple `X2N ???' + else: + n1 = 0; n2 = 0; n3 = 0 + mbs = np.zeros(0) + for mbtype in mbtypes: + if len(mbtype) == 1: + mbsi = get_boa(mbtype[0], zs) + if alchemy: + n1 = 1 + n1_0 = mbs.shape[0] + if n1_0 == 0: + mbs = np.concatenate( (mbs, [sum(mbsi)] ), axis=0 ) + elif n1_0 == 1: + mbs += sum(mbsi ) + else: + raise '#ERROR' + else: + n1 += len(mbsi) + mbs = np.concatenate( (mbs, mbsi), axis=0 ) + elif len(mbtype) == 2: + mbsi = get_sbop(mbtype, obj, sigma=sigmas[0], \ + dgrid=dgrids[0], rcut=rcut, rpower=rpower) + + + if alchemy: + n2 = len(mbsi) + n2_0 = mbs.shape[0] + if n2_0 == n1: + mbs = np.concatenate( (mbs, mbsi), axis=0 ) + elif n2_0 == n1 + n2: + t = mbs[n1:n1+n2] + mbsi + mbs[n1:n1+n2] = t + else: + raise '#ERROR' + else: + n2 += len(mbsi) + mbs = np.concatenate( (mbs, mbsi), axis=0 ) + else: # len(mbtype) == 3: + mbsi = get_sbot(mbtype, obj, sigma=sigmas[1], \ + dgrid=dgrids[1], rcut=rcut) + + if alchemy: + n3 = len(mbsi) + n3_0 = mbs.shape[0] + if n3_0 == n1 + n2: + mbs = np.concatenate( (mbs, mbsi), axis=0 ) + elif n3_0 == n1 + n2 + n3: + t = mbs[n1+n2:n1+n2+n3] + mbsi + mbs[n1+n2:n1+n2+n3] = t + else: + raise '#ERROR' + else: + n3 += len(mbsi) + mbs = np.concatenate( (mbs, mbsi), axis=0 ) + + return mbs + +def generate_acsf(nuclear_charges, coordinates, elements = [1,6,7,8,16], nRs2 = 3, nRs3 = 3, nTs = 3, eta2 = 1, + eta3 = 1, zeta = 1, rcut = 5, acut = 5, gradients = False): + """ + Generate the variant of atom-centered symmetry functions used in https://doi.org/10.1039/C7SC04934J + + :param nuclear_charges: List of nuclear charges. + :type nuclear_charges: numpy array + :param coordinates: Input coordinates + :type coordinates: numpy array + :param elements: list of unique nuclear charges (atom types) + :type elements: numpy array + :param nRs2: Number of gaussian basis functions in the two-body terms + :type nRs2: integer + :param nRs3: Number of gaussian basis functions in the three-body radial part + :type nRs3: integer + :param nTs: Number of basis functions in the three-body angular part + :type nTs: integer + :param eta2: Precision in the gaussian basis functions in the two-body terms + :type eta2: float + :param eta3: Precision in the gaussian basis functions in the three-body radial part + :type eta3: float + :param zeta: Precision parameter of basis functions in the three-body angular part + :type zeta: float + :param rcut: Cut-off radius of the two-body terms + :type rcut: float + :param acut: Cut-off radius of the three-body terms + :type acut: float + :param gradients: To return gradients or not + :type gradients: boolean + :return: Atom-centered symmetry functions representation + :rtype: numpy array + """ + + Rs2 = np.linspace(0, rcut, nRs2) + Rs3 = np.linspace(0, acut, nRs3) + Ts = np.linspace(0, np.pi, nTs) + n_elements = len(elements) + natoms = len(coordinates) + + descr_size = n_elements * nRs2 + (n_elements * (n_elements + 1)) // 2 * nRs3*nTs + + if gradients: + return fgenerate_acsf_and_gradients(coordinates, nuclear_charges, elements, Rs2, Rs3, + Ts, eta2, eta3, zeta, rcut, acut, natoms, descr_size) + else: + return fgenerate_acsf(coordinates, nuclear_charges, elements, Rs2, Rs3, + Ts, eta2, eta3, zeta, rcut, acut, natoms, descr_size) class BaseRepresentation(object): @@ -575,47 +582,67 @@ class BaseRepresentation(object): Base class for representations """ - def fit(self, X, y = None): + def fit(self, X, y=None): """ - The fit routine is needed for scikit-learn pipelines + The fit routine is needed for scikit-learn pipelines. """ - return self + print("rep, fit") + raise NotImplementedError - def transform(self, **params): + def transform(self, X): """ The transform routine is needed for scikit-learn pipelines. - Just returns the generate routine, since generate is a more - natural name. """ - return self.generate(**params) + raise NotImplementedError - def generate(self, **params): + def generate(self, X, y=None): raise NotImplementedError + def _preprocess_input(self, X): + """ + Convenience function that processes X in a way such that + X can both be a data object or an array of indices. + + :param X: Data object or array of indices + :type X: Data object or array + :return: array of indices + :rtype: array + """ + + if isinstance(X, Data): + self._set_data(X) + self.data.indices = np.arange(self.data.natoms.size) + elif self.data and is_positive_integer_or_zero_array(X) \ + and X.max() <= self.data.natoms.size: + self.data.indices = X + else: + print("Expected X to be array of indices or Data object. Got %s" % str(X)) + raise SystemExit + + def _set_data(self, data): + if data and is_none(data.natoms): + print("Error: Empty Data object passed to %s representation" % self.__class__.__name__) + raise SystemExit + # Shallow copy should be fine + self.data = copy.copy(data) + + def _set_representation_type(self, X): + X.representation_type = self._representation_type + + class MolecularRepresentation(BaseRepresentation): """ Base class for molecular representations """ - representation_type = "molecular" + _representation_type = "molecular" class AtomicRepresentations(BaseRepresentation): """ Base class for molecular representations """ - representation_type = "atomic" - -#class AtomCenteredSymmetryFunctions(AtomicRepresentation): -# """ -# Variant of atom-centered symmetry functions used in https://doi.org/10.1039/C7SC04934J -# """ -# -# def __init__(self): -# pass -# -# def generate(self): -# pass + _representation_type = "atomic" class CoulombMatrix(MolecularRepresentation): """ Coulomb Matrix representation of a molecule. @@ -646,7 +673,7 @@ class CoulombMatrix(MolecularRepresentation): """ - short_name = "cm" + _representation_short_name = "cm" def __init__(self, size=23, sorting="row-norm", data=False): """ @@ -661,106 +688,80 @@ def __init__(self, size=23, sorting="row-norm", data=False): self.size = size self.sorting = sorting - self.data = data + self._set_data(data) - # TODO Make it possible to pass data in other ways as well - def generate(self, X, y=None): + def fit(self, X, y=None): """ :param X: Data object or indices to use from the \ Data object passed at initialization :type X: Data object or array of indices :param y: Dummy argument for scikit-learn :type y: NoneType - :return: representations of shape (N, size(size+1)/2) - :rtype: numpy array + :return: self + :rtype: object """ - if isinstance(X, Data): - # TODO Replace with isnone - if isinstance(X.natoms, type(None)): - print("Error: Empty Data object passed to CoulombMatrix representation") - raise SystemExit - if self.size > max(X.natoms): - print("Warning: Maximum size of system increased from %d to %d" - % (self.size, max(X.natoms))) - self.size = max(X.natoms) - - nuclear_charges = X.nuclear_charges - coordinates = X.coordinates - energies = X.energies - elif data: - - else: - error - - - # Raise warning is size too small and increase it + self._preprocess_input(X) + self._set_representation_type(X) - if (sorting == "row-norm"): - return fgenerate_coulomb_matrix(nuclear_charges, \ - coordinates, size) + natoms = self.data.natoms[self.data.indices] - elif (sorting == "unsorted"): - return fgenerate_unsorted_coulomb_matrix(nuclear_charges, \ - coordinates, size) + if self.size < max(natoms): + print("Warning: Maximum size of system increased from %d to %d" + % (self.size, max(natoms))) + self.size = max(natoms) - else: - print("ERROR: Unknown sorting scheme requested") - raise SystemExit - -class Data(object): - """ - Temporary data class which should be replaced at some point by the ASE-interface - """ + return self - def __init__(self, property_type = "energy", filenames = None): + def transform(self, X): """ - :param property_type: What kind of property will be predicted ('energy') - :type property_type: string + :param X: Data object or indices to use from the \ + Data object passed at initialization + :type X: Data object or array of indices + :return: Data object + :rtype: Data object """ - self.property_type = property_type - - self.ncompounds = 0 - self.coordinates = None - self.nuclear_charges = None - self.energies = None - self.properties = None - self.natoms = None - - if isinstance(filenames, list): - self.parse_xyz_files(filenames) + self._preprocess_input(X) - def parse_xyz_files(self, filenames): - """ - Parse a list of xyz files. - """ + nuclear_charges = self.data.nuclear_charges[self.data.indices] + coordinates = self.data.coordinates[self.data.indices] + natoms = self.data.natoms[self.data.indices] - self.ncompounds = len(filenames) - self.coordinates = np.empty(self.ncompounds, dtype=object) - self.nuclear_charges = np.empty(self.ncompounds, dtype=object) - self.natoms = np.empty(self.ncompounds, dtype = int) + if max(natoms) > self.size: + print("Error: Representation can't be generated for molecule of size \ + %d when variable 'size' is %d" % (max(natoms), self.size)) + raise SystemExit - if self.property_type == "energy": - self.energies = np.empty(self.ncompounds, dtype=float) + representations = [] + if (self.sorting == "row-norm"): + for charge, xyz in zip(nuclear_charges, coordinates): + representations.append( + fgenerate_coulomb_matrix(charge, xyz, self.size)) - for i, filename in enumerate(filenames:) - with open(filename, "r") as f: - lines = f.readlines() + elif (self.sorting == "unsorted"): + for charge, xyz in zip(nuclear_charges, coordinates): + representations.append( + fgenerate_unsorted_coulomb_matrix(charge, xyz, self.size)) + else: + print("ERROR: Unknown sorting scheme requested") + raise SystemExit - self.natoms[i] = int(lines[0]) - self.nuclear_charges[i] = np.empty(natoms, dtype=int) - self.coordinates[i] = np.empty((natoms, 3), dtype=float) + self.data.representations = np.asarray(representations) - if self.property_type == "energy": - self.energies[i] = float(lines[1]) + return self.data - for j, line in enumerate(lines[2:self.natoms+2]): - tokens = line.split() + # TODO Make it possible to pass data in other ways as well + # e.g. dictionary + def generate(self, X): + """ + :param X: Data object or indices to use from the \ + Data object passed at initialization + :type X: Data object or array of indices + :return: Representations of shape (n_samples, representation_size) + :rtype: array + """ - if len(tokens) < 4: - break + return self.fit(X).transform(X).representations - nuclear_charges[i,j] = NUCLEAR_CHARGE[tokens[0]] - coordinates[i,j] = np.asarray(tokens[1:4], dtype=float) diff --git a/qml/models/__init__.py b/qml/models/__init__.py index e375525cb..1047b011d 100644 --- a/qml/models/__init__.py +++ b/qml/models/__init__.py @@ -20,4 +20,5 @@ # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE # SOFTWARE. -from .kernelridge import GenericKRR +from .basemodel import BaseModel +from .kernelridge import KernelRidgeRegression diff --git a/qml/models/mlmodel.py b/qml/models/basemodel.py similarity index 64% rename from qml/models/mlmodel.py rename to qml/models/basemodel.py index 2241a7a43..4a3069b5b 100644 --- a/qml/models/mlmodel.py +++ b/qml/models/basemodel.py @@ -1,6 +1,6 @@ # MIT License # -# Copyright (c) 2017 Anders S. Christensen, Kristof T. Schutt, Stefan Chmiela +# Copyright (c) 2018 Lars A. Bratholm # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal @@ -22,31 +22,12 @@ from __future__ import division, absolute_import, print_function -class MLModel(object): +class BaseModel(object): - def __init__(self, targets): - self._targes = targets + _estimator_type = "regressor" - def train(self, data): - return self._train(data) - - def predict(self, data): - return self._predict(data) - - def save(self, path): - return self._save(path) - - def restore(self, path): - return self._restore(path) - - def _train(self, data): - raise NotImplementedError - - def _predict(self, data): - raise NotImplementedError - - def _save(self, path): + def fit(self, X): raise NotImplementedError - def _restore(self, path): - raise NotImplementedError + def predict(self, X): + return NotImplementedError diff --git a/qml/models/kernelridge.py b/qml/models/kernelridge.py index 9097241dd..51ab6b2b3 100644 --- a/qml/models/kernelridge.py +++ b/qml/models/kernelridge.py @@ -1,6 +1,6 @@ # MIT License # -# Copyright (c) 2017 Anders S. Christensen +# Copyright (c) 2017-2018 Anders S. Christensen, Lars A. Bratholm # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal @@ -22,93 +22,127 @@ from __future__ import division, absolute_import, print_function -from time import time - import numpy as np -import ase - -from ..ml.representations import generate_coulomb_matrix from ..ml.kernels import gaussian_kernel from ..ml.math import cho_solve -from .mlmodel import MLModel - -class GenericKRR(MLModel): +from . import BaseModel +from ..data import Data +from ..utils import is_numeric_array, is_none - def __init__(self, targets, representation="coulomb-matrix", - kernel="gaussian", sigma=4000.0, llambda=1e-10): +#class KernelRidgeRegression(BaseModel): +# +# def __init__(self, targets, representation="coulomb-matrix", +# kernel="gaussian", sigma=4000.0, llambda=1e-10): +# +# super(GenericKRR,self).__init__(targets) +# +# self.representation = representation +# self.kernel = kernel +# self.sigma = sigma +# self.llambda = llambda +# +# def _train(self, data): +# +# start = time() +# X = [] +# +# (properties, compounds) = data.get_properties() +# +# for row in compounds.select(): +# # for key in row: +# # print('{0:22}: {1}'.format(key, row[key])) +# X.append(generate_coulomb_matrix(row["numbers"], +# row["positions"])) +# +# self.Y = np.array(properties) +# self.X = np.array(X) +# self.K = gaussian_kernel(self.X, self.X, self.sigma) +# self.K[np.diag_indices_from(self.K)] += self.llambda +# self.alpha = cho_solve(self.K, self.Y) +# +# Yss = np.dot(self.K, self.alpha) +# mae_training_error = np.mean(np.abs(Yss - self.Y)) +# rmse_training_error = np.sqrt(np.mean(np.square(Yss - self.Y))) +# +# result = { +# "MAETrainingError": mae_training_error, +# "RMSETrainingError": rmse_training_error, +# "ElapsedTime": time() - start, +# } +# +# return result +# +# def _predict(self, data): +# +# Xs = [] +# +# (properties, compounds) = data.get_properties() +# +# for row in compounds.select(): +# # for key in row: +# # print('{0:22}: {1}'.format(key, row[key])) +# Xs.append(generate_coulomb_matrix(row["numbers"], +# row["positions"])) +# +# Xs = np.array(Xs) +# Ks = gaussian_kernel(Xs, self.X, self.sigma) +# +# return np.dot(Ks, self.alpha) +# +# +# def _restore(self, path, load_kernel=False): +# self.Y = np.load(path + "/Y.npy") +# self.X = np.load(path + "/X.npy") +# self.alpha = np.load(path + "/alpha.npy") +# if load_kernel: +# self.K = np.load(path + "/K.npy") +# +# +# def _save(self, path, save_kernel=False): +# np.save(path + "/X.npy", self.X) +# np.save(path + "/Y.npy", self.Y) +# np.save(path + "/alpha.npy", self.alpha) +# +# if save_kernel: +# np.save(path + "/K.npy") - super(GenericKRR,self).__init__(targets) +class KernelRidgeRegression(BaseModel): - self.representation = representation - self.kernel = kernel - self.sigma = sigma + def __init__(self, llambda=1e-10): self.llambda = llambda - def _train(self, data): - - start = time() - X = [] - - (properties, compounds) = data.get_properties() - - for row in compounds.select(): - # for key in row: - # print('{0:22}: {1}'.format(key, row[key])) - X.append(generate_coulomb_matrix(row["numbers"], - row["positions"])) - - self.Y = np.array(properties) - self.X = np.array(X) - self.K = gaussian_kernel(self.X, self.X, self.sigma) - self.K[np.diag_indices_from(self.K)] += self.llambda - self.alpha = cho_solve(self.K, self.Y) - - Yss = np.dot(self.K, self.alpha) - mae_training_error = np.mean(np.abs(Yss - self.Y)) - rmse_training_error = np.sqrt(np.mean(np.square(Yss - self.Y))) - - result = { - "MAETrainingError": mae_training_error, - "RMSETrainingError": rmse_training_error, - "ElapsedTime": time() - start, - } - - return result - - def _predict(self, data): - - Xs = [] - - (properties, compounds) = data.get_properties() - - for row in compounds.select(): - # for key in row: - # print('{0:22}: {1}'.format(key, row[key])) - Xs.append(generate_coulomb_matrix(row["numbers"], - row["positions"])) - - Xs = np.array(Xs) - Ks = gaussian_kernel(Xs, self.X, self.sigma) - - return np.dot(Ks, self.alpha) + def fit(self, X, y=None): + if isinstance(X, Data): + try: + K, y = X.kernel, X.energies[X.indices] + except: + print("No kernel matrix and/or energies found in data object in module %s" % self.__class__.__name__) + raise SystemExit + elif is_numeric_array(X) and X.ndim == 2 and X.shape[0] == X.shape[1] and not is_none(y): + K = X + else: + print("Expected variable 'X' to be kernel matrix or Data object. Got %s" % str(X)) + raise SystemExit - def _restore(self, path, load_kernel=False): - self.Y = np.load(path + "/Y.npy") - self.X = np.load(path + "/X.npy") - self.alpha = np.load(path + "/alpha.npy") - if load_kernel: - self.K = np.load(path + "/K.npy") + K[np.diag_indices_from(K)] += self.llambda - def _save(self, path, save_kernel=False): - np.save(path + "/X.npy", self.X) - np.save(path + "/Y.npy", self.Y) - np.save(path + "/alpha.npy", self.alpha) + self.alpha = cho_solve(K, y) - if save_kernel: - np.save(path + "/K.npy") + def predict(self, X): - + if isinstance(X, Data): + try: + K = X.kernel + except: + print("No kernel matrix found in data object in module %s" % self.__class__.__name__) + raise SystemExit + elif is_numeric_array(X) and X.ndim == 2 and X.shape[1] == self.alpha.size: + K = X + else: + print("Expected variable 'X' to be kernel matrix or Data object. Got %s" % str(X)) + raise SystemExit + return np.dot(K, self.alpha) diff --git a/qml/utils/__init__.py b/qml/utils/__init__.py new file mode 100644 index 000000000..931e7807f --- /dev/null +++ b/qml/utils/__init__.py @@ -0,0 +1,23 @@ +# MIT License +# +# Copyright (c) 2017-2018 Silvia Amabilino, Lars Andersen Bratholm +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +from .utils import * diff --git a/qml/ml/representations/alchemy.py b/qml/utils/alchemy.py similarity index 100% rename from qml/ml/representations/alchemy.py rename to qml/utils/alchemy.py diff --git a/qml/aglaia/utils.py b/qml/utils/utils.py similarity index 100% rename from qml/aglaia/utils.py rename to qml/utils/utils.py diff --git a/setup.py b/setup.py index 2ead66e22..d7d0df7eb 100755 --- a/setup.py +++ b/setup.py @@ -162,6 +162,7 @@ def setup_qml(): 'qml.ml.math', 'qml.ml.representations', 'qml.models', + 'qml.utils' ], # metadata diff --git a/test/test_armp.py b/test/test_armp.py index b7bc3ddad..4c6fb10ab 100644 --- a/test/test_armp.py +++ b/test/test_armp.py @@ -27,7 +27,7 @@ import numpy as np from qml.aglaia.aglaia import ARMP -from qml.aglaia.utils import InputError +from qml.utils import InputError import glob import os diff --git a/test/test_mrmp.py b/test/test_mrmp.py index ffe32e585..1f98d7602 100644 --- a/test/test_mrmp.py +++ b/test/test_mrmp.py @@ -27,7 +27,7 @@ import numpy as np from qml.aglaia.aglaia import MRMP -from qml.aglaia.utils import InputError +from qml.utils import InputError import glob import os import shutil diff --git a/test/test_neural_network.py b/test/test_neural_network.py index 775335f12..21ec5318b 100644 --- a/test/test_neural_network.py +++ b/test/test_neural_network.py @@ -29,7 +29,7 @@ # TODO relative imports from qml.aglaia.aglaia import MRMP -from qml.aglaia.utils import InputError +from qml.utils import InputError # ------------ ** All functions to test the inputs to the classes ** --------------- diff --git a/test/test_qmlearn.py b/test/test_qmlearn.py new file mode 100644 index 000000000..4194be8c3 --- /dev/null +++ b/test/test_qmlearn.py @@ -0,0 +1,76 @@ +import glob +from qml.ml.representations import * +from qml.ml.kernels import * +from qml.models import * +import numpy as np +from sklearn.pipeline import make_pipeline +import copy + +class Test(object): + + def fit(self, X, y=None): + print(X) + print("test,fit") + + def predict(self, X): + print("test, predict") + +class Rep(object): + + def fit(self, X, y=None): + print("rep, fit") + return self + + def transform(self, X): + print("rep, trans") + + def fit_transform(self, X, y=None): + print("rep, fit_trans") + return 1, 2, 3 + +if __name__ == "__main__": + train_data = Data(filenames="qm7/0*.xyz") + train_energies = np.loadtxt('data/hof_qm7.txt', usecols=1)[:train_data.ncompounds] + train_data.set_energies(train_energies) + test_data = Data(filenames="qm7/1*.xyz") + test_energies = np.loadtxt('data/hof_qm7.txt', usecols=1)[:test_data.ncompounds] + test_data.set_energies(test_energies) + + # Generate the representations + rep = CoulombMatrix().generate(train_data) + print(rep.shape) + # Generate a symmetric kernel + kernel = GaussianKernel().generate(rep) + print(kernel.shape) + # Generate an asymmetric kernel + kernel = GaussianKernel().generate(rep, rep[:10]) + print(kernel.shape) + # Generate symmetric kernel from pipeline + model = make_pipeline(CoulombMatrix(size=23), GaussianKernel(sigma=30)) + train_kernel = model.fit_transform(train_data).kernel + print(train_kernel.shape) + # Generate asymmetric kernel from pipeline after fit + test_kernel = model.transform(test_data).kernel + print(test_kernel.shape) + # Fit and predict KRR from kernels + model = KernelRidgeRegression() + model.fit(train_kernel, train_energies) + predictions = model.predict(test_kernel) + print(predictions.shape) + + # Fit and predict KRR from pipeline + model = make_pipeline(CoulombMatrix(size=23), GaussianKernel(sigma=30), KernelRidgeRegression(llambda=1e-6)) + model.fit(train_data) + predictions = model.predict(test_data) + print(predictions.shape) + + + #model = make_pipeline(Rep(), Test()) + #print("fit") + #model.fit(data) + #print("predict") + #model.predict(data) + #quit() + #rep = model.fit_transform(data) + #print(rep.shape) + From 35b53325ff895c119c833f992e07ede23d984b50 Mon Sep 17 00:00:00 2001 From: larsbratholm Date: Wed, 15 Aug 2018 17:58:43 +0100 Subject: [PATCH 05/25] Mostly hacked the searchcv routines to work --- qml/data/data.py | 38 +++++++++++++ qml/ml/kernels/kernels.py | 4 +- qml/ml/representations/representations.py | 16 +++--- qml/models/basemodel.py | 25 ++++++++- test/test_qmlearn.py | 67 ++++++++++++++--------- 5 files changed, 114 insertions(+), 36 deletions(-) diff --git a/qml/data/data.py b/qml/data/data.py index 9175ae048..afcedf278 100644 --- a/qml/data/data.py +++ b/qml/data/data.py @@ -5,6 +5,7 @@ import numpy as np from ..utils.alchemy import NUCLEAR_CHARGE + class Data(object): """ Temporary data class which should be replaced at some point by the ASE-interface. @@ -32,6 +33,43 @@ def __init__(self, filenames=None):#, property_type = "energy"): if isinstance(filenames, list): self._parse_xyz_files(filenames) + # Hack for sklearn CV + self.shape = (self.ncompounds,) + # # Hack for sklearn CV + # self.iloc = self.iloc_override(self) + + #class iloc_override(object): + # """ + # Hack for sklearn CV. + # """ + + # def __init__(self, parent): + # self.parent = parent + + # def __getitem__(self, i): + # self.parent.indices = i + # return self.parent + + def take(self, i, axis=None): + self.indices = i + return self + + + + def __getitem__(self, i): + """ + Hack for sklearn CV + """ + return i + + def __len__(self): + """ + Hack for sklearn CV + """ + return self.ncompounds + + + def set_energies(self, energies): self.energies = energies diff --git a/qml/ml/kernels/kernels.py b/qml/ml/kernels/kernels.py index 533197a1a..90a081a6b 100644 --- a/qml/ml/kernels/kernels.py +++ b/qml/ml/kernels/kernels.py @@ -24,6 +24,8 @@ import numpy as np +from sklearn.base import BaseEstimator + from qml.utils import is_none from qml.data import Data @@ -301,7 +303,7 @@ def get_local_kernels_laplacian(A, B, na, nb, sigmas): return fget_local_kernels_laplacian(A.T, B.T, na, nb, sigmas, nma, nmb, nsigmas) -class BaseKernel(object): +class BaseKernel(BaseEstimator): """ Base class for kernels """ diff --git a/qml/ml/representations/representations.py b/qml/ml/representations/representations.py index 8223a3749..75fbd1a48 100644 --- a/qml/ml/representations/representations.py +++ b/qml/ml/representations/representations.py @@ -28,6 +28,8 @@ import glob import copy +from sklearn.base import BaseEstimator + from qml.data import Data from qml.utils import is_none, is_positive_integer_or_zero_array @@ -577,7 +579,7 @@ def generate_acsf(nuclear_charges, coordinates, elements = [1,6,7,8,16], nRs2 = Ts, eta2, eta3, zeta, rcut, acut, natoms, descr_size) -class BaseRepresentation(object): +class BaseRepresentation(BaseEstimator): """ Base class for representations """ @@ -586,7 +588,6 @@ def fit(self, X, y=None): """ The fit routine is needed for scikit-learn pipelines. """ - print("rep, fit") raise NotImplementedError def transform(self, X): @@ -613,7 +614,7 @@ def _preprocess_input(self, X): self._set_data(X) self.data.indices = np.arange(self.data.natoms.size) elif self.data and is_positive_integer_or_zero_array(X) \ - and X.max() <= self.data.natoms.size: + and max(X) <= self.data.natoms.size: self.data.indices = X else: print("Expected X to be array of indices or Data object. Got %s" % str(X)) @@ -624,11 +625,11 @@ def _set_data(self, data): print("Error: Empty Data object passed to %s representation" % self.__class__.__name__) raise SystemExit # Shallow copy should be fine + #TODO solve clone issue self.data = copy.copy(data) - def _set_representation_type(self, X): - X.representation_type = self._representation_type - + def _set_representation_type(self): + self.data.representation_type = self._representation_type class MolecularRepresentation(BaseRepresentation): """ @@ -690,7 +691,6 @@ def __init__(self, size=23, sorting="row-norm", data=False): self.sorting = sorting self._set_data(data) - def fit(self, X, y=None): """ :param X: Data object or indices to use from the \ @@ -703,7 +703,7 @@ def fit(self, X, y=None): """ self._preprocess_input(X) - self._set_representation_type(X) + self._set_representation_type() natoms = self.data.natoms[self.data.indices] diff --git a/qml/models/basemodel.py b/qml/models/basemodel.py index 4a3069b5b..c8a5c00fd 100644 --- a/qml/models/basemodel.py +++ b/qml/models/basemodel.py @@ -22,12 +22,35 @@ from __future__ import division, absolute_import, print_function -class BaseModel(object): +from sklearn.base import BaseEstimator + +class BaseModel(BaseEstimator): _estimator_type = "regressor" + def __init__(self, scoring='mae'): + self.scoring = 'mae' + def fit(self, X): raise NotImplementedError def predict(self, X): return NotImplementedError + + def score(self, X, y=None): + + return 1 + + y_pred = self.predict(X) + + if isinstance(X, Data): + try: + K, y = X.kernel, X.energies[X.indices] + except: + print("No kernel matrix and/or energies found in data object in module %s" % self.__class__.__name__) + raise SystemExit + elif is_numeric_array(X) and X.ndim == 2 and X.shape[0] == X.shape[1] and not is_none(y): + K = X + else: + print("Expected variable 'X' to be kernel matrix or Data object. Got %s" % str(X)) + raise SystemExit diff --git a/test/test_qmlearn.py b/test/test_qmlearn.py index 4194be8c3..a6317807f 100644 --- a/test/test_qmlearn.py +++ b/test/test_qmlearn.py @@ -4,6 +4,7 @@ from qml.models import * import numpy as np from sklearn.pipeline import make_pipeline +from sklearn.model_selection import GridSearchCV import copy class Test(object): @@ -29,40 +30,54 @@ def fit_transform(self, X, y=None): return 1, 2, 3 if __name__ == "__main__": - train_data = Data(filenames="qm7/0*.xyz") + train_data = Data(filenames="qm7/00*.xyz") train_energies = np.loadtxt('data/hof_qm7.txt', usecols=1)[:train_data.ncompounds] train_data.set_energies(train_energies) test_data = Data(filenames="qm7/1*.xyz") test_energies = np.loadtxt('data/hof_qm7.txt', usecols=1)[:test_data.ncompounds] test_data.set_energies(test_energies) - # Generate the representations - rep = CoulombMatrix().generate(train_data) - print(rep.shape) - # Generate a symmetric kernel - kernel = GaussianKernel().generate(rep) - print(kernel.shape) - # Generate an asymmetric kernel - kernel = GaussianKernel().generate(rep, rep[:10]) - print(kernel.shape) - # Generate symmetric kernel from pipeline - model = make_pipeline(CoulombMatrix(size=23), GaussianKernel(sigma=30)) - train_kernel = model.fit_transform(train_data).kernel - print(train_kernel.shape) - # Generate asymmetric kernel from pipeline after fit - test_kernel = model.transform(test_data).kernel - print(test_kernel.shape) - # Fit and predict KRR from kernels - model = KernelRidgeRegression() - model.fit(train_kernel, train_energies) - predictions = model.predict(test_kernel) - print(predictions.shape) + ## Generate the representations + #rep = CoulombMatrix().generate(train_data) + #print(rep.shape) + ## Generate a symmetric kernel + #kernel = GaussianKernel().generate(rep) + #print(kernel.shape) + ## Generate an asymmetric kernel + #kernel = GaussianKernel().generate(rep, rep[:10]) + #print(kernel.shape) + ## Generate symmetric kernel from pipeline + #model = make_pipeline(CoulombMatrix(size=23), GaussianKernel(sigma=30)) + #train_kernel = model.fit_transform(train_data).kernel + #print(train_kernel.shape) + ## Generate asymmetric kernel from pipeline after fit + #test_kernel = model.transform(test_data).kernel + #print(test_kernel.shape) + ## Fit and predict KRR from kernels + #model = KernelRidgeRegression() + #model.fit(train_kernel, train_energies) + #predictions = model.predict(test_kernel) + #print(predictions.shape) # Fit and predict KRR from pipeline - model = make_pipeline(CoulombMatrix(size=23), GaussianKernel(sigma=30), KernelRidgeRegression(llambda=1e-6)) - model.fit(train_data) - predictions = model.predict(test_data) - print(predictions.shape) + model = make_pipeline(CoulombMatrix(size=23, data=train_data), GaussianKernel(sigma=30), KernelRidgeRegression(llambda=1e-6)) + #model.fit(train_data) + #predictions = model.predict(test_data) + #print(predictions.shape) + + # Gridsearch CV of hyperparams + params = {'gaussiankernel__sigma': [1, 3, 10, 30, 100, 300], + 'kernelridgeregression__llambda': [1e-10,1e-8,1e-6,1e-4] + } + + grid = GridSearchCV(model, cv=3, refit=False, param_grid = params) + #grid.fit(train_data) + #print(grid.cv_results_) + #print(grid.best_params_) + grid.fit(range(20)) + print(grid.cv_results_) + print(grid.best_params_) + #model = make_pipeline(Rep(), Test()) From 3f39602f4677e1bb2f7a6d6bb33796e15ea8ea4a Mon Sep 17 00:00:00 2001 From: larsbratholm Date: Fri, 17 Aug 2018 17:25:10 +0100 Subject: [PATCH 06/25] Implementing atomic gaussian kernel --- examples/qm7_hyperparam_search/config.yaml | 35 +- examples/qm7_hyperparam_search/idx.csv | 1716 +++++++---------- examples/qm7_hyperparam_search/idx.npy | Bin 8080 -> 0 bytes .../make_model_pickle.py | 23 +- examples/qm7_hyperparam_search/model.pickle | Bin 1039473 -> 4444353 bytes .../qm7_hyperparam_search/osprey-trials.db | Bin 217088 -> 40960 bytes qml/data/__init__.py | 2 +- qml/data/compound.py | 684 +++---- qml/data/data.py | 90 +- qml/ml/kernels/fkernels.f90 | 73 +- qml/ml/kernels/kernels.py | 128 +- qml/ml/kernels/wrappers.py | 62 +- qml/ml/representations/representations.py | 305 ++- qml/models/__init__.py | 2 +- qml/models/basemodel.py | 62 +- qml/models/kernelridge.py | 45 +- test/test_qmlearn.py | 71 +- 17 files changed, 1695 insertions(+), 1603 deletions(-) delete mode 100644 examples/qm7_hyperparam_search/idx.npy diff --git a/examples/qm7_hyperparam_search/config.yaml b/examples/qm7_hyperparam_search/config.yaml index 79a0f6a08..f805eccd6 100644 --- a/examples/qm7_hyperparam_search/config.yaml +++ b/examples/qm7_hyperparam_search/config.yaml @@ -2,44 +2,23 @@ estimator: pickle: model.pickle strategy: - name: random + name: sobol #params: - # seeds: 5 + # seeds: 30 search_space: - hl1: + gaussiankernel__sigma: min: 1 - max: 100 - type: int - - hl2: - min: 0 - max: 100 - type: int - - hl3: - min: 0 - max: 100 - type: int - - l1_reg: - min: 1e-9 - max: 1e0 + max: 1000 type: float warp: log - l2_reg: - min: 1e-9 - max: 1e0 + kernelridgeregression__l2_reg: + min: 1e-10 + max: 1e-2 type: float warp: log - learning_rate: - min: 1e-9 - max: 1e0 - type: float - warp: log - cv: name: kfold params: diff --git a/examples/qm7_hyperparam_search/idx.csv b/examples/qm7_hyperparam_search/idx.csv index ac96ab9d3..9b039b475 100644 --- a/examples/qm7_hyperparam_search/idx.csv +++ b/examples/qm7_hyperparam_search/idx.csv @@ -1,1000 +1,750 @@ -0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -10 -11 -12 -13 -14 -15 -16 -17 -18 -19 -20 -21 -22 -23 -24 -25 -26 -27 -28 -29 +5707 +3988 +5940 +367 +2832 +2540 +5733 +7001 +3233 +6120 +4154 +2365 +6840 +4445 +5498 +7036 +7054 +5473 +4761 +1360 +4622 +2639 +1269 +435 +1423 +1439 +5539 +469 +5100 +5189 +7004 +4973 +2188 +6364 +282 +1532 +6775 +3168 +4337 +3063 +4193 +4364 +5369 +6409 +5644 +1378 +5982 +4299 +4482 +4272 +5831 +1595 +7075 +2058 +6937 +420 +5096 +386 +3282 +4816 +5736 +5136 +6165 +6473 +854 +5901 +1726 +2599 +2406 +2688 +5990 +2420 +1718 +4664 +2439 +4353 +4541 +6841 +6685 +6375 +2730 +1456 +6990 +6234 +5344 +3570 +5725 +3778 +3544 +584 30 -31 -32 -33 -34 -35 -36 -37 -38 -39 -40 -41 -42 -43 -44 -45 -46 -47 -48 -49 -50 -51 -52 +6567 +1312 +370 +5602 +717 +3173 +5205 +4106 +5134 +4505 +1199 +5491 +1596 +3737 +6299 +5529 +705 +1000 +3752 +6224 +6715 +3892 +364 +5082 +3305 +5192 +2354 +2222 +319 +4667 +5088 +4805 +6744 +610 +4099 +417 +2825 +5413 +5887 +1990 +1293 +2586 +5355 +5124 +6703 +1435 +5771 +5267 +5713 +4820 +1605 +4184 +4127 +1167 +979 +4249 +4681 +4264 +2549 +5572 +3643 +6427 +6711 +5805 +1145 +4426 +5537 +3738 +3873 +2481 +5008 +3457 +1482 +6341 +1665 +2636 +4394 +6163 +3667 +3211 +3965 +3454 +1625 +5404 +5098 +5419 +975 +1398 +2644 +6336 +600 +5196 +3184 +3401 +3670 +3159 +3444 +1957 +4076 +3238 +3320 +1920 +2145 +6086 +6870 +6895 +778 +895 +2883 +5800 +5937 +1182 +3219 +1716 +2180 +607 +5616 +5770 +1084 +2732 +6192 +6109 +1654 +564 +3040 +709 +873 +3068 +2002 +5858 +6756 +4226 +3256 +3781 +4657 +1031 +5969 +5297 +5567 +4969 +5279 +255 +4211 53 -54 -55 -56 -57 -58 -59 -60 -61 -62 -63 -64 -65 -66 -67 -68 -69 -70 -71 -72 -73 -74 -75 -76 -77 -78 -79 -80 -81 -82 -83 -84 -85 -86 -87 -88 -89 -90 -91 -92 -93 -94 -95 -96 -97 -98 -99 -100 -101 -102 -103 -104 -105 -106 -107 -108 -109 -110 -111 -112 -113 -114 -115 -116 -117 -118 -119 -120 -121 -122 -123 -124 -125 -126 -127 -128 -129 -130 +2392 +762 +5850 +5716 +4153 +4315 +5790 +3757 +6229 +3036 131 -132 -133 -134 -135 -136 -137 -138 -139 -140 -141 +2746 +6202 +2585 +614 +742 +6363 +4520 +310 +5701 +1071 +1832 +3276 +2372 +1228 +5180 +6147 +3813 +2778 +2625 +1467 +5878 142 -143 -144 -145 -146 -147 -148 -149 -150 -151 -152 -153 -154 -155 -156 -157 -158 -159 -160 -161 -162 -163 -164 -165 -166 -167 -168 -169 -170 -171 -172 -173 -174 -175 -176 -177 -178 +6868 +1797 +3200 +3973 +935 +2094 +3895 +4366 +1241 +5246 +2262 +1362 +446 +738 +1569 +325 +2125 +5507 +6286 +2964 +4634 +1723 +2770 +718 +3902 179 -180 -181 -182 -183 -184 -185 -186 -187 -188 -189 -190 -191 -192 -193 -194 -195 -196 -197 -198 -199 -200 -201 -202 -203 -204 -205 -206 -207 -208 -209 -210 -211 -212 -213 -214 -215 -216 -217 -218 -219 -220 -221 +4608 +6908 +3327 +88 +6727 +4167 +670 +52 +2062 +6492 +4934 +6098 +101 +6188 +4091 +4285 +5542 +5243 +1348 +3299 +3124 +3084 +6408 +2049 +2749 +6013 +6802 +5746 +4899 +6643 +568 +3608 +829 +5157 +4605 +6044 +2721 +3592 +1644 +6047 +3319 +720 +5663 +632 +3507 +6106 +6016 +728 +3916 +2736 +5118 +4972 +6770 +808 +1571 +3583 +5483 +1940 +1243 +1230 +4143 +3466 +280 +6232 +2134 +6809 +5208 +4429 +624 +1866 +2304 +4557 +305 +2056 +2573 +2651 +998 +5595 +2617 +6696 +817 +1888 +1613 +1115 +3222 +1794 +1282 +5762 +3721 +5226 +3120 +6845 +4983 +6677 +2498 +7087 +5385 +529 +4551 +425 +4117 +4173 +5545 +5499 +6656 +4784 +835 +5612 +3158 +2903 +3409 +6926 +2474 +3657 +2402 +2922 +5844 +1260 +759 +3227 +823 +4469 +1091 +3617 +1131 +4433 +6534 +4919 +6951 +734 +4384 +2454 +6973 +5977 +4411 +3712 +3028 +5647 +1366 +2880 +5842 +407 +6759 +2706 +3918 +6469 +1911 +98 +696 +4441 +5214 +2097 +6414 +5376 +4659 +4174 +17 +6787 +175 +4107 +5702 222 -223 -224 -225 -226 -227 -228 -229 -230 -231 -232 +5760 +4247 +2980 +3179 +924 +5345 +3591 +2799 +4447 +5256 +5904 +3472 +5560 +2528 +394 +1564 +6879 +4570 +2346 +2028 +3307 +3023 +2361 +3945 +2659 +3954 +5600 +5227 +855 +5563 +6815 +767 +6021 233 +1387 +2671 +1839 234 -235 -236 -237 -238 -239 -240 -241 -242 -243 -244 -245 -246 -247 -248 +5374 +6626 +4945 +3610 +1406 +2384 +6142 +2890 +1056 +475 +1487 +4195 +796 +6634 +5548 +6005 +412 +543 +3244 +2360 +5788 +1426 +1736 +1805 +6521 +4877 +5229 +1492 +5097 +3743 +2514 +6207 +6938 +3018 +1941 +1266 +5952 +4976 +3234 +124 +1052 +7077 +6805 +595 +2698 +6360 +1261 +1862 +2153 +2842 +6612 +641 +5206 +5005 +4845 +5528 +3041 +2115 +1849 +7092 +2319 +6091 +6407 +4690 +6801 +4163 249 -250 -251 -252 -253 -254 -255 -256 -257 -258 -259 -260 -261 -262 -263 -264 -265 -266 +1937 +6562 +994 +5919 +2 267 -268 -269 -270 -271 -272 -273 -274 -275 -276 -277 -278 -279 -280 -281 -282 -283 -284 -285 -286 -287 -288 -289 -290 -291 -292 -293 -294 -295 -296 -297 -298 -299 -300 -301 -302 -303 -304 -305 -306 -307 -308 -309 -310 -311 -312 -313 -314 -315 -316 -317 -318 -319 -320 -321 -322 -323 -324 -325 -326 -327 -328 -329 -330 -331 -332 -333 -334 -335 -336 -337 -338 -339 -340 -341 +5181 +5481 +4735 +6671 +2138 +1523 +4164 +3883 +5466 +1062 +5554 +1909 +3979 +5025 +1720 +3037 +29 +1104 +1753 +250 +4013 +6376 +1663 +5678 +2064 +1474 +3278 +6423 +5908 +6933 +6351 +5216 +1466 +1517 +6506 +3429 +966 +5747 +2884 +1899 +1265 +5361 +2200 +6639 +3450 +1816 +6810 +3017 +7013 +1725 +659 +4746 342 -343 -344 -345 -346 -347 -348 -349 -350 -351 -352 -353 -354 -355 -356 -357 -358 -359 -360 -361 -362 -363 -364 -365 -366 -367 -368 -369 -370 -371 -372 -373 -374 -375 376 -377 -378 -379 -380 -381 -382 -383 -384 -385 -386 -387 -388 -389 -390 -391 -392 -393 -394 -395 -396 -397 -398 -399 -400 -401 -402 -403 -404 -405 -406 -407 -408 -409 -410 -411 -412 -413 -414 -415 -416 -417 -418 -419 -420 -421 -422 -423 -424 -425 -426 -427 -428 -429 -430 -431 -432 -433 -434 -435 -436 -437 -438 -439 -440 -441 -442 -443 -444 -445 -446 -447 -448 -449 -450 -451 -452 -453 -454 -455 -456 -457 -458 -459 -460 -461 -462 -463 -464 -465 -466 -467 -468 -469 -470 -471 -472 -473 -474 -475 -476 -477 -478 -479 -480 -481 -482 -483 -484 -485 -486 -487 -488 -489 -490 -491 -492 -493 -494 -495 -496 -497 -498 -499 -500 -501 -502 -503 -504 -505 -506 -507 -508 -509 -510 -511 -512 -513 -514 -515 -516 -517 -518 -519 -520 -521 -522 -523 -524 -525 -526 -527 -528 -529 -530 -531 -532 -533 -534 -535 -536 -537 -538 -539 -540 +3980 +5461 +2045 +6968 +2339 +1335 +2208 +2182 +6539 +1029 +5244 +3160 +5522 +1733 +5173 +6316 +6015 +4946 +5914 +811 +944 +3144 541 -542 -543 -544 -545 -546 -547 -548 -549 -550 -551 -552 -553 -554 -555 -556 -557 -558 -559 -560 -561 -562 -563 -564 -565 -566 -567 -568 -569 -570 -571 -572 -573 -574 -575 -576 -577 -578 -579 -580 -581 -582 -583 -584 -585 -586 -587 -588 -589 -590 -591 -592 -593 -594 -595 -596 -597 -598 -599 -600 -601 -602 -603 -604 -605 -606 -607 -608 -609 -610 -611 -612 -613 -614 -615 -616 -617 -618 -619 -620 -621 -622 -623 -624 -625 -626 -627 -628 -629 -630 -631 -632 -633 -634 -635 -636 -637 +2357 +6978 +1320 +4743 +1346 +2516 +4600 +6083 +5808 +6818 +5927 +2509 +1356 +1643 +4627 638 -639 -640 -641 -642 -643 -644 -645 -646 -647 -648 -649 -650 -651 -652 -653 -654 -655 -656 -657 -658 -659 -660 -661 -662 -663 -664 -665 -666 -667 -668 -669 -670 -671 -672 -673 -674 -675 -676 -677 -678 -679 -680 -681 -682 -683 -684 -685 -686 -687 -688 -689 -690 -691 -692 -693 -694 -695 -696 -697 -698 -699 -700 -701 -702 -703 -704 -705 -706 -707 -708 -709 -710 -711 -712 -713 -714 -715 -716 -717 -718 -719 -720 -721 -722 -723 -724 -725 -726 -727 -728 -729 -730 -731 +3565 +5197 +6253 +5202 +6568 +4615 +5223 +2232 +6362 +2612 +7082 +4279 +4671 +5635 +2242 +1552 +1518 +5717 +6328 +2316 +6606 +5722 +329 +6824 +6251 +840 +2868 +2568 +4493 +1215 +4572 +1025 +1353 +5474 +4079 +61 +5979 +3293 +421 +1111 +5579 +2364 +2254 +36 +3890 +2238 732 -733 -734 -735 -736 -737 -738 -739 -740 -741 -742 -743 -744 -745 -746 -747 -748 -749 -750 -751 -752 -753 -754 -755 -756 -757 -758 -759 -760 -761 -762 -763 -764 -765 -766 -767 -768 -769 -770 -771 -772 -773 -774 -775 -776 -777 -778 -779 -780 -781 -782 -783 -784 -785 -786 -787 -788 -789 -790 -791 -792 -793 -794 -795 -796 -797 -798 -799 -800 -801 -802 -803 -804 -805 -806 -807 -808 -809 -810 -811 +4150 +3468 812 -813 -814 -815 -816 -817 -818 -819 -820 -821 -822 -823 -824 -825 -826 -827 -828 -829 -830 -831 -832 -833 -834 -835 -836 -837 -838 -839 -840 -841 -842 -843 -844 -845 -846 -847 -848 -849 -850 -851 -852 -853 -854 -855 -856 -857 -858 -859 -860 -861 -862 -863 -864 -865 -866 -867 -868 -869 -870 -871 -872 -873 -874 -875 -876 -877 -878 -879 -880 -881 -882 -883 -884 -885 -886 -887 -888 -889 -890 -891 -892 -893 -894 -895 -896 -897 -898 -899 -900 -901 -902 -903 -904 -905 -906 -907 -908 -909 -910 -911 -912 -913 -914 -915 -916 -917 -918 -919 -920 -921 -922 -923 -924 -925 -926 -927 -928 -929 -930 -931 -932 -933 -934 -935 -936 -937 -938 -939 -940 -941 -942 -943 -944 -945 -946 -947 -948 -949 -950 -951 -952 -953 -954 +1433 +464 +7086 +2324 +215 +2142 +2173 +745 955 -956 -957 -958 -959 -960 -961 -962 -963 -964 -965 -966 -967 -968 -969 -970 -971 -972 -973 -974 -975 -976 -977 -978 -979 -980 -981 -982 -983 -984 -985 -986 -987 -988 -989 -990 -991 -992 -993 -994 -995 -996 -997 -998 -999 +5071 +4516 +3424 +5141 +3595 +929 +4546 +3437 +1750 +6701 +5057 +5781 +2936 +4053 +4243 +667 +5159 +4084 +652 +1964 +2657 +6676 +5807 +4940 +3102 +6141 +1773 +6644 +1113 +6317 +4712 diff --git a/examples/qm7_hyperparam_search/idx.npy b/examples/qm7_hyperparam_search/idx.npy deleted file mode 100644 index 9cb5d111677091a42e650562df09cfb1b69e2989..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 8080 zcmXxmQ(zbf5(ePfwwp3V+NQQ`n`xVxNp0JIhd0{%*EWy!(irRKIUfu7Gxn7W)T);F&1YDmSicGW*L@cIhJPyR%9hsW))Ut zHCAU0)?_W#W*ydLJ=SLfHe@3fCD**gE@plIgGI<=2mj<>{F@OXx&MsB$c)0MjK=7U!I+H2*o?!tjK}y)z=TZ1#7x4Z zOvdC)!IVtJ)J(&)Ovm&LU%?Od*e@0?tMqyM& zV|2z~OvYkt#$jB>V|*rHLMCEjCSg)0V{)coN~U6JreRv9V|oTK12ZxcGcyabGLYGr zojI73LCnS6%)?;jWj^L-0TyH-7G@C^Wib|K36^9jmS!22WjU5-1y*DwR%R7eWi?i3 z4c25W)@B{nWj)qs12$wMHf9qxWivKs3$|n{wq_fKur1rMJww@n9odPU*@a!%josOU zJ=u%B8OCt-VPE!Re-7Y44&q=A;ZP3aaE{84j-r{ZE;a%S2 eeLmnrKH_6O;Zr{2bH3n9zT#`X;ak4rd;SMIg#vN_ diff --git a/examples/qm7_hyperparam_search/make_model_pickle.py b/examples/qm7_hyperparam_search/make_model_pickle.py index 2269037bb..446c498eb 100644 --- a/examples/qm7_hyperparam_search/make_model_pickle.py +++ b/examples/qm7_hyperparam_search/make_model_pickle.py @@ -22,24 +22,19 @@ import pickle -import glob import numpy as np -from aglaia.wrappers import OSPMRMP +from qml.ml.representations import CoulombMatrix +from qml.ml.kernels import GaussianKernel +from qml.models import KernelRidgeRegression +from qml.data import Data +from sklearn.pipeline import make_pipeline -estimator = OSPMRMP(batch_size = 100, representation = "sorted_coulomb_matrix") -filenames = glob.glob("qm7/*.xyz")[:1000] -energies = np.loadtxt('qm7/hof_qm7.txt', usecols=[1])[:1000] -estimator.generate_compounds(filenames) -estimator.set_properties(energies) -#real_estimator = OSPMRMP(batch_size = 100, representation = "sorted_coulomb_matrix", compounds = estimator.compounds, -# properties = energies) -#estimator.set_properties(energies) -#print(estimator.properties.size, estimator.compounds.size) +data = Data(filenames="../../test/qm7/*.xyz") +energies = np.loadtxt('../../test/data/hof_qm7.txt', usecols=[1]) +data.set_energies(energies) +estimator = make_pipeline(CoulombMatrix(size=23, data=data), GaussianKernel(), KernelRidgeRegression()) pickle.dump(estimator, open('model.pickle', 'wb')) with open('idx.csv', 'w') as f: for i in range(energies.size): f.write('%s\n' % i) - - -#np.save('idx.npy', np.arange(0,energies.size)[:,None]) diff --git a/examples/qm7_hyperparam_search/model.pickle b/examples/qm7_hyperparam_search/model.pickle index a014c8808ec3f91c06993bfdc1dfefa6276603bf..9f9218643688ac3d2aa5485b687b111879b8451f 100644 GIT binary patch literal 4444353 zcmZU530zFy`+xhc-4-FCY-I^qpJU&*5hZIS%2qQ?D#}_yQACnx5n9kDBz2~RvPG-X zzV9{DzV|<8@|oA`|NVJ+J;$6o?%e17oO|xGyq|kn)Xe78X>$|SSsiOjYxC2VXU%1H z{r3;1(9A1LVJD`Dwvo7ykdO`A+}egIdWb2eZ8VB}HnXxlZFR=&0{8P6Q&&hN9jLxsE}Dl1#=WRrMweUVdOHc zSytBOEVjcDwu7}fQ_)C*oXh;IIqQT2Q^`n%d^l@nb;jDt_Uthmrt`|j_H5u z%%$q+Fqpbj9Vx2LJPYlcYdKQtY-!h$R;UUxUZ?8i3PPSLGu)5I;Tf-PN2yE_L-Nn2 z*gIZuh1)UpQe|ZRUA$yn;>1d@+Vo6B>p(3S=M9=^ENc+V+sEB#O?{OK%K3iRPONEy z;{{%v8uESF$HF@4V_NaiRr}{3RR4qJ6F#S0*dd9QI%eGz%Wj3IG9PMtvs2-GYF(TQ znOEzsS+Y&~U$r+~mB|)zkg)$HLAm)j$dKzu*6F{(ryET*7AR3$A)8 zz3OEQ=5l(~V)Uw4(5pV`bYE%g}6^OxHO|L?lUd0-E75Sseeav@?JjR<{ zqFko9HbT%Lm;3$(;Sltpa{EQMT72uC=cdA*BHX^>!}3$_qp-F8>UAfL>mbiyjLRX_ zzc}4w-X)hGQTTqB#1Wru&+#Hh*~TyNo#5DN_51KjE|h3|7P>hj1c#YD$Xe*s1LfKd zws)@7<0%`YJQioP<1hWp{%WskJW!>LROPssbEvCv$NP3HKF%9T$28)+CG+BTSJi{X zYVFoaiw-F5pYYajawUG|^*pbC6c@k$@a07Gr4D%fEhy40s16%gjXM3HtR9$a+01nU z?blQ7H!zqRsrIr|`%OgqHBw}^5S6|)*{p>NWgR*t9)&Nl_Emiat$9v(y@l9>KmKCq zNy5UvVbk8>25$?xvt>Gu;ASwm9DtD;gRy1fhU9+Q}u_q-U(dN+t1uBnDY&#I9^ zY%^5qts6ZgmjnHpa-!ZvRj||UseOibJ$|(m&D-cv4ZD~V?PJM&XF!1J>!M%SN@wVu zmJ}BU3{U!WE>sddk8eDk4k@5K_VR{$Mk}=Lovm__%pdP^aBtgI19N}AR@_%vjg@t( zl!`Vfp4)Qa=z>+J$hm{e!D`UCF5@~KA85QqJn&X z)5*YNtcE^pnfhau$K7ts+{|Wf5dgN80=A98+%8xyJNV1xi0pYQxg;(SBMKKZKF1-c zw(hfaqOs3khNQzGF(e&YSYS9Q5_?AlJ|4Vs6}$C&Y;%>DL2udOg*UYRVlDqkg$v5K z5dA59(maJ4XdSL}mVTBA_XB6W$oj>_1v0+gM&qlXwJO`}!1XHl9;a&fz^NLR9==ky znaoGpWFPmt&c#=LpYK_)n*(Dv$7iPbNufu_KmOV5^b=NnZ!c?^(~Qfr{XV6ylR(i0 zcN>=D8oXci-0lv}Ff1%?3{E_#gyNp(dVWYNz;iVwcbxX-0&^#uxl5q^ZmK zU;K06ybamOIxz}8y|`=aO!EETv!9+1PI?aC3zc@r9Ucb{kpY{4eauigqfQUz4U{eQV|N8DxgqYZz6N1id;I+fA&^!mnL z-&nAXvV5+X{v4Ql+01 zvNlbARxEZ^bD0*itP8gHJil0?N-zzJohRFT8Vr^;&sFa{g5c=!1oRgt?g@W!0HjGBlt%hlvysQJTtuw?ml@R zRdBZ)_eX8Mn5x@=nFrX+g90rNQ7sQMm`A9VGSp^Asg~Rq8|qFpbhd!w+VU7%>qyY% z?pU?*T@-wm=P4ZcD~cw!d52Ga91JZInZ}tlZ^7q!xu)qiNu(4G$@9wEaq6#23Wwth z@MSsmeqIDm8hFUJd?hZKT$dSqr@Qda`3w zfZ9$m=A<3T*m=~>i2K#eST|XO8Bna3$Ao;s6CY>qB%1_x?l(=a&3M#^2y!SKYw zP@}q51f^fd{pNfm3^XEH7I(>f`t&<16vxXUzT3_?|Yugm8KgUxK`T%$GRAi6He zeWp<^j{7`N^jozV^G{Fa_Xa2DuTkVda;fg|t^3k|y>`#JzSZT>)fBgF`_&5k%3#;- z&e@H?Ji%sK2(H?aUiC=^^OV4rr}?(*)B0ev`&u(xZu$PGS>_3R(lm?hz4RO}J!!t>vp@8grrh>+Z|iJ{4da zIrPiqA{Tbnb!zeAvhgq7Ex(V{m4J_Bzu`XL4E(NrFvq^W87w4tSJW?+L1D2d`^b)J zJboWY{w=&*G? zUU)Kb%c^q)_&|MZfzizuXW(CoTe-L${KL+aJBK<%jA)4!)CGt+SpQU&M}yFR2#lSVo$Zn zm#bbEFB8fIlVRJlf8F0>&+f1X(%2d6&Dn71l$96?Z``$gb6zYaPsJ=e?K!5Nru-`@ zB;T>C@Oowo1YWlyu5|;p+T!rq`D``(nj^hy&(M!dy z>11NjL;WJ>HMsb%=cVI{ldIv%5zWntZyJHP^|N-l<$>91lU*;ElIXf^W$|4e7dKAu zDbkM4$7TME>EE78qV(0$3(pKk!?v=#;uX0q5I5NN^_H*%BF==a|EorbIB4}a)UzKm z9oWqC0_`tQ?Hw6RCxI<5@@;vJ{7BT_7FgGNb&T@7Fnr4>H8RELEgV~3y7TsXQPg>L zr%l_`h1hjG-B&(>-U+8~^W&t@^3P$beF~dlUB6M1RZ1Crwzd^YEv|sr8`?9hF1KR+ z1s(%QoLszB6v9k!1vKp-?#PaM*nnI#H%%|Yr+@A{+(zb?#vBS%@Qwm&nc1Ta?(;z2 z_4$>_osuY&^L4ScYZF$_%s&^q;-=*6+XeQ6DVQ)UPZJmN$o`T-C8)jabon9w_H4ZmEF}r zLvggxu)by4gh6;dua)`YD;K-%@HTEG^9E;wyleL5!P6++0()yA4%3Crye`n*m1=*3 z!MsVem!nUfTSWUcQ;2i9<03spE=1I?UOiKd3ro~E2d1n}!R9X~=d9ntg*s&#DdplI zKf3RxfaJqFOMm6(V(S?%WTVcLMws#%{I&f??N~2SS7&d1CCrgOadDC45R6Zc`Tl-Y5zH|;{?P$zEZLhH1Wk7qsEkArW>2-E&%d21@aDq=^+5p^FJUPSzFMy>c*hS^Qd6L>CjQyR_j&DbH4)&TQ9ep-6)R6Aa9H|z8y$toCaB(z6DJ^&Adr#Xe& z)?-7l;R9dqG{MpRmJFF?t+@Kbo@GXts^Qedf#>U;8UWKxVs$;tzpv6fYg!$6-IhGi zG)4wV?H{vkV7D_)M^6tWcGklJAp>Tmw-_=g|HVCiyA98vw}kUBp1LGuR8Ns2YU?PDZb!`s@h+O6d!RewZLzz8_tOfy018=8c$oJ zng5|h39U>08G1~p1M3GW3cK#A!B6Cmt&BZi0Ygq-Ttox8u%5;#c@<#w(cHC1vJChx z^sHK}79?}&`jSdmlv(@CUB3)me&r=4l_bOXDVbuo@|$se-l_AuA|=qz6>%@n@%OOu zsrTd$GCX+Ur`vV-lLT5l(>QLdPYXPYnDWy|uNd>gF)uOXIHSMT618I{@uW-I`>^}9 zt!aL)V(7CNiA>#uIH$~-oFA5@;&xF(2a!M(v=gbAq-5t~_v`K(L!H}kfo6_eC_@E_ zbRN%fTqMj{GnH;4?YKE#Dz(v91=-V;G@}DArQQ|^@0btSOdkR8k0|iI4CZ42@P7XR z&-gIOrEzyR{La3-I_qUM%$T6?F2DI1Zmg2~tZpTVektDB-!<(oFfKBhvw0aXmgXqF zNFm}exW_ftV!dN4woG?y2c5l^Zl?z1(6zI=mhRqtu!-&rdUcRX9=1m}b;1>i$ICpv zb;Hk$AiG&d>oE2e6|*79k#67^$2PY7Td7T8q4{U~%!n4)$v)JR)cXw=L^~e3K3s?m zLJLG+BUv;Od0_~o;i+2XuQD$zdBK*mMVE} zU~CBVjZzuK;3%Wq{c|JJQoG@$xNe_JRx?Db3+R(%N~4SXTX7i9e@Q%6xw;>iPuR?- z0^t2A@Bs{FAO)W9{Ra`?*K&)sqw{v1Y=rqg*(BD=#e-XFQ)0w!;?ZG`-S8P9bi{Wm zr~Xg|#G7ugzo7UIYf^VUR2)qjAO3SVya}F4$lY6Yh6@e$Kalwf1!QMQyqKtdNSzO#%2iD+_+3M zMGOVKx3~XwzM!8ce*@UeN^NNtK21*M}jT_&zzMO6Z;uqh(cE zH=JJkS^tuXFYq6DXOvLNwc~#h_71}a>Z6746S6`=AZ)gOwL_!`daSWmH6r>QY*@|h>3S20Kas4UpMw;=`Nzu-9V8L)s_(+MNG!5B|Gn3Z zMi@A=(=%UQ3He9t*O~C75`5uBtfibZhZ)Ldh6%I`r&_*ZFke$G`7zELs-@=;@!e-c zwSqzA^T(-EB5~cx*r#r~~pLtR0w#%#hpZSDkX_uQi zqw)Wq=T5}Rm}$Gl{-5tNBiPJH!Bt1mtBz(cW9U`$y@t2ss`*|+pO1nR`$a3Ht{n%r zFFb=MiBAat@mPGpa<`ItQ51iH20%}MI8hFJ?ndMBZ{#n3_%4l}%!w_rdi#4M>5|b3 zHFmx~>b_cfFP2K z*{NR#J~yWZpFL3tDsr`I$H_c@QOs_toxSLvY+bP5YO-Dk8DheDPC`c*Vjb>o+EEAn#7a zMB8Jks3F=X-*I9GUL%$MaF=i={LT#XW|Mip1xreAv=qW?y_w?0hsh$3Wi#IifP7DZ z{J>zw34r|gACN;bk7i#x*$M-!+n1|`gFN8nr_@?Z+x*8YxCk%nxGE**H8bl zV6Xzhjv6qu0;)lN{^9npjB1$9e`;31rJMGR1BF~{Vr9B%$BaDanb7>9+^z!RDk?&X z*2Tf1=-G2MZCmlZNGmO+)ne#ggc!rrum^NUPiS8j)CIYc_w3}Wweyo>}YV=}q7p{o% z|IFx6M9uGHf7+^XG0t7aTv}3xnelAqCjpQN6v#vd6H_4hdnkwAM_fpcK3d_?1Pa=Z zSuS>gaPHJ2xwQ8E;FxnHB9bG5tPlJ3oZJ`)>u!3lcwQVM%d#BwI7RAaNXOZr8tYw`PEvlc$@ECr(}=NU6>i4(2#tai-DAK3U} z*{q$OqfxojULAvr)gYpi?p=MV9?K|-FML7fqg1~dip*@qz7JDcmS0c7{GTQBfg9bQ z9;hn-i#1-U_gF2!Oky*W1*oJ@R8kqt&w}Ok<-g@MiTJq7bh)5>iZ@f}OE@lymvE|Y zP6L|P+9fB3_8s5SXeAv0Z8Gb2ehdDH4L%$6NuQ8JAG@u;GG8|0=uIlELLTjyZo;m@ zawzP09s`jS?X`i^RUf~KaBiYbg{hBfVc)v#9laqHa69M1td~w)ShoLULTOzMST6QR zIl796y@Ug9O?Y1jar(2C?A}@k3;FrS-*BpZI7~CG8Xu9&7{B|x6nZ|obXU*mU)X2# zEjGNYhsb*##o}8ekmBmxx?CY1bc|7GT^gB-qtvzvlEKjyvEyme)dexgzi{z0Z&9TL{;qZtzGC3k-n@(dh(*;0& zr9ggTFux0y*AM>kDyEwUNx|BtnLakq`v}*kNKXIp_BLLmrMp1*k{EhFO|t1^+beh* z?%C#bG6H){xVWVFo)nts^g~kYV-p_7FCKC6#MtDHgfFt_nv0j-CAoh5^I%z3;e{TH5XhUBuvsD<#l1Lkb9^={MxBpoa~(W>h$DPXnx#JqLv)IBX?zV zTtf?<%HOY}aO#W&Hdns4;*Z30e-)MkMz!NS6G~FR>S=@LTJUr4@@^!{;mWGs%gz2?1)Q_KbTqycbz=Vlbk!@Ys6&JYc2(NveKTZHz$$c&ap{n z*S!H?m7`5#BCBz9+xMgHo5;3b^dz=Wrvgj-|(QvjUWRb)M^t@5Cd^ zt`ttaKVLR?II4zL7@zF1IeW@~Z(DUTNQ~%*96fo7>a$%AmqVEB}&u4X%+` ztG@3;74%egeb)*v2mbRpGz8BYFBJBH*NV4H%lxkGT>*x_Hn&bF<${3$bwBe#xBT~k z6%!@Vsu>L5pOv40+rLD2#|o0x_q?(q-cJmDs^zepw{_sLt!d2PAQ_x ze@!Sw$pw$(PLfC@sip!p^RED~LJC+BgIP=g z8%2Y`5(=0;Sz-KS?4f4hhaxe)FleD~el{@*-^<{ne?KpZ3P{Y;C;kiv4=P0NYmUZm z&lS4VkYF%_cz-d<-FWPksKhH@lYlsR#tIQiNFp)lK*Y>$ENnR}`fBfIz^$d9^m*m< za>c%Up@tKe2w$g}g<^WiV}Ky@`cw4*vf41*tVO!PR|7(;N1c0xc;%SY3%P z?5uwK@>nCDtJ(SV+gAy+n_sM|1O-|*zEuNOmmjs3xk{jy4Kz&7gC8G4uhiIZ@fhc} zRd+Usp=BZ;_J7J3;;ak$9o;gy3cD>KX|QrJL~>(?5@ZH3EmjLf@gSr4&BGVkD(LJD zp0A;88}6_c$+kRHi_7Lr|Gs&g3S#YO3Urba_X840yJm3mTsU5bFFZGr7-@vpEvDFm2Hi6CtS+1M8j5iqK3u4h0` zELO?k&VAJ{iB8Pc?e!mc4Hs#Z{6_~q@L*ln>Ml9dxw$3s=8spfBspm64lNPR4yP#& zwYv)7f}%*TH_2f6U#QHw3($QzL6f3XEGn-;rkL7W7NB}FOM z?UC>~{oYbiIr_Ee; z4?0dYz8r26LjOBIr>)U7zfAW3%s)BNv~}?v<^MDPSxtwv(=Y}9?|GPwY-W?->YC}* zwJ?~ig2l}JZ!tes4^3>}$_0&8Bw>Z-;o28J^+)ZFg}L?Gk};db&^Sk5ud$B;@Z$;D zD&F>=AV(`MWd=(MidHu?1t-pFR+T~$$0F=5tsv#9y>l+^G2-IP zA6}qn%%uM!BBr!TL{s zT|OVv2a0vu{u&T(ddFItJmSD_>isMc8S(0U6cInU+E9k44w-EActOf+TbD0-ct9M9 zRoD(*Jv4y#t^e>fWI+q2k5fEG5rOpsEyb=zJbQPSH)|CaGkI)gn*iW;3UCL5*+~KB zr{%k-`SWSn|L3}7AH1{csD(tf7E>Gp+bcf-|T608@V$C_-?yohlyCCvKw01s- zI=*+P3~s@vE@$m}=2->#vf`&DH7lS|xZQK6e?90J!s4$URY3E28%e5`l{9`}(DWTh z&a{__%aTMSS+XL~umu{^?~=_d2J@XeGXJ?^iOU0-YOLzBqj9lzJ!W>ZnLPsSd#Uz) z3}(Mzc@6xxyfpclp;k!jG)dS!`xTZ^jk;mtdIj4kvfqY@i=n$0`PIetM&B)m*8AZMT2_%86G(W9D5sV7IUnlPX2ay$^LG6)}-2 z7zvIV!0O4R_$I@09N?jq zw0YK^AFFqD;_*us{dUawi9Hr+UMnKC=!uUu4jnBM=Gbj3esCzV1DJzs=8yp3VFEBq zh&*~2EMX!k%+%>&iI5KzY3lckQI#-P`^klYJLM#934``kcWTnVCEmW?MfKjmT*>STsZ!cT2uLV zR~we@Kc}oj;)yBK&z-R%aS5rpFN*Qb9?4{tiz5EYRqDi%cPd@CXb$2VdAAm)%~V9i zvCfXa+uL!ASbv1go<@wujvDX55=BdMjeT4ex52)+UA=!;LLBD;v2os|3aHP-nQNih z32ldIeT0Wuq67>!ON_|>@MDQnoCt&xWMr@;iTrP+vM6-_E!M>e~&s zN?4=iEdSA&3%g!qR+-9ipgMP{#4+;yXLVORl5SPuRIMkQ=H=GF{vB)TWlh9U^6s{V ziPtLNK4Fus@Dk$8u$j{uW~Yd*rk?8-A`Xf(FZ;t1%)^t;x?cYrB93MW&zBzCHVm`W zPF*_D-->PgCyL)N6-VO^+H4&Do`(I-oZYr_oDdo*`iY}@M|M~-1528K$Y#k9DS%v-;z}TwBO`+~ib%n#m8XaPt=hlt z*@}Dhd0<6voBoSP(dS+JujOpZ%w>MiR}?+dOBcT@3K^^x928 zDrS3+vBi-byIHnQybC8Rxh_;@G>Dn^YCJ3nU&z+($*A7icD(M?`p1}cBhq;zJo_2uOp5-YyW5`CVKtN%$6p0i-s6=rh z5Gs?A!BQbo0HG>9^bbNl8&JL#7LLDDyD;DjR=V9X<=l%@eARQ#QIWBtNRs4YUk4=t z`>@%=kdpI|;!8+hrzB8eAIZihw!_UMQPf=kDr*{fXe7e0v5Vs1SwDQi!WDlfv(0VH2p$HN>4W<~fR zDU^I`&uf?cwb-tROGq7`u|`IU^tKaX$m!nl(^D4~!>Z`GlTV*~!p&Xvim2W+zU?_fazPf?GMH&*)XcrVY5}{$(F$w<<90x#6R{ zITKVIvpWie>tTgbs)N3XBs%xpJh&vg8AAP+j1Ham0;ZNe&zYVgfszPq!5WofwvVnZ zC$SX^iR6VHY9ECjm=ui)>J;KskDAx6zCr@ACJ+tTtcgSl@KC2p6Fer7k-?fwB*lYf ztyr}1O7Y-^1=-wluQUPL8$*I@;@DURBi)mV!YInYZy!Md$hmvq#s;pB~31LI4$u(xnK zce_^~HcyZ0lq=5!{cU=ilCrBYOM__6W=$nhfRiQ#gWxocj11OvA_X|jpob_v}iFGZ&=dP0b@!NBrY{aoH{N<(4xL*MS_}k}kMrpeh(Y(10UPXZ& z_?1xUnUW86n59KDXS1}46yP+IfCfR<-dzeQpjO6ei|#NaCO1f0ghe{Zj<5fy0w zl}B-$YMv6RH&^{W_cae2uC17~b#fD&p_}{YY8YUE)vPS?d1BNUf00^D^Uh_nsz9y0 zd9_q}CY+o%DrIA8GkED!E@B0IVtb7f%WJ@QzHQIRDk{V@GdM(&RXguHulxI@6-IJ0 zSy1xzW!k%KB#FgmFvej|k^q!B&HI;{6d$S9Z5QG^G&r*)X-NsZd-86BSfw~}#WPOj zq%`9mhW+oR^MgQYlh;GHt03LMt&B%O?RfD@ZSD25T`_Ae0g27hAyNRXF2#>PJCBSE z)_fuT0jqxCCG7U9US#{b~D&?^dxM*^bY)|L=~?@4>9G|Zm1JREm0kTvJDr2 zPc-1&#*K?Z2J636484dcNg; z!t-1pRA(FRwrm)+Q17W!MI-PD?|ymkETd>`a`zjYlgT7W`4-H0>aeuwsRYu#YHlarj>KC7%40FtZW5M%1CC3|AQ(Y6V_I?;sYVhn#S_)>oUdBXBysS z4&Vc8LlbVCZi;9=WnAvA#eM6g`(y`sII52@rQeC8D8}mbxvc~6e*L3S3pO`k z)=~l%o3)Hc0kq30egxVTWMr^b5-Hf#R?$OrR|_6G?;>v31e^Ky+&y4f*ZZ`)KZj6f zt&YBhq9}zVHSJ8V!}TvJtNmQVM`)^A(#V{3@`R6ZHS}jtHXsk*aZqQbOGqKJoY<~a zd&@An(e}x|Zh_@X&U~@FQ49&mSNnxM%3!9K+Fzm54Ve9)g)EYC2-$b!&}i*4u%LB^ zt1VzCkt?A!wH)HyAFElAxcm5yho4hxh+Qd@rCAh3k+`=hhzH*Wo>^tF2k{!3FEW-# znGUxP1)YzXy(V-Luq3p>u2_NaBk{4BS>gJER%2ckKfg@_d3w3>oJ za9Tq~25T*m0-V;-Llh?m%KsyU_LtCqcj|&C(hDK6d4UU9*8tJf!o6= zbKM`<-us>duiNDG?q%Q|=I&?^6hOL@1}U z4(3W#J5BtV4aTPCdTsk9k@)DnIIFiAPnh+k_wlZ5{6or$kuG(I7642%pEIUWSy3RI5&C{u#oo^+aG$Y5dt$_QDQvSk1H(%eg#XY5(4&k`rkhI8E(r;u89_E#;^Y_(s)f$OasEy zmB_?ScC)v#$M8T6bGO{vD2X;Ssi64_nxJq1H>BUq#3T}MSw2=7@sl}=g*X<|U-w+U z+KV5UMSmFXQAVr0SFu;G?8c|wwcQ*38UvIvGh18*?Z&bPVjg7Sh(m<+Y)l;2tpqGK zYa5XQXtz`R2(&xM$YAXxQm_OmABn#NjnWT}Hz?%76_SZ-dJ+l=leXRR*vWw@zaPB~ zH6{h`54R<_V>o#G(`er}T3_(1#lsQHSdvKbV70WfdL6`2>PiP;o=&e#o%D_LJZy+$ z4lb{S@oAg7q!xAKDYBUw*3;^NzqH2H;y-70^}kkcgoLmdFGD>Fwg&N04Kf?4jkGced8tqusJ+bb|T63sj6_5=r)&8*s@E{c5Z{DL&& zg7LIV^SDJ9La_Xp9a#fRDHK6DG)j&5?bmXz9aTImqp~3;`iT^pI@o1fHnIBOx+3f| z#|wV~BC24-u_TUY-1xJc+ANs!>HB)eS_wqp8sRZfuc ze|?SQvW!9&{aikXZ4XQx>zAZ}HvF2)o-Wdb#m^T1+@s!vcc0V^YPJ_e^WwD>WET-e zf2Q!}{3Sz}GJV{|712nWSld7DZ3|`@63y8xBO(Pj?WbT6oQ%oHU>zWm;zVg5tb_Cr z!O1gd+O(|2Lrq|(^g{dLmlSZ8t175I;Exk%(@B94TI(@x<@?{UV8v&{oy0!7lFG%S zBoLqPJ(GuNN=|=vG5jz!ZD%&iB7G7f=c@>DHe@6}J!jtllY7R!((e>Ui)7y}5B6z> z<>k}$JFn#8Sl5JIoW+VL=beeF*e?lAAZg{>eXk8zhsc$&S%--f=y8NzEYagA85t}S zA_ev~rH81!ztf(kjl|)nY@C_N@cY>Cic#Y-AcqttdwWDs$-Rk@;a@RsEx$*Yv2St5 zjZ7O^Ct@cTQH}kGb~mqoWQ*s>A({3`Cw#3%Ie%@f)|-m7;E{lA z6rvBn%8$~Vkw-}_u*`_&Y}PR%1vr^gFbGb^$;e=xAd=$5_Z}_iA%YX9bm*You}@7f zQc&{*KD|id8MaZ7F73Y{f<$SYxz!aYd#QhK2oBcYsIG1(f$rlA8*fWj;^+C9`Z-n| zm^9af#(k1Te1?M!5BDn6zmHfe!f|QVGIR8)1f_`&nR`SF!L*}jjp5XKFy}`$;%I37 zZ=cOM<+y{BcJy;7;VJVA^5Te+O(QRpz;o(&*6EdApg-3>(=d|Y(b=KYbFBv(=V!}b zIwQ?tSrSdytdm3v^gczEBYK}EBZGB@NUAqK;6FkQdV$L$F_=^9A!s(kTDhW*+Zm3-(Yqs95%; z1PbgUKc7_~eMF>kzOAty8&&;Lh)OB|JAPzX4OzFWv>o5I07D_OPQIcXNBr{ZoVl@T zgpBDJIFdfw+yY|iuTO3!^B1-&$ZW~1$7J(xRa#St$6hK_nU*Jo0tjWkB&z{#1m=Dn zf2$OqpmE+o6~y=c!+UXSI4iI%q5)V|L_;>qnn(d2OsX`&gGELL%Z5k+9&CDu;vsVJ z_O#U>n*VX}??S}U3A|%V!tvIxEafeiMA5iaQ(vr8e*KR~bPxZ$C$uQmT?+mDxr9_C z8*uo&^i}B#8Zf`@)mRFZoBUnU{fFe@{ZwmK5?;?`zj-El;bs3ik!tam!YvU$H`HRh z$%4=&D*#biusBJe08RDg;WvLeNOxZa+|UWJdX)VczbW{(b?9Xa@X4~CVu(88j^DcR z!uf}8~9FYP%?5NTN4|_5)SPn!A z@HkHoQ9SrOJK<&!CULL+Oh2r%(zoI3fETdWC2kGoyDRo!+BC4%{Uw}Jvkd)2 zo(oqbY0|GAjx;6o1g&g@RlocfhiNJy8`_JU+yEYgr2ddp8plvR$&(5oi%&1#pbjgd zlPX%>HJDOV9cPpR3A;^H#6qE2PEK zMji>Y-}PBI+ph8-zrO;f%Di`+<|2wTUVTx{ljsIMgQ1}tSQm(fY?dRD0z90k(gcr- zWMr@|5h=jqGCf4`*xdQ#$o(l?IK@qzHFj1!w&z_ey4d>zNE=I2_cT%T+ErI&{iZlP z(#7}+!#-}IqNX^i>Xi#w@~jrrCpK-D`^v+C`J_#;M;_g|+oW_-q8@G^t?aXN>cal) zsYwP4n}B2|NOiCtkI|bOwLp#dT_;uV{b(ozetKNC1b_1OR5-u53R6DE5fbf7&~`qr zCfKK)HRvZ&286C*_LG$MXd38VS{G!b#Ri%d33Esz=V8Sw%4p=ezX9dMyL@WItSdxA zHtQ;p0z9r!r3oI+WMr^hh@^P%i!Imbp?^77LUfClRc(c@Ucskklt#jYWB2aHL`ICb zHpfKK<*RO^7x{Pt+rj3;z6+7?X8}phW=oE6eQev|FlpDdn^c2KNM}@ryAmQky2Ix^ z?J!PuZB>VD9ghA>eas374YGd_`?Cr+rjeu}RD-|R;@YxDmGE-;&J`zQKjFWh^|BK( z8nOTKq@l(l3B>Q1@@m3?f6LZ&oce?B{JdjOI$i?xYb+)-j2Iv+aZ-=w!KQ&*iw>re zQoqq~?Yzmt$cbOe&%pc?r-(QrO+&4FMz@14pL5nf!tZWVMiu*3%dzKnf_d3x?>9#= zuv`f!Y}O4T1rXk(I1vbMk&(f2Ba%X>KwS-YdWb+sA3^*MK$AM$Jw>c)D+zEMFJ0?u zsQiVuZyRgcDcO!^u4y>e^&lUrRMo#$!*3Y!p?SXRIQvvgl(J0+(6FVyhJ&-Or(|({ zwZpPQl-HLJdE-mRJ`8h)4bdAz#%I-lMxDjpt$kb^O3U}hk|8O)I%8X8BUU~5DkO=` z#d{hDH!MHfiQA6S#A7G!*074};|*Z@$4hp|w5P(przn4B9_fKNpW;-Nmk00H>F8g( z)Cv2G8x|-*#6+HCWU%fMDKODJdgz~tZk&)_Aaari|61KvouHko5Qb+= zJn!cEUJSkI_Qh5oX*3uZ#G?<<0gY- zO%{1aF0KR8(@3i8ZI~a3tW-uEJA(DM9w=BgmH?xCtl$FKIS$_=;$)u>*g;D>>DTiIg5DERABSSSL{93pgpJ?l~By1wG3GoF*L#sD90=>0y>*@pGbj; z9#A_F6M2)7!Foufz(hXu&_5IT+ursT9qa;`z-v0?zvJ=t6I$1vRTbbF^i3Ea#nBeo zhkNtoDuJuYdj3>A6mO)hHLD~M-45Da+kw#V(Wbl(I6b5Ev<)FCMr5)AT>pB5yW*rX z6X!RB#xk)xV#`Porqg$lYEB1Eu+d(-*S#0`oS)8C&@kGoi?b{7k4u|KUbz?)9;g2h z*Qo)*K&tXHss#G-g-2m^*pH;Dr?~XOG2_Qyx{S%)7jsyyTCDwF_oXx%O*{oXlQx*L z{z;sxSQCDHj`SY13L#U4W;5|&E8zT-kcn4`oh8BV*|0-Dc6??qlxL=!J`9enYk z6Vnb^%W-1pgv^sa`PHJR`uEhr-izJ%&ew_|O_d>FJtCm9S-wOHO!S!AftbjTj11Nj zA_XRTN)P=r(bA?wGyA4`n0r*@eM?OosCq5o{IL6k_tRgzD1vD7o!-nixN+r;lW@Xi zxWC1IO@u!o9Nnh0j*sx0w)muWLC(cfT3eE(P)y5Y5|~%x9hBO+hl@W&E;``wRSJ#E zSoGs}51}sn7Vn%(=v;?sOJZOtSZkWQ1s8vSI~lz*?_6ucTUP{}{x&)v;tZb8eZ#B6 zg!C+l2-D_z%`#GS&4=@pIsCL0_V7DW|G?SggauZ)qKF?2{HenpZYPy?yAOj$$gbjp zMWV?0TEb0@xK`-graRx;au~L}t>V2ssDQ%H_dc1i@FQ3sIADKxR`)+%p%`*jAcXJ_ zZJ4URs$u}A|0XHZ0R`k#J9c_qL>s<2oswW`Fni{kS^a0lkOSq(k!G$fw7+=FAeNzB ziDHU~A%9%us81(!k`(uWni|aVC!n)g0YnN+6iDqrOcX>$2J0D-0uu$(L;p;~ZyfyC z2C?H`YJ_ng;@ySUcPp;=0TMcdpFupC3rR1RI^2zdIXS*H>MxTq>48757b}S_oAkDy zNlFA-y;$7dizmb^e<+~&A1NzsgvFB0g^j{NN581P z!bLiIe=MsyAnOic_a2f)!w*kois;p2-!-+1&fUxc?xK8J0%+y((u$#7>O26CQe=33lYvJ$`Bn8NB7T^ z$>-GbaADF7Cy@|YPU;b_MGfU`uz)qD>+Q7;SaWdcg%f%roaf!QBUWaUh?i6chOT$P zXrapSdU~=P)^h?noArW7fr(yHI}j6vkdeU(C6by*iMk76^bj!--SGC(7pD*wcN_G^ z+vWIR&&@Qga2;yC1^g{LCWc<#qtW9_9Faw{EB<&in>01%NFigh-_n8Cim@1N=g{iH zYY4Ar%td9CrFLzHlzJ%+*+arV*>245U<^@3w`rQWwh7xVojzrTi3F#evT{{Pvn;>k zD-(vJ6xBu2T0vZl38=+% z(M}@`Yf^k~xfZ*}os#Yis0H#ws9k(T5{-9Py18a+3)Fb;BSo#BF!Ej-ag?Jy?NRRX zA0zGKNSH%XQg`Lz#1Ub7Nd8#cjc-pUFD2U3hx@WO*e(-NMYI*j{!JUUq1hXA@{$Zn zgpyW8d|uqbGN28MmZfdjl{W5eR6#3c{yejG8p3=U&@CxWemsVovqDJoax?B` z6NdQ(2@XFf-=%`K@w2Nx@E}RI>xIa1SmDGtY}P9x1txt>K2VdsAtQqoL8QQ>k@V0% zlPEi3YF2No$HnhjcV8+d)Wghkgbno_4{q09$IZ|!0cm7C zjb=&;@QW*xb0bSzfd48Q<32dQ=dO&IY%2sWl6a!Hj5xn^@Aevnro`U#>Q@hg7}~%; zsT&{J5LJ{VBZdYSjMG&9EW)96e(l>0c%-Dz3H&z$v~}VvE7Ah6e*izTH6X3lT$p@Z zw{XXVeoR@q@iQW^Xx5OVQKKlj?Rc1%qgjdXoFhz!(0-s<9)BBoM6=C-+NEHkH@nbv zMF*_*PxlITQ$Xbprv2(1Y{w%6-xi1_?a`5++hF-C>+c(t+x{`SJMm2Vx|scqFsPb# z#qfK*fD#$Ab|A(3urO+iibtDQ4*{jTl&^k|pVAk7_$ncb!}9U1ArVYw@huvVz%F&P0)}*G@ahhPKTYZ_EcZEF3Mp=5hJl#(L@SN z8bj?yO!}6L3|1_W0+YU@hyIz=%UZ*uY(p#XBia7XAfjm3Icq@y@cUnX2%|OA-T`<*=;w zrM~>9A*c(|sNDTcjHA->z#&|x8-|`VXHICTgxI1-(%0J?VCA3g@+lexc-jSLUya;0 zATK{saO$puC$!jdu>-!UFTVRoHWhSfzw^>EJg_HhUBQ7W5K5O;*+Am7(>{T`b&4sl zR&BI?fdr|}Q36%tFf64xLAiE()^gFbmZ)aDZQRdqT_#fK#uie~>utcMs$w2iRzvWL zWK17TRYd%^8f_QnwAs-H*gS0ev+D7gizKuvMZ z;kkClua$4^%^#O5;&bw4gb4HWb02zu^0_#075}WmbNP$V?sNIN3Ka!EbiYJnqlppGTqKAk{ zrN}sU?Zn24#eIqIy}w$oC~XJlS&J?Nj5Y$xFs*iv#5P#>!jUj#a=>7(#T4(K+raX* zg^eG130iV+PmlBI7W~P9zOG<4o?h?!b;*HR9690L+;3$5QJY4HTEa3+sc?@-93)Aa zfPSsd0Z!D3U*;#*fOA~&cPExA2=^+_(svw!J*2)F(b0%U4+cKc`ZEnrY_Kyqb#4d( z122(ay9sFf#>uC`sK%&duYyk%nDHJbZxP=MS>s+m3%l6~!Fvo>eJ71A{I;IITX9Jk zc_qfpPQ0VVPFsdJ{~8gJq#Jq<=06jk3UelOJbw=ls{kfWBd?X%0IUQ82Ah>gqyR!p zaUu|M$jD$N5lJDW-2$v+dWb-%PewUmqLfDmw%Fh8;ZipRKOCY(X^i>4s%JiKkhHzK z=21HyZeW(4eDZ0e!Lx-tV?1>Oudc1eBd>3Y#P2q*H*0Qfg->EHS9z|e20k%ADH5+* zLpe72Sc+DZeA*#~lzTqve1dc2K{)=HFz4!>eX+sMR>0S~%hQ?;lmT62W|DpQ7~i|! zf_)c2 zh?jW7q!DSo`2nd}1&Yp#s(iQ=MVaYlHNaoD>4XWOEaB_p3*X}W^EWZLc@1n7Eh0o2 zlZ(m*Mz+c7Z8a$gMLlY!2HT!Lbla=biMT^Zt>RHUbge5u9PtbDK?QSSM718HkKHCu`lxSpw^n~l0D)0OiN%vw<%rPd7Xos_cZt2d92S))oKmgZSoH) zn|#}QX_ISWSJdu^I`%ES|{E|&Wg1t`j^XDyh3 zKWtBdh<{7ly-zAYp;2+WhX<*8qtNk1V?(<+u?>CEFlcwnY3t3IF@23EoqmiTDOEOZ zWj2h;&o_*kVa9fd4ynxp_7&rfTlz$!6>* z9PPK^R{;oM&yN~BT-tE@$u2gmpFzKwVU~pum#}I1^yAqm{i2}ArC=J>Aq?Shq7eS3 zStfLh=}^fE&mEU~bcdZM*hm!KnEg%&xc9W(QCB?(E+6 zHS9c9k!54+2I8EMX~0CK>QEZ;*Nb_dJY4->%GQ6gE!m}>K8k5Y>G+g7gl9XX;VRzB zZDq?|j=0rf`;Krt$2%E~%CTYFSuJwtv#rzT^|8{VFc573^~=E#pJn_SPykea+S}kw zF&cb}4F0h*3wmofm^VGC%XUuOTQlesRV25!JZ<^%5uCnFf}p+nEVy3t{lgu?@gPxe z4!*5LOY!~EWxUZTL$eBnc&0-xHvL-P?=)3_XFG;F_-hs7wFAy{KmUZio|PJR^lVfi zujd2Ak*d`lq@xfoo-qyB`Q=}dR+L9!>9kXoQ`;A#6h=XmOW`!CLmt87N+FM=Stdo% zsJ?6e@=JAht*2eu+Cd9SfN)~P^JS5+XD+qtFIGY$PQ_WP#a>Sxw0Wf46g2D9Xv1}% zx4=tTKkgZ@{ygE@BNwJ?zWMiReHt1+4gR^p){==@?znWQAwKDwzf{CA+rL;d)RxKp zZ7RP}Y&~ zW+AwrUOY9giw@HpWz=DJh7t+Cl<@Pw0!Z&n<;!z=Y-I^+-{np%jC~W8`mrb$cRtvd zD9zSmqRDJt1|(I{E_!zrc5KzlF+$IQbrn4&(U4J(vqJu%^pAofm!fG@hcJf6i9#4l zvrLMkk%v$e)8hFhUK}{oolqJ3^T7V=kxe!-UyQlJU-l2dqlnT0n=~2izo-s-AH8Tr zd#w+E)kF*E+r*Rw70P_NK6b~I#bHzS%zOn7&s|3S_rsit9=Mm(&|ykcopD62q%KQ; zUrRd1<9zQIn#`|h?e``F{P0G9ekk`x6;b6Tgl)x;lg8%EiBaes3Y^U-|bW8Ef_MEQGUpXJn1@Sqa-Vf-tllkrUby%$z z*ksm4Z8rP+FVnG6{y2N;g6;}kAMDOw!p#ktNT6h7B8OAWAVD@s70-{*(qS`JwehbX zQiO$%Ti(Yo&^6!dFY-4sau5@h6*!VaVNP{H z%6q3$eMI0Pl)2a-$nEA>vk>5S$SeAL>A9{0zRwK@-U)m#J{!Z%o|L*c7%)LJ)F$fO zkcb104U?hM^!vtJE|-W-s7<3ParC2w`Bw*#4&k6g{8K$Qp}{Dr(Q?&6#}2s$(pICA0P17-U+dvLf+$rj%PQyJ6f|CKT5;2 zOF&zj_5Hp~t-&9s(9e_pe(kUQ-2kHb^jMyLBYt5TN@)~qxs*<$dWkalI?xhj(kzp* zXjH#wv-zdE7w!C6TuhV&-O%I~2KQ;rQ? zra!v`a<`m|yM3-0RUoE#_)n*|EMUY(EWXhbYGgUEte&6N;HgnCTILn8sFG<<=8@eOB^!&2G(q4`JM;;_e7%OaGtODgU?T{|L~!yH47NNHt^JJEp~1L zC0wZ)XZaJaLeDaLCeP|g64fHu$SZVI9FcM;h;k{HMs>*Zcw8yuN}6R-K8@-ZZ2`Yj z_oCf?wI+MGQ8tM8#7P>ms5fHHEyzdgM_cILpJ1qWsgD9~Y4lO|<$ED8$7 zHVwsm9{Y9-QW1^n5Ek<|Q3y+DmPw^Fs_)t|eyQ%R zJvohfua_3U>A&y4RGZOGYHqt;sT+TSFQuf0Q`L^$iF;$__7*uGyqLHbc)tWQUlT$|h@+-iUYlH3YCuU*V zr+G8~w6n9D??Q8eWhvt+)4zZj*8<{at z-&v*0R69@F)Rt2aX29RxKhZgXqHLUG^+2~8h)wnuKn}udxr0HDuEpfx;*cP z7pho%Zm44a=+v=m7l4Yi=8xK3zpG=!WP8q0t~3KCf5cyb${d{R>C&mDml2!eS+(lF zluVr6@WNG3k4z9r?^9+>r0>kj@wXSjgv!fdqRpg?)^5pgw#jVG0~JXiK66XcQPoJA z2Am^SGv|jY{iw`^OI@y9$Qhib>iRRvS7|rqeaCvl6|4A4nk$Z<{ZIn4c;>~i0!7_p zss%Ioc4Mva#SCaq-LcR1=i==h596N=v|us+CH9?lQ&T~;PMmjKi%GQ<6uG29_;_^) zHR(U_5NgpZleB5%ArzlO9e#;I$Ulc`8vI)|FfSLlN}^H1L!hnWBe&OYK|rQ5wV9}o zZoCA#hcmexbC}eYs!>XTk*bUe#b8KyoNz z*xlKlh3ilJT3q|87BYQ0YX6k<*b2G6Pj}rM*pU^q)jhlxgiJ8boRxk~z5m}yPl%{! zcB|vPV))Ie?elfnKbwmA_70i~lXG5sCyxAs?GHZCZqlh9+j}A_uK4m7>~iCw-J_pc zih-PsH?SU?)W6gHp^G&X=4W-s#P_1Koh}7kF6q&zULt+I4zxrDG|MDI8r3frBYvsw z#q!;?{(Q$qdGPks%E-JU&*808c~RuWTgW?U0yK2ki*BC1`dIh?f5Yw^=Ysu#AeRRUzcF6+&Af9*6BuM@q7 zH;7I`{@-G>MVGZ+!Bqv7NE!|5lrDg*8p`hPHD|+2oQ{vgVw~mM?%mj7AAnlE-$Xjm ziv?V8dpBxXq9nprXCsAVOhJ@OCN!!;Zp!0IAvdF0Ce@=+{bDiam+D?D&Abj!N+J)U zn(_v8UvR!U<;DZ=Koq3C-P-KzFHXJof<}FR+uQv|aEV#xTt>GrWnVdix9bJiqO{Tv z&U44E7k0jja;6$|>Xd6?T373!O>QVC3B&8Nysw?8y)y^XD+V-d@jf10>R-v~@jD4P z@%!PNJP=~4qlswoaePwjb|MhA-c@?Cfjd}#zqGgKxE#2>h#;0`Ns#M2TZ(R*f#?5~ zIA=aiR-JcR*|NkLyl*%gOxK3T?;xH3SjFQrEe_kT!AVD&{M=WCzd6#fpAv{3y+m=@ zfd#ZLdfTm67S0u_mS5uMj;TATVq%!5 zA|u}SfRr2qV&C0|H=H{!(`GVJ@%DhfX9b!OzRF1FQAv2E4!i{J0o27os^|OQ#M2P@ zR*Yw^o_lzJev1WhQ$L7iFI{05jv%&9lT3v&mFv0VEP%Uug0yj zUL!UuXF+G{CmtX(dhz;B*ZUyVr(Y}XWS!=C?Wibn&QWs(Dp>KAQ8eyQ$7tI|Pf zu@-p^zfK4%!B53`72jusgJ5V6)nmNBu3*MpJb&QvP{T9vW5-lOcJhhm;n0z3 zkVGe+zSRZLig(JzSTL?e-K9^;8jjv~EP>XY=tp{6F$exQJJ&+kj}LJr)R-AAEo)^* zN)>8wf8sf!7Q;WR>9MaytkqvBY$5~mUYk~jAt+Ccj)NkyGJI{-#cap$_RQIl{&6$nkjnmhjX)sba z_}I=fVq76`sIv{@WucXaNgMG!cuTY86i1H9F1BUvdZAsm>y$$EYip{b($TOH1zj#R zrcu2_j(iN~lr3MpOgk|v{NNNYlz|KHxaeAN*TmR{5 zmF+b*c+rV#Ti(}6Wepe?-ZpnC!8PleRrnVN!L?VPb*f%yFlzaTu(+qKP?2S0VUj@2 zTCn5fYS1XQ0!4nt?k)NRT`DU5T?M-DTg?8PR0hLxcG|?hEQK?{7dNzZ(PudXRgW7L z4cydWW-~2j=yGTE8F>*1M%j8_jCCSN;eRzyL(MJS&5W7<`>)HB3sa!2LF2gCo?lTh zi0ZHns(_#lEgPj!=2&&UjQV{;ZvG; ztRwLqH)*MXL(70oiKHI00UM6$jO zzf|`v{AkaQ$9b^doM+$N;qIT7Np@czp=cK^)Mm?rgW>Im0N_B_*~ee2 zc%oG)2tvVgNF{6BO`&Y+wSjzJ^_4Ez@!V5lit;dvf4_7a#$EuEeJ=oKR za7VmSvn#ny3y_Cvm(Ba{^!1O}dqVuR?OTbk-6H$n$RT7JBM>sVIQvBQsBa@;VcgaA zM1_roz4t>NXnxI5(d?r1*@l;{mC0=)@!BZ8-R*MI@Xi!Mx(ufyW8W>*+#Zh!+2JHg z%!7f;e&5Y}W59%f=&B_f%Ttf}|KM2llG0=!OSWy*{0BN$%~|!0pzZ%@epGScCRnhk z^tstK#+<2o3(G)iOF@)N?PyeoygiRAg}ejJGN~ht>ZhU;zf|{BxE;eBLs=Pc?wbPA zA$aGV!zK?rvBS`U3ye$!hV~~Ob<`5l4cL9a@h=TYQO>ou06K{@BuByZKhvBzd zLjwxw-L2_j+u8PDa$N~22c+HvsUSf&mjk=6eGf~ipCB+G9g8^czLg2^!Gqo`p%+1HM#zTtWb@nUW2jZ$X{x?Ji)qk4(D@^zpk>PE9n>Q1Bjsp!Ek)jbtD zq*Gm8UkTzandJse1cmQ014^CY~*{WLI}}z+-K7|4NK;c4AP)JyS$7yQIYUeTQ~L#{hfoKP3I;-N8aT& zKNFJ(gniLH36hB&*8X4;)W4xMzyH3As7hnzqS&<4bq#@x%?0NIeeGXFMe{A$959?X zF2~+;qHzv-pzFJROCwRp9eU}r8f9VJzKFZ9gTE11}yt>K20NWr?gr|H-i@*Nt*uV_$SUV|=m@io5vr>pPfoR6Qw(a;X=M>X7&5 zaix&=p;;#NrICkRG~xH-muP#V92eZ`(D;gcBwsG^ukVX~d~ULNAaD+S>+mqV?t1jZ ze=+%>D&h5ms;OD;Q_6DdNZw?GmE(W2^rPa<|NA@g{$xTa;s5*R)cw69+#kxdmYRyX z|5}ZQ{dsr!D3t^>C>@{oyyG)>=qbJY7L@qF*kT`qQY6lQRzn#c2S|G6`tpaT3ktj9;Q# zvk9Fq^Pf2nT=+5q-R}M}3(PHmq@Is6AKC^2|8gWeF2og1cXPjX>I1Vy>*GzNe#$m; zDk4yDRjUecO8NCO1}FT06N1jIR|2n`1&^nyE8!RyDm{@aBNuH|g5Q&SjcSM|5MEU_ z`dGys44JXi;%LWsy|PXgLL z@#C)Nx*7-_)%5h}=ps}xz5OxXoLfiKz=!_rSA2L@1j277zz3y)6fC(kh(>j22lMz* zXot`&lZMjBLu8g@ivwZw~>E*@KkBwLd(<;w{)v0hbiL-AD zuv2Z*>XM%H`Tkj1+b5v}T3$T=@>hNq_NJQbvy~)Ua3O(5c^0M}NK$D) z3J5ZS#!lE26fpGij{yUO*WH5+#EyU(i!S-@dZ4nxE6VBdjP8}77VsuHIbDB9@<$cWy>kM~NLFI8!kLh*c8veUbU-0qDRB-AJ zCTwiD=b`0A*=V8iA%QNDrvGf0B*5~gnH&6VGQfG{i=_oW4#MSqH{GY?DB;YaEc;Q{C7erC4J_meT)QIhqEbfW8fr=%?lAn}&`MDJ3S8Hu#WV1!B)*lS z3Cp3)&+B`>Smudy7c`|r+yyo&dz3ZnYxZhpF5&)oLLm4<1@O$~Fi!(UM-Thn#g%w% ziS4b(qgf!0rQph?aWtxzVmyyGEyV$o#s4<2AWnY>yse$;yU0CR)kih-1U}dt*r>tOA{_l$pRx)u5WNI*Vx%EK)tr>Z+#)KD>xT(uZcM>Jyhro&k|S1W0-uMLXLPvn}`vE;vF zL=-oDsRYp;JVBQU4f^XM>`?Fbq~zdAG;KhRLk8x|i}b-&IVBjd$7#a7{_*I+4R^Ze zv8ZM%dRI(W!Xh^`JV2TXX%YojE={IUy%bY;ylE+>(kzpv(WqVuZejVqOL6mf{(=%o z34=+Uad^jVNEr6jcI3&*`;0hP^8vWl(LI~l5a z&TwzsEECCY#JM@i%}wh+zxUu}5;P-eQ-hoP!7lU3yN-@INPI8)I;5yj?5lIZDEh_K zdoL2fqQTpb#hne<{c8)z(CR!;rw3We&{K7V^zVOknHezw7GEty6)?X_g~LlQXG-FJ zJor^mmyJooB2R+oP<_XZ7o9^W6W-?M#Nmg^u;!W{xS%yUg})M2C?ZhOk*Axvss)V=tH5B#}gzG0+| z5v#W8MvRAi5M7x`*)VUH^p^mW-&|2NdHSX#Hk zsf1n|8{5>nnzHv>-KSpKSqq*YNB2H5FcRi*+`qmhTe>=Z>8_iB(2D?nL;Gm4l?NwJ ze3hxgN_C6gUHevq|Ml-kmDDQL4tkv_BPB!Ddyi`5xbHRbF(A#O;L4@>G^&?k0gpE= z#X_28(jpqwOR<<=s#}Ui_uGEy=aviN75^v@_%VWdsnFHCc<7+3+H6I2N6oT(A3<>H z3VeX)(|e3DVk#q`Bovam3z_)em)hjK8bemtsOr5YnN$2|LXc3%!%MN$Ic#mo7Djz< z>NT_+OsDWHWD1JCDFGLF+w3|@*ZS@@r2G6Ii1VFgG66(Q$NMQsT$I;ztdwPB^UM;Ha6i*seBN24)=7~@@A{IY)zpm>4u?gn_eoz zH`Jf7aXA?d-P0z0Qz4wVJmht}RW-h@8h6f~3dg+PcacT`aPsK<)DIvUE~yWP*JJ`{ zY0Jy@XTLfTjX#Hmvv%=kSXH!= zP@cswh5wf4hK&2;{JETpR_*${YH~XldvKS5kA_T(^H`T;!0x~Y)6Rq_(SQ7>MSAUG zVBLJl{qtXvr{s@#)&E)qRPYVV;3W7+jHmp;9v}>3CKMt`IY)-xjDgME@?(+`cs97| zU@+*7-l*A`l>F;!&twE=Bx2sl;{+hAfliO|FK5T;ux(36?!Iz`Fyve=f2#`0-95XP zp4Mg7^swF-P=sSG$RYH4B{rO&>QuJEoHaPSr2C=d0vy(?2h}waaQQ_YPcIK$#(V4T zI21v*CPX`KUW2)G<9B3!Sz{Cua^M(q8$ znD6PssXHtnK>5l#AMbcw{|*(&P_mGdvdPnj9B9HV@8c8W3-V`KSu(@k_EQe74+5QdQXA2ckXG>0 zVog>t^Z(-BDa8d&WB&r-0&ME*TeEdyZEoYw?jOPeWJa%l^V>d^kj<42+0O0!JbMk5cc$Qo?tmw47-CZ`ECrz5P` zZ@h#>CJ&R!Zb`^}&iw9cvQ2_Ad4VyxpgK2q=eKAh#hYS7F@ zo)AcejK8;8FrhQNn2OI_dyqUri6eLh?Y;#Q+lY^8P-Yu`V)O7^?A(*X)gNF~ZsVJk z(tEfY$hVf4u-){e)%|4VOth3WUW7ticSou6&FpB#2A|1Iv%XV+YlNLBefucwx23y{ z>Kjb-PE7%$%|J><`#Ct9u_W1rfWU(VAoye{X?`vE2T`w7J57Zzxv)N_#z^1g{haI@rXrzj@;>fOgX9$?&3ToedI2W=MwcGD zpRv7pxM$NLYwZ?^Fri1i`DfztAPzZR^BdIT&80de=Ap15=ox~d%5ALAL|1cJ zD2(@6Ns#Cqm6sf`Ufn;>9PHIdgNcUodO0}$%;5Jw+7fhQ9Vcz4fuM!H(@>aBnw8n) zu@zI&3Ab0nO{(1G%3_S(M|ttLR!qq7rxk<1gYEc)CS1}KZplhDVqCjlt%PY@HhIB7 z!M`Wx8rZ3LT>UiI&`4#`s2K$F*++{Zm-f@BUb+K(eQ4(Wt(~5A#d>7I!2& zu)yoNaIZPXbzVT#Bl!!s2>0wyYPaSo-@TuJHz8Ije6YXboNcFL1a%^_fbISXC~~2< z=$4!Drc=;76ZWMk2Wo|b*sV3xR0uZ$doq0?oZAzleApF}@zbiazCzVu5(~GM#DMK> zvWg_cZcKy{gb$ZNcwHhg4aKgmbtpEj$enueFgdYtF+sDOIo)obd*MQoR0BeP@Im_ER#;q z$m1m7E+_dVij#ua_UpF1Rf5ow#2m$e9j^9m;r|5r+0m#purk)Yt>XM*6DcB?M-+a zK6_H~FQyTUMZ99+IW9XaFFup*!6^G}QTqQ`>oZR@msU!agRck5AtMm1_C}zMewT zK5}88KEKA5_Jq$4y-_{lOa}f_H1eGE#)$oR(C>k6Nj8R0;;upjVK?bt+PD)n2=j*d zMOq4hj&lfsW3gVlw>6H3IWqqkow#`LA+_t>+36~@Q(hAAs!RN$P&~1%(4o~8G9xC3 zj~KlyaCHp)YBiz|&A1kFdMcjoy18T18`|tQDF6;ns)93xH83u+7==^qZgV#8Pwq6$ znnKVHIoHkWRysbSmt5~y0~SGcl@pB8kZVJa9*K|JMYQc;m&4{*_f~k*MChxl* z51(ET%Pv@xef~A<(bWPS1wDs{jiw|(#-9tn3bvauVG`n?uh4q_o>&g)Al;(i%B9;h zs+Yo@$D5Yo4$U&jgGTi&dzWA0w=8F8G$5AqF78TjSi)D+yvLMz1fUCN;7!+IBPbd9 z;ocMUu0ObM$fNrx-e%Q?jK9oV991IUf3`1*hX~B-zAVO&`F$NjhixVBO@m@qDyq~n zW^8*%PCo-KPQQLZhlyT|?a$G+)bX;}H5UvZ4CsNS2JGa<^$w;t3sE~@A90Lw@R7ya zK7rJzt=dZ$fsnoB#-K=%^|NFn?hd|}7E=Th?^F`H^)Cv}Nw*Z>hmwjVO6(*IT%tiV zsU1xM!JvJE&~AOBq@%~O(Z3h3g;6g(Kluiw!fL;U>9+#XA>8F<@a#|eOo%~zLNJL- z!)%ps)|HxOO38CpxZM$SgkV)4kJhLlAVs?DsCIg3^;J!Uu(~`Mk8XBn&mAE!tH4<= zFT}pve;?7YQX=;U{gY3gmE0C^x(55UugJRZ5Css<^!B`2ddz(QKl_&eHzB*^SccR% zn%7{z6(d1h*w^<(&`PvJC)ByBTd_^NyK0CgOW>KwbfQ~%(_+Y_M>MLJ?lE5GbznduD?xfniy@c%XjCtqKVKhOx&WGG(lZ+Q z(utFEAiqRQ$4}0@&bG@u4}^aIZ6NXv&-eOwq0;G}d+#c3#>tJHU%f#c5wnHZpL;+| zGhq8i9Ber2R4Gy?1|c6}P>5Ko^;u9GQt-^BJ&jw(CVtnm(Blw?%o#K8Uwq|w72I0M z6HsNS+cns~^UWj>g3Yijq!wU1-+pNjGBYH#4?)Klneism-zebsqp1$hDa=ii;9m?U zTx9`=AgicuirabJ z##NVvt{f4)$WTksVeuoA`=;MjQjo*u%zNyjZ5@dnFhNjEVo`cd!IDcaXjF&xC66D4 z_7%-C={1e&(7xf9>Ow0>yKQpoa`w+r=H`C#F>Sgn8_{dp zj}4T(h~`4dX<@)?n5Mq*Almmck70!?ef-@BS37n1Nx2@K{E|oh?s+Z^g!uIWNgv!l z*atOC1g=&aFasmbLQ;(3KCMV zXS#s2r7 zBZLg=yzu~msn1=qtLq$lggm3TWR-WFX>k%9{E8Q+J&K07x8oZ!$W@_;bdzAx9Bv03 z1N>pXwLTL>m$hCJzy@N8Ps+!Nkab_9w?9DX1N{=Y^pQsO@BYMp9R2RkG|QwfG^$^; zU-_lF7wx_-&$c}uug`?4zUd=$h`Mm`4AF`j{;XPQc3OuC4e&KDe7M1cykTD8!WEpU zXTW+v&F|Yag?O{J@vBW4St>HdAsu$go)?0O(Bj^?)f;DMDQvin(kM&DrMZo=3vnIE zr#!x9pzsu&pu<$f+Dv$3bHs3ESv7pR;kkNxVLi50Fo|B_Bb@fUY7z0g4SugCx24&Yr#*P z#9O&YW0r49xzT}|ib~!$z9<~$T01MdM(VL9zdFTLJc&e6OMa=R5Pj$uN}>6!j%hMC z+Ah8r-fP+#EybTjv-I3N^qByu_R>_S?9}xXybMwK#;UI0H5v;fbrUVtZ)isbgR_n4{uEzjb1I}w@ zl#L?!>l%+_zt!&k&X~F)Z#M9oHh)MIx^ef)2Qko$cM6xhz)b>lo)11{V;6X6>9L(c ziFF8y390ElI}gPR(?E}P{ZxeWOLb30Ko{d2w^Rzz z{9c*|Q@(+?+XrOf?+NZ%!<^9np3{Vs)a6Y#8 zo04b?c5J+6m+J$+;aJPjUt^Y$viAgqa zL%+T~vw{2s4==h_QRhb-2Vt2H9%%VTf}k2TJ&Ih+=Y3WQsz#FxY|OE*Hih`;xqEF& zg)Vd87KBa8AnCqu-qHBC*q#4tk`a5Aaj{XN;~z*d;vQN>D4xDZeHL|5b6MM{Qj~9v zJ)wIb0!|EBuH(>Kn_YDdHZqy<2RZIF;AI*3a5?r;OLmOo($Ay=?@M18RtU?;LGR}c zD<=GjZ;;zWW3oEkSq)#;hJ2AZSug_|4Mq@F&~6)Md?{2 zE)SVtzeGW$zp2Aq*Yd6#1(wbX-nPi#GzgWxmjM$5y~*kDpJqvj*8}oy8Ak3l$rfzS zmN_`TcLB6>>hr{a$p zH>5Eav+$vu`!1!!k8cae%7U_^gnsn<5QX!9XzuMx#ow!Ghqt@Rexa?{+GfGV8RWX_ zV!*D4a+6bEI53y0%`bCNbSFgVspv^A1sHJq$9fm5P~dW<>us#r4kD$?_oV?h6PVWT z8T_l9Y}3Hbnz@p7v*VUR;K0)ISJgn4WuZ$6QI!1h{ENk}wAiCp$+V}8MnUxMX~^2Y zOzPWF8ltX1v@hj-s>Gp!$pI>x12Jo04@~FG>q}LV(T9?(R&e4idVMAH0x5L#U6T* z>$q3i=s zPx;_tr?~MaT=iJn^pdk5mlNSHl&e`EWBdW{sIL$7SlLmR*LU?}@n?Kq{ihyzD0o!C zM9T*Iq@dO1b@n-hFu2dssB;9SB8`pLUzgu-HtX*Bh@%BCuKfYlX@M1M=VcWmdu4Sa&(p0fg5BGrt;SA!JaZJJ4Rjpoum4cA(0}t4ezMU7PTUN z(hy6PFN$;8Dr2_sp=|5jY1P;!VDH_Q;bnLu{e^Ph zhD5NPNJ!nIO6X5#hUnqn;TaVwC$CN*P1K(6KW``F9-R-Tx)LpEwjc=q1C<2*HyU5e z)SA0A7kNLoAALh4*`?AxlPF6hRKYp;pcQ>Qn&>kDAa40d1!gv`fy9)I*%Qez+wFL} zpunrepn3rotHAz`bs2Bz|2sufA(SCQBY|wn8boHAu#?{0;;8^v8UGsnAyrEe=|)Za zPfeM~_2=p-PAeRD)@-XpDT9J5mojNoFGUuQH!Vdr%`z#6M!po{JeA8Y(Ngg96c;!) z`Hy-;=w?};><`q}O4yH^kmWd}Tk;kic0$QTT@Qd;BJUpZtnTP&z_^Si)A=p%;?cU2 z5*60rm;tjt)5hU)tnvO3+W6^eG+f1tdQKNEK87|D}Lq(dwx?9zBB2@2pdP zmyx*^a7d3n<7OUjNZH7H4X0j;!HNQ{Z1+MvrokIMTcx6kJzWa|!5Xg1hDOc$TAoNY z+}t)bfGA4V-qFUZ&6vn}I~XY9*7_wKjH3!$9tBG-DQQ%PHlN3jLR&zyOe&;N{emju zm+D?n4}w{TtHd52E4=PL!(HXgsEwUVRO#RtZRVUi%iGE|2)>(e3)vT##*KW^^qI~S zu1-?oDl?LME=tD=vz>_=f3?{Zj;?E4fHlv~HgZVNQ2dkW|K}WN!3J`JIHJsw7QfeI zlXSG=Ro6;Q7SV-z&~7Th^)oeoZmxz!$}TxqubMG&YHuG3XAKBI{i0muHwW5vooQ(p zRPX1#lVR21VR=0`i;k);6;B=8Hp_=*g29uBD{n0(ubUF+!5ss&K7wj_*&>@56~L*PQ;huCBE`H#2c4xGXGXV1fXQ|A?-BWDw7 zTCtDRy|u$xLoxVx#H;&-)M#Eziy@avXjCstwv&&-0HMO10j>rj4R+S&foieTxOw+)tcH)N_K85MyDGo^7X3w}@GE#BnTIF`R+ zh?IWtb4FHq5{4e&Aelta9Lsyz?m&6Tqk;1(a^U+F^3A0z#P1EpbB?fVT#?{wbahW6 zil>0mwE~htUC$&A=zVx<1&GhCwh2qKuRMF=Qw<0%=k4zxmO_&>Td%x1(ikoU(9r(Br6BKgzg~3rIR5=)v~w1H?i@H7 zB@*((%b$H~tXbD;dkY8jAAQcTsbGGjY#~W252j&s51qsyrah#nN6MU>`8jQD>$4cJFt#lqm znX#WIC&!iR7D1I){_n91cLT|Lj%a_@V{16={d+MCecoed!jw{ooqXqA~r3lil@7@6(-{Q zUw~W)>Ci`)_2u1#O^g*>@X`Bl7zhr`PYaeG)vr&aLm`&M#*v?14%#)l_c(2p729m? zG%20xy1t!U+%@|rf}}~ol}lPQs+U5W$D5Wyhh~|iOQZUQrpGVUz0j^Y?piUCz`QR; zZH|ZvE&{61kaI%}Zu(gI!DA)uOn7hG?1$LN;acIOhWC(X16@8FurOaP9gN0nd(OD_ z%r6AaQaZ4ayiLg-Xoqn%IwhVqz4NUcTD_kw^WR~~a`wLMxw+^shH0lZoiv~r%%RqK zi?aUS` z^z>|s>)OBMX((LLr>2p4I;jAwTr#6ky%hC$ylE-SX_iSAH1eeoJ1$FpiI#%zxXQ|T(*)t7 z>89!MJOHLd2Y#rE-VI@MqICDI(q`m~J9n|x3n($>-5sydKbCJr^qFvxFeKrI_&2E- z-n+~5iRUet$f><3z(BS)@zbRgoc(QapX>8%Sqr0%Sv3T4AMkif=9z=ZxQWYvGn0X; zKKIU0o#$}m{Dc-i;y{&iO9b9(@*gGhBP{v-mI!D2bB1UFsL19|R4=A0dDw+l_!9Kk z=awh?o6iV_CeP*#n9(%@N_M(8>itTG1*SA_o6)oo#kb>TB?`^yHXZgho@nu)p$H_P z$fd}Sm^;wmv+K9f$Fs5 zmngIb6x!BZ4IVc9n*&^jZ@&I1IP*j1Rs;fXNvN&OIvR5O@}oFcoE7~+&(az`HYC$? z;_7Z+orij+#MW4?1ge|4ON%U+EKFx&=|CllhJj&9_`*%$D{a_>`fbl2pemy9m)@m> zFh;l}C&E~7Xum&fdKLtm1wOiDk%&Q~xa(XzaIkmhqTjHNW08vT;UfQ(6vx4ay}hVf zmW5$+hy3fh@f(V2E-}lvaZWW#N&3a#gpFxz4oN0snq-j zb(sBR2e-pD8MvcGOJ|c48cedMpva{LG^#`Bz~e+AY)G?AYDA;@MbVgFs(Vq)EE=a9 zT}UeS_MCX{hnvK*ox(Q!Nz2n_Fm09(W}n2+aX(V@20w>xDT62U{AUr#-q`}0w_AF^B!xf+*LcwRH- zM&DH_FpjL2Cs-%psHgvtC3X^A%;E=zJt}40?QGC>Kif}ptUep9%{%#A(ciCz$BTlQrb9eb*hEUML7GIZ46}JdD9_BZpo%S z{YeDLhxno2lsC(gwOErK=k5+%V9k;XkKSv%;t7hD`|4y!{q;u9-%PPP>(t>Xb((`}XtpM4{N!HW@vVJT^}=3hu#Vd)Y}IGNkN#~T9=gLp zFgZY8HH$weiM;)=;mpPVv=tpMm61sKC77or=-h2?&Ia&|$deKjEUw?axPfQ}Km9dX z^1qjQJ60rv3;&*O%g2Xj=B%)9n+)grUfdj}9Sy?PuwMoyQEU41Pz^Ti+`Wf)SEWKD zZQqEOPz0)kN)8Am;D96)@76AbAne>C4H?fX^jKDmLzi71+GbNNj209{GD&?&K>y1+ ziq{j{kBQ4lsk6f~Vq>V@-ExqzVy8Uw>oRgC74jTqJ$6aq{qnzdSr}^Bo|Gzi_-lWw zBQ0B6Gi}bGzeA`|evV4|4N^-Au3TzGqk1V?^LW!zw4qrhwWX0Sg%Fjt9rCaizC{1pf(~=6!_r| z4juKkdafbAaPDtFSpXF`v=D`Py?+&MAP;;0NNXlqSoakoc_!JG8KvWuU(f#)`QS)%H&WcSr(Z6YI?$+&Q%4>Qic=?= zWm0Dvd7K21qzk`9af%c+SM8~5T+s2KfU}!|6y07pbrd0v2Wv7CW$%3Jau^2*sKHj? zu6~zZn6Npb;OB=wDdU)X!3VOViGN23Chi?jD9b>BU>mB$yWy25Hj!5l=eI}4B}1X+ zo}YKd<{(cpA1SZLR$u%!EUn59g{-Y>J}lwRGtW$!ke7uP;)idc&u{;!!sr`>gij;V zJvm~WYE+DFYx!$A9p9!l>~GbNs_Lg-ZD_MEA8L6y?gXU(y3#L{OWkNx$D=#{(G-s! zG|Qx(G^!urUi?zsBP_gq+T=s3H<$0u+ zg3-V6AIqL-FqIBD71r@2#@JxkB(|Wrb%5hk;Ii@NB)YZVWOuG6HW70FtMQw-L!VtA z!n=p^A+Mi9!7r0=+{A0Q_4Qdd+H6fXN=B81b_NP!dJkPDth9BCupr>~s;BGA;nI<2 zcU;nRS(1ynx6PVpATN{sc{?dAw;U`q3)QIOX+6L*N4+Pw=0yDSHiE9OAE(@`^ue#V2?$WG72{ZTr07vV8ZCMaxACd|FPb=CZ1x@=a~vHbxJ^I&A@2BU(f;lR%lS9NvR(eq15 zBo_uk4b(wH(TuAMTalp)2M!U;pUbtoyOd+P&3NbE>l!c>6ebFF`km`F!`s; z;6J`yfQtr`Ynz-s14}q}+|-b9o1U<~X_(ZTDtq@c@uQ+28Sq*%;gvWx73s509Q6zo z^}^4#D&(U68On}@Kf4G1zBeTnquX(^ za{|^Hl4CH_XQGQGHyM@7@olKro5Ew3eub#w4V{D?`mkQ3`VB4+t$LA2Zhu2} zSd!l*CtSVthM(itH{8-w3H^9gC?y5G<9G|92h0^Q%>m&$LLXigyz&+9I&_-=8GlfC z7ugonK`y_ciTAjl6K2<-%4Plw@YYY;&S}8G@$v!H7}S<`&bZ*Ol?2~DtHkqT#?g6* zdXP5po3<8Jjr(cXTBs&!7KxrfeMJ|n^Q60k3JWdqCF*M7hu{5f0sd4-gzojDT z8?Nc=8ji#aIo`fgKC0Lc>oLq{RL14T1$b^+n%(Kw*P*2Cr^V}Q7oe&|Ap*xcQ9J6Q z0w~mZH_m6}1mI^?jqQ1`oe&m#x28gLVbu!veq(`}=R7(D<-oxFMc?OqOU8yfzU%As zTZFv4IR1>LqI2mNf{eXGZj;x4?<%~}Z)#iTNohEp=##CcXtGLYpK%xRRsiq$n59z# zn{8w}WjRG4Y%KdOhatrL>6KFf#wO&E_Phj(*1A@|2wj8Hc>0BMX#$PvcueF!n&Kg& zStd=Qk;g+oc_#BqbW3r9{>}}&jFo{;xA=4!Vy49Hp+&@?{DvW!$tE8f6{0|eYecCC z_y2wYNB@PE9);P-*B|9_YB2Y^Oa6$yrqO%)wJd<+Z}#?B6cPqpo^etu5IvT2hD<7q zESMd4NLz{2Gj4FImoJp|^Xfd>DG!Z@jQ@J)V?GMgxw}Qk4?vGQ7b8z{T}ZCN*Y-W< zwR8+sQO=@@Q00{N7p}b~92srynh!N)@XYwf#s(`2(TAX7uYYC3IMW_CHuNq=OO7XV z%ZIzJcHfFd7eO^GaHA0wn6$aUs!@|I&zNjcSy)uHD~4i2)yI4 zV4}r6O4yOiKPLoLmL1!baA|e_$c4RAupw;*N9A0odod}&>E4W+F-=oo4FN=Y+aJKs znEH9!c#!yYstXhPiGWunCMc2&=v_1+)5-aU2Kx@W%bba!p} zeZ0e>C>dP|KfPU1fpcm+O=@{I6qs&y_h|?vB|f!+v-2Y2ks3r z`tGEx6YuVJOP39vUGh;(aRr*OYoj(P8$BXlyt-1SD3*ufJif(_HE6Pa#VV?zFYebo zH*s+)+N^h-u%~?-9OJb3=oEb95jJc6=upiAAPKrYov`~$d;l0EtgGBTBM-$F@q`o@ zoF0gc&!~YsMEbRKiG#;9Ei^N+am%7_VXx(S8wu0GG1RATM0mLC#q>{#61 z{swjQr~_>uw;=OJX%2~~l;%>@N0LKRx z52J?-6IJuk5AkbB`z_&Ji$FmCPnfab)AQ>dYpcd4zl*kY_e%!qJV`B~Mg`u8^{6@p z8t_zcQW%O!xEty4s_&z2OZR<1qOoQ?AcBan+|^5mY6It1OP{LoRqv1^7E&lST{Cgu z$j;OtK)I~9XP<%?$H`O6_>9qpOkn2EmVmG(d|wLFwym3RZM_u}qtgefKuobG1cDwn zOG>t30uSl&9rN~5rvnwa+J4i0>;BGx2@3iLPn2uopkC>oU8ysv>CK)RVQ3_}+k?Nm9;N=TIWAxB z;doO5lZ9`65Nx4>>HBl#FmyOEB%2jrW0(1pzg(<@0r6hzR?)L2kcdiYB}J{|tGKQt`D$7f(i)0d$tAv}k=!x9DW?PhrSG5o z>hK(G$BrN{@o^MQWLNZA6SuSl4&S};o`Q2PuA>|IJ&ibN%EoQ{llA6h4g`f3bori~ z4yol-z45SS+z<2c`v(|#3OTmRrD*);Sc90qVd&h>dH+)hqqf?yO74CV%PXnq*b??2vy9~G=*5wNUu$jO2#v8DM+<^FUDGDpVg{44rgYk#ltXOAafj%2e-ac&q zvAW|!L0U^9DkUXFt>o*tt|a+-S{2d;idxAx@-2(sQG)4-R4&zVh47}+{B(g?CKyv`vMw zZPh7X$W**7vG0!)NnT`=))vctZ#$DQ`&Imgzim>0Tn)(MPmN23Y$V!W)A0)_D$ z!1@52oEq2bQ;-^8Q)YDS$~c)W^*I$d!*eiAdf+H>iix zWFd6E7IVb?H@P6*Z~ww3DHEp%z)cSBYUtDF%#(`cv0@JTWrA(e^Lla9G&gX2JV~{7X6re67z6 z`a!QrPB9!tUoKQ1s93qV0{%vHei3?uh$sYW@03uC#xp)%fBC8$10U)R8<1_yrt>GE zTMSlTMjJ3gCrUYf%7aTZ#^1e36&_Nz+x^wTMs6~r*Z;-)Hpf4ICq}{! z5>YAbq^OmA7uS^}-%YDR+Cxz*`Ch)Ik-V+P{?C2RX2Whiu{d=%P7LY#CGo~*ij%;W>$(V69D)XMRLLa zuMc?N>^`_|r5O{iXWwd>>8%@w!hawt5`NqN4%hpbjfZQ-XZJISMKPy#F%1?!{ySuN zaw_oQ2DhYa94HJBB7jdmhKvb@L1qi$Gc!}clUt!u3GyGO1GY*;0bV>pr)!R<+dz`% zPwv#wFBkcoeB#JVAcS&xw_HOeS}(W!gPEq?Rzq*WmuqR1r|r)5{ZMUwN=a#tz( zbUyVGtfc(@NwP(aro4Zey)RxaCPReUx|)%*dvJHgh)&U&pCO?Ex#n~>VX-92By?}ad%+WE|z1%l(3hmThHo3ZF^$BsK4&}Cf$Qg0r4 z^&I#8T9#!up%N~Q`=&F-g6oPfY5m)Bb`UeUCT{)cLsFr3ob1Lu5PJn_Zx zJdjPFgMlaCAKkwNm7VsKsAwv6?_BbjN+<${M#EDk&#VApZ(Glv0dM=WNC9Vc$4`U8 z1dIlgLi$?2V+4kJijolOnsL1A~@cKP&|D7dcuEtsJ=^W8z$ zdL=a;b@xQw6%VH@BJ<}hHhMHWor=!llPn)clXuIgN*fk$);#F1Z!+Xhqpm&j6Zh$! z7kn(?vuuYu;xjnESlFsvsTZE0@6PMNH1s1JiJfo+AAmcS4gIJ6y8X)ElDQE_WfKi} z_nbyZ>E6{?H{IgMyDSGLxUjYbSp3o;?y7=pq)yT=R7!3XwR)W5KbrJ7O{+pWLy_wt zi2G;x7Tr>uxG!ex!~ekonwDDU{{@B$&R#e)YIEjgqj$#am0zD3(VhpP6yckOlq z-9NyY$|U$`I`kt{v{TFze~BR zty}xA-|om7P=lm8;x+PVgMeev*rUHueipgJb7W1OY5AC#a4%}$ihqzfo+y|H|G^sW z-|0{3MDY$P1KGHb+H}ZZ#FE~1)MiEB9pTt%>HA@Q*PF1L2AsBC1@dszPf+$xyEdGT zdM#^*jCcA%mLFTf*GEh(M%}UP=b5EAu)P=`x9;l&(a=1C5P$|nAGIt zwiC%j&^PYgey2JBM3>O>9GqrR@6NVU#w_=tbUpvcZ!9?+?lytWQ9t`M-kST{h-K$4 zr-8~CXs2jyKC6B;T9qy`%A=C3@CcoE4rMzgA6QI1vReWjSY-#A?dX*cd}uvnZh>rS zu+N@PeBZ2kQ;NlOnE$n)45gbSqEfm=Q7idvt}97?hgOAjm!kHry~npSckQOB+v8e} zrkP)nT)qXo8ndL+W7Ol7s3LuK@b#^wZ%yxj^6$b+`|tUqZo5v)PY0W_2QMzhUA&_P zL5|;+f?p<*fK!laIyu>4~$L6lJTRq4EvRc>$$Pj*aSGy-eR)km@ALR z+eVJ6zb!6Xs8i;5xP63_MU3sKEnNesZo>!r#Z-$gBADC zd_YoR(bhbKpgb%q+B;g}k6_N)*u=*VIJr1^(TNBXH@Qt$IJ-?Sqs!-4y?no`Nya(Fo46f!;!P{ z5AP+(#@uLf{BArr35sM6@G?o;Iw5eJP8-qB5~203(_!9yQgHHwcKCOG4$vHY*!gqh z;p2C1%m*^mY;PAKwKddHsmW^K>N6Xb)H-X)0Sv@FBe*Ge9tt?aGwM4uKeMAvyCN{G zs5NwV`!YE8lQ$6Buv0Zxwesm20i9M8Vx)Z@aGUJht)^@tXI6eIMnOz+r_$FG5>Y8V zrKpwMlj};7KciJ4c~R8fwa@vM=B}k=dcuQeYN+ndMSPB;wXpsLFwPG49RLPwArV2% z+!taUZcKZ^Q{$>=%7leaP#*9~)hNSXKo(rhid@WCy;@E;A9pE~{gdZ-BQK;}Y-hx3 z&m|KMA3c@mg8BFh%TK(VUKy4E8Uvm}nr1u@*q?jvw#^5@5rVb-&}*@vT3QFPHJ zrh_n1`AtW?*4aVRSH|J;sr>1mjl$Nbr4b`kXAqYvc|LSPZMCD`z@qAG2k3X54<8xjY}`Tj35$-Xw}rdO=Yu;Y+R)N%)FZh4h*t zmr#6Byy06UA^)P_rhT&s|1EN|m3Kj~&MLi7u*aw~r{m*S-CK5BgWUAvlj%F0Hm7Xn zfH5X)GaoQ7DaN1|bI2L_Es~)^@X2{5?8%R5>*_v8gCgp-{ajIqQHym}?FzDG^+&iK zAMHo|Z=428wy7r~srT@FA~1mnj6t4uf3hwec5TQ{_w+~rk%oI5iMa`tuPXC%P#?k}Q~pV)mDCvuC?9yIRUnt$W`k8_a@!~2-`-&Eg_c@r?*&bJzUj&+}} zpPqvjyfrDxka1gr)hnuTc&{$!Y9FCz{g))m_IqsDh@*dg2Q4Ejar4cqEC`A}nnT>; z4A?gAgSox}mR+#$EvQ+9hd;J;trup`gMi#MsZTtLCa&@lqyf#CdybZD2zh%I zmPJE!mi6h#^$Eb=YKgD(*!>O!QEgcPS2mKhWJE45{FgN6dPpeD3h)Rmy!ILBub!q@ z;;?(i;}suNG@nHFk}M??n%1^{5tf|`0ev!SwHaIhA~(D+6C)<_)}v%t&Vk|umTY|e zAu&Bp2jV8~q)9$%yJ#%O%UFlesUl{=@i_c&@75@`R!{Z%>!N2n$d6418kb-n&IIw# zm!Xi8?U?W)?VljK%&AuD6ZUTI zjecfVcB|p+)ZA_ngt7a3%R#xt)P#weZ+;?a78|kIb$QBz%*i;!5U5=QiJd)w8{m-t<#8xA-mC>QY`b!aZo>X zNij{JOTHwcQu3pymHZpml_dA4RUv(+sD1YT;9Ht!|JdtqUYAVE0lClaIYtd$Nr-+= z2&G<{C@RmVq~&e7SoI1yi?x}}W9UoYUaS6CuwkP3`3C>#dG$Y;^8{8?6*2j*YX zwO33M;(>b~$vD&5buC!HICxY<7Np^|N6oqj$cxMn$X{UeVl*? z%4USQx1pha!VORE`WDt1toGVS%{uRg-yi<9Lp%tZqAhCV^Va4W0l1{cp8+c_ILUUMZlFG8(B#JmbK2P32#g!Jr_#!qD55 zDE%Z+l+rJXS_yx1ok+qzv?`-A4Do zm?rv6psCGX$>*2P+ko?ur}|m2c|qR9d+|ZJ*Zn4bkGQjOdSeIHRZNb)#d1&1!Mlb8 zH@$fx$%!Q%>sH$44Oxvx&{U647KSdlbvwT}9v)COGv`Gfrsv+K@u4%|JA?>0pW-2M zOov~wdn18h7Vodt%0#oF3yZ9}CBSfxf3>q;MS%Q;m@T)RA+#!_P>R}jVHn@iJZJ^#bdC}0l4qs9WtavrRlXtd+QzQeeE!kwt%OC_%%|VipEWeo)?oXM^}yp8-of z(mdPpiY?nqfb}*@(m?*=D}ZjDs7N(TpEYh~|I%#k4-mfBy>(TSdC#YnHA|ts<2j2n zU>?VeniM1NboXtl&xFU2O9j09Msvb}g=p5was2qYMr{4~F8zWkb1~J=zule=ZYUlW zJ&WKvX8^2I;<*#no4Y=IBm3tO6-J+Pv}#w#oj08d;jn^#${I1=3{Z2ZTCPP;H((sD z(&I@neEc}n9qxF`hk-N$R_Q=#`@_bJ=Z9i~!?!dKT9Fz3LFUW?w9pNr09zAfs1u{uZlkkm!WtQCVb|e7lFvgZVZ7p3l5YVpf1~-C-we$(EL!@ zcCwtYCoaO0sEe`ftytJag8$stVK#igqE84OHM_K2e?vTM>od!1`im^|nw@a`#;*|c z`P!C-xX7S}2Dh4T$%N8{JVzafh9ieRI`TCQY6?Sxr+7S}&N^jOKA2M~{+*E_J5pph zVx@x*xNWdnaL=_6hujVe=pIDVe!T2L)`xPh`gh2gN}gIhOjxfPo}+cXq(W>L@|YmB zh8Pu_Kp+}ETvJDfiCdK{3RlyZXzxM`ruwl#ZBJ?4p&ZXrOP5KpB%)G^qo|cUp6g1I zC(x>p5-DmWPvTn|$&D7&8rgPGHW1Zs>bY94(J`gZug-HmK>4JJ?x*MK%MH(;AK@gk zUDFf%Zla*x5)$G&H@lb`iQZir)-ZPdhpPiSH+P?8EAu5o!{s*lKV^9G@=MR~>*N3> zrt;rooCC4ZnrEU{KVHE|f`D}HMDi6lV`SwzVm((JzSiMX92`EtZQAp2hH&hBiUvBq zSNfcb&`?wo!Rm%L&3eaU(A$l$5@?oImy?uJ(M4cOlA-T<=jJV_3i$mU73FK?fM64? z`UW#?D^@k3A)D1=DE1c!1ji9wmG1}|cQiF%P5wLiA}_TV`LvU{YdJ=IfAKfuydB%y z+PrDiy=3IC!mkeuRZ=pEsFYGDY9&wQx{~B+v?`=@irNQJ2H(;=h~yGdCM@Fc+~Ow? zxq`zJcf;J$f)yX~_1Hvpqz##&0QuZeb>|p5#f~XiKEs#|5IuGou%iZ{#5*Cb*|S>Q!~DZP$cnC?JK?*oYW44eyGJvk;W8nE>b;Le zk+$>71oIu7)E5cG7V+m7j1I@I6;XdqU(P}U;a3|aqa$7Np!e-PMLVp#z$5Zn@w9S7 z_Kd&kB6E=Ug}6NTK%)_HA<=bcD95xNod)B8W07ALXMr6j)*Lrr4TwEwQeFvU*|^;J zWDJ?i^s**BEbfuI&p8AdZnSe`3ix%POH@d!LMo!jFOjIf6!R^*M7;ju!0GqOd>}vW zh@*XTWpl6Vv-M^DX51R{49A@wPX89ZIOmqLcX_Z0+g-)e2?-FM zcIj$D_Y4sBu^-5Kul0p_nZoPGtV=P=EtCVc0z$x~(K{is!^y7lxdF;3`3oCR{UNHqTnq07DAS_RE;~ z`&&AS#$CDvWd2@|1|<6(%7Um}@S7b=pG6a)1*I_f@kPbVOao@Q=F#EEuo^5@@ymHj zyEJ$(d%U9A>}rtyktj;3grZi$Qmzw8SVpTtDyOKGu!3)CB$RVi*|1D_G#rNw1>+Ay zS$u(Go*W*d&qN{C_8dBF;p8_D`Fo!z9OB#fAsN!w8f~zqUM%!~0*8KAny{-4h$TTy zJdJ1a27m-4{AWw&VFt|W`M9;Ocj?N{(8csoC>)X!z}GJ1YuKSI6m0=DlYz4`Uw4d# z2|OEk$QugD+N&+q-u=uu10#V9KLf^(UYG^_#K;8oz;cU<-!|-3Jq{Pnh0z=6aCC|=Ii zXP;+RnO$h*i)6i(u_pEb{Eaf3)1}yi=^G8bI`O0$EH-c&?pNH_h<{ERG2y~-J{>pG zAjO$SMId~yCz`Xb%c#rPqZ*q$o;Q7s=Rdfp$HB%)kVqe5(=Fp5Z&1ofpw1PkBXLk@ z0=*5-Zr-7Oue#iNG97BaUO8-&N0>}eS(S~&P4_3wPk0Po)Oj}}jS@)Q7&9UR-J24( zWL_G+B-UHU`!s`FzV2?X-(?_%8NbCqv6*>}ag+gjI{A}JPP#t^4O-XpM3FArx9NDN zk=-rX3f`*o%s}PVi|{4q65%jGBd>0;WS#jz$xdH2rsv7o4=xbf$CN}=N@f(blACi~ zNpcHX6_O=I?O(AK-_raQpGqQ!$^%*8PC0emJx_7<29C&l1m47lSv69h9T~|X8Q1aU zKyC`Y1B!N!UYk#u=7O|AlgED!Mc>yWUrkzAh|AsT=e1f-h(FO1QwZF~vf+VfIK(Ny zFU;7HsK3E=AL=naUilvxR(?Fl-Jj#Z#K@-aX#&IupMk^xobb_pbLG$&5cb_?QZy9E zM7&@dQ*W#+wlZw%BLjz!QC~-sjJ|@f@NfA_*eElc=k>1Tpv^hUNB~5BsohtUv5p~ zwI$GC2G2E<9W5W_nn`z|J{>|=q+uo@$!$?Bk16jqVmyzYYnA{VlfuYJ;T1N2HF@lw z`$h1uPu4!Y^*Snj^V$9?LURc8g_#Ze9|pLAXi-nbiM+4%hb8OzYsE*Ojs_}DdphB( zqmmp+M5W|JQ7gGK*Oer%L90TlNl`0#Exx6Z{J{lphZE}!*f#=2xF@S2;8{+6^ML#K zZavSs>u79ye&QhWKQG+Z1!B{2e00)?4d}DYplWpi%880lSjLOubSmEe)xO`Q0|u%# z|8dmHZyc=C?XYug8`hf)a9<~trmbEcpm*;WQ^Flb+?Z7hc!IV)K$Y&##%I&)e9CQAr zZ|;hoWH$VMM?u18T_&!K!cbN4=$2#bpL{fpEwBCjr6D`N(r%Axsk!Q){rORuR6h|t zxZP|S|2``(Da6Ls1D%%KD8Mqp(zRPbaBwPZ3>~DSnkeRzLqMucA}Xai6t$Aq<+_sO z^=MT{^(ks4Z@{-Sl8fZ}^E?nlUn4gR{dDe7!^TPAv7BriGWA(nC*E-71!C^-Ul{O` zaQ7e+rc2HMk!SLdPsW_E@d0a|a=lCgZCz-HN2dV&Hl7}~vHu6uYcu0R*SQX?rSRz1 zSBW`-9#Jv^uPPOqCg-f(`aBnw&*BXR8Nl7q%#LJ%{}*m*lY~o`4Qy|fL$&Pt!=_E2 zmjEqK4s3UGas;+kI}!%$HwNm&=$SsyXLGs}B8J9^Kg1i_-e2+-#1|C-Cc}U7gG2xZ zq*)PJ!3X;{9&pTqpwu3!6Tu$a%~WzpD-z$>IytU!b!KmaOe#lgGgNgj?l<`1)Id1o zIRXccaAH!ab^6QW`l?AA&mAuLMVaIONJOR7kfK)dMqF2tyfLi`sR>0cxv0N1BFt>m9ogrdpjS~5;^ zpGM3mUB-)itM_GL{Tl&wKOg&t^B;Asep1VjRauZtyiXC_2_|nmw?ZsjmU*|D5EI>R zraVmgL?&6h17mZkreVFDT6HX5j_^1voK)wgUJ4t%$+DmxYT|0k7gexHjRZx*(i99hZL z=Qj}sMuw*wud~g95&~LWB-eEAK(nxSIc%UZbB&ChvfII9Gp7C+vOKv-G`JLns3w{c zqN4V|WiA&>X*B4GG5dKOpWn0ARn7mNzpmyLLlsfqTrQVZ zz>){e3XbY9&rKpv>9oLm4Z%?a2>OfH|c}Kpbk$mkty?|5C%3v8UYU|Sj=u-dzgzi9Z zvKxG}##sAT;nd~){`f{gFd=+{i~6bw8{78itc_$ZvsE}Se8;qyPfdKu@SamopBa~9 zJ`rtfCSQll?%VHmwxbRT-m05M{XQkP>`j{{K;r6@No!ddh*1^2GT>j|lMlu~nIKx9)x_Px~HkCsbJcrfYgaZR7#yGY9;T&btTEW z(yEZUQRI?~({gvdMUwN=vRvlRgKS&wXZI9VT~qZcEj);Q_`}fFfXx-O;af5?-1t2H zl<0h1ZOr&dq75Ol$qeO&;o^8?cf6M@|82`Qts`p8!+czN=uCj$qIlfm@OW^%yAxB0 zueSmW>{>N&!KZlaL=JEF?xaJL48yA_^VMK`Y)WzCr{U1$&(>ejZByXgS;|(wNX2Y( z^8KUlV*!|KoC3xfrw4cYcmvlRU;fyAks94Sx$ll4yB@x7l}p1S6eYKorC3Vi9qZ`< z#ozL`TI#4o!+632;I6AzqH3`g8_ixTMi<~_G3^kKQV$YEDfOhNm9Q7ri6rbzt3v8S zQG3_+}F)dn2w2 ztwwWDE1@KSL{v(QqE_;LTvw93KdlOB07dOxJCJW_?poekUrHr~ly#hrI|+x2iK$d@ z;=G>!JEeb z6UN_KV`Hcba0o%=@-t-9-KLIgzUU=yh+i5!l@_1;2_}7wT8;W*$7EM;b>A0@bt5Jg z6yq0q1iv=6WDCg3_i8LsD$Jr{_+M+kj>6Sp97SkGM4@yOA{;FSH0{atq zl$jaxYM5>vGFO#Xf$+w8Ipv%Fp0Fk(~ziX2F68s1!w zX6kKvOt3KT8mNLBS*~10?bEe*OU5xP8Gtg>O!cQpoy~`h6ToG`CwKoFRK4T<04hBd zb(frIwJsN5a)jmxBPL#pn@m)X3;5J4k&=d!h)QV$MXlr`xvnJnC|VWLXo}jqb`0Or z+_gqK@O+!_d`;du8Ra3+3-H)#=c|V==_9Lh?u+>op5xe?C+1Jgy$3lo1<}*hm>u#i zanMQp3!%}yga7MR2!u~8+eXjY6E}N2%l{q>dO<AB>K<{je-UAgzyd$LaB?9Z^zuq%W9pIZee zj!tLE1bF8xDjmE4abyh6l_`TmoY<+ai|bApc>>6Gy_-W{%JG=~^MB`Uz}`?kB*h~K z8_eo&c)&SXGg@QFqT0|9-qd{f#OtH$N;SAvLnikEB;fY@SGzf;S7QtE@bz6+i!~R` zw6~F)Q)yL5(&>ZjP;knO%#u5+tvAn;1Baa|q8 zTW?e14cSP3HY>{qF$?_V6^h1EuR;yZ^fS1f^^de2ngW#Ho<{8o>0&xUra*0LRowiNo| ztS0rzyWUiFhkBf9Y^otUptwb1&`VD1(o;z@NJOPHlcHAgSzK3=d^W8LX%0p0U-4YN zrTHt?o1YfeuSqr#q=gJ4Zs3lO;T{V*ZpEuaJGMTmr}18#2mM4(o41UL4mM%}=t|6y zlivFkTWkx2>*(C>^j9pbmXdI+IJOo*yE3!DBZf{)cvY&BP$oSU!1sSg1}@%9B2bCz5a>tqN%oMJ}Ph4=v_f^jEAdKRxf+D=3uD zMG0$AU{Kqs(EWVd&l$(Vahp-l#j%7t z-q}tl%>EVV-h>Y6HP!GY!nJb22BJyiZ6`t%Wks6!ncQ3b5u_#b%azhnidvnPaWP1z z<+LiK6%@Hn!e(P7-y)q<3vaI7?wS~ioT|MtP>td(F=>trj~6GH&=nQihhYwT1`>=yVvj^N55O5RUxgVsQukazQupH_&VJeg5=AkQ-jz# zZ&b;i3+J!Kg-h^gF^!px;t;hj7;8@H=$rpI4`h4X!Em}H%CvY^ zE{x!(!;)l}arY!Wbo$AbY|&XT?aaX9eol$lv1bF*^7V63ygY8}s{FRHYn_&SK>fz< z%=YkN?MrFHa6 zl+t>N+TXo_|2X>H8);QYn<(<%Eto!=`4&A|IMXNmlq#i8g<76$NVu)xh`U4kHW{a) zc9Cf_#|x}|3LS+8cc2x&78}hNH-NEUtA>bbon^sovtd3TnQmsvDwlDLTN*~L${sk> zF9jDQesFqQ*N!z@-T2S!zsX2EvSlV#NqAuxcV0?@gO^(mQIe zz^xcBYuN(+a2be0;pV^K{n>wU68-a0L4j26iG_r;meDirmO_}d=ZB9&?U^8c`4=Id zjW!A(0`e9TODS!osFijb*N>##POC!NK~XF1PQIm)w!`zDJqCCZJX(2rPq@KXpfu94 zFI_#s$olX5*96>5Ta>Zxufz{mw-!(0MFA!g#{6~YQ+~NDh&;m~TC zC6i&p1drUDuq=8iO@|s2n&tRaGN5uGOFQnQmL^q&5Il6)G0?33s1kSygDHz16PUeFHd-1#$BR`6HsPS5ufUdjd>T= zWD_P>GRdJR`wJm>C)75(I#Q2``jTrNOx8cWzt+AK*s`bHh&vh9YyxL7^w4GD4VQE) zJD-3J$vH!xOpt5x!H$D^Y%zDmnLytf2S_ZXbdaJ}+CyAFlGc@0g>;xAmsY$ekMJ%2 zqP)|rH6aGF;LS4A(rb%uqRdD11E2BBe`Ma%WlM!+?^{ge71O0h(PenO(ch0#DRSA| zT~E?}p$pF|<`TAb3c>SJtdQl|SDj928 zUe`di?CLRk1;#-_Uujw33Qu^mOD{LczCc6y@RhrS(Hr-Zufy=FL~QmxV(g^>#%y*y z4$t(^RGl-x=kw${v-I?s%tr{svN?Ozi*M?v?B4x%)GW)6`3~=PsNd}0Xt1`h#AjqI zN=NCJE2U!;wK^T=VvtTJXjMojDRP~xxcRdi-=bTP+nDIj?!9M!RVH+*c`UWlwJS2+ zWBmizlvxqEpuRqHf3ai1?V%r#+Zp)=T!NpcLt^bzG=#(7^^gx@{HVq!cuw=_tw$L% z&jZDtyRWK|`{R8l2iMuUBWiYMRN1(0)hg;)9;WE)T^<^l2ON+QJTDHC_L|)xaRXn_|I8YcF8V+&8)-N|m$c5nVxt%ZMh0AkK?bUc;IG?0?h~JJK zB?7x43msiOrGWgHW)?c$s(CmYDJ{8k_Ax^?HK^I^Pfg1~6uaJJ0q>(~+>S`p>+Ol0 zR0fSkzSHl$Bn2rU{PZIgrNn^ca$?h_&A(sNm)N-EL0~c4fEm9S_OyGQWEo2wp~J|J ztFPWHGX2%#-TwYa^LHZ_+tJhD6p5>pPE*ugiZfhqx)f(=RY>P3YM%z@`IhEs@OLry zT+acEA7<1z`xf69|8Udk;*GobS?jSrTf1+=>w~!;VTV1>H+f(HnF{pyY{bYKhgyq@ zfX|EkicCit>2AtI+e^4U8!W76eqzRdwW9kSF=OI$Z zS&WEMvguGP)D_0ei6GqFZv4S{n?KrF-^)V*PA!iB`Kgc!1s>K^lL&zxl+`|IWWv_p zw=KI+mV&KDc{?q<9)*f+M7>>XB)h#03<9bm()3skh?&`@M(oCR8W_$jgl*LL^tH4C z`E<4Q1#`A2M8P8Jo3Yfoyv9Nmi3=o_QgWxLmG&ankEFdst3tX=QG1tP;ai%!T(E0? z=YWCmW4w$Uit&u*pMLr^*kJpWaKJNLTb%rR4hCg!cZ_@Mjmrp&rkHBX=Fuzc?!!;u zW=M^Amy7X7e_}K>Fkz!JyWXn|$%oYAOIFm1Oa(E?9b(ID>#eBc=~M{9Y&1C@#Dm-| z5yEJgZotk=B)_YQ}4tBB!1j=;t}~T4ZK!nd@>f>^fgp%q}~eAPEz$5 zyNRT3w&t_Oi0O)v@@#zbW}fHBmLYH>cjbaZl$uG&>{xNPb{ueivia;H^uO2f$lL3N zjCZ9}>Ey{6WX}pX`#Zn)tBVFKn7m?phv~4Zlen*CC^Vpa=>oDc7ogYlct|-raX>;$ z9*82;Fni5mSeJ39k2B-w1#y+cRZ7<=YA?lgt~XtZ8?-83E*`D;K6NAZVOZ9Q5+b z{BWz}Je+i2PX_J;#isd4(_zN!92r}^D2xW~9;;rF2O>>A*O-}I=}%rY$?*B>h0yoboBe}wCH1^`k>2ey8b9b0LWg__{V@d73Qisrk%rrHU*$`o5WH|cPMJ5 zz036@Y46dhknU5|-sKPYmgX+EX`EZfUv17x3F`aLrW6Ep8q}yuMP661(!8#MK-j_F z3`S>y0|RW7v?{1@lUL!q< z=+8C6j&=QZyu|Uaj!ICEh8U>)ysyu_1kTJr%o{|b)4?kQGD<+{A&IM$9#Pa@3Jlx}JdN8on=(PL4Nt`bjw?%tn5l8ltbXpEex|J9k*i;QHl>5GO+A(Y6-mvR19ibB z{CJY0^C1ITa{eKeItR@R`K%Y43}FlBHA$J0fN}$KI!4gp)-*>yy@UaGL%fUe#Y_=5QAl{_Wr3^dOpbk4HYbM;r zHvPjRk$hY+7DY;g{&|n!spt2te~%T%;yte#mZ4~|xNg9X>_0O2zCkgF>BrF}AbFBl zO6eIzt+ZZTKa%!2tqRGTqE^}$d`lzk*XRaQTScV9Rm1D=j^Fgwz^^Z1&E6xZ1x}ZB zTg@w7Z$VhzcRhw8CwbI_sWSWs+OCFKnTke*UU6`i@_CM{^x0(YYaOF2PeQaU0FL_( zAo_)DVVe#0;)a@^n-C3$Rg+^?Nii^6pMa{JGVp`nclXJ&Taw=h`Qg z0?M*<#?Qx6t&X6<69bO&rrQSEjSan-M|G)H{egW#Se(Lkazif5f$PyU>Olh zy>Vz^N7v-DFW~df`}1|^%}&njXMPzMW?RC2^8hF15MS|>K zJT6*5gTNN#2g|!x=jf?=3j3=FJU#OD$G6iHV5v{+VqZ%djvqd`_Pl1%@J*lFEr-I* z>9tVppNeuAmc#2<%ANka-O}a_ERo4WObUWNA>OV5v)V1d;FkT-SG3F zwg-Obu{fWbU-yQUq8NtxRfq!3PS#exjWQ6vP?&6)7*%qYVubwTW&N%0W2z2nyd zrTWY&oV^JvRO170Zbg}nPd%^n{gl66Pe3`1fW*pF2z3mZY0h}M--t!Ed*q|*?spEa>mDq#i97zf5ei?qYg~&wFtwHOIy!bafDo&&53*# zX~Rs{BxhX?FjL7CkK0IQg(2N+oSDNJ%VjN>8L5QVO}iu*>-7G82&K#Gcc5-&b0byk zxE|RpUgx3ojl@z){uH&+e&_m;v_EK7NIxlRrTxXXG|~oB`Eq<$-_@ z&Mrfd;k;|hI8cJ@Z((v$C?)0|p!A!>Qc8a)YNZX}`jNDOv?`>(6t&U@@hy$C`;{9_ z5}W0M7>)>f0F8IFS4V9*2vc~LzMsCv^Xx6qzxcB)`6u8>#v`Y&o2HDI-v9hk@}P&A zjbW!dWGu!FQS5El<}tjR!4odb=bQ?9Gq+k$Gogkl%OpJSTt9N9iZ@ZTD9!{y=^LB? z3$OY%|J^$c#xcs3W|Fm%NP#B8<|J#QD_eX)Fjbo505ymF-fd&Vrcdq`Q@7?zsLQ== z2dPmtjTmm92v8W&^~mmPrC{56pz#mg`bKIrJ#ItPbN04n9RPe)PYqN!QjhO$DmuF^@dy^=s zQ3Gue&9|AIBJc-S#5xw9?KF*gbFmB!`3SU)DNBH8wrfL}-1XU>JfV%6c5`i=2sIQVDNi3xlMNunlG}n)$jiFT`#Zu(b ziu6((-y&&wdTEm3(heI-OqqBF9w-K#6ZTDyETqOz-WI;ZkexX2{_UA2&tSUlzRhJO z?|^2#&KQ3qcUX~MPeef&9h#3%SmCV^8Fn%^H#i;n-)|kSJnY6b2y$d)7x|z@xP0pR zmxt30zLptLE@JhsgMl&EDnO)6919?me2ncZj2X9zyfHKBO<5>XBeei*gWw}s75Wr?sAiwd4H(9R%0ZQ+r zqvCVN>jDkgo>3-sm8R586*lAAEmnwPu&plq@3@zS+*;E+WJG*FuT1)Ayx4(!>yg&d{Y;VZ78Z|1H zzJE3TcAI$lbzDmgWZCo5u$3L-{ltm%5xt4`DL=`U5z<{9+0!Y!&YT6?llyIUy+of0 zG!9f+lmLF&=SEK?8wuXsBQ4EA@s(2?1sqPIoD~eB6c|p4k}Ff}t_?C~t|yzUZB0q} zLe6Qgn~h=a9`BNfPVzfyz{h{vvw;8?X2D0R!1YvlUDuEHZ2wxzPMxp31A$G=zKT*R z{c@$0Mp3I%Iv0a<%Ai#tWm44Yl*PBWPSlBhz~=@v? z+f0wm=-M;&Zucv2EHk#~g4YjN=`dwkb4Me#&FkmU^ly1^?~~cft&OAL`_)zEXDLhY zmA_xw7vkFiuDW(36JQ>ZH@ueFFhPY`U5H}7Y9h_(erg%x*E=2t3V>WHv|s$6?#=In zZz9Wdzn_WljbrMPX6+~uB$p`xUb)Z{ly?<(=S&_ITW=z<#I7dCpE1K zDUTx8Nt9Ug`4;J<>dk#iR#FL6%mD7g^JMqedHy5hxs8M+dMs+_z!K+vUYZG(y)bHi zir%@w#_R`;{vM1Y199ONl?MOF@VIWUF?lN(91giDwIL+>>o|+5l{Tz; z!@zwRcj$0fuAls38&wtx=$9y^LWi7p>E&emdS}?g2vK;6g zzP9CsG`OYUhNQ-fENg$5eW`?SFRJYlDCtOha-BRz#(`A8Svc$=whyz6bXxER8XY)z zaoiL;#-~-5)OZ2HU{P6$qEbONq$-+sHeFMN|FeCOO6iv?r80_Ioyxfwq*Dd03aOGJ z*GZIQtN51Yt?I}-+#8U!@{RuM*5*CI#upzF8ubfwOYA(*v9CURwZXD1F!DH-5eKQ! zWgjR=jof{^rzxvPW-87j)G$S5urhjNE$7+ca&BCV7S9%xzBy6>`_5CTqBhk%d6j!uHCkWv`(fBN2Bc~d zODX9X{y%AT=|6C3^=MT{`V_UdpaI|Fw;*rxeH@w(F4q@CgVQtQoYr5ilnb_?vf6t+ zmOPh3@;z~&KwkU6gxhT%=k+jWbE$@uXrhLzPe1GX-+70@e~f19D{R?!JUi*}OHYhi zbml@}SsV(hyq@OlGC5428kmps#AlEX3J;ROmQ29CCzgQh$ddt28s5}*TW!h&CONlU zW5|$&f-g^|QhP6R%D1P6JaGfhv`A@?HpG$@bD+b8GR)b$#8szR25|f3;x(o$Xx=oV zOT_dIiQ6`z?RcVp8qzOUN=6j5x1cc>gLE>XRUw&D)ZT(-d`oi+t}41z_&&gljTNWH zr$APF^G=Zsh_ErH4<>`&+$Q%6%9hsiU{p)$g7PKD9hwt(=Mf6Gmyk7VTp?OHCRU70 ztAyqJRJ}7B3XY9aHg?uiwaB`WwriCUyFZdV>i$%~(HUgT)1VYZx$~9@d-yweVph9Q zoO|<8?1&!~F!$QBS&VofS=a1BJ3OYku96cN1K_`DXPZ>D8oMlX}VkL&Bpr@y`Y^x=C!VNxV;IxDJc6Cx~Ue447|- zmL0B`sTy186SP{7&2g2MRP-k&+rvK2v!17;WKLo!B@2pLX)U>aB&`*#3dx$H_7=3^ zTl^M0RsAj2CsBv-X`n?WUf7A7c+@OFyPK&_Z|wABKn(3&jG^aeT_6TWe39v z{V-9jIa{OvH@K~3Y|KV&KRoK}8-3M}96sC=30clr`&aa`VUt%w+%DZpEco$#&cC|! zHlb2LFI72+Tk~@=EU9__-}q5UavA2WIpebd5d^ILcqonLq07OstN!#8$s|xu8+f^L!`wP>~61xYwWILV0Sm$V|R<)-GKoriXE(j z0Tz;ih)8$Cp`_k5vw6<@&&=m>-^Y8;tQ|9J{nq!#kf_r)R3_?^{}KMfs{oo`bMuND zT9b*c1-(+dMyvam1VoY*{c@FLO`~?5Z1`f(I@!{!l9kL0Blf!<5kGzOne1>JgtEd{!GwQb*D&Ceh+YO7eD6S> z?4(sxb7GFO2Ck;Q&fT>-AKwz(s(!2u<2Boyr#fsq_1n4h`UT-%i1D;l!Pa)_K5WK2 zefBtpZmd`ND4v5^hV1E@RbKxEn=w*~o0&M!8LSrla+TyjqjsHY^TnWbszb91Rh&ZD0`tG5x zYdpc%c0BhI3bA_~9&MO!%>;Pf-b5Xrls0Ybun6Fsp+Ed(MXX<|9F> zPrqCxHK0+uP7V2D&^k4uSt&Wv$k$2ix}5kWS|_#OpU%vHv;jNsEeMz3R#(HzmOcGY zlrZcd^VDnW6KcnL3pvvl2&2b#*vfp2B-BpM`xt0SOt$<}pD?PF(|?ko4i`YE6l@6f zl{k%x0J*?Yg{_CH&P=wg!+58-;p?wJ4djP2hmi`YG5r#i)PzRu?{3O}9R2QQG%Ka% zH1gjq`Calid;pu~hn098HY=7LTzj@M)B5RT=RL!)5*WL$so zzfqAN}8=93; zTN<_Z0l`J(`+($=?|=0RM*$)=EW&HGKN_tqcfp8%PGk^UfQOC!d~eTvSS=`=iI;t@ zJeu6@7an+M2Hj)wKyIF}1j2X^jy;+U5k+f9%HsH{{KwJn?nJXv>P(~dcX#2J`2BD3?$d`& z7JpIv`^}57tssOf!P(qPcC-e^R}4$};OzyXAQ&P(IzeduLMsdV(~0|Gp}Cb$PWs&b8g2;uJ?W6Hs_ju znT*n{&K>ff#3MEPdNeH4VbrWd)$1{^o{*ev37)%&b~bkiaK5BNjSquAQoOs;FHuR| zXw?4h|M-ui-`$;NrPPB){<|%?PNyfoL|d8dlt>{^S`hdK%SqxT>2Q_Iyq2B(j=WR9 z{)J>bFx7wvI^WPO?CHV9?%&|xs*_{qJEQ=wpgR4F$2o06+gq)?jXFs=f5#@2!t8tP zV&_v@PLy!SmvWoXei^YL`|~O2QFyCz;~s+-1t40cyepx_Z1=30lYQ}b-OaN*5sg7W zcVEQ=~+IFris3NXs9G zt}cAv7=thC6GC=HF3Jb%&$xSciADPb;_7&I zd<^=<@nSCR654mtMZ*G6$`}UVjZaydssAVH&$%mJ+3?1MdrM@?K9C-NErfs{5(tpb?!GrU0#d1z6{tMiqZ$o?}F4>38mQwxWreSD@%|t3XOInWx z-?Bg$Ss1(H`I^S6X3LY&e5&{9v(zrjJ^yC;f50PNclU?76``;)P0d4`O)$;|x{Lq{ z#dq-`+SS=VzwT)C-2}ZwmZoe=3MLtZA*|%$VY8jqjyGCB#R|aJ3q$yogBDbI! z=fRbyt0wOW_>Hphvp$n$$}e$W^g8pSTQafd_p`NT#gxHwuKr(6nQH!3YF~nL+7VMq zkpT@ks^)JY%J%0~D9ZawuYv!z9gWytHAsCa7%HhBjoQ%l=kcM?4WL;m4Wv;Ux