From caae62a88758964a506b1ee1d96a3f780a29c210 Mon Sep 17 00:00:00 2001 From: Michael Anderson Date: Thu, 27 Apr 2017 04:45:34 -0700 Subject: [PATCH 01/13] switched to distributed srm --- brainiak/funcalign/srm.py | 181 ++++++++++++++++++++++-------------- tests/funcalign/test_srm.py | 80 +++++++++++----- 2 files changed, 166 insertions(+), 95 deletions(-) diff --git a/brainiak/funcalign/srm.py b/brainiak/funcalign/srm.py index bf4ac0d66..8aa4fb56b 100644 --- a/brainiak/funcalign/srm.py +++ b/brainiak/funcalign/srm.py @@ -38,7 +38,8 @@ import scipy from sklearn.base import BaseEstimator, TransformerMixin from sklearn.utils import assert_all_finite -from sklearn.utils.validation import NotFittedError +from sklearn.exceptions import NotFittedError +from mpi4py import MPI __all__ = [ "SRM", "DetSRM" @@ -47,7 +48,7 @@ logger = logging.getLogger(__name__) -def _init_w_transforms(data, features): +def _init_w_transforms(data, features, comm=MPI.COMM_SELF): """Initialize the mappings (Wi) for the SRM with random orthogonal matrices. Parameters @@ -85,10 +86,16 @@ def _init_w_transforms(data, features): # Set Wi to a random orthogonal voxels by features matrix for subject in range(subjects): - voxels[subject] = data[subject].shape[0] - rnd_matrix = np.random.random((voxels[subject], features)) - q, r = np.linalg.qr(rnd_matrix) - w.append(q) + np.random.seed(subject) + if data[subject] is not None: + voxels[subject] = data[subject].shape[0] + rnd_matrix = np.random.random((voxels[subject], features)) + q, r = np.linalg.qr(rnd_matrix) + w.append(q) + else: + voxels[subject] = 0 + w.append(None) + voxels = comm.allreduce(voxels, op=MPI.SUM) return w, voxels @@ -151,10 +158,11 @@ class SRM(BaseEstimator, TransformerMixin): K - the number of features (typically, :math:`V \\gg T \\gg K`). """ - def __init__(self, n_iter=10, features=50, rand_seed=0): + def __init__(self, n_iter=10, features=50, rand_seed=0, comm=MPI.COMM_SELF): self.n_iter = n_iter self.features = features self.rand_seed = rand_seed + self.comm = comm return def fit(self, X, y=None): @@ -175,19 +183,25 @@ def fit(self, X, y=None): "({0:d}) to train the model.".format(len(X))) # Check for input data sizes - if X[0].shape[1] < self.features: - raise ValueError( - "There are not enough samples to train the model with " - "{0:d} features.".format(self.features)) + number_subjects = len(X) + for subject in range(number_subjects): + if X[subject] is not None: + assert_all_finite(X[subject]) + if X[subject].shape[1] < self.features: + raise ValueError( + "There are not enough samples to train the model with " + "{0:d} features.".format(self.features)) + assert_all_finite(X[subject]) + # Check if all subjects have same number of TRs - number_trs = X[0].shape[1] number_subjects = len(X) + number_trs = self.comm.allreduce(min([x.shape[1] for x in X if x is not None]), op=MPI.MIN) for subject in range(number_subjects): - assert_all_finite(X[subject]) - if X[subject].shape[1] != number_trs: - raise ValueError("Different number of samples between subjects" - ".") + if X[subject] is not None: + if X[subject].shape[1] != number_trs: + raise ValueError("Different number of samples between subjects" + ".") # Run SRM self.sigma_s_, self.w_, self.mu_, self.rho2_, self.s_ = self._srm(X) @@ -222,7 +236,8 @@ def transform(self, X, y=None): s = [None] * len(X) for subject in range(len(X)): - s[subject] = self.w_[subject].T.dot(X[subject]) + if X[subject] is not None: + s[subject] = self.w_[subject].T.dot(X[subject]) return s @@ -258,11 +273,18 @@ def _init_structures(self, data, subjects): rho2 = np.zeros(subjects) trace_xtx = np.zeros(subjects) + for subject in range(subjects): - mu.append(np.mean(data[subject], 1)) + np.random.seed(subject) rho2[subject] = 1 - trace_xtx[subject] = np.sum(data[subject] ** 2) - x.append(data[subject] - mu[subject][:, np.newaxis]) + if data[subject] is not None: + mu.append(np.mean(data[subject], 1)) + trace_xtx[subject] = np.sum(data[subject] ** 2) + x.append(data[subject] - mu[subject][:, np.newaxis]) + else: + mu.append(None) + trace_xtx[subject] = 0 + x.append(None) return x, mu, rho2, trace_xtx @@ -343,7 +365,7 @@ def _srm(self, data): The shared response. """ - samples = data[0].shape[1] + samples = self.comm.allreduce(min([d.shape[1] for d in data if d is not None]), op=MPI.MIN) subjects = len(data) np.random.seed(self.rand_seed) @@ -351,11 +373,13 @@ def _srm(self, data): # Initialization step: initialize the outputs with initial values, # voxels with the number of voxels in each subject, and trace_xtx with # the ||X_i||_F^2 of each subject. - w, voxels = _init_w_transforms(data, self.features) + w, voxels = _init_w_transforms(data, self.features, self.comm) x, mu, rho2, trace_xtx = self._init_structures(data, subjects) shared_response = np.zeros((self.features, samples)) sigma_s = np.identity(self.features) + rank = self.comm.Get_rank() + # Main loop of the algorithm (run for iteration in range(self.n_iter): logger.info('Iteration %d' % (iteration + 1)) @@ -363,68 +387,85 @@ def _srm(self, data): # E-step: # Sum the inverted the rho2 elements for computing W^T * Psi^-1 * W - rho0 = (1 / rho2).sum() - - # Invert Sigma_s using Cholesky factorization - (chol_sigma_s, lower_sigma_s) = scipy.linalg.cho_factor( - sigma_s, check_finite=False) - inv_sigma_s = scipy.linalg.cho_solve( - (chol_sigma_s, lower_sigma_s), np.identity(self.features), - check_finite=False) - - # Invert (Sigma_s + rho_0 * I) using Cholesky factorization - sigma_s_rhos = inv_sigma_s + np.identity(self.features) * rho0 - (chol_sigma_s_rhos, lower_sigma_s_rhos) = scipy.linalg.cho_factor( - sigma_s_rhos, check_finite=False) - inv_sigma_s_rhos = scipy.linalg.cho_solve( - (chol_sigma_s_rhos, lower_sigma_s_rhos), - np.identity(self.features), check_finite=False) + if rank == 0: + rho0 = (1 / rho2).sum() + + # Invert Sigma_s using Cholesky factorization + (chol_sigma_s, lower_sigma_s) = scipy.linalg.cho_factor( + sigma_s, check_finite=False) + inv_sigma_s = scipy.linalg.cho_solve( + (chol_sigma_s, lower_sigma_s), np.identity(self.features), + check_finite=False) + + # Invert (Sigma_s + rho_0 * I) using Cholesky factorization + sigma_s_rhos = inv_sigma_s + np.identity(self.features) * rho0 + (chol_sigma_s_rhos, lower_sigma_s_rhos) = scipy.linalg.cho_factor( + sigma_s_rhos, check_finite=False) + inv_sigma_s_rhos = scipy.linalg.cho_solve( + (chol_sigma_s_rhos, lower_sigma_s_rhos), + np.identity(self.features), check_finite=False) # Compute the sum of W_i^T * rho_i^-2 * X_i, and the sum of traces # of X_i^T * rho_i^-2 * X_i wt_invpsi_x = np.zeros((self.features, samples)) trace_xt_invsigma2_x = 0.0 for subject in range(subjects): - wt_invpsi_x += (w[subject].T.dot(x[subject])) / rho2[subject] - trace_xt_invsigma2_x += trace_xtx[subject] / rho2[subject] + if data[subject] is not None: + wt_invpsi_x += (w[subject].T.dot(x[subject])) / rho2[subject] + trace_xt_invsigma2_x += trace_xtx[subject] / rho2[subject] + + wt_invpsi_x = self.comm.reduce(wt_invpsi_x, op=MPI.SUM) + trace_xt_invsigma2_x = self.comm.reduce(trace_xt_invsigma2_x, op=MPI.SUM) - log_det_psi = np.sum(np.log(rho2) * voxels) + trace_sigma_s = None + if rank == 0: + log_det_psi = np.sum(np.log(rho2) * voxels) - # Update the shared response - shared_response = sigma_s.dot( - np.identity(self.features) - rho0 * inv_sigma_s_rhos).dot( - wt_invpsi_x) + # Update the shared response + shared_response = sigma_s.dot( + np.identity(self.features) - rho0 * inv_sigma_s_rhos).dot( + wt_invpsi_x) - # M-step + # M-step - # Update Sigma_s and compute its trace - sigma_s = (inv_sigma_s_rhos - + shared_response.dot(shared_response.T) / samples) - trace_sigma_s = samples * np.trace(sigma_s) + # Update Sigma_s and compute its trace + sigma_s = (inv_sigma_s_rhos + + shared_response.dot(shared_response.T) / samples) + trace_sigma_s = samples * np.trace(sigma_s) + + shared_response = self.comm.bcast(shared_response) + trace_sigma_s = self.comm.bcast(trace_sigma_s) # Update each subject's mapping transform W_i and error variance # rho_i^2 for subject in range(subjects): - a_subject = x[subject].dot(shared_response.T) - perturbation = np.zeros(a_subject.shape) - np.fill_diagonal(perturbation, 0.001) - u_subject, s_subject, v_subject = np.linalg.svd( - a_subject + perturbation, full_matrices=False) - w[subject] = u_subject.dot(v_subject) - rho2[subject] = trace_xtx[subject] - rho2[subject] += -2 * np.sum(w[subject] * a_subject).sum() - rho2[subject] += trace_sigma_s - rho2[subject] /= samples * voxels[subject] - - if logger.isEnabledFor(logging.INFO): - # Calculate and log the current log-likelihood for checking - # convergence - loglike = self._likelihood( - chol_sigma_s_rhos, log_det_psi, chol_sigma_s, - trace_xt_invsigma2_x, inv_sigma_s_rhos, wt_invpsi_x, - samples) - logger.info('Objective function %f' % loglike) - + if x[subject] is not None: + a_subject = x[subject].dot(shared_response.T) + perturbation = np.zeros(a_subject.shape) + np.fill_diagonal(perturbation, 0.001) + u_subject, s_subject, v_subject = np.linalg.svd( + a_subject + perturbation, full_matrices=False) + w[subject] = u_subject.dot(v_subject) + rho2[subject] = trace_xtx[subject] + rho2[subject] += -2 * np.sum(w[subject] * a_subject).sum() + rho2[subject] += trace_sigma_s + rho2[subject] /= samples * voxels[subject] + else: + rho2[subject] = 0 + + rho2 = self.comm.allreduce(rho2, op=MPI.SUM) + + if rank == 0: + if logger.isEnabledFor(logging.INFO): + # Calculate and log the current log-likelihood for checking + # convergence + loglike = self._likelihood( + chol_sigma_s_rhos, log_det_psi, chol_sigma_s, + trace_xt_invsigma2_x, inv_sigma_s_rhos, wt_invpsi_x, + samples) + logger.info('Objective function %f' % loglike) + + sigma_s = self.comm.bcast(sigma_s) return sigma_s, w, mu, rho2, shared_response diff --git a/tests/funcalign/test_srm.py b/tests/funcalign/test_srm.py index 5f9694260..e2b0391be 100644 --- a/tests/funcalign/test_srm.py +++ b/tests/funcalign/test_srm.py @@ -11,8 +11,9 @@ # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. -from sklearn.utils.validation import NotFittedError +from sklearn.exceptions import NotFittedError import pytest +from mpi4py import MPI def test_can_instantiate(): import brainiak.funcalign.srm @@ -20,13 +21,16 @@ def test_can_instantiate(): assert s, "Invalid SRM instance!" import numpy as np + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + nrank = comm.Get_size() voxels = 100 samples = 500 subjects = 2 features = 3 - s = brainiak.funcalign.srm.SRM(n_iter=5, features=features) + s = brainiak.funcalign.srm.SRM(n_iter=5, features=features, comm=comm) assert s, "Invalid SRM instance!" @@ -39,39 +43,58 @@ def test_can_instantiate(): S = np.vstack((x, y, z)) + # DSRM: broadcast S + S = comm.bcast(S) + X = [] W = [] - Q, R = np.linalg.qr(np.random.random((voxels, features))) - W.append(Q) - X.append(Q.dot(S) + 0.1*np.random.random((voxels, samples))) + # DSRM: only append on rank 0 + if rank == 0: + np.random.seed(0) + Q, R = np.linalg.qr(np.random.random((voxels, features))) + W.append(Q) + X.append(Q.dot(S) + 0.1*np.random.random((voxels, samples))) + else: + W.append(None) + X.append(None) # Check that transform does NOT run before fitting the model with pytest.raises(NotFittedError) as excinfo: s.transform(X) - print("Test: transforming before fitting the model") + if rank == 0: + print("Test: transforming before fitting the model") # Check that it does NOT run with 1 subject with pytest.raises(ValueError) as excinfo: s.fit(X) - print("Test: running SRM with 1 subject") + if rank == 0: + print("Test: running SRM with 1 subject") + # DSRM: cyclic distribution of subject data, otherwise None for subject in range(1, subjects): - Q, R = np.linalg.qr(np.random.random((voxels, features))) - W.append(Q) - X.append(Q.dot(S) + 0.1*np.random.random((voxels, samples))) + np.random.seed(subject) + if subject % nrank == rank: + Q, R = np.linalg.qr(np.random.random((voxels, features))) + W.append(Q) + X.append(Q.dot(S) + 0.1*np.random.random((voxels, samples))) + else: + W.append(None) + X.append(None) # Check that runs with 2 subject s.fit(X) assert len(s.w_) == subjects, "Invalid computation of SRM! (wrong # subjects in W)" for subject in range(subjects): - assert s.w_[subject].shape[0] == voxels, "Invalid computation of SRM! (wrong # voxels in W)" - assert s.w_[subject].shape[1] == features, "Invalid computation of SRM! (wrong # features in W)" - ortho = np.linalg.norm(s.w_[subject].T.dot(s.w_[subject]) - np.eye(s.w_[subject].shape[1]), 'fro') - assert ortho < 1e-7, "A Wi mapping is not orthonormal in SRM." - difference = np.linalg.norm(X[subject] - s.w_[subject].dot(s.s_), 'fro') - datanorm = np.linalg.norm(X[subject], 'fro') - assert difference/datanorm < 1.0, "Model seems incorrectly computed." + if s.w_[subject] is not None: + assert s.w_[subject].shape[0] == voxels, "Invalid computation of SRM! (wrong # voxels in W)" + assert s.w_[subject].shape[1] == features, "Invalid computation of SRM! (wrong # features in W)" + ortho = np.linalg.norm(s.w_[subject].T.dot(s.w_[subject]) - np.eye(s.w_[subject].shape[1]), 'fro') + assert ortho < 1e-7, "A Wi mapping is not orthonormal in SRM." + difference = np.linalg.norm(X[subject] - s.w_[subject].dot(s.s_), 'fro') + datanorm = np.linalg.norm(X[subject], 'fro') + print('Difference/datanorm: ' + str(difference/datanorm)) + assert difference/datanorm < 1.0, "Model seems incorrectly computed." assert s.s_.shape[0] == features, "Invalid computation of SRM! (wrong # features in S)" assert s.s_.shape[1] == samples, "Invalid computation of SRM! (wrong # samples in S)" @@ -80,26 +103,32 @@ def test_can_instantiate(): assert len(new_s) == subjects, "Invalid computation of SRM! (wrong # subjects after transform)" for subject in range(subjects): - assert new_s[subject].shape[0] == features, "Invalid computation of SRM! (wrong # features after transform)" - assert new_s[subject].shape[1] == samples, "Invalid computation of SRM! (wrong # samples after transform)" + if new_s[subject] is not None: + assert new_s[subject].shape[0] == features, "Invalid computation of SRM! (wrong # features after transform)" + assert new_s[subject].shape[1] == samples, "Invalid computation of SRM! (wrong # samples after transform)" # Check that it does NOT run with non-matching number of subjects with pytest.raises(ValueError) as excinfo: - s.transform(X[1]) - print("Test: transforming with non-matching number of subjects") + s.transform([X[1]]) + if rank == 0: + print("Test: transforming with non-matching number of subjects") # Check that it does not run without enough samples (TRs). with pytest.raises(ValueError) as excinfo: s.set_params(features=(samples+1)) s.fit(X) - print("Test: not enough samples") + if rank == 0: + print("Test: not enough samples") # Check that it does not run with different number of samples (TRs) - S2 = S[:,:-2] - X.append(Q.dot(S2)) + if rank == 0: + S2 = S[:,:-2] + X.append(Q.dot(S2)) + with pytest.raises(ValueError) as excinfo: s.fit(X) - print("Test: different number of samples per subject") + if rank == 0: + print("Test: different number of samples per subject") @@ -191,3 +220,4 @@ def test_det_srm(): print("Test: different number of samples per subject") +test_can_instantiate() From 65fe066994e2347733c8fa516a33ff44f3ff5612 Mon Sep 17 00:00:00 2001 From: Michael Anderson Date: Thu, 27 Apr 2017 07:29:44 -0700 Subject: [PATCH 02/13] passes test for -n 3 --- brainiak/funcalign/srm.py | 35 +++++++++++++++++++++++------------ tests/funcalign/test_srm.py | 10 +++++++++- 2 files changed, 32 insertions(+), 13 deletions(-) diff --git a/brainiak/funcalign/srm.py b/brainiak/funcalign/srm.py index 8aa4fb56b..95674bb3f 100644 --- a/brainiak/funcalign/srm.py +++ b/brainiak/funcalign/srm.py @@ -40,6 +40,7 @@ from sklearn.utils import assert_all_finite from sklearn.exceptions import NotFittedError from mpi4py import MPI +import sys __all__ = [ "SRM", "DetSRM" @@ -184,25 +185,35 @@ def fit(self, X, y=None): # Check for input data sizes number_subjects = len(X) + number_subjects_vec = self.comm.allgather(number_subjects) + for rank in range(self.comm.Get_size()): + if number_subjects_vec[rank] != number_subjects: + raise ValueError( + "Not all ranks have same number of subjects") + + # Collect size information + shape0 = np.zeros((number_subjects,), dtype=np.int) + shape1 = np.zeros((number_subjects,), dtype=np.int) + for subject in range(number_subjects): if X[subject] is not None: assert_all_finite(X[subject]) - if X[subject].shape[1] < self.features: - raise ValueError( - "There are not enough samples to train the model with " - "{0:d} features.".format(self.features)) - assert_all_finite(X[subject]) + shape0[subject] = X[subject].shape[0] + shape1[subject] = X[subject].shape[1] + shape0 = self.comm.allreduce(shape0, op=MPI.SUM) + shape1 = self.comm.allreduce(shape1, op=MPI.SUM) # Check if all subjects have same number of TRs - number_subjects = len(X) - number_trs = self.comm.allreduce(min([x.shape[1] for x in X if x is not None]), op=MPI.MIN) + number_trs = np.min(shape1) for subject in range(number_subjects): - if X[subject] is not None: - if X[subject].shape[1] != number_trs: - raise ValueError("Different number of samples between subjects" + if shape1[subject] < self.features: + raise ValueError( + "There are not enough samples to train the model with " + "{0:d} features.".format(self.features)) + if shape1[subject] != number_trs: + raise ValueError("Different number of samples between subjects" ".") - # Run SRM self.sigma_s_, self.w_, self.mu_, self.rho2_, self.s_ = self._srm(X) @@ -365,7 +376,7 @@ def _srm(self, data): The shared response. """ - samples = self.comm.allreduce(min([d.shape[1] for d in data if d is not None]), op=MPI.MIN) + samples = self.comm.allreduce(min([d.shape[1] for d in data if d is not None],default=sys.maxsize), op=MPI.MIN) subjects = len(data) np.random.seed(self.rand_seed) diff --git a/tests/funcalign/test_srm.py b/tests/funcalign/test_srm.py index e2b0391be..2aeaa08f5 100644 --- a/tests/funcalign/test_srm.py +++ b/tests/funcalign/test_srm.py @@ -33,7 +33,6 @@ def test_can_instantiate(): s = brainiak.funcalign.srm.SRM(n_iter=5, features=features, comm=comm) assert s, "Invalid SRM instance!" - # Create a Shared response S with K = 3 theta = np.linspace(-4 * np.pi, 4 * np.pi, samples) z = np.linspace(-2, 2, samples) @@ -95,6 +94,7 @@ def test_can_instantiate(): datanorm = np.linalg.norm(X[subject], 'fro') print('Difference/datanorm: ' + str(difference/datanorm)) assert difference/datanorm < 1.0, "Model seems incorrectly computed." + assert s.s_.shape[0] == features, "Invalid computation of SRM! (wrong # features in S)" assert s.s_.shape[1] == samples, "Invalid computation of SRM! (wrong # samples in S)" @@ -120,16 +120,24 @@ def test_can_instantiate(): if rank == 0: print("Test: not enough samples") + # Check that it does not run with different number of samples (TRs) if rank == 0: S2 = S[:,:-2] X.append(Q.dot(S2)) + else: + X.append(None) + with pytest.raises(ValueError) as excinfo: s.fit(X) if rank == 0: print("Test: different number of samples per subject") + import sys + sys.exit(0) + + def test_det_srm(): From 8ff2ef7d7d5107c0adbc245bd7c521dbb0e03f0a Mon Sep 17 00:00:00 2001 From: Michael Anderson Date: Thu, 27 Apr 2017 11:02:39 -0700 Subject: [PATCH 03/13] passes pr-check --- brainiak/funcalign/srm.py | 26 +++++++++++++++----------- tests/funcalign/test_srm.py | 5 ----- 2 files changed, 15 insertions(+), 16 deletions(-) diff --git a/brainiak/funcalign/srm.py b/brainiak/funcalign/srm.py index 95674bb3f..b24d1a0d2 100644 --- a/brainiak/funcalign/srm.py +++ b/brainiak/funcalign/srm.py @@ -159,7 +159,8 @@ class SRM(BaseEstimator, TransformerMixin): K - the number of features (typically, :math:`V \\gg T \\gg K`). """ - def __init__(self, n_iter=10, features=50, rand_seed=0, comm=MPI.COMM_SELF): + def __init__(self, n_iter=10, features=50, rand_seed=0, + comm=MPI.COMM_SELF): self.n_iter = n_iter self.features = features self.rand_seed = rand_seed @@ -213,7 +214,7 @@ def fit(self, X, y=None): "{0:d} features.".format(self.features)) if shape1[subject] != number_trs: raise ValueError("Different number of samples between subjects" - ".") + ".") # Run SRM self.sigma_s_, self.w_, self.mu_, self.rho2_, self.s_ = self._srm(X) @@ -376,9 +377,10 @@ def _srm(self, data): The shared response. """ - samples = self.comm.allreduce(min([d.shape[1] for d in data if d is not None],default=sys.maxsize), op=MPI.MIN) + local_min = min([d.shape[1] for d in data if d is not None], + default=sys.maxsize) + samples = self.comm.allreduce(local_min, op=MPI.MIN) subjects = len(data) - np.random.seed(self.rand_seed) # Initialization step: initialize the outputs with initial values, @@ -407,11 +409,12 @@ def _srm(self, data): inv_sigma_s = scipy.linalg.cho_solve( (chol_sigma_s, lower_sigma_s), np.identity(self.features), check_finite=False) - + # Invert (Sigma_s + rho_0 * I) using Cholesky factorization sigma_s_rhos = inv_sigma_s + np.identity(self.features) * rho0 - (chol_sigma_s_rhos, lower_sigma_s_rhos) = scipy.linalg.cho_factor( - sigma_s_rhos, check_finite=False) + chol_sigma_s_rhos, lower_sigma_s_rhos = \ + scipy.linalg.cho_factor(sigma_s_rhos, + check_finite=False) inv_sigma_s_rhos = scipy.linalg.cho_solve( (chol_sigma_s_rhos, lower_sigma_s_rhos), np.identity(self.features), check_finite=False) @@ -422,12 +425,13 @@ def _srm(self, data): trace_xt_invsigma2_x = 0.0 for subject in range(subjects): if data[subject] is not None: - wt_invpsi_x += (w[subject].T.dot(x[subject])) / rho2[subject] + wt_invpsi_x += (w[subject].T.dot(x[subject])) \ + / rho2[subject] trace_xt_invsigma2_x += trace_xtx[subject] / rho2[subject] wt_invpsi_x = self.comm.reduce(wt_invpsi_x, op=MPI.SUM) - trace_xt_invsigma2_x = self.comm.reduce(trace_xt_invsigma2_x, op=MPI.SUM) - + trace_xt_invsigma2_x = self.comm.reduce(trace_xt_invsigma2_x, + op=MPI.SUM) trace_sigma_s = None if rank == 0: log_det_psi = np.sum(np.log(rho2) * voxels) @@ -475,7 +479,7 @@ def _srm(self, data): trace_xt_invsigma2_x, inv_sigma_s_rhos, wt_invpsi_x, samples) logger.info('Objective function %f' % loglike) - + sigma_s = self.comm.bcast(sigma_s) return sigma_s, w, mu, rho2, shared_response diff --git a/tests/funcalign/test_srm.py b/tests/funcalign/test_srm.py index 2aeaa08f5..2a1333b78 100644 --- a/tests/funcalign/test_srm.py +++ b/tests/funcalign/test_srm.py @@ -134,11 +134,6 @@ def test_can_instantiate(): if rank == 0: print("Test: different number of samples per subject") - import sys - sys.exit(0) - - - def test_det_srm(): import brainiak.funcalign.srm From f8d05bd1db170dc6c754412a7ce3ba2db1f9b23f Mon Sep 17 00:00:00 2001 From: Michael Anderson Date: Fri, 28 Apr 2017 04:50:21 -0700 Subject: [PATCH 04/13] Added distributed example --- brainiak/funcalign/srm.py | 36 +++-- ...rm_image_prediction_example_distributed.py | 147 ++++++++++++++++++ 2 files changed, 168 insertions(+), 15 deletions(-) create mode 100644 examples/funcalign/srm_image_prediction_example_distributed.py diff --git a/brainiak/funcalign/srm.py b/brainiak/funcalign/srm.py index b24d1a0d2..ccb81f1d9 100644 --- a/brainiak/funcalign/srm.py +++ b/brainiak/funcalign/srm.py @@ -61,6 +61,8 @@ def _init_w_transforms(data, features, comm=MPI.COMM_SELF): features : int The number of features in the model. + comm : mpi4py.MPI.Intracomm + The MPI communicator containing the data Returns ------- @@ -121,6 +123,8 @@ class SRM(BaseEstimator, TransformerMixin): rand_seed : int, default: 0 Seed for initializing the random number generator. + comm : mpi4py.MPI.Intracomm + The MPI communicator containing the data Attributes ---------- @@ -140,6 +144,8 @@ class SRM(BaseEstimator, TransformerMixin): rho2_ : array, shape=[subjects] The estimated noise variance :math:`\\rho_i^2` for each subject + comm_ : mpi4py.MPI.Intracomm + The MPI communicator containing the data Note ---- @@ -164,7 +170,7 @@ def __init__(self, n_iter=10, features=50, rand_seed=0, self.n_iter = n_iter self.features = features self.rand_seed = rand_seed - self.comm = comm + self.comm_ = comm return def fit(self, X, y=None): @@ -186,8 +192,8 @@ def fit(self, X, y=None): # Check for input data sizes number_subjects = len(X) - number_subjects_vec = self.comm.allgather(number_subjects) - for rank in range(self.comm.Get_size()): + number_subjects_vec = self.comm_.allgather(number_subjects) + for rank in range(self.comm_.Get_size()): if number_subjects_vec[rank] != number_subjects: raise ValueError( "Not all ranks have same number of subjects") @@ -202,8 +208,8 @@ def fit(self, X, y=None): shape0[subject] = X[subject].shape[0] shape1[subject] = X[subject].shape[1] - shape0 = self.comm.allreduce(shape0, op=MPI.SUM) - shape1 = self.comm.allreduce(shape1, op=MPI.SUM) + shape0 = self.comm_.allreduce(shape0, op=MPI.SUM) + shape1 = self.comm_.allreduce(shape1, op=MPI.SUM) # Check if all subjects have same number of TRs number_trs = np.min(shape1) @@ -379,19 +385,19 @@ def _srm(self, data): local_min = min([d.shape[1] for d in data if d is not None], default=sys.maxsize) - samples = self.comm.allreduce(local_min, op=MPI.MIN) + samples = self.comm_.allreduce(local_min, op=MPI.MIN) subjects = len(data) np.random.seed(self.rand_seed) # Initialization step: initialize the outputs with initial values, # voxels with the number of voxels in each subject, and trace_xtx with # the ||X_i||_F^2 of each subject. - w, voxels = _init_w_transforms(data, self.features, self.comm) + w, voxels = _init_w_transforms(data, self.features, self.comm_) x, mu, rho2, trace_xtx = self._init_structures(data, subjects) shared_response = np.zeros((self.features, samples)) sigma_s = np.identity(self.features) - rank = self.comm.Get_rank() + rank = self.comm_.Get_rank() # Main loop of the algorithm (run for iteration in range(self.n_iter): @@ -429,9 +435,9 @@ def _srm(self, data): / rho2[subject] trace_xt_invsigma2_x += trace_xtx[subject] / rho2[subject] - wt_invpsi_x = self.comm.reduce(wt_invpsi_x, op=MPI.SUM) - trace_xt_invsigma2_x = self.comm.reduce(trace_xt_invsigma2_x, - op=MPI.SUM) + wt_invpsi_x = self.comm_.reduce(wt_invpsi_x, op=MPI.SUM) + trace_xt_invsigma2_x = self.comm_.reduce(trace_xt_invsigma2_x, + op=MPI.SUM) trace_sigma_s = None if rank == 0: log_det_psi = np.sum(np.log(rho2) * voxels) @@ -448,8 +454,8 @@ def _srm(self, data): + shared_response.dot(shared_response.T) / samples) trace_sigma_s = samples * np.trace(sigma_s) - shared_response = self.comm.bcast(shared_response) - trace_sigma_s = self.comm.bcast(trace_sigma_s) + shared_response = self.comm_.bcast(shared_response) + trace_sigma_s = self.comm_.bcast(trace_sigma_s) # Update each subject's mapping transform W_i and error variance # rho_i^2 @@ -468,7 +474,7 @@ def _srm(self, data): else: rho2[subject] = 0 - rho2 = self.comm.allreduce(rho2, op=MPI.SUM) + rho2 = self.comm_.allreduce(rho2, op=MPI.SUM) if rank == 0: if logger.isEnabledFor(logging.INFO): @@ -480,7 +486,7 @@ def _srm(self, data): samples) logger.info('Objective function %f' % loglike) - sigma_s = self.comm.bcast(sigma_s) + sigma_s = self.comm_.bcast(sigma_s) return sigma_s, w, mu, rho2, shared_response diff --git a/examples/funcalign/srm_image_prediction_example_distributed.py b/examples/funcalign/srm_image_prediction_example_distributed.py new file mode 100644 index 000000000..2c164ba4f --- /dev/null +++ b/examples/funcalign/srm_image_prediction_example_distributed.py @@ -0,0 +1,147 @@ +# Copyright 2016 Intel Corporation +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +import scipy.io +from scipy.stats import stats +from sklearn.metrics import confusion_matrix +from sklearn.svm import NuSVC +import numpy as np +import brainiak.funcalign.srm +from mpi4py import MPI + +comm = MPI.COMM_WORLD +rank = comm.Get_rank() +size = comm.Get_size() + +# Load the input data that contains the movie stimuli for unsupervised training with SRM +movie_file = scipy.io.loadmat('data/movie_data.mat') +movie_data_left = movie_file['movie_data_lh'] +movie_data_right = movie_file['movie_data_rh'] +subjects = movie_data_left.shape[2] + +# Convert data to a list of arrays matching SRM input. +# Each element is a matrix of voxels by TRs. +# Also, concatenate data from both hemispheres in the brain. +movie_data = [] +for s in range(subjects): + if s % size == rank: + movie_data.append(np.concatenate([movie_data_left[:, :, s], movie_data_right[:, :, s]], axis=0)) + else: + movie_data.append(None) + +# Z-score the data +for subject in range(subjects): + if movie_data[subject] is not None: + movie_data[subject] = stats.zscore(movie_data[subject], axis=1, ddof=1) + +# Run SRM with the movie data +srm = brainiak.funcalign.srm.SRM(n_iter=10, features=50, comm=comm) +srm.fit(movie_data) + +# We define a function to present the output of the experiment. +def plot_confusion_matrix(cm, title="Confusion Matrix"): + """Plots a confusion matrix for each subject + """ + import matplotlib.pyplot as plt + import math + plt.figure() + subjects = len(cm) + root_subjects = math.sqrt(subjects) + cols = math.ceil(root_subjects) + rows = math.ceil(subjects/cols) + classes = cm[0].shape[0] + for subject in range(subjects): + plt.subplot(rows, cols, subject+1) + plt.imshow(cm[subject], interpolation='nearest', cmap=plt.cm.bone) + plt.xticks(np.arange(classes), range(1, classes+1)) + plt.yticks(np.arange(classes), range(1, classes+1)) + cbar = plt.colorbar(ticks=[0.0, 1.0], shrink=0.6) + cbar.set_clim(0.0, 1.0) + plt.xlabel("Predicted") + plt.ylabel("True label") + plt.title("{0:d}".format(subject + 1)) + plt.suptitle(title) + plt.tight_layout() + plt.show() + +# Load the input data that contains the image stimuli and its labels for training a classifier +image_file = scipy.io.loadmat('data/image_data.mat') +image_data_left = image_file['image_data_lh'] +image_data_right = image_file['image_data_rh'] + +# Convert data to a list of arrays matching SRM input. +# Each element is a matrix of voxels by TRs. +# Also, concatenate data from both hemispheres in the brain. +image_data = [] +for s in range(subjects): + if s % size == rank: + image_data.append(np.concatenate([image_data_left[:, :, s], image_data_right[:, :, s]], axis=0)) + else: + image_data.append(None) + +# Z-score the image data +for subject in range(subjects): + if subject % size == rank: + image_data[subject] = stats.zscore(image_data[subject], axis=1, ddof=1) + else: + image_data[subject] = None + +# Z-score the shared response data +image_data_shared = srm.transform(image_data) + +for subject in range(subjects): + image_data_shared[subject] = comm.bcast(image_data_shared[subject], + root=(subject % size)) + +if rank == 0: + for subject in range(subjects): + image_data_shared[subject] = stats.zscore(image_data_shared[subject], axis=1, ddof=1) + + # Read the labels of the image data for training the classifier. + labels = scipy.io.loadmat('data/label.mat') + labels = np.squeeze(labels['label']) + + # Run a leave-one-out cross validation with the subjects + train_labels = np.tile(labels, subjects-1) + test_labels = labels + accuracy = np.zeros((subjects,)) + cm = [None] * subjects + for subject in range(subjects): + # Concatenate the subjects' data for training into one matrix + train_subjects = list(range(subjects)) + train_subjects.remove(subject) + TRs = image_data_shared[0].shape[1] + train_data = np.zeros((image_data_shared[0].shape[0], len(train_labels))) + for train_subject in range(len(train_subjects)): + start_index = train_subject*TRs + end_index = start_index+TRs + train_data[:, start_index:end_index] = image_data_shared[train_subjects[train_subject]] + + # Train a Nu-SVM classifier using scikit learn + classifier = NuSVC(nu=0.5, kernel='linear') + classifier = classifier.fit(train_data.T, train_labels) + + # Predict on the test data + predicted_labels = classifier.predict(image_data_shared[subject].T) + accuracy[subject] = sum(predicted_labels == test_labels)/float(len(predicted_labels)) + + # Create a confusion matrix to see the accuracy of each class + cm[subject] = confusion_matrix(test_labels, predicted_labels) + + # Normalize the confusion matrix + cm[subject] = cm[subject].astype('float') / cm[subject].sum(axis=1)[:, np.newaxis] + + + # Plot and print the results + plot_confusion_matrix(cm, title="Confusion matrices for different test subjects with Probabilistic SRM") + print("SRM: The average accuracy among all subjects is {0:f} +/- {1:f}".format(np.mean(accuracy), np.std(accuracy))) From cf8d818e9b12d3f9f782ae4bb86601b787610330 Mon Sep 17 00:00:00 2001 From: Michael Anderson Date: Fri, 28 Apr 2017 07:22:01 -0700 Subject: [PATCH 05/13] pr changes --- brainiak/funcalign/srm.py | 31 +++++++++++++++---------------- tests/funcalign/test_srm.py | 1 - 2 files changed, 15 insertions(+), 17 deletions(-) diff --git a/brainiak/funcalign/srm.py b/brainiak/funcalign/srm.py index ccb81f1d9..cd7d654e9 100644 --- a/brainiak/funcalign/srm.py +++ b/brainiak/funcalign/srm.py @@ -144,7 +144,7 @@ class SRM(BaseEstimator, TransformerMixin): rho2_ : array, shape=[subjects] The estimated noise variance :math:`\\rho_i^2` for each subject - comm_ : mpi4py.MPI.Intracomm + comm : mpi4py.MPI.Intracomm The MPI communicator containing the data Note @@ -170,7 +170,7 @@ def __init__(self, n_iter=10, features=50, rand_seed=0, self.n_iter = n_iter self.features = features self.rand_seed = rand_seed - self.comm_ = comm + self.comm = comm return def fit(self, X, y=None): @@ -192,8 +192,8 @@ def fit(self, X, y=None): # Check for input data sizes number_subjects = len(X) - number_subjects_vec = self.comm_.allgather(number_subjects) - for rank in range(self.comm_.Get_size()): + number_subjects_vec = self.comm.allgather(number_subjects) + for rank in range(self.comm.Get_size()): if number_subjects_vec[rank] != number_subjects: raise ValueError( "Not all ranks have same number of subjects") @@ -208,8 +208,8 @@ def fit(self, X, y=None): shape0[subject] = X[subject].shape[0] shape1[subject] = X[subject].shape[1] - shape0 = self.comm_.allreduce(shape0, op=MPI.SUM) - shape1 = self.comm_.allreduce(shape1, op=MPI.SUM) + shape0 = self.comm.allreduce(shape0, op=MPI.SUM) + shape1 = self.comm.allreduce(shape1, op=MPI.SUM) # Check if all subjects have same number of TRs number_trs = np.min(shape1) @@ -293,7 +293,6 @@ def _init_structures(self, data, subjects): trace_xtx = np.zeros(subjects) for subject in range(subjects): - np.random.seed(subject) rho2[subject] = 1 if data[subject] is not None: mu.append(np.mean(data[subject], 1)) @@ -385,19 +384,19 @@ def _srm(self, data): local_min = min([d.shape[1] for d in data if d is not None], default=sys.maxsize) - samples = self.comm_.allreduce(local_min, op=MPI.MIN) + samples = self.comm.allreduce(local_min, op=MPI.MIN) subjects = len(data) np.random.seed(self.rand_seed) # Initialization step: initialize the outputs with initial values, # voxels with the number of voxels in each subject, and trace_xtx with # the ||X_i||_F^2 of each subject. - w, voxels = _init_w_transforms(data, self.features, self.comm_) + w, voxels = _init_w_transforms(data, self.features, self.comm) x, mu, rho2, trace_xtx = self._init_structures(data, subjects) shared_response = np.zeros((self.features, samples)) sigma_s = np.identity(self.features) - rank = self.comm_.Get_rank() + rank = self.comm.Get_rank() # Main loop of the algorithm (run for iteration in range(self.n_iter): @@ -435,8 +434,8 @@ def _srm(self, data): / rho2[subject] trace_xt_invsigma2_x += trace_xtx[subject] / rho2[subject] - wt_invpsi_x = self.comm_.reduce(wt_invpsi_x, op=MPI.SUM) - trace_xt_invsigma2_x = self.comm_.reduce(trace_xt_invsigma2_x, + wt_invpsi_x = self.comm.reduce(wt_invpsi_x, op=MPI.SUM) + trace_xt_invsigma2_x = self.comm.reduce(trace_xt_invsigma2_x, op=MPI.SUM) trace_sigma_s = None if rank == 0: @@ -454,8 +453,8 @@ def _srm(self, data): + shared_response.dot(shared_response.T) / samples) trace_sigma_s = samples * np.trace(sigma_s) - shared_response = self.comm_.bcast(shared_response) - trace_sigma_s = self.comm_.bcast(trace_sigma_s) + shared_response = self.comm.bcast(shared_response) + trace_sigma_s = self.comm.bcast(trace_sigma_s) # Update each subject's mapping transform W_i and error variance # rho_i^2 @@ -474,7 +473,7 @@ def _srm(self, data): else: rho2[subject] = 0 - rho2 = self.comm_.allreduce(rho2, op=MPI.SUM) + rho2 = self.comm.allreduce(rho2, op=MPI.SUM) if rank == 0: if logger.isEnabledFor(logging.INFO): @@ -486,7 +485,7 @@ def _srm(self, data): samples) logger.info('Objective function %f' % loglike) - sigma_s = self.comm_.bcast(sigma_s) + sigma_s = self.comm.bcast(sigma_s) return sigma_s, w, mu, rho2, shared_response diff --git a/tests/funcalign/test_srm.py b/tests/funcalign/test_srm.py index 2a1333b78..b07cbe000 100644 --- a/tests/funcalign/test_srm.py +++ b/tests/funcalign/test_srm.py @@ -223,4 +223,3 @@ def test_det_srm(): print("Test: different number of samples per subject") -test_can_instantiate() From cb4e6dd68ac32e4d2a4b1ca4d31b65606eb0aa91 Mon Sep 17 00:00:00 2001 From: Michael Anderson Date: Mon, 1 May 2017 04:17:15 -0700 Subject: [PATCH 06/13] one random seed per rank --- brainiak/funcalign/srm.py | 1 - 1 file changed, 1 deletion(-) diff --git a/brainiak/funcalign/srm.py b/brainiak/funcalign/srm.py index cd7d654e9..9d832ccb9 100644 --- a/brainiak/funcalign/srm.py +++ b/brainiak/funcalign/srm.py @@ -89,7 +89,6 @@ def _init_w_transforms(data, features, comm=MPI.COMM_SELF): # Set Wi to a random orthogonal voxels by features matrix for subject in range(subjects): - np.random.seed(subject) if data[subject] is not None: voxels[subject] = data[subject].shape[0] rnd_matrix = np.random.random((voxels[subject], features)) From 6a950a36ec5727821cf77366062cd325dde483d8 Mon Sep 17 00:00:00 2001 From: Michael Anderson Date: Mon, 1 May 2017 04:50:48 -0700 Subject: [PATCH 07/13] fixed over indent --- brainiak/funcalign/srm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/brainiak/funcalign/srm.py b/brainiak/funcalign/srm.py index 9d832ccb9..5d44e195d 100644 --- a/brainiak/funcalign/srm.py +++ b/brainiak/funcalign/srm.py @@ -435,7 +435,7 @@ def _srm(self, data): wt_invpsi_x = self.comm.reduce(wt_invpsi_x, op=MPI.SUM) trace_xt_invsigma2_x = self.comm.reduce(trace_xt_invsigma2_x, - op=MPI.SUM) + op=MPI.SUM) trace_sigma_s = None if rank == 0: log_det_psi = np.sum(np.log(rho2) * voxels) From 2fea955c921ff6d8364574555536321d1f975967 Mon Sep 17 00:00:00 2001 From: Michael Anderson Date: Wed, 3 May 2017 06:39:52 -0700 Subject: [PATCH 08/13] Added old srm test --- tests/funcalign/test_srm.py | 84 +++++---------- tests/funcalign/test_srm_distributed.py | 136 ++++++++++++++++++++++++ 2 files changed, 162 insertions(+), 58 deletions(-) create mode 100644 tests/funcalign/test_srm_distributed.py diff --git a/tests/funcalign/test_srm.py b/tests/funcalign/test_srm.py index b07cbe000..caeae8cea 100644 --- a/tests/funcalign/test_srm.py +++ b/tests/funcalign/test_srm.py @@ -13,7 +13,6 @@ # limitations under the License. from sklearn.exceptions import NotFittedError import pytest -from mpi4py import MPI def test_can_instantiate(): import brainiak.funcalign.srm @@ -21,18 +20,16 @@ def test_can_instantiate(): assert s, "Invalid SRM instance!" import numpy as np - comm = MPI.COMM_WORLD - rank = comm.Get_rank() - nrank = comm.Get_size() voxels = 100 samples = 500 subjects = 2 features = 3 - s = brainiak.funcalign.srm.SRM(n_iter=5, features=features, comm=comm) + s = brainiak.funcalign.srm.SRM(n_iter=5, features=features) assert s, "Invalid SRM instance!" + # Create a Shared response S with K = 3 theta = np.linspace(-4 * np.pi, 4 * np.pi, samples) z = np.linspace(-2, 2, samples) @@ -42,59 +39,39 @@ def test_can_instantiate(): S = np.vstack((x, y, z)) - # DSRM: broadcast S - S = comm.bcast(S) - X = [] W = [] - # DSRM: only append on rank 0 - if rank == 0: - np.random.seed(0) - Q, R = np.linalg.qr(np.random.random((voxels, features))) - W.append(Q) - X.append(Q.dot(S) + 0.1*np.random.random((voxels, samples))) - else: - W.append(None) - X.append(None) + Q, R = np.linalg.qr(np.random.random((voxels, features))) + W.append(Q) + X.append(Q.dot(S) + 0.1*np.random.random((voxels, samples))) # Check that transform does NOT run before fitting the model with pytest.raises(NotFittedError) as excinfo: s.transform(X) - if rank == 0: - print("Test: transforming before fitting the model") + print("Test: transforming before fitting the model") # Check that it does NOT run with 1 subject with pytest.raises(ValueError) as excinfo: s.fit(X) - if rank == 0: - print("Test: running SRM with 1 subject") + print("Test: running SRM with 1 subject") - # DSRM: cyclic distribution of subject data, otherwise None for subject in range(1, subjects): - np.random.seed(subject) - if subject % nrank == rank: - Q, R = np.linalg.qr(np.random.random((voxels, features))) - W.append(Q) - X.append(Q.dot(S) + 0.1*np.random.random((voxels, samples))) - else: - W.append(None) - X.append(None) + Q, R = np.linalg.qr(np.random.random((voxels, features))) + W.append(Q) + X.append(Q.dot(S) + 0.1*np.random.random((voxels, samples))) # Check that runs with 2 subject s.fit(X) assert len(s.w_) == subjects, "Invalid computation of SRM! (wrong # subjects in W)" for subject in range(subjects): - if s.w_[subject] is not None: - assert s.w_[subject].shape[0] == voxels, "Invalid computation of SRM! (wrong # voxels in W)" - assert s.w_[subject].shape[1] == features, "Invalid computation of SRM! (wrong # features in W)" - ortho = np.linalg.norm(s.w_[subject].T.dot(s.w_[subject]) - np.eye(s.w_[subject].shape[1]), 'fro') - assert ortho < 1e-7, "A Wi mapping is not orthonormal in SRM." - difference = np.linalg.norm(X[subject] - s.w_[subject].dot(s.s_), 'fro') - datanorm = np.linalg.norm(X[subject], 'fro') - print('Difference/datanorm: ' + str(difference/datanorm)) - assert difference/datanorm < 1.0, "Model seems incorrectly computed." - + assert s.w_[subject].shape[0] == voxels, "Invalid computation of SRM! (wrong # voxels in W)" + assert s.w_[subject].shape[1] == features, "Invalid computation of SRM! (wrong # features in W)" + ortho = np.linalg.norm(s.w_[subject].T.dot(s.w_[subject]) - np.eye(s.w_[subject].shape[1]), 'fro') + assert ortho < 1e-7, "A Wi mapping is not orthonormal in SRM." + difference = np.linalg.norm(X[subject] - s.w_[subject].dot(s.s_), 'fro') + datanorm = np.linalg.norm(X[subject], 'fro') + assert difference/datanorm < 1.0, "Model seems incorrectly computed." assert s.s_.shape[0] == features, "Invalid computation of SRM! (wrong # features in S)" assert s.s_.shape[1] == samples, "Invalid computation of SRM! (wrong # samples in S)" @@ -103,36 +80,27 @@ def test_can_instantiate(): assert len(new_s) == subjects, "Invalid computation of SRM! (wrong # subjects after transform)" for subject in range(subjects): - if new_s[subject] is not None: - assert new_s[subject].shape[0] == features, "Invalid computation of SRM! (wrong # features after transform)" - assert new_s[subject].shape[1] == samples, "Invalid computation of SRM! (wrong # samples after transform)" + assert new_s[subject].shape[0] == features, "Invalid computation of SRM! (wrong # features after transform)" + assert new_s[subject].shape[1] == samples, "Invalid computation of SRM! (wrong # samples after transform)" # Check that it does NOT run with non-matching number of subjects with pytest.raises(ValueError) as excinfo: - s.transform([X[1]]) - if rank == 0: - print("Test: transforming with non-matching number of subjects") + s.transform(X[1]) + print("Test: transforming with non-matching number of subjects") # Check that it does not run without enough samples (TRs). with pytest.raises(ValueError) as excinfo: s.set_params(features=(samples+1)) s.fit(X) - if rank == 0: - print("Test: not enough samples") - + print("Test: not enough samples") # Check that it does not run with different number of samples (TRs) - if rank == 0: - S2 = S[:,:-2] - X.append(Q.dot(S2)) - else: - X.append(None) - - + S2 = S[:,:-2] + X.append(Q.dot(S2)) with pytest.raises(ValueError) as excinfo: s.fit(X) - if rank == 0: - print("Test: different number of samples per subject") + print("Test: different number of samples per subject") + def test_det_srm(): diff --git a/tests/funcalign/test_srm_distributed.py b/tests/funcalign/test_srm_distributed.py new file mode 100644 index 000000000..45eb36f1b --- /dev/null +++ b/tests/funcalign/test_srm_distributed.py @@ -0,0 +1,136 @@ +# Copyright 2016 Intel Corporation +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +from sklearn.exceptions import NotFittedError +import pytest +from mpi4py import MPI + +def test_distributed_srm(): + import brainiak.funcalign.srm + s = brainiak.funcalign.srm.SRM() + assert s, "Invalid SRM instance!" + + import numpy as np + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + nrank = comm.Get_size() + + voxels = 100 + samples = 500 + subjects = 2 + features = 3 + + s = brainiak.funcalign.srm.SRM(n_iter=5, features=features, comm=comm) + assert s, "Invalid SRM instance!" + + # Create a Shared response S with K = 3 + theta = np.linspace(-4 * np.pi, 4 * np.pi, samples) + z = np.linspace(-2, 2, samples) + r = z**2 + 1 + x = r * np.sin(theta) + y = r * np.cos(theta) + + S = np.vstack((x, y, z)) + + # DSRM: broadcast S + S = comm.bcast(S) + + X = [] + W = [] + # DSRM: only append on rank 0 + if rank == 0: + np.random.seed(0) + Q, R = np.linalg.qr(np.random.random((voxels, features))) + W.append(Q) + X.append(Q.dot(S) + 0.1*np.random.random((voxels, samples))) + else: + W.append(None) + X.append(None) + + # Check that transform does NOT run before fitting the model + with pytest.raises(NotFittedError) as excinfo: + s.transform(X) + if rank == 0: + print("Test: transforming before fitting the model") + + # Check that it does NOT run with 1 subject + with pytest.raises(ValueError) as excinfo: + s.fit(X) + if rank == 0: + print("Test: running SRM with 1 subject") + + # DSRM: cyclic distribution of subject data, otherwise None + for subject in range(1, subjects): + np.random.seed(subject) + if subject % nrank == rank: + Q, R = np.linalg.qr(np.random.random((voxels, features))) + W.append(Q) + X.append(Q.dot(S) + 0.1*np.random.random((voxels, samples))) + else: + W.append(None) + X.append(None) + + # Check that runs with 2 subject + s.fit(X) + + assert len(s.w_) == subjects, "Invalid computation of SRM! (wrong # subjects in W)" + for subject in range(subjects): + if s.w_[subject] is not None: + assert s.w_[subject].shape[0] == voxels, "Invalid computation of SRM! (wrong # voxels in W)" + assert s.w_[subject].shape[1] == features, "Invalid computation of SRM! (wrong # features in W)" + ortho = np.linalg.norm(s.w_[subject].T.dot(s.w_[subject]) - np.eye(s.w_[subject].shape[1]), 'fro') + assert ortho < 1e-7, "A Wi mapping is not orthonormal in SRM." + difference = np.linalg.norm(X[subject] - s.w_[subject].dot(s.s_), 'fro') + datanorm = np.linalg.norm(X[subject], 'fro') + print('Difference/datanorm: ' + str(difference/datanorm)) + assert difference/datanorm < 1.0, "Model seems incorrectly computed." + + assert s.s_.shape[0] == features, "Invalid computation of SRM! (wrong # features in S)" + assert s.s_.shape[1] == samples, "Invalid computation of SRM! (wrong # samples in S)" + + # Check that it does run to compute the shared response after the model computation + new_s = s.transform(X) + + assert len(new_s) == subjects, "Invalid computation of SRM! (wrong # subjects after transform)" + for subject in range(subjects): + if new_s[subject] is not None: + assert new_s[subject].shape[0] == features, "Invalid computation of SRM! (wrong # features after transform)" + assert new_s[subject].shape[1] == samples, "Invalid computation of SRM! (wrong # samples after transform)" + + # Check that it does NOT run with non-matching number of subjects + with pytest.raises(ValueError) as excinfo: + s.transform([X[1]]) + if rank == 0: + print("Test: transforming with non-matching number of subjects") + + # Check that it does not run without enough samples (TRs). + with pytest.raises(ValueError) as excinfo: + s.set_params(features=(samples+1)) + s.fit(X) + if rank == 0: + print("Test: not enough samples") + + + # Check that it does not run with different number of samples (TRs) + if rank == 0: + S2 = S[:,:-2] + X.append(Q.dot(S2)) + else: + X.append(None) + + + with pytest.raises(ValueError) as excinfo: + s.fit(X) + if rank == 0: + print("Test: different number of samples per subject") + From 120fc09ddb1689a6b4f0a849acf417b0366664f7 Mon Sep 17 00:00:00 2001 From: Michael Anderson Date: Mon, 8 May 2017 04:17:18 -0700 Subject: [PATCH 09/13] Gets same answer as v0_4 --- tests/funcalign/sr_v0_4.npy | Bin 0 -> 12080 bytes tests/funcalign/test_srm.py | 5 +++++ 2 files changed, 5 insertions(+) create mode 100644 tests/funcalign/sr_v0_4.npy diff --git a/tests/funcalign/sr_v0_4.npy b/tests/funcalign/sr_v0_4.npy new file mode 100644 index 0000000000000000000000000000000000000000..6472c68b9cd414a67fb07f3b32af7e838492746d GIT binary patch literal 12080 zcmXY1c{r5c+qdt#7DXs!7fM>WC!)RzAw`Q4Q3+)&C8R7R(S|lEWJ@F><$i2ssmQ(? z`_7md%Z%+e?{&S;Kj%94%yYImbMDV_o|svi9<^m>GiJM~1(`ULN*tXWcyqAQyHL>&CBF#Oa{oGt@!w8 zW)Vtlw?Dfe!A7!5aLIh5&Oy>0%2)XAz)4!PcvF_Qz(q3s(=>EnfrnHtcHg+RiI>!M zgmdMw2tTQ1x9)IlfB?yP;M{4e5J6H}ca{3K1R>HynPvT*wZbIVuO;g5HAF~4{fCBy zghff^S1JtWvP4M-w+EV0%*04_>)J1D^%WyMtPu*_lrKiQ(acPKS|>(&eO9;X!KfI? zB+SoMSX7)O@TM+qey2F;$diSn8$9A9Q|`&tPU7Mu^Y+M>W!uC__R7{f7bVd=SAy%`Y_-!;d|#fu`O zM|JwN{$0W(`|)QQ!cjsbPqR?}XZnJq%C*Zo719Mr*_5mw+E)A|kpRZl4PktwqKyLL zmKnUH-&#($b;5W^FCwI*QjNJueW?Zc8}D+GhMc*Sd|Nn3S4kABGs5g731&p{Hd!{3 zWw^;jx9TNG%GWy-9LR*F=B>6HKh1%)yF+0@_6%g6*GScuqQa)E%;2!zNoe`z{p35_ zD5N?@?~Ib}ho^1bD~dy1@Nn>Oru=9VOx_!RK(Q!;+!y(dXMKFJJ<@iuxUK~IdE>(~ z0~>Mug;&jtqaEn&DL88TupcS+PnGFr51~TGxt%dd<4BJSJQ#I~f^kU=5vsyeocaAT zejCv7{-McGt?RRRswa%#-Z_uv_vlrQ>zSzSQ|hj;djWOQq_XL^7E$|E{`SZ%%jhn$ znwSk)!IsVPTQs%Vh~S`Y7w^xp5!*NW^xrwnPW)Nde_S_=owzA0I&fnr2LYj~L8?(4 zL}IMey-^(Fd0kOiDmCf){d95`2@h)JW9X+vrq!+Zd8GFV}41#>x^efLN!@!=mN7;U4 z6pV7aGA4Z{;ryxS-P)@tFdRDnt98dT_%BszE9Ox_(D_2uV^12`N2%`~*-D4VoGSJ= zt7c#&sONyqwHX-X2;~{*o`naHnRhE^4jPW{Ird}qJj52Bt@Qjo4`J2fGnU&JV5xC# z-Rh$ZV4S=>ed#g-D6&}vT@Fl0xP99oBZ3K6v`9mf$xMjX+j}!4p9x!}?}!AhS%ATl z_g}cCFk#?_Mb{0%1!!J6Kd$|k36^-|uzo!Q%I|VK+_8oA{@Z$!Q?U#X*to%|Thei*Q2qESQ?fV%xwhB=g)Fb~-u>J6EZG++#cg z;p~6;wGYq%f<~*Dd^9LFH1e;$PX#r12Q8ZUG~^8N$Tjs*Ks0}OV?Nsyh&bQ#Zn`rG z+IE}$1{TKPu;=%{n1B)3!YffI_k9R1Cm4~|Y#D@m>ZJ6wTl;`k?K-!Y)&n(Jn)5PS z$)Nn#xYE!1R_P@8t3@OGnOAGIv{i-qYhm*JQiCQ!lb~LWeX~Y## z0h-pQW;{Q!ar54-?U)y+dHvVdE-dB$*|a9H2Sc9omWZ$J#Y*Rd4V&EtaDGFXLrPFT zCe{~^yk9$n@=f1+kMIuT{+7?`R`-Sx7Z{uEj7Bh~zW3X5!YE2stloX$@HmF)eB_*7 z8b?aJ8FRPfB!->vw)%c>0#|41>|Z8NpgecM(+J;59Q7}%2z*FEvnv&~8>mxA`n#SX zAVb03*BR_@3Z_suXyRPXQwq9wcHVJ#Lcw7AR|gqa3f5K3OdlSe!U5|mNzJkpjJD!_ zbk2Z+?52dc@X!ZdxA2N!a$mfSj63X;zPd#KZYz^+w?QuAFrRbIHvq^}nNWA#5e*!#IN)mkhr@-;HVaib_ z3Jk4He(cUZ1%;XqWgkURVEnu<Yjr<{Ub1QOX112&=J^J zmcFmtX9!>-mWzF80Ls%}JA}#%LU4G2nh(7XRMs~bTw(74b?E~`)x6zMBq2h`pX`L; zR=zkj^LA(p91`^Z&;-t>QPko@9n^%FyLCSO2OdKqIy9?auxs|&UonF+_&W2D_3zw*&QAOHR9^EevpX67Qazk`x{=d`IGIU{tNjY?5W*k^atM@irp*oqz*mwE;BJoY3Qk+a$(TQBIdRIN4>M-J9T6<($BO0?_pUAs;w=C6jCA@idT;$D zX<9Xc!3Vr@>paKt&z>0n!SqSgUD(LpRY$?jo>I;;rqjrK<8qAG?`i!0SFWP1go>ko z-rQ77rJ>H>N`vpibkvX592HTUL92!T7Ou6=VE!d9fn#EGI2p-aR}nLZf4)z5m%7d4 zor)+Mz4Z*t9{X{uIg)|l*{0;%=Nb6OEzh&hh>84^Z=z@D40PCGRXjP)z|L(ILj#pe z+;cB1uqJ|uFQlFbRNP|Xd7lUJ0>>8cmE63OdkYh5zjr_KjAml=L({cw;tQBPcY2kN zI}^)2Exg=W_+Y7>n62S3v zW?%+y?rdT@vgYCK+>qzwNyh_1i%vWzX!v;g`(}zP4c|$fF0Kxr#?JR#+pS+vFlpzW zpFR60@zaxsdjkz7@P3v38#S3RED~-}tot;K)9D8m)Bune&-(D+-ltLWszj8S%i_o1M(`)-grsd;_iG~0dTK=#~E?1 z8q_?vk{`C$!@~I;T(@4f!;KWb({kr~!1%~=h}_#g$nn|BlVUUky7ab`QKwOOEqGV+ zx$`8%Ugd0(Po#inm)@x8@-)1^CAhi1iUz-UgazGR%mA+ckmAHS3m<(H4{-;~gKl5v zY{PK|5bM*b?l^*>U{Tg_Op|?#m9fc7j{zjzO%dd)i_AcRnNpU0taa-`$NQ;=Nu#<+pFGf zKR8Iu2jB9TGB`*%=?aRR@|+}zp>Vk~a-5_?8sDBPwsVkHN0B_d&v25o{0yEA_OaR- zTD66}n}ZbLvAEdX%0U`3DZbX1!a=&W#wfu|iG$Sn+O%_{4?8Jo$i)4rK0B%BzwSkX zjh*zuxA0|D%?b#gIeWrw`wFlHJy=pGT7tec+?(_6FT$1IlV?qzF2Gi(QuLE#fMUwn zyTG71*v7o^`m)yyoby-}$FQY=+5Yl#>qDp@d-CwQKOd(cBH)^hT>lv8ft1q4pg~yS zRM>4(+Y7VH6Wz_b+Ti3kU$lNa>wmu|?eCYeVmO!dNWT!0FzCs2a`C#q$d=!r;py3o z%D4RvYJ4Q)vA5By8+H$(GO;gnzxEhTt+9;d$f4k{$AP{=9c zX5htXVHqw=DZlOjy?Hp zahe!WqIGsSqf?A{vLL$CDPN2@ZW0!)Y9LNry6|1}{c&-^a$EnyybIz4!RQhn0&$|e zCG^V!J#k{*KsAFyS)6!hB;^&t#)|(aa(u}rPCOT#^(I*D{A@e@qt#W6xF040ZTrNC zr(^mX5+_87pU>8e4ois=@$2&q`|Lyr?YT?mCzOSWrjD8eDndd8qaY~Idr*LoZe??O z%ppMddq|1Lr}7bQn+Ix}^>_(Bwq`qpT5dvE_I0SzD=xy6T<2uZ#z{D8zxAmwWG6la zZ&cbFxq`Pcu1Wv*egR!=9~*M3GH`d^j;i7E8JuGbxMe0xqs+aJENDH7V+W|Ny@Gn{q+0S?pOGM$yt~(NkSR_&)+zD+{sD2)qOqPX2a6|Xw)pq941BISbpZX!AGCFKRstY{K?-%vMKTy6XkzpB{ z2$BQuZs`mcW6Mf2sXDG1$4-U?TC4S7#Rrnr@4G{ICNqBd?!HO1m^IjwHaU%7!-ta# zZqm{Cd-bIgm*!C7^{Cox7z5X|ym+2|eG$uSb!?^eSMY~dTdHa%8H>Eo2MQ2#M^~G#dM7}*inv5xNf02uEr&G^UKb$xUe`s0=n4>*y5ENn?GPZ8 z{BO8V_VN>2Tdftg`LX&DUY047!kxC4S8}hkb11 zA^KWwG}&|U5MSQ}hD`f&6V+2XV}-6ee$L31mo-E zms(DCVpChVrp+2QLOi_v&+WKn)H?R?@P`A7n9NUll{L&j^Cg?H#&2`zu0K$ox_K5q z@GGlI#nSMVdDeIBkZDX^`AGH?pTt+sUtHEo9mN<+?S%br1~G~+Rlho;2UWzIhf|b0 zF{skv>bkQHs5IoRQeIPr++LNR_H*T8#8Fw_UJN|ap0OfXD3CSf&&JEvpb3@@byt* zRP6*63x8I^4FMz+bi(0_N3d1{1xafzMgn?WCgxwWzb`zR$wi?TIfRL z3Ot{5nOt*k85G0le<5xelB)PM4y7&uDN9x}L3Rn^onC&j3|#`R*EhKms0+Y%{`sdh z#}?q7n(w!Tr%dq5d$Comm;o9a-ZsvS%)?&xgXVUY^Pus-a-AcKZyz4E$tszi0ond~ zJoT3jZT$`UCr4;d;J9$KK7tDQ0wT>Jr>4OvprB7~a0={=PObQCo`kJN)rPMgjKT2x ztMQ$x!yv^kCT&zQ01XFKTK2u`174YuTJ4c8h#NS5g?D`y9QLbkSbwJ(KD|E`==$~_ zENE#(zK*GdRjYSh%u6Le_>DyY$F&qNI(JUw%DYVT{`WA=frO|ODr?pJ`#1ir5mqbv z@ej9ktp40`tN}l+(^_@#MGGqHGWdQQkkRJJ<_o#bJ;?vLCP<6Xhrt{N(&Uf!qg(B2 zjq_~7=p&0c;`1Xo-C}Vsdh;0me&g@*H4J=;i zZvJajkc!J*+l4-?rlPPW*UpchsF8BC7(GnfBi2A}QT=ir$} zM~fK0scBX`lP}XZfV!FMC>W`r$Md)z0;|^--|&Na)U~aMrwPD)&9( zn!;q)+dm5|Cb50sY^c861bXg!_Tb0gG4u_u=^u6;!;L#@-bl}l;3&~hkavFw`*O>E zdHVNb-Ttp08j^je=AfY&CD?<4JT9D~SG&=qkgdmp+=WC*cGTW;?U?wVja|@#cC<>j zv>H9ph}@I|pH@X4)Ic-6Z}eSyBeX>ogcrNF zz&5gU$acYY=rp14ce~XI?2#KUvYB;*UC3pt`T{b%D0-aeRoV;fXNiLO2mO$8?N5w5 zwI96yZTHGKJO~Ey;^$y?07}E>aHL&n<#H;5|w8r7<$9?C=LdIb@cJ`vtw+Z-oalppt*aUPYJX;)%oPg<% zFFnn#Ov0U3p&csHlh8WKeKzyPB&?rqZ!W$&2|w_jL-yj~Z5rDrAZtWmCwu86 zK#G<3+Nv>lWh$gHy>=WfN6h%1Uo{S%S0n{O+sELD@4Nl{cg8?+i=f{%^HF$MH+|T9 z%?JqU1{a+@Is)0$ytM8oLqL*0*OzX>(vzGzd6jtv;j!4;a`IRo1cY38{ZhXdG~HxR zR2%ex<^if*xjGpb$AX=jM9DBlc*Nda(*e^~-_I5LH^GSpWkZR;M(8TM+%+uS2#uVg zs>_!3V5yz{rbeg|{C9cxoR=Q>k zyS-U7oVi<)T`@*QdG+r6kS$cyU#1JOD^u~3bz%}(m4$~Qf1BBlQSjYW%8hR<{p{yl zoxzz86PR5)^IY`h7@9e6vAVTw4Eat|Y}lWVpcCbQdEnO}Jo=Ft_4@Jv{=IX~{CH&_ ze#%&Qrgxl-hPu4=lNp`ZAYw$Gr8Z%^V%;H$)CL?8V@khOsztU9{yN%$<=Cu|wv+$K zH?-V)orZ57;ju>{4Vu!Q;gz;z%OOJo?1NM;7|Z{Ld!PA#e@Utblj|p_KMdR8z!B~= z=gf9k@9@R@HjA%bWe<0dYZCG= zAG{WPlLC(WDY+vZ({P^8-8d$N3Ky$ae+l122kxnN+GN`qSY?p?>d?zskY61sO0AlO z<_Hg;PqB009d__Wcq;==D=p>;3b1&r$D-?{e@sBvcjd^p2n@E;9~&79pipH08{HN` z>0ZC>pMoW*c=c4y^vW{irQO#$lDP~;8mi9Ug;pRYTtnLK_c8=k@H-AjYIB>KGK6_3Dvokr_2`oLWFeuY(`-e#o z=kJnl^&W?Voj1&>;Y09mSj~m>wjYumDSW#Od!cGTfAo+?8$9~@>3Qg%8hBDHZ?C$( z6k;`tqEv3bhndGB$W9cZ-&VI(@$0M6G_**(C9w&|ZvC)X*i6Qv3;k;y>-sR^7{<>v z4&%8E4~3FC$B=8wnJ~i?3c3wzPv(_VF)`X~{6h~NONLwPZ{KUZk^G#BU0>s2uzb%eV0>rcAhBUsj0z`9* zNvnOP0CDZne3VD30I{VgL5oRa;r$(id}xpWAs-fGulrPh(7ZOI`oLO%C_L@o7stZ; zkG65{c*^7__C&~VMu+kf@aJpiU1@$oN-Ff1YbYOaEk&qC@&PX)maYEqY#=YuXdlK@ zy23-8$vq^L^o5&vGx4CTDu9bP+w%GT`rn+y@_laVu0EEp6;T^~W{90o+-rD7>I)lD z86BTl_F)-YciQ+1Zdt^vu2**Z9x?Ec9#7Q!&{>xMbTi}AS=KpplV7h$m&Hp{w!dzx z8^`fUb607VK`fDxmphu$jpY`3)wzrRFw5kq{nyBR)Oy7fwQVhg^k7=&W}5~`>J&QT zTto&=PK{yR#X+ciL{2BO^4l6WfN_sRmHZ9{tS-)c!lS$d8yKxJ zDbj2tjtjrS5(#z^mE7Qzpv_5=seNjv(!fRP5}V#|;XDti>wAFaLv=nu^_ZyZ1A zIK5@DAVPo?yT!}DTSthbluM@X3l%1*gjTyIw}_C!tF_<#CoM)QsLsl%e=0_@^>wdT zUSavX@lCwpdE%r=_UM6tX>rov%CE0p7bQrP8s~i(T@s|F2Y0Q@_ehdXz;-V7Tau(- zffYV&o{}U*ts3sqEJ>1e+pI`o;0oS$?A3E$wSuW~rHa|c%V^0w+eV;=G3QeJR%doidHFeeJh^b?z47o>fe|vaV}n7OP?(CFZdB*R$-w(AcfB_*#pd z>scEeU%8Z*eK|UVUTJ$0ON!~pv)SEll$D2sX3v^w8q)E+`@dKr3JrJon`%h-(D2~N zN17AztoR2L{e5Rx`Ne{It?C*LO9GuOeI#hubXP9T{4f>Y$}%q2=Tq@Xsg3fNJSqk` zBo1%=$->ENPMsR)rec=BM%fLQsJPFIFqWRC;_HRoS{-UMY`@t`Qr<$tF!$UJ+aFYv zUZwq8;~^D)%Q~8#WA!H@adrnUE3cibxUIsgL&F}sI)ffA8uFC+Qp4(Lm{=8Htm#3+ zWb3OBV&v#Jv03j{QacS<-j?XgaT;!Fev|f!VCC7bK4+vq(Q$ZOJyNoSj)TLQJ80@N zSnRs0WqI`s?l{-Dd*8bmH2vd2G|tT6?yk(8K4CNX{(|F){GJ(9A@eEPUz^3q`rIp1 zC+09Ve@oel=Nu}yyGOO=&Eeha?Lsvwb2tnkb7>;;c+d13zrnvbY-JMs)z-S77<7tpZ>={2HMYB+dYQ#Q%Uu+0UJ)ib1Xjuk9cQ9SQ*CV!ORuh7Sf9ErkBQ0! zV(UzjnYhLD_`%;>m?&P&S1*2>rIQoVfpY0g%&#>M{BfR%)jtg5ukB>wM!CGBb@fc# z{Vh`QG3)zMHPJE;CmFc*%6pDy*$njF?5=T7fr*?yf5pX>Gw^Wl@@4;P3>44uxO$mY zzevBKI^V~@>nk>YWKJ@$E>NgPf`fr{#lVStPX)S7Xtc}%wv%ewC{ zk85u&`_(h%Fym8e#}y;idhtrRoOWcL>jz&SQEi*Wl)a62(+O6&XeIgTz%1@%&y3VK zI*UPSLlytHXOZu_$^uxFfmHQTuQ+me7A%|pH87w@tIp)yQfe^VrVy6cM|8L zbq`gjPhdI+*Ny|bCUCz=?8G0oaRkojxzVsObmDl+l_@uhk*)R`kNiiG{WQrVKxh=B z6jN-w(njzz{lIlqn-N^+GE_O@G>q-#FZS`;viwB5I<`{7L2NgPKJvJ45Mx@Gn^UU? za9g+hrQi(%7{4^lJUY{dcI=wV3v=dW^pM( zw-0HLKP;}K^&q$L$EvNY_#w+C&wU=f7%z6Erl63FZf+K*d3bwq%r1)0*S81b3I>H(&%_-(T#aB_sEy|d(fHB@xM7_eJ&EXLgpajVSm}$4aa-X_soUI z6{%z-yT^2$cOql$Y^jpts}2kb8xZBK?L^I$O&`*}b>r~`acu{^PLvzebtp{jM%x16 z3b)&xXzX^gDN~<}?=q|#77lfyuWD*crd1b;nys1IKGT6S`D}#=N?n++sXuOyWjCw6 z>hX>H$gJ}e%7-qH@q?SZN$e>yu2L5$na=9Ql*?q{A3MmX>1EmxB-4e;?FILYeso}` zTsrlGM;GF`Y$37_8LxcH6;3=!#&Bk4*rvzaU0AjLd!4pZ zH`1k74V^pNg(~uQ^lyoF;kq;Wu9>X)J}au}3;)@Tk9w?+we9XikMFyKWaTUF9X4LCXzr(@TikB3|M|S75AaPVxSN(Z2RvH_qiRd?>d?N4EG`SIF_pu#XqBmi} zwp9oFbepjEn#r01s}?M{-KrRn(S#G)luv<%P1wyDOg6M_zzauC1vT;1C)Fo00?YGd z<;M?XK-|l1pD+GMg(H%>?gk$+Ax5pY|IAP}xCAV{`IeLoJtAzz`M{_82_Blv@Y;OzyoM@UDmOy?Gm| zRO`T*^IG9eomx=)p=EmPM+4kGQ%Y0V+ypx}oKN0W*$5WFtv_#{Zi4dcgZH=du=x9y z?dGRsn__g&m^iuA0)EQ}jb>N3!}pY(+RvR@U~A*irjo{1P>s%|&l$D@ z^M&K;@Pu|q+_<^@v2h0+G|%tk{n!rjd4mVexwJvsn2vIYAM2h##fisGtp091NpF-sO@?fTSI?B$ z$nae)uyl1g8OEi>5_-ON!`1}TsCL%4eU+rRZ++{8o(R$9%q?UvRxw+6Fx&-UVfn|e zu5?1ZWOHk`4;d(-YEh{aGE~Is)V57_LubzTqthW>;PlO_QR-wDL;C_CUV& z?d^O9J6Xl3}&MK<(< zqOXihN_ii69SX{~*MqP4Z`iInxbg) zL13&kcm1<{5Xu*29NmV8;LV+UPfoF6p!k2Bh)y2{9^PA1=BF8 zdEriQ?=&m9lq@UzQNeh&`PipADwOLt9m~<7!NwQN8p%u=g#EdG?UE=RL?mwB5M|O} zv}EA2$0Irj-MP=W_kj)@F1gO|U!nsg1KDSOw3W13I3^PwdEa3cBtKI9ym!xmN#R6Jv*tWRd^;6?xML3DTpV-> z&UsLx-~4@orT2gOmNo0g!2rF#!jGpo8KBL6%sgC?0V~hbGI>`qz%X%JobY=FXssWz z+a<_^nAaEAf8NgohoNc{??X&@SL7hn@`C}x&EXqXcQ7EQH{x!HHWPLShm;G@8St!` zhmm-P39H)&7|n`I==if{SwofyEsFn1_4zYlY@BQBpCJaw1o;^pT+M{s3t|*~D<)`% zGp~Hv!-Q|P%?V1=3{c5&l5l3#PySq`KwB$>KvN=VeJfPo?|0zo0rA32%Z@QcmGiZ?CXny~501a9h|D_yf-6L$u%n>~_O@;ldL*t$*)1bFh z#dkAH*AkT#FzK_Rg2ZxYcxM;ue!!-p( zOOKcP0u|gY9o=5WL4&QVpUpeM@aXx z^ef>GO3WT{Cj46xInI{O1bY+9hxflOKy?2mrO}23&~Xv&o!Gnx>B$!cU3ixuBHDgD f?avbM+Lbw9d$SBYg6g6p>C5m?d_?M0`wILIIls)W literal 0 HcmV?d00001 diff --git a/tests/funcalign/test_srm.py b/tests/funcalign/test_srm.py index caeae8cea..7c5a04464 100644 --- a/tests/funcalign/test_srm.py +++ b/tests/funcalign/test_srm.py @@ -20,6 +20,7 @@ def test_can_instantiate(): assert s, "Invalid SRM instance!" import numpy as np + np.random.seed(0) voxels = 100 samples = 500 @@ -62,6 +63,9 @@ def test_can_instantiate(): # Check that runs with 2 subject s.fit(X) + sr_v0_4 = np.load('sr_v0_4.npy') + diff_v0_4 = np.linalg.norm(sr_v0_4 - s.s_) + assert(diff_v0_4 < 1e-5) assert len(s.w_) == subjects, "Invalid computation of SRM! (wrong # subjects in W)" for subject in range(subjects): @@ -191,3 +195,4 @@ def test_det_srm(): print("Test: different number of samples per subject") + From af5a756bdc9616bd2a0bb53e6b59f6dd9738664a Mon Sep 17 00:00:00 2001 From: Michael Anderson Date: Wed, 10 May 2017 02:13:26 -0700 Subject: [PATCH 10/13] np.allclose and fixed file path --- tests/funcalign/test_srm.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/tests/funcalign/test_srm.py b/tests/funcalign/test_srm.py index 7c5a04464..00a3445ab 100644 --- a/tests/funcalign/test_srm.py +++ b/tests/funcalign/test_srm.py @@ -63,9 +63,9 @@ def test_can_instantiate(): # Check that runs with 2 subject s.fit(X) - sr_v0_4 = np.load('sr_v0_4.npy') - diff_v0_4 = np.linalg.norm(sr_v0_4 - s.s_) - assert(diff_v0_4 < 1e-5) + from pathlib import Path + sr_v0_4 = np.load(Path(__file__).parent / "sr_v0_4.npy") + assert(np.allclose(sr_v0_4, s.s_)) assert len(s.w_) == subjects, "Invalid computation of SRM! (wrong # subjects in W)" for subject in range(subjects): @@ -195,4 +195,3 @@ def test_det_srm(): print("Test: different number of samples per subject") - From 742d1b279f43f23568a43da84829b72b9b7988d3 Mon Sep 17 00:00:00 2001 From: Michael Anderson Date: Thu, 11 May 2017 04:19:33 -0700 Subject: [PATCH 11/13] new npy file --- tests/funcalign/sr_v0_4.npy | Bin 12080 -> 12080 bytes tests/funcalign/test_srm.py | 1 - 2 files changed, 1 deletion(-) diff --git a/tests/funcalign/sr_v0_4.npy b/tests/funcalign/sr_v0_4.npy index 6472c68b9cd414a67fb07f3b32af7e838492746d..c9426d92ab953bdad7dc8794a6196c125c875f7b 100644 GIT binary patch literal 12080 zcmXY1cQ{w?`?q=9JE2b~MMNZR_fgW45;7vALYhQHLWqWx?24ok6-ideee9jR_uk|E zw)fZXx_;-M`#Sd@=Umr$p67M$*9g$oKXcBAo=%hQ;r1(6?_Rdwen4XTA*+4c}VkYoL8$5kX}T+N5|JU|_UBJ>63$Xv4xjemB!aXRJtkgr z5(f9!jeq>+B-)3?zg%48BsMG6C+v;pB!p(%4lt ztwS6{p1Hu~*a>!m*|unvonR+4e@o1Xe4!8}0gof1c2NjQ(vDIzVIwY{D1Gmx&q`z* zqFmp9jfJpj;plmJkD1tYYx8R6V{8E6!SUD>X`R$b808|K z{Z1K$_qW>5=V=VUnTz$&V)gCdBtEGzrQQgqk0|Eu?97MC)@u5E@lZ(WcDZ<&Ascqe$YHzkAVW=h z@?zf6@QLdTWXy$A@{h+D$m3jR*p?d@NxI7Dn)W+PBsZV*k0(h?G_aD(8?dj&R9?rOrP zZI4;W>|V%&a$ zbEe)dKFUm{|2DdhM8<*2GBiy|D9<`Nc?%bKOzcSj0dwkN#zBo~0*eA}&?D zJVZ|h_{a6N`q7ceRo8cvZdyY(`iY&2@+(;P@p>W0{1T>|jbtytMcimOIo@?|9&0%E z^0#o$p_-!bGdJxS)O6d|^!+ytmFVy040cdaa#&Eh-)IsKJ}LCgfD!cexK}scKZNBG zUw<)p_Ttp6rrD;IE)-Ag?){L^j+avd-hX+~fX|FE#ZtEd>l-aZg_%fv^YqA(S-6IU zmm1DXd`f`rX+*KKK@kLfbWD61TnSDDeL-?X9cb+0wTTRAffKSKU$2FCK)OSro0okz z?74UT0WWti&`B3gZn5l#gGbzlf1MqMHy43&CutP^8i|aSDvv|V-t!YljuYV2d_SH_ zM+Lsb{mz>nQb8)8)w)`f1{%*p|8TlZfs4Is1MzSQ>O;ZvvdA>NY6z(dxH1j1zoRx? zx19zDg;6^Db2AWdl;d`F+6+8m|B=q|W(HW5U!3;KpMjiY73#U<8MqN=aWlwm20~}Q z-;Lv%fjyfFzkj(l15^4NM$Vm?0r`~$cdycENY0-w7LA#Pq#%=qvUk&Po9jdFXxJ2p zovLm>nmPqKEjFYT!xZ=%S2(B?QQ^SLn}LHnsBm8A&qk`lBs{RF`DlH40xmKKJu{3P zgJ=9b6JPU2!F<_*&-dUk#I_yEG5c=_Di?Un9*quy8|$Nu?YsK{X=gXO9_oeh`*c!X zkzJ5PG=$&S)B)Db)e+V)&2UFH;olp*dQiDwndqH{qHEpsf{nIDh;23KEf=t9qt-p4;vghN^o68Pv2_?uRv2U|lbvfSt zFXqWYUL~3w+k5yrPb2a>jZ(rZ+E8S)zBx&z6H^}^)VBWCiRqfXiqS$nSh!PAODees z%c#O{XG8i>!gRUEF`*xYPud-5XdA>Uf4NEJ{lmD~Ote(^D9_9PAj3bZOt*REzNnD!yPo+4061O}vF^_sTiQ^G^0ZMUH6x@`RTA)kC{rkG2 zRJdt4H0)P+*pP-r-R%;zTP69L&H~ZwM67=Xh^85mSl#}aC4_M zN6L2^wy>K!4A)R`@k_jo)Cv{bM&*e`DH;ksk70Asox^=O>Be)YAW%mzz2&|k z)}my-NUzKF8g%}k867iOiB3V@?ELQv@K-}z<5n|7idlQd+GH|{?D$c8{Z#@|Qg=J~ zS9zhJ(e7hGvl$@1Nmys3m;~u#Tkg-EL6DybDfOEzhk|PxWtsac;p27X)hv-(NdDMJ z|9ZS0cA4LN8D7&2(zZ&*T?bk~_==Eqt7r%4owwXN{iqw3J{1XV^6P@~y;uAb(|Uk< z|8SMgy*`L!*qkQn(+B4&yi07l`k^xC?@){T0BDFEW0RQ~hUngy1iPRSFuJJp_nGww zEcO-HpQsrD9oIDDm-eGD?mux*T3{UX?4IPO7>$D6w`3{7o(XX0^?CZ3Iu3Gen*<(8 zPJ-JtQ97|(lb{f0AJ=|j0*)RFqN8?ALMh}kwI)u2$S0vsoN`onA|xVqY8w^A21er> z7bn0mesS<}Y|n9Tg#J^>82C2}LLC&8o7=gYIS3AoS}XG5`{fXaKKY{I)I;IYvg zm716_P<%Od{P59nP!-SNmp31SfVdJDp|KIzU3LDbLE#W+otl0ne|i{B{r32Luyg?W z{+uPvT!ujRi@7gvZ$DH?UFp6T(hIeR6b}E|-wXdqwz}Oc>IUui1|?>+F7O%QXlZ6@ zgHm7p6%*4|p!@q|8~kjBOo5HKL#7`7a(i{WcvlJ0pVTv*CX3h`?n+mi7AJL%bL)9+iAQ9Yu`0HhJ91#Ga9yt8RwH$jjAx z>E~)Csu$}G424vq0Q08;Cb33LYSnR5&}lznyY|-)<3Fydw_2x0P`yaJHZXn^V|BF`LT5%% zwlQZ_IBy(R)+8?_DooYh`udT=*;2Tc=Y%b+P|ip=ir^j7dnHVHW^LhNx=y%^V`!X)i%wm z_;eagzs7ygK0J-N@^o2+!qX_{s&vMibs8^k`*Uia=@j-n{CSdG#}A7xgK_sk8iv0* zt~hN#!z7bZ6VGKTZk1$`e6B#l=bOwIuIyRI3AG^=9mPr9wBdF}*xgBN>S8E#IW~@U zk^c4ELZf(LQqWa?-!S@L`Pf=^dk~evDvj=EOY!5-b(M#t_pU^X;QD*xRDYhI#nVsrc8uhhL~ zb_N503=-RzSBHTz>3Y3Mdkhwb>PJt@O~R3aZ%(?$X|NpocgKFqDG1-h3Jq#Aa1s|H zZ#T`tys9|iraKSrYTShr8y2DZu-NW&#U-dH_I(qI%W%-{Yd&Yj3M2_|1xWo}gUjAO zp2X?V5w9Lz(;h9PBNn@Y%EHa*2@#&J@4r^l6QNA40z3H`2-f~Xg=qqegqw`acLNF| z!KPp@c$jn`B#*B$*of9h^&0!)U<5fP%@h}s4Ru^}- zXfYF7oBM#807>p#VHo^suJ&V0g~ zfhbpto?n@vC*;+)y0UoF6A3-+91)-B2xZnEN#|78AS#7xoZfmB4rFQ1Kdf5@d;VvY z%BG7TN$&wl3JVYwdhW9<>m2B3ruY1NFats+AJmu*PQ%#EO=rXHsPHdPRh9a10^BB# z{+aU~f*(txof0FxVDi~~>X2qP{M_ec;_#>$DkJV}zbIM_7v{=dOgm)54X3rW*Rz3` z+!vXD)wURGJ3NO?>}xSN+%W88XB%=1?A7o*y?#zxMUq+H3}aryu04h-o!t@g`Ry7`bspEF z^QR*-*@8v`1L(n&bEFOD&z;O@3dsgM8S09L+ zu8T2f83lExnq4Jlsi5*#;9}IS893%?X&m%;0X}~bH@Ez=0+|{&Y)K;HEVL z@q&j!Va(+tSbiLBr_bUej))`1w^3XjtWGB-3vV)50Da4*yv6V4f zHX_bsUvZN#E5V&UdgVw63sIr_Ex+EEna~Vbwf}vGiFnliFmcs^kzjY;HQ|uLKumYl zju_sgCuqy>%-SB)5o|2YQCt44g3}n^)%%XiaJyl@wb1G!JZn!J3YeLPdT;MbsTy<8 z!R0sdQF;cfO?GeU=A8lwN*t{qlnQ#K$9O}YPJnY0>vTf(2&{2UCQXbCKwieh!o|5B zpyRt(bwH&X-fN5RvE14Q)9TObuE*Ab*X(PBEehozZ&P?L#2^p2%HOF<7DPj2kK>Wt z=RxQn&mc1un1M#kOL~J2g;>IFnIbh=j_+SN7)$HdVA?X##?(`f$uT=`_TFwrlhyR6 zn&%WhWcuu2bcAm zaeH9x2JXuZm@1`v{=#YvUZh{WwOF$gHKKSIV zbpLB~Kj^p$w?2aIL7HERA%b^n$UxR$G3f3VOnET#4`#=DHV@I)LipJW_M`gEV7&K& z{^5x>ctZZRnS9#`J+(tixA%3!)?cb;G+g`OYwF-tvHX6h*7W+USv?3JU0Qc5>x@D+ zbC9b8(>Sc2eJrdTH3|GbALVc{P@(nnJI$L0G^ncHD{b?E1{LyCLZL-d;O@OzEVF)& z!YAC^gu2)Hq5G+)UXHUcUL}(Ft7#6#3*$qg7Un?oLZGh5?ggNIt|-10xd0!w3FV$? zUxfMg=+>7XmcZ9sFQ!Os88SGlGUIM6gXQtCMOxr8)ZFO05q)X}OfJxKe3@JUIfG9e z9OEl6H1VEC)Nd8oPSj3%M6AN}z&@jYuIu%vQVILQD#S~7Fx~1~g#bNS>%i1CxE=d% z>%Q}A@bKiGjL=7Gu=Zmq@j&Pr_<5|ou1j2lQ2l2|b_lKkYvS6N#;a8rx70uV-E0-= z+BrXIzgYpBMc0R?AFsfn^18?s_Z2wY?Kjc z=`3jN{b985b_Qzq2yySVpN2ymrGlPX)8MzTYQ4jC3WltyWA7ztFmS5pz&AxI<_-4ni~aKQ_TRgks&DADcq@+KMb43#2Q`d2jD-u zr65tgKJX=z=bx?5(T1LrO%6}HU^(^28NKFCcthJJ;Lh3s(j9aGVM8rYx~KQni3cs9 zCNkHsYqkz{{Q9{x?p6a{uj|W8eip;Vr5zn1b0v`KRWE*IH4i>a6|gUSNr4A8$F6k^ zq=5bMIp%x$5ilKfoqo|V3>MysB;T9=iJu!!WL6)H!MvC~S)!3S*!9p?hEpmV5B$=- z+-RALJG+cE>+a=ZlKeKJd{+@NX9*qFn<_+>KmplBr3$2$kk6$BD$$Zff&103Doo$N z|GrkH2BXHRJ_H2UBJpGgz4P%JEXYV93+8GuEH=6|?ROo5#;g-a)T7@DvRLxg;by_{ z`zfY1_{c>-V}M-mw<^Z;-Kj&7+T+W19QAm+bofe+K`nlCxx+*MtPYL(uW#YLQHy_@ znm6oHt;QW86K@X9SFhs?lZ1nDHM$E*zBkxdgBf2$a~;lB;osGQ){ph&IF{x9u=@8u z+%{F9n6y-k36h`FMb=92iQMNGnGr?Ed&^nQtpU+L)g#h(F&Y26WbyEfcN#uXe)PuvHBQCV!?SEc6`(0i47SUCL`^onK2 zjSc>Rn{WO*!TC7>vTZ{3sbwUnKTbTYYFGebZ>Pg#-seH!+wkj?C-UK4kik)MUnx{Q zjbpPqUJB`&#$01|Rd7b~+vL~qT6mn38F=_wJq&$qV(JxZgdj%8lm3B?V9W71D<-oM zQXWkm+QqMfiVFFE|l zYJ)BQ`f1x8+F^kE?CTuo4xpMUW*g>p!dVu&$QNdvaQmc#+P&MIpf$C-52`xhNT;qq zxl%WzsChh{QtE=|59({^_`6_Bg;K|JpDrlo-ORga*98|#^^Y>NcZ0d%_lJGlT`*l7 z6sO$W1&8Qz&xfn8Kj**g>MhJ2a3_avp!#hGh+TMDNh$AyC1xA{jI-?!R<)Y^PN)^~ z3K)!-x3)vO9_Yu?w}Vc#HNSsID_nK0v){U<8JgcKylRbags4p>CK+Y*@Slz@2e(Tt z_`k|LAJ9nFpsiQXaES z<-$SMc7A*39Eew=@W-6Zf}tS!QKbt%A?@;_|EAd}xYs2!#rHlEf-mp<#vdJlGEev0 z4eyRYN6mv#c{ftgZncG~kFp&}+`NM(i?n+I}~)8BeAj^KN+9hR0%Oq%Q__;-PUh3FppE zEPYZ}?|i)*BWZ7Rw_5e$7Tp_S1#kOsp4ou=;eBqNj zX%ooRzpU@zI*F~mUjmz)Cvn-r!@@pm68X>jF#pJ*;=a8~e%6Up>}R=Y*$_`kTWmHiU+f3dy3Azi7DWLdF+rEDif=?DrLx)3DrvU%$>^3K!mP zrJmV5gW7u4p?AkKS+(r~ACFgY_CeqhF-q-|tNKAEi$rD>j9MeNScl1UvQb5Nd z*>)786nK)j8iz5%Y2IL`&Jad21ur=^4B@EoK@*nB0n|k@q{Evx&e1|J~+K!rv^JzT0I|ERiebxL-p@V%CLfkv(ST#6o^M~>m z=WPfvp2}p5-25x&C*K=9`gM2WaI&jP*tXjD**Y@ zXCC;^ z*Cf)h*&Br5x(>vn{NK%eQ&{R0%zOXWG=8wi81UOMgAwU{`;OGiB3IW(3G2BzRQ@jg zB1L}z7mA2ZM$RQ<`?vM`(#B;RXbsNvj9$TA?;7H7POf53Yt>f9g*A+pHl|IHbmXmA z4z{{S^yDdhPx%%#2GTsE`lPEb1Nm^gNYhw`k^KC|S(T5*NRI4}>CZaIM83#vz5DSs z6S>XQhNh6rO#1CSVlK1pSE(rszlE*)Nx8pviF6%dC66%>zXzVMlDjQy z3bxO(lG^5@OwzB|NUJE?^DiVD`CoeNXYgVpSDp%fDUD_$S&vP=I4(^gg{z&E9V{qh z5`+KOP!|f>t98e!(w;&_wy#96x>LxX;jd>>Z&FCnE8|-G11Th@yFBT>kwTtjyvEFD zN+Ij{N7$%h6!Oz=HSdxNHuCn-i609-Y-HHO4c!vQ*vL&M_6hUUvy%7A!oy@vtov*B z{XFWqEabW@D|Duqh5YwbaiW%)h0HwFKYTZYnM{1Mv$b&#GkJ=^^3n#Mb$)*++wqkr zBe|W}(eF=TBxB?Z;xC&qkgoYl8~(=7lWChH$Ez~w$V$x?eMxmX(rGpOsFKPmQvdR7 zwBEOZgux4YueXbswk1BY_{##`R{kz;*gJHpWp*8B!nr^6Phkj8K68ja4SOaQ zI0J1a(QHd+%-H8KOrxew+>jl{ogY~|;eVSVB#n9Sy)3%_s3ep1WyK?!O z;rP>QkKIJu;ewBXOucS57#Es)sh;hF>=G&q|Ly^pyvZHdBRB-vzrwV3IFCU3zY5yO z^)YC2RP}OL*O5g=#8iX6O#*%Q>x}oesPNOV%hPn63S{;~bA$Fe|8Oa&K5~l&Z`SO+ zx7}kVe>{1g`*!`i+MZ*zy|SH=lshGx5hlz)`ro6y(yXK-b*C?B+-Y0EGe4XShF&kD zmd?GN#DoP*&w=^LeX|%)ZS?hMAPrT=ZpWrrPU2;Q6A?xoL&!yG9NkFk#hUA%Z5iLT zp{ay?h|i@;{IcJM^v2}0GHgwXUtlcA^Yb8ZPznv;8yXx z#%cpS5p1`%Yt4~?usQvcFMcy4VQC$vcyOK9zxA}+967>FM8%wmKVHE?Oivy6PMBpS z(!J*`1)i`G>$0j%OB99JC42HzVg@@Q+A7fys>?x$BtNQIY2hG>oER>KzU3rlcZ68& z80REBn#J!Xedi+b+BAg)Xk0|0{FhJ;3vS|oS7h9L1~i#)^`@$#LF&v*#q@Q=2C`nU=4xgb~HO+3WjR^QdT zRBnRu&PtelCpU3b`&_`a_4g#2nz-?Va1jApl`gCEa1nC1ouW3K;3NcseJ&MoauN*- zo)#}u*K-P#{f!RM>_n*amEOy}6r$hToXwt#LO9C4yep~AMtHRNPc8?t5}Ay<5Aj`L zA?l(V$2L-#iR1dL8aLPVN<1g^gK>$G*tq2j^Ennq;@nZm|Msuzk^M4yvSrzHMESpj zJsC}F;I71${EuY?QWfjnx;dAi-g0F|XL=s`Ty4zymS!R6Vx`02f9rf#EobJG^b`n5 zz46mNLF=xILD>B`ayXyT2dmTY$K`MZf-{#i@m%Zd}tnt+V z?{T)@n^+6LDm$Q3Iwutlr`<{uQHw_5#};bu{uH2raBzK-X9aFh3+LUvw;oeVXGZ9G zTd{wA#<_E;13x*Wwo7vL;o+4i-4ynI^k&z7-fTUH1}&C(@zNs*{39ko??#cYd5U!J z9mhWDh+|RqlX&csQbeR&?LF>SjKL{f%}S z9;erH3%9o~b&}IKXwqPn)-{ciL81rZd}i=H#|hVE`56p2k?ouCb_VI&&07qzrxCLC zZycMSMk-Z$B28-sNAv5S3yV(UKOI;9EdFV{XL0bMe9IJWdTh5Ann%Nzo}mYo{?hQO zQHA({{`|s!(33magSv30g{}aaU;Mpz;m9HzpaMuww1^;09GM7|-edF8((KF0BYl z3Nz-yfum$-d-p$ZI&?rFP`4Hk7<)~P8{xOtW^N1JHt;Ov)0GnHg6d>(vGVwCm^IjL zAsgKb52HnLF9{BU(cy0&xCaNpr1o&C&$@rtoMJLE(>Ds+9{Ki#?i+{Md-W9J_#~8X z^H$-Op~5-B*>;x<4VcUqdmde*fhpzfQ6Hgci0aT6)ia!ieRan-FCLnPY&EXO4|mJ} zo_Q1#7cdL=7K4tLAD)Gf80Vj!wzE(slpweJ&n#rPpMQREU5_c2`BX)8%|VoB*6#|+ zJp6cS?Hm_62SE`LM+_Y1A7pSY&q!~A-FULj)DO?*SzLoo99hJ@x?qe39o8< z+08>_er!z`-8`JzbKU;-k$LE^OlZpaI0rTlWzH>`&4F-iq-BoD9Grh6w(#uQEYOy( zl7j6s5WD5IX-Ult2uHGIe%77=y?;Bs1uxG))ARc&pYx`mVW&KCEollk9{7s7EYg5# zY*F)>LW9Ll1#c0pHc^pn!9-C!%xfB5hI z4!AMk)OfkS63E8Uk$ln&f{FI z%aRI%tip$0O8$Tsg2K`#oui(J=&#=`&7WC|9A6s)v@bSdO5+QW;q~0fzax3nl|@yT^?e?H_n6@kuJzo3$ptIQ=^ngGyU2Rxa}RQMX7Ged z_u+}5ISIDEefYHHPxOVq{rE~{i`#vr0X#D8@L$1?A>8^sMQCzx0EdmEJg3G6@TS%V z2QI$>bUjmEEtfQi$+6mzohJv;-T8qEbMp{>=Cw(btsTUeV=E_4F%RQNtl{Ms9)ma= z$ai|gb`V$lZ{K|McK|S^~Wob4&KZ!NeYn{bXQ-c3_q3^Rcz1b~J66IbB}RhT5NR6a|&Fpxy(w^Jk^%(Pq&}e@93I>fH4x zV{NEKEirlRymyr-;H_@a$5VmrncuGZ=oFyCrTBOL*K<%-Vv9yz(I0HkRNv^bISKs~ zoVgzB#9{Sqjqc>SaPVo;H{AR%9ljQAzec&41slE>>4_Xo2gN;V8b&vA;q_mQsg!RO zP~;Tz!~a($sPg}<^HHmXbL576eX#~e`SI1)zM~0}ujdtM1vEmEd!F$#r4C4B-qP{a zwgY7Q8!OM6cfpQ?|2#CMy21B2ILdK#Lr3yzwuwv+_!dQm8Lmqrg>luQ!`b~{d&+xh zM57-@K5cVO)foWSL#AhSy9Xe?uh3xj)F9A{`h7Lc9Ryw0bIF|-haiY|w`IhhVR(N< zR=8ng1hkGAOXTnk!s)yl`{Qg#bmk}@-jIHL#9D*=&=Q$z&AvmmHaPaPfL8whoDQd_X z1cS?bNIn{bTjJ8=YdZaK|C2f8eN;dA>vi5Oor3WV>vx+1j z!c?a+ERo-PH(9U>S=Cto7(*RyYkPLIzo`j@zY5XyPPF4SKFZ^7Bb_L-`>MgW5{ky>~M%4N11Kc?hL9-;x}8){OgT9)yqXU!tL-P0=5({m`T(q!ds zHIMygDdyjz7jQ{`wN6ug5r-|s*#|lokwLf3G+THXedn8>3kohHr(^T=pYF?;BOfe% z=GF>+GU9VOzqW#&^BlwR=T`AEvvTO$1z zPs-Kk$gx27WaX1|K{G8rq4h} zZYV4KwoQ(Xbh&#k@UizA9t)m2+BOwDLA1$Hrx3@jNVU z?Y)fLvfQ$P%*(iVK1%-1mqomDs_~W<<06LdKQ-GDvVe*oxhr?Pnn&wDx~1FuXHnWp zaN~B9S#&OJ$`|LG!I<~opKbmzg%0AI>2uffP{F6yN@BYwu%)Op@ub}ZN*dhVcz$jK zpX|BvR_^B@hV|XFs9wh(NdZR7=Exr0&_A(nb{%KA**RG=Z`NZ-SBa>`QWdh(kNw~v z%Q07V>gGeGZ}`;gZeO}sF!X=eS?JT81WRkHx4xH@fuc{#m)Lvtu#aEn%-Fg0Jj$}} zrnlR>z|?Nxo6p-`7@RM#99`=NWm`#AN9keE?!6F1Y+TP3>=f52Sm!QaS-8G2|F7GTInc=2zxw(7JS^^NGrP5U0ZvDhJm8&Q zgshVrelkcc!G0N~>gB`BAW&ldC2HMAd1kqwW$|SN8e~DE6^1bb5Yan!uTzSV%N4PQG<<90=gKv*1O0CD%;I|u9I%;VJUToWT zr}o_{FnWS-<;)70?&Y#Pd29u+Ypdw4#br=ax;GMjdl?ROy-4M&S%d=91#Om<1rW&F z0_|_sbIvb%ZUvXm!R;q{7k(7aK(nHri-OZM{2bNj)%#6@RKtx)(@&_tut)VH``-z8 zf8Qhd*Z435JJDBWT^#`CGm6&B>D>^%Gkvl1OA|cns(SxJw;r-y+l%we=7Y$OZuy-5 zzQT-Hb+`YqS1SwXr&#w>h;s!S6!nY_aL|$JgVc%aD$h)5F#3@*StO}_xWL3y>QxS?q|a79juL_uzB-AV7YXXYAZ!C_sKRle3F&<|kR| Le;XEm;wS$Hy-|0U literal 12080 zcmXY1c{r5c+qdt#7DXs!7fM>WC!)RzAw`Q4Q3+)&C8R7R(S|lEWJ@F><$i2ssmQ(? z`_7md%Z%+e?{&S;Kj%94%yYImbMDV_o|svi9<^m>GiJM~1(`ULN*tXWcyqAQyHL>&CBF#Oa{oGt@!w8 zW)Vtlw?Dfe!A7!5aLIh5&Oy>0%2)XAz)4!PcvF_Qz(q3s(=>EnfrnHtcHg+RiI>!M zgmdMw2tTQ1x9)IlfB?yP;M{4e5J6H}ca{3K1R>HynPvT*wZbIVuO;g5HAF~4{fCBy zghff^S1JtWvP4M-w+EV0%*04_>)J1D^%WyMtPu*_lrKiQ(acPKS|>(&eO9;X!KfI? zB+SoMSX7)O@TM+qey2F;$diSn8$9A9Q|`&tPU7Mu^Y+M>W!uC__R7{f7bVd=SAy%`Y_-!;d|#fu`O zM|JwN{$0W(`|)QQ!cjsbPqR?}XZnJq%C*Zo719Mr*_5mw+E)A|kpRZl4PktwqKyLL zmKnUH-&#($b;5W^FCwI*QjNJueW?Zc8}D+GhMc*Sd|Nn3S4kABGs5g731&p{Hd!{3 zWw^;jx9TNG%GWy-9LR*F=B>6HKh1%)yF+0@_6%g6*GScuqQa)E%;2!zNoe`z{p35_ zD5N?@?~Ib}ho^1bD~dy1@Nn>Oru=9VOx_!RK(Q!;+!y(dXMKFJJ<@iuxUK~IdE>(~ z0~>Mug;&jtqaEn&DL88TupcS+PnGFr51~TGxt%dd<4BJSJQ#I~f^kU=5vsyeocaAT zejCv7{-McGt?RRRswa%#-Z_uv_vlrQ>zSzSQ|hj;djWOQq_XL^7E$|E{`SZ%%jhn$ znwSk)!IsVPTQs%Vh~S`Y7w^xp5!*NW^xrwnPW)Nde_S_=owzA0I&fnr2LYj~L8?(4 zL}IMey-^(Fd0kOiDmCf){d95`2@h)JW9X+vrq!+Zd8GFV}41#>x^efLN!@!=mN7;U4 z6pV7aGA4Z{;ryxS-P)@tFdRDnt98dT_%BszE9Ox_(D_2uV^12`N2%`~*-D4VoGSJ= zt7c#&sONyqwHX-X2;~{*o`naHnRhE^4jPW{Ird}qJj52Bt@Qjo4`J2fGnU&JV5xC# z-Rh$ZV4S=>ed#g-D6&}vT@Fl0xP99oBZ3K6v`9mf$xMjX+j}!4p9x!}?}!AhS%ATl z_g}cCFk#?_Mb{0%1!!J6Kd$|k36^-|uzo!Q%I|VK+_8oA{@Z$!Q?U#X*to%|Thei*Q2qESQ?fV%xwhB=g)Fb~-u>J6EZG++#cg z;p~6;wGYq%f<~*Dd^9LFH1e;$PX#r12Q8ZUG~^8N$Tjs*Ks0}OV?Nsyh&bQ#Zn`rG z+IE}$1{TKPu;=%{n1B)3!YffI_k9R1Cm4~|Y#D@m>ZJ6wTl;`k?K-!Y)&n(Jn)5PS z$)Nn#xYE!1R_P@8t3@OGnOAGIv{i-qYhm*JQiCQ!lb~LWeX~Y## z0h-pQW;{Q!ar54-?U)y+dHvVdE-dB$*|a9H2Sc9omWZ$J#Y*Rd4V&EtaDGFXLrPFT zCe{~^yk9$n@=f1+kMIuT{+7?`R`-Sx7Z{uEj7Bh~zW3X5!YE2stloX$@HmF)eB_*7 z8b?aJ8FRPfB!->vw)%c>0#|41>|Z8NpgecM(+J;59Q7}%2z*FEvnv&~8>mxA`n#SX zAVb03*BR_@3Z_suXyRPXQwq9wcHVJ#Lcw7AR|gqa3f5K3OdlSe!U5|mNzJkpjJD!_ zbk2Z+?52dc@X!ZdxA2N!a$mfSj63X;zPd#KZYz^+w?QuAFrRbIHvq^}nNWA#5e*!#IN)mkhr@-;HVaib_ z3Jk4He(cUZ1%;XqWgkURVEnu<Yjr<{Ub1QOX112&=J^J zmcFmtX9!>-mWzF80Ls%}JA}#%LU4G2nh(7XRMs~bTw(74b?E~`)x6zMBq2h`pX`L; zR=zkj^LA(p91`^Z&;-t>QPko@9n^%FyLCSO2OdKqIy9?auxs|&UonF+_&W2D_3zw*&QAOHR9^EevpX67Qazk`x{=d`IGIU{tNjY?5W*k^atM@irp*oqz*mwE;BJoY3Qk+a$(TQBIdRIN4>M-J9T6<($BO0?_pUAs;w=C6jCA@idT;$D zX<9Xc!3Vr@>paKt&z>0n!SqSgUD(LpRY$?jo>I;;rqjrK<8qAG?`i!0SFWP1go>ko z-rQ77rJ>H>N`vpibkvX592HTUL92!T7Ou6=VE!d9fn#EGI2p-aR}nLZf4)z5m%7d4 zor)+Mz4Z*t9{X{uIg)|l*{0;%=Nb6OEzh&hh>84^Z=z@D40PCGRXjP)z|L(ILj#pe z+;cB1uqJ|uFQlFbRNP|Xd7lUJ0>>8cmE63OdkYh5zjr_KjAml=L({cw;tQBPcY2kN zI}^)2Exg=W_+Y7>n62S3v zW?%+y?rdT@vgYCK+>qzwNyh_1i%vWzX!v;g`(}zP4c|$fF0Kxr#?JR#+pS+vFlpzW zpFR60@zaxsdjkz7@P3v38#S3RED~-}tot;K)9D8m)Bune&-(D+-ltLWszj8S%i_o1M(`)-grsd;_iG~0dTK=#~E?1 z8q_?vk{`C$!@~I;T(@4f!;KWb({kr~!1%~=h}_#g$nn|BlVUUky7ab`QKwOOEqGV+ zx$`8%Ugd0(Po#inm)@x8@-)1^CAhi1iUz-UgazGR%mA+ckmAHS3m<(H4{-;~gKl5v zY{PK|5bM*b?l^*>U{Tg_Op|?#m9fc7j{zjzO%dd)i_AcRnNpU0taa-`$NQ;=Nu#<+pFGf zKR8Iu2jB9TGB`*%=?aRR@|+}zp>Vk~a-5_?8sDBPwsVkHN0B_d&v25o{0yEA_OaR- zTD66}n}ZbLvAEdX%0U`3DZbX1!a=&W#wfu|iG$Sn+O%_{4?8Jo$i)4rK0B%BzwSkX zjh*zuxA0|D%?b#gIeWrw`wFlHJy=pGT7tec+?(_6FT$1IlV?qzF2Gi(QuLE#fMUwn zyTG71*v7o^`m)yyoby-}$FQY=+5Yl#>qDp@d-CwQKOd(cBH)^hT>lv8ft1q4pg~yS zRM>4(+Y7VH6Wz_b+Ti3kU$lNa>wmu|?eCYeVmO!dNWT!0FzCs2a`C#q$d=!r;py3o z%D4RvYJ4Q)vA5By8+H$(GO;gnzxEhTt+9;d$f4k{$AP{=9c zX5htXVHqw=DZlOjy?Hp zahe!WqIGsSqf?A{vLL$CDPN2@ZW0!)Y9LNry6|1}{c&-^a$EnyybIz4!RQhn0&$|e zCG^V!J#k{*KsAFyS)6!hB;^&t#)|(aa(u}rPCOT#^(I*D{A@e@qt#W6xF040ZTrNC zr(^mX5+_87pU>8e4ois=@$2&q`|Lyr?YT?mCzOSWrjD8eDndd8qaY~Idr*LoZe??O z%ppMddq|1Lr}7bQn+Ix}^>_(Bwq`qpT5dvE_I0SzD=xy6T<2uZ#z{D8zxAmwWG6la zZ&cbFxq`Pcu1Wv*egR!=9~*M3GH`d^j;i7E8JuGbxMe0xqs+aJENDH7V+W|Ny@Gn{q+0S?pOGM$yt~(NkSR_&)+zD+{sD2)qOqPX2a6|Xw)pq941BISbpZX!AGCFKRstY{K?-%vMKTy6XkzpB{ z2$BQuZs`mcW6Mf2sXDG1$4-U?TC4S7#Rrnr@4G{ICNqBd?!HO1m^IjwHaU%7!-ta# zZqm{Cd-bIgm*!C7^{Cox7z5X|ym+2|eG$uSb!?^eSMY~dTdHa%8H>Eo2MQ2#M^~G#dM7}*inv5xNf02uEr&G^UKb$xUe`s0=n4>*y5ENn?GPZ8 z{BO8V_VN>2Tdftg`LX&DUY047!kxC4S8}hkb11 zA^KWwG}&|U5MSQ}hD`f&6V+2XV}-6ee$L31mo-E zms(DCVpChVrp+2QLOi_v&+WKn)H?R?@P`A7n9NUll{L&j^Cg?H#&2`zu0K$ox_K5q z@GGlI#nSMVdDeIBkZDX^`AGH?pTt+sUtHEo9mN<+?S%br1~G~+Rlho;2UWzIhf|b0 zF{skv>bkQHs5IoRQeIPr++LNR_H*T8#8Fw_UJN|ap0OfXD3CSf&&JEvpb3@@byt* zRP6*63x8I^4FMz+bi(0_N3d1{1xafzMgn?WCgxwWzb`zR$wi?TIfRL z3Ot{5nOt*k85G0le<5xelB)PM4y7&uDN9x}L3Rn^onC&j3|#`R*EhKms0+Y%{`sdh z#}?q7n(w!Tr%dq5d$Comm;o9a-ZsvS%)?&xgXVUY^Pus-a-AcKZyz4E$tszi0ond~ zJoT3jZT$`UCr4;d;J9$KK7tDQ0wT>Jr>4OvprB7~a0={=PObQCo`kJN)rPMgjKT2x ztMQ$x!yv^kCT&zQ01XFKTK2u`174YuTJ4c8h#NS5g?D`y9QLbkSbwJ(KD|E`==$~_ zENE#(zK*GdRjYSh%u6Le_>DyY$F&qNI(JUw%DYVT{`WA=frO|ODr?pJ`#1ir5mqbv z@ej9ktp40`tN}l+(^_@#MGGqHGWdQQkkRJJ<_o#bJ;?vLCP<6Xhrt{N(&Uf!qg(B2 zjq_~7=p&0c;`1Xo-C}Vsdh;0me&g@*H4J=;i zZvJajkc!J*+l4-?rlPPW*UpchsF8BC7(GnfBi2A}QT=ir$} zM~fK0scBX`lP}XZfV!FMC>W`r$Md)z0;|^--|&Na)U~aMrwPD)&9( zn!;q)+dm5|Cb50sY^c861bXg!_Tb0gG4u_u=^u6;!;L#@-bl}l;3&~hkavFw`*O>E zdHVNb-Ttp08j^je=AfY&CD?<4JT9D~SG&=qkgdmp+=WC*cGTW;?U?wVja|@#cC<>j zv>H9ph}@I|pH@X4)Ic-6Z}eSyBeX>ogcrNF zz&5gU$acYY=rp14ce~XI?2#KUvYB;*UC3pt`T{b%D0-aeRoV;fXNiLO2mO$8?N5w5 zwI96yZTHGKJO~Ey;^$y?07}E>aHL&n<#H;5|w8r7<$9?C=LdIb@cJ`vtw+Z-oalppt*aUPYJX;)%oPg<% zFFnn#Ov0U3p&csHlh8WKeKzyPB&?rqZ!W$&2|w_jL-yj~Z5rDrAZtWmCwu86 zK#G<3+Nv>lWh$gHy>=WfN6h%1Uo{S%S0n{O+sELD@4Nl{cg8?+i=f{%^HF$MH+|T9 z%?JqU1{a+@Is)0$ytM8oLqL*0*OzX>(vzGzd6jtv;j!4;a`IRo1cY38{ZhXdG~HxR zR2%ex<^if*xjGpb$AX=jM9DBlc*Nda(*e^~-_I5LH^GSpWkZR;M(8TM+%+uS2#uVg zs>_!3V5yz{rbeg|{C9cxoR=Q>k zyS-U7oVi<)T`@*QdG+r6kS$cyU#1JOD^u~3bz%}(m4$~Qf1BBlQSjYW%8hR<{p{yl zoxzz86PR5)^IY`h7@9e6vAVTw4Eat|Y}lWVpcCbQdEnO}Jo=Ft_4@Jv{=IX~{CH&_ ze#%&Qrgxl-hPu4=lNp`ZAYw$Gr8Z%^V%;H$)CL?8V@khOsztU9{yN%$<=Cu|wv+$K zH?-V)orZ57;ju>{4Vu!Q;gz;z%OOJo?1NM;7|Z{Ld!PA#e@Utblj|p_KMdR8z!B~= z=gf9k@9@R@HjA%bWe<0dYZCG= zAG{WPlLC(WDY+vZ({P^8-8d$N3Ky$ae+l122kxnN+GN`qSY?p?>d?zskY61sO0AlO z<_Hg;PqB009d__Wcq;==D=p>;3b1&r$D-?{e@sBvcjd^p2n@E;9~&79pipH08{HN` z>0ZC>pMoW*c=c4y^vW{irQO#$lDP~;8mi9Ug;pRYTtnLK_c8=k@H-AjYIB>KGK6_3Dvokr_2`oLWFeuY(`-e#o z=kJnl^&W?Voj1&>;Y09mSj~m>wjYumDSW#Od!cGTfAo+?8$9~@>3Qg%8hBDHZ?C$( z6k;`tqEv3bhndGB$W9cZ-&VI(@$0M6G_**(C9w&|ZvC)X*i6Qv3;k;y>-sR^7{<>v z4&%8E4~3FC$B=8wnJ~i?3c3wzPv(_VF)`X~{6h~NONLwPZ{KUZk^G#BU0>s2uzb%eV0>rcAhBUsj0z`9* zNvnOP0CDZne3VD30I{VgL5oRa;r$(id}xpWAs-fGulrPh(7ZOI`oLO%C_L@o7stZ; zkG65{c*^7__C&~VMu+kf@aJpiU1@$oN-Ff1YbYOaEk&qC@&PX)maYEqY#=YuXdlK@ zy23-8$vq^L^o5&vGx4CTDu9bP+w%GT`rn+y@_laVu0EEp6;T^~W{90o+-rD7>I)lD z86BTl_F)-YciQ+1Zdt^vu2**Z9x?Ec9#7Q!&{>xMbTi}AS=KpplV7h$m&Hp{w!dzx z8^`fUb607VK`fDxmphu$jpY`3)wzrRFw5kq{nyBR)Oy7fwQVhg^k7=&W}5~`>J&QT zTto&=PK{yR#X+ciL{2BO^4l6WfN_sRmHZ9{tS-)c!lS$d8yKxJ zDbj2tjtjrS5(#z^mE7Qzpv_5=seNjv(!fRP5}V#|;XDti>wAFaLv=nu^_ZyZ1A zIK5@DAVPo?yT!}DTSthbluM@X3l%1*gjTyIw}_C!tF_<#CoM)QsLsl%e=0_@^>wdT zUSavX@lCwpdE%r=_UM6tX>rov%CE0p7bQrP8s~i(T@s|F2Y0Q@_ehdXz;-V7Tau(- zffYV&o{}U*ts3sqEJ>1e+pI`o;0oS$?A3E$wSuW~rHa|c%V^0w+eV;=G3QeJR%doidHFeeJh^b?z47o>fe|vaV}n7OP?(CFZdB*R$-w(AcfB_*#pd z>scEeU%8Z*eK|UVUTJ$0ON!~pv)SEll$D2sX3v^w8q)E+`@dKr3JrJon`%h-(D2~N zN17AztoR2L{e5Rx`Ne{It?C*LO9GuOeI#hubXP9T{4f>Y$}%q2=Tq@Xsg3fNJSqk` zBo1%=$->ENPMsR)rec=BM%fLQsJPFIFqWRC;_HRoS{-UMY`@t`Qr<$tF!$UJ+aFYv zUZwq8;~^D)%Q~8#WA!H@adrnUE3cibxUIsgL&F}sI)ffA8uFC+Qp4(Lm{=8Htm#3+ zWb3OBV&v#Jv03j{QacS<-j?XgaT;!Fev|f!VCC7bK4+vq(Q$ZOJyNoSj)TLQJ80@N zSnRs0WqI`s?l{-Dd*8bmH2vd2G|tT6?yk(8K4CNX{(|F){GJ(9A@eEPUz^3q`rIp1 zC+09Ve@oel=Nu}yyGOO=&Eeha?Lsvwb2tnkb7>;;c+d13zrnvbY-JMs)z-S77<7tpZ>={2HMYB+dYQ#Q%Uu+0UJ)ib1Xjuk9cQ9SQ*CV!ORuh7Sf9ErkBQ0! zV(UzjnYhLD_`%;>m?&P&S1*2>rIQoVfpY0g%&#>M{BfR%)jtg5ukB>wM!CGBb@fc# z{Vh`QG3)zMHPJE;CmFc*%6pDy*$njF?5=T7fr*?yf5pX>Gw^Wl@@4;P3>44uxO$mY zzevBKI^V~@>nk>YWKJ@$E>NgPf`fr{#lVStPX)S7Xtc}%wv%ewC{ zk85u&`_(h%Fym8e#}y;idhtrRoOWcL>jz&SQEi*Wl)a62(+O6&XeIgTz%1@%&y3VK zI*UPSLlytHXOZu_$^uxFfmHQTuQ+me7A%|pH87w@tIp)yQfe^VrVy6cM|8L zbq`gjPhdI+*Ny|bCUCz=?8G0oaRkojxzVsObmDl+l_@uhk*)R`kNiiG{WQrVKxh=B z6jN-w(njzz{lIlqn-N^+GE_O@G>q-#FZS`;viwB5I<`{7L2NgPKJvJ45Mx@Gn^UU? za9g+hrQi(%7{4^lJUY{dcI=wV3v=dW^pM( zw-0HLKP;}K^&q$L$EvNY_#w+C&wU=f7%z6Erl63FZf+K*d3bwq%r1)0*S81b3I>H(&%_-(T#aB_sEy|d(fHB@xM7_eJ&EXLgpajVSm}$4aa-X_soUI z6{%z-yT^2$cOql$Y^jpts}2kb8xZBK?L^I$O&`*}b>r~`acu{^PLvzebtp{jM%x16 z3b)&xXzX^gDN~<}?=q|#77lfyuWD*crd1b;nys1IKGT6S`D}#=N?n++sXuOyWjCw6 z>hX>H$gJ}e%7-qH@q?SZN$e>yu2L5$na=9Ql*?q{A3MmX>1EmxB-4e;?FILYeso}` zTsrlGM;GF`Y$37_8LxcH6;3=!#&Bk4*rvzaU0AjLd!4pZ zH`1k74V^pNg(~uQ^lyoF;kq;Wu9>X)J}au}3;)@Tk9w?+we9XikMFyKWaTUF9X4LCXzr(@TikB3|M|S75AaPVxSN(Z2RvH_qiRd?>d?N4EG`SIF_pu#XqBmi} zwp9oFbepjEn#r01s}?M{-KrRn(S#G)luv<%P1wyDOg6M_zzauC1vT;1C)Fo00?YGd z<;M?XK-|l1pD+GMg(H%>?gk$+Ax5pY|IAP}xCAV{`IeLoJtAzz`M{_82_Blv@Y;OzyoM@UDmOy?Gm| zRO`T*^IG9eomx=)p=EmPM+4kGQ%Y0V+ypx}oKN0W*$5WFtv_#{Zi4dcgZH=du=x9y z?dGRsn__g&m^iuA0)EQ}jb>N3!}pY(+RvR@U~A*irjo{1P>s%|&l$D@ z^M&K;@Pu|q+_<^@v2h0+G|%tk{n!rjd4mVexwJvsn2vIYAM2h##fisGtp091NpF-sO@?fTSI?B$ z$nae)uyl1g8OEi>5_-ON!`1}TsCL%4eU+rRZ++{8o(R$9%q?UvRxw+6Fx&-UVfn|e zu5?1ZWOHk`4;d(-YEh{aGE~Is)V57_LubzTqthW>;PlO_QR-wDL;C_CUV& z?d^O9J6Xl3}&MK<(< zqOXihN_ii69SX{~*MqP4Z`iInxbg) zL13&kcm1<{5Xu*29NmV8;LV+UPfoF6p!k2Bh)y2{9^PA1=BF8 zdEriQ?=&m9lq@UzQNeh&`PipADwOLt9m~<7!NwQN8p%u=g#EdG?UE=RL?mwB5M|O} zv}EA2$0Irj-MP=W_kj)@F1gO|U!nsg1KDSOw3W13I3^PwdEa3cBtKI9ym!xmN#R6Jv*tWRd^;6?xML3DTpV-> z&UsLx-~4@orT2gOmNo0g!2rF#!jGpo8KBL6%sgC?0V~hbGI>`qz%X%JobY=FXssWz z+a<_^nAaEAf8NgohoNc{??X&@SL7hn@`C}x&EXqXcQ7EQH{x!HHWPLShm;G@8St!` zhmm-P39H)&7|n`I==if{SwofyEsFn1_4zYlY@BQBpCJaw1o;^pT+M{s3t|*~D<)`% zGp~Hv!-Q|P%?V1=3{c5&l5l3#PySq`KwB$>KvN=VeJfPo?|0zo0rA32%Z@QcmGiZ?CXny~501a9h|D_yf-6L$u%n>~_O@;ldL*t$*)1bFh z#dkAH*AkT#FzK_Rg2ZxYcxM;ue!!-p( zOOKcP0u|gY9o=5WL4&QVpUpeM@aXx z^ef>GO3WT{Cj46xInI{O1bY+9hxflOKy?2mrO}23&~Xv&o!Gnx>B$!cU3ixuBHDgD f?avbM+Lbw9d$SBYg6g6p>C5m?d_?M0`wILIIls)W diff --git a/tests/funcalign/test_srm.py b/tests/funcalign/test_srm.py index 2e7c750aa..0ac44f352 100644 --- a/tests/funcalign/test_srm.py +++ b/tests/funcalign/test_srm.py @@ -194,4 +194,3 @@ def test_det_srm(): model.fit(X) print("Test: different number of samples per subject") -test_can_instantiate() From fc595fb1d07491bc62ec62873b9fb03601679f9f Mon Sep 17 00:00:00 2001 From: Michael Anderson Date: Thu, 11 May 2017 04:31:06 -0700 Subject: [PATCH 12/13] changed random number generation --- tests/funcalign/test_srm_distributed.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/tests/funcalign/test_srm_distributed.py b/tests/funcalign/test_srm_distributed.py index 45eb36f1b..061681aef 100644 --- a/tests/funcalign/test_srm_distributed.py +++ b/tests/funcalign/test_srm_distributed.py @@ -21,6 +21,7 @@ def test_distributed_srm(): assert s, "Invalid SRM instance!" import numpy as np + np.random.seed(0) comm = MPI.COMM_WORLD rank = comm.Get_rank() nrank = comm.Get_size() @@ -48,11 +49,11 @@ def test_distributed_srm(): X = [] W = [] # DSRM: only append on rank 0 + Q, R = np.linalg.qr(np.random.random((voxels, features))) + tmp_noise = 0.1*np.random.random((voxels, samples)) if rank == 0: - np.random.seed(0) - Q, R = np.linalg.qr(np.random.random((voxels, features))) W.append(Q) - X.append(Q.dot(S) + 0.1*np.random.random((voxels, samples))) + X.append(Q.dot(S) + tmp_noise) else: W.append(None) X.append(None) @@ -71,17 +72,20 @@ def test_distributed_srm(): # DSRM: cyclic distribution of subject data, otherwise None for subject in range(1, subjects): - np.random.seed(subject) + Q, R = np.linalg.qr(np.random.random((voxels, features))) + tmp_noise = 0.1*np.random.random((voxels, samples)) if subject % nrank == rank: - Q, R = np.linalg.qr(np.random.random((voxels, features))) W.append(Q) - X.append(Q.dot(S) + 0.1*np.random.random((voxels, samples))) + X.append(Q.dot(S) + tmp_noise) else: W.append(None) X.append(None) # Check that runs with 2 subject s.fit(X) + from pathlib import Path + sr_v0_4 = np.load(Path(__file__).parent / "sr_v0_4.npy") + assert(np.allclose(sr_v0_4, s.s_)) assert len(s.w_) == subjects, "Invalid computation of SRM! (wrong # subjects in W)" for subject in range(subjects): @@ -92,7 +96,6 @@ def test_distributed_srm(): assert ortho < 1e-7, "A Wi mapping is not orthonormal in SRM." difference = np.linalg.norm(X[subject] - s.w_[subject].dot(s.s_), 'fro') datanorm = np.linalg.norm(X[subject], 'fro') - print('Difference/datanorm: ' + str(difference/datanorm)) assert difference/datanorm < 1.0, "Model seems incorrectly computed." assert s.s_.shape[0] == features, "Invalid computation of SRM! (wrong # features in S)" @@ -134,3 +137,4 @@ def test_distributed_srm(): if rank == 0: print("Test: different number of samples per subject") +test_distributed_srm() From 1ee3e5b64fa2b2400e4ee7a94b684a2a16ff470b Mon Sep 17 00:00:00 2001 From: Michael Anderson Date: Thu, 11 May 2017 04:47:39 -0700 Subject: [PATCH 13/13] switched to npz --- tests/funcalign/sr_v0_4.npy | Bin 12080 -> 0 bytes tests/funcalign/sr_v0_4.npz | Bin 0 -> 11809 bytes tests/funcalign/test_srm.py | 2 +- tests/funcalign/test_srm_distributed.py | 2 +- 4 files changed, 2 insertions(+), 2 deletions(-) delete mode 100644 tests/funcalign/sr_v0_4.npy create mode 100644 tests/funcalign/sr_v0_4.npz diff --git a/tests/funcalign/sr_v0_4.npy b/tests/funcalign/sr_v0_4.npy deleted file mode 100644 index c9426d92ab953bdad7dc8794a6196c125c875f7b..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 12080 zcmXY1cQ{w?`?q=9JE2b~MMNZR_fgW45;7vALYhQHLWqWx?24ok6-ideee9jR_uk|E zw)fZXx_;-M`#Sd@=Umr$p67M$*9g$oKXcBAo=%hQ;r1(6?_Rdwen4XTA*+4c}VkYoL8$5kX}T+N5|JU|_UBJ>63$Xv4xjemB!aXRJtkgr z5(f9!jeq>+B-)3?zg%48BsMG6C+v;pB!p(%4lt ztwS6{p1Hu~*a>!m*|unvonR+4e@o1Xe4!8}0gof1c2NjQ(vDIzVIwY{D1Gmx&q`z* zqFmp9jfJpj;plmJkD1tYYx8R6V{8E6!SUD>X`R$b808|K z{Z1K$_qW>5=V=VUnTz$&V)gCdBtEGzrQQgqk0|Eu?97MC)@u5E@lZ(WcDZ<&Ascqe$YHzkAVW=h z@?zf6@QLdTWXy$A@{h+D$m3jR*p?d@NxI7Dn)W+PBsZV*k0(h?G_aD(8?dj&R9?rOrP zZI4;W>|V%&a$ zbEe)dKFUm{|2DdhM8<*2GBiy|D9<`Nc?%bKOzcSj0dwkN#zBo~0*eA}&?D zJVZ|h_{a6N`q7ceRo8cvZdyY(`iY&2@+(;P@p>W0{1T>|jbtytMcimOIo@?|9&0%E z^0#o$p_-!bGdJxS)O6d|^!+ytmFVy040cdaa#&Eh-)IsKJ}LCgfD!cexK}scKZNBG zUw<)p_Ttp6rrD;IE)-Ag?){L^j+avd-hX+~fX|FE#ZtEd>l-aZg_%fv^YqA(S-6IU zmm1DXd`f`rX+*KKK@kLfbWD61TnSDDeL-?X9cb+0wTTRAffKSKU$2FCK)OSro0okz z?74UT0WWti&`B3gZn5l#gGbzlf1MqMHy43&CutP^8i|aSDvv|V-t!YljuYV2d_SH_ zM+Lsb{mz>nQb8)8)w)`f1{%*p|8TlZfs4Is1MzSQ>O;ZvvdA>NY6z(dxH1j1zoRx? zx19zDg;6^Db2AWdl;d`F+6+8m|B=q|W(HW5U!3;KpMjiY73#U<8MqN=aWlwm20~}Q z-;Lv%fjyfFzkj(l15^4NM$Vm?0r`~$cdycENY0-w7LA#Pq#%=qvUk&Po9jdFXxJ2p zovLm>nmPqKEjFYT!xZ=%S2(B?QQ^SLn}LHnsBm8A&qk`lBs{RF`DlH40xmKKJu{3P zgJ=9b6JPU2!F<_*&-dUk#I_yEG5c=_Di?Un9*quy8|$Nu?YsK{X=gXO9_oeh`*c!X zkzJ5PG=$&S)B)Db)e+V)&2UFH;olp*dQiDwndqH{qHEpsf{nIDh;23KEf=t9qt-p4;vghN^o68Pv2_?uRv2U|lbvfSt zFXqWYUL~3w+k5yrPb2a>jZ(rZ+E8S)zBx&z6H^}^)VBWCiRqfXiqS$nSh!PAODees z%c#O{XG8i>!gRUEF`*xYPud-5XdA>Uf4NEJ{lmD~Ote(^D9_9PAj3bZOt*REzNnD!yPo+4061O}vF^_sTiQ^G^0ZMUH6x@`RTA)kC{rkG2 zRJdt4H0)P+*pP-r-R%;zTP69L&H~ZwM67=Xh^85mSl#}aC4_M zN6L2^wy>K!4A)R`@k_jo)Cv{bM&*e`DH;ksk70Asox^=O>Be)YAW%mzz2&|k z)}my-NUzKF8g%}k867iOiB3V@?ELQv@K-}z<5n|7idlQd+GH|{?D$c8{Z#@|Qg=J~ zS9zhJ(e7hGvl$@1Nmys3m;~u#Tkg-EL6DybDfOEzhk|PxWtsac;p27X)hv-(NdDMJ z|9ZS0cA4LN8D7&2(zZ&*T?bk~_==Eqt7r%4owwXN{iqw3J{1XV^6P@~y;uAb(|Uk< z|8SMgy*`L!*qkQn(+B4&yi07l`k^xC?@){T0BDFEW0RQ~hUngy1iPRSFuJJp_nGww zEcO-HpQsrD9oIDDm-eGD?mux*T3{UX?4IPO7>$D6w`3{7o(XX0^?CZ3Iu3Gen*<(8 zPJ-JtQ97|(lb{f0AJ=|j0*)RFqN8?ALMh}kwI)u2$S0vsoN`onA|xVqY8w^A21er> z7bn0mesS<}Y|n9Tg#J^>82C2}LLC&8o7=gYIS3AoS}XG5`{fXaKKY{I)I;IYvg zm716_P<%Od{P59nP!-SNmp31SfVdJDp|KIzU3LDbLE#W+otl0ne|i{B{r32Luyg?W z{+uPvT!ujRi@7gvZ$DH?UFp6T(hIeR6b}E|-wXdqwz}Oc>IUui1|?>+F7O%QXlZ6@ zgHm7p6%*4|p!@q|8~kjBOo5HKL#7`7a(i{WcvlJ0pVTv*CX3h`?n+mi7AJL%bL)9+iAQ9Yu`0HhJ91#Ga9yt8RwH$jjAx z>E~)Csu$}G424vq0Q08;Cb33LYSnR5&}lznyY|-)<3Fydw_2x0P`yaJHZXn^V|BF`LT5%% zwlQZ_IBy(R)+8?_DooYh`udT=*;2Tc=Y%b+P|ip=ir^j7dnHVHW^LhNx=y%^V`!X)i%wm z_;eagzs7ygK0J-N@^o2+!qX_{s&vMibs8^k`*Uia=@j-n{CSdG#}A7xgK_sk8iv0* zt~hN#!z7bZ6VGKTZk1$`e6B#l=bOwIuIyRI3AG^=9mPr9wBdF}*xgBN>S8E#IW~@U zk^c4ELZf(LQqWa?-!S@L`Pf=^dk~evDvj=EOY!5-b(M#t_pU^X;QD*xRDYhI#nVsrc8uhhL~ zb_N503=-RzSBHTz>3Y3Mdkhwb>PJt@O~R3aZ%(?$X|NpocgKFqDG1-h3Jq#Aa1s|H zZ#T`tys9|iraKSrYTShr8y2DZu-NW&#U-dH_I(qI%W%-{Yd&Yj3M2_|1xWo}gUjAO zp2X?V5w9Lz(;h9PBNn@Y%EHa*2@#&J@4r^l6QNA40z3H`2-f~Xg=qqegqw`acLNF| z!KPp@c$jn`B#*B$*of9h^&0!)U<5fP%@h}s4Ru^}- zXfYF7oBM#807>p#VHo^suJ&V0g~ zfhbpto?n@vC*;+)y0UoF6A3-+91)-B2xZnEN#|78AS#7xoZfmB4rFQ1Kdf5@d;VvY z%BG7TN$&wl3JVYwdhW9<>m2B3ruY1NFats+AJmu*PQ%#EO=rXHsPHdPRh9a10^BB# z{+aU~f*(txof0FxVDi~~>X2qP{M_ec;_#>$DkJV}zbIM_7v{=dOgm)54X3rW*Rz3` z+!vXD)wURGJ3NO?>}xSN+%W88XB%=1?A7o*y?#zxMUq+H3}aryu04h-o!t@g`Ry7`bspEF z^QR*-*@8v`1L(n&bEFOD&z;O@3dsgM8S09L+ zu8T2f83lExnq4Jlsi5*#;9}IS893%?X&m%;0X}~bH@Ez=0+|{&Y)K;HEVL z@q&j!Va(+tSbiLBr_bUej))`1w^3XjtWGB-3vV)50Da4*yv6V4f zHX_bsUvZN#E5V&UdgVw63sIr_Ex+EEna~Vbwf}vGiFnliFmcs^kzjY;HQ|uLKumYl zju_sgCuqy>%-SB)5o|2YQCt44g3}n^)%%XiaJyl@wb1G!JZn!J3YeLPdT;MbsTy<8 z!R0sdQF;cfO?GeU=A8lwN*t{qlnQ#K$9O}YPJnY0>vTf(2&{2UCQXbCKwieh!o|5B zpyRt(bwH&X-fN5RvE14Q)9TObuE*Ab*X(PBEehozZ&P?L#2^p2%HOF<7DPj2kK>Wt z=RxQn&mc1un1M#kOL~J2g;>IFnIbh=j_+SN7)$HdVA?X##?(`f$uT=`_TFwrlhyR6 zn&%WhWcuu2bcAm zaeH9x2JXuZm@1`v{=#YvUZh{WwOF$gHKKSIV zbpLB~Kj^p$w?2aIL7HERA%b^n$UxR$G3f3VOnET#4`#=DHV@I)LipJW_M`gEV7&K& z{^5x>ctZZRnS9#`J+(tixA%3!)?cb;G+g`OYwF-tvHX6h*7W+USv?3JU0Qc5>x@D+ zbC9b8(>Sc2eJrdTH3|GbALVc{P@(nnJI$L0G^ncHD{b?E1{LyCLZL-d;O@OzEVF)& z!YAC^gu2)Hq5G+)UXHUcUL}(Ft7#6#3*$qg7Un?oLZGh5?ggNIt|-10xd0!w3FV$? zUxfMg=+>7XmcZ9sFQ!Os88SGlGUIM6gXQtCMOxr8)ZFO05q)X}OfJxKe3@JUIfG9e z9OEl6H1VEC)Nd8oPSj3%M6AN}z&@jYuIu%vQVILQD#S~7Fx~1~g#bNS>%i1CxE=d% z>%Q}A@bKiGjL=7Gu=Zmq@j&Pr_<5|ou1j2lQ2l2|b_lKkYvS6N#;a8rx70uV-E0-= z+BrXIzgYpBMc0R?AFsfn^18?s_Z2wY?Kjc z=`3jN{b985b_Qzq2yySVpN2ymrGlPX)8MzTYQ4jC3WltyWA7ztFmS5pz&AxI<_-4ni~aKQ_TRgks&DADcq@+KMb43#2Q`d2jD-u zr65tgKJX=z=bx?5(T1LrO%6}HU^(^28NKFCcthJJ;Lh3s(j9aGVM8rYx~KQni3cs9 zCNkHsYqkz{{Q9{x?p6a{uj|W8eip;Vr5zn1b0v`KRWE*IH4i>a6|gUSNr4A8$F6k^ zq=5bMIp%x$5ilKfoqo|V3>MysB;T9=iJu!!WL6)H!MvC~S)!3S*!9p?hEpmV5B$=- z+-RALJG+cE>+a=ZlKeKJd{+@NX9*qFn<_+>KmplBr3$2$kk6$BD$$Zff&103Doo$N z|GrkH2BXHRJ_H2UBJpGgz4P%JEXYV93+8GuEH=6|?ROo5#;g-a)T7@DvRLxg;by_{ z`zfY1_{c>-V}M-mw<^Z;-Kj&7+T+W19QAm+bofe+K`nlCxx+*MtPYL(uW#YLQHy_@ znm6oHt;QW86K@X9SFhs?lZ1nDHM$E*zBkxdgBf2$a~;lB;osGQ){ph&IF{x9u=@8u z+%{F9n6y-k36h`FMb=92iQMNGnGr?Ed&^nQtpU+L)g#h(F&Y26WbyEfcN#uXe)PuvHBQCV!?SEc6`(0i47SUCL`^onK2 zjSc>Rn{WO*!TC7>vTZ{3sbwUnKTbTYYFGebZ>Pg#-seH!+wkj?C-UK4kik)MUnx{Q zjbpPqUJB`&#$01|Rd7b~+vL~qT6mn38F=_wJq&$qV(JxZgdj%8lm3B?V9W71D<-oM zQXWkm+QqMfiVFFE|l zYJ)BQ`f1x8+F^kE?CTuo4xpMUW*g>p!dVu&$QNdvaQmc#+P&MIpf$C-52`xhNT;qq zxl%WzsChh{QtE=|59({^_`6_Bg;K|JpDrlo-ORga*98|#^^Y>NcZ0d%_lJGlT`*l7 z6sO$W1&8Qz&xfn8Kj**g>MhJ2a3_avp!#hGh+TMDNh$AyC1xA{jI-?!R<)Y^PN)^~ z3K)!-x3)vO9_Yu?w}Vc#HNSsID_nK0v){U<8JgcKylRbags4p>CK+Y*@Slz@2e(Tt z_`k|LAJ9nFpsiQXaES z<-$SMc7A*39Eew=@W-6Zf}tS!QKbt%A?@;_|EAd}xYs2!#rHlEf-mp<#vdJlGEev0 z4eyRYN6mv#c{ftgZncG~kFp&}+`NM(i?n+I}~)8BeAj^KN+9hR0%Oq%Q__;-PUh3FppE zEPYZ}?|i)*BWZ7Rw_5e$7Tp_S1#kOsp4ou=;eBqNj zX%ooRzpU@zI*F~mUjmz)Cvn-r!@@pm68X>jF#pJ*;=a8~e%6Up>}R=Y*$_`kTWmHiU+f3dy3Azi7DWLdF+rEDif=?DrLx)3DrvU%$>^3K!mP zrJmV5gW7u4p?AkKS+(r~ACFgY_CeqhF-q-|tNKAEi$rD>j9MeNScl1UvQb5Nd z*>)786nK)j8iz5%Y2IL`&Jad21ur=^4B@EoK@*nB0n|k@q{Evx&e1|J~+K!rv^JzT0I|ERiebxL-p@V%CLfkv(ST#6o^M~>m z=WPfvp2}p5-25x&C*K=9`gM2WaI&jP*tXjD**Y@ zXCC;^ z*Cf)h*&Br5x(>vn{NK%eQ&{R0%zOXWG=8wi81UOMgAwU{`;OGiB3IW(3G2BzRQ@jg zB1L}z7mA2ZM$RQ<`?vM`(#B;RXbsNvj9$TA?;7H7POf53Yt>f9g*A+pHl|IHbmXmA z4z{{S^yDdhPx%%#2GTsE`lPEb1Nm^gNYhw`k^KC|S(T5*NRI4}>CZaIM83#vz5DSs z6S>XQhNh6rO#1CSVlK1pSE(rszlE*)Nx8pviF6%dC66%>zXzVMlDjQy z3bxO(lG^5@OwzB|NUJE?^DiVD`CoeNXYgVpSDp%fDUD_$S&vP=I4(^gg{z&E9V{qh z5`+KOP!|f>t98e!(w;&_wy#96x>LxX;jd>>Z&FCnE8|-G11Th@yFBT>kwTtjyvEFD zN+Ij{N7$%h6!Oz=HSdxNHuCn-i609-Y-HHO4c!vQ*vL&M_6hUUvy%7A!oy@vtov*B z{XFWqEabW@D|Duqh5YwbaiW%)h0HwFKYTZYnM{1Mv$b&#GkJ=^^3n#Mb$)*++wqkr zBe|W}(eF=TBxB?Z;xC&qkgoYl8~(=7lWChH$Ez~w$V$x?eMxmX(rGpOsFKPmQvdR7 zwBEOZgux4YueXbswk1BY_{##`R{kz;*gJHpWp*8B!nr^6Phkj8K68ja4SOaQ zI0J1a(QHd+%-H8KOrxew+>jl{ogY~|;eVSVB#n9Sy)3%_s3ep1WyK?!O z;rP>QkKIJu;ewBXOucS57#Es)sh;hF>=G&q|Ly^pyvZHdBRB-vzrwV3IFCU3zY5yO z^)YC2RP}OL*O5g=#8iX6O#*%Q>x}oesPNOV%hPn63S{;~bA$Fe|8Oa&K5~l&Z`SO+ zx7}kVe>{1g`*!`i+MZ*zy|SH=lshGx5hlz)`ro6y(yXK-b*C?B+-Y0EGe4XShF&kD zmd?GN#DoP*&w=^LeX|%)ZS?hMAPrT=ZpWrrPU2;Q6A?xoL&!yG9NkFk#hUA%Z5iLT zp{ay?h|i@;{IcJM^v2}0GHgwXUtlcA^Yb8ZPznv;8yXx z#%cpS5p1`%Yt4~?usQvcFMcy4VQC$vcyOK9zxA}+967>FM8%wmKVHE?Oivy6PMBpS z(!J*`1)i`G>$0j%OB99JC42HzVg@@Q+A7fys>?x$BtNQIY2hG>oER>KzU3rlcZ68& z80REBn#J!Xedi+b+BAg)Xk0|0{FhJ;3vS|oS7h9L1~i#)^`@$#LF&v*#q@Q=2C`nU=4xgb~HO+3WjR^QdT zRBnRu&PtelCpU3b`&_`a_4g#2nz-?Va1jApl`gCEa1nC1ouW3K;3NcseJ&MoauN*- zo)#}u*K-P#{f!RM>_n*amEOy}6r$hToXwt#LO9C4yep~AMtHRNPc8?t5}Ay<5Aj`L zA?l(V$2L-#iR1dL8aLPVN<1g^gK>$G*tq2j^Ennq;@nZm|Msuzk^M4yvSrzHMESpj zJsC}F;I71${EuY?QWfjnx;dAi-g0F|XL=s`Ty4zymS!R6Vx`02f9rf#EobJG^b`n5 zz46mNLF=xILD>B`ayXyT2dmTY$K`MZf-{#i@m%Zd}tnt+V z?{T)@n^+6LDm$Q3Iwutlr`<{uQHw_5#};bu{uH2raBzK-X9aFh3+LUvw;oeVXGZ9G zTd{wA#<_E;13x*Wwo7vL;o+4i-4ynI^k&z7-fTUH1}&C(@zNs*{39ko??#cYd5U!J z9mhWDh+|RqlX&csQbeR&?LF>SjKL{f%}S z9;erH3%9o~b&}IKXwqPn)-{ciL81rZd}i=H#|hVE`56p2k?ouCb_VI&&07qzrxCLC zZycMSMk-Z$B28-sNAv5S3yV(UKOI;9EdFV{XL0bMe9IJWdTh5Ann%Nzo}mYo{?hQO zQHA({{`|s!(33magSv30g{}aaU;Mpz;m9HzpaMuww1^;09GM7|-edF8((KF0BYl z3Nz-yfum$-d-p$ZI&?rFP`4Hk7<)~P8{xOtW^N1JHt;Ov)0GnHg6d>(vGVwCm^IjL zAsgKb52HnLF9{BU(cy0&xCaNpr1o&C&$@rtoMJLE(>Ds+9{Ki#?i+{Md-W9J_#~8X z^H$-Op~5-B*>;x<4VcUqdmde*fhpzfQ6Hgci0aT6)ia!ieRan-FCLnPY&EXO4|mJ} zo_Q1#7cdL=7K4tLAD)Gf80Vj!wzE(slpweJ&n#rPpMQREU5_c2`BX)8%|VoB*6#|+ zJp6cS?Hm_62SE`LM+_Y1A7pSY&q!~A-FULj)DO?*SzLoo99hJ@x?qe39o8< z+08>_er!z`-8`JzbKU;-k$LE^OlZpaI0rTlWzH>`&4F-iq-BoD9Grh6w(#uQEYOy( zl7j6s5WD5IX-Ult2uHGIe%77=y?;Bs1uxG))ARc&pYx`mVW&KCEollk9{7s7EYg5# zY*F)>LW9Ll1#c0pHc^pn!9-C!%xfB5hI z4!AMk)OfkS63E8Uk$ln&f{FI z%aRI%tip$0O8$Tsg2K`#oui(J=&#=`&7WC|9A6s)v@bSdO5+QW;q~0fzax3nl|@yT^?e?H_n6@kuJzo3$ptIQ=^ngGyU2Rxa}RQMX7Ged z_u+}5ISIDEefYHHPxOVq{rE~{i`#vr0X#D8@L$1?A>8^sMQCzx0EdmEJg3G6@TS%V z2QI$>bUjmEEtfQi$+6mzohJv;-T8qEbMp{>=Cw(btsTUeV=E_4F%RQNtl{Ms9)ma= z$ai|gb`V$lZ{K|McK|S^~Wob4&KZ!NeYn{bXQ-c3_q3^Rcz1b~J66IbB}RhT5NR6a|&Fpxy(w^Jk^%(Pq&}e@93I>fH4x zV{NEKEirlRymyr-;H_@a$5VmrncuGZ=oFyCrTBOL*K<%-Vv9yz(I0HkRNv^bISKs~ zoVgzB#9{Sqjqc>SaPVo;H{AR%9ljQAzec&41slE>>4_Xo2gN;V8b&vA;q_mQsg!RO zP~;Tz!~a($sPg}<^HHmXbL576eX#~e`SI1)zM~0}ujdtM1vEmEd!F$#r4C4B-qP{a zwgY7Q8!OM6cfpQ?|2#CMy21B2ILdK#Lr3yzwuwv+_!dQm8Lmqrg>luQ!`b~{d&+xh zM57-@K5cVO)foWSL#AhSy9Xe?uh3xj)F9A{`h7Lc9Ryw0bIF|-haiY|w`IhhVR(N< zR=8ng1hkGAOXTnk!s)yl`{Qg#bmk}@-jIHL#9D*=&=Q$z&AvmmHaPaPfL8whoDQd_X z1cS?bNIn{bTjJ8=YdZaK|C2f8eN;dA>vi5Oor3WV>vx+1j z!c?a+ERo-PH(9U>S=Cto7(*RyYkPLIzo`j@zY5XyPPF4SKFZ^7Bb_L-`>MgW5{ky>~M%4N11Kc?hL9-;x}8){OgT9)yqXU!tL-P0=5({m`T(q!ds zHIMygDdyjz7jQ{`wN6ug5r-|s*#|lokwLf3G+THXedn8>3kohHr(^T=pYF?;BOfe% z=GF>+GU9VOzqW#&^BlwR=T`AEvvTO$1z zPs-Kk$gx27WaX1|K{G8rq4h} zZYV4KwoQ(Xbh&#k@UizA9t)m2+BOwDLA1$Hrx3@jNVU z?Y)fLvfQ$P%*(iVK1%-1mqomDs_~W<<06LdKQ-GDvVe*oxhr?Pnn&wDx~1FuXHnWp zaN~B9S#&OJ$`|LG!I<~opKbmzg%0AI>2uffP{F6yN@BYwu%)Op@ub}ZN*dhVcz$jK zpX|BvR_^B@hV|XFs9wh(NdZR7=Exr0&_A(nb{%KA**RG=Z`NZ-SBa>`QWdh(kNw~v z%Q07V>gGeGZ}`;gZeO}sF!X=eS?JT81WRkHx4xH@fuc{#m)Lvtu#aEn%-Fg0Jj$}} zrnlR>z|?Nxo6p-`7@RM#99`=NWm`#AN9keE?!6F1Y+TP3>=f52Sm!QaS-8G2|F7GTInc=2zxw(7JS^^NGrP5U0ZvDhJm8&Q zgshVrelkcc!G0N~>gB`BAW&ldC2HMAd1kqwW$|SN8e~DE6^1bb5Yan!uTzSV%N4PQG<<90=gKv*1O0CD%;I|u9I%;VJUToWT zr}o_{FnWS-<;)70?&Y#Pd29u+Ypdw4#br=ax;GMjdl?ROy-4M&S%d=91#Om<1rW&F z0_|_sbIvb%ZUvXm!R;q{7k(7aK(nHri-OZM{2bNj)%#6@RKtx)(@&_tut)VH``-z8 zf8Qhd*Z435JJDBWT^#`CGm6&B>D>^%Gkvl1OA|cns(SxJw;r-y+l%we=7Y$OZuy-5 zzQT-Hb+`YqS1SwXr&#w>h;s!S6!nY_aL|$JgVc%aD$h)5F#3@*StO}_xWL3y>QxS?q|a79juL_uzB-AV7YXXYAZ!C_sKRle3F&<|kR| Le;XEm;wS$Hy-|0U diff --git a/tests/funcalign/sr_v0_4.npz b/tests/funcalign/sr_v0_4.npz new file mode 100644 index 0000000000000000000000000000000000000000..6990802212ca4cc5681a304ca5e493c3b614e7c6 GIT binary patch literal 11809 zcmV++F5b~lO9KQH000080C**f00~w_Uw$?`Us%%lx+e zDT?&%BckYx^kbEfl|AYO=o4E?S4cY)r7=sQAM~ zX!bokH@csNaL(GQ3EZs2!F?*Cc64lnj*x``JF*dmc4ww6`zgdL-e0^=ZP|&-dp2m+ zoaP{61Ugc;z2_j*^g^^flsJj#FbVs!$2kcX>CFdEdU6sWx=NmtFFA?xciBz8{o*7# zMkGF6SmGo$DK{kUiRC1QD_Ei@BRL5J8AVcRg@bq_hI2u&97OW?6*ip%97KVIz@_*} zc7oZqWR0C*C$xV_&Wn7a5T!wnqGEPZ2ukwyGBjl)E*vX+>utzNWFDYg-+PUPuxaJ! zeSVjj*m-NyTGk^bV#?bv;HoYoafI`5P`oVzp?X-%YQT}6xL$err%Kg2ytt?7Ch}$# zECwv*Zr)u2zbN|?>7@&hZp9ZRd3F|5f|l*lPEEtn_*5Cav?&Z^}S5)Ku zx>x6sbv6I9Q~eT-3O%bAp{}5s+2h{{f7fuhE4e2A79BZaR}o^QO;27ZI2<{7oq>!y zcS7OOI0Jc<>lE8cBO^&y6&&lO(>R+5{uF~Mb=mGoMcH#w}zMjm=HTlAKVjkFZw zHVD>ZBg5$2yb_PGk%|UFwLF4s|GXS zB~MOV+7Yvbm9!Xb(~eVSB{NHpH{4HPA#=UU75@`sA!TWELVdo>q`C|B_MgMdWX3P! z+sc=j$SscV6lu+jB%hz1uk24oa)RrQqG$;N$vl?L*fK{?&PH9Ve13qQ3<~_--xfee zrc__wUbb-^-RUQHC@HLB{rl_191F{sdODiD2$xXQXlkPS?gG|w?BQ?ao<|KO;ivAp zv#9N^-u(3!4VCE~@`gI8C^aG|GhjT0`yUtkXTm7@c;2m_7#PNisLwwcy!vo@PTPFr zYBx%x_4K_5Cc0z_@vAef@5A42s_C7CnAJEAZ zPi?jufc=L&Mt+_ifmauRawmBVvW-Q?%2XyGZqM1tWT#1RZn>90rK190(q5O152zqr z$ZB1qO#`iG;lDZEr@__Ut&w;z4GrPobxC9fUNnZ)2VI_lxnD6GuiMUmqv9By{h3(^ zI?QpqCVdtjvVY6qcr^>GDt;#e3TGiNMU{FcWfpFHx4aqRJ`3SOdaX9( zRfcH@G^un{Euq4`=Qo3gwo~D(-fvN=;}qPttbK2NX%a3lhdeck9*3v=y_25{#=v65 zlFxtt2*kG^$TRFH#mjc`;kd_vF+zhH}*pOadXAdJ+Ty-B-J-*)KuW@|Kc7m7F40>kv#{m z^E4s9^B5(vvK>Xn8d{QNyD;s+eqHM?U6`TWrxYvHi^V$xb)-{zv79RWdM>OVCCyfP zoe~F7__*D^#`Yn+oXt(D>>a^P=3;gFWuy4MM)Z!K-Z*uFJ}mTYoj@M( zTh*rkptq|fzwExNqZ zj*XkDLgx@4cK$a-__OhQ(-w0?ig`!p`cw*vZ2wkw{Y4^D(snrqR(qqM@vb8wbD1Eq zQCM&E4+%0yHs6~&g`hASRu(W<0Y%qD<(LPm;Qe)#wH%Q;NO|8x|8k-Mc3Rwg9$DK0 zGPcSl-TPWW__C07n^-3poVD6A^RNe&Ka>b<4Csc6J(mNM(tCk=?_jmw-F}E>*px2j z+Ye_deM)V*2cRl1d$`qO5VXXPu*uGjKx|)JqFu-+7++A%eri1mOZ`Rm$7)AG&n?~L zx&0VS1WxXk5tsl2yT^s8#$#ajB}H1WcM?2!eV;s{PJn#-Mu7)XQ{aA0j86R46evd6 zfA2Up35Sn_&{4ampbQF`+LER~rke93P;*Yf|)cj*GI2=}k z0?Z$Zn8cefxlPYqQLh~-o!=6o>^d;uxNhPXn=bTf3uip{uU_fohbF2|^yBfr&)Z(k z^`n~lG8(B5qrkV{JNGt>;9sul*E%OgQL{v%E;wNfY9xwCtkAOheP?OYL%#RE$3oxa*+<4c{LtqeOd6V{w>p zRR88_bm4GoI(&2*?O#&Pa`4WepWe`ijm9&0TyRpy;`R(mx6kk@J()qX&)?ta9-P5^ z1-hJK;Te>7Q$FRxI)j(C{yw2@HjTXxejF$N@x!wF{Dj944I^J1Rhl_Z!(`JkQ?C^& zZjoY=dZtLjXB#aRFYo@x35{V@J*6q!xZ!qYgu@gzcQX{b9+|-U=)eYUp)ou+CFrK0 zK7xUl-?x?D9zx~FD#@q0{TN^+=Pq}y52gJz9+LyzC?q$oX}G5y9~evhOzmw&_kz(M zq;Ne(oa<+`u`b8pN7o0^SBvm_jCxtquN0gOp);4O3dMK65&J^4Nm$t|UFw`(2ni7$ z{a=2T1E<4-F}Ip}*sWjJFP$!`MmrDL7R0#aaI-4OZf_x9_!@hRBVq(5NvB$8jCk&x?{z3 z#8P)id87qBA;R<+dY-~auqoOL`6461>;6HD z?h+HRcW$9;={yrLZ0?~kYtBTp&5M_fM#>JnjGpX` zWg={t*&yl&6Y==Y$AUv|83`46c1zAMM#A-v&VNe(JmtE}yv3vs15u$7yRbS-Pbg?^ zabxkJClY(vIifz$5h|?TlFz8ELrf~y1iked?90(zcu>Cr_WVz)RLqt@iry2H6&E2U z{LDu=)_E|@%IN)je-?yH-)S)IpMmk48&5~tQQ>c}ni}=NB)Cr<{ypzM4BwW=x+F*Y z!1SZV^a1T2_@VA>>iDn)s-o^}yC7Br=jO}(W*l?jhV%OR%ei1o>5neFV*3Z{I=x0r z?dvc!(kSA6S37bH?$Po(`R|>!iKMW;8o`3box6=xCs3FCJMSV)#jEU>*|NH((fG{q z%@HB9$he~+T*z(#&UH+#l5VS*e0qD>$JgsP-F4J}E|8AQVhb4y4x%R; zj1$y{N*Txn8T;nFZj7Y0EECi4TPCva^46`d6q!l0RUM`q5-j9pmQMe!?<}PFHd6o5 zPge5vpkE?CFB|#!SY0Z{vynj!%TD%XY^0wHzs3_Q3aOD8C~@fflhofH`p4YvOm)0McLHV;rp58rH=R(T33 zVHsqjd!0h=;1QWBonRx6JZF0-aE*;DcD|>w=`|~PDT4BJ`v431t9S3TDmM$+5^clz zV+S*t<+^Ra{wWhFM-5f}WXDMQt}u+M%+Qm*!;}b333^h>J2|UiWgUYK*l4^@S;gVv z&rcg!mvP7URZhe3dGy$nmV99M3>tnDigTu;;smew(b&i1=n}@RonSJ6^rH$>!d~rY zb@J1%Xx>UZ=4O5SXljDg)FTbwxlr~q2TLZbkmxF@Z+IS81wlEmT!kU z=yUjpLn2ATYft%z%gx^_ZU^xZd%`r-+&uXR`!mdZo>F`SqvR`^>I5&bCs^*%SS>Gc zGGWEx=4D>u$`C*KU7DA8k$!>StdED#d*@eFEW$%b*NCPY9^fYGotT+#BykZZi7-iB z3KtQ7{Ci5kDhII;%A)-(i=Cj7JamtY*oh3j+>k%?6k>Os`0BVV8}Z##{ZF$nE5V&H zcKJ{k3sI^6rLe)Dna~bfv;TF6iFi2hAZg8!kzjY(Iq8_nK+JU4jT+sgCul2g%-bK) z5o|0iF`NIcf%7=um3vMraJzA@wb0rUJncvt4w_wn1|Oe`XI*>}4IP1iD>Y?A84oM9P*9^OuT z@DFgaa@}YvBYRO!trvxE{NM>r=|kJ;13z?zd-47kJyEmFM68UBNknS9;JwZo|cV zzaGrl!}6F{zZ+*W1{UItyYLFjRn8~$9e6w^TX4v?9St?Y_pcbX;I`np4cwO+F-=a6cKR>10L#1bud1!`|}--tc23a01iIX+FTg8cap92^s?Fg*E|M=W3s z*pAgrc}A_l%%HmQU$=kpuyQH;;u<8#bTZxQUxOe6IqTrGb+{e>cZ>Sjb$D=mcV_s* zby)wloU|`|9RfVpU)CqBL%8A7L)!({fi-D;TdoE>Bt8T<2W`q3C4x;ji;xKI40R@6&lW-+wOjnb{oZ?D=NA_<9!V zb_;Ruv7dnh9A$!DIx`TkxMscGZ5oEHspD^@XfSx9ci$H!D&!q`e}>XE1syta!P=G+ zpw-*2Dc>~#QtbCXpO_y5IWz4b^U+}_-67niQ8)q{$Hkjm8wTM&yX6ougMRQQQx=~7 z>!XdlrkWj}bi+#8w^If!UGR#wRltL_6J$E+f+B`np=@{Gtz-9FK|^G|apzn;Z2$RV zdBVLGyk9m{l>YbwqRZPm!{$pN&AUP3&{_e!n=WEs{FDm!ZH`>)9!v%Ml{3tD3!`8r z<~se7Qv@u&7D>6g@B=?K9m}fOABP2TyK}^%^RWAYzbvP8F7ErOf2qkTA9r+{XxHD( z$7F@AM8(b$WX=&fXfR!jEWrYDOUji_NKhfTgZ zuu$`lGfa|>CN=0GDE0QdXf0-b63cfyU5(jmMQ!gJDsVi<<3Y`@zqoa}NGW;w4<<@| z%n(^G#mDj={j#D;koT60yn7>JV47#N|56J6chT~}DW7zF{8UXep)e5(EP`s2tdg)K z#4L$!EeUmX*SBr22}BK-$`9Q8LQzFr{j>5jOX#~oJt&;<6Z*t+zmE_7hMTYcJI47j z5pr$94XNcMXg*3hsb*9J;;&~SWZxD*@axFyQ^yM7P00Dfqiv~x&Sx?q{vCNT4K2Sik_rMwYpgMuOkW9BU#&|v_E z@$?;_S7Xf|7}f??-0JPOY;J*;w~8;?5}F`pqp4|Tc?0~Xr_aIdS_gqIvd#u|*1>C! zZILD3b&z2%Wqjg$6_k9YDmd{}f?(wnm3F}z@F|pA$%-z3lN_m!*rxMgKWhiSy-OY> zXi)g$PUpaIh{BlixgU^zX(@2yTnyapmYwE%8x5hCc6{NFjY8Qcd+kPc#i5h-{+NOr zY3O>t!r|xr0?e`X^$<8wj3R%ou{<2CK*_BkwW+m#@$yA$&M5zCEc#5bn22n^b)oxO zXS5n|R3PLfXL}QNn>cTC$ZEmkX-9k-AGG6<_*t0?AzgT2LPOG}s|(8>*EhIa@4;x= zEB!53`*5@V4e_Ga{V2L^$5rOWe&p_F5!Ln`#QgBa=J~TjNGE568)k+u;NcyA6SYwc zX4@snh@)s)#j!6uc?`p*I2*ZI$I%N4BQ-=PvFDz@jLp6&97-sDd?$Spxdv7Y9o?p| z&Hqzyv&$5&SbAF8=S(60SzqRFc~n&2qa0wJM8yGydoNeE)39Z0a`B;SG?Y)8-gHKI z8mYfHyyL@YD5aPpHuaN+8_#8aqQ=v(ztxPpcimi&hG=cjS;^%m->P1CqJs8CS& z$u!z}uIr0+(D3jLk4Uc z&I{*v=nZ2uQ|Pi&<1mg1?>A+s8btl%pQ7#xgV_7-u|mM{0ZbRNUU|Xahre>&HrI^w zpb}Z3zoy%befjR9bgw(ng!_c_nshrRn1^*om$qR=%d7L(9h>pBj_t(Rla08&>;B2T zdbQZ8+UE7Bx(X$q9B6o3T8@=0oNYc^3eifZKAxHaIDU|pRG(K@C0#%*`^*gj;pkFSSP}Q%5<;of=Lsu2-iT|f0vz$lCr@+R!teN?2S$&-pr-SD^yKh6Y~(9F-L#T{Hz|J% zZPZGzUk-!jj#OYo<^*F(ZZ(?6E(s1NH{dTd>qVvUR(x@nJ!_9gCq6m8>2*y|HyWd6 z$8G68eB^jCN^x`mkFZZnhxQNS@U=#v*EyrOO|)W~QAiWI}$@mk422GXLjJk?upR=VgG{+Chr5xm*R)hB|V%3|3_)+1AXM`b9aaE-Hyqa}q*W(fQo?n)v1bnaZO zvZs*I9jj5S9u)FNXtC}acwC>vFrLVozA z;Zs`4M&3R=`EAjcjf{A(p-1uv8@chAx-d@zD|xRxGD7wkD=B9`z@wSZLS9hL7CQBZ zh5Y+jX|j%)h0HoJFyau#OeVeB(blw^nLNQ@b#a65Kfgba>-55lk=#aX9|)u{l5z6q z6E2xEkZy&`8?xi*$@ERp6V+LCWR-TSp_C>a>AaSESXp%qso6ZD*6OQBocFW$e!Yb0 zn-ik{d|JfYDqj_h`sOh#^c_i|%>Cm`B=@_4X$<2jWDX0YVejN3XRyr_ns4rk8~-?t z>D08z8*(GK<2{R~{FQ#RbyaBpncIVTR|8-8I<#Zvcf+RliVfJkr`7xDlTtK|r1R?1 z%)+RqjNUgazWAskVYkqhEO_)%U0bK)543gPv@L3@hVz04-4K;}-iH0u8I57(laL$_$~YTe#v>s@B@+vB(Sum7E^ z?HN|v%i9=9`4e)P5yA{);9c4a?J7D_f99gro%U5c_08q{@XHm{(YxE5l(>i)d9X00 zK8Hay#-9%d(@<^vc6_SU6ka-iEXufZ7`Z4-W1_S^tiArxmhp8vnn@~z`ChETFMIoJ zhu(e1Ej-gdts4u#Q-_CJzr7CJS$L#>-)V=6ng!W2vI9Vy3}THh8wU>#gIc3)RPa!j zXiL{bvv~V?v>AKuWq0xLhaUft~)Uh zHYb1ZC2V3OtgIuH_W$$x*Iss;qDPsDn7C63M=M!~ndzfGiF2$(hR=ePz+*NdIwsH_9&I}jAUvm<3+rzGIpWq}sTO=Hkzj6@; z?b^ZuG%li8;Zrz=B{#9pJNo-VCO5(KV?D{5pNDvPGBk<%0}pX$Z{03eUmoIt9i(*G2+!8Qsg+<>B8zd?0lv#DM15@2xG0sGIBLkM zbu)*Fz%$a{7?&9d(aoQj&#*8OXAVpKw|ADF7?3rPE6=4PD*h(!&TL)>4`sHLzbvbe zrqtlx!?_F%R;#mmGYin~W@FyJJO_CfsvL*@`{%nFd9x>Era?&hReBNl}pS!%rbU4-X_LmQgCDsh8GB=4?04VYFoJ4(;nh6De4oI4jg@q=Sp zhZI*o9$bylPh}rKA9memE!IPLzSXKAL1q+zf7CSO%^32vOp_jc6WA{kbtJ}q3XeQ@ zO*?doijvl+tu<<>_<>hvbD9zjUATXD{|uSN%APw{+%2YYpvg|l^W+Rx-QKp`Mb6-m zY2(%O?irK{5!?6OcNX7r9CKSyn8l!Dx&Dc-XOX_cqV;_43_`BqjUx*)NTte5rt8e& zSYg96VX+zftLGM&!#{&}E%!fAXr0E5kL=dN3uySD0*9`bOg}y-pTaK$tqK}R)t4rbp7R%HZ&Bc$R<`<;|m4E$@8f9SYS>_NvAAa7( zp*(<_g65=B&-LT#RUVdh(QahQ=h4{N+KG9oJq|m1TTsy3x|n4{Jw}Po*fDRfMg`?M zy4LrlXeFJ_U8$9aDmM&XnPz^#%JpAE8`QjEB5&ZUM0P4%ToslQX3U3uhsp4cp13M#hxsB+6v;SAwoyHl11 zOcqPM53kX{jPm-huh0y{bQ+2o7|nos{n1TJ2WB8wgX__Q?X!TV9>#qSnuEJbAxA3? z&Oun5%MUNxIVcxOl;8Dx4l+H?KHI-I4}WAoR1@9v5aX5etCF$+-=0{zd=H<8kf^9b z=N%WIz;O?6mFfa)KJFGJxI7O|fdu0}-yJ@Jmu8{)*}c?{1=G;DLxH%KJPjQ8{l#3DXh1cwto=x#LEGK6^1gT~ zBnJ=NdQC@#Bs(cW!ezu5idJktcb_C7z#YFYd#vOh(=7R-~ON~BD8Vs=tA9ybH4g3U! zWsbYVKsSTY!*}1$f|c6In>7|OSfbo_)z&T(53&6j_cSX)n&a0fL-9XY|4yvUP@@_H z4SQtxv+9uJb5oG+g(ggG@)H^HYs9~Y3aF{qTQEHTPORIuPFx>KpIy&xL&g`2n^muN zAxHW9h-d%mJi+0J(IKv0EH*uNm2$Ed9cUL=Pkrn~&aORkAxw$a zjqW-=gdQ&URhe6c@guKIl3d*o#vNHbc7k~XN8^ny`FRfETrl6sQQIM08@PS*VfG+? z2>WWMw`mx^UX(Z8Y&L|AuXaZS(FTyh8mj%{=>YEfXg3j)I)Jz2?@kN!_T#35(hUJ> z16bR}^SKoIFz%Mzx&2AK82e>)1Yh-{n#rqM^+&owPuq z;0(E;&``V)Qonunx9@C*l``FyRJjw9m^XKRw(SJDfu^d{7TvHt@jp*( z=^pSu3QqD|J7 z{hmQc=r2A$cVYzS0U3&Rk?yUQwS_XxbbEGOJJItn_6OeFL8Mqu;( zPZ&Nm3Lg9PxAFcOg&NDRC1hc%Zo%g+B|?^-Xr*H9pFtFQZ*3yX&A}rcNJL?v3%<={bNUMh~~j1P$Zg ztL1|_Wy9FIa!JQEaLK~ATZ>pwhJFi#;==G3iK{9w%Ie0F^my%sn| z63(pQM`o4q*DGrnn0SX5YS;1avanE%<~n|vQ)7Im_0RXX`5#wk(2?W8>?tb8>Bt4U zsOHcA(UGTky%&t`(vepZE{_U|(UGyyv>J}Jbu>J@kIj&Qj@(dQ{AH^=9qHZi0ZZ4GIo5i@d zU!QLJHjR!Fo9OeoX&8EPy)?dO5?f2kl8)O=qSSc@(X;cT_;~l_*YZDxFrxpaWz9eS zNC_}nwM6&ghJi`-xqqDHX6Izhy4iqX-KAn$%hkwEKmLt_tiXJ=>6;Ifzu*&dhyD!l zP#Ac(qu9448J5@AZhb8+2PNOuPw{sfK%HOq)cBb;pj*-3_}6g0q<2` zmbwjda9y-8+x^fyXyxr)`*?N%mUgzA-`cbYC!qS}p2Wep-b_IklsYCDtJ9a{xz}>Kc5Veuuh3YY^@^rXwS?4qm8ocboY- z2wfvr->}mW?u-uHxjgIe?Z`U(a;M6~EU$v!)~$Ey-mC$m7x-7ru7cSfF0125 zRsp-Wi0xcj0S)E5qmj2);6S%u8ei=a6qzmRvaBwGK*46{c(nkWe!aIsE9T+$V}o(y)l7XwD^9=kL)9N$$MT0aW(d3!ORAAVx_MSa^65ihPO!+x60-?_IRXJA%f%%k@ z^-4w$MDEB~>iX0SPrIw%e$#J&oR{_zJadI0@~uZ9@4wG*V|VDv+jn`m9I|J_l43nF z8SF`~scl8(J(tc#ZSKX3r=ld|@`f<*(qw*--XwlD-V>KoN=3i&HsSxGXR!L(BL(@Y zd8F^ zD}OSR^4Hl48rNCK%E(*ZwkO$0)u?Uzq;)8yq`&mFH_zG0EJKMFl~xW?;_7`pwN6e_ z<}0(%=mBo>Md7aIQ(`=%Nq-vCVJ%+L&THlNt7KmCmEn*58p?d+yC$2K z>jWQ3YD$qo1V4Fn@xn&8ul%I-e;clrP4ko4j3R#tK>;#nm+Ea+QvuS!>z_D<3Xs*2 zbQKX90%TsgVB)!X0rH{7U5x%0zl^tr)pZGwSEuWFlQ#*H=i6@v)%_75e`wf?Y-tc6 zMK4-zkEszLv+sHGIlmJizbY_xZ8j1h-