diff --git a/.idea/.gitignore b/.idea/.gitignore new file mode 100644 index 0000000..e7e9d11 --- /dev/null +++ b/.idea/.gitignore @@ -0,0 +1,2 @@ +# Default ignored files +/workspace.xml diff --git a/.idea/chromo.iml b/.idea/chromo.iml new file mode 100644 index 0000000..4f2c9af --- /dev/null +++ b/.idea/chromo.iml @@ -0,0 +1,15 @@ + + + + + + + + + + + + \ No newline at end of file diff --git a/.idea/inspectionProfiles/profiles_settings.xml b/.idea/inspectionProfiles/profiles_settings.xml new file mode 100644 index 0000000..105ce2d --- /dev/null +++ b/.idea/inspectionProfiles/profiles_settings.xml @@ -0,0 +1,6 @@ + + + + \ No newline at end of file diff --git a/.idea/misc.xml b/.idea/misc.xml new file mode 100644 index 0000000..7a85af6 --- /dev/null +++ b/.idea/misc.xml @@ -0,0 +1,7 @@ + + + + + + \ No newline at end of file diff --git a/.idea/modules.xml b/.idea/modules.xml new file mode 100644 index 0000000..ea896ee --- /dev/null +++ b/.idea/modules.xml @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file diff --git a/.idea/vcs.xml b/.idea/vcs.xml new file mode 100644 index 0000000..94a25f7 --- /dev/null +++ b/.idea/vcs.xml @@ -0,0 +1,6 @@ + + + + + + \ No newline at end of file diff --git a/chromo/Polymer.py b/chromo/Polymer.py deleted file mode 100644 index a75c112..0000000 --- a/chromo/Polymer.py +++ /dev/null @@ -1,243 +0,0 @@ -""" -Polymer class. - -Creates a polymer with a defined length and discretization. -""" -from pathlib import Path -from enum import Enum -import dataclasses -import io - -import numpy as np -import pandas as pd -from .util import dss_params, combine_repeat - - -@dataclasses.dataclass -class Epigenmark: - """Information about the chemical properties of an epigenetic mark.""" - - name: str - bind_energy: float - interaction_energy: float - chemical_potential: float - - -class State: - """Defines the epigenetic state of a given bead for a single mark.""" - - pass - - -class HP1_State(State, Enum): - """HP1 can be bound to none, either, or both nucleosome tails.""" - - NONE = 0 - LBOUND = 1 - RBOUND = 2 - BOTH = 3 - - -def make_epi_collection(epigenmarks): - """ - Takes a sequence of epigenetic marks and returns a summary DataFrame. - - Parameters - ---------- - epigenmarks : Sequence[Epigenmark] - The epigenetic marks to be summarized. - - Returns - ------- - pd.DataFrame - Columns are the properties of Epigenmark. - """ - df = pd.DataFrame(columns=['name', 'bind_energy', 'interaction_energy', - 'chemical_potential']) - for mark in epigenmarks: - df = df.append(dataclasses.asdict(mark), ignore_index=True) - return df - - -class Polymer: - """ - The positions and chemical state of a discrete polymer. - - The polymer carries around a set of coordinates ``r`` of shape - ``(num_beads, 3)``, a triad of material normals ``t_i`` for ``i`` in - ``{1,2,3}``, and some number of chemical states per bead. - - Its material properties are completely defined by these positions, chemical - states, and the length of polymer (in Kuhn lengths) simulated by each bead. - TODO: allow differential discretization, and decay to constant - discretization naturally. - - Since this codebase is primarily used to simulate DNA, information about - the chemical properties of each epigenetic mark are stored in `Epigenmark` - objects. - """ - - def __init__(self, name, r, t_3, t_2, epigenmarks, states, - length_per_bead): - """ - Parameters - ---------- - name : str - A name for convenient repr. Should be a valid filename. - r : (N, 3) array_like of float - The positions of each bead. - t_3 : (N, 3) array_like of float - The tangent vector to each bead in global coordinates. - t_2 : (N, 3) array_like of float - A material normal to each bead in global coordinates. - epigenmarks : (M, ) Sequence[Epigenmark] - Information about the chemical properties of each of the epigenetic - marks being tracked. - states : (N, M) array_like of State - State of each of the M epigenetic marks being tracked for each - bead. - length_per_bead : float - How many Kuhn lengths of polymer are simulated by each bead. - """ - self.name = name - self.r = r - self.t_3 = t_3 - self.t_2 = t_2 - self.epigenmarks = make_epi_collection(epigenmarks) - self.states = states - self.length_per_bead = length_per_bead - # Load the polymer parameters for the conformational energy - self.sim_type, self.eps_bend, self.eps_par, self.eps_perp, \ - self.gamma, self.eta = self._find_parameters(length_per_bead) - - @classmethod - def from_files(cls, *, r=None, t_3=None, t_2=None, marks=None, states=None, - lengths=None): - """ - Instantiate a Polymer from input files. - - Each of the r, t_3, t_2, epigenmarks, states, and length_per_bead - arrays requires a separate input file. - """ - r = pd.read_csv(r, header=None) if r is not None else None - t_3 = pd.read_csv(t_3, header=None) if t_3 is not None else None - t_2 = pd.read_csv(t_2, header=None) if t_2 is not None else None - states = pd.read_csv(states) if states is not None else None - marks = pd.read_csv(r, header=None) if r is not None else None - return cls(r, t_3, t_2, marks, states, lengths) - - @classmethod - def straight_line_in_x(cls, name, marks, states, num_beads, - length_per_bead): - r = np.zeros((num_beads, 3)) - r[:, 0] = length_per_bead * np.arange(num_beads) - t_3 = np.zeros((num_beads, 3)) - t_3[:, 0] = 1 - t_2 = np.zeros((num_beads, 3)) - t_2[:, 1] = 1 - return cls(name, r, t_3, t_2, marks, states, length_per_bead) - - def write_marks(self, marks_file, **kwargs): - self.epigenmarks.to_csv(marks_file, **kwargs) - - def __repr__(self): - s = io.Bytes() - self.write_repr(s) - return s.decode() - - def write_repr(self, path): - np.savez(path, r=self.r, t_2=self.t_2, t_3=self.t_3, - states=self.states, marks=self.epigenmarks.values, - length_per_bead=self.length_per_bead) - - @classmethod - def from_repr(cls, path): - npz = np.load(path) - mark_info_arr = npz['marks'] - marks = [Epigenmark(*info) for info in mark_info_arr] - return cls(r=npz['r'], t_3=npz['t_3'], t_2=npz['t_2'], - states=npz['states'], marks=marks, - length_per_bead=npz['length_per_bead']) - - @property - def num_marks(self): - return len(self.epigenmarks) - - @property - def num_beads(self): - return self.r.shape[0] - - def __str__(self): - return f"Polymer<{self.name}, nbeads={self.num_beads}, " \ - f"nmarks={self.num_marks}>" - - @staticmethod - def _find_parameters(length_bead): - """Determine the parameters for the elastic forces based on ssWLC with twist.""" - sim_type = "sswlc" - lp_dna = 53 # Persistence length of DNA in nm - length_dim = length_bead * 0.34 / lp_dna # Non-dimensionalized length per bead - - # Determine the parameter values using linear interpolation of the parameter table - eps_bend = np.interp(length_dim, dss_params[:, 0], dss_params[:, 1]) / length_dim - gamma = np.interp(length_dim, dss_params[:, 0], dss_params[:, 2]) * length_dim * lp_dna - eps_par = np.interp(length_dim, dss_params[:, 0], dss_params[:, 3]) / (length_dim * lp_dna**2) - eps_perp = np.interp(length_dim, dss_params[:, 0], dss_params[:, 4]) / (length_dim * lp_dna**2) - eta = np.interp(length_dim, dss_params[:, 0], dss_params[:, 5]) / lp_dna - - return sim_type, eps_bend, eps_par, eps_perp, gamma, eta - - def compute_dE(self, ind0, indf, r_poly_trial, t3_poly_trial, t2_poly_trial, - states_trial): - """ - Compute change in energy of polymer state from proposed new state. - """ - delta_energy_poly = 0 - # Calculate contribution to polymer energy at the ind0 position - if ind0 != 0: - delta_r_trial = r_poly_trial[0, :] - self.r[ind0 - 1, :] - delta_r_par_trial = np.dot(delta_r_trial, self.t_3[ind0 - 1, :]) - delta_r_perp_trial = delta_r_trial - delta_r_par_trial * self.t_3[ind0 - 1, :] - - delta_r = self.r[ind0, :] - self.r[ind0 - 1, :] - delta_r_par = np.dot(delta_r, self.t_3[ind0 - 1, :]) - delta_r_perp = delta_r - delta_r_par * self.t_3[ind0 - 1, :] - - bend_vec_trial = (t3_poly_trial[0, :] - self.t_3[ind0 - 1, :] - - self.eta * delta_r_perp_trial) - bend_vec = (self.t_3[ind0, :] - self.t_3[ind0 - 1, :] - - self.eta * delta_r_perp) - - delta_energy_poly += (0.5 * self.eps_bend * np.dot(bend_vec_trial, bend_vec_trial) - + 0.5 * self.eps_par * (delta_r_par_trial - self.gamma) ** 2 - + 0.5 * self.eps_perp * np.dot(delta_r_perp_trial, delta_r_perp_trial)) - delta_energy_poly -= (0.5 * self.eps_bend * np.dot(bend_vec, bend_vec) - + 0.5 * self.eps_par * (delta_r_par - self.gamma) ** 2 - + 0.5 * self.eps_perp * np.dot(delta_r_perp, delta_r_perp)) - - # Calculate contribution to polymer energy at the indf position - if indf != self.num_beads: - - delta_r_trial = self.r[indf, :] - r_poly_trial[indf - ind0 - 1, :] - delta_r_par_trial = np.dot(delta_r_trial, t3_poly_trial[indf - ind0 - 1, :]) - delta_r_perp_trial = delta_r_trial - delta_r_par_trial * t3_poly_trial[indf - ind0 - 1, :] - - delta_r = self.r[indf, :] - self.r[indf - 1, :] - delta_r_par = np.dot(delta_r, self.t_3[indf - 1, :]) - delta_r_perp = delta_r - delta_r_par * self.t_3[indf - 1, :] - - bend_vec_trial = (self.t_3[indf, :] - t3_poly_trial[indf - ind0 - 1, :] - - self.eta * delta_r_perp_trial) - bend_vec = (self.t_3[indf, :] - self.t_3[indf - 1, :] - - self.eta * delta_r_perp) - - delta_energy_poly += (0.5 * self.eps_bend * np.dot(bend_vec_trial, bend_vec_trial) - + 0.5 * self.eps_par * (delta_r_par_trial - self.gamma) ** 2 - + 0.5 * self.eps_perp * np.dot(delta_r_perp_trial, delta_r_perp_trial)) - delta_energy_poly -= (0.5 * self.eps_bend * np.dot(bend_vec, bend_vec) - + 0.5 * self.eps_par * (delta_r_par - self.gamma) ** 2 - + 0.5 * self.eps_perp * np.dot(delta_r_perp, delta_r_perp)) - - return delta_energy_poly - - diff --git a/chromo/components.py b/chromo/components.py new file mode 100644 index 0000000..05161ce --- /dev/null +++ b/chromo/components.py @@ -0,0 +1,286 @@ +""" +Components that will make up our simulation. + +Various types of polymers, solvents, and other simulation components should be +defined here. +""" +from enum import Enum + +import numpy as np +import pandas as pd + +from .util import dss_params + + +class State(Enum): + """Defines the chemical state of a given bead for a single mark.""" + + pass + + +class HP1_State(State): + """HP1 can be bound to none, either, or both nucleosome tails.""" + + NONE = 0 + LBOUND = 1 + RBOUND = 2 + BOTH = 3 + + +class Polymer: + """ + The positions and chemical state of a discrete polymer. + + The polymer carries around a set of coordinates ``r`` of shape + ``(num_beads, 3)``, a triad of material normals ``t_i`` for ``i`` in + ``{1,2,3}``, and some number of chemical states per bead. + + Its material properties are completely defined by these positions, chemical + states, and the length of polymer (in Kuhn lengths) simulated by each bead. + TODO: allow differential discretization, and decay to constant + discretization naturally. + + Since this codebase is primarily used to simulate DNA, information about + the chemical properties of each epigenetic mark are stored in `Epigenmark` + objects. + """ + + _arrays = 'r', 't3', 't2', 'states', 'bead_length' + """Track which arrays will be saved to file.""" + _3d_arrays = 'r', 't3', 't2' + """Track which arrays need to have multi-indexed values (x, y, z).""" + + def __init__(self, name, r, *, bead_length, t3=None, t2=None, + states=None, mark_names=None): + """ + Construct a polymer. + + Parameters + ---------- + name : str + A name for convenient repr. Should be a valid filename. + r : (N, 3) array_like of float + The positions of each bead. + t3 : (N, 3) array_like of float + The tangent vector to each bead in global coordinates. + t2 : (N, 3) array_like of float + A material normal to each bead in global coordinates. + states : (N, M) array_like of int + State of each of the M epigenetic marks being tracked for each + bead. + mark_names : (M, ) str or Sequence[str] + The name of each chemical modification tracked in `states`, for + each of tracking which mark is which. + bead_length : float or (N,) array_like of float + The amount of polymer path length between this bead and the next + bead. For now, a constant value is assumed (the first value if an + array is passed). + """ + self.name = name + self.r = r + self.t3 = t3 + self.t2 = t2 + self.states = states + self.mark_names = mark_names + if states is not None: + # does what atleast_2d(axis=1) should do (but doesn't) + self.states = self.states.reshape(states.shape[0], -1) + num_beads, num_marks = self.states.shape + if num_marks != len(mark_names): + raise ValueError("Each chemical state must be given a name.") + if num_beads != len(self.r): + raise ValueError("Initial epigenetic state of wrong length.") + # first if makes sure len will work in second if + if not issubclass(type(bead_length), np.ndarray): + bead_length = np.atleast_1d(bead_length) + if len(bead_length) == 1: + bead_length = np.broadcast_to(bead_length, (self.num_beads, 1)) + self.bead_length = bead_length + # Load the polymer parameters for the conformational energy + # for now, the array-values of bead_length are ignored + self.eps_bend, self.eps_par, self.eps_perp, self.gamma, self.eta \ + = self._find_parameters(self.bead_length[0]) + + @classmethod + def from_csv(cls, csv_file): + """Construct Polymer from CSV file.""" + df = pd.read_csv(csv_file) + return cls.from_dataframe(df) + + @classmethod + def straight_line_in_x(cls, name, num_beads, bead_length, **kwargs): + """Construct polymer initialized uniformly along the positve x-axis.""" + r = np.zeros((num_beads, 3)) + r[:, 0] = bead_length * np.arange(num_beads) + t3 = np.zeros((num_beads, 3)) + t3[:, 0] = 1 + t2 = np.zeros((num_beads, 3)) + t2[:, 1] = 1 + return cls(name, r, t3=t3, t2=t2, bead_length=bead_length, **kwargs) + + def to_dataframe(self): + """ + Write canonical CSV representation of the Polymer to file. + + The top-level multiindex values for the columns should correspond to + the kwargs, to simplify unpacking the structure (and to make for easy + accessing of the dataframe. + """ + arrays = {name: self.__dict__[name] for name in self._arrays + if self.__dict__[name] is not None} + vector_arrs = {} + regular_arrs = {} + for name, arr in arrays.items(): + if name in self._3d_arrays: + vector_arrs[name] = arr + # special-cased below to get correct epigenmark names + elif name != 'states': + regular_arrs[name] = arr + # first construct the parts of the DataFrame that need a multi-index + # TODO: remove "list" after numpy fixes + # https://github.com/numpy/numpy/issues/17305 + vector_arr = np.concatenate(list(vector_arrs.values()), axis=1) + # vector_arr = np.concatenate(vector_arrs.values(), axis=1) + vector_index = pd.MultiIndex.from_product([vector_arrs.keys(), + ('x', 'y', 'z')]) + vector_df = pd.DataFrame(vector_arr, columns=vector_index) + states_index = pd.MultiIndex.from_tuples( + [('states', name) for name in self.mark_names]) + states_df = pd.DataFrame(self.states, columns=states_index) + df = pd.concat([vector_df, states_df], axis=1) + # now throw in the remaining columns one-by-one + for name, arr in regular_arrs.items(): + df[name] = arr + return df + + @classmethod + def from_dataframe(cls, df, name=None): + """Construct Polymer object from DataFrame. Inverts `.to_dataframe`.""" + # top-level multiindex values correspond to kwargs + kwnames = np.unique(df.columns.get_level_values(0)) + kwargs = {name: df[name].to_numpy() for name in kwnames} + # extract names of each epigenetic state from multi-index + if 'states' in df: + mark_names = df['states'].columns.to_numpy() + kwargs['mark_names'] = mark_names + return cls(name, **kwargs) + + @classmethod + def from_file(cls, path, name=None): + """Construct Polymer object from string representation.""" + if name is None: + name = path.name + return cls.from_dataframe( + pd.read_csv(path, header=[0, 1], index_col=0), + name + ) + + def to_csv(self, path): + """Save Polymer object to CSV file as DataFrame.""" + return self.to_dataframe().to_csv(path) + + def to_file(self, path): + """Synonym for *to_csv* to conform to `make_reproducible` spec.""" + return self.to_csv(path) + + @property + def num_marks(self): + """Return number of states tracked per bead.""" + return self.states.shape[1] + + @property + def num_beads(self): + """Return number of beads in the polymer.""" + return self.r.shape[0] + + def __str__(self): + """Return string representation of the Polymer.""" + return f"Polymer<{self.name}, nbeads={self.num_beads}, " \ + f"nmarks={self.num_marks}>" + + @staticmethod + def _find_parameters(length_bead): + """Look up elastic parameters of ssWLC for each bead_length.""" + lp_dna = 53 # Persistence length of DNA in nm + # Non-dimensionalized length per bead + length_dim = length_bead * 0.34 / lp_dna + + # Determine the parameter values using linear interpolation of the + # parameter table + eps_bend = np.interp(length_dim, dss_params[:, 0], dss_params[:, 1]) \ + / length_dim + gamma = np.interp(length_dim, dss_params[:, 0], dss_params[:, 2]) \ + * length_dim * lp_dna + eps_par = np.interp(length_dim, dss_params[:, 0], dss_params[:, 3]) \ + / (length_dim * lp_dna**2) + eps_perp = np.interp(length_dim, dss_params[:, 0], dss_params[:, 4]) \ + / (length_dim * lp_dna**2) + eta = np.interp(length_dim, dss_params[:, 0], dss_params[:, 5]) \ + / lp_dna + + return eps_bend, eps_par, eps_perp, gamma, eta + + def compute_dE(self, ind0, indf, r_poly_trial, t3_poly_trial, + t2_poly_trial, states_trial): + """Compute change in polymer energy moving to proposed new state.""" + delta_energy_poly = 0 + # Calculate contribution to polymer energy at the ind0 position + if ind0 != 0: + delta_r_trial = r_poly_trial[0, :] - self.r[ind0 - 1, :] + delta_r_par_trial = np.dot(delta_r_trial, self.t3[ind0 - 1, :]) + delta_r_perp_trial = delta_r_trial \ + - delta_r_par_trial * self.t3[ind0 - 1, :] + + delta_r = self.r[ind0, :] - self.r[ind0 - 1, :] + delta_r_par = np.dot(delta_r, self.t3[ind0 - 1, :]) + delta_r_perp = delta_r - delta_r_par * self.t3[ind0 - 1, :] + + bend_vec_trial = (t3_poly_trial[0, :] - self.t3[ind0 - 1, :] + - self.eta * delta_r_perp_trial) + bend_vec = (self.t3[ind0, :] - self.t3[ind0 - 1, :] + - self.eta * delta_r_perp) + + delta_energy_poly += ( + 0.5 * self.eps_bend * np.dot(bend_vec_trial, bend_vec_trial) + + 0.5 * self.eps_par * (delta_r_par_trial - self.gamma) ** 2 + + 0.5 * self.eps_perp * np.dot(delta_r_perp_trial, + delta_r_perp_trial) + ) + delta_energy_poly -= ( + 0.5 * self.eps_bend * np.dot(bend_vec, bend_vec) + + 0.5 * self.eps_par * (delta_r_par - self.gamma) ** 2 + + 0.5 * self.eps_perp * np.dot(delta_r_perp, delta_r_perp) + ) + + # Calculate contribution to polymer energy at the indf position + if indf != self.num_beads: + + delta_r_trial = self.r[indf, :] - r_poly_trial[indf - ind0 - 1, :] + delta_r_par_trial = np.dot(delta_r_trial, + t3_poly_trial[indf - ind0 - 1, :]) + delta_r_perp_trial = delta_r_trial \ + - delta_r_par_trial * t3_poly_trial[indf - ind0 - 1, :] + + delta_r = self.r[indf, :] - self.r[indf - 1, :] + delta_r_par = np.dot(delta_r, self.t3[indf - 1, :]) + delta_r_perp = delta_r - delta_r_par * self.t3[indf - 1, :] + + bend_vec_trial = (self.t3[indf, :] + - t3_poly_trial[indf - ind0 - 1, :] + - self.eta * delta_r_perp_trial) + bend_vec = (self.t3[indf, :] - self.t3[indf - 1, :] + - self.eta * delta_r_perp) + + delta_energy_poly += ( + 0.5 * self.eps_bend * np.dot(bend_vec_trial, bend_vec_trial) + + 0.5 * self.eps_par * (delta_r_par_trial - self.gamma) ** 2 + + 0.5 * self.eps_perp * np.dot(delta_r_perp_trial, + delta_r_perp_trial) + ) + delta_energy_poly -= ( + 0.5 * self.eps_bend * np.dot(bend_vec, bend_vec) + + 0.5 * self.eps_par * (delta_r_par - self.gamma) ** 2 + + 0.5 * self.eps_perp * np.dot(delta_r_perp, delta_r_perp) + ) + + return delta_energy_poly diff --git a/chromo/Field.py b/chromo/fields.py similarity index 50% rename from chromo/Field.py rename to chromo/fields.py index d027e0e..1ddaa12 100644 --- a/chromo/Field.py +++ b/chromo/fields.py @@ -1,12 +1,12 @@ """ -Field class +Fields discretize space to efficiently track changes in Mark energy. Creates a field object that contains parameters for the field calculations and functions to generate the densities - """ - import numpy as np +import pandas as pd + from .util import combine_repeat @@ -16,10 +16,21 @@ class FieldBase: Must be subclassed to be useful. """ + + _field_descriptors = [] + + @property + def name(self): + """For now, there's only one field per sim, so classname works.""" + return self.__class__.__name__ + def __init__(self): + """Construct a field holding no Polymers, tracking no Marks.""" self.polymers = [] + self.marks = pd.DataFrame() def __str__(self): + """Print representation of empty field.""" return "Field<>" def __contains__(self, poly): @@ -62,9 +73,80 @@ def compute_dE(self, poly, ind0, indf, r, t3, t2, states): pass +class Reconstructor: + """ + Defer construction of `Field` until after `Polymer`/`Mark` instances. + + Constructs a kwargs object that can be re-passed to the appropriate `Field` + constructor when they become available. + """ + + def __init__(self, cls, **kwargs): + """Construct our Reconstructor.""" + self.field_constructor = cls + self.kwargs = kwargs + + def finalize(self, polymers, marks): + """Finish construction of appropriate `Field` object.""" + return self.field_constructor(polymers=polymers, marks=marks, + **self.kwargs) + + @classmethod + def from_file(cls, path): + """Assume class name is encoded in file name.""" + constructor = globals()[path.name] + # read info as a series + kwargs = pd.read_csv(path).iloc[0].to_dict() + return cls(constructor, **kwargs) + + def __call__(self, polymers, marks): + """Synonym for `Reconstructor.finalize`.""" + return self.finalize(polymers, marks) + + class UniformDensityField(FieldBase): - def __init__(self, polymers, x_width, nx, y_width, ny, z_width, nz): + """ + Rectilinear discretization of a rectangular box. + + Computes field energies at the corners of each box. The bead in each box + contributes "mass" to each vertex linearly based on its position in the + box. This is much more stable numerically than just using the boxes as + "bins" as would be done in a e.g. finite-differences discretization. + """ + + _field_descriptors = ['x_width', 'nx', 'y_width', 'ny', 'z_width', 'nz'] + + def __init__(self, polymers, marks, x_width, nx, y_width, ny, z_width, + nz): + """ + Construct a UniformDensityField. + + Parameters + ---------- + polymers : Sequence[Polymer] + Name of (or object representing) each polymer in the field. + marks : pd.DataFrame + Output of `chromo.marks.make_mark_collection` for the `Mark`s to be + used. + x_width : float + Width of the box containing the field in the x-direction. + nx : int + Number of bins in the x-direction. + y_width : float + Width of the box containing the field in the y-direction. + ny : int + Number of bins in the y-direction. + z_width : float + Width of the box containing the field in the z-direction. + nz : int + Number of bins in the z-direction. + """ self.polymers = polymers + self.marks = marks + for poly in polymers: + if poly.num_marks != len(marks): + raise NotImplementedError("For now, all polymers must use all" + " of the same marks.") self.x_width = x_width self.y_width = y_width self.z_width = z_width @@ -78,18 +160,58 @@ def __init__(self, polymers, x_width, nx, y_width, ny, z_width, nz): self.vol_bin = x_width * y_width * z_width / self.n_bins self.bin_index = UniformDensityField._get_corner_bin_index( self.nx, self.ny, self.nz) - num_epigenmarks = np.array([len(p.epigenmarks) for p in polymers]) - if np.any(num_epigenmarks != num_epigenmarks[0]): - raise ValueError("All Polymers must have the same number of " - "epigenetic marks for now.") # one column of density for each state and one for the actual bead # density of the polymer - self.density = np.zeros((self.n_bins, num_epigenmarks[0] + 1), 'd') + self.num_marks = len(marks) + self.density = np.zeros((self.n_bins, self.num_marks + 1), 'd') self._recompute_field() # initialize the density field + def to_file(self, path): + """Save Field description and Polymer/Mark names to CSV as Series.""" + rows = {name: self.__dict__[name] for name in self._field_descriptors} + for i, polymer in enumerate(self.polymers): + rows[polymer.name] = 'polymer' + for i, mark in self.marks.iterrows(): + # careful! mark.name is the Series.name attribute + rows[mark['name']] = 'mark' + # prints just key,value for each key in rows + return pd.Series(rows).to_csv(path, header=None) + + @classmethod + def from_file(cls, path, polymers, marks): + """Recover field saved with `.to_file`. Requires Poly/Mark inputs.""" + # the 0th column will be the index, the 1th column holds the data... + field_series = pd.read_csv(path, header=None, index_col=0)[1] + kwargs = pd.to_numeric(field_series[cls._field_descriptors]).to_dict() + polymer_names = field_series[field_series == 'polymer'].index.values + mark_names = field_series[field_series == 'mark'].index.values + + err_prefix = f"Tried to instantiate class:{cls} from file:{path} with " + if len(polymers) != len(polymer_names): + raise ValueError(err_prefix + f"{len(polymers)} polymers, but " + f" there are {len(polymer_names)} listed.") + for polymer in polymers: + if polymer.name not in polymer_names: + raise ValueError(err_prefix + f"polymer:{polymer.name}, but " + " this polymer was not present in file.") + if len(marks) != len(mark_names): + raise ValueError(err_prefix + f"{len(marks)} marks, but " + f" there are {len(mark_names)} listed.") + for i, mark in marks.iterrows(): + if mark['name'] not in mark_names: + raise ValueError(err_prefix + f"mark:{mark}, but " + " this mark was not present in file.") + + return cls(polymers=polymers, marks=marks, **kwargs) + @staticmethod def _get_corner_bin_index(nx, ny, nz): - """Setup the index array for density bins corners.""" + """ + Set up the index array for density bins corners. + + TODO: this should just be a higher-dimensional array to avoid having to + do any of this math. + """ bin_index = np.zeros((nx * ny * nz, 8), 'd') count = 0 for index_z in range(nz): @@ -121,9 +243,10 @@ def _get_corner_bin_index(nx, ny, nz): return bin_index def __str__(self): + """Print summary of the UniformDensityField.""" n_poly = len(self.polymers) - return "UniformRectField" + return f"UniformDensityField" def _recompute_field(self, check_consistency=False): new_density = self.density.copy() @@ -138,49 +261,63 @@ def _recompute_field(self, check_consistency=False): " bug in this code.") def compute_dE(self, poly, ind0, indf, r, t3, t2, states): + """ + Compute change in energy based on proposed new polymer location. + + Requires the polymer to be moved in order to compare to old state. + """ # even if the move does not change the states, we cannot ignore them # because they're needed to compute the energy at any point along the # polymer if states is None: states = poly.states density_poly, index_xyz = self._calc_density( - poly.r[ind0:indf, :], poly.states, ind0, indf) - density_poly_trial, index_xyz_trial = self._calc_density(r, states, ind0, indf) - delta_density_poly_total = np.concatenate((density_poly_trial, -density_poly)) - delta_index_xyz_total = np.concatenate((index_xyz_trial, index_xyz)).astype(int) - delta_density, delta_index_xyz = combine_repeat(delta_density_poly_total, delta_index_xyz_total) - - delta_energy_epigen = 0 - for i, (_, epi_info) in enumerate(poly.epigenmarks.iterrows()): - delta_energy_epigen += 0.5 * epi_info.interaction_energy * np.sum( - (delta_density[:, i + 1] + self.density[delta_index_xyz, i + 1]) ** 2 - - self.density[delta_index_xyz, i + 1] ** 2) - return delta_energy_epigen + poly.r[ind0:indf, :], poly.states, ind0, indf) + density_poly_trial, index_xyz_trial = self._calc_density( + r, states, ind0, indf) + delta_density_poly_total = np.concatenate(( + density_poly_trial, -density_poly)) + delta_index_xyz_total = np.concatenate( + (index_xyz_trial, index_xyz)).astype(int) + delta_density, delta_index_xyz = combine_repeat( + delta_density_poly_total, delta_index_xyz_total) + + dE_marks = 0 + for i, (_, mark_info) in enumerate(self.marks.iterrows()): + dE_marks += 0.5 * mark_info.interaction_energy * np.sum( + ( + delta_density[:, i + 1] + + self.density[delta_index_xyz, i + 1] + ) ** 2 + - self.density[delta_index_xyz, i + 1] ** 2 + ) + return dE_marks def _calc_density(self, r_poly, states, ind0, indf): - num_epigenmark = states.shape[1] + num_marks = states.shape[1] # Find the (0,0,0) bins for the beads and the associated weights - x_poly_box = (r_poly[:, 0] - 0.5 * self.dx - - self.x_width * np.floor((r_poly[:, 0] - 0.5 * self.dx) / self.x_width)) + x_poly_box = (r_poly[:, 0] - 0.5 * self.dx - self.x_width * np.floor( + (r_poly[:, 0] - 0.5 * self.dx) / self.x_width)) index_x = np.floor(self.nx * x_poly_box / self.x_width).astype(int) weight_x = 1 - (x_poly_box / self.dx - index_x) - y_poly_box = (r_poly[:, 1] - 0.5 * self.dy - - self.y_width * np.floor((r_poly[:, 1] - 0.5 * self.dy) / self.y_width)) + y_poly_box = (r_poly[:, 1] - 0.5 * self.dy - self.y_width * np.floor( + (r_poly[:, 1] - 0.5 * self.dy) / self.y_width)) index_y = np.floor(self.ny * y_poly_box / self.y_width).astype(int) weight_y = 1 - (y_poly_box / self.dy - index_y) - z_poly_box = (r_poly[:, 2] - 0.5 * self.dz - - self.z_width * np.floor((r_poly[:, 2] - 0.5 * self.dz) / self.z_width)) + z_poly_box = (r_poly[:, 2] - 0.5 * self.dz - self.z_width * np.floor( + (r_poly[:, 2] - 0.5 * self.dz) / self.z_width)) index_z = np.floor(self.nz * z_poly_box / self.z_width).astype(int) weight_z = 1 - (z_poly_box / self.dz - index_z) - index_x0y0z0 = index_x + self.nx * index_y + self.nx * self.ny * index_z + index_x0y0z0 = index_x + self.nx*index_y + self.nx*self.ny*index_z # Generate the weight array for all of the bins containing the beads + #TODO if you clean this up, please use append instead weight = weight_x * weight_y * weight_z weight = np.concatenate((weight, (1 - weight_x) * weight_y * weight_z)) weight = np.concatenate((weight, weight_x * (1 - weight_y) * weight_z)) @@ -190,6 +327,7 @@ def _calc_density(self, r_poly, states, ind0, indf): weight = np.concatenate((weight, weight_x * (1 - weight_y) * (1 - weight_z))) weight = np.concatenate((weight, (1 - weight_x) * (1 - weight_y) * (1 - weight_z))) + #TODO this can be done in one line index_xyz_total = index_x0y0z0 index_xyz_total = np.concatenate((index_xyz_total, self.bin_index[index_x0y0z0, 1])).astype(int) index_xyz_total = np.concatenate((index_xyz_total, self.bin_index[index_x0y0z0, 2])).astype(int) @@ -199,14 +337,15 @@ def _calc_density(self, r_poly, states, ind0, indf): index_xyz_total = np.concatenate((index_xyz_total, self.bin_index[index_x0y0z0, 6])).astype(int) index_xyz_total = np.concatenate((index_xyz_total, self.bin_index[index_x0y0z0, 7])).astype(int) - density_total = np.zeros((8 * (indf - ind0), num_epigenmark + 1), 'd') + density_total = np.zeros((8 * (indf - ind0), num_marks + 1), 'd') density_total[:, 0] = weight - for ind_epigen in range(num_epigenmark): + for ind_mark in range(num_marks): for ind_corner in range(8): ind_corner0 = ind_corner * (indf - ind0) ind_cornerf = ind_corner0 + indf - ind0 - density_total[ind_corner0:ind_cornerf, ind_epigen + 1] = \ - weight[ind_corner0:ind_cornerf] * states[ind0:indf, ind_epigen] + density_total[ind_corner0:ind_cornerf, ind_mark + 1] = \ + weight[ind_corner0:ind_cornerf] \ + * states[ind0:indf, ind_mark] # Combine repeat entries in the density array density, index_xyz = combine_repeat(density_total, index_xyz_total) diff --git a/chromo/image.py b/chromo/image.py new file mode 100644 index 0000000..948932f --- /dev/null +++ b/chromo/image.py @@ -0,0 +1,189 @@ +r""" +Imaging module - Generate image files from simulations for PyMol imaging + +Notes +----- + +""" + + +import numpy as np + + +def gen_pymol_file(r_poly, meth_seq = np.array([]), hp1_seq = np.array([]), limit_n = False, n_max = 100000, + max_method = 'mid_slice', ind_save = np.array([]), + add_com = False, filename='r_poly.pdb', ring=False): + r""" + + Parameters + ---------- + r_poly : (num_beads, 3) float + Conformation of the chain subjected to active-Brownian forces + meth_seq : (num_beads) int + Epigenetic sequence (0 = No tails methylated, 1 = 1 tail methylated, 2 = 2 tails methylated) + hp1_seq : (num_beads) int + Number of tails with HP1 bound + limit_n : Boolean + + filename : str + File name to write the pdb file + ring : bool + Boolean to close the polymer into a ring + + Returns + ------- + none + + """ + + # Open the file + f = open(filename, 'w') + + atomname1 = "A1" # Chain atom type + resname = "SSN" # Type of residue (UNKnown/Single Stranded Nucleotide) + chain = "A" # Chain identifier + resnum = 1 + numresidues = len(r_poly[:, 0]) + descrip = "Pseudo atom representation of DNA" + chemicalname = "Body and ribbon spatial coordinates" + if len(meth_seq) == 0: + image_meth_seq = False + else: + image_meth_seq = True + + # Determine the bead indices to reduce the total represented to n_max if limit_n True + + if not ind_save.size == 0: + connect_left, connect_right = find_connect(ind_save, ring) + else: + if n_max >= numresidues: + limit_n = False + if limit_n: + if max_method == 'mid_slice': + ind_save, connect_left, connect_right = find_ind_mid_slice(r_poly, n_max, ring) + else: + ind_save, connect_left, connect_right = find_ind_total(r_poly, ring) + + # Write the preamble to the pymol file + + f.write('HET %3s %1s%4d %5d %-38s\n' % (resname, chain, resnum, numresidues, descrip)) + f.write('HETNAM %3s %-50s\n' % (resname, chemicalname)) + f.write('FORMUL 1 %3s C20 N20 P21\n' % (resname)) + + # Write the conformation to the pymol file + + count = 0 + for ind in range(numresidues): + if ind_save[ind] == 1: + if image_meth_seq: + atomname = 'A' + str(int(meth_seq[ind])) + else: + atomname = atomname1 + f.write('ATOM%7d %4s %3s %1s %8.3f%8.3f%8.3f%6.2f%6.2f C\n' % + (count + 1, atomname, resname, chain, r_poly[ind, 0], r_poly[ind, 1], r_poly[ind, 2], 1.00, 1.00)) + count += 1 + numresidues_save = count + + # Add a nucleus bead to the pymol file + + atomname = 'AN' + chain = 'B' + r_com = np.mean(r_poly, axis = 0) + if add_com: + f.write('ATOM%7d %4s %3s %1s %8.3f%8.3f%8.3f%6.2f%6.2f C\n' % + (count + 1, atomname, resname, chain, r_com[0], r_com[1], r_com[2], 1.00, 1.00)) + + # Define the connectivity in the chain + + count = 0 + for ind in range(numresidues): + if ind_save[ind] == 1: + if ind == 0: + ind_left = numresidues_save - 1 + else: + ind_left = count - 1 + if ind == numresidues - 1: + ind_right = 0 + else: + ind_right = count + 1 + + if connect_left[ind] == 1 and connect_right[ind] == 1: + f.write('CONECT%5d%5d%5d\n' % (count + 1, ind_left + 1, ind_right + 1)) + elif connect_left[ind] == 1 and connect_right[ind] == 0: + f.write('CONECT%5d%5d\n' % (count + 1, ind_left + 1)) + elif connect_left[ind] == 0 and connect_right[ind] == 1: + f.write('CONECT%5d%5d\n' % (count + 1, ind_right + 1)) + elif connect_left[ind] == 0 and connect_right[ind] == 0: + f.write('CONECT%5d\n' % (count + 1)) + count += 1 + + # Close the file + f.write('END') + f.close() + + return + + +def find_ind_mid_slice(r_poly, n_max, ring): + r""" + + """ + + # Find the value of the x-coordinate and sort the beads according to distance from mean + x_ave = np.mean(r_poly[:, 0]) + delta_x = r_poly[:, 0] - x_ave +# ind_sort = np.argsort(np.abs(delta_x)) + ind_sort = np.argsort(-delta_x) + + # Select the first n_max beads based on distance from the mean + ind_save = np.zeros(len(r_poly[:, 0])) + + for ind in range(n_max): + ind_save[ind_sort[ind]] = 1 + + # Determine whether adjacent beads are saved to define connections + connect_left, connect_right = find_connect(ind_save, ring) + + return ind_save, connect_left, connect_right + + +def find_ind_total(r_poly, ring): + r""" + + """ + + ind_save = np.ones(len(r_poly[:, 0])) + connect_left, connect_right = find_connect(ind_save, ring) + + return ind_save, connect_left, connect_right + + +def find_connect(ind_save, ring): + r""" + + + """ + connect_left = np.zeros(len(ind_save)) + connect_right = np.zeros(len(ind_save)) + + for ind in range(len(ind_save)): + if ind_save[ind] * ind_save[ind - 1] == 1: + connect_left[ind] = 1 + + if ind == (len(ind_save) - 1): + if ind_save[ind] * ind_save[0] == 1: + connect_right[ind] = 1 + else: + if ind_save[ind] * ind_save[ind + 1] == 1: + connect_right[ind] = 1 + + # Remove end connection if ring False + + if not ring: + connect_left[0] = 0 + + if not ring: + connect_right[-1] = 0 + + return connect_left, connect_right + diff --git a/chromo/marks.py b/chromo/marks.py new file mode 100644 index 0000000..25c583f --- /dev/null +++ b/chromo/marks.py @@ -0,0 +1,83 @@ +""" +Modifications to the polymer's chemical state that affect the energy. + +Each simulation type will typically require a particular type of Mark, for +which there will be a class defined here. In addition, for common, physically +derived marks, such as known epigenetic modifications to DNA or +well-characterized modifications to other real polymers, the best-known +parameters for that mark should be documented here, as instances of the +appropriate `.Mark` subclass. +""" + +import dataclasses +import inspect +import sys + +import pandas as pd + + +@dataclasses.dataclass +class Mark: + """In order for our code to work, all marks need a string name.""" + + name: str + + +@dataclasses.dataclass +class Epigenmark(Mark): + """Information about the chemical properties of an epigenetic mark.""" + + bind_energy: float + interaction_energy: float + chemical_potential: float + + +hp1 = Epigenmark('HP1', bind_energy=1, interaction_energy=1, + chemical_potential=1) +""" +Heterochromatin Protein 1, binds H3K9me marks. + +For now, we just use filler values for the energies. In the future, this +documentation string should go through the process of explaining exactly how we +arrive at the values that we actually use. +""" + + +def get_by_name(name): + """Look up saved mark by name.""" + all_marks = [obj for name, obj in inspect.getmembers(sys.modules[__name__]) + if isinstance(obj, Mark)] + matching_marks = [mark for mark in all_marks if mark.name == name] + if not matching_marks: + raise ValueError(f"No marks found in {__name__} with name: {name}") + if len(matching_marks) > 1: + raise ValueError(f"More than one mark has the name requested: {name}") + return matching_marks[0] + + +def make_mark_collection(marks): + """ + Construct summary DataFrame from sequence of marks. + + Parameters + ---------- + marks : Mark or str or Sequence[Mark] or Sequence[str] + The marks to be summarized by the DataFrame. + + Returns + ------- + pd.DataFrame + Columns are the properties of each Mark. + """ + df = pd.DataFrame(columns=['name', 'bind_energy', 'interaction_energy', + 'chemical_potential']) + if marks is None: + return None + input_type = type(marks) + if input_type is str or issubclass(input_type, Mark): + marks = [marks] # allow the "one mark" case + for mark in marks: + if type(mark) is str: + mark = get_by_name(mark) + df = df.append(dataclasses.asdict(mark), ignore_index=True) + return df diff --git a/chromo/mc/__init__.py b/chromo/mc/__init__.py index 407ac78..8f1af93 100644 --- a/chromo/mc/__init__.py +++ b/chromo/mc/__init__.py @@ -1,27 +1,51 @@ """Monte Carlo simulations of a discrete wormlike chain.""" from pathlib import Path +import warnings import numpy as np from .mc_sim import mc_sim from .moves import all_moves +from ..util.reproducibility import make_reproducible +from ..components import Polymer +from ..marks import get_by_name, make_mark_collection +from ..fields import UniformDensityField -def polymer_in_field(polymers, epigenmarks, field, num_mc_steps, num_save_mc, - mc_moves=None, output_dir='.'): +@make_reproducible +def simple_mc(num_polymers, num_beads, bead_length, num_marks, num_save_mc, + num_saves, x_width, nx, y_width, ny, z_width, nz, random_seed=0, + output_dir='.'): + polymers = [ + Polymer.straight_line_in_x( + f'Polymer-{i}', num_beads, bead_length, + states=np.zeros((num_beads, num_marks)), + mark_names=num_marks*['HP1'] + ) for i in range(num_polymers) + ] + marks = [get_by_name('HP1') for i in range(num_marks)] + marks = make_mark_collection(marks) + field = UniformDensityField(polymers, marks, x_width, nx, y_width, ny, + z_width, nz) + return _polymer_in_field(polymers, marks, field, num_save_mc, num_saves, + random_seed=random_seed, output_dir=output_dir) + + +def _polymer_in_field(polymers, marks, field, num_save_mc, num_saves, + mc_moves=None, random_seed=0, output_dir='.'): """ Monte Carlo simulation of a tssWLC in a field. This example code can be used to understand how to call the codebase, or - run directly for simple simulations. See the documentation for the + run directly for simple simulations. Parameters ---------- polymers : Sequence[Polymer] - The polymers to be simulated. Epigenetic marks, discretization, and - other parameters should be specified before input. - epigenmarks : Sequence[Epigenmark] - DEPRECATED. Should be integrated into the "Polymer" class. + The polymers to be simulated. + epigenmarks : `pd.DataFrame` + Output of `chromo.marks.make_mark_collection`. Summarizes the energetic + properties of each chemical modification. field : Field The discretization of space in which to simulate the polymers. num_save_mc : int @@ -31,11 +55,15 @@ def polymer_in_field(polymers, epigenmarks, field, num_mc_steps, num_save_mc, output_dir : Optional[Path], default: '.' Directory in which to save the simulation output. """ + warnings.warn("The random seed is currently ignored.", UserWarning) if mc_moves is None: mc_moves = all_moves # Perform Monte Carlo simulation for each save file - for mc_count in range(num_save_mc): - mc_sim(polymers, epigenmarks, num_mc_steps, mc_moves, field) + for mc_count in range(num_saves): + mc_sim(polymers, marks, num_save_mc, mc_moves, field) for poly in polymers: - poly.write_repr(output_dir / Path(f"{poly.name}-{mc_count}.npz")) + poly.to_csv(output_dir / Path(f"{poly.name}-{mc_count}.csv")) print("Save point " + str(mc_count) + " completed") + + +polymer_in_field = make_reproducible(_polymer_in_field) diff --git a/chromo/mc/mc_sim.py b/chromo/mc/mc_sim.py index 2dafdb7..e4a06b5 100644 --- a/chromo/mc/mc_sim.py +++ b/chromo/mc/mc_sim.py @@ -23,7 +23,7 @@ def mc_step(adaptible_move, poly, epigenmarks, field): if np.random.rand() < np.exp(-dE): adaptible_move.accept() poly.r[ind0:indf, :] = r - poly.t_3[ind0:indf, :] = t3 - poly.t_2[ind0:indf, :] = t2 + poly.t3[ind0:indf, :] = t3 + poly.t2[ind0:indf, :] = t2 poly.states[ind0:indf, :] = states return diff --git a/chromo/mc/moves.py b/chromo/mc/moves.py index 6290ef3..14154b7 100644 --- a/chromo/mc/moves.py +++ b/chromo/mc/moves.py @@ -7,6 +7,7 @@ import numpy as np + class MCAdapter: """ Track success rate and adjust parameters for a Monte Carlo move. @@ -95,7 +96,7 @@ def crank_shaft_move(polymer, amp_move, amp_bead): rot_angle = amp_move * (np.random.rand() - 0.5) if ind0 == (indf + 1): - delta_t3 = polymer.t_3[ind0, :] + delta_t3 = polymer.t3[ind0, :] else: delta_t3 = polymer.r[indf - 1, :] - polymer.r[ind0, :] delta_t3 /= np.linalg.norm(delta_t3) @@ -129,7 +130,8 @@ def crank_shaft_move(polymer, amp_move, amp_bead): t3_poly_trial = np.zeros((indf - ind0, 3), 'd') for i_bead in range(ind0, indf): r_poly_trial[i_bead - ind0, :] = rot_vector + np.matmul(rot_matrix, polymer.r[i_bead, :]) - t3_poly_trial[i_bead - ind0, :] = np.matmul(rot_matrix, polymer.t_3[i_bead, :]) + t3_poly_trial[i_bead - ind0, :] = np.matmul(rot_matrix, polymer.t3[i_bead, :]) return ind0, indf, r_poly_trial, t3_poly_trial, None, None all_moves = [MCAdapter(move) for move in [crank_shaft_move]] + diff --git a/chromo/util/reproducibility.py b/chromo/util/reproducibility.py index 14ed934..adb30ba 100644 --- a/chromo/util/reproducibility.py +++ b/chromo/util/reproducibility.py @@ -16,8 +16,10 @@ """ from pathlib import Path import inspect +import collections import pandas as pd +import numpy as np # get version number (includes commit hash) from versioneer from .._version import get_versions @@ -37,9 +39,11 @@ def make_reproducible(sim): f""" Decorator for automagically making simulations reproducible. - Requires every input to the function to have a __repr__ that is actually - useful, and that the wrapped function has a keyword argument with name - {output_dir_param_name}. Also requires you to pass a "random_seed" kwarg. + Requires every non-builtin input to the function to implement a + ``.to_file`` method and a ``.name`` attribute. These will be used to save + the state of the simulation's input to the output folder. The wrapped + function must have a keyword argument with name {output_dir_param_name}. + Also requires you to pass a "random_seed" kwarg. The versioneer version (including a commit hash if necessary), the date/time, a unique folder name where the simulation inputs and output will @@ -73,11 +77,11 @@ def make_reproducible(sim): ├── paramN └── (Any output of simulator goes here) """ - params = inspect.signature(sim) + params = inspect.signature(sim).parameters if "random_seed" not in params: raise ValueError("Requires the random seed to be passed in explicitly," - " even if the random generator is passed to the " - "function as well.") + " even if the random generator is passed to the" + " function as well.") if output_dir_param_name not in params: raise ValueError(f"Could not wrap simulation function: {sim}\n" f"Missing required kwarg: {output_dir_param_name}.") @@ -90,19 +94,28 @@ def wrapper(*args, **kwargs): "no default value can be found.") kwargs[output_dir_param_name] = outdir_param.default # get a new folder to save simulation into in a thread-safe way - outdir = kwargs[output_dir_param_name] - output_subfolder = get_unique_subfolder(Path(outdir)/sim_folder_prefix) + outdir = Path(kwargs[output_dir_param_name]) + # make the requested base output folder first if it doesn't exist tho.. + outdir.mkdir(parents=True, exist_ok=True) + output_subfolder = get_unique_subfolder(outdir / sim_folder_prefix) kwargs[output_dir_param_name] = output_subfolder # get the inputs that can be simply saved in our CSV file simple_params, hard_params = split_builtin_params(sim, *args, **kwargs) + # add the extra outputs that we promise users + simple_params['version'] = __version__ + simple_params['start_time'] = pd.Timestamp.now() # append to the CSV file, unless it doesn't exist - # TODO make thread-safe - sim_tracker = output_subfolder / sim_tracker_name + sim_tracker = outdir / sim_tracker_name + # TODO make writing header thread-safe header = not sim_tracker.exists() simple_params = pd.DataFrame(pd.Series(simple_params)).T simple_params.to_csv(sim_tracker, header=header, index=False, mode='a') - # TODO add support for non-simple params - return sim(*args, **kwargs) + # the remaining parameters should implement `.to_file`/`.from_file`. + to_file_params(hard_params, output_subfolder) + sim_out = sim(*args, **kwargs) + with open(output_subfolder / Path('__completion_time__'), 'w') as f: + f.write(str(pd.Timestamp.now())) + return sim_out return wrapper @@ -111,8 +124,15 @@ def wrapper(*args, **kwargs): # builtin_types = tuple(getattr(builtins, t) for t in dir(builtins) # if isinstance(getattr(builtins, t), type)) builtin_types = [bool, bytes, bytearray, complex, float, int, str] -# the additional types are not allowed by default because they make loading in -# the CSV a little trickier +""" +Types to be saved as CSV entries by default. + +These are basically just the built-ins that permit simple save/load cycling via +CSV entries with Pandas. Some extras are included in the split_builtin_params +function itself. Please see that code for a full description. +""" +# these additional types are not allowed by default because they make loading +# in the CSV a little trickier # ..., frozenset, tuple] @@ -121,7 +141,8 @@ def split_builtin_params(sim, *args, **kwargs): builtin_args = {} non_builtin = {} for arg_name, value in sig.arguments.items(): - if type(value) in builtin_types: + dtype = type(value) + if dtype in builtin_types or issubclass(dtype, Path): builtin_args[arg_name] = value else: non_builtin[arg_name] = value @@ -157,3 +178,36 @@ def get_unique_subfolder(root): except: i += 1 return folder_name + + +def to_file_params(non_builtins_kwargs, folder, suffix=''): + """ + Call ``.to_file`` for each parameter, handling sequences. + + Parameters + ---------- + non_builtins_kwargs : Dict[str, RoundTrippable] + A RoundTrippable must implement the + ``.name``/``.to_file``/``.from_file`` interface, or be a standard type, + such as a numpy array or a DataFrame. + folder : Path + The folder to save the output to. + """ + for arg_name, value in non_builtins_kwargs.items(): + dtype = type(value) + if isinstance(value, collections.Sequence): + for i, data in enumerate(value): + to_file_params({arg_name: data}, folder, str(i)) + elif hasattr(value, 'to_file'): + value.to_file(folder / Path(value.name + suffix)) + elif issubclass(dtype, pd.DataFrame) or issubclass(dtype, pd.Series): + value.to_csv(folder / Path(arg_name + suffix)) + elif issubclass(dtype, np.ndarray): + np.savetxt(folder / Path(arg_name + suffix), value) + else: + raise ValueError(f"Argument not understood: {arg_name}={value}.\n" + "This simulation cannot be made reproducible.\n" + f"Please implement `.to_file` for type: {dtype}.") + + + diff --git a/doc/Makefile b/doc/Makefile new file mode 100644 index 0000000..d0c3cbf --- /dev/null +++ b/doc/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = source +BUILDDIR = build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/doc/make.bat b/doc/make.bat new file mode 100644 index 0000000..6247f7e --- /dev/null +++ b/doc/make.bat @@ -0,0 +1,35 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=source +set BUILDDIR=build + +if "%1" == "" goto help + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.http://sphinx-doc.org/ + exit /b 1 +) + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% + +:end +popd diff --git a/doc/source/bd_sim.rst b/doc/source/bd_sim.rst new file mode 100644 index 0000000..01a69ab --- /dev/null +++ b/doc/source/bd_sim.rst @@ -0,0 +1,169 @@ +.. _bd_sim: + +Brownian Dynamics Simulations +============================= + +Langevin equation for translation +--------------------------------- + +We consider a discrete chain of :math:`n_{b}` beads with positions :math:`\vec{r}^{(n)}`, +where the bead index :math:`n` runs from 0 to :math:`n_{b}-1`. +We consider a general potential energy :math:`E = E( \{ \vec{r} \} )`, where :math:`\{ \vec{r} \}` indicates +the full set of positions. +We define the potential forces on the beads to be + +.. math:: + \vec{f}_{E}^{(n)} = - \frac{\partial E}{\partial \vec{r}^{(n)}} + +The Langevin equation for the bead positions is given by + +.. math:: + \xi_{r} \frac{\partial \vec{r}^{(n)}}{\partial t} = \vec{f}_{E}^{(n)} + \vec{f}_{B}^{(n)} + +The Brownian force :math:`\vec{f}_{B}^{(n)}` is governed by the fluctuation dissipation theorem, +which states that :math:`\vec{f}_{B}^{(n)}` is a +Gaussian-distributed random force +that satisfies + +.. math:: + \langle \vec{f}_{B}^{(n)}(t) \rangle & = & 0, \\ + \langle \vec{f}_{B}^{(n)}(t) \vec{f}_{B}^{(n')}(t') \rangle + & = & + 2 k_{B} T \xi \delta (t - t') \delta_{nn'} \mathbf{I}. + +This stochastic differential equation is numerically integrated by +discretizing time in discrete time steps :math:`\delta t`. +Thus, the discrete time at the :math:`j`th time step is +given by :math:`t_{j} = \delta t (j-1)`, where :math:`t_{0} = 0`. +Based on this time discretization, we numerically integrate the Langevin +equation as + +.. math:: + \vec{r}^{(n)}_{j+1} = + \vec{r}^{(n)}_{j} + + \frac{\delta t}{\xi_{r}} \left( + \vec{f}_{E_{j}}^{(n)} + \vec{f}_{B_{j}}^{(n)} + \right) + +which suggests an explicit integration algorithm. + +The fluctuation dissipation for the discrete time + +.. math:: + \langle \vec{f}_{B_{j}}^{(n)} \rangle & = & 0, \\ + \langle \vec{f}_{B_{j}}^{(n)} \vec{f}_{B_{j'}}^{(n')} \rangle + & = & + \frac{2 k_{B} T \xi}{\delta t} \delta_{jj'} \delta_{nn'} \mathbf{I}. + +Thus, the value of :math:`\vec{f}_{B_{j}}^{(n)}` is selected from a Gaussian distribution +with variance :math:`\frac{2 k_{B} T \xi}{\delta t}`. + +Langevin equation for translation and rotation +---------------------------------------------- + +We consider a discrete chain of :math:`n_{b}` beads with positions :math:`\vec{r}^{(n)}` and orientation triad :math:`\vec{t}_{i}^{(n)}`, +where the bead index :math:`n` runs from 0 to :math:`n_{b}-1`. +We consider a general potential energy :math:`E = E( \{ \vec{r}, \vec{t}_{i} \} )`, where :math:`\{ \vec{r}, \vec{t}_{i} \}` indicates +the full set of positions and orientations. +Since the orientation triad forms an orthonormal basis, the dynamics must maintain the following conditions: + +.. math:: + \vec{t}_{i}^{(n)} \cdot \vec{t}_{j}^{(n)} = \delta_{ij} + +.. \label{eq:constraint} + +for all :math:`i,j` pairs. +These conditions are enforced through Lagrange constraints. Thus, we define the constrained energy + +.. math:: + E_{\lambda} ( \{ \vec{r}, \vec{t}_{i} \} ) = E( \{ \vec{r}, \vec{t}_{i} \} ) - + \sum_{n=0}^{n_{b}-1} \sum_{i,j = 1}^{3} + \frac{\lambda_{ij}}{2} + \left( + \vec{t}_{i}^{(n)} \cdot \vec{t}_{j}^{(n)} - \delta_{ij} + \right) + +We define the potential forces and torques on the beads to be + +.. math:: + \vec{f}_{E}^{(n)} & = & - \frac{\partial E}{\partial \vec{r}^{(n)}} \\ + \vec{\tau}_{E,i}^{(n)} & = & - \frac{\partial E}{\partial \vec{t}_{i}^{(n)}} + +The Langevin equation for the bead positions is given by + +.. math:: + \xi_{r} \frac{\partial \vec{r}^{(n)}}{\partial t} = \vec{f}_{E}^{(n)} + \vec{f}_{B}^{(n)} + +where :math:`\vec{f}_{B}^{(n)}` is a Brownian force (discussed below). + +The Langevin equation for the orientation triad must resolve the orthonormal constraints. This leads to the three Langevin equations + +.. math:: + \xi_{t,i} \frac{\partial \vec{t}_{i}^{(n)}}{\partial t} = \vec{\tau}_{E,i}^{(n)} + \vec{\tau}_{B,i}^{(n)} + + \sum_{j=1}^{3} \lambda_{ij} \vec{t}_{j}^{(n)} + +The Lagrange constraints are satisfied by setting the time derivative of Eq.~\ref{eq:constraint} = 0, +or we have + +.. math:: + \vec{t}_{i}^{(n)} \cdot \frac{\partial \vec{t}_{j}^{(n)}}{\partial t} + + \vec{t}_{j}^{(n)} \cdot \frac{\partial \vec{t}_{i}^{(n)}}{\partial t} = 0 + +This leads to the solution to the Lagrange multipliers to be + +.. math:: + \lambda_{ij} = - \left( \xi_{t,i} + \xi_{t,j} \right)^{-1} + \left[ + \xi_{t,i} \vec{t}_{i}^{(n)} \cdot \left( \vec{\tau}_{E,j}^{(n)} + \vec{\tau}_{B,j}^{(n)} \right) + + \xi_{t,j} \vec{t}_{j}^{(n)} \cdot \left( \vec{\tau}_{E,i}^{(n)} + \vec{\tau}_{B,i}^{(n)} \right) + \right] + +With this development, we now write the Langevin equations as + +.. math:: + \xi_{t,1} \frac{\partial \vec{t}_{1}^{(n)}}{\partial t} & = & + \frac{\xi_{t,1}}{\xi_{t,1}+\xi_{t,2}} + \left[ + \vec{t}_{2}^{(n)} \cdot \left( \vec{\tau}_{E,1}^{(n)} + \vec{\tau}_{B,1}^{(n)} \right) - + \vec{t}_{1}^{(n)} \cdot \left( \vec{\tau}_{E,2}^{(n)} + \vec{\tau}_{B,2}^{(n)} \right) + \right] \vec{t}_{2}^{(n)} + + \nonumber \\ + & & + \frac{\xi_{t,1}}{\xi_{t,1}+\xi_{t,3}} + \left[ + \vec{t}_{3}^{(n)} \cdot \left( \vec{\tau}_{E,1}^{(n)} + \vec{\tau}_{B,1}^{(n)} \right) - + \vec{t}_{1}^{(n)} \cdot \left( \vec{\tau}_{E,3}^{(n)} + \vec{\tau}_{B,3}^{(n)} \right) + \right] \vec{t}_{3}^{(n)} \\ + \xi_{t,2} \frac{\partial \vec{t}_{2}^{(n)}}{\partial t} & = & + \frac{\xi_{t,2}}{\xi_{t,1}+\xi_{t,2}} + \left[ + \vec{t}_{1}^{(n)} \cdot \left( \vec{\tau}_{E,2}^{(n)} + \vec{\tau}_{B,2}^{(n)} \right) - + \vec{t}_{2}^{(n)} \cdot \left( \vec{\tau}_{E,1}^{(n)} + \vec{\tau}_{B,1}^{(n)} \right) + \right] \vec{t}_{1}^{(n)} + + \nonumber \\ + & & + \frac{\xi_{t,2}}{\xi_{t,2}+\xi_{t,3}} + \left[ + \vec{t}_{3}^{(n)} \cdot \left( \vec{\tau}_{E,2}^{(n)} + \vec{\tau}_{B,2}^{(n)} \right) - + \vec{t}_{2}^{(n)} \cdot \left( \vec{\tau}_{E,3}^{(n)} + \vec{\tau}_{B,3}^{(n)} \right) + \right] \vec{t}_{3}^{(n)} \\ + \xi_{t,3} \frac{\partial \vec{t}_{3}^{(n)}}{\partial t} & = & + \frac{\xi_{t,3}}{\xi_{t,1}+\xi_{t,3}} + \left[ + \vec{t}_{1}^{(n)} \cdot \left( \vec{\tau}_{E,3}^{(n)} + \vec{\tau}_{B,3}^{(n)} \right) - + \vec{t}_{3}^{(n)} \cdot \left( \vec{\tau}_{E,1}^{(n)} + \vec{\tau}_{B,1}^{(n)} \right) + \right] \vec{t}_{1}^{(n)} + + \nonumber \\ + & & + \frac{\xi_{t,3}}{\xi_{t,2}+\xi_{t,3}} + \left[ + \vec{t}_{2}^{(n)} \cdot \left( \vec{\tau}_{E,3}^{(n)} + \vec{\tau}_{B,3}^{(n)} \right) - + \vec{t}_{3}^{(n)} \cdot \left( \vec{\tau}_{E,2}^{(n)} + \vec{\tau}_{B,2}^{(n)} \right) + \right] \vec{t}_{2}^{(n)} + + +Brownian dynamics simulations with constrained bond lengths +----------------------------------------------------------- + +Overview of bond constraints +Ref [Hinch1994]_ \ No newline at end of file diff --git a/doc/source/code_structure.rst b/doc/source/code_structure.rst new file mode 100644 index 0000000..fc3c3c9 --- /dev/null +++ b/doc/source/code_structure.rst @@ -0,0 +1,9 @@ +.. _code_structure: + +Code Overview and Structure +=========================== + +We provide an overview of the 'chromo' code structure and the individual classes and modules that +are responsible for the individual simulation capabilities. + +9/18/20 - Bruno will update this page with a code map and discussion \ No newline at end of file diff --git a/doc/source/conf.py b/doc/source/conf.py new file mode 100644 index 0000000..59bb2d3 --- /dev/null +++ b/doc/source/conf.py @@ -0,0 +1,71 @@ +# Configuration file for the Sphinx documentation builder. +# +# This file only contains a selection of the most common options. For a full +# list see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +# -- Path setup -------------------------------------------------------------- + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +import os +import sys +sys.path.insert(0, os.path.abspath('../')) +sys.path.insert(1, os.path.abspath('../../')) + +# -- Project information ----------------------------------------------------- + +project = 'chromo' +copyright = '2020, Andy Spakowitz, Bruno Beltran, Joe Wakim' +author = 'Andy Spakowitz, Bruno Beltran, Joe Wakim' + + +# -- General configuration --------------------------------------------------- + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# ones. +extensions = [ + 'sphinx.ext.autodoc', + 'sphinx.ext.autosummary', + 'nbsphinx', + 'sphinx.ext.doctest', + 'sphinx.ext.todo', + 'sphinx.ext.mathjax', + 'sphinx.ext.ifconfig', + 'sphinx.ext.viewcode', + 'sphinx.ext.napoleon', + 'sphinx.ext.intersphinx', + 'sphinx.ext.githubpages', + 'matplotlib.sphinxext.plot_directive' +] + +# Add any paths that contain templates here, relative to this directory. +templates_path = ['_templates'] + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +# This pattern also affects html_static_path and html_extra_path. +exclude_patterns = [] + + +# -- Options for HTML output ------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# +html_theme = 'nature' + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ['_static'] + +# Autodoc settings +autosummary_generate = True + +# Napoleon settings +napoleon_google_docstring = False +napoleon_numpy_docstring = True diff --git a/doc/source/field.rst b/doc/source/field.rst new file mode 100644 index 0000000..aa682da --- /dev/null +++ b/doc/source/field.rst @@ -0,0 +1,69 @@ +.. _field: + + +Field-Theoretic Treatment of Interactions +========================================= + +Our treatment of interactions uses a field-theoretic treatment of the densities to determine the interactions between polymer +segments. Following work by Pike, et al. (Refs. [Pike2009a]_, [Pike2009b]_), +we define + +The simulation has a fixed volume with sides lengths :math:`L_{x}`, :math:`L_{y}`, and :math:`L_{z}`. +These lengths are discretize into :math:`M_{x}`, :math:`M_{y}`, and :math:`M_{z}` bins of length +:math:`\Delta_{x} = L_{x}/M_{x}`, +:math:`\Delta_{y} = L_{y}/M_{y}`, and +:math:`\Delta_{z} = L_{z}/M_{z}`. +The bins are defined by the three indices +:math:`i_{x}`, +:math:`i_{y}`, and +:math:`i_{z}` that run from zero to +:math:`M_{x}-1`, +:math:`M_{y}-1`, and +:math:`M_{z}-1`, respectively. + + +We consider the :math:`n`th bead located at position :math:`\vec{r}^{(n)}`. +We define a weight function :math:`w_{I}(\vec{r}^{(n)})` within the :math:`I`th bin. +The :math:`I`th index is defined to be a superindex that combines +:math:`i_{x}`, +:math:`i_{y}`, and +:math:`i_{z}` into a single unique index :math:`I= i_{x} + M_{x} i_{y} + M_{x}M_{z} i_{z}` that +runs from zero to :math:`M_{x}M_{y}M_{z}-1` (total of :math:`M_{x}M_{y}M_{z}` unique indices) +The total weight on the :math:`I`th bin is given by the contributions from the three cartesian +directions, \emph{i.e.} +:math:`w_{I}(\vec{r}^{(n)}) = +w_{i_{x}}^{(x)}(x^{(n)}) +w_{i_{y}}^{(y)}(y^{(n)}) +w_{i_{z}}^{(z)}(z^{(n)})`. +Figure~\ref{fig:weight} shows a schematic of the :math:`x`-direction weight function (same method for :math:`y` and :math:`z`). +This shows a linear interpolation weighting method, consistent with Refs. [Pike2009a]_, [Pike2009b]_. + +.. figure:: figures/weight.pdf + :width: 600 + :align: center + :alt: Schematic of the weight function :math:`w_{i_{x}}^{(x)}` that gives the weighting of the particle in the :math:`i_{x}` site in the + :math:`x`-direction based on a linear interpolation method + + Schematic of the weight function :math:`w_{i_{x}}^{(x)}` that gives the weighting of the particle in the :math:`i_{x}` site in the + :math:`x`-direction based on a linear interpolation method + + +The number of epigenetic proteins (\emph{e.g.} HP1) to the :math:`n`th site is given by :math:`N_{I}^{(\alpha)}`, where :math:`\alpha` determines +the type of epigenetic mark. +The :math:`\alpha`-protein density within the :math:`I`th bin is given by + +.. math:: + \rho_{I}^{(\alpha)} = \frac{1}{v_{\mathrm{bin}}} \sum_{n=0}^{n_{b} - 1} w_{I}(\vec{r}^{(n)}) N_{I}^{(\alpha)} + +where :math:`v_{\mathrm{bin}} = \Delta_{x} \Delta_{y} \Delta_{z}` is the volume of a bin. +The maximum number of epigenetic proteins bound :math:`N_{\mathrm{max}}^{(\alpha)}` gives an upper bound on the +number of proteins that can bind to a bead, accounting for coarse graining of a bead to represent multiple nucleosomes. +For discretization of one nucleosome per bead, the maximum :math:`N_{\mathrm{max}}^{(\alpha)} = 2` implies binding +of a protein to the two histone tail proteins for the :math:`\alpha` epigenetic mark. +We define the number of :math:`\alpha` marks on the :math:`I`th bead as :math:`M_{I}^{(\alpha)}`, which can take values from zero +to :math:`N_{\mathrm{max}}^{(\alpha)}`. + +Protein binding to a marked tail results in energy :math:`-\beta \epsilon_{m}` [non-dimensionalized by :math:`\beta = 1/(k_{B}T)`], and protein binding to an unmarked tail is associated with +energy :math:`-\beta \epsilon_{u}`. The chemical potential of the :math:`\alpha` protein is defined as :math:`\beta \mu^{(\alpha)}`. +The binding of :math:`N_{I}^{(\alpha)}` proteins to a bead with :math:`M_{I}^{(\alpha)}` marks results in a free energy that +accounts for all of the combinatoric ways of binding. diff --git a/doc/source/figures/weight.ai b/doc/source/figures/weight.ai new file mode 100644 index 0000000..840f5d3 --- /dev/null +++ b/doc/source/figures/weight.ai @@ -0,0 +1,1119 @@ +%PDF-1.5 % +1 0 obj <>/OCGs[23 0 R]>>/Pages 3 0 R/Type/Catalog>> endobj 2 0 obj <>stream + + + + + application/pdf + + + weight + + + Adobe Illustrator 24.1 (Macintosh) + 2020-08-02T19:05:18-07:00 + 2020-08-02T19:05:18-07:00 + 2020-08-02T19:05:18-07:00 + + + + 256 + 52 + JPEG + /9j/4AAQSkZJRgABAgEASABIAAD/7QAsUGhvdG9zaG9wIDMuMAA4QklNA+0AAAAAABAASAAAAAEA AQBIAAAAAQAB/+4ADkFkb2JlAGTAAAAAAf/bAIQABgQEBAUEBgUFBgkGBQYJCwgGBggLDAoKCwoK DBAMDAwMDAwQDA4PEA8ODBMTFBQTExwbGxscHx8fHx8fHx8fHwEHBwcNDA0YEBAYGhURFRofHx8f Hx8fHx8fHx8fHx8fHx8fHx8fHx8fHx8fHx8fHx8fHx8fHx8fHx8fHx8fHx8f/8AAEQgANAEAAwER AAIRAQMRAf/EAaIAAAAHAQEBAQEAAAAAAAAAAAQFAwIGAQAHCAkKCwEAAgIDAQEBAQEAAAAAAAAA AQACAwQFBgcICQoLEAACAQMDAgQCBgcDBAIGAnMBAgMRBAAFIRIxQVEGE2EicYEUMpGhBxWxQiPB UtHhMxZi8CRygvElQzRTkqKyY3PCNUQnk6OzNhdUZHTD0uIIJoMJChgZhJRFRqS0VtNVKBry4/PE 1OT0ZXWFlaW1xdXl9WZ2hpamtsbW5vY3R1dnd4eXp7fH1+f3OEhYaHiImKi4yNjo+Ck5SVlpeYmZ qbnJ2en5KjpKWmp6ipqqusra6voRAAICAQIDBQUEBQYECAMDbQEAAhEDBCESMUEFURNhIgZxgZEy obHwFMHR4SNCFVJicvEzJDRDghaSUyWiY7LCB3PSNeJEgxdUkwgJChgZJjZFGidkdFU38qOzwygp 0+PzhJSktMTU5PRldYWVpbXF1eX1RlZmdoaWprbG1ub2R1dnd4eXp7fH1+f3OEhYaHiImKi4yNjo +DlJWWl5iZmpucnZ6fkqOkpaanqKmqq6ytrq+v/aAAwDAQACEQMRAD8A9U4q7FXYqlv+ItLXV20m aQwXuxhWZSizAgH9yx+F6VoQN8VTLFXYq7FXYq7FWJeY/wA0fKfl9i169xLaRSCO/v7S3kubazqD RrmSINwWooTvQ9aYqyXT9QsNRsob7T7iO7srhQ8FzA6yRup6FWUkEYqiMVY/5284W3lPTbbVL2Ll pz3cNrfXBfgLeOeqLLTi3L97wSlR9qtdqFV5Za/mx5i8/aro+l+WFfQdXNrNcXCyyGWOC8jmMM0V wvpUkSCKM1VgOTTR04kc1Vep+SIdfh0iZdbMhuTd3DQ+tIJZPSaQkVcKnwluRQcVohUcVpxCq3zd 5/8ALHlRYE1S5LX94wTT9KtlM97cuTQLDbpV23/a+yO5xVMdH8w6XqyuLWQrPFtcWkoMc8R8Hjb4 h+rFUxxV2KuxV2KqdxcQW0Lz3EiwwxjlJI5Cqo8STiqSWfnjy/dXSQCSSBJzS0uZ42ignPQ+lI4A b+PbFU/xV2KuxV2KuxVLf8RaW2rrpMMhnvdzMsKl1hABP75h8KVpQA74qmWKuxV2KuxV2KuxVCan pOnapata39ulxA37LDcHxVhup9xiqR+j5k8v7wGTXNHXrA5rfQr/AJDGgmA8D8WKpzpGt6Zq1uZr CcShTxlj3WSNv5XQ0ZT88VR2KpbrPmHTNJRPrUha4l2t7SIGSeU+EcY3Pz6YqlX6N8weYPi1d20v Sm6aXbv+/lU9riZfsg90T6Tiqf2enWFnaLZ2tvHDaqOIhRQFoetR3r3xVgGoflbf6Fey6z+Wl9Ho V5KxkvNAnDPo14x68oV3tnP+/Iae64qjfLP5q2N5qieXfM9lJ5W82GoTTL1gYbqlKvY3QpFcLv0F G/ycVZV5g8v6N5h0e50bWbVb3TLsBbm2ckK4Vg67qVOzKDscVed+f9C03yZcHzt5eksNP1oG9P1P UJFit7271J7P1pHmmuIBH6dvYE8VP+UFJ+F1XW3nzzz58t44PIdp+i9KdFF35z1OArGWKjl+jrIs 5nYMdnd/TFOrbHFWU+Tfy18u+WJpdQj9XU/MN0KX/mDUH9e9m8R6jf3aeCRgL02xVNtZ8tadqjpc NzttQiH7jULc+nOntyH2l/yWqMVS0a5rOhEReYo/rFiNl1u2Q8QP+XmEVMf+stV+WKsjt7m3uYEn t5FmgkHKOVCGVh4gjbFVTFUi1PzVBDdHTtLhbVNWH2rWEgJFXvPKfhjH4+2KqNv5WuL+dL3zNOt/ Mp5Q6fGCtnCfZD/et/lP92Kp5e6fY3tq1pdwJPbOKNE6grt029sVY/8Ao3zB5f8Ai0h21TSl66Xc P+/iUdreZvtAdkf6Diqa6N5h0zVkf6rIVuItri0lBjniPhJGdx8+mKpliqB1fW9M0m3E1/OIgx4x R7tJI38qIKsx+WKpN6PmTzBvOZND0dukCGl9Mv8AlsKiEHwHxYqnmmaTp2l2q2thbpbwL+yo3J8W Y7sfc4qi8VdirsVdirsVdirsVdiqRa/5f0uUtq31g6TqEC1/SsTKhAHaUN8EiezYqxf/ABx5ilto 4ZBFZ2ckhh/xQYpTbsv7LpEyjizeLHhXFWXaH5e0rT63kLG8vLgBpdSmb1ZZAelH6BfALtiqb4q7 FXYqlPmfyp5d80aVJpWv2EWoWMm/pyjdWHR43FHjcdmUg4q8bvfPPmzyD5gl8seVbiX8zLeGJ5H0 dmkk1PS+IJRZr2KOWORDSgSX970Ar3VTf8v9HsfzEB80eY/Mg1u+Cm3uPL1iDa2VjG0iyvZTwSKL lyZIkLmWnLjSnHqqzDWvzW8p6DriaFqAmt7k3trp6HjGIh9chaWKavOohBQxk0rz2ApvirFY/Pv5 jarb3Or6PYyJodwYp7H17eMXCCaKOSzgRfUq8d1yX1X4sYy5Cn7RhVeu4q4gEEEVB2IOKsN1+ytP LYbU9Ivk0yaVqnTHBktrp/5VhWrK58Y8VQkOuatr+oJpmqyN5ZjdFZbMFlubqoBISZlVVX2X4sVZ lpmk6dpdqtpYQLbwLvxUbk+LE7sfc4qi8VdirsVSjXPL2lahS8mY2d5bgtFqULelLGB1q/Qr4hts VYj/AI48xRW0sMYivLOKQRf4oEUot1X9pniVTyZfFW4YqyjQNA0uIrqv1k6tfzrX9KSsHJB7RBfg jT2XFU9xV2KuxV2KuxVgH5A/+Sa8p/8AMCv/ABJsVZ/irsVdiqRal5qgium07S4G1TVh9q2hICRe 88p+GMfj7Yqo2/la4vp0vfMs6386HlDYICtnCfZDvIw/mf7sVSCGKJ/z51KJ0Vom8p2SshAKlTqN 2KEeGKp7J5b1DSHa48szLFETyl0eck2z16+k32oWPt8PtiqM0jzRZ31wbG5jfT9WQVk0+4oHP+VG 32ZF26riqc4qkXm7zv5Z8pWC3uuXq24lbha26gyXFxIekcEKVkkck9FHz2xVhv1L8yPzBqdQa48j +T3+zYQsBrd4h/3/ACjktmjD9hKydQSMVUfzI8qeXfK/5eabpWgWEWn2Meu6KfTiG7N+kYAXkc1e RzTdmJJxVPPOv5T6R5h1OPXNPmGh+ZYkaNdatoInnpVWRiWAPNGjADV+wWX9rZV5t5t0ixbVrj/l cVlJaG+sF0WHztpJ/wBxrxC5FxDJPGySvZXKyCgd6pu1DTFWWaFJ5t0fzDqOjaHpsSaX5k1KS+0j zAUlvbC2tItNtVj5xWxjUiZ7eRErcR0O56osirJ9C0C90zz95o126jjisNTsNJX64GVVkuLT62Lh ihZmQKskf2u3c0OKomTzJqGru1v5YhWSIHjLrE4Itkp19JdmmYe3w++Ko3SPK9nY3BvrmR9Q1ZxS TULihcf5Ma/ZjXf7K4qw78pNK07VPI2oWt/AlxA2ua18LDcH9JT/ABKRup9xvirI/R8yeXt7cya5 o69bdyDfQr/xWxoJgPBvixVOdI1vTNXtzNYzCQKeMsZqskbd1kQ0ZT88VR2KpNq/mizsbgWNtG+o as4rHp9vQuP8qRvsxrv1bFUHH5b1DV3W48zTLLEDyi0eAkWyU6eq32pmHv8AD7YqkWrxRJ+dnlaJ EVYl0HVlVAAFCiezFAPDFU/uPK1xYzve+Wp1sJ3PKawcFrOY+6DeNj/Mn3YqraZ5qhlul07VIG0v Vj9m2mIKS+8Eo+GQfj7YqnuKuxV2KsA/P7/yTXmz/mBb/iS4qwv8mfMf5pW/5W+W4NL8mWmoafHa KLa8k1lbd5E5N8RhNrJw+XI4qzT/ABX+cv8A1IFj/wBx5f8AsixV3+K/zl/6kCx/7jy/9kWKvM/K nnL/AJyY8y3uqWer+VfQsbaWjxmQ6OxBJ/dRXMkdz6qUH2kXp+1uMVeg6Zq35qaXara2H5c6dbwL vxTXl3Pix+pVY+5xVF/4r/OX/qQLH/uPL/2RYqwuDzH+aQ/Oa+uB5MtDqh8uWkb2H6ZUItuL65Kz ev8AVdyzll4cNqVrvTFWaf4r/OX/AKkCx/7jy/8AZFiqA1e//NHV7cQX35c2EgU1ikGvqskbfzRu LLkp+WKvPfNXnj/nKDRta0vQ9L8uI0d6CkFw1dUAZyUAuLyNLdE9IUeroPE1GKp95Q8s/mXoV82t 6h5HtvMHmyZaXHmHUdfjeen8lugsvTtotzRIwNtiTirM/wDFf5y/9SBY/wDceX/sixVhf5s+Y/zS n8s2iaj5MtLGAavpTpMmsrOTKt/E0UfAWsdBI4Clq/DWtD0xVkNx5t8ySX1zF5xvR5AeG3STR7a1 nt75byRi/qMJZrceu6cVX6tGobep5ck4qomLzP8AnDeaYkd5+XlhMlxCBcRTayiBw6/ErwtaSca1 3QsadKnFXm2r6V+efkmC81fyF5WTR9MKs935bi1Aazbc3NPVs7T0IJInBbkVjcqf5NsVTnyrqv56 +ZdGstS8z+TYb5GWsdhPqH6JRirECSe0a3nkqaVAdqf5NDirOI/M/wCcMaLHH+X1gkaAKiLryAAD oABZYqu/xX+cv/UgWP8A3Hl/7IsVYX+U3mP80oPLN2mneTLS+gOr6q7zPrKwEStfytLHwNrJURuS oavxUrQdMVZp/iv85f8AqQLH/uPL/wBkWKpLq7fmxfTi+t/Ilpp+rIKR6hb6+iyf6si/UuMi+zYq wzQvPX/OTWs+Z9Y0O88srDaWhIaQE6eAqtwAgv5EnjkMi/FyRD4jiMVZ9pF/+aOkW5gsfy5sIwxr LIdfVpJG/mkc2XJj88VR/wDiv85f+pAsf+48v/ZFirC9V8x/mk35ueX7iTyZaJqSaRqSW9iNZVkk iaa1Mkhm+qjgUIUBeJrXqKYqzT/Ff5y/9SBY/wDceX/sixVCanq35qapata3/wCXOnXEDb8X15dj 4qfqVVPuMVefebfOP/OS/lq50u10fyt69lcTUWMSHWWABUCKa5jjtvSQ1PxOOn7W2KvSx5r/ADlp /wAoBYj2/T6f9kWKt/4r/OX/AKkCx/7jy/8AZFirC/zm8x/mlcflb5kg1TyZaafp8lowubyPWVuH jTkvxCEWsfP5chirNPyB/wDJNeU/+YFf+JNirP8AFXYq7FXYq7FWAW3/AJP3Uf8AwFLH/uo3eKs/ xVD6lJLFp11JEeMqQyNGwAJDBSQaGo64qwr8vPzG/wAR6FokKQXWo6y2n2EvmK6iijgis7i7tEnD Sidrfn6nIsot1egoSArKSqmX5YarrGqeT4rrWLk3moJe6nbSXLRpEXS01G4toiUjVEB9OJRsMVZX irAPzu/5RCx/7bmi/wDdSgxVn+KuxVjP5jHWx5XP6FFwb83+mClrz9X0DqNuLn+7+Lj6HPn/AJNa 7YqifJXnLTPN+gx61pscsNrI7xiOf0vUBjNDUQyTLQ9VPL4hRhVSCVU9xV2KsA/JH/lEL7/tua1/ 3Up8VZ/irsVdirsVdirANa/8nj5X/wC2Hq3/ACfs8VZ/irsVdirsVdirAPz+/wDJNebP+YFv+JLi rI/JH+Fv8J6Z/hTh/h30R+jPT9Tj6VTSnq/H1r9rFU8xV2KuxV2KuxVI4/8AC3+N5+HD/Ff6Mh9f +85/o/6xL6X/ABXT1vU/yvoxVPMVdiqU6V/hT9Lan+ifqH6X/d/pj6r6P1mvKT0vrPp/H9r1ePPv yp3xVNsVdiqR+cP8LfouH/E3D9H/AFy09Hn6lPrf1hPqv93vX1+NO3jtiqeYq7FXYq7FXYq7FUj8 n/4W/Rc3+GeH6P8Arl363D1KfW/rD/Wv7zevr8q9vDbFU8xV2KuxV2KuxVI7v/C3+MdO+s8P8T/U 7r9HV9Tn9U5xfWaU/d05+nWu/h3xVPMVdirsVdirsVSPzv8A4W/wnqf+K+H+HfRP6T9T1OPpVFa+ l8fWn2cVf//Z + + + + uuid:9E3E5C9A8C81DB118734DB58FDDE4BA7 + xmp.did:fe646f74-ee7e-4b6e-92b8-1740a48566ef + uuid:b5a937fd-0226-ef46-b8bd-32e27ebcdfcb + proof:pdf + + uuid:eab21be5-e766-664f-af13-7ab4a523d93d + xmp.did:cc4dc041-9a60-4315-bc5d-7e2d77be13d2 + uuid:9E3E5C9A8C81DB118734DB58FDDE4BA7 + proof:pdf + + + + + saved + xmp.iid:fe646f74-ee7e-4b6e-92b8-1740a48566ef + 2020-08-02T16:40:33-07:00 + Adobe Illustrator 24.1 (Macintosh) + / + + + + Basic RGB + Document + AIRobin + 1 + False + False + + 6.000000 + 3.000000 + Inches + + + + + ChalkboardSE-Regular + Chalkboard SE + Regular + TrueType + 13.0d1e2 + False + ChalkboardSE.ttc + + + MyriadPro-Regular + Myriad Pro + Regular + Open Type + Version 2.106;PS 2.000;hotconv 1.0.70;makeotf.lib2.5.58329 + False + MyriadPro-Regular.otf + + + + + + Cyan + Magenta + Yellow + Black + + + + + + Default Swatch Group + 0 + + + + White + RGB + PROCESS + 255 + 255 + 255 + + + Black + RGB + PROCESS + 0 + 0 + 0 + + + RGB Red + RGB + PROCESS + 255 + 0 + 0 + + + RGB Yellow + RGB + PROCESS + 255 + 255 + 0 + + + RGB Green + RGB + PROCESS + 0 + 255 + 0 + + + RGB Cyan + RGB + PROCESS + 0 + 255 + 255 + + + RGB Blue + RGB + PROCESS + 0 + 0 + 255 + + + RGB Magenta + RGB + PROCESS + 255 + 0 + 255 + + + R=193 G=39 B=45 + RGB + PROCESS + 193 + 39 + 45 + + + R=237 G=28 B=36 + RGB + PROCESS + 237 + 28 + 36 + + + R=241 G=90 B=36 + RGB + PROCESS + 241 + 90 + 36 + + + R=247 G=147 B=30 + RGB + PROCESS + 247 + 147 + 30 + + + R=251 G=176 B=59 + RGB + PROCESS + 251 + 176 + 59 + + + R=252 G=238 B=33 + RGB + PROCESS + 252 + 238 + 33 + + + R=217 G=224 B=33 + RGB + PROCESS + 217 + 224 + 33 + + + R=140 G=198 B=63 + RGB + PROCESS + 140 + 198 + 63 + + + R=57 G=181 B=74 + RGB + PROCESS + 57 + 181 + 74 + + + R=0 G=146 B=69 + RGB + PROCESS + 0 + 146 + 69 + + + R=0 G=104 B=55 + RGB + PROCESS + 0 + 104 + 55 + + + R=34 G=181 B=115 + RGB + PROCESS + 34 + 181 + 115 + + + R=0 G=169 B=157 + RGB + PROCESS + 0 + 169 + 157 + + + R=41 G=171 B=226 + RGB + PROCESS + 41 + 171 + 226 + + + R=0 G=113 B=188 + RGB + PROCESS + 0 + 113 + 188 + + + R=46 G=49 B=146 + RGB + PROCESS + 46 + 49 + 146 + + + R=27 G=20 B=100 + RGB + PROCESS + 27 + 20 + 100 + + + R=102 G=45 B=145 + RGB + PROCESS + 102 + 45 + 145 + + + R=147 G=39 B=143 + RGB + PROCESS + 147 + 39 + 143 + + + R=158 G=0 B=93 + RGB + PROCESS + 158 + 0 + 93 + + + R=212 G=20 B=90 + RGB + PROCESS + 212 + 20 + 90 + + + R=237 G=30 B=121 + RGB + PROCESS + 237 + 30 + 121 + + + R=199 G=178 B=153 + RGB + PROCESS + 199 + 178 + 153 + + + R=153 G=134 B=117 + RGB + PROCESS + 153 + 134 + 117 + + + R=115 G=99 B=87 + RGB + PROCESS + 115 + 99 + 87 + + + R=83 G=71 B=65 + RGB + PROCESS + 83 + 71 + 65 + + + R=198 G=156 B=109 + RGB + PROCESS + 198 + 156 + 109 + + + R=166 G=124 B=82 + RGB + PROCESS + 166 + 124 + 82 + + + R=140 G=98 B=57 + RGB + PROCESS + 140 + 98 + 57 + + + R=117 G=76 B=36 + RGB + PROCESS + 117 + 76 + 36 + + + R=96 G=56 B=19 + RGB + PROCESS + 96 + 56 + 19 + + + R=66 G=33 B=11 + RGB + PROCESS + 66 + 33 + 11 + + + + + + Cold + 1 + + + + C=56 M=0 Y=20 K=0 + RGB + PROCESS + 101 + 200 + 208 + + + C=51 M=43 Y=0 K=0 + RGB + PROCESS + 131 + 139 + 197 + + + C=26 M=41 Y=0 K=0 + RGB + PROCESS + 186 + 155 + 201 + + + + + + Grays + 1 + + + + R=0 G=0 B=0 + RGB + PROCESS + 0 + 0 + 0 + + + R=26 G=26 B=26 + RGB + PROCESS + 26 + 26 + 26 + + + R=51 G=51 B=51 + RGB + PROCESS + 51 + 51 + 51 + + + R=77 G=77 B=77 + RGB + PROCESS + 77 + 77 + 77 + + + R=102 G=102 B=102 + RGB + PROCESS + 102 + 102 + 102 + + + R=128 G=128 B=128 + RGB + PROCESS + 128 + 128 + 128 + + + R=153 G=153 B=153 + RGB + PROCESS + 153 + 153 + 153 + + + R=179 G=179 B=179 + RGB + PROCESS + 179 + 179 + 179 + + + R=204 G=204 B=204 + RGB + PROCESS + 204 + 204 + 204 + + + R=230 G=230 B=230 + RGB + PROCESS + 230 + 230 + 230 + + + R=242 G=242 B=242 + RGB + PROCESS + 242 + 242 + 242 + + + + + + + Adobe PDF library 15.00 + + + + + + + + + + + + + + + + + + + + + + + + + endstream endobj 3 0 obj <> endobj 5 0 obj <>/Resources<>/ExtGState<>/Font<>/ProcSet[/PDF/Text]/Properties<>>>/Thumb 28 0 R/TrimBox[0.0 0.0 432.0 216.0]/Type/Page>> endobj 25 0 obj <>stream +HWKo0 WǵYQ`@71:쐭{`-v&`SGESZ.נެZ7 7wG_ۃטgP:FVĭ"/_jBA!.1<#d$x`j}gL6ՎLi +ہĺp;ĩF0X<(T{ gj',zVZ] +Ca{?Xjq?.#:y:ϭyC "Iu";hmo"lI`c*"7f)z.>QRRuG3urmp#tڀsD;#hBSR{SI*ʾ)dh2ik Faᓕ8!6Jn?7ʃL:KRkUL\!w(b}|,$0&Z0HF5b#I1 8Kzj J%jMffl*Ch vW{>u~?^󹌿t%W9tx,inڲ6=K,[l~4kMK7!|=_ g1v⸠х +F p kQF^2o[8o@>stream +8;Z\ob726@#X^bOpL-\FBqU%>K6lg^MBro+OFa>O*PG)opC0Bq"d2^G8Uo&q'VK=_ +L^36cdBd!21_;+uTU2Apgt/&fpoK7mAJhk?s3"LWZU/5G%3eh!OhO+='QIO endstream endobj 29 0 obj [/Indexed/DeviceRGB 255 30 0 R] endobj 30 0 obj <>stream +8;X]O>EqN@%''O_@%e@?J;%+8(9e>X=MR6S?i^YgA3=].HDXF.R$lIL@"pJ+EP(%0 +b]6ajmNZn*!='OQZeQ^Y*,=]?C.B+\Ulg9dhD*"iC[;*=3`oP1[!S^)?1)IZ4dup` +E1r!/,*0[*9.aFIR2&b-C#soRZ7Dl%MLY\.?d>Mn +6%Q2oYfNRF$$+ON<+]RUJmC0InDZ4OTs0S!saG>GGKUlQ*Q?45:CI&4J'_2j$XKrcYp0n+Xl_nU*O( +l[$6Nn+Z_Nq0]s7hs]`XX1nZ8&94a\~> endstream endobj 23 0 obj <> endobj 31 0 obj [/View/Design] endobj 32 0 obj <>>> endobj 21 0 obj <> endobj 22 0 obj <> endobj 34 0 obj <> endobj 35 0 obj <>stream +H| Xd) .Y.Ⓕ4}AE\ĭ\Ps ܻԫ.7-擅=뽏<3ͼsw~3i)mx8U-щP]*#t*)&Kh*M},=SܵRYcOHLO!89ql]9mF:gII>Y;OU1Dssͅ&=Q QU^BuԀj[FE=G4 fYVOcMWd-hh耎ī:ʾ@zz#}  śC8" !b`C,! H@"0ɘ" LTiLŷ1s00 KX +`%Va5:>6"%؄؂؆!J>vNn^'؏8Cy8ϐ#8&s'pp_gp_y\E\e\U|︆op7psvZ@ 0xUI3SV/\SU2+_NuT}TeSTZrWR'-UGʠTԖF:CָWg/ܒw șx_µhm-Vjmmnh=H{aZga*11dzlgbGorE^O7VOoAN~–*gifi%Zj6[#6s1equtr++ ]+9*:_WՔ*AuV!*Tũ$j:SʮTꐉD|*سBF܂[*4+,-(-^KѲlmmvkeV~ʼM0nJMwMK?to@]t v&[';|DNv.QpIa9E7dq&'{W +$M' +ěũBc'3b(4.܇7wn4!>!܏pQMCF=DAJx? 1- ,{%yጐc}D6Agoݧ{/O}J-yyjxXȃj=Ujmc;$I|É&j-&UQըcm&4!6uXhT46lAFcu@W$ni ck5}.˵D==߽|]B+%^RTbK䣅xx| Ըqx!.fټV*˻ +~1 ;&ig [ [+E%*:_:F +'V#V()sS6>[೹i+ka& +ӿZqb|\G=*. ({U֨Va[U'W/cpvsv[D!\;sBuc%^g#~p!\7pÍVs@8TFxH-lVJk o9ͳ!f v  o11" !Ô-{W0@ Z +bu"~Z5 Z7I:'ܱ;:Ⱦ|/* AQ>,j^V~|vIu85i8h8i4:\O H#>sO0 }b /\z^j}M<F07Y.soz#+Άz3Mgcs2Dc+>NW/mpZ-6[U8Zn"͸ +Q(huL"! +l&R,MLf!2YoU;Qj zjsfʩB~}jX:LƻGsf;ir(VfԺ~~IY}(H+_GM{1-cA Xa }79' Wp82x0geZ2 +l55z}![VgS+O:T{,Fl-zO+3 t|Vwbtth!X޹G3-W;>l`Z!)% a'fL쬨AjfFN@|YqAzlK8STOP_m"{)‡ 54%vfs7֏P_}SzxY(l|-_N^AXGY=cZrW(ӈNQ; ~Fmdn>v 4E=l=\4ğCAAays46 ~w#ۘ%1T'lJ6|Zq)6(busU_ 2U3Tԡ7AZ~ٔdCoXC @)Odm(MkD?2rr4d 8Ŕ1~} +}CvL)ȉԠg*ߏIr}GԾ7lЦk N?C84Sz]&D=B5cv&k+tio'䤱蟢Uy#-47_ig#XV+;o]\nCc +<\9\ЭЯOO>c5FM~1ODO/k@+t]~.Zm)o2үJ5Y?f挜u79Y mHE̦J\%WN0\g R-eOϳ=Yo釡{}p=&Vt=9z2{J<*W ^#dTBﲿS zFேRt Ŀ^vz߽|~9 zxz~ZCƬF^㾼3߫Y>oϚ]_W8u ~Jb j~Qo9ނĕ*ӈwWA +*G?N$FaJK~B݋~CJ(AGnK}8@x5~:2gK<>a3) c2Xoېׁ3AJr lB=uo,`! D]<4|wa;}`g~gq`s75X4.Mxqf'{q{s >87b?ohUeǿ fws9-LsP٤:umk7wl +_"zQ(7z + AA 2 ,հ!Y>9޶Eo|<<<9*/>SʝYig;<Of眹5a63RE]ͺMd΂~V(=жQ@;!E\2I EViN$=C=-G~+82؝<5w8ޝ_ګvvKQO݊k@G41Wbڧjfty7W9ݗT26߸6FIlV$jMh QhLV~Tٱ|"Ŋ +"'="3L\gvZ[Y[qt5Jx#2oLw'ک z endstream endobj 33 0 obj <> endobj 36 0 obj <>stream +H|P]HSa>zYEb}7\ F~Ƙ˹9΍?4;h)P4q^Tw]IAHE̫HBQUBw淠sMW}/2(8{.D܍{uq")C>5B2s {ҾGf=X"mPnS(yfZSXl !^`xTke9?˹Ɲ^ P!pȥ =cH\nu"VC +y|HBg|9#]-lO< _P<3 6Jy:JOQ<Q T urAy`|!GdH* ̈UIT+¡ԨVDe:_(ӣсŘF'ߞȱb쵄 +ȡ:L..n?o׆&iCi1sު[?QFM v&Jh8m( V]ˡ5Ya! E: 싅V{VSĕKJ` 08al9x4$JNNE L]֔'REDzYC؎{sg.OfjM>IJ8_) 6齤L\Y76nh`dB7 endstream endobj 27 0 obj <> endobj 26 0 obj [/ICCBased 37 0 R] endobj 37 0 obj <>stream +HyTSwoɞc [5laQIBHADED2mtFOE.c}08׎8GNg9w߽'0 ֠Jb  + 2y.-;!KZ ^i"L0- @8(r;q7Ly&Qq4j|9 +V)gB0iW8#8wթ8_٥ʨQQj@&A)/g>'Kt;\ ӥ$պFZUn(4T%)뫔0C&Zi8bxEB;Pӓ̹A om?W= +x-[0}y)7ta>jT7@tܛ`q2ʀ&6ZLĄ?_yxg)˔zçLU*uSkSeO4?׸c. R ߁-25 S>ӣVd`rn~Y&+`;A4 A9=-tl`;~p Gp| [`L`< "A YA+Cb(R,*T2B- +ꇆnQt}MA0alSx k&^>0|>_',G!"F$H:R!zFQd?r 9\A&G rQ hE]a4zBgE#H *B=0HIpp0MxJ$D1D, VĭKĻYdE"EI2EBGt4MzNr!YK ?%_&#(0J:EAiQ(()ӔWT6U@P+!~mD eԴ!hӦh/']B/ҏӿ?a0nhF!X8܌kc&5S6lIa2cKMA!E#ƒdV(kel }}Cq9 +N')].uJr + wG xR^[oƜchg`>b$*~ :Eb~,m,-ݖ,Y¬*6X[ݱF=3뭷Y~dó ti zf6~`{v.Ng#{}}jc1X6fm;'_9 r:8q:˜O:ϸ8uJqnv=MmR 4 +n3ܣkGݯz=[==<=GTB(/S,]6*-W:#7*e^YDY}UjAyT`#D="b{ų+ʯ:!kJ4Gmt}uC%K7YVfFY .=b?SƕƩȺy چ k5%4m7lqlioZlG+Zz͹mzy]?uuw|"űNwW&e֥ﺱ*|j5kyݭǯg^ykEklD_p߶7Dmo꿻1ml{Mś nLl<9O[$h՛BdҞ@iءG&vVǥ8nRĩ7u\ЭD-u`ֲK³8%yhYѹJº;.! +zpg_XQKFAǿ=ȼ:ɹ8ʷ6˶5̵5͵6ζ7ϸ9к<Ѿ?DINU\dlvۀ܊ݖޢ)߯6DScs 2F[p(@Xr4Pm8Ww)Km endstream endobj 7 0 obj <> endobj 16 0 obj <> endobj 17 0 obj <>stream +%!PS-Adobe-3.0 %%Creator: Adobe Illustrator(R) 24.0 %%AI8_CreatorVersion: 24.1.2 %%For: (Andrew Spakowitz) () %%Title: (weight.ai) %%CreationDate: 8/2/20 7:05 PM %%Canvassize: 16383 %%BoundingBox: 32 -145 395 -72 %%HiResBoundingBox: 32.588231612538 -144.40625 394.588231612539 -72.14794921875 %%DocumentProcessColors: Cyan Magenta Yellow Black %AI5_FileFormat 14.0 %AI12_BuildNumber: 408 %AI3_ColorUsage: Color %AI7_ImageSettings: 0 %%RGBProcessColor: 0 0 0 ([Registration]) %AI3_Cropmarks: 0 -216 432 0 %AI3_TemplateBox: 216.5 -108.5 216.5 -108.5 %AI3_TileBox: -162 -396 572 180 %AI3_DocumentPreview: None %AI5_ArtSize: 14400 14400 %AI5_RulerUnits: 0 %AI9_ColorModel: 1 %AI5_ArtFlags: 0 0 0 1 0 0 1 0 0 %AI5_TargetResolution: 800 %AI5_NumLayers: 1 %AI9_OpenToView: -82.2887423922857 86.0357908581646 2.75291966104253 1780 1011 18 0 0 168 140 0 0 0 1 1 0 1 1 0 1 %AI5_OpenViewLayers: 7 %%PageOrigin:-184 -408 %AI7_GridSettings: 72 8 72 8 1 0 0.800000011920929 0.800000011920929 0.800000011920929 0.899999976158142 0.899999976158142 0.899999976158142 %AI9_Flatten: 1 %AI12_CMSettings: 00.MS %%EndComments endstream endobj 18 0 obj <>stream +%AI24_ZStandard_Data(/X\h^ *r8reoi#H9vT R 3n.JB(qq2"eDfݩѸ-Nt1xO&zqsnqXYM-vZ;vUU=QumsFf3bNkfFX %}X-StRoԤdCcZF$OƜ rg*.cj%\O[xijpn|10:I5 vAEt܄~ uQQ0I_Q4BqV"5HdEqs}Z[b~bȋy\_f")~uJ(B ;$[2Ng ?*w;llfAÚ&>]yZBojA!_F:иZƆ޴32/lEQ?sI(=vv"~<5ń^5tUG#\UkdSG/Ĭf۱C/fCoB ƿ]:F_"vsuWFg!rBqgr)k Du2D%-ceDdr3"=4V7,?Vg c 7踰=fTj.%~!r6>Ue\lUР>eT\UJҽ;zM(kδ2*:vArXQ;wk9C2+\%{ OH7" k7 ^;JK T.-F efjTqQZqQ+>z4\E_LU{Fօvt'mo2Xk&9rZcUn +wT$RVϨD57ŹGݣ_Pu7pHF#.hw-HFцx+̘[ +~)kYdRՆ Sq[-QޢFvªFݢWxqlAI u-^.jvg:``158UDR4\KbT$f{Ƞ u|BIN<\unL%bA"(B'8]6k2d1rqVr^$"$ $$$$.q&Q$aHHDDD@D b Q5@0HQPx`K4H8H<Ї΃@PepCCÃÃCCăăă@(@ ! + hhppxxDL#AARB$!B BArrbbER@0@ uK\=5ϮnǓ0shbQP j$tdDd" I8vpﳹliZb"aanleHCw: GFe(D!Ї:ԡ὞ZR]N۲R)C,Ta Hp` pp`!…ТbXȂ 3CK;! $[O_UUC(m#c,.d  8 &(,@0  P\0 +, H0 &` +*E .h  )ń \` .h`"ă +ɰ 4qaB$T4 Lp\ a \ (,LPPA +P`BB$T4("PX ,v tX , HҠA,P +Hp0\` 0H@aAPHP.LD +$\A&@h@Q + +h a4hPQE  \$Ѐ4`h@"""0@0Q„X@qEHxaB D hP (\ BE * D (  (&4@CCThP"PXE ,B*(L $DxPQLDxpL0 *P\ *"\xPa&4P! +@0 $P + *\pC*LD  @@ (LX@aBPA\0A((4`$LX@!`B&.h`"B h aBT " + X@ +( 4@QA.4`pB &LЀABEEEE N@\x@Dh "TX@ + LD4ŃHРBE + +\( +*PQ„ $ha,PT &P +(0``ALx ! +@`@A)‚5Xғr5 PQA PT01 @q( 9bV/\<:Ȩ[]g L $(P`\YaB + +X„8AD0o@ + @D +X@!aB0EE ((LศѢ>XhqB,hKaa"PX,:^k,Pa&,`x@„ +((Ё  8@ +$TX@Q +&,hk@(*pӀ '^,*.88( (LP!E&P + +( 0Pp + + (LL"Laq +$L Lb (*B +@aB *P<@a$PXp Ȏo`qXPYѨ@A +0`P3> 2 ,p H  榸8E +**T(  0l!m8  *P11, + +`q*,@ +,T`&. Hh&H0`. aBE +pp h`A@a +*\(L0AA1&@h0@&'$6,oax&L2S̈3. "WUFEF9YLfvWrDE훹W$Uh E{#G7 G,#GD}>kA-bТby|zƤ#۾V0HYR7ESR~^LQTFG2%5S3ʼ!{x]KCc #ߨYTIiјE&Tc_3<S5zU$>8f\E[*W.Z<~\^碓 x<|&k{:סỏvns G03trټn|QECGk+]z;LhZCzVYQg.\V,N|eTU$.,ӳtAZ ΂|;,fdįΫ--R׭?wіV5ςl<b~qgn&~>&M w3$MYTH_#0̏[y,:nATDŽ3~\q6¼"7EF3601_KEdQ-Ë[UYd\b(*/Z04$Sjfad*F/I)Fbs>*MDRutECsEd*y-*3Q/J*-5hrRG!|%k,-L7.qq^]sqf3Z6θ(ӇxC *Cijƽ MgO:E-}{aOC|FLƜkC5#9hiA'!iΡ߼hQIn,Y b63fCSTZꦨFJ&Dzhar˞QBQ·q=Zr5rs"7˵C!*ʪVɫεUT,QQȧ$2΋aKdb0``j,KsݕfFd"CfL1Jz"GcjkT왓yWSxLgRf?ӘZ&{%]b\Rf~dЕs!5*!/NdZTl|f֢P>0d~]\tD5,YLӜR1 'g,"Ylh~47jAf85.HRis76b!Z}aRgËJb<-q협F E7KkUnsEEjɭ2lx6Dv߫?6DCx%i#iN*816oEdclaZ^E?_nk1SkTťA_5guFaRʓkF0|esU~aƳUCX;H?..7&Bd̍K_nB9)fU}A82~JSUb|ld^MO.GU'vث/Rսd2"yHv ulhֈ|D"7u si,4[IنFJiǔΪDWbmE&vy nj3Z|KQtQ}8c}bċ4cqJ+Qd)fd8kpH] >"+"Ϊn3fPfhGCrΆ I}g]eZ2dϜq^.jU+|ti|ER.JTjկ.ȮߚOŷ7M P VWkf+HIFϰi*Saخ\br[j:M>D,NCYe\ DVRI ͭXͯA GbyƎʡndyK55KV*!f\ƒͧYaQGGC1\mLoX\txT+G48H]Mh j66;Uoʃ6.\$qAt-#aQZa:f֗ ^0/2Mʏzq2!^.vF6wA/E3sva)jd(8 +]n iQ󂬗ꕎTü)=Bs25uúQ]pD*ew1V+[,512VlFɤ<"st>8.l9èxCÈ*W # Ɉh>\U;בp >Ƭ-+:o[hcb*BV}LPZhXd-X/'*Y7QS6_Fw^Qy ])! +iXX&ush^X&zz mjε%Z+k9c +[gjn#-!"-!uWP̮Cctņu-uWSo"*(d-:*5rH!KJ_dQL bœlڮl:)'I2\js܌9ҺA6SDŽ6DxMDN~T֨g4yK#&FY:yՒ~ +f|ZQ9͜:T73T5! Yʆ* -s sg74.وhql lF=off:uT=5l8ѡR'uH)OM>f͸TlbHPd̍!%%cdH+C 02DCCo972>nReHe>xHǙ!Bbt3CFq ar;34hfFM̍hyx!#Ӻəq46h[|aa2&0 YjQM;C,$Q}Q*QKtqm2̣TDDd%qu.t::6vi%!6 _ni'ac GDGG%pk&k'uY6B)2V_D*jIm6~cu6ŝ<mAeۙ0*Df cEzG("qnxG=f<8?f3^"uYnHuqQF6SQ ,)48@aa + +8Ѐ`\0 8Ѐ}16+))e$g nj'v|bGcU"xTolNUZ{B{yk͍\)2JGh6}QVTUbBX*ɬDduΊz3ݕ}XVUfzf.v,YkxbD;I*]!:;ݝwԷʮU=\ ˊ>͊feoo3t~.ҏnQbM)#⪑XL[7cΟ񊤳t+ۢu+4SSsI%"NjH[ef݂NCsEqnA/3:IrNttfh},YYc ko\ +IN!?*tOaER4*Cg[|kA TOm.EEUF;A-~ OE\j'ɍ І{S3N㺠j kXǹ0c_4öv ]kji);񺠪΢؏IFR!9RRXQ1YKw5̫]FfSlfq#Crڡ}~rSUZH#SS5TIhzTHN5TSTvYY]b!d5ZTTG4G*to9VтFV $$ 4Ӻow/YT8uS9!ۉ,|O,h<EL3 +iabl~ULgQu{PY(yJhvto1zb:!~Sr"~77LrcrSQѺ; &z #.1".11H ⢴nUwVS)[uqjf!Xg:dH m3Y=6u$u,':.$5ewfȤA jy Seq :xZgj$#-l}S w ˪84.|4f*s W,XFrOQG'n58#-*-]tը1,C;#juev&lv*N$"V+!|XФŲHsjnŢoHsSCXFNTnTUoWi'&sT8?̨/HEԘX f/Hf*4VhfR3e9gl5,:*NyPMyf>3V>qG*g3,d1Mǵ}7E̍=fr kԱ>7's3fi͍37vOfMCd(Soa#? )5Cv3Sj>vX4w tEez!:X3t#㙩1zk["S_I<macR}*Ղ2:a-'Njx.3ނ&5`[!_Ԅ2殮!M VzG#kɑ֢LBsi94Jg\n&.s CþOj_Sc8*|hԾ 싺gF_&ٝ9~){_TL*&j'+Kl峊jNѩLSBF63\67JJGU.0eD"65;Z"q+doʞ=攔oPvfO>Ko)٪fT9umAR9q)L8>K|.UjШ~;oΏ^B:})LYȬ&;ff*0#^929.6Je5^b2ޝĊ7椫HTƣ\_t,f0He.Yhv-RҹH睦l,I0V|z<ʈ2X4|;_j-\"eӍmjrg]q4zscZԏn߾b+Σ͉O OrБEƣ"_%] ņJj;8*gh+z22r4"kwd<ľUu{h&ȞzM˥:Uù6rp(Vb,9!C_)w>)FԨ)V\}PwsH6C7=D҈Bl_+/D6e9w U$hEl؈gTڛ߭>Vnn޳b2WNGQ8䖋,idǐb֥t֡!bDxpk5)y49)d 4]4gfצhCd3<݈şi!D0eTljDr;Ӊ{fsx2xgvqҥxnX%d gu:]ͅE1#ͺKFcGGxt:]È!gCPniXސW*/)՚9!kdql$ۑh +X9j3K&Zb#Cz%yjm/ |4#7F2ԑy7smyW5ȭR͕EfMb&tq>{XM^C6uֆH9G)HB9FCv׺|# B|^#+1߰xvZZI˟[X]qj[ٷeع_n2NvD4"b5oz&3u3S{7=Kqz:#Wsq)2VyZ~z25uZ|:ߠMsokkKמp!s 5fB.YK[Nm>}^e3nzStnZ]v:6n][ΤU )UmX[Uo.V83kjbWvUT2YތYT9Lf,éƇaf>6\\}jf6Hӎ~:YhG^REttw#Α ;Olw20`xpv7ӏG"*}Q*yqZP8D8N:dLc:ZNh\3qW%5g 6l6 +ʓ}r#rDzjIsW#Nb͑ՐJ˯ +;e<׫xab/StR*/!!19Θr)骗g=c$j0g6")8(BrB^AK2Pp<'Jreӥ!K:Hjq|VBv(Bʐ*pt(&g.V! +Rb`(aDQЋ!lV G@{_DN2TT&JD iC (\􋳆@"v(K2i>9؍!ʯq1DE3eF\}A JFa@ϣaLrÄ!` %`t`09ДrN+.Zaå^ZT`p2^8LZ,6܅8v!:tz + ~-BTJ! BTI> +tg b*bȩU- T!0&޸ThhaJbU,_Z1"XUhf, ++,=,w6XJi +NVjW\[!° +V: +'ɹҏžu<]U/ +h+OW29,,cq +%WʯTmbiXucA*7k ,+ʢJx&B2bSYƔ)v) +8gI @A -C}ɴ`( `·l-] -b'0r9!;iKw*޶/S v { ɫ &xp$3 V -=+ZDh.%xz.'.->]JXuy|ND:JH]>Ҡ&f!}# 9   //r@/@2VmGK +'(F|!!k}}iBa59؊Tq'«j!" %<Lb<*oE6g*MÕe{0r0h"Ɖnabu & 8<>[AbHIhU&9,:?@_}0\>SK[.T&̺2I292`w3x1.]d&h]33a3aI˚9oLG$Zs@7q*Qh'g ԜL 3(Y*kgtj0>s Ip4*- -j4Mh ѱgR9& %S:{g4l2<T!`줹1ri1 KCH81H08TOs xP;`ʨ!DJ͑ OMzʩic5hiWvklɩP@SF +g ~i ZPXM,Nk8[֌3,8_W ?RY+&h9jט KhU5v O} +n@*,u [#"55 +r5" +@+؜(ش'l`'8#es+!'!&hh>@isLRK׆Ll=,o/y%XYwHܰ]678@J㺩*pxua@ʽ빔Si1}{rD9UMx֧C /sS{%}AOx(W')[tR,ZVi)8d.lV8G`H3 rB9@ .9|NsF\e9vw9>sɈ Ӣ5loΫVsCZ90:WdԻ*:tZr/N:H΋ݽt#M)"xRY>v{P"(a-p|u>5`|ؾ=+>@%>Z!G=j,dh9AB, 8@D|x`rilblNkM9@K)`hPTl->'>"(o7BAki2(2& ƺ4" #Y\y`E[YW *Zy?/ ITvm߭V]]W0I]5d߱.S@LRUǬQ;Q06EuHGBʠjWy^=@x)!RQn76iߞځE Rа +O6ՅсU怿;p?q䀄YDvlGG'{UТpm0 +Q;TcAP@Ә'}4%1Na,ԑqqXHqhƔ8sP]qP6,P,P7|3+H4Cgx}}s^(yLxXY0e#mՀt MN= nw O*qM3`1Uy@ @TF/6 uQ"e\IHL@y˱y&:7jp{@6/mdj=0H?f "e̹ȡjl؀Y`΄/Ob7v6LG3 g . K>f s ƨ#^f]a' h:pW^3IӀ).]cEd ̓au: +Kp7F{P6fmO;;L_UVCy{ 0AӐD!1nE0p!6kCs|E0ޔ x hИdZP5QI{:TL`#'RF:!  5^yʃ70| y\`LAx9`B㦆 B` .nA0zea GHS돛B/ 2$ (о_DwK,t}l+L2аypg52@||y< +<#;?v  a`ᒍ60) W*DZMYP{ҝd)x #Z/٦̺pGfi V"]' ?\0_u\BWT0{wĻgv?%b-#,LmE/hM($mHC>'FХ-R 95. |2Cq)UOJ)*5Sq]W Eyp=>ؘZX9#Ẃ}~tdHf}4rh.nAþB^HC}PZr<}߼ +I>s#H>NLV@BU?N_V7~׶L_71su1H]h)0`ўxSu>k\171Ji?QL4_i[L=(rSǼ)E8?GBwл( 7 bFB Or(AA/nPj˧XKgʴSSӜߌ/u¿W/XdX'h +$Bw>co~tϰ{xads %e@w3~+meX"aw]֍[XQ:K/V  + +U"PSwhHQ7((%~J\?,SM0Mb1+6n3Q\W{/qא;I5ʂ?M[ ~k-^4?Fzpd@! ^.|3+ #y\H‡r:޴bL/:('oًD*Gf E.# qheJyøOH>P6"!EF\p\D# &&>嗿P^ z@ƟE^4!?GߢGhIdPCݼ7,Έ{kHSJdc}`&Ѷ@# +דocooqCPIWR",h# dBjRʫctPe6-30P +=24AKV4 3a D0g{/b(`3z7jc!Ǵ6u;΄:e"^`I\P)J}S q-+$fhgcތzƪ|״}E[%he-& MCɦ9Lqo᎕>5`6ֳX]z#r|YKU_g*#R `?6mf "sgXԧ'A$:=3)k-C'mM|j!.۠wT)85҇a+d=F߫ UFW[{/ڃ#ι8 C/+M227řځ +?LiYDry|O΃.ufu7t^ -jx,kz2~XndH#'EA?-u9Ɵ)x(2~XJ [µ'"\>?Ŧl.xZvZ +4rSS-e{&&`KbC!OU (7CamUx4|/ɋw@ `W«,34Jݻ3|}nF%nwNT)x&dE8_`صUDrׅ9yEaNoii/^g䉰u?n=^ <;\{F 4_,΀O$gQATzݍfZ#^pS 3]`93T=x'{Kt;pɂ|ػX2\fj2>!p8(I udz=׸svԺq͠!3|tSkh[-gz{bLdnXvZ՝ZoR$~ˢVu+<'sr569_7x]/ʺhzT6gc~|DmbH,eLUܮE"1 oՅʎ.rѧ}#l -p@l9*BRG珸RQ/lD]ﯣ6 _BF\TTs\ш;B +=n=[q@ǵ+[nh5>GDY10]Q{|\Pm\S: gn dtΉ^g$]ǎPD*"tVnjD@Ʈmyi7*Y^m%ry)42_!?F3ä"n/Z k ۹ +Lb0@%{ _(Z';՛:҇I#֟v+9p nG8w673cp[LW}hH؛ 8 I>+T,&=76 xzM+UF9 xO +hxݧ~L=1+!18nb?m=QU~03 p˘~/y)IG獊 @HhB>k˴栂EK:ߴ'?|W'G]+6W}#a ZfQ{gP9bVRJƒԊT'HDG&BB]θy/@O7qWpF)Ue|qv7YVxk&h{[T + p /CG0j@[OW?ns%On%Ew/ZjCT]?t(S_'tGnka|-O Y{AL9 s5[Y"u]p,9;eȆp1n6l5Uw7 w|}PzQ6nRk^b zk,0\3ZpS$xۖϭdُ.PE8p@ywAX.q8Il[մ%nh$dvE^{܁ YMpo@1+tL ʶ[Vf ׏./to;;)E( xM3LԄ |f+(Lb1Zq MimqPC*vٞ(~=0#)!t A#+< ?#Ku^{kN@6*BG@9K=PC´s~c(aB$8'*k0'm)M$hE-ɐv-yf=0Ֆl< *=Dq 2gmQKsMfybjvF!>p +.4, G5To.-n%X2Q=%>ٸJu !#[\/2dGւCv0*2c7mѱ S5++$fSd1b\23Wb4*"`0&;@73_ Ⱋ2D?QVa7L&xN˞ {^9m 4],b2O׵wMhѤ~DH_W9/Jk=k EH8:5 +"]+݉|4^@4^l4զv׸,<dJvzy50گ`/;;GF|$56T^ć+G^&7o=!;QhȏR0 +1BZ?KXH!?}mͣY]hֈ} XYpAڳ*FDZ)T>CDq!b+9WjÌ[dX0VO_j^6\UC߭A"zXQuoxi~(T_e*sQ$q?5" j_N=6hS4<dԶxbFvښ J*I-hQgH;#`r٢yonD} Ivj"]1ƈ& 2/Uwbd7maO AUi^v9vp4ptvNB֜VmYLdvoxSYAN4/t,uck|{n:o :e(2=ȉ1 ;~L[vu4B kt,-<ܕX?ҋG\تF+Ta V揆c'ޘN*'ʂI%-$ݻV{O>tfҺJPH{x裝Ƭ Ge+B :Mou^PrcwN܁Md + U&3+F&fM1061wW`\5"d*s Za189_B8̀3/;b31ƒ0Oa|ţ;a@ [لYr YI^|'I6UKaT0aVKB +S,xI5ꌂyI.8J >uo:0&:̽gL +AEfYMKW| vΎ FP"Np#Wh$Dy 3 +\ r6U=Œ]tV{0X Q`ff2:1>de/(F;3_VMc+D}yL2 eGJp0L-M)^n+b7#9|.lUl"޲;,s͊ejK5MomfH򱬙%cw I>6Fl]yP]syLZkAUuI 16 &+dkn `f7)k־.Opr=N2^'p.6-lj,R+:9%<2Ks{s*fs\!3DڗٔGkp:fN՗v+v^@v(hF!yEw2%4MRʸ5IN9]Fc]:}ΎgHkVӝ04F(ơvG;Vl|>lbqt;<^F1e 7(-1"BK?iBDzVB_(ߠ9T +ġ%i\ Xඬx:襔::Ed5,J[1!ѯr8Q6ءPܨ0U}zr:Ot?jHՏƱ2 Y~~Hpi(A&NQ;I688:(}M)&JRMPf5\GpK2񹦔v0](@NLoqA*|ں)Ҧ A8_{0NP u\!`m"O,e)8Y6 +XIUBPаtVQgPGD}(: +Qa`J;|RV*Lj>*URLpP<aA,#~Rȥsj[!2+^jn¥sT< Uhʨ`ê֕tZ`NAXUNX +!b1Uxՠ:)|Lm,<Bvd:-c+F4qVr̠uUIPҚP:z{Il cVkMpxNߏƸʮ\˄JHĜ* r=x9S^wiI^%^ђ4G}]·_9'`Wݤ 1ٵv!_a'dacvaIS?'vLbGRA+*^~͓~A]U]Ar.GMlU{}Nƾ __0֨N瘛%RI9 +{6rK>IHhG!(-k +I[Ŗ:հ9l?`Ԣ-[NjqYN&v`]2.^ǘkL: b;Q%|ŸM[a$mv).p[*v: w3+ )K8mr6 JumBo.Hl.Ldž:B?L>K⪪öbXAr|(7bF oRl/ +kT}~:dH8!]knL`Խ~q0u=ib}n&HKL&4.=p7]^3-TrĽ{%ގgD..xP;_ 7`a/8 t}PU0e:2'=S2e˾J$4[Tq멯ťqc6vA0[P}\pDB}s$c&Vsu5A Pu:yʻW䱒 [^?v +]wMo5rFbo{;=O_aBesP_#@TR ?C-*@b 2pg~%,WʥB7̲;wJ5ÿd{oƚk50|GJgp8qo,^o–zoaή+_MQ | s/v [ޡ|VϘ?i/442}-6*WR'SN +[a7sglڑbW2AR|5x~$*sGS|C!@VUsڐRDbdIR_UNJDorƗ %f!/'ZXߩ!W1ι~j}zs7Xn;$ėU,rN촉2dzs3,!0F5\m؆#E|7NB| _+ +ķ 7W+oX V[l+9pY⋷3CGϾ>*b% j0;fy:0|F;{W;N㪢F| ,;pҚ bS[N|Q"{H]p{SUEj'שj-ӭlJψ@QmG0oBUQ|VALm^ZݕCJݞɗ˩oA-IJXKR|92PqpB KX/&MkPS|7. cߛhWy?¬Lҽ$ !Mt8_1}l?[r_XABHWC7.&͜/a윯33">|ACݽWRB/mz2+<*qٟˍ󕹉 9ss/T[9@ +#mҧx{$`DRy#7 +y}#&"E_ꝷIWj)7jIz\5{7QFoYp!%+-1yOUs]o;صhگ59os]+[-IFdWiL/E,TeOe<O.:(lnELeFJp2js%k + a/T8 + qv%kb +c7\&]Pos* gLU`hŢ3pը\`N̬ª̈qzl9҈ ~#O?4=U Yl^5nI7;-gMͬL] ÄNH O6 !FQ9a01}>RGD +AP_11;9?dCUoB6s keWfgv[AcտO@%4Ε_(F٘nE$;q+~MEs)b}|XML|sky~%EGhTLBN?bJ3z zm6>!΀"/w!7@?dVRэGƴCdO` g`-ş#*J!`EٕEKUrٍqQ{/pZ G腁#sϤ(q>Bko a rn|&`L3ŀ 侧}EodO_ 'ЅfdD1Eؒ] +`Kr_#,@^{A%Of fI\=ȡ$Z(XJJx[Jm~0F}ESrDJy[JV=m#x&%[O ׈er(%'ȽrCBIe=ᇷh LɀϠ4KIHQJW\7G-XEJF\ՉUa %jlU_,&J"8kJdr[A:5C$>lGftKd2!R; !Hye6ܝ,z$ye6rp& )iʘMy?E*(hm4 Iܯ rj<-C&ü!b2Fඕ?2d}_l6<idΘEdܕ*`̗t{,Qo2f +G5PJd"hffSZrs-_fnF9uU~MgFJHdsy? ]3vP=|_=Io?Syje&+}+6h`Yh=z*G͇Y#:^0$+ )X,$u1lMg% 5Px$j`WAxC FP@TaBP99A=0P+MP(/t⒠^03Ps^;Cͣj3Բ H jy;{+PÚ X3k +­=A@= 1oQ8&uz=)- L}pꎥ`GW(Nͨ* +1R|Rk, Yʗz0#Vb](*0UWH꣙ 2꞊ɭU3[Zvnmzc(񥿚 X'/F/Y'h֍hF{aaqm"m`kk+HYk)dL ƴ5O+Bm)el+"g Tml2-U8XCIQTZO&.ZfN"uMb3.^ztru{J/1eAO5>bwQkz)BGWt!nխ <[;!4o]ԿOҏB +1z5Bo)2z?䈮 Ƥ7z /b:;X[c.ehf'۲V1Zz%aK<0\䘸\>|/_.嚼0.8uhU3k1FɯM2.{^Pع>/ +k및Ç ːH9Wئy 9a;LQC\l*=e=&f? kvȦeu%?Ř_1s@U)V="`\]m ml;z̟l~?=)7n=],ǽ{2WCe e =thQQk&u+̌v_8GջS.7Q`&{p e?wd/0)jo@>2|<ewpnajٷߏo[ww3pk +Ghd~BfopJ@ꗄ[(Sis4=W p8ϛywZZt Q|O +U+F V_,\[q/缎'Ǒue~\Dl|#ka'GV[MfYUa9ŚhV6A//f؅L~ oB3'1Bs>0Ӵr(y,q~4" :@ow^E-Tp^86j[cz2T14Cgz}R]stΐ,>*HEKu { A/ +`-E]ԯkUtԬ▌Mx:-yqJ\։{u: ʚ@OV="vNcvi@5xv*>Չ0]D 4%j&u{0?x!tǁqt`6T莎s; >6cN;ٝh»/{7a}Y:Fy+sT735j3|_S^<£uyqK HܽzP(V xxt)G*Y6 IE! +:%L:d֍^,},-ů")1KF +A'XS,"Yxu㐕NȻ!97 e>0b5Tk<;L jEig/$Lّ蒷Pne6jLX@`a~d`I%gix0y_ŒzT$@EOS4KTBϐVD?P}F`~f8]ȃr[R]顧1G e=w|[,3.g`h@ 1Zg:һĴkSbt6ܞ s;pR ]L~Mȉ#X%A~^{vHq>i!~8yӼ{J5Kϙ'AP+2k_/9|d' _K+} > +}q :,sҿ;JQ-`˾}rJߵ!?.dή(9 )*S[kc3"th9j +wcf#ϨLX*} M=mNhݻ>iV06 +j\r z7 ߳+t 57C\ )DjL'ob(:GSKJr>IZ:[;-;a;dHv < (=Pmwm'y0`[`ܲr,j jE$Ee%+?U8rVLbu.uf!7N~k(5y_'?tl D ʤw2틇 +/ԇaYcX؂|myXy\[aZ2˛Tj:`VT٘o-E_V +,$I|>iebr7\U"D$E7B')̌]2s<\P֕*l*ҳNNpN,/-@@Q:? +R/k/jnH𮬦0v(ݶ,н`]^5$ Â8A8A. b=Y;lC!P6Z?3?~Bt\68X!ӅG(x*KsΌ\殂ϊZe,(ͳA5M2xBo]Cq}^yAƃyhm<]v۔݌TOtMnɒJ{Ѝx̻!BT#sxӝ6w/Kj(`U؇:k+ZUX\Ŕ('LPϽ67]1"2T 5yP +OrPc1Xq}ԒF&fN +Uf;e1r(">gwX?{amԬyF" ^w3sEG'.uȒfWлHSVO9esQI Y ӽ~ax-I*b+uJסN!JX%m@;1~767"S= Tg+UIg7F3+zcAMN<|E%XsB"+ImRZMPV`VpSCIpo>eC ͦ%2) ްc#;Zqי2VF|c=N{F_2 "YSS`?g&<ȣjB&H#ORQhXq:G`640q)@ύzBDC5"u@ StZdpPes~6J1IZ``^I҃!pѨ(;R;; +HݻP\Ƅ)^&Ksx{X)t<)\ԊW):ñVPAge .{a쟟5,F.j|"1Ē^7 54TGzcAD4ZwS(|BЃ9N35 E.z̳RB|IĒ j_?lU% Fm`H(rWي %BBNk< (/mV{i\2p3yB veo,2gaxH\6K{C<5;zzqgn#=6ҡ*O@ouu)Ksω*WN|FD)*>bѳ@W;?qu>H[-{ ~%WU1|m?G9 ʪǻ*BuWL麐u_Q&Ia؅]趪l?Ow+ sR!#*gJӟϔYP{R) +L{; j^P&똖4O4LHJ_T⏅/LW] K裒o +#oM35p;Ɨ'2G?L` gvm:I%.M) ǐ +XZJyĒo4ǡ,@+` +}IY38hw]?2Km+E q +p +9f+F-I1֫dL$820~y $*,s" "F hƫCy`A8`.ՀF^HxMҎ.\RI>wL 'ǘn$D(252BL?,TߠSi=!E=,>Nӣ^l7G'@t745"+Gs|'OJʮfWV} 8P#m jSҋPcdBU4्6 |ِqp@3Tˇ`vqo dXDQ"4ﵬ[t-=#Sn g٥ +w7(H9lz7%0s,ZyMsmʹ.y*<>0@0#k]Mw9R5@84ND/ԢyƊh;g3R/}j"D̶ckrrުi,T=9>yɳ9l+Zx( i7:ȥ\[c8Վu^ỷΗUSZ $ă߷h!x1_apa¿y&*QHmӯרgE+n1sZIցb]S&b ^1:Hp duGx1DQRb `^; `"l]|ޒ@Gc?2hwd<J,\y萊ƾD9C'Agiƣm\7" Mbb3g܅yK߰ e49Nח\6 T^fv494%?>f2ܶ!_Wa<И61i0 6˷78rҘ'.@˒Sʗ[ј bg w9[!Cx6%fXu. <vk)[w{%п3%;gI|wW~JwVKTB؜î69Qi 2xWZEA@~Z2N.X.dr}癸=gi|yL.5rm mL`B!PܢRLJ3YB RbȤ,J0S5 A@LNWؙm +iFc$ "W*qP YI5)1`#YP2e$W *>؁,'X@ {Wc+9Ii7̗dj$ {CA9l:I1`$KTPR}1#IP<.2*\ + sD~ޤXwP0u_D\h9+WYoo̢nTi\/3R٘);Ǹ8 FRʱ湁iQ"/gj8*Z^ 0&{ݲrbJa+ =0RәC:4k6ĠNAn٧#1bM +代/E/@.to.HYVފÓ]Pg]HQ +[(`4V@4jA թȻm;<]`b-x*tLlAFFQȠk3J9Sշ@uL`s +P^R}>mmϑ MtVAZaxޖyZz+ۦز72~̏.6-Ha8aFQz-Ȉ {ۛpqF԰(ƴR{#ɶӻ4/cMfZ>!"{^1` eH7]PրKS<7 AWk/hn5GX[RUg ͲwDZx.=KO{'%3ʐ(3ZIMO;xS%J} Z!# -W +laGoI1'6" VbBڂ[eescAb:>f|`t|Ռ)ȣ9VhgE  @Bھ'UU>0XXfX&̀䝵 +3c/qv[~95I|tDR>TG .ZaC:z/`yn#\C@7N5N%3ei??'(/@P^SMU9c8{aȡ>G0/ɯn!zG0Ug( IPD̲99#"!gDKVUGa:BOL摫șzbC z c8s,_1Yr ,n2o*bwiP f2ׇ61p`5(F! D(^=՝α#uAdކ(Z@x8_8-ء#nhځ_/uo/I3)vAƹ4aVgOͅdn%[N\X%Id FYu&A.ƹ 0[uSf)Uea(D p6%*wD,ie *dsRi6%z +{}㱓W-\Luݎ=ZO ˠG,1`)i&TSb`27Փ4C0%u#+[۰}W)-peJ ֪%b!:JA%U܀puU%@Y,LTj67OGm nU&y+^/Kn# Ts/Wv4I7{Fb؍,xFV?Cji͕>f*+<*2Fp?آ:i9 '}ftgW^6>LӛYs I( R1b>uCZ[4[tVjx<řMYŴHdt,\Fl6X)˖NDHnb^mcb,6w(DDMQCW +n\䶗 裍NfwXP:! g!я1iBUn,{: }~~2_=zyC3e奚l"$">VhsWB} ,X żgO/[[?<͡7/e@Wwsw;u_y)`5[mH%h$֍u6@7w\3M@3fۨ)8 ,`N\ $Ea`(BQs{-wxuJHrN7SU + P7/z7xZ{=rt Ѷ]r:>[*0}d cP,NUǥqy`2œd`h-՘eDre2@9l^}L)T^Hk[;ZX.S+6( !ɍ#MS,7W,Hp)6V +1 WްyCk/s#%yȿ3-Y9t4@ ̴ F̎l^2M zQs͔aȴ3*Ya SxkT Ĝ}&O$ir-uh{);Qǥc3 C(zQ¼&ݑ-}'K&zjP՛DҮ+]'-z{5=y[c@uSFЕ .b6q q j)K71>e+dB?zW㆏Oe֏Sz&OvMix'ܾ.fPKD3`gF۵, +.'7 v0/ް`.z̈[Jm[2|#GҁV~dPb!OBo 镦2a5h0jQ0zX 4cJ3<;zW#hMҵcgy|l@wb s0qXC ڤfieEk͔6@~b\$Sw;ìRsh644^ ti#J-Jv/``gѶc2kU0pXfgQG X%kC]~1"BrfY*ˎGD=JIYˮ3y(G|^|TaVr@f4.%[dL}'x1 w,$4X+WJ9dJ8Kk+ñCi\*K9-(2ŀ0-Ua_HfP[oZˋ0mL̛~2XNUčre"P[3%l+ cQ$T$~7&OU.u3ALM hlhwkrsc"?thJ;})d.M\cyio\Ӭ(=+Y81ƳABA>0ofm+|Xx14 J!~B +4pBȾCM]N(>lh'Pl ]8H%~e!- eLWzWNa<'DsXMԤ8A?9Ty/J.dB Bh?. +֌_ܞ_ >/qXJfg{:n6sHC: .~O娘H~%Q=o(pֳ 1}I_%Qͥ65F.b޵xCbƿ#}oKk/Í-PzVunh*nזF8_wQycQ4.}Z~ɺv&oVRb**+xVOc_mŷI_ wv +˱OȓCSyγ0@m܏@ +L M 3;hZ45/7}K-{<d^H`^_vr EK0( @z/3*~1$P 0P:+Җa@UP +3 +49H~z s1˫E@zo8]"}ږI\ۀ4 +-K:PJqJLٷG" ݢ$ +4ueV(ehmN,sSnfr]_n)m0\xXE+:zzgZgɴC=Ts?߁x\} h*֏yD A`(O, +ˣ +K]AV%TJUJTT=8g}/4˗L b@T)"rt= p4Ɯ2)]-v +/܀{=n + + lSpi]a# #'--M"gqj,Y^);ҺumJ< +xfj'),Dj9ʱ_5:.ͿA 9t[}T5stIn=,--2,?Џ/' 0БEeN.E~dx7 2Y/4(.FA`5?^w׮G"M5E!)Z<>6"صiY UϤ7Ƣk+83яFI4, g MCU yLĊE!n܁Hj 眜-D tq|V*|je(sCCk%]_a!E~pCI1d ECQ8˟"(wȀiCP`΋̎(S̈́9+;-t@hA#G/2ڟDAQ.W$'n]=?g=--ЊH Ζ L0{: D3%ezJZC|è^fpD/2oF4"ꗾ5*.OOfi_sg4d4%( bֿeE2C=E=mRџZQ#~L3;KApg) ec)*lw4CgS֎qGsHΥڪ"a~Hkm}kwle|%ѧUdIUE.C1ݪ=5]wf]tۢU.LM~97R9|w^śNN{3R2M4#-dg2oLVkC*\Esl&l&qwg펇|;"ifiK[3C>}G:=T濦5-223Yլ3gPYji_rʰY3cAPX[KUO2Dۍmn/iVݛa⎯2u|6~wT s3 C5sApdV_X"__sSր`A4 #_jsjMNG@Ijh Ed + X' Xb@  eT Qqx0df#PVnNsѨnM=ڠE^KVu;u^6XV\g5V >kO}_TE9>J$)^2"W5W֟L^EjXN3<ўGTQO[z:](:SCkMi.j)swNwvƻk^<#$o:Ԙjƛ33YnihYGq׋73o ׄsbke <(ͩkZKh$;4̛IiIJ}ٖ,Kn-gvph2Y7YxVJ5TI*i))U8^u>4g㫛Z:<ޫ^W3׎ݯZMD'މzΪ7M+5>wT}f6u6g>۽u;`jM/}Kt7RއwxFx$̻rGv7tRV2$&}wb-FfOV^P<(LJvHMd97/uz?i^Y6XPٰʆXz#ˋE2ozʋc398I2T+3MIIydKKK,[^:f}͝3\b朱Ls&Cf^}jմX779NM$Φ)IJ:YJҵID_h5*jWIk]ړn5w8uȺ3Ty~ ^Df"˰EKcs˷~VMo+G լ:giKK#դG1)4j7tiӛZM5RG]ygi9guejΰ3.ͼs5qp84#-<պgw<r·Wttvtsnq/|/ٽ6TzOm[>oY_'~[ys{w^;O|jehmQ}ý/gVd?3Qt'!l &Lp͌eg>mYUUiٚ3CB  &" ,$ Da(T!APX& $B&Aaj1P L-fJh1@ D!р8H G 0ps o@, + `@$r8r8u8CCҘL8 DBiT$ʄFHaRT CdByPi8Ypph,=(FcTDCd,$Bр$8Y@0I<80$ dA84 CԘ`Y(@&0Y`@Y ɣȬ"",ʃҨ< Ap@ cM@$c f +jBE"( yT < hp%B.UF TS/k0WDAWD_*Ss*[]S7Uw u;;EEUAl+¿z}\߻ȸBng3d]Gg4lx-!XZ!'w 9r5eb#wSYE2(f$ ȍz9ێ1_^62Y~tl?&v~dx9˸6ă;LDݐlQ5X[DJPG W:?SEVc(2CРNaY +}î2n-7AМ*%() , }*4C +xM21MbTK 8O^zCPv?i ۔(O0եj Aw%13{9#%t6h"V籴 X9hKF{jLu27bjDF>1#<1@3nQku՝`c?e`^Y <8d'8 틐U싏"D:֪]Rk[uϱJki} }v2%KiJJSHe:Nb#@p?Ugy\Ga;DC3Hpo"fT3Gb2vB>݆SŕRrv 4EYz _iVuYhնQΟb%B`d,.*7gE'5ĽOY҂Ȑ9&d٥c{}l:XЭ`D{#r^ nQٯKopۋ}C<Vos^= #ھu)]v4id_c$ 6@%]QEC`#0}jpWf3bF~ n{_5ȂJs) }[1.zgUWܭKbp?,#vFnkW+ 1Ը-Pq/׷]6g$%!-R5KGk~ndM`=U/.}f±`ȸL0(n]2^?V>myPtܢn>Pr2yt^X8ׄBe3겿MuAP+*\C,j72h(#u42Al4[㌾ٗ0#0X~%[RSg/Yť; @d?lG'^2c|-bX{rQ4DМqj  GS SAIP. ZK"o\jW)<$c`4}{ ;DG]:FddA1{*!aDd=WY~%fԍLd\_l#bRHl2pbH#hAAg+'8ՉZvPn=IǜܹeP_pZ"{3YiͅuޑW4boLNg=ALzer=p-5|+ fEBt'jGV[q)_Ixp9]be%n(`<[9om3,lɎLOM|%,hO SbI%JQs?ԌS;(0@qC{M. R8N3Wz:nS @2(-Cx#Ռ?^wEZXrj&<[Ȳ32ADP-ޥ9]-d} !w!vwRpa;_`{fYl"S^9p&#}s<N@]i=9%h;w + = Qq*2Mf+eP7]6*ddMiFaixc|ȩk+c91V#8t1aV3Wϊ4@0(g{.A p 'c)TPsI%>ycp >6>0hԙYpTaDkPenՑ>]Tq"OI1/yxDR̗fenrhrdmCn01DMѺAySW#Ҝ&j&κQV;t@. '9\6R,Q&L?J ujMV|2q14f&ӛD +#ɆKky=QxS1Xێ܇{%{&dBJhr3\Cqe3J%X2FZ t0ru0@}0vfd9>UJ:*Mxo)V*-=ĜFzEraJYRēY rc TLGϑj｝_L p~T+|&l cd\>3@`KmBD<{_2jƷ+Iig7q2ǥ1Ùf` {zu%Ưö b,<qR2B[Z3݀՚st37rȸ;,?ߕ}iqR\,xB,7RMLFE.ՂzNrsuGUEedvpl0F=mk2Z:#B ] Z8=@xGѢX=gvˠehoMGٿ<+@sB˂>62rqO >JFpM_i0"(19gn'zN 3' JB5l Bl=Zbu<$t5-L ؐyEonM7A>4Y a"r;FXU]ߛF"ָsn/3~|FhAʤ!(@"ZK 6-d uFmz9x隔\1>odjH169%0'ԶDdLX # 駹7~7Z[$ VB;˞zN%HJAL-#GnN~)r4.xL.1z1!.GHjRPR=T2Ǭ Eꊴk6ĕj΃ZɅԠ5ODRɖN6EȾ C$[jCT.eMcp^gL}~nO>J00)- +QN*ɚ.~W%|%9cm*[M֐odྲQT)f[0_9 +@mAOR,d-U?~+OʛիT9D@}!+>=-4|:.}$Ea< qodtZSeh/pO˸˓X6.5I=Tk|D`ȉ+0th4%Ӹڇët'2lָ >7rO@Q,'Ru$|L8ȟ1!{`?wQ]`3 c=DCG+2̘7b0%iV}ͯF+LM[7LI_  Թ H(b0ka(gCa?oX]isսRP+ -pX&RWCLT_I _=G 2_h$6ɍe'~kl~i !l,[$\몙cF(QY9L/O3a,^)izlk TO@&uwpKuh+і% |7.N Rc~n!=X><Ў☃NЫ`>pN|wBdw>ڔcL>Xp8F6.,ގC^z_>_AGEn4;zdpޮ IxyW4@Um}p"gVnl49- >؋l0nw`|L~2[tJ'W!c3"-AWt>tSwX:0L^8'&8`Š}9Dc%}c$ѦY\3Q„_b'PYkTgH@ 6;Fϻ*Xv05p"=ÆݟQpZ0HhbUڶ;Іڄsì`mq4}1Za"Ds< 3DD26G[5$V)N{nHhfCPa䨧0Kz8ɭ_bU ǧtR=f IR*l?GɠB+B쪌V-0Ѱ]zULx,fN/mW*hdDY-갰_9ыh䪮M,Oϟ`c(+@*Y tFV`AdՉk4r Sf$cR\PR,lL-l|圐 &\4"Õ{a]M}x窎mE|iE˿E1qKʄ9g~NJ8\:@A3W2Fne:avGݶӗ;z$xt'|mabc]zvriZi_j= cQm@ 4Vd}eR?l$\XkOR.V]ĆltȺIL9-xP4nb Я*incy +LǻɫnwԸ|̒XâU1+܆u;r_h$fv3j.s~!}A{3y&5%h~ܸ҇}􆳖W?yZ6+p> u8kH@5ف^!+Ϸ O6ˣx%LPg7}gc1Zz4*NEԭ|1&*o:70s3~ʩx44ӛ9oi3EPS0C5ĥh^M:b܅BMS`uCzlr@j#2{x4f8vmDœф/h^Q\Ҟ@-݌ׁR+ؕM$P^ fl\$$C"#ܛ̥gU,0껕"͕*#D F,_m1QoEf< o!ũ**HO qǶt=; ܚS;~lAQ9으!X?+F!,4$B1h$Y7(1^YpgD25 +(d7W/q /11SFy !Yʺ 7%f6n5$YdM+ 뚔_$R EV:k,NiE$,JufmdkqV/呥 ByMSksW5_L^M-Y.h>m@G/uwUBMUC_if}E B*$fEZCaULR`kg ! cb '/3 /MliKVR +JiNmPG|$1ֿ(tNnXG2t %Uuk|i&0vIpezz72{9۳NIXնVUryV%{ )\Qy^P(DYZr +tԥ&XGsH|IBx`1 +*b#֚gipm;^G CD&#)-<J& ki5OymABE*m^qڧ+4sR${Aϥɀ6C tO#kIYηAylk0XxFDы$&=i90K>z$r1Oz LgmIYZ +\;_;&]c 7p2zUI+s'׹0clt6uc`M[c= Lld8J9 &e!I4r x.H>*|wYq ƎPӪf k#0QM3v}XWznY\NcXٯ?{+Jn1hUGI5jsBEn +̺r66M4to.4?e9'$v&PBZ* tɆ s *D<Ԃ.T=<v=4ۈƣV[}iErC_-G&ܿ !ҭ"g.6kԲ@'9rq5. s{Ml}G +z5y+1ŧ8OP#ċ2{}լ$\w؊|+HẚWn +j9̅eb{ ["N̈́Rc/Ny~U,A%;SR +Uk[ŀ:;Įŧѐxq`D2JK^[g]SDuv}ljb"XsOQ<m`[&I+?܅z]؆QTH<ä4hs9֧v(a/gx IPʇtY~ iy55d'c#i`錎"UK7O +d;UܾEO傛6!~>g#x`,GEJFr}猉1&"zg8AcAG z-9a,16EQRX 2e #rj +:EI\58 nY##q)?& 1woKLrr0Mb_"tNY7Xw^Mb ͬIŃ^aF!/JLI~rK Ed +]#ŕw i$6Կfz5-n W%EjTd0X3h- v2~[;Sej YKzՅ5);&"4>nS9&83xJ( +,2Hi?/]V>kZ E +*v@PR|~.J<؀+,֎_3Zhv'!2PSy0U"昃9vF M ~EUs:`C<%3#c.hS]N\{REyC6HxOWF@P XEյ#U*fW؏{2Lu2GP5tfn(JKG +c5`^5 i3Jƾ@\Ndžnn]B3qgehK,N*1~8 +Aw &WK\ʪa|z@ysLKm <-YZ/k\՚v|Ʋ̦qnmܱD2mM &i f|6+AAlNcm( ܎Z҈=of%1g"+:-Ҍm E~=i5>i(Y…(2b#ܖRFՀ|mͩ2@=!e0Cr%p [oc`i$'dž&PTXa1 d%}I@ B–pid2gXUa휑i0a8Nfdwo/()KH|xj1{)c]T{mj`8Ʀ6#E\/< Qln,u$E[_.#\UF2P8$HWŬeN6tZ=g-YE=XB6mz;H ~ALD4*m P}g $FhنN buܝ ˞xm4^ 8rJK91 + BZKR_C7CC{G}Ii8*t*[nDo`V(}3QQDJ2`>f>B -&&&vL̎2pRsG/{78Aan$:6GAvHm$Me h+ŗß2I5ܒxZУ[/KlI4| ,@+lUp5HOGBRN*GKu֬;kEK{6qJ(ŋ P*53AFH! +qncG¹@thS-x&M}WC'dǛpB+*u4j^"~Vb,z -Z=؈ \嫟Z pq`CEsJ68g( >.Dm+m 5_fHAHEA$tc'zZZ=W ͅKiH8"A@ܕM %M0iȘqr0J*%̻(y/>#ruz?lR)FƖ8NޑPR/ |;(s{`Aw__ +D6h`\k}58c[ˀ}o Y\xehKV5an+Bm)VP{낾ռ qsXr~-ʦ10JTVh֘g:$I*P3/??2_z4#$sbFce2azA%X_(W$ڭGW6΅<wnj˪M#xռu0jsFƌ? ×-%&XW6B)5Fj(dt +ٔcd.A8PVd@!ޣ?3 d?X}rPu diG:?kz݈+Ηj&|XEz0B 6ûebeXS˂sY nb|%/O@/{ 7?&:ns֍nb$i䋉T/^FH}ܙq's/E?Xޡ^쥘̬xrE<ɡ?^&\j D"R=}q.*BO!sZ.@b@ a!҉X8DELkY`Dc"D-9n!A h]H@|:z e6HD0h_0o"TB b>"k8qT !1\<0G"02`̘Yైئ3}bh菈L?La5ԅ!U $֐t bf@D` BopQG7"CBox/2p`!DTxZ= ^0&7p!( IDuwC]46dhh;0?C>1di5 q> .+!B`C,﹇Dɂ(ȣRӇ'dkB)?`ܸ\Bq5 3Pa(d D)sAeE ] Oì 0uAX@hTT5 +A58 ȇWj0DHGF$YbGA@F36fT#@T r@@$d|YH a!azLG[,L9FGŒ$3Y?l$ ?$P_'P_c)eXp} w%V8KpW@zX-s&Om5ZjNq>8a5Ze =儆l:aȚq Yk6[A|ƇscG{ + +|0cQsĹ,=Ę {Xr)NLb=@rSj)?E^AQ0zp@*@&py`FΪXyHV_"ôv+(xϰk PcX;zF Bx;8 РcPZ};S P1 +v-` rp(ŜH\R)I..o.d DB[ ["vaou@RC^9nr>z[鵽(7ҁE`K:/v:쎀@0Hi0#/ 9Gͱ b jY%Z)̱ c11ɱ6 D$]$7Jh9o[42a20: pw8m}1U8Yg9 =8)8 xÞCE*iٕ}5 |EB|8j,7Oxd5wt58Սŵ1tkcn(% :> [ڷmL YJ HjchVk#kԆm6w1d M ٍILo$mc :.=mr(5!5$q!Q2EEAw eF*3Kؾ[Cs**% #Bq4=2[c 5l1<~qL58*>ʱ˼zǁϪ1ތj8q8550Rcְ:h^~MD{Ħ8F1 0r^=iD5qi 'YDpMh+F(uAIg`qA8vI8b"q8=fКBklc߀Aǁ&1iN,C" ؘ1f8^f=4 p͗1f ZF0,ÜM22aÓQpdC2'CmÁ c֏11Оp16G1~SIEWh1ίbHD1Lbˮ&uNhh-tnJ5 l71 t7 c!`>({ Өc_8cJ'ppH`l c_r8Bo8_ G\8a +(_Vc²yAǽA {ARz^ Tv7w ^$¯8fqE!/BB8r̍@6)ϼZ`#pxE8r%ڔ]|œ2 -t9Leǎi0tGQ<:f`Y upMt]ux,fkxdH 0=qV,%b2,wy`10<d}+OhbЮp]c2=or1м֊Tiq嬘= +/@b+F|*@\(UcWPUx +1TDL TPA*ևDqPd & KOo ^NB)DMh1`d R)xD얂WxkL)ؽ xIU "M R_Q5q:0A!"(d!l"Y QՅC 2Y'*I,.+*XB9(Tb'!|[v'[f= =v<h>'N 0! i/9y&Pvc:cӁքr b2&0p&0U$$2AL=DŽSܐ6GzWn/gy@l7+70Z % +%rH WQv +!s*:r̞fCJF \BO$: 8Dmr8&AZ:Jy@)PIBzC4v$đuzH DA{xBBB I|KD2H"9- Ꮀt#1r#5]R`0Kfg?eD&2Q0E2dH],"hYN̤ԣ8D#t<;rDL #VR y*DLI ADL$>^୉CP"s8| :_IRC$-%ABIpa +&$%P.BVA#b) Q 9`a|/E^+A2Alr,A=)@b".S F uɼ@K@Duw0tm{cZ3 " 8Kn,x%A/\e | Ai.>&&&7Srr=01'?! ny@Ӱ<8\|LuDxXp'wc4C2t_J3dAZJ C$&V Ў߸*R^k5ODIl/>stream +ckV{s3VŢkVk߶~~=z^yrOu9oǶbޯ~]q:goO/hN[oWOsV=qUWq;}3V鶸z/ݛW{3kuX[Qi]j1>og]᫑\r'X(WX㸐ez3`o;M)P{kX⫵8q:Yqʔ:#Qن'Er>tb +/U/mLp:t8۶fj8۶mm$,m##۶M2ᴓm[3mmȌ3T@ Vc6HE`o۶Ŷmܶm+gN۶mmPm۶$h{N[b+[f-dd&Lrd+rS⪒JPi['hV%^f5WW)[}.PVPUMR&U3*IUOYZcʱ% AG+P-ujU $?bxV+ʕH^UY Ŗ ˱ʑ4 WWWTB +MUU# aTy^dzt}//I$I)1" yAQWqj,Us%PUU#D(%d(j/34[(6zU(-M'hod3+4g+ pAj''h*7Tf0/@ʔ5).QRk p4CWkg[ţV1m/w}wko=f=}nߺ7}Ww7z{nY{oǘߍ_}~{onhl1ޛgn3|k_wǞͷţw[<1cs--ݣޫ-#ʱE1qǙcK2oww}g[(M~qjG̓yXSyy[=w;-?n;}uWuj7/~loc=gousyƞs9^kV-Yo{u[m~﫫oϛg^mg;nq|5z㞷\kmy;os{Gx][?nn|-znvo{{sbܽqw5㋳b}{wc~GC꛻z^/=gs볾8_s}9zޭZ X0Wš~ *MC4F0CBbpΥU.,8ﰠre8M<B.Z$6 +6p* +L +Ѓ@" +:҉s' b;I@ dpKP<*SdL-Q .i!p$VĆpJ *1.9` +dx a+! tF 0g^.9<N#^FBF(qDEP ^ !%d> zDN& Ġ| +f@2K4Hs>hB~J3ZqN@@ mRb=2TgIx:uX߹FïCT)M?sTlVB(lNE_$Ȧ6@} +DB`%9͒ڻD@Gd Yf@!Fp&ExLXaT#.t_R r.*q,tη1PȖ xXD-%E3P6)B Qf ֏B"T]xB:*Tu +F ]بnln eJ$`'٢[bNXMr)a0WI19#:x1R )~D!')F q39FZ|$cI <3"HcCP "0D*6ć4`=Vģad`p<*of8)HsHx1z$#\p`.%D9 +9XtF]t5thʃаL 0jPG&,h1T!C4-\q.8 K"d +Z͔Z, b TXĄH BAItfQP,>e! ,0PF7/]Lpfdz, y$c1@'V<4n v~k]XMd`mj!"Hl]b㴗%q&#ö6ȴ`̬ m  HT:@l)ƃfD*&eAa6es$6ѝQc7޸,o'#o`7P(SՍm( d"h@<8m,m{unSIv0JDc5m۸n۶1 $bX?>:m۶6`m_l m,mN̶iJζͶQ|X-0BۀdYgTM+oOf(2m0 $*3eRzWq[X z@yShef>&h`pb]KZ! +RuZ1M@:,WLLE `V3:"0W4H*KAB@52H † &80 vH$*\ih50LC"{ ǶH۶ T!l6cJmi)ɶmcr) 6`8Mr(^qCCwЀRl pQ H8-KYh}CE@ʼnx*v-n‚* < JBP Ges5 :ʔ) eDy$I}#6DJBSdN ! &#TB<3/9Еl3xTAlfp1ec y`@ qbDž?ؼ ̉9/,<|>IRP#YY$X>*Ge9[;Z|vFhŭ\4hg]H2p4s7p"<(\Me#}x E2( d"72\MQ#jDy\D.npbˉZ'2񱚛@+yb%%V\L ՘`L|CB"CiN7 +&DB, +"=i!*`r#!2q +|*yfp1#sT|mtBa䋑 \ XTr6ss^f4ܡ"<0Q8qx +Ͽt2/gEwy,ۇqڄ3w=,Ł +@zpÛT`Eq.9"p{LFvRH:XE^N3U Em9zN+Si3 < &QH,p洮dqEY xD.De )AKFP3 :d~DyI#3.N&M"7! Ql NNdS r_FVC,hLb8hcK Hu\$$D%1V3؉5}/8M$Rd2R|*B"yx ,4X?4>x@HHS'Rp8, [)Y l$J=?X2GUcs&d |A1 +`QD (/LE эf貲0[jƸC o5**sNN[-O+cUfTBL1ތֻC!T+"ۓtZ|L%auH+#_Ȥ,-RVn1]єĔk +1:nLx&kWC;ncJFIQQ=.)}lq ( /"" H.]jVjv x3.,g5=4[R>F?iVi^=J~GW'&6mE^TuBr0?'}Ç~潨6lMsA-"NeF[_џTUWA9<@0sۨ7l|,TlϰeQ8ySٹlЪHSs?(?ZjQô"Fpf_^x^Xy,/ +hd2y%g-A^>beձb'tX&3Ek.WnC,4[VBl*rA=~tJϮ)]u̱'ʢ(]|02'" d9THx:yAq阸Óq~miܦ^g>ho3tG6L+7=w㒸^)KO?j|_ɓx@|ҝ6Sb{EY :=Plh |&xJ7M5@pOy +Ni0e[º~;amY%{]@5g{ju.&.,-[tܝs7WjѭrJ( +GઠtT8fNM⢋(>8` jstJُ=Xt:\q(@>@)t[tw4:KcDpO͚%Β <S+a'=(}Ua•SAGZd d} u5gN954 8%Ê`UMԬp{I eU;Qhł3 S1k-!KFsy8M}̠d1(4x|< CA vVpٴ\S 3Sne}kBv32-rab @BI]_s3уeM(6Ê%W{WR. 2F}d=* x{ &{U+@ +_&)EFP4 Jnc>"cH,c#LR\U;U^%|AJSĪq$| B[>~ߏX}K0~#&792 -߆CNDbZ^M8)`T #䲷.T +^c.|W!`bCO\"f+KM'~6u5 c>kvr7> adwɘ~uTA3Rk@|^Gz@omIFe6*pXIϪ0_k O-V6IE-XuzZHd@P:6Hy5`Ί?W O;LcH`KِVdSh @-ᚷwȐ8 +1(cDr>fsh>iD6D |ј/`G1uL  [vi@?ڎU | +Ƙ;\>*m`zCD@0I\àn"ſ퀛:ܨla:r> Sxہ)0Hqk Pc*7ۊ^5dVꣷFp=fD̫[6w@_s` $`F@0FࣳN˭m{mi39?x ĤbDϩ]h[?2S#@%+ +JK֐F%&X;2|ݐcGDSatG,)B~TBiizx`%Ы9cJQ[;Yf.C\>0AqPf3pa.!Q 4\er?tyB^91svd< {> ƒRpP]Rzئ0Ne8`]Aw4_b Pau,͘)>9bp舙C={L v_4QL&zK] ? &T©dPZ3޺Y 4f%{})V&pw%l`6 E}ԼXʒ3 {yIY0n~2 +┮:{l3qo {R/VpUt8b @$+~ng;=ëy N@ܲrzVez ܠ̵..2b-V<`ex:g/Fi[#P3w24 W%l\ԋ?>+TȖ]I` {uvS-pA'`ۚFN蕹?׬.Z0󺇚 7Ù)Ζ{@R%/F} {7ʁA@vt _bgTl\>%r.x"wfҔ- jo\\XwE̸ԫާ4FN}xࣶVA}l٩4ͷma:;y~m߿}^u_5%jڥ)*{q]耔\'K9~.A+ZWQLF81"19iYܴF\u=5TP+[ ?15`0Y%bW у]`y̳G u6m庘'~ 4V.Շoh)Ɋl+E]iaN->H) s 9.XZk?3@NFDo1q-Z:>1þ~ߓ΂K;lD{DDy1mn41s-9犉l!dGCϯp 4WV +_Q.dEq{E3Uz0Y[{tW- bTp^^5K9T1}b\o:q.nN[\C߇pճ${z+K||=!j`\6UX  Tb%-TT=eX9Băqzˋo̧BrPF]ap7aJNk'UH$;-1j46PfU {|FJc^Ih|qwc`b,5f L|`Y0Di>s($xfH'NY`0c>TdZҦ!ŢĸP9AJPC%%ePt&[魁mZ'§L{8?6DȽ3l乮: vl,Z/H +?[g A^Y1PY45$P'/'.UL ݃M9[bl@n[w~ b~jxޱ(іB0m$ Hk\KAIԏ-NjqPeǕnٌk@g$9"9X4 ui% +dRx6$o!W艍EhP )ؙoR@ KG\"n`$;p)*yq&9€j1yqdn\8*Q b',?TИ5z|]r&?#a6kG$3w wss_9 ̔p0 %مYQ]G}[\轧\&4bYI6KC2}pn=MHjN+)r +f=4t=d;9C4m7dIC.}4Z%/.N(*ǍJ#FE#}VROz!TD.(0fUКiuW}&LN+w\=hR&6*f@[\ǤԜ&ǑBްJ_%\1|'k_h'«G\>}1'W u&R (oMb +e`~4:SwKZ +Ep0宔54J|gY¥$AS\Hy;~Gu1+xSl\<ȓ3WҺi\D2*F]Mc67 A8m:6Cru@AUvm2(Nsvak rm-|M:P:JR@s3ek-/T6[_:TN*>"od}_DR]wlVfQ;UDOpЄo1~O?Z.W啤 pp_1!\RVO8\[iE T؄_8^0G|^ܝcjb>;ο|UwV*d+7"OV"\ßW;2àlNTfWU+x*rB]C=}\ wjJ#Ȱ2_cU\4!l?F*,t׎栜ɅV|֩"{?BRb7(9U~QQNM; :tdxT[q;LRӮqA' ]ZW9>I#fԊveĺKknD|I()8J1,O%:~>6I+T c{( +ڍ&l*)pnf#82RG)nÃL,v?!Iej\=2 zfc[z6HAEx}I)W0Wxpldx;Ӗ.]zTcH@^lqnm#{%S$g"٧c#u @v` ېLz>dxy"f[#DF@hZpt@wM7t+dazD#1%|cdσA0J֙^4(k|a^-]꟰G ?D,Z!$ qͭ\<w^bO0ul.9$0xg2fBCaf~ ce߫Zq}'*_`tqP 6$CVwKkP ;a)A˟ChA<|XPԖ,:jl#L'f1Ԯ 6k)tO| ` >ݏl2aE5E-UA'enC72GcCHVq\{.8#]sUsh❫)RY%30FrZe2vgtv\.bs8's7jBX>0sv 7(De3[W 0ngk]ŅhGl @™M nŠ!lq9Z|ܙ osfB> _6'.gZ& yj /}[OAr Fo~E_wBc!*[ +SfXw5 ֎`_@.$&+Fw\91Cy ,b9/9 +Z(*pdM i+HV7^򡟖: +Aq4RXOSɧ@7 dwq@Ⱥ[5Q )nѕ`֥IR'&P#-L Gǝng,n#&~? 2P̛J`,ͮ2SEB>pS@2Ƶf +E Pu1u3OH_R*<ֶg˪,G1CFGhb);H3;6fT7RpZ"įq_CdpƮUIw|Q2 +?yu|V3Fmc|u=8gfgK>[5FPnXm«! O(jzz>C0@;VZ_:R[p0ڃM[BkFe`x hT#%a5CY` CȡZF0RaHsIy;fSOZqaR|U6NGJ>U]$6hدi1* #S#__ib0}$PkIPР~LQse 쳿Ty>.3k‹?yg3aa@{Xo6FP9syN٧GDgrc?c߂ 3*^'–aggՕ@p Tr㍹qe[DQoz>4za !H> /rUT&)\i!x.qT[b #B$mV*vA:gW[R p.`Ѹq9,QaoS˽۱*L4aJkށ<j:t΄+:r$@[T~U_ZK\9 +YƘ6x +JfRN57mc_A9<JkmW[Q^Iހ95o+ކ ]P*,Nvqn-R\I٣b</[v=3 ٓ׀w; Hs껡B`z0":h>|F9Xh Xv"yL># ٶ#[%2^Ro/5)_Sq>~gFdjhъ42}, + u[l8!qൌC$vPua|`˘j-&TiO},L:#655 Wã +bSBod:f% +/v^\ +iJLUN9SHCf4o~Go1t|3}VU82С)7uߣpY ɅU9BC(EٍçcQ|kj-Ah/uD|2}esM{gLwzPmZG;/<??n@\ldB,#*vz5&qrx{&2F?W˿Lȿ :@!y4N~T8r/W WH'\.q;a)w=..l:U9m ab'MJ[}@LL!"+1xIOXtiqqzo,=J:::$ژH4iTVXvàq@@M%gdƮ}-N>$pևA._hY2d<[3"ÔުN) +xïO#{jNF=݂=b+yL9m(J䥮>m4a*B3A0`0$n03D@OjDP'z@&)1jDd!-j}#\56v}ptT96Z0LP`M4ʳ`U/DË0WH7kcTH${ 4pݵ/Lcޟ$kI S{܆(,^ءN|i'IG3 .9_PZh߯:~_3R / gۙRL -l5u]xiE0[mJp#d%H/<t/͑r_.a=Ÿ ^4셤H"SS8r_ƶq/V"Xp6U +yXat$dƯ2?IaZDP 8 +J77A:ERsj[j\_(!2zNy2)S&ʚȮ굪u`Y x?MO+"lo{un &|;J?dg.H>90 y\lMK/J4QS)r<VPbe ?uicݦ>W= X p<0E#PD2 hƍH{[<1 K!lߑv;-kg(T˟^CPքy}{%csUQڎZz $4VH hhpMD]&ѻ(5@0F~+}iSM+Iˇ ?z6 r1 +~=o8XksڝU)-u ;ۏY+(!CƬvɐNfЂ.Dܘw^E/a!pvԼ@+/YVg 'ICaYM8қYPZIWDTeP,rkHN=?uV,?QíL N|fc,hsj*WFψH򝴹 F +Λ]xNALR'dGAaY4U֮<7BÍޡ +jԶ @H;q`YxvNO FH{DH±@D:-Tn,D4.㰲4/kC:h|wiШ2_TA' ~gƏ> d|]gFy[au(b2ѣaL-@,ѯصqUWTJKW,Ε?:rnZUT$sJ:mD`=qߎ S1.X Rjdbr|y9|EZk_$.2lٸF ]0KZ 8/xi4 H=]_$^w>`_2JZٚtS%rִ\Sܱyn|aJ JuX6|H( SOG`:^ 'X?;=Fi߇](Z^ViVgE>/՘\* )8 H@6Ml׮[9fE6H a> "&7ńJǣ[TdJօT'N/(0O}apd!U% ǂMCΫj17h1A7 i9[G(tpPt˝~S#vC9Y/ +bjbH9PG @QC Ξ{C!5 3JI +<i {d)MO,vd6a=`A]ZabFE /S;~S}U,ާ\#2vv ݶHJP9 $>n2%.%L|6&a'Vȴy5I2 x9ºļ3w NjI -91r0FهfJ( +Fpt1O~+]yW"vֿ288Pgɷ-4PuzpTd[T_ZO>eQ J ?VBţnK*N⭮GKvFYtt{X=_(WA ] +Vd8:%jX1{}Ybec?c+U>vg0 VԺ9}?p&u75qLpZh9ݤЗIƺÅ4},0%{z +}ӬV{p\c[ت!g?Ռq|x9Q ^K̳+9ӍtP1/͝c,TFTC[t7 bpIW59yTd8]dLݓOC9PPuMJK) i7J;fxbIJ2$@6#@H}s4zr{6R4zѤ2RAèݽ?@3Vq]$D. + v 4Ee\,qsz%!h-^]J!cZMz1i"KtPel^aCj_ I֤@`NZ c(=v_VuK6usZb4lm0k5"Lc"|0MՊ䈭 A&4a%P,dڱ⨲aN +p։jbY9잝C#9RbN7╉*Y̪!"r:RCN2ZHQVga(+K '"F|$wo~Ɲb՚Zעi<05٪ljgr[N/FxF"[ka%zs0;]wgf}Zڢ0 @0<)Wxg00usPsp$G3+,D/-*s|Czd# +i +7J T2JFP{.%B\b aU1=qe0=i..xQQ ぽ˾<>d>dqt܎ۉcx &c/6hiXUpX9jw=X;N[`~$M-0b$iJ$#9)ju'xsbԦ—}]CD:I^PyA-S@hfF1䝦i<;m䍑)ĆiyGn;ZE46C%5ˉ͒Jкwj-o  c:8Mh[rWfٿ"N3Q@_Jykz,W[-RMs| 7G8q + #:nn&`jϥו9^w.0ʰ sZ  S%>+O-wpX!YJJ#7d<崇CoyF.+䊉g9^i^ª!eizQ-T!DCrf(5 q~og)V$Xͤ[I}L @/F7/}.n *T&C%5ۯkk4q,+p~X@Nc_#Io3܎z=dzMۜ&!Cof=D5x#.\H"l)l7@.'XJz Ԏ2xf +f8dŒYuLj;S`R>+Y$sX* +`ʃ@P0{M^֫N"3LQ.^˳ A ހ\5Ҽbc`\A ܯy!D6(?x̂iTTxޅI>> {i P-Đ8lY qPg +HooٓpkRMx؊j8 Q [Rkg0*;; h,B!w /ě 4/"_j]# ɝG_h]Ldq`% r%_*N|tį>Lk')9;RKN?ֿw>i(up#v.,ƯaaEAxsay#{#_Z>ZX5܄VQ/pe2'8 wyW:v>27pwM?&9xPg[EqEZϩP'`,|5Ib V3ޜ_dG` KoMQq9T/ ^-qIzU d="~FIڜ_ E{WC&ś4!h UsqFoW vw;e`1Rm9FUjhخQc`PU`r$|!ݧda4r0j˦jA9Y+$- ڵ_PH/qZ'{&|D lke1 @@S qQp0j7)"&҈h]$+uxpr"#EyD;9ި=dVSA}XzF)GhxW|CP@p,nʇH/}tyCрއ F/h̡\HI ++FX}=R5i-Z 9S +$$jH2ZI4u}}XM)t>0eG5eE'\WslC8 z@W)0m~5񽉘&røj8rXo8Qc;MRrr9e[4E EC }BӶ^a#IO74j8]i hia3I\Jg -׭XՕ@,p^/c,5bKOHK= Hl` ދZ~1[:ĕ^:|+MZDxȽˠ-J"5[EnTՈ$DQ|ޕ^oFIƯW4_ lt."x~e^z1TNEt׻K6&'k&{sv^7xN[+6'^1CDԢE2R"w/l TJx+zדTrwNJew^aW.<}:# +9nQO4qޣ7,=wQ@◂Ť{Jtr*NW/Mş$]œVb,R^O5\$P۰9L$Vhx:1›iCb$#A931EaϘI^px' V w./h!'[ُ Fms)q/-V*Tv?J?}fPJgd}M 'ls@ZV.Ym', "I2k4lup0AiyIv͔n4MDMdr(=ES4,M-g(z%lwn[_SkӧYO-?Os'7:PSc &jzw_e?ύ'k{?'{\G]9H"r|prD"%X<Kzoz>z;iL~q H(K'\ V6a&6\N 9}x՛88m)K!BmD<:(#-Z,g<V A_ E.,i&q!)MۈƑi)b,Df&~L"qa*RVZ;iCERuXΒ 9M3R"'O^DBl'2`J`~14 i)%ѸAY*d$,pI\ǚ-cb%m ;84$.2+-1WzQWᬻr L ShKK6 j/rWCLeREze(({ږS$U8Js{cD uN^W+-6#^Ĉ+ &U2Rnkԃշ 4>qk+|K nG\ql KKjCuG K5N X|.N/."&ي{AALA)bSm+`fVy߮*EYB.kb1"1f +%7!'4hPI\4_Jj$wHD_1| + l1 jg`gA#_s>q%?vq4ON,6h-*t9 K߶DV!WfۚK7 QDB?Eh)V `\JV VCARYHdIXfdCx+)1vz ڳcI,x d 9 +H9#l@, XlU +(6E*{*t%iXDVY-Wx[ 4vV:pƲAk8ЖvUOpNom1YU(!H`)@jFGart܈kSa8Sb|ű{tMb**0OB.Xw؃Hku(͑t֞yZ[f3FEgԎu+8, +{%Hj13)e{W̲/ϘuGfoQ ̳6Q]0dHP֐GSft=aɍ*ʬgf[Xe&b=$\IOCQ{z1N&Ǹ&)eכ_V2+qaa#/1H9 _QSi#h!|mᬥ<ت60P$e徂56_QԔv[axtC"$(a2^e +wNVICf|D$FhtO욡בEL.4Wn r.&eڀV_M@,?v[%el!`MӈZ{*?ʬ%2`GR aչ A~(]J%08R9zr8ik-~`6" B5#M*YxŠ0֖ +(``kUWdxxj`ˆ%ȜIWuο5 gQiA3e]HN@׿ɘ$G PATej@30<OASkM@('Dɞ "r>o7}\nr Ko=r+h~Ip@G#R|V +=NgdN}S'ǘTgRErz<#~;b! Ɇ#= csYOt=F{Mg\H4]Ji~)l%2fɨ +w ,$|B@blwnAA$vGf5ՈG> +< +r4)V 9GOg8<=ZIHMDn*eeG$VH;̮rP<'-^?W eAm +%-\!ȥ̺̋"BG1sXz!-yJtPA\QcX>/lLp4XHR2hRw)nԇʋ8&bNM!~]QcB(8*VT,+ΐtrE <#F"T7QTX$jUZyq\2ȧ'Ouc0:R!^;%|~s(ת__;Z8~ :%Yik +wD.o_}(mP. h. *!XO6CE۲vr_sk$U74-YT>TBCVgl(!#uX*`$ҳK4miF_u +j>Û󷑹X;w9paw~)l7)7Z#TOPE`[.dW?FY1p%Q#md/9 +ͥ;O:P9 +ldJ.*k3 eb8^8fRǒ1A;lx;>iqU]-Qha>48Xմ ?[X8$&( T, S~0@ އz!ޱL 7]eii@]۶Chnuм?riI/lWWzzK@&MUyjlV--kß15nmݔ'7Z(vP]x%`JJiGjTDd dgK_'bW4w5W!!^]#z>@{9Ct GLYLQXdg +C6g'Bc?Tvxv>>_6w(;+84 gMI9ȷ!A2IA#bD(hJq' s#}ПK&܏2~7x *2нBh̓"W|$c.&,ƭ6~EHr}9 `4jH#x -X&+7nUǠ1' =dFM; ' :.EO'az;ZqF-b*B)/Py&oOјހ[uL QemVoVQ,;$.t%UOYZ#3;k|8P?=Qb^r"DD +XyE. [ +՟9㝏O'_cuXt9pZTI=HX1۾[홫4"\d6X%QM+(aL' ( #ф& ,rB;X&?;U]Ni]ί +ĪP+WE51dۃjHCDg^BO9|; 5TB8{h!.]&&P'IW hrYReDB;.$aܜ +SI/&FnXob F>ШT^:& &EU+~w3]+[I)G@5d{4K%!XQ?'JVsg F]|powEn(qe{m'$x Q@Vd?_dжmdUHj(j{ ,;9k}3T).}ŪM `O~t`j2HhA8aU9)IKa9 -75o .ehp{}=e{lPDݭ3u]}0P~G#ys^S}Up3o 'cݶ2ar IDCrU80Q1)F 2Yknq:^ .VΜ@ yOl&]ȽrBoFZHՒ Y%$$$|ieph%6W/H'tdq\J:e#[(.i*~_B$Qu%GÈO +IŴ{R +B8k <KAA8}Jq6VDv{"{zSUsRv:uf"+~o."Nh?^5I0%f8]gl]"?3ՆApƘ *\v"BmDq +S@_B`8=#£/G*c/܁^Um732vr0Q.2'(ӷ:gL<ח˄ؚnmQh[U~.jV 4: ?<^qoFݸXWqv0W]qyڶ4V.K:Ʊ!BAXdrXLP%ϴ_/hpO}\f9ϲUd&g`nv;Jm0.t*Q&TհRƍ؁Xde9E1 *Ak);q%g(9.F?2=St֑(| 8F@|KȚ7!5OP3v"`QZo6P`w@ZY5ۯڹX49~Ik~ȁwzs(FK.mܑYz479Cy e!+ Yχ0*d~)rD[{/*]4~*"+y+UdM]Y NYP;iFYpXAP-l!Ieb_8WA2k[EnP}B9"a˼fA|@z65N~ P9+nxvX[@UmCVD4/(z1+Tr ~Fq"yNopq\6bpX8_zV +( hWp\ueg,HJoz%((/U^~s|QڇS7lO4Q' +G+}8O_i~Y5lD +[uS_ Oo{Ͳ-X&Yp#̺iG?ɷ@(2R=њ۫5%!Io1\/M=i)k+M0<?fYQDXw%#P4/PwYrla2Yn4@q6b&VEHQjrgrӞzWlcZ>{Pַkl{l"Y13KI(<`7$7t⼊b*.SMU8D KRor!6*/cN.A LdVJsR ޚAzwnZyTEH(3"6T?+d)>,^a p9s{ֻ=K> P!W`*睍$uD1qd‚hKw+?X(FM{=3?pFy8 *d^z +SNJ#)`p&\Ю~`QGyJ E +JDAAIJ[ҽCR@we" 1ӹtvB:˿vql-RjΓ[9K2ݜk +)S\/9X,}.|.ATDpf#{WU j Ǩ悶٧GŞ(f6{jTMhzb` 22NggW͒I#m x- >4ePg$H+Zи[/FkBF^u0_[T7_p:HiDJ(>˄TVr)E@ŖnFN.rC3/ƽ&va9]7${zF(YNe~LF,&yhZ+rvEf;YGNq4qcY/wA)3Lg 1uԌA`? +u|1~_ @t!9RJ2'{Xz$;a,.R +5c4Z\}MN$d#r|@,<$gCujaCm7:ad~i@e7L~[Q05&!V"/Û;b%F2eUشEϜW9uW3PoCZigx6^Gvo +rZ;눇X%Oc=# +AFe#TVe捥siD/˞_W+i=!!oW<hJ5P+$#cڅ; v֥`:% ^::Z:P ǥ@x&BnM`Q9fBke#p!uli2kEݫ-mm#|\ߣ#W W#* Rߖ J9`ES584nBF<68xs/ע/lu ɰtIX@ͳfRs7ucE>WĉALx<<$QW1f +U{>1x?ݑ2^ `kP}Ԝ.{][پSUv#Ac5@mVQrЌo =;=Q#Mnju_^5QWm;F%E-%޶ +:XAglMk<an>Re*/fDԨ; 0\>: d2{OpX ]i}'t % xJ™oR\!Fd+L b3r"pce/|tXnx#!udǮ1}aR +M.atnr)1 rdu6'̃"'vifweIAډ&Xbdy~[<J&}RU˒pZ&&p` z][<xl >bM>sR{ ޵Rv#̫̈́I=^\F0뎴vTp d7Bۥttn{"?/L;mtc98̡az%%ɖboLBge cM8H +1˜޲X^m@P  ݕ0d{8sm`NΪt얬Rᆍi%U#K)~tՃ%4 +w +1e,x'rDI+˞q}̞;FEhD5 endstream endobj 24 0 obj [23 0 R] endobj 38 0 obj <> endobj xref +0 39 +0000000004 65535 f +0000000016 00000 n +0000000147 00000 n +0000039285 00000 n +0000000000 00000 f +0000039336 00000 n +0000000000 00000 f +0000052193 00000 n +0000000000 00000 f +0000000000 00000 f +0000000000 00000 f +0000000000 00000 f +0000000000 00000 f +0000000000 00000 f +0000000000 00000 f +0000000000 00000 f +0000052266 00000 n +0000052440 00000 n +0000053596 00000 n +0000119184 00000 n +0000000000 00000 f +0000041595 00000 n +0000041917 00000 n +0000041408 00000 n +0000155045 00000 n +0000039778 00000 n +0000049510 00000 n +0000049397 00000 n +0000040562 00000 n +0000040847 00000 n +0000040895 00000 n +0000041479 00000 n +0000041510 00000 n +0000048341 00000 n +0000042251 00000 n +0000042494 00000 n +0000048628 00000 n +0000049545 00000 n +0000155070 00000 n +trailer <<127DBCABF56046DFB7D3ECC9A5BC80AE>]>> startxref 155256 %%EOF \ No newline at end of file diff --git a/doc/source/figures/weight.pdf b/doc/source/figures/weight.pdf new file mode 100644 index 0000000..c8310c8 Binary files /dev/null and b/doc/source/figures/weight.pdf differ diff --git a/doc/source/index.rst b/doc/source/index.rst new file mode 100644 index 0000000..b75edab --- /dev/null +++ b/doc/source/index.rst @@ -0,0 +1,31 @@ +.. chromo documentation master file, created by + sphinx-quickstart on Thu Sep 17 12:55:50 2020. + You can adapt this file completely to your liking, but it should at least + contain the root `toctree` directive. + +Welcome to chromo's documentation! +================================== + +.. toctree:: + :maxdepth: 2 + :caption: Contents: + + Code Overview and Structure + + Polymer Models + + Field Theoretic Treatment of Interactions + + Monte Carlo Simulations + + Brownian Dynamics Simulations + + References + + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` diff --git a/doc/source/mc_sim.rst b/doc/source/mc_sim.rst new file mode 100644 index 0000000..91932cf --- /dev/null +++ b/doc/source/mc_sim.rst @@ -0,0 +1,7 @@ +.. _mc_sim: + +Monte Carlo Simulations +======================= + +9/18/20 - Joe will update this page with an overview of the Monte Methodology (with schematic of moves and methods) + diff --git a/doc/source/poly_models.rst b/doc/source/poly_models.rst new file mode 100644 index 0000000..9077869 --- /dev/null +++ b/doc/source/poly_models.rst @@ -0,0 +1,334 @@ +.. _poly_models: + +Polymer Models +============== + +We introduce the polymer models that can be used in the chromo simulation. +These models are designed to capture different physical effects. +They are all based on a common representation of the +chain as a discrete set of beads. +However, each model uses a subset of the geometric variables in their +calculation of energy, force, and torque, and the simulation automatically +adjusts the underlying calculations based on the model. +Here, we define the general chain representation that includes the most +complete geometric variables for the chain, and each polymer model utilizes +a subset of these variables (defined for each model below). +Each model also has the option of creating a linear or ring polymer, +which simply requires the model to add a bond between the first bead and the +last bead of the chain. + +We consider a polymer with :math:`n_{b}` number of beads in a single chain. +The polymer chain is represented by the +bead positions +:math:`\vec{r}^{(0)}, \vec{r}^{(1)}, \ldots, \vec{r}^{(n_{b}-1)}` +and orientations +:math:`\vec{t}_{i}^{(0)}, \vec{t}_{i}^{(1)}, \ldots, \vec{t}_{i}^{(n_{b}-1)}` +where :math:`i = 1, 2, 3`. +Since the orientation triad forms an orthonormal basis, the dynamics must maintain the following conditions: + +.. math:: + \vec{t}_{i}^{(n)} \cdot \vec{t}_{j}^{(n)} = \delta_{ij} + +Each model represents a chain whose total length is :math:`L`. +For linear polymers, the chain is discretized into segments of length +:math:`\Delta = L/(n_{b}-1)`. This implies that the end-to-end +vector is given by :math:`\vec{R} = \vec{r}^{(n_{b}-1)} - \vec{r}^{(0)}`. +For ring polymers, the chain is discretized into segments of length +:math:`\Delta = L/n_{b}`, which includes the bond between the last +bead at :math:`\vec{r}^{(n_{b} - 1)}` and the first bead at +:math:`\vec{r}^{(0)}`. + + +.. table:: + :widths: 30 30 30 30 30 + :align: center + + +-----------------+--------------+------------------+-------------------+-------------------+ + | model code | 'track_norm' | 'bond_contraint' | 'ring' | 'link_constraint' | + +=================+==============+==================+===================+===================+ + | 'gauss_chain' | False | False | True/False | False | + +-----------------+--------------+------------------+-------------------+-------------------+ + | 'wlc' | False | True | True/False | False | + +-----------------+--------------+------------------+-------------------+-------------------+ + | 'sswlc' | False | False | True/False | False | + +-----------------+--------------+------------------+-------------------+-------------------+ + | 'twlc' | True | True | True/False | True/False | + +-----------------+--------------+------------------+-------------------+-------------------+ + | 'tsswlc' | True | False | True/False | True/False | + +-----------------+--------------+------------------+-------------------+-------------------+ + + +Polymer chain model A: flexible Gaussian chain +---------------------------------------------- + +The flexible Gaussian chain captures the behavior of a flexible polymer chain that +is representative of a Gaussian random walk in the absence of additional interactions. +The chain is defined by the bead positions :math:`\vec{r}^{(n)}`, and the +bead orientations are not utilized in this model. + +We define the polymer energy function + +.. math:: + \beta E_{\mathrm{poly}} = \sum_{n=0}^{n_{b}-2} + \frac{3}{2 \Delta b} \left( \vec{r}^{(n+1)} - \vec{r}^{(n)} \right)^{2} + +where we define :math:`\beta = 1/(k_{B}T)`, and the Kuhn length +:math:`b` defines the statistical segment length of the polymer chain. +This energy definition is valid for a linear chain. +However, the ring representation only requires the upper limit of the +summation to be changed to :math:`n_{b} - 1`, and we note that +the condition :math:`\vec{r}^{(n_{b})} = \vec{r}^{(0)}` connects the +chain ends into a ring. +In this model, the bead orientations :math:`\vec{t}_{i}^{(n)}` do not +contribute to the energy and are not evolved in the simulation. + +From the polymer energy, we determine the force on the nth bead +is determined from :math:`\vec{f}^{(n)} = - \frac{\partial \beta E}{\partial \vec{r}^{(n)}}`, +where we scale the force by the thermal energy :math:`k_{B}T`. +This results in the force on the nth bead to be given as + +.. math:: + \vec{f}^{(n)} = \frac{3}{\Delta b} \left( \vec{r}^{(n+1)} + - 2 \vec{r}^{(n)} + \vec{r}^{(n-1)} \right) + +for all beads within the interior of the chain. +This expression also applies to the end beads :math:`n=0` and +:math:`n=n_{b}-1` for a ring polymer, if we note that +:math:`\vec{r}^{(-1)} = \vec{r}^{(n_{b}-1)}` and +:math:`\vec{r}^{(n_{b})} = \vec{r}^{(0)}`. + +For the ends of a linear chain, the end-bead forces are given by + +.. math:: + \vec{f}^{(0)} & = & \frac{3}{\Delta b} \left( + \vec{r}^{(1)} + - \vec{r}^{(0)} + \right) \\ + \vec{f}^{(n_{b}-1)} & = & - \frac{3}{\Delta b} \left( + \vec{r}^{(n_{b}-1)} + - \vec{r}^{(n_{b}-2)} + \right) + +Initialization of the flexible Gaussian chain in the absence of additional interactions +(e.g. excluded-volume interactions or confinement) +is performed by selecting the bond vectors from a Gaussian distribution with +a variance that is given by + +.. math:: + \langle + \left( + \vec{r}^{(n+1)} - \vec{r}^{(n)} + \right) + \left( + \vec{r}^{(n'+1)} - \vec{r}^{(n')} + \right) + \rangle = \frac{\Delta b}{3} \delta_{nn'} \mathbf{I} + +which leads to a mean-square end-to-end distance + +.. math:: + \langle + \vec{R}^{2} + \rangle = + \sum_{n=0}^{n_{b}-2} + \sum_{n'=0}^{n_{b}-2} + \langle + \left( + \vec{r}^{(n+1)} - \vec{r}^{(n)} + \right) \cdot + \left( + \vec{r}^{(n'+1)} - \vec{r}^{(n')} + \right) + \rangle + = \sum_{n=0}^{n_{b}-2} + \sum_{n'=0}^{n_{b}-2} + \Delta b \delta_{nn'} + = \Delta b (n_{b} - 1) + = L b + +which is consistent with the solution for a Gaussian random +walk polymer with length :math:`L = N b`, +where :math:`N` is the number of Kuhn lengths in the +chain. + +Polymer chain model B: shearable, stretchable wormlike chain +------------------------------------------------------------ + +We consider the shearable, stretchable wormlike chain potential, given by + +.. math:: + \beta E_{\mathrm{poly}} = \sum_{n=0}^{n_{b}-2} + \left[ + \frac{\epsilon_{\mathrm{b}}}{2 \Delta} \left| \vec{t}_{3}^{(n+1)} - \vec{t}_{3}^{(n)} - \eta \Delta \vec{r}_{\perp}^{(n)} \right|^{2} + + \frac{\epsilon_{\mathrm{\parallel}}}{2 \Delta} \left( \Delta \vec{r}^{(n)} \cdot \vec{t}_{3}^{(n)} - \Delta \gamma \right)^{2} + + \frac{\epsilon_{\mathrm{\perp}}}{2 \Delta} \left| \Delta \vec{r}_{\perp}^{(n)} \right|^{2} + \right], + +where :math:`\Delta \vec{r}^{(n)} = \vec{r}^{(n+1)} - \vec{r}^{(n)}` is the bond vector, +:math:`\Delta \vec{r}_{\perp}^{(n)} = \Delta \vec{r}^{(n)} - (\Delta \vec{r}^{(n)} \cdot \vec{t}_{3}^{(n)}) \vec{t}_{3}^{(n)}` +is the perpendicular component of the bond vector to the tangent vector. + + + + +Polymer chain model C: shearable, stretchable wormlike chain with twist +----------------------------------------------------------------------- + +We consider a closed ring polymer, where the :math:`n_{b}` bead is the same as the zeroth bead. +We consider the shearable, stretchable wormlike chain potential with twist, given by + +.. math:: + \beta E_{\mathrm{poly}} = \sum_{n=0}^{n_{b}-1} + \left[ + \frac{\epsilon_{\mathrm{b}}}{2 \Delta} \left| \vec{t}_{3}^{(n+1)} - \vec{t}_{3}^{(n)} - \eta \Delta \vec{r}_{\perp}^{(n)} \right|^{2} + + \frac{\epsilon_{\mathrm{\parallel}}}{2 \Delta} \left( \Delta \vec{r}^{(n)} \cdot \vec{t}_{3}^{(n)} - \Delta \gamma \right)^{2} + + \frac{\epsilon_{\mathrm{\perp}}}{2 \Delta} \left| \Delta \vec{r}_{\perp}^{(n)} \right|^{2} + + \frac{\epsilon_{\mathrm{t}}}{2 \Delta} \left( \omega^{(n)} \right)^{2} + \right], + +where :math:`\Delta \vec{r}^{(n)} = \vec{r}^{(n+1)} - \vec{r}^{(n)}` is the bond vector, +:math:`\Delta \vec{r}_{\perp}^{(n)} = \Delta \vec{r}^{(n)} - (\Delta \vec{r}^{(n)} \cdot \vec{t}_{3}^{(n)}) \vec{t}_{3}^{(n)}` is the +perpendicular component of the bond vector to the tangent vector. + +The twist angle :math:`\omega^{(n)}` gives the +local twist deformation of the chain. +Geometrically, this is defined by the relationship + +.. math:: + \left( 1 + \vec{t}_{3}^{(n)} \cdot \vec{t}_{3}^{(n+1)} \right) \cos \Omega^{(n)} & = & + \vec{t}_{1}^{(n)} \cdot \vec{t}_{1}^{(n+1)} + + \vec{t}_{2}^{(n)} \cdot \vec{t}_{2}^{(n+1)} \\ + \left( 1 + \vec{t}_{3}^{(n)} \cdot \vec{t}_{3}^{(n+1)} \right) \sin \Omega^{(n)} & = & + \vec{t}_{2}^{(n)} \cdot \vec{t}_{1}^{(n+1)} - + \vec{t}_{1}^{(n)} \cdot \vec{t}_{2}^{(n+1)} + +where :math:`\omega^{(n)} = \Omega^{(n)} + 2 \pi m^{(n)}`, +and :math:`m^{(n)}` gives the number of additional integer turns +of twist within the +nth segment. +We write a differential change in :math:`\omega^{(n)}` as + +.. math:: + \delta \omega^{(n)} & = & + \frac{\vec{t}_{1}^{(n+1)} \cdot \delta \vec{t}_{2}^{(n)}}{ + \vec{t}_{1}^{(n)} \cdot \vec{t}_{1}^{(n+1)} + + \vec{t}_{2}^{(n)} \cdot \vec{t}_{2}^{(n+1)} + } - + \frac{\vec{t}_{2}^{(n+1)} \cdot \delta \vec{t}_{1}^{(n)}}{ + \vec{t}_{1}^{(n)} \cdot \vec{t}_{1}^{(n+1)} + + \vec{t}_{2}^{(n)} \cdot \vec{t}_{2}^{(n+1)} + } + \nonumber \\ + & & + + \frac{\vec{t}_{2}^{(n)} \cdot \delta \vec{t}_{1}^{(n+1)}}{ + \vec{t}_{1}^{(n)} \cdot \vec{t}_{1}^{(n+1)} + + \vec{t}_{2}^{(n)} \cdot \vec{t}_{2}^{(n+1)} + } - + \frac{\vec{t}_{1}^{(n)} \cdot \delta \vec{t}_{2}^{(n+1)}}{ + \vec{t}_{1}^{(n)} \cdot \vec{t}_{1}^{(n+1)} + + \vec{t}_{2}^{(n)} \cdot \vec{t}_{2}^{(n+1)} + } \nonumber \\ + & & + - \left( + \frac{\vec{t}_{2}^{(n)} \cdot \vec{t}_{1}^{(n+1)} - \vec{t}_{1}^{(n)} \cdot \vec{t}_{2}^{(n+1)} } + {\vec{t}_{1}^{(n)} \cdot \vec{t}_{1}^{(n+1)} + + \vec{t}_{2}^{(n)} \cdot \vec{t}_{2}^{(n+1)} } + \right) + \left( + \frac{ + \vec{t}_{3}^{(n)} \cdot \delta \vec{t}_{3}^{(n+1)} + + \vec{t}_{3}^{(n+1)} \cdot \delta \vec{t}_{3}^{(n)} + }{1 + \vec{t}_{3}^{(n)} \cdot \vec{t}_{3}^{(n+1)} } + \right) + + +With this development, we write the torque vectors as + +.. math:: + \vec{\tau}_{1}^{(n)} & = & \frac{\epsilon_{t}}{\Delta} \omega^{(n)} + \left( + \frac{\vec{t}_{2}^{(n+1)}}{ + \vec{t}_{1}^{(n)} \cdot \vec{t}_{1}^{(n+1)} + + \vec{t}_{2}^{(n)} \cdot \vec{t}_{2}^{(n+1)}} + \right) + - + \frac{\epsilon_{t}}{\Delta} \omega^{(n-1)} + \left( + \frac{\vec{t}_{2}^{(n-1)}}{ + \vec{t}_{1}^{(n-1)} \cdot \vec{t}_{1}^{(n)} + + \vec{t}_{2}^{(n-1)} \cdot \vec{t}_{2}^{(n)}} + \right) + \\ + \vec{\tau}_{2}^{(n)} & = & - \frac{\epsilon_{t}}{\Delta} \omega^{(n)} + \left( + \frac{\vec{t}_{1}^{(n+1)}}{ + \vec{t}_{1}^{(n)} \cdot \vec{t}_{1}^{(n+1)} + + \vec{t}_{2}^{(n)} \cdot \vec{t}_{2}^{(n+1)}} + \right) + + + \frac{\epsilon_{t}}{\Delta} \omega^{(n-1)} + \left( + \frac{\vec{t}_{1}^{(n-1)} }{ + \vec{t}_{1}^{(n-1)} \cdot \vec{t}_{1}^{(n)} + + \vec{t}_{2}^{(n-1)} \cdot \vec{t}_{2}^{(n)}} + \right) \\ + \vec{\tau}_{3}^{(n)} & = & + \vec{\tau}_{b}^{(n)} - + \vec{\tau}_{b}^{(n-1)} - \eta \left[ + (\Delta \vec{r}^{(n)} \cdot \vec{t}_{3}^{(n)}) \vec{\tau}_{b}^{(n)} + + ( \vec{\tau}_{b}^{(n)} \cdot \vec{t}_{3}^{(n)} ) \Delta \vec{r}^{(n)} + \right] + \nonumber \\ + & & + - \frac{\epsilon_{\parallel}}{\Delta} + \left( \Delta \vec{r}^{(n)} \cdot \vec{t}_{3}^{(n)} - \Delta \gamma \right) \Delta \vec{r}^{(n)} + + \frac{\epsilon_{\perp}}{\Delta} + (\Delta \vec{r}^{(n)} \cdot \vec{t}_{3}^{(n)} ) \Delta \vec{r}_{\perp}^{(n)} + \nonumber \\ + & & + +\frac{\epsilon_{t}}{\Delta} \omega^{(n)} + \left( + \frac{\vec{t}_{2}^{(n)} \cdot \vec{t}_{1}^{(n+1)} - \vec{t}_{1}^{(n)} \cdot \vec{t}_{2}^{(n+1)} } + {\vec{t}_{1}^{(n)} \cdot \vec{t}_{1}^{(n+1)} + + \vec{t}_{2}^{(n)} \cdot \vec{t}_{2}^{(n+1)} } + \right) + \frac{ + \vec{t}_{3}^{(n+1)}}{1 + \vec{t}_{3}^{(n)} \cdot \vec{t}_{3}^{(n+1)} } \nonumber \\ + & & + + + \frac{\epsilon_{t}}{\Delta} \omega^{(n-1)} + \left( + \frac{\vec{t}_{2}^{(n-1)} \cdot \vec{t}_{1}^{(n)} - \vec{t}_{1}^{(n-1)} \cdot \vec{t}_{2}^{(n)} } + {\vec{t}_{1}^{(n-1)} \cdot \vec{t}_{1}^{(n)} + + \vec{t}_{2}^{(n-1)} \cdot \vec{t}_{2}^{(n)} } + \right) + \frac{ + \vec{t}_{3}^{(n-1)}}{1 + \vec{t}_{3}^{(n-1)} \cdot \vec{t}_{3}^{(n)} } + +where + +.. math:: + \vec{\tau}_{b}^{(n)} = + \frac{\epsilon_{b}}{\Delta} \left( + \vec{t}_{3}^{(n+1)} - \vec{t}_{3}^{(n)} - \eta \Delta \vec{r}_{\perp}^{(n)} + \right) + +The force on the nth bead is given by + +.. math:: + \vec{f}^{(n)} & = & + -\eta \vec{\tau}_{b}^{(n)} + \eta ( \vec{\tau}_{b}^{(n)} \cdot \vec{t}_{3}^{(n)} ) \vec{t}_{3}^{(n)} + +\eta \vec{\tau}_{b}^{(n-1)} - \eta ( \vec{\tau}_{b}^{(n-1)} \cdot \vec{t}_{3}^{(n-1)} ) \vec{t}_{3}^{(n-1)} + \nonumber \\ + & & + + \frac{\epsilon_{\parallel}}{\Delta} + \left( \Delta \vec{r}^{(n)} \cdot \vec{t}_{3}^{(n)} - \Delta \gamma \right) \vec{t}_{3}^{(n)} + - \frac{\epsilon_{\parallel}}{\Delta} + \left( \Delta \vec{r}^{(n-1)} \cdot \vec{t}_{3}^{(n-1)} - \Delta \gamma \right) \vec{t}_{3}^{(n-1)} + \nonumber \\ + & & + + \frac{\epsilon_{\perp}}{\Delta} + \Delta \vec{r}_{\perp}^{(n)} + - \frac{\epsilon_{\perp}}{\Delta} + \Delta \vec{r}_{\perp}^{(n-1)} + diff --git a/doc/source/references.rst b/doc/source/references.rst new file mode 100644 index 0000000..154c251 --- /dev/null +++ b/doc/source/references.rst @@ -0,0 +1,43 @@ +.. _references: + + +References +========== + +.. Spakowitz Lab references + +.. [Spakowitz2004] + Spakowitz, Andrew J., and Zhen-Gang Wang. "Exact results for a semiflexible polymer chain in an aligning field." Macromolecules 37.15 (2004): 5814-5823. + +.. [Spakowitz2005] + Spakowitz, Andrew J., and Zhen-Gang Wang. "End-to-end distance vector distribution with fixed end orientations for the wormlike chain model." Physical Review E 72.4 (2005): 041802. + +.. [Spakowitz2006] + Spakowitz, A. J. "Wormlike chain statistics with twist and fixed ends." EPL (Europhysics Letters) 73.5 (2006): 684. + +.. [Mehraeen2008] + Mehraeen, Shafigh, et al. "End-to-end distribution for a wormlike chain in arbitrary dimensions." Physical Review E 77.6 (2008): 061803. + +.. Simulation methodologies + +.. [Pike2009a] + Detcheverry, François A., et al. "Monte Carlo simulation of coarse grain polymeric systems." Physical review letters 102.19 (2009): 197801. + +.. [Pike2009b] + Pike, Darin Q., et al. "Theoretically informed coarse grain simulations of polymeric systems." The Journal of chemical physics 131.8 (2009): 084903. + +.. [Hinch1994] + Hinch, E. J. "Brownian motion with stiff bonds and rigid constraints." Journal of Fluid Mechanics 271 (1994): 219-234. + +.. Wormlike Chain references + +.. [Kratky1949] + Kratky, Otto, and Günther Porod. "Röntgenuntersuchung gelöster fadenmoleküle." Recueil des Travaux Chimiques des Pays‐Bas 68.12 (1949): 1106-1122. + +.. [Yamakawa1997] + Yamakawa, Hiromi, and Takenao Yoshizaki. Helical wormlike chains in polymer solutions. Vol. 1. Berlin: Springer, 1997. + +.. Mathematics references + +.. [Arfken1999] + Arfken, George B., and Hans J. Weber. "Mathematical methods for physicists." (1999): 165-169. diff --git a/notes/chromo-derive.log b/notes/chromo-derive.log index 801f624..63c5952 100644 --- a/notes/chromo-derive.log +++ b/notes/chromo-derive.log @@ -1,4 +1,4 @@ -This is pdfTeX, Version 3.14159265-2.6-1.40.20 (TeX Live 2019) (preloaded format=pdflatex 2019.5.8) 3 AUG 2020 16:42 +This is pdfTeX, Version 3.14159265-2.6-1.40.20 (TeX Live 2019) (preloaded format=pdflatex 2019.5.8) 17 SEP 2020 13:17 entering extended mode restricted \write18 enabled. file:line:error style messages enabled. @@ -357,7 +357,7 @@ File: figures/weight.pdf Graphic file (type pdf) Package pdftex.def Info: figures/weight.pdf used on input line 60. (pdftex.def) Requested size: 281.85587pt x 140.93474pt. LaTeX Font Info: Font shape `OT1/ptm/bx/n' in size <14.4> not available -(Font) Font shape `OT1/ptm/b/n' tried instead on input line 88. +(Font) Font shape `OT1/ptm/b/n' tried instead on input line 90. [2 <./figures/weight.pdf>] (./chromo-derive.bbl) [3] (./chromo-derive.aux) ) @@ -379,7 +379,7 @@ texmf-dist/fonts/type1/urw/symbol/usyr.pfb> -Output written on chromo-derive.pdf (3 pages, 188127 bytes). +Output written on chromo-derive.pdf (3 pages, 189721 bytes). PDF statistics: 64 PDF objects out of 1000 (max. 8388607) 44 compressed objects within 1 object stream diff --git a/notes/chromo-derive.pdf b/notes/chromo-derive.pdf index 3a59cb5..2bf7974 100644 Binary files a/notes/chromo-derive.pdf and b/notes/chromo-derive.pdf differ diff --git a/notes/chromo-derive.synctex.gz b/notes/chromo-derive.synctex.gz index d3714b8..6ac7a59 100644 Binary files a/notes/chromo-derive.synctex.gz and b/notes/chromo-derive.synctex.gz differ diff --git a/source/mc.ipynb b/source/mc.ipynb new file mode 100644 index 0000000..a915c00 --- /dev/null +++ b/source/mc.ipynb @@ -0,0 +1,223 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# How to use MC simulation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The main entry point for MC simulation in the current code base is currently found in `chromo.mc.polymer_in_field`. This file shows the basic usage of that function." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "sys.path.insert(1, '../..')\n", + "\n", + "import importlib\n", + "from io import StringIO\n", + "\n", + "import numpy as np\n", + "import pandas as pd\n", + "\n", + "import chromo\n", + "import chromo.mc as mc\n", + "from chromo.components import Polymer\n", + "import chromo.marks\n", + "from chromo.fields import UniformDensityField" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First specify what epigenetic marks should be used." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "epimarks = [chromo.marks.get_by_name('HP1')]\n", + "#or just e.g. epimarks = [chromo.marks.hp1]\n", + "# making it a list is optional since there's only one element" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[Epigenmark(name='HP1', bind_energy=1, interaction_energy=1, chemical_potential=1)]" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "epimarks" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "They are expected in a DataFrame format." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "marks = chromo.marks.make_mark_collection(epimarks)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now specify the initial state of the polymer(s). We will use just one here." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "p = Polymer.straight_line_in_x('Chr-1', 10, 1, states=np.zeros((10,)), mark_names=['HP1'])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally create a field object to contain the polymers." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "udf = UniformDensityField([p], marks, 10, 20, 30, 40, 50, 60)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As per the docs of `chromo.util.reproducibility.make_reproducible`, check the `simulations.csv` file in the requested output folder (i.e. `{current working directory of this Notebook}/{output_dir}`) to see information about all simulations that have been run so far.\n", + "\n", + "All inputs are saved in the output folder as well, for exact reproducibility." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Save point 0 completed\n", + "Save point 1 completed\n", + "Save point 2 completed\n", + "Save point 3 completed\n", + "Save point 4 completed\n", + "Save point 5 completed\n", + "Save point 6 completed\n", + "Save point 7 completed\n", + "Save point 8 completed\n", + "Save point 9 completed\n" + ] + } + ], + "source": [ + "mc.polymer_in_field([p], marks, udf, 1000, 10, output_dir='output')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here's an alternative implementation that takes only \"simple\" inputs:" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Save point 0 completed\n", + "Save point 1 completed\n", + "Save point 2 completed\n", + "Save point 3 completed\n", + "Save point 4 completed\n", + "Save point 5 completed\n", + "Save point 6 completed\n", + "Save point 7 completed\n", + "Save point 8 completed\n", + "Save point 9 completed\n" + ] + } + ], + "source": [ + "mc.simple_mc(num_polymers=5, num_beads=10, bead_length=1, num_marks=2, num_save_mc=1000, num_saves=10,\n", + " x_width=10, nx=1, y_width=10, ny=1, z_width=10, nz=1, random_seed=0, output_dir='output_simple')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/tests/test_components.py b/tests/test_components.py new file mode 100644 index 0000000..2ce3d56 --- /dev/null +++ b/tests/test_components.py @@ -0,0 +1,14 @@ +from io import StringIO + +from pandas.testing import assert_frame_equal +import numpy as np + +from chromo.components import Polymer + + +def test_polymer_saving(): + p = Polymer.straight_line_in_x('Chr-1', 10, 1, states=np.zeros((10,)), + mark_names=['HP1']) + df = p.from_file(StringIO(p.to_file(None)), 'Chr-1').to_dataframe() + df_pre_round_trip = p.to_dataframe() + assert_frame_equal(df, df_pre_round_trip) diff --git a/tests/test_fields.py b/tests/test_fields.py new file mode 100644 index 0000000..04b39aa --- /dev/null +++ b/tests/test_fields.py @@ -0,0 +1,49 @@ +from io import StringIO + +import numpy as np +import pandas as pd +import pytest + +from chromo.components import Polymer +from chromo.fields import UniformDensityField +import chromo.marks + + +def test_uniform_density_field_roundtrip(): + marks = [chromo.marks.get_by_name('HP1')] + marks = chromo.marks.make_mark_collection(marks) + + p = Polymer.straight_line_in_x('Chr-1', 10, 1, states=np.zeros((10,)), + mark_names=['HP1']) + + udf = UniformDensityField([p], marks, 10, 20, 30, 40, 50, 60) + udf_round_trip = UniformDensityField.from_file( + StringIO(udf.to_file(None)), [p], marks + ) + assert udf == udf_round_trip + + # ensure the code errors if the wrong number of polymers or marks is passed + with pytest.raises(ValueError): + UniformDensityField.from_file( + StringIO(udf.to_file(None)), [], marks + ) + with pytest.raises(ValueError): + UniformDensityField.from_file( + StringIO(udf.to_file(None)), [p, p], marks + ) + with pytest.raises(ValueError): + UniformDensityField.from_file( + StringIO(udf.to_file(None)), [p], pd.concat([marks, marks]) + ) + + # should also error if the poly/mark has the wrong name + with pytest.raises(ValueError): + p.name = 'different' + UniformDensityField.from_file( + StringIO(udf.to_file(None)), [p], marks + ) + with pytest.raises(ValueError): + marks.loc[0, 'name'] = 'different' + UniformDensityField.from_file( + StringIO(udf.to_file(None)), [p], marks + )