From 517d880ce6d5cb95862a077458b772bc8dd4d487 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sat, 25 Jul 2026 18:21:02 +0900 Subject: [PATCH 1/2] feat(mokken): Mokken scale analysis with Loevinger H and AISP Implements sample Loevinger scalability coefficients H_ij, H_i, H, their Mokken Z statistics, and the automated item selection procedure (AISP, 'search normal') in Rust (mlsirm_core::mokken) and Python (fast_mlsirm.mokken_analysis, MokkenResult). Formula verified line-by-line against the mokken R package (van der Ark, 2007, doi:10.18637/jss.v020.i11): H_ij = S_ij / Smax_ij where Smax uses the comonotone (sorted-column) coupling. References: - Mokken, R. J. (1971). A Theory and Procedure of Scale Analysis. - van der Ark, L. A. (2007). doi:10.18637/jss.v020.i11 - Straat, J. H. et al. (2013). doi:10.1007/s00357-013-9133-6 Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com> --- CHANGELOG.md | 20 ++ crates/fast-mlsirm-py/src/lib.rs | 46 ++++ crates/mlsirm-core/src/lib.rs | 1 + crates/mlsirm-core/src/mokken.rs | 377 +++++++++++++++++++++++++++++ mokken.patch | Bin 0 -> 93802 bytes python/fast_mlsirm/__init__.py | 3 + python/fast_mlsirm/mokken.py | 109 +++++++++ tests/test_paper_features.py | 76 ++++++ tests/unit/mokken_tests.rs | 391 +++++++++++++++++++++++++++++++ 9 files changed, 1023 insertions(+) create mode 100644 crates/mlsirm-core/src/mokken.rs create mode 100644 mokken.patch create mode 100644 python/fast_mlsirm/mokken.py create mode 100644 tests/unit/mokken_tests.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index bf870b754..27b107a91 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -114,6 +114,26 @@ ### Added +- **Mokken scale analysis** (`fast_mlsirm.mokken_analysis`; new + `mlsirm_core::mokken`; Mokken, 1971, as cited in van der Ark, 2007). + Computes the Loevinger scalability coefficients `Hij`, `Hi`, `H` and their + Mokken Z statistics, and partitions items into Mokken scales with the + automated item selection procedure (AISP, "search normal"), with sample + statistics and selection mechanics verified line-by-line against the mokken + R package source (van der Ark, 2007; Straat et al., 2013): `Hij = + S_ij/Smax_ij` with `Smax` from the comonotone (sorted-column) coupling, + `Hi`/`H` as ratios of pairwise sums, and per-scale Bonferroni-adjusted Z + gates. For LLM-as-a-Judge item-quality management this flags evaluation + items that fail to scale (label 0) and detects multidimensional item pools + before parametric calibration. Complete integer data required (dichotomous + or polytomous). Rust-only numerics; the Python wrapper validates and + marshals. Tests include a brute-force covariance oracle, an exact Guttman + `H = 1` anchor, a hand-computed Z fixture, a Z-gate anchor whose deletion + seeds a spurious scale (this test caught a real sign error in the normal + quantile during development), a hand-constructed Criterion-1 design whose + negative-`Hij` exclusion is the only active gate (mutation-verified), a + two-cluster AISP recovery, and an `#[ignore]` 500-replicate Monte Carlo + (normal + skewed traits). - **Many-Facet Rasch Model (MFRM) rater-severity calibration** (`fast_mlsirm.fit_facets`; new `mlsirm_core::facets`; Linacre, 1989; Eckes, 2015). Fits `ln[P(k)/P(k-1)] = theta_p - d_i - c_j - f_k` — the rating scale model diff --git a/crates/fast-mlsirm-py/src/lib.rs b/crates/fast-mlsirm-py/src/lib.rs index 437e33ea8..5c8d98852 100644 --- a/crates/fast-mlsirm-py/src/lib.rs +++ b/crates/fast-mlsirm-py/src/lib.rs @@ -60,6 +60,7 @@ use mlsirm_core::rasch_cml::{ andersen_lr_test as core_andersen_lr, fit_rasch_cml as core_fit_rasch_cml, }; use mlsirm_core::facets::fit_facets as core_fit_facets; +use mlsirm_core::mokken::{aisp as core_mokken_aisp, coef_h as core_mokken_coef_h}; use mlsirm_core::rsm::fit_rsm as core_fit_rsm; use mlsirm_core::rt::{ fit_rt_lognormal as core_fit_rt, rt_person_fit as core_rt_person_fit, RtConfig, @@ -1470,6 +1471,49 @@ fn fit_facets( Ok(out.into()) } +/// Mokken scalability coefficients (`mlsirm_core::mokken::coef_h`). +/// `x` is a row-major complete `n_persons * n_items` integer score matrix. +/// Returns a dict with `hij`/`zij` (flattened `J*J`, NaN diagonal), `hi`, +/// `zi` (`J`), and scalars `h`, `z`. Sample statistics follow the mokken R +/// package (van der Ark, 2007, https://doi.org/10.18637/jss.v020.i11). +#[pyfunction] +#[pyo3(signature = (x, n_persons, n_items))] +fn mokken_coef_h( + py: Python<'_>, + x: PyReadonlyArray1<'_, i64>, + n_persons: usize, + n_items: usize, +) -> PyResult> { + let res = core_mokken_coef_h(x.as_slice()?, n_persons, n_items) + .map_err(PyValueError::new_err)?; + let out = pyo3::types::PyDict::new(py); + out.set_item("hij", res.hij)?; + out.set_item("hi", res.hi)?; + out.set_item("h", res.h)?; + out.set_item("zij", res.zij)?; + out.set_item("zi", res.zi)?; + out.set_item("z", res.z)?; + Ok(out.into()) +} + +/// Mokken automated item selection procedure (`mlsirm_core::mokken::aisp`, +/// the "search normal" algorithm of the mokken R package). Returns per-item +/// scale labels: 0 = unscalable, 1, 2, ... in formation order. `c` is the +/// scalability lower bound (rule of thumb 0.3), `alpha` the nominal +/// significance level. +#[pyfunction] +#[pyo3(signature = (x, n_persons, n_items, c = 0.3, alpha = 0.05))] +fn mokken_aisp( + x: PyReadonlyArray1<'_, i64>, + n_persons: usize, + n_items: usize, + c: f64, + alpha: f64, +) -> PyResult> { + core_mokken_aisp(x.as_slice()?, n_persons, n_items, c, alpha).map_err(PyValueError::new_err) +} + + /// Marginal-EM fit of a mixed Rasch / mixture-IRT model (`mlsirm_core::mixture`, Rost, /// 1990). `y`/`observed` are row-major `n_persons * n_items`; `model` is "rasch" or /// "2pl". `n_classes` latent classes each get their own item parameters. Returns a dict @@ -5066,6 +5110,8 @@ fn fast_mlsirm_core(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(fit_crm, m)?)?; m.add_function(wrap_pyfunction!(fit_rsm, m)?)?; m.add_function(wrap_pyfunction!(fit_facets, m)?)?; + m.add_function(wrap_pyfunction!(mokken_coef_h, m)?)?; + m.add_function(wrap_pyfunction!(mokken_aisp, m)?)?; m.add_function(wrap_pyfunction!(fit_mixture, m)?)?; m.add_function(wrap_pyfunction!(fit_lltm, m)?)?; m.add_function(wrap_pyfunction!(fit_testlet, m)?)?; diff --git a/crates/mlsirm-core/src/lib.rs b/crates/mlsirm-core/src/lib.rs index e55b47ff6..748091f05 100644 --- a/crates/mlsirm-core/src/lib.rs +++ b/crates/mlsirm-core/src/lib.rs @@ -14,6 +14,7 @@ pub mod mhrm; pub mod mixed; pub mod mixture; pub mod mmle; +pub mod mokken; pub mod nodes; pub mod nominal; pub mod oakes; diff --git a/crates/mlsirm-core/src/mokken.rs b/crates/mlsirm-core/src/mokken.rs new file mode 100644 index 000000000..934d96746 --- /dev/null +++ b/crates/mlsirm-core/src/mokken.rs @@ -0,0 +1,377 @@ +//! Mokken scale analysis: Loevinger scalability coefficients and the +//! automated item selection procedure (AISP). +//! +//! Implements, for a complete integer response matrix (dichotomous or +//! polytomous), the sample scalability coefficients +//! +//! ```text +//! Hij = S_ij / Smax_ij +//! Hi = sum_{j != i} S_ij / sum_{j != i} Smax_ij +//! H = sum_{i < j} S_ij / sum_{i < j} Smax_ij +//! ``` +//! +//! where `S` is the sample covariance matrix (denominator N-1) and +//! `Smax_ij = cov(sort(X_i), sort(X_j))` is the maximum covariance +//! attainable given the two items' marginal score distributions — the +//! comonotone (sorted-sorted) coupling maximizes `sum x_p y_p` by the +//! rearrangement inequality, and the means are marginal-fixed, so it +//! maximizes the covariance; the N-1 denominators cancel in every ratio. +//! +//! `Hi` uses the ratio of PAIRWISE sums, exactly as the mokken R package +//! computes it (`coefHTiny`). Verified caveat: this is NOT generally equal +//! to a "max Cov(X_j, R_-j) holding the realized rest-score marginal fixed" +//! reading of van der Ark (2007, Eq. 2) — counterexample: X1=X2=[0,0,1,1], +//! X3=[0,1,0,1] gives fixed-marginal max 1/3 for item 1 but pairwise-sum +//! denominator 2/3. The pairwise-sum form is the de-facto MSA standard and +//! is what this module implements. +//! +//! Mokken's Z statistics (null hypothesis of inter-item independence) follow +//! the mokken package's `coefZ` (`type.z = "Z"`): +//! +//! ```text +//! Zij = S_ij * sqrt(N-1) / sqrt(s_ii * s_jj) +//! Zi = (sum_{j != i} S_ij) * sqrt(N-1) / sqrt(sum_{j != i} s_ii * s_jj) +//! Z = (sum_{i < j} S_ij) * sqrt(N-1) / sqrt(sum_{i < j} s_ii * s_jj) +//! ``` +//! +//! The AISP ("search normal") partitions items into Mokken scales: a start +//! pair maximizing `Hij` among pairs significantly positive (`|Zij| >= Z_c`) +//! with pair `H >= c`, then repeatedly adds the free item that (1) has no +//! negative `Hij` with any selected item (nonnegative allowed), (2) has +//! within-augmented-set `Hi >= c`, (3) has `Zi >= Z_c`, and (4) maximizes the +//! augmented set's total `H`; the scale closes when the best augmented-set +//! `H < c`, and further scales are formed from leftover items. The +//! significance level is Bonferroni-adjusted per scale as +//! `alpha / (K1*(K1-1)/2 + sum of later step candidate counts)`, with the +//! candidate-count vector resetting at each new scale, matching +//! `search.normal.R` (`adjusted.alpha`). +//! +//! Verification status: the coefficient definitions, rules of thumb, and the +//! Mokken-scale definition (all inter-item covariances nonnegative in the +//! selection sense and `Hi >= c > 0`) were read in van der Ark (2007) and +//! Straat et al. (2013); the exact sample statistics, Z forms, tie-breaking, +//! and AISP mechanics were verified line-by-line against the mokken R package +//! source (CRAN, `R/internalFunctions.R::coefHTiny`, `R/coefZ.R`, +//! `R/search.normal.R`). Mokken (1971) and Sijtsma & Molenaar (2002) were NOT +//! read directly; claims from them are relayed via the above sources. No +//! primary-source derivation of the Z normal approximation was verified; +//! it is implementation-verified only. +//! +//! References (APA 7th ed.): +//! - van der Ark, L. A. (2007). Mokken scale analysis in R. *Journal of +//! Statistical Software, 20*(11), 1-19. https://doi.org/10.18637/jss.v020.i11 +//! - Straat, J. H., van der Ark, L. A., & Sijtsma, K. (2013). Comparing +//! optimization algorithms for item selection in Mokken scale analysis. +//! *Journal of Classification, 30*(1), 75-99. +//! https://doi.org/10.1007/s00357-013-9122-y +//! - Mokken, R. J. (1971). *A theory and procedure of scale analysis*. +//! De Gruyter. (as cited in van der Ark, 2007, and Straat et al., 2013) +//! - Sijtsma, K., & Molenaar, I. W. (2002). *Introduction to nonparametric +//! item response theory*. Sage. (as cited in Straat et al., 2013) + +/// Scalability coefficients and Mokken Z statistics for one item set. +#[derive(Debug, Clone)] +pub struct MokkenH { + /// Row-major `n_items x n_items`; `hij[i*J + j] = Hij`, diagonal = NaN. + pub hij: Vec, + /// Per-item scalability `Hi`. + pub hi: Vec, + /// Total scale coefficient `H`. + pub h: f64, + /// Row-major `n_items x n_items` Mokken Z; diagonal = NaN. + pub zij: Vec, + /// Per-item Z. + pub zi: Vec, + /// Total Z. + pub z: f64, +} + +fn validate(x: &[i64], n_persons: usize, n_items: usize) -> Result<(), String> { + if n_persons < 3 { + return Err("mokken requires at least 3 persons".to_string()); + } + if n_items < 2 { + return Err("mokken requires at least 2 items".to_string()); + } + let expected = crate::checked_mul_usize(n_persons, n_items, "n_persons * n_items overflows usize")?; + if x.len() != expected { + return Err(format!( + "responses length {} != n_persons*n_items {}", + x.len(), + expected + )); + } + if x.iter().any(|&v| v < 0) { + return Err("scores must be nonnegative integers".to_string()); + } + Ok(()) +} + +/// Pairwise machinery shared by `coef_h` and `aisp`: covariance matrix `s`, +/// sorted-column max-covariance matrix `smax`, and per-column variances +/// (diagonal of `s`). All with denominator N-1. +fn pairwise(x: &[i64], n_persons: usize, n_items: usize) -> Result<(Vec, Vec), String> { + let n = n_persons as f64; + let j = n_items; + // column means and centered columns; sorted centered columns for smax + let mut cols: Vec> = Vec::with_capacity(j); + let mut sorted: Vec> = Vec::with_capacity(j); + for it in 0..j { + let mut c: Vec = (0..n_persons).map(|p| x[p * j + it] as f64).collect(); + let mean = c.iter().sum::() / n; + for v in c.iter_mut() { + *v -= mean; + } + let mut s = c.clone(); + s.sort_by(|a, b| a.partial_cmp(b).expect("finite")); + cols.push(c); + sorted.push(s); + } + let denom = n - 1.0; + let mut s = vec![0.0; j * j]; + let mut smax = vec![0.0; j * j]; + for a in 0..j { + for b in a..j { + let cov = cols[a] + .iter() + .zip(cols[b].iter()) + .map(|(u, v)| u * v) + .sum::() + / denom; + let cmx = sorted[a] + .iter() + .zip(sorted[b].iter()) + .map(|(u, v)| u * v) + .sum::() + / denom; + s[a * j + b] = cov; + s[b * j + a] = cov; + smax[a * j + b] = cmx; + smax[b * j + a] = cmx; + } + if s[a * j + a] <= 0.0 { + return Err(format!("item {a} has zero variance")); + } + } + Ok((s, smax)) +} + +/// H and Z coefficients for the item subset `idx` (crate-internal; `idx` +/// indexes into the full `s`/`smax` matrices of width `j_full`). +fn h_subset(s: &[f64], smax: &[f64], j_full: usize, idx: &[usize], n_persons: usize) -> (Vec, f64, Vec, f64) { + let k = idx.len(); + let sqn = ((n_persons - 1) as f64).sqrt(); + let mut hi = vec![0.0; k]; + let mut zi = vec![0.0; k]; + let (mut num, mut den, mut vsum) = (0.0, 0.0, 0.0); + for (a, &ia) in idx.iter().enumerate() { + let (mut na, mut da, mut va) = (0.0, 0.0, 0.0); + for &ib in idx.iter() { + if ia == ib { + continue; + } + na += s[ia * j_full + ib]; + da += smax[ia * j_full + ib]; + va += s[ia * j_full + ia] * s[ib * j_full + ib]; + } + hi[a] = na / da; + zi[a] = na * sqn / va.sqrt(); + num += na; + den += da; + vsum += va; + } + // each unordered pair counted twice in the row sums + (hi, num / den, zi, (num / 2.0) * sqn / (vsum / 2.0).sqrt()) +} + +/// Compute `Hij`, `Hi`, `H` and the Mokken Z statistics for a complete +/// `n_persons x n_items` row-major integer score matrix. +pub fn coef_h(x: &[i64], n_persons: usize, n_items: usize) -> Result { + validate(x, n_persons, n_items)?; + let (s, smax) = pairwise(x, n_persons, n_items)?; + let j = n_items; + let sqn = ((n_persons - 1) as f64).sqrt(); + let mut hij = vec![f64::NAN; j * j]; + let mut zij = vec![f64::NAN; j * j]; + for a in 0..j { + for b in 0..j { + if a != b { + hij[a * j + b] = s[a * j + b] / smax[a * j + b]; + zij[a * j + b] = s[a * j + b] * sqn / (s[a * j + a] * s[b * j + b]).sqrt(); + } + } + } + let all: Vec = (0..j).collect(); + let (hi, h, zi, z) = h_subset(&s, &smax, j, &all, n_persons); + Ok(MokkenH { hij, hi, h, zij, zi, z }) +} + +/// Standard-normal upper quantile via inverse complementary error function +/// (Acklam-style rational approximation; |error| < 1.15e-9, sufficient for +/// an alpha cut-off). Returns z such that P(N(0,1) > z) = p. +fn normal_upper_quantile(p: f64) -> f64 { + // invert the CDF at 1 - p using Peter Acklam's approximation + let q = 1.0 - p; + debug_assert!(q > 0.0 && q < 1.0); + const A: [f64; 6] = [ + -3.969683028665376e+01, + 2.209460984245205e+02, + -2.759285104469687e+02, + 1.383577518672690e+02, + -3.066479806614716e+01, + 2.506628277459239e+00, + ]; + const B: [f64; 5] = [ + -5.447609879822406e+01, + 1.615858368580409e+02, + -1.556989798598866e+02, + 6.680131188771972e+01, + -1.328068155288572e+01, + ]; + const C: [f64; 6] = [ + -7.784894002430293e-03, + -3.223964580411365e-01, + -2.400758277161838e+00, + -2.549732539343734e+00, + 4.374664141464968e+00, + 2.938163982698783e+00, + ]; + const D: [f64; 4] = [ + 7.784695709041462e-03, + 3.224671290700398e-01, + 2.445134137142996e+00, + 3.754408661907416e+00, + ]; + let plow = 0.02425; + // standard Acklam sign convention: lower branch yields negative values, + // central passes through, upper is the negated lower expression. + if q < plow { + let r = (-2.0 * q.ln()).sqrt(); + (((((C[0] * r + C[1]) * r + C[2]) * r + C[3]) * r + C[4]) * r + C[5]) + / ((((D[0] * r + D[1]) * r + D[2]) * r + D[3]) * r + 1.0) + } else if q <= 1.0 - plow { + let r = q - 0.5; + let t = r * r; + (((((A[0] * t + A[1]) * t + A[2]) * t + A[3]) * t + A[4]) * t + A[5]) * r + / (((((B[0] * t + B[1]) * t + B[2]) * t + B[3]) * t + B[4]) * t + 1.0) + } else { + let r = (-2.0 * (1.0 - q).ln()).sqrt(); + -((((((C[0] * r + C[1]) * r + C[2]) * r + C[3]) * r + C[4]) * r + C[5]) + / ((((D[0] * r + D[1]) * r + D[2]) * r + D[3]) * r + 1.0)) + } +} + +/// Automated item selection procedure (Mokken's "search normal" AISP). +/// +/// Returns a per-item scale label: 0 = unscalable, 1, 2, ... in formation +/// order. `c` is the scalability lower bound (rule of thumb 0.3); `alpha` the +/// nominal significance level (default 0.05 in the literature). +pub fn aisp( + x: &[i64], + n_persons: usize, + n_items: usize, + c: f64, + alpha: f64, +) -> Result, String> { + validate(x, n_persons, n_items)?; + if !(0.0..1.0).contains(&c) { + return Err("lower bound c must be in [0, 1)".to_string()); + } + if !(alpha > 0.0 && alpha < 1.0) { + return Err("alpha must be in (0, 1)".to_string()); + } + let (s, smax) = pairwise(x, n_persons, n_items)?; + let j = n_items; + let sqn = ((n_persons - 1) as f64).sqrt(); + let hij = |a: usize, b: usize| s[a * j + b] / smax[a * j + b]; + let zij = |a: usize, b: usize| s[a * j + b] * sqn / (s[a * j + a] * s[b * j + b]).sqrt(); + + let mut in_set = vec![0u32; j]; + let mut scale = 0u32; + loop { + scale += 1; + let free: Vec = (0..j).filter(|&i| in_set[i] == 0).collect(); + if free.len() < 2 { + break; + } + // Bonferroni accumulation: k_counts[0] = K1 = #free at scale start; + // later entries are candidate counts of each add step (resets per scale). + let k1 = free.len() as f64; + let mut k_rest = 0.0f64; + let z_c = |k_rest: f64| { + let adj = alpha / (k1 * (k1 - 1.0) * 0.5 + k_rest); + normal_upper_quantile(adj) + }; + // start pair: max Hij among free pairs with |Zij| >= Z_c. Ties mirror + // mokken's eps rule (search.normal.R subtracts row*1e-10 where row is + // the LARGER member index): smaller larger-member index wins, then + // smaller smaller-member index. + let zc0 = z_c(0.0); + let mut best: Option<(usize, usize, f64)> = None; + for (ai, &a) in free.iter().enumerate() { + for &b in free.iter().skip(ai + 1) { + if zij(a, b).abs() < zc0 { + continue; + } + let h = hij(a, b); + let better = match best { + None => true, + Some((ba, bb, bh)) => h > bh || (h == bh && (b, a) < (bb, ba)), + }; + if better { + best = Some((a, b, h)); + } + } + } + let Some((a0, b0, h0)) = best else { break }; + // pair Hi == Hij for both members; require >= c + if h0 < c { + break; + } + let mut selected = vec![a0, b0]; + in_set[a0] = scale; + in_set[b0] = scale; + // add loop + loop { + let candidates: Vec = (0..j) + .filter(|&i| in_set[i] == 0) + .filter(|&i| selected.iter().all(|&sj| hij(i, sj) >= 0.0)) + .collect(); + if candidates.is_empty() { + break; + } + k_rest += candidates.len() as f64; + let zc = z_c(k_rest); + let mut best_h = f64::NEG_INFINITY; + let mut best_item = None; + for &cand in &candidates { + let mut aug = selected.clone(); + aug.push(cand); + let (hi, h_total, zi, _) = h_subset(&s, &smax, j, &aug, n_persons); + // candidate is last in aug + if hi[aug.len() - 1] < c { + continue; + } + if zi[aug.len() - 1] < zc { + continue; + } + if h_total > best_h { + best_h = h_total; + best_item = Some(cand); + } + } + match best_item { + Some(it) if best_h >= c => { + in_set[it] = scale; + selected.push(it); + } + _ => break, + } + } + } + Ok(in_set) +} + +#[cfg(test)] +#[path = "../../../tests/unit/mokken_tests.rs"] +mod tests; diff --git a/mokken.patch b/mokken.patch new file mode 100644 index 0000000000000000000000000000000000000000..e7c2883cbd1ef8a2c951123eac3d7b5be4915ce8 GIT binary patch literal 93802 zcmeI5`RD59ncz?Kc03k~X=p$K@zBr+fsp*C z@9zD?x4wNNA~UP9XvEGQ`*7&0%8ZP-ao_jGjmZD=zi%z>F7Cv?t;J6K`!;?(iSIWS zCl?yh zM)1Skzy`ecJgzg3+l4Fd7MQ^AlhWH>@Z{aay_h*r>sNl?i)&yW9@V&aW1NJp#=zF% zn`i|kZp43R#NYk+6h73;4{`T#xx?K@(KmPk&f76heFwHjakcH;a}ySRD$jg(2v@>% zqu>gx0f%v~{(kz}kE9%lx;s_@Khom5^q*@?`V{mY*$u!;gvQy~Dd)zY}9Z^_(YM&s?Dv z90it`1+#(&!QW2I?NQ+D$7uO5V3Mr-Sg0q-db@(l@^P}jE$M(~nXUMlUylM7;Sy;9 zUZ@K;8Ksv8=(!VDQ`h~mzdhA;ms` zGhj;10zDqZKP*H-$0&E=s$^ce3E5>7q-;NW5MAzubS0G^7A#vgBj-THXkhAg!5jF$ zRpyILW(2SdM~Hgc1!m-1+zF>b-Gojw1n&DG*4;R>eH^GN3mXL)n_pmSpkjD9sD@SmmqS#U=e&uzL~RlK^&N{aUZIKBXDvM zSZ55{LGxY$pqL~LZj$zf4}T0ez6Z7GD$vMi#$TJ zTT96BD&(hNbM6PEKNQZQtt=ihk##agWoe){5-L9bsLYdJn(xEF=ILe633uW#VMAXF z?{|yc+c~rr&=;N%uUX=mF}O6pr&PdE%#fLEE&Jr#cFYLgOwF}Z+5o97tmm3=JqfBm zidlGc%?)2p9)qM^{LQSNA96B$CXU_h6Bm>=tPkT_j>7U|x`EW=VI{6@IC_fWVeyjM?m*U>TfOaSB2N;v( z-i`0*XZ)T^F`l>_++B)xyYUyRDc#6yv1Ql<(cAtf5cN8Ay$=r!@a(Y{g6ANx^bcju zTQRCV_|=|!|1@?H*`QDGjGW`UzZ3t_m!HISVh~ws`77rOzObyq9#|Iydoilgn>#`J0Pt(F&a;o`ln`2VI5pTg9u<_x<>OZ+ZPreE+_1!{2WCNn;qek;~d* zViiOcpqh?9nXP!~AfN;*XqKIjfsC7W0=nmM_aNHt$ET=vB7PxR#HmN^1sItH5EIXy z3t#YJd_NoiZU;}ESe%bB-iiBAoWgG5R)hi&CAZ!La3l+*QdttFIMy7&6}Gz@hs>4vvLtFOZ^QsZ2MFpFlyG7m?p z+k!#017~0aAM$$9k3E*!=jyhbU18heFKHPw)QGH%ch$yc^U zeYY_WI!jst?S{Rz_i&|$mq)Y>Q?UziLPQTNZw@Aw~zLZ1|9@D=fJ z_9DVP9{-iYx*X$QUiM2g3-{oO$fw&2nY~+Dwlvfz%nfSLmwiCxuxxu33o+(*r5C>a z9RE+mlUHIK;}mOcACqs+K~T6B-=haYzdc=d0}2e*dGSIuzAk|6M7mb@0M))u%F?afB=1bH#h~)Q`7}M z{E$`|)1z z4UWFfu!+65=ZSS|VYBsEy;1?x_lq6GCi4vWa^F251C_SchN~?cqM_@5VxZRQi^WXdfJ)to#>F!mDy6}Rz8L-4|91O+8RAUyS>m| zd!;?POfwvZtoP9Hb;+9GAKi@aK!<;^9aqq{%p~zf&h=(o!G^N30o;29*FOY~Zxk4e zU-CT<$_$`I4~=Z3vQK1woEK$8<>Bn(*}A?y>&vY6<3Bb*Gvx}i&@A)Lsfgj%*=qB0 zZ>Rs7ZLclJe&O6Z#p+uUeJ}9=ceia;+3r6Dh4;gT{}>!}X2Lp8S>nXy(XFtk+)@^q zKZ*?&m$zD#Yg@$^RosL(O|^tP^qsKkmkR#SQX&X#NH~> z&+O0{tmpOd*ZgRU>?HQWvq{}{EP6IBjzuf>#`tJgbjhi`N~ECl(=JE1X7#OJxtXQ#rJoQ^A-taoe*$Mj0a zraf6hw-ta?1LhhBYutFnvR8b5S2UY-+3>zC`|1|9@6)ZA&F#b3d^k>KF{b*!ODB$? zfJe~U+#X~ODsANQc z7xym()x9RvuCWq{;n}kbGj2!wJS0?Ai>~?7 z3b^8IB0c2--3ooISe$RJK}D0YiblrNYjC*=fd8#5e_!{MKC4z!pvQRa6Pn-AbX{wi z6zijTuan?mQIqTX*Lwc-cf4%=F)!+Y+Iycx53JP_2fI>H@kd*&ID87L&r_-dydEAT z{rs+6yBk@jzQ0xNh{s$P*^eaeq8;gO^!|SQzqfLnzSZ<7!caYO{auw&c>M4Pm706; z+1BT@w1Y;&dN5^ow1Onwfp(ceAp=M3aXE70tWChBWNKPi7{8v0sGt$E9&~&0=kU(& zMlW!ft3-A6UR-^1WDe?;+Cy?cl8=5CEiA|at3ir%mFFU!QgssXlX+WpMj1~X#A*e) zxQ&#>^#K}C zp5SSqaP-~dj*3M49HzPO-NKVp&Y5|fF}|&f-!Ntm&ZK-5D_LYm{;kxkeG;RHPSW^UZ>g$MVrx@M-auN(VTQQH zowy6{3ErU}Wpr3uj6t-K*0spm$^xw=E0J<>;U?P+j7y5L^X1iSHVC%&Uyuv+(vl&KCYF$`qz zBb-leN_LR?Dl$}-&z=M2WhtL^tkYVO0wh{vDc_8CfkTw#_6`!{)jLSwt)9bv-?JMYFUN1uH`-0}qS=M!pUj zqTg6aAM4L>pv_x*Mc-%)a$w?|OJP?&jCN>T*&Q?u)Yq>THuXlOD?D6tFxuC-Cq^}0 zto6v6)qU09Td(!+i-di12orh(=v;F^ji~%>c$u|fbUqkiO&Xes%iuY2%1&J8oA9A( zGgf=_P3sd#MqSO&!e<$&eRSbJ_nHXq)Ck`cTZ@kH3igSREZOa2QIUPDXo~k%U`753 zkPkh$Z(>eE5AJxu+Fyc#PePi~qad?<5b~t!tZ?6rf2xP!*S(?}+Bq|Vy*(@`oVv>2 z1A}wNpXzHaWa%B&kOA+|KYSUu+Meb?wa%Ng)H;$u->qKj1aj@@fcyEch(z{hlig-$ zM`9mIIb6`v+M2`~%^r~Tw)Z&U^PqEqgBo5xht#<`re66Y(}ndyM!C-%uz<{cZ}Ce= z0bKk`{Qg6UAifC>+b?p_!k0UyJR*J>yCw7-&tNkor`Q(k12T@D!kXC&Ah!*)mT2}| zv9lp@s(aFJ+e?oeeaBPiVRyaPuT01kner&b$Sq_N$s|US)@Ftnk2H^31SJX!fKqEA{X_uukEjg-x7= zCNW=mmaKUoJ&ZxjmT&;GeRQwti@v$554zupr}578ylJmBrLRlBEi7YFrF};9Mn7Ok zc5m*kJ5@T=&;z^oZP`ECE7hP@~QF5gf)Oy%D&I`R#75VKmJNxIPvS_?$ts9;Qb z!rSojlhQ)ZwNUJpUZFXCoCu!7_v0D%iSEWsq~-cukF2n{-k`5|{avSB{uFRIhG5pV z0YIvHOe*ov?v@O&3755bBKg`}V<)D2QKaN*5IDH62dO1oBeXcU6i+ygQk1~lkV-tw zemwBOeR1rVcLwx*?!kI@OfNDY&%K4L*0Zjx;hHmgmgrvpp$VxCfimcC@Xvhs3nh^? z##|r6@VrK8A-fXzC$cjSN;APT;3xMM*AaL|vrs;xp5%ktM}}(n&(=etGgQ|eYt7~R zG80#!XU5W)Nk0Y40grRXUU*_gx~z_Ki6N1wm)p-hRH+Sm2H=-b_po2);t<9G0+A=L`LTBQe+52$4YO8cyZ zK8w}+IarIM@@ zFWIK?Bs%?hd?E#$+5m*sJLIcn*Wk&dk1X5exYEMKi3N;H|LoVRpK)dYU$NGCshyBW zIK0-mz{*}=>Qj`%tDpTT?lQ-oTeYT@+>ds&Jsn1|+2e4RCE8OJ)~h1j-N2vr#Y0gz z7^pa_fqd3htm=SK@;$eL=Qv@Y#oJEY8`nK>?H5Ihd|3WwHf4-V`pn98RT=NuuK*Wk4=Jo@;!$Ga7z_N znM-ljX4%s{o_AwF`5F&EDDD>hYi$h8ju%;!ZUfUq zZnleKXRSyxt-|f`(RqIya?W`G_+XrmfVUw(K-%hMPTyf}hPCE)bV~Xfv`Sy2rHJID zwnVgsyZIcmNnHQDK>lw__?5}ULJ6Or2j}DWC`*V9vR9@&B%1u;Qc7HHhrH}}l!RLx zHQBeXS88jTfpg?kr+}=+>t>ird=VF2S?apEYqm+BnkQw^(4DOot>qmY!7+&eVdS6j^CY$EBX=RWDuj`J?m;UpSRhd^mqY=Q*i&;v*s{YLivX zJly9*)C(=YD%h(t#+K@yx{4gwD<|BKjlIdpUk4_g@$i_?+CFF>(J}t2(HXm~-8O`t z>&w%w51ca}}<*~Pbd+>sfbV73_*SlPmI3A^M{;R!gI zd_39)O*dTWKuh5Vs}}E|?DcqAE%0SI0YMr=__ChRFWbSdV(#@ix-!3HynqIsq?H>o zV9ba7hbn*e${dO4IXz=9uG{n8D}BOOdqtar3+w@&8EuriQihppcu0GJA$*+{I`koH z4vgEL_2<<(U7gcA<8ht>5^hMPb^D(C-s;r>tpyCNy`mt^LVo{;AxZaRtjy=*okDve zUGfRgmYkN>s%^(*KfTt+)A+fs-+m6g#-S%U&&FDhKa9DexrcfLilBSsv-MK!=e_3N zYT>@Md|7zT=R(3OSX)<T=sdjmx_(Ay@74T=#aaD?v(97z#ci#o@;r2EnY>%h zZo;E6w*6^e9iIi9={lCq{KRT8SwGfs&}Qf_s+@X!^rO~u`fW#+HIg^LU+B-y(-Q&b zv?tQTtKFJCUWAqHWq(7CLI$D14@Im0P;~mU@I6_Hqf!G1@B+acbx$?_@vL|OSnxW~ zqdi`2MH61}m>}DaK3W^2DiSD|k(g{hI_2SH|HdK~NjR+gmZu{Kr*UdnkB1({UxuSB z6I1iDz3{XC%=uCi2(RaiwH}S#ScNh2;^xBZ#hPGwZxqdZvry_rv{2TNTnUi`+)gHO zJK95^+9#lvB-h(+$FDJ8!0?N*z(AzTX_(v5_cVPCsj|}f3@f%h?s}XVSE>1hDQ!&W1^P5()?+TFr!IpWP2vF%G|No6eA_u}8RCfW|3 z!vlF5+(YfpkMSS+@p~?8Q<`Mn^Jtso$DV@?vp5oieUxZN{jW)lgOGjX zueM-dBBkT$Kns|X{-h7&vF-zD+N=+0lGV__G3)rtp3-WVR`t~DT+dT+_87)A+xLq` zktSVlImfd-(r@5G@%FLM8bm7(0-G=5>aX#s8nvxbSMcI%mGyqs72o9xKT&B$YVX3Qu( z!?F2H!|lGL<-rU%7SGaW&)II@f1>x$PT-W&=z8Ai+gqFV*Zms}M`Q^k*d9S9Tbmlk zu_v~B-H5E9WRJGHNADLpX20tomCa-xnpJubo6Twpe@L}_Ny&=vC|~u)B#w~e;ZuTA z`iLHwDG{glfDC0Vsi!(#(GM7=s3U1N{-DuROKSnQ1)9EXj}_&x$S0^qmkK_tFrzvP z>kOf6xA%sTpQYYN)=%~`Z72SRYJHif^O~OhcFFppYotT8D?0TW*ygTwA>reI-7|q( zWyyQIbGoTwm9}L|i`9P0M*%T$zxE$uxjzdi=D+c1xLOMbuc|XZD$okv#f1)SeU!bx zPrs5hzRu73->^nNK8^GIz{7aZwf0_ErO>!u(U)&RTW8&0cUHKr+4!vR!j-@!z833~ zW9zGMd%q4W4(l=KtM;JDuMAhu(q8rcuJnebqRyMU=XQm{U+Rh_U@>$krfBj7Q1?!)&{*o`iXj-q_#1Fv!P>aLU z=rLwbW+$WTu}D&P(SnWy#(U`Xy&>*Dh4=A3#DCLJ0JZrJh=^H+J+*c`t$5Cs?CKm_ zY1n()HFdpjXIVi%@OFfZRWcg#>WB0EtZ&T8&AackXx|lWWefK6p$9vR^xlZsqZ{~? zx5Y2*C4Zx6v-{=#&;|_s3dQEF*2Fg@D*(UIGvD;dUvYrH>65oHthSP9VMX(-va(hx ze*yUKN3WdpAnmsmUJ&p4lFtjC@Q>t8BBSUK>QG#LWh+QT!87#FiZHFPZLD7sN3*Ky zy2jR8Qr}Ho1~sz}0y3=c^Z3OGd9r~r14M>Xb?oS{U#fA*p0Ywu4*hgg!?Ne%-GWJc zj^QpKX_gi?@W7bflZ55R)7Xib(F5E@wQN@Ukdyd1{xgPrGx8EdTGY?6?-*WX=RGkJ z{F}V?L-eNcWyxzAQK@1YT_4bh^Qw^$>2-Qm&69o)qaFE=Nl&12{^X;mi~KZx3ugA{ zvl79&nB-2oTfZuEhRU2#ko40I#c3*Oe%1k&!@oo?SWm5Yc*I^S;%Ta0lecu}TuWB} zsbrxCu<3l=8zG|Y3u<0PsJ0rOjNl~WYjxw{YiV=98*jhZ*!Bg1=f&c&5;(F51vTzcS)hzWHc7^@^#a;MRe81w7tcp4ldPr}NOUlR^w_k4^G zWx+6dLnrH@KR97;T_5M72ihX<%mjJkKi?()>v{kN#)PiCAq3u%+`}i#WU7bN)ctrs zUFv$J>!);T(iXh~Zy_hhi{=EUt)?wHSu)DGgzS6adHR7Pn45Z-=6CRc+rT?@?z)4679t*^@KHUv9OvpmV*H z&T`U;89uFA{s| zZnrl@@P44hc5oy-b1J?QXZA<-==hdFR)dKYXVcq^2B{ z&@=WREiZHF2nXQpB3nb5GA2xYyk^FPLZV17g?_@6CDy5FOApYNt2M^z>ovye>ovye z>oslb>!Ph8ZEXtPotz2YX2p8 zY;$?tIh0qZ``4GdPLp>4Q9#F)K+ryx8M*V@5L;g z8)018;}*iYqWylfvCo0#)=qdl3gIkI?cp3@Rxvo2h#gz@ch1K*VaB#!^T9(QW1)A? zB{kc%5ue%L?|d2Cd!|`fBM!59&^0rWy@%7T1FUX)^!Lr|2Cm3BX>Z$~;#d8y8Mx|o zufxbLlbMm93Qt=%C_}4kA+$V*R&NBvcqCIA9sg}?H{BEu??+!$DBuBmUyb+%TLwnR zz#V1xS(Q!rYyYXOa*!` z$8c7AjjZ$7y2{J01ep@`ircGz1b_9nTn%cqr)S}sAtaFsI_=+tqli5BBa@4C>K$3w z7k0b%K6T%hts54Ep7Qw~hDYHYkjI|VF?CO>@Hjmf-S)sx_>4rjAoVUbj+27N*Qc<1 zXliVh_7vVNwca&WsJFRRW~4b`(TLGm3tPYEM7qha)%AJGho4!21<#OZp9+)jCvjk3 z&QJJLKHb=fd#oz#$3Ln|uvS=Xc{Skg```)Ys9F`INKr5t;B-CF_xIry02%bw8bwVv zWjt6F==r9FPGc*tp_~U*3aTk3S3>W@v-DK?Vfb{~57DUWaP)pCBn#@oYrIjAUg30n zGOw~F%;~EbLVudazCWD1XWmk@u9lpdwfnLJf@?vBeH4!o5oyI}T`xZuN?9+~eL+X{ zfhL@iw0?#<2T=XhR!(uK=8wkO4jzFIh}!H;D56!wt{9|-XSn9$vHiKpDx2W!V;lGM z>l&U*-b@P$zO;(GwX9r%s~?0EXr0Xz&nMs_@9|^B^PAw?;W@EH179CTE~2pSo}Aevk|Cz3QBf2i3^EvOdKubXF*`ru)k4 zt%A^A{}&Ot?Z)XqpTu7*kf}O~ zX=578vfhY?mbhFvbGi=J^*K}}t8Rgs2>eEQ-{$9brxg43HmgfEjy&EQLH%`O90YG` z*A}ojU$H9c@F@P2U6n89x@+w>B$_9Vp*OsI(`=mvbkF`&9vF|Qo4P-B-Q3ker^x+gXhT6r94+cW z8~iN1nOq%s5x>%prA}&UKC*v~Yh{Z~_0K{Y??q46HRQ^$2d&-Fi2(Rid7cm41de@D z+c-79E7~0n+K>P1sG0Mbn#($z+9uR_RpnfPH6yVbOFV<`VI1bUv{4JsT#H`lDSg~2 zNZV&M&fI53F?-zPfr#$L<;{23W7m8)^!?j&7}92r#dzaT7-RlC5ghe0k>F8#@9Ey# zVm$Esq0Kn$=3wbt+D^!(a}~K0cPqn}I|@C&sn1O!S@$8mH%1kp=PeX?ariZya{+AY z}jvC@5h?Jx|5T7 z)URSVRfvBQe)jcPC%GJTf&BeEB8z{Y)THaj$dw&1WII2lBix_0FPUQD)c#&xsj)nt z1B_K;L+?Vc9?G>Ro$xU43=2c;PxaZJDWO@crH&w5>!bJAf6MYg3J3fBx29 zc*YTPW*}*wtIuBV9N@0xKzq;f#H39&g?p%r+A-GY3D4i3*eNZ$S-pbtGjltGSVeF{8JmJ`@j9-}9a8 zWY5TEP$fd%+9w6xk1Kj31N9(KipT=_?^}7TVC*h_1KwEnx?A`r&;R_i^a>wf!-s9b z31cUe-^O>~9f!_$Tc~suXgK}-L3o|!K1Qvtjg9w1NYttkeZL=5*EOoX_v6-Aha)2i z)UhbvJwD&IO$$MN4=6sZ`fDK>-Vra>c(3KZ$9{c}oiy58m2bEZT&r6F;icdeuNH_? zT&trVrtS}MhfbQNneVpv(pBKtEl;;Ns;>>lgTwV!)*d%^3a@c$o?E0{hWlC=##;$D ze*V+IlHPswq~qjEhhq%ew0(RYT!m-wRgBN73i(W~fCIyny{&3J#09bblNH#{`|aF) z9Z0`e_R-Ryy~@beajN9EmgCpO(0T7fPile1P!@u3@-?m&t7vPEEv#wmGllE>re!Tz zW9^3W;I%?N{-KpFFDlueHbKwtx6E z+(VT=r&E9{Vqx;HcrSeGt$ghxukV6} zA5*Gst5sI{;?n02qo=xk&c=0`$Qy68ZRc?*D9BsXzl!zbFCriKb$nim|Nk>ql&(iT z{D*Oc%;A@zr$3KxSDH~Y;y(q(iEru|>Q2sh+RyO&zpO`BNzLXGupRi#b#hsw=OxBA zAE$i5T~LTz13cpuea*%2?M1z`^5}d6pU~3$%r#3D*E#h(f9nb54UiP>>fAc7S@}eh zgsPo$3(K?yV2?WQjTyePg!Rm=fKsy%M23D?4L>zP_EGoHLRV*$GRGiEAL=btqVJX5 zCnx{xh4$tzxzPHT#zUH+mm$wwj5KGq|D>q0r(+SP({S3(!*7=*J-ytTyAlvcnD@0f2P;cO{1tgU6NfQ*mc?ZrM!D(|E-^laOMBo*ACA&Ixq zRkj4fwsTw3PVWPOi`srvSu(FD+RkmcYRu%I)xFg?f}Z4tYMKmv;(nI}a@W31WNO4d zFDG7rhd2*h6{5Gw>=ToU5?G5Qf70R=?X9;px0uLhIA_66kw?|}8JA+Lu?Tk)%LAYH z{u;F+=tm?0Z-F-QZtT^d)9TUTTfY^R@3P(fJMh|B&y#?fuGTCcHL%rI>eYCe=&H|R z#?scd2wdS#{^nOZr>uOGUZAeA#nkk>;)FMOd00&{Wms5gs(O}bcSRTMQ9DWtRmu>( z^e~Qj(P`V!g1N?MtvSBG&H7uG7;S4kY1<{N^w6!>x%T_l?UzH6zKwD5K^Tjyjr@?` zQcokb$up8#yiC>QP$hgNAmRNjx<=c!zLYChN*`dprgiyahn3beTr(gxP?8K`jMFSOQsNAqo^_7U5ee;%_W zL-0l5bSIuWD1Ue2t&%?mcaSYeeY8_>MBmJuH5>Q~&Y~TeM|SR^ORbIM&u|}Xkb~f9 zP7$+i&E6g9nv>xa&FNOqLKIkE^vnJHa+lCu&u5^NA}2dQjLHl zG83N&0mLoNlI}#>Xj*1RgeNZ-oXLtfAL|Up=i!}x9d}t*28&b}B`(G(1f6R+QKZjR z>nrvgBTx38;u7Z&TVM(C`dF(4 z-zUQ>p+~GRcp|@2w~=NctI5e8mVI-GvKiJY%8EdHo%E6VYisehTP;kzZkS`?fB;E@ zazjiy69G>F-}=}=5Vdga7XH9fg}bmA;A#Ck%q~ z;~rL@YAW8v(ntTs2=_;m1{JTv$Xz1E)bV0@TP^DqG%${DSmz4WkJzPkD5==k4bb8OFg@Elv2z}tg5 zqP95RexP?IXLXNuI`R%2GS;`_Nj;Z2Hck}uNG&#>hgKxBg53Lj?K86^Ws2hPa+@t2 zVgjG|OvExD$LDG<7%F+@DIuQ2u+kYl(4OHUOfpZc4KiMf?J>TBj_jqFk~dbA=4q6p z;0X=Z#tlG3ytT?R=X{iLYMY<*R%E!$Z-6FKE`E`=)BI9_@$+T^Y>L0BoY2ekvlNmQS1zHa*RxM}&* z-0sgiLEkal$#!d{9BWM*t(Pe(dhrr-_%(VU9u}tQQ>QI>CEYOjv*$ZvS~MBBMUPoW zrDTD7=_hJj(H*KEz^bw1yzfzM)T^yObjystd=FcG-lc(7FqMZoLs{ASzI|}sRnd=M zXE>PI6bw^mRvA-AgNhlYexJ7geY>uuPOM_O%9FYOhmJ1$R*JrBb9@J)~PK|?Qbu`odI(@<({aXp-BHy_eh ziAlzu2jP79rZD+)j8t=eV)EzUr&X9t$(@2VTa0s~bcD2351-@*PinmW9u12Aqu$@& zc~=g#Qd%Ls6Q5{WorOnyp0)RWV$-MLe4r+$iZ13viW{MCZ-)Q%Zj8YB(KMdp#)FXkw?@PU-i!$6LeX?OIg)c& zWy#2Qoho`N+l#YVU%<|c!Hv&yv-F3~lb!w`C=2(>lW`UFJ< zsFQJDYY+HsnWt$n3MTB^=>2Q$efu5x{ugB-TBv>sOu9<*ac~8>K-TW0i_k>Son0Hi zO`W-XNJIE2)*tABs4?SQuY+<9gSZy&8O?*NJqfxJJ8&8lnGjZHTBt;2ajyC$GlhqU zCk9uGNBBghSz%?3wB{GqA@GFpPtc}bnNKKL3u2~3mfUw6mKszFjh%_w8dE;Q&El0tXk~%GO`_dq1L*6FRcy2D&yfvMznTG z55PVpr>lEpU7jKV=5gxT4 ztURxV#vFjXK+{vZLVtpeN}{WAjd4X6L)ucJ#aZ2FW+D!XuHx+vx{oi9pI47XKEUv{ zBOo8P3ii0Mrt-b9VPyS*LnmytRtwKSlHxvUd(w{dIQldv2~sRE6DFbhg$~juG2>U;BYuo}BnuB#iJ}wt zoK;q1Q}B^16<5k%%eNTdzY2@!#vkHo-cBG)Aa&qc9EXP{%d8y#c9D5@owrzg7SFs6 z2*K>jL8ur*Fs{>s!p(Pq1?)dHsTKwIk3l)&r36R;1gk2?S1P6 zTgDKy?LFRl(PF*TF7OB$V`}ro8&Gn1jdp&c>G)WyMzQp1onj2m_d<5B#TYsRPv_kk zp7l?ki}*F4z8Rc02gcvWco_-lTpP|{I|wWAApVh4R3rtYb8_|~4237%e=YD*1 zZO>dh)b_y`QC`e8f_V5=aL!iDZp=zabJcmsv{KlY77}aF)t@aSCnAFGp~<)RsDJ(4 zJho<}by*B-lcLe-IQrj;R`Yu+i6rsieHxG=Dk5Y0vIk0C09gRI#H%*k=bCEObcbJx z%ZTUe)@RCkplNQ|X)!MUb*EVTodP9V7;kq!;^@q7!$rBOKeh%(90~`@8n+URHpbH> zKY_P_#q(Nz8#%Yo;kOZ0I*xP{FYB-RJ@=DRLY-(GEoLt}(UT%a$Ed@o_^&-yT&?3< z{H=>a|2rO(rul1#2Zj@*Y6?=F&S%@@ngO9ffL|jqvfmVkVah7SSs|70;P(c(ht!M{F2fAwp`Xn*bySLKN!C9p8DC%Sxg)-l6z56xi8uDS?jFuUdt&GJ2zAg4;33yg$xA~C<^5UBE&3BBz|GL^c)Ot!ygqz~(LSyj2}J)}F<&CwaHF zpFJDx^yY8g(awK8%h`u54F6hihRlEytk<^z@yVEna;oATe0$#h0mm?3^CfGqWCY>C z2XR%a*;p*?*R9unITJ=27^=7ZvNndCwXh-UdV@7QC}~RGzaG}{dZ7TEO~pN42hdQ* zAgX9KR0FZb3oXdSv+4{jMO7+&XiYXB{fP|G(+@Et)kUPP;8an|k@Jst3O-t1VyCP| z6Gy6+MNteo4j6h~QSC+Y3Ngy?yxCeC@rFop-I<|ZEf?N*ARe=KOeIzFnBH8^yW2Te zjyJ(8>#fx&AW|lswWoLDp3jYP{PZf|I~Aid5)ckK`BX>(?|{D;<7*{_to-YP5A4mu zNqb%`dJ--4ac~tp#HtE>hqgd-?Z>@mg_dX!cowd=9YivR^%RZRhH})zUM7Qs;4tGUB!&flw`bE%BJaTuE)1xT?oHIr;4~QuccjyJ5+HE zdD(p#3ipx~;96ahR&Z-$)xQ%G=1yf8bK$M1vPft(K_H%DRX}Z~AyOXIc|M&dg-&)) zXjkLPaQGAd>pzMH=iB>Bc{w|?{U6Gz#l@Me-TQR=8e-X%XOn$XrOgTn7-H8a_(Gl& zU+^z0XQHuZhjr@PuBMq~;Yv_U{)98OwKmJ@sXvBXJPUnrDMn_k1?e$I8}B_I)IGf# zZul`~E*vISv1l#dp&L{^1*W`4Kh%_l+r0MU;Y#W1-lzDke@{Ema<)Hf>C#xdvy^Hp zSr%4ok=!RmZr@wV1b!HvmwXz$g7?BN|1eIb(;6>OGDF(von7OUwv@IO{Ys<7@BnuU zhshi8ssoH?&$4F9EMctN&g-X0+55 z(naz&xr3*wHE*tAWz_!a5?{)LCM3L>)iDY!rD;<0{3#&j{nB@fmO|2c&1!yN#g2J1 zX8n$ap8jrcZHEl=oMHjJ6;`9M#)0n!<;UI*Yg_X}+?SW9*pX*?D0}V&0~+^nu^3bi zSTb5k1Q$Rsj=h;<_;xzXB|BWPT@!y5%!6~2l#F}u>hK56~{d)kKj+^lT(B= zqv^T<*9y2o2beUrDhc>CH6xGrEPP5ZE~;q1liy#<{IHxuoZ8BP13WT8`y#%fqkfJ5 z*gP}|w9`6x#zh|6pO)6C-_Mp(<#w%9|1sdBRtjAwI_b0@{8FM%Fp4a=Ruuo6m|8I& z(Po_i!sd|4z}nIS{f?i)(?98zVOtlk6$mzZw7bJSiB&k7RY!f%!y z8A_AkW=FiLTlSY)iWz91Jo{Ox0s#Z+m6f3}J(0ZuIqf;;2z;rx*6vzz8JSZ%6a`2V z>f{zI`TI-!|6w@?_z^y<=W8i)kN#vXueRr*ZuFeN%+te}uA-uFjRw6I?THw{7Fa@` za+m!c^p*$!Y@Y|1LYyh9Vnboq4 zAA6&tbIk=n{yMloGX8Wa-9YcvnQ_b6XC&L4FcM29fRi4fq|Xla_r190YAxAJYEQ9T z!?UEYJ9t6XOR8AM@+wbB&KQ1gEjwAjb(URv3{FZrt6ZII7#^Ck5L$Dw&4U}UueMUQ z+HFRXKFFfOvD6qTFX_mfyJSh^1@09NC$op0$M?}`#{FB|z}t4=b+G5WA6}|5GQfo_ zSrZyE@`0KYx}$zex9lx+v2=|2EcPCNR<+?|klDkN^;n`OquU=aex`GeXiHXghP%1s zClFns6--NVB#eNUcD;PFmnO?LWsj|$R$Zrsq1I_*<*ot7H+{OhF+(lM>1?vW4j;2K!s8b1JqgvgZnTn_}9oe z_n6Z>6#O%=X-}1>@KLjXVQdyt8%=hC@-9q5;wBW}FyB-cDyUsG=XO zL1a&KZC~8)e>payJ>0(~S{VVqbsfe7t3Hygu25#u=yVUujczzve%Dh zIFVm)GCq>2**|Sq=RIwB@i^zg>yafF@^pcdKIZbbasb47lsSa=4OEql$Cw!$fVz8TUJn7r{ch9T7?oVTV%;yu;(3-Zab;RBCG`7#X zA%+&v2 zGUMu4`O9u!28dqxU>6+iK1UIbyH;0$zLOsOBYK*U^U@hH#(=`~M(xRSgE!TdI z7-QXYl2T+7JA_=#l~F@#t>;naK-)vz_b~_ZDH)u)-pfYI^}H__UL#B3nV}=$^Th2s zUgwO%&q5EGi?7A+ADoNVYG@2a8`bTG_ER+vnZLbehw5kH$rE+% z#T{ZDvZndw9mkvSkaD7%>_c`W<9Mo#^dx8a==j~T*oyD)4E}aY@DK}PGjt-PEtb|2 zR!3j%S>v=h=N3OiN??delFY&4b)e0${_Ld2-jQopy{9aT&ip3JiT0AL5f`y$r8t}H zuWCG?DK%;2SmjF**0W+sZL07y58H59=k{?_rdu`MRAyy?Y(J={vBExAmU(`EBbb*7fUM z@brGPw$}!>!&s386-U55KJ6qua^CqueGVKl-y>hI2c2PE?Quv83qJe>d;`@fvReAx zh(@kuvi5V~p(%UlXL{dw)idUo^o_m0bb1^=WBXKVJF~Lma>uqfd%Qw-NSZ+dw7#tM`?^-rz|*aCf1Z_?*Vk05qiArQf~I<& z`5K&A^JYHr7apYZ;s^1a-C`%pU#d;WFyk5FtH@KK&Ya32+O~X1WsP{;x<(^WW1&50 z`R#bAl|t|3we4jq=|lg@`~7wGsh;rV_$>{pZ}hHj+$95ZeiS~ur)R~>%62D9U;d4k z2sZF{p2a`uB7QwTd~Q2(s>frDGx6k`fk9yL^JgcJ_Y-Hz?5T67!dLa`Ckx)OamK~P z$m>#Xcs~BAP0KUZ4K+VgDV{w*4f0CoL#vU;X63Z5i-Ws?X;`tj-YV-XLE=m2sOe2L z%=4nR#P^t2Biqj?Lkm6)oljp#=z8qTfvu+U5Vm@ZA#9RlMX|&tde5A_E6y>oB}Xsl z1b^D`OXIdwUNu+A@Mbeii0An+5{H{yz z2uGO~{~6!mRY`>^!LYq(SuomG{7Z^CqX!wnUq(}~Pwva|-F1svDI?Oh-M7W}_^tOr z?nH8j=fFx^&L1Cu3@=`XYNz~;!QtKvRgLLn3UIa^xMS|h0syV5%)v7mh^(}1eIw2m zWG>!iWK6>ojDaPhkNuD$&urS4F|Cij-{il&8T2B2x4AcyyWF50oNT&eL~4JA5v0NW zjW?8z#i{P2&Y9W!sCCTs;F-THSf!}pp!a%hSfeRQI2JOx6S9CW3opaT_UOG&sg@JV zUJc1nbrJ6j>g#>R8we~``?B6N!QZ=}%3+2hZ6F7c$=n+^m;Awc3g?bHk zzoALc5mZvMy9~*e@8^BDTJ_UgUsPeu`&yqxJ2*xERWrxFg!(!=vcWsMVSQck9oJjE z;5Hf2`C0X;=NQk|1%t0M7JeL->3)0;>j~^9!fRWLKbDybBhcRXQ?*D7+h+kko{9Y& z{4rJMOMl`mOMeQs&!RPU)9tO%>BrY~_M6_$oADr=MojIKgy2qFU)~7P!o{fC|BBpY zMiO3;mM71dd=9@u&4lc9$e!&T=id>V8fqZ*`yQbpcFk?DZTN)jMWa;@erT-z!#hnk1pZn1M{p|=+wE$u0BjHn04xwgx^FZ&{ybF1&8 zSM9(Omcbh~&34DnwicCWiG1zT7_GIa&c%B+r;6N;(Nwaf*o(1|JI|HuKbd_r5j!2> zHQj5m+}dBYpdX&6R!V!dYwm#--Yv?zmIQsRHGu-JL#KjDw()0;kvWz)F)h7<$J2B5 z^-=`}%&Zr;v?cG1_FHX{TOB36bTyv28qcZ%2H#%#9zCmlBGM5T%1RS_ezl9O4M{^TjO?;)J`l%Mrj(yQyn z>LDD%#&p(09Gg4#R2j~imTUjTE9f=70i0S{Y?-1$^4a}*s^b(^6CXuCU=Xb?dxbue zboKl)o$qyz*AJgXN2$(Sy-q)C>O70~$N>_LKBuY=_#VW!gOF_2nOIeMV+r+eUhc_$ z#?)`e)a%y-z$4h%bX5Q?C|?TwdZ*Y;>#8aHNqf$Qg%(*mK(uL>P801=f-rwXNzfGMA>i z8tt;TAusf>v@7SF?^{2w7@ORhAX>dA8c$`@5gCtJlS^R?WhqzHpm<&0{sX#q92@F; zQ;5JpW>A1hF~uC1ur$yucQcMd<}rq_Sq586${ZcTWE_-IC|yrRPVs(HOS+@3qH{lW z1Y(U)bH|)m$o8VPHOD5rs5&-k?Rx8%)^J&4NuKOO^{Y*u(OzJw#apg0c~2TUEEJ*P zvNJ_fy-DE5`?SXkATO>ZV_fe#z&Ra%1A@?`d(7i zbH@v2Un3Tin1b=Zev7lr?`UlwkCRx$UL}~(JIOhB%Mr9!ft`8U%(HB9OpadI#cQZJ zM0=lyR@!oEKV)^Se7<|W_c`1T-ZoN>1ap?tR68m{)N+&azx&=~>DrmXp~P>tq>O_e zIT3KdpSOaBikN^^J3qvkP=H!xEOd=KBICYC^6>OXwRGQyV^V7)~a{n{Jd4mu*c zuV-(yaZ|69knhRbe*BB}B=*gBSDO>iGxy72TFw}ymhcslaLekxWWvod^^bKrY@Tl-aHm+R44eKtlt@|94Wg{--vP8 zZ+$=h>lFNiSGmk%!B2IsFAIxo$Y!wAn3@ZXxwiFIdhfH{%efb`XWq!LR-deOa~I7{ za2=5>z9IOqA0}_Juj9djDvxJ>sNcPrk+WCOlzFSOxH-CgGTNrKul-L|k7CWWI-oP( zX-|fq(T2MaTO4_>L|c7}N5K^`(c>7P((7UPiTNCF-B|6_3g+K+0ZqyAi-JWWRsp@F2W--sAmIfFVL<-X!+X9}sl7-}NJ>iOn5 zY0CGbbIFR~Z4#G^@3ptwk@bcq_GbC^qRaf)PPHTDc`f+ASW@K$v^#>j5!L0^QiKgp zS;In56{1~=K!L?kRla#sK6d|pF6$mcW}J`*WB%s9vE6fnUO6e1J?=!;d%>I7F6EZ-iM8^FKgvq4-Xz1B1K^nz zq1AdBeJsxM8d>{J4rOq?I_Wlfns)X-t%cF0^*dN}s=suw&M8X;E8Z*yDH zUE|l)((4&&v)%{nB^XF}R&SJ59#p#^MKf#z_m2hk*=5Mp+$F);P-jPXDD@mTHtyA+ zrTjZaIJblYJ{*Stx{?nf{>OIGida-9pVk??DLLq&@Cp>WzJ)>l3a7pzt0Dh?HGcZ+ ze`@5HgRQ2pp{b#q3)%E}=(Pczix z26~)__G4_Vk=qi=C&t@+6|ew79V<@b1CQ~C;_J?XXpWxU3u$29{n_2jHkUX+<526p{ zTUpa`j*b3oQ_}*tCWTyfy_VKzrt7nr>$#DYZef!jjvYdWD+-c#r#Gl@MV>nIBo5(h zS@OSBJCUPgT@Y+?N7WDP76*sQR0A`+)6kP(+V4x|8P#jjli#iFTsheBiBzn%oru=N z20oi!=d!udLlon!YRrP2gv1!cnnWBlzB|;u;dQ&LKe^$=T^jj)&|4VU)+mhGe?TG_w?Im z+#Lle%SwcY|3geWCuixXZA`0(zwNd7)3bJjL|x<7;_u3s$QeEsvgc|ct*?x0CZfh# z|41W04ZV?PLx@7mChbG&P($ULWcX8=Gt5uQpp;`_6OL!^gEDTHx4HO2>3mt*3c+(9_UJ=<-@qr>+%IsM8TU=33{ zEo4n(clwNK>mf_mq3@rRoR95rTlxKc^sIT~pP&z=6VSJNF}gVQ=(`K*ciPOcC*T&I z0v<)i_y;YuxsIpwptPx}1D0MzR+1W;bven^E8n(B=O_Dlbv6kP5Uy38QMh>=*wYB( z_^`G*9=H|<&s_tA)2RdBbhXX0?~(ui(UWAJwGJ=cSnDL=Zrw_e=icg-OGPum$y-r{ zgeOOa7>`@k=i)GFCdCoxAUrNK8uBWgf`5S&5I^;X^YevcOYP^$u@vQp~})Iw8fLvLpuFfm`j%qDN_Wa&wx$K3iOYe?0fZC2=VJ ztSyNmGoU77(paJ@Z&UiaxaPGUw43T;Yuff$_#*mHZF}9yF#tTI_ne^Z!4KF6a<7i5 z>b>^RXdyoGUbchAUVzW~Y*TvE88O3mM%jjVrhBFD;kp|-Ju`3fPXA1c+~v^Q?5Q4n_U?li z)!y51eX)NUs_AWG(^&H}#-5%rrm&rQr>gpyTJv!)rj%#CdM#D&UL+1<_0JVX`<=b$#&Q%hvVv zu;*y$n$e(N?JS3~A`qn8b*7-E*AXn?pc68N3vZhC@tZK=& z1vy|Sdofx}?T|FR+FHb4N1_^J;zf90i2*Wi;}8#X8y;t~87F0HY)aFOW7?4|^Ck{@ nnR$EEokGnv&dPXG_|OS_<5=*V-CxTSt7b9Y-w> MokkenResult: + """Mokken scale analysis (compute in Rust; Mokken, 1971, as cited in + van der Ark, 2007). + + Computes the Loevinger scalability coefficients ``Hij``, ``Hi``, ``H`` + with their Mokken Z statistics, and partitions the items into Mokken + scales with the automated item selection procedure (AISP), following the + sample statistics and "search normal" algorithm of the mokken R package + (van der Ark, 2007): ``Hij = S_ij / Smax_ij`` where ``S`` is the sample + covariance matrix and ``Smax_ij`` the maximum covariance given the two + items' marginals (sorted-column coupling); ``Hi`` and ``H`` are ratios of + the corresponding pairwise sums. A Mokken scale at lower bound ``c`` + requires nonnegative inter-item covariances and ``Hi >= c`` (rule of + thumb ``c = 0.3``; Straat et al., 2013). + + In LLM-as-a-Judge item-quality management, AISP flags evaluation items + that do not scale with the rest (label 0) and detects multidimensional + item pools before parametric IRT calibration. + + ``responses`` is a complete ``persons x items`` array of integer scores + (dichotomous 0/1 or polytomous); missing values are not supported — + Mokken sample statistics assume complete data (van der Ark, 2007). + + References (APA 7th ed.): + van der Ark, L. A. (2007). Mokken scale analysis in R. *Journal of + Statistical Software, 20*(11), 1-19. + https://doi.org/10.18637/jss.v020.i11 + Straat, J. H., van der Ark, L. A., & Sijtsma, K. (2013). Comparing + optimization algorithms for item selection in Mokken scale + analysis. *Journal of Classification, 30*(1), 75-99. + https://doi.org/10.1007/s00357-013-9122-y + Mokken, R. J. (1971). *A theory and procedure of scale analysis*. + De Gruyter. (as cited in van der Ark, 2007) + """ + from .fitstats import _core_module + + core = _core_module() + if core is None or not hasattr(core, "mokken_coef_h"): + raise RuntimeError("mokken_analysis requires the compiled Rust core") + + if not np.isfinite(lower_bound) or not (0.0 <= lower_bound < 1.0): + raise ValueError("lower_bound must be in [0, 1)") + if not np.isfinite(alpha) or not (0.0 < alpha < 1.0): + raise ValueError("alpha must be in (0, 1)") + + y = np.asarray(responses, dtype=np.float64) + if y.ndim != 2: + raise ValueError("responses must be a 2-D persons x items array") + n_persons, n_items = y.shape + if not np.all(np.isfinite(y)): + raise ValueError("responses must be complete (no missing values)") + if np.any(y != np.floor(y)) or np.any(y < 0): + raise ValueError("responses must be non-negative integer scores") + if y.size and int(y.max()) + 1 > MAX_POLYTOMOUS_CATEGORIES: + raise ValueError( + f"responses imply more than {MAX_POLYTOMOUS_CATEGORIES} categories" + ) + x = y.astype(np.int64).reshape(-1) + res = core.mokken_coef_h(x, int(n_persons), int(n_items)) + scale = core.mokken_aisp( + x, int(n_persons), int(n_items), float(lower_bound), float(alpha) + ) + return MokkenResult( + hij=np.asarray(res["hij"], dtype=np.float64).reshape(n_items, n_items), + hi=np.asarray(res["hi"], dtype=np.float64), + h=float(res["h"]), + zij=np.asarray(res["zij"], dtype=np.float64).reshape(n_items, n_items), + zi=np.asarray(res["zi"], dtype=np.float64), + z=float(res["z"]), + scale=np.asarray(scale, dtype=np.int64), + ) diff --git a/tests/test_paper_features.py b/tests/test_paper_features.py index c0a6b2964..9835d9c0f 100644 --- a/tests/test_paper_features.py +++ b/tests/test_paper_features.py @@ -5008,3 +5008,79 @@ def test_fit_facets_rejects_malformed_and_flags_disconnected(): yb[6, 1, 1] = 0.0 resb = fit_facets(yb, n_cat=2, max_iter=50) assert resb.connected is True +def test_mokken_analysis_coefficients_and_cluster_recovery(): + """Mokken scale analysis (van der Ark, 2007): every assert reads crate + outputs (hij/hi/h/zij/scale from mokken_coef_h / mokken_aisp via the + wrapper). A perfect Guttman scalogram must give H = 1 exactly (kills + covmax mutants), and a two-cluster simulation must be partitioned into + exactly two AISP scales (kills selection-logic mutants).""" + import numpy as np + import pytest + from fast_mlsirm import MokkenResult, mokken_analysis + from fast_mlsirm.fitstats import _core_module + + core = _core_module() + if core is None or not hasattr(core, "mokken_coef_h"): + pytest.skip("compiled core built without mokken_coef_h") + + # perfect Guttman scalogram -> H exactly 1 + guttman = np.array( + [[0, 0, 0], [1, 0, 0], [1, 1, 0], [1, 1, 1], [1, 1, 1]], dtype=float + ) + g = mokken_analysis(guttman) + assert isinstance(g, MokkenResult) + assert abs(g.h - 1.0) < 1e-12 + off = ~np.eye(3, dtype=bool) + assert np.allclose(g.hij[off], 1.0) + assert np.all(np.isnan(np.diag(g.hij))) + + # two independent Rasch clusters -> two scales + rng = np.random.default_rng(2013) + n, per = 1500, 4 + bs = np.array([-0.8, -0.3, 0.3, 0.8]) + t1 = rng.normal(size=(n, 1)) * 1.6 + t2 = rng.normal(size=(n, 1)) * 1.6 + xa = (rng.random((n, per)) < 1.0 / (1.0 + np.exp(-(t1 - bs)))).astype(int) + xb = (rng.random((n, per)) < 1.0 / (1.0 + np.exp(-(t2 - bs)))).astype(int) + res = mokken_analysis(np.hstack([xa, xb]), lower_bound=0.3) + a, b = res.scale[0], res.scale[4] + assert a > 0 and b > 0 and a != b, res.scale + assert np.all(res.scale[:4] == a) and np.all(res.scale[4:] == b), res.scale + # crate zij symmetry read-back on real data + assert np.allclose(res.zij[off_idx := ~np.eye(8, dtype=bool)], + res.zij.T[off_idx]) + + +def test_mokken_analysis_rejects_malformed_input(): + """Wrapper validation: incomplete, non-integer, negative, non-2D data and + out-of-range c/alpha raise ValueError (each assert exercises the wrapper + guard in front of the crate; kills guard-deletion mutants).""" + import numpy as np + import pytest + from fast_mlsirm import mokken_analysis + from fast_mlsirm.fitstats import _core_module + + core = _core_module() + if core is None or not hasattr(core, "mokken_coef_h"): + pytest.skip("compiled core built without mokken_coef_h") + + ok = np.array([[0, 1], [1, 0], [1, 1], [0, 0], [1, 0], [0, 1]], dtype=float) + with pytest.raises(ValueError, match="complete"): + bad = ok.copy() + bad[0, 0] = np.nan + mokken_analysis(bad) + with pytest.raises(ValueError, match="integer"): + mokken_analysis(ok + 0.5) + with pytest.raises(ValueError, match="non-negative"): + mokken_analysis(ok - 1.0) + with pytest.raises(ValueError, match="2-D"): + mokken_analysis(ok.reshape(-1)) + with pytest.raises(ValueError, match="lower_bound"): + mokken_analysis(ok, lower_bound=1.0) + with pytest.raises(ValueError, match="alpha"): + mokken_analysis(ok, alpha=0.0) + # crate-side guard surfaces as ValueError too (zero-variance item) + cst = ok.copy() + cst[:, 0] = 1.0 + with pytest.raises(ValueError, match="zero variance"): + mokken_analysis(cst) \ No newline at end of file diff --git a/tests/unit/mokken_tests.rs b/tests/unit/mokken_tests.rs new file mode 100644 index 000000000..f93aeef35 --- /dev/null +++ b/tests/unit/mokken_tests.rs @@ -0,0 +1,391 @@ +//! Tests for Mokken scale analysis (`mlsirm_core::mokken`). +//! +//! Every assert reads values returned by the crate (`MokkenH` fields or the +//! `aisp` label vector). Each test names the crate value it reads and a +//! mutant it kills. + +use super::{aisp, coef_h, normal_upper_quantile}; + +/// Reads: `normal_upper_quantile` directly against published anchors +/// Phi^-1(0.95) = 1.6448536..., Phi^-1(0.999) = 3.0902323... . +/// Kills: sign/branch flips in the Acklam approximation (one such flip was +/// caught by `aisp_z_gate_blocks_insignificant_pair` during development). +#[test] +fn normal_quantile_matches_published_anchors() { + assert!((normal_upper_quantile(0.05) - 1.6448536269514722).abs() < 1e-8); + assert!((normal_upper_quantile(0.001) - 3.090232306167813).abs() < 1e-8); + assert!((normal_upper_quantile(0.5)).abs() < 1e-8); +} + +/// Deterministic xorshift for simulation without external deps. +struct Rng(u64); +impl Rng { + fn new(seed: u64) -> Self { + Rng(seed.max(1)) + } + fn next_f64(&mut self) -> f64 { + let mut x = self.0; + x ^= x << 13; + x ^= x >> 7; + x ^= x << 17; + self.0 = x; + (x >> 11) as f64 / (1u64 << 53) as f64 + } + /// Standard normal via Box-Muller. + fn next_normal(&mut self) -> f64 { + let u1 = self.next_f64().max(1e-12); + let u2 = self.next_f64(); + (-2.0 * u1.ln()).sqrt() * (std::f64::consts::TAU * u2).cos() + } +} + +/// Simulate Rasch data: P(X=1) = logistic(theta - b). +fn simulate_rasch(rng: &mut Rng, n: usize, bs: &[f64], theta_scale: f64) -> Vec { + let j = bs.len(); + let mut x = vec![0i64; n * j]; + for p in 0..n { + let th = rng.next_normal() * theta_scale; + for (i, &b) in bs.iter().enumerate() { + let pr = 1.0 / (1.0 + (-(th - b)).exp()); + x[p * j + i] = if rng.next_f64() < pr { 1 } else { 0 }; + } + } + x +} + +/// Brute-force oracle: two-pass covariance and sorted-column covariance, +/// computed with a DIFFERENT code path (per-pair, f64 accumulation in a +/// different order) than the crate's matrix construction. +fn oracle_pair(x: &[i64], n: usize, j: usize, a: usize, b: usize) -> (f64, f64) { + let col = |it: usize| -> Vec { (0..n).map(|p| x[p * j + it] as f64).collect() }; + let (ca, cb) = (col(a), col(b)); + let (ma, mb) = ( + ca.iter().sum::() / n as f64, + cb.iter().sum::() / n as f64, + ); + let cov = (0..n).map(|p| (ca[p] - ma) * (cb[p] - mb)).sum::() / (n as f64 - 1.0); + let mut sa = ca.clone(); + let mut sb = cb.clone(); + sa.sort_by(|u, v| u.partial_cmp(v).unwrap()); + sb.sort_by(|u, v| u.partial_cmp(v).unwrap()); + let cmx = (0..n).map(|p| (sa[p] - ma) * (sb[p] - mb)).sum::() / (n as f64 - 1.0); + (cov, cmx) +} + +/// Reads: `MokkenH::hij`, `hi`, `h` on random polytomous data, compared to a +/// brute-force oracle computed by an independent path. +/// Kills: any algebra mutant in `pairwise`/`h_subset` (wrong denominator, +/// mean subtraction, sorted-column pairing, aggregation order). +#[test] +fn coefficients_match_brute_force_oracle() { + let mut rng = Rng::new(42); + let (n, j) = (60, 4); + // polytomous 0..=3 with item-varying marginals + let mut x = vec![0i64; n * j]; + for p in 0..n { + let th = rng.next_normal(); + for i in 0..j { + let mut score = 0i64; + for k in 0..3 { + let cut = -1.0 + i as f64 * 0.4 + k as f64 * 0.8; + if th + rng.next_normal() * 0.7 > cut { + score += 1; + } + } + x[p * j + i] = score; + } + } + let res = coef_h(&x, n, j).expect("fit"); + let mut num_tot = 0.0; + let mut den_tot = 0.0; + for a in 0..j { + let mut num_i = 0.0; + let mut den_i = 0.0; + for b in 0..j { + if a == b { + assert!(res.hij[a * j + b].is_nan()); + continue; + } + let (cov, cmx) = oracle_pair(&x, n, j, a, b); + assert!( + (res.hij[a * j + b] - cov / cmx).abs() < 1e-12, + "Hij[{a},{b}] crate {} oracle {}", + res.hij[a * j + b], + cov / cmx + ); + num_i += cov; + den_i += cmx; + if b > a { + num_tot += cov; + den_tot += cmx; + } + } + assert!((res.hi[a] - num_i / den_i).abs() < 1e-12, "Hi[{a}]"); + } + assert!((res.h - num_tot / den_tot).abs() < 1e-12, "H"); +} + +/// Reads: `MokkenH::h` and `hij` on a perfect Guttman scalogram. +/// Kills: covmax mutants — any error in the sorted-column max covariance +/// breaks the exact H = 1 identity (for a nested dichotomous scalogram every +/// observed pair is already comonotone, so S_ij = Smax_ij). +#[test] +fn perfect_guttman_scalogram_has_h_one() { + // 5 persons x 3 items, nested pattern + let x = vec![ + 0, 0, 0, // + 1, 0, 0, // + 1, 1, 0, // + 1, 1, 1, // + 1, 1, 1, // + ]; + let res = coef_h(&x, 5, 3).expect("fit"); + assert!((res.h - 1.0).abs() < 1e-12, "H = {}", res.h); + for a in 0..3 { + for b in 0..3 { + if a != b { + assert!((res.hij[a * 3 + b] - 1.0).abs() < 1e-12); + } + } + } +} + +/// Reads: `MokkenH::zij`, `z` on a hand-computed 2-item fixture +/// (X = [0,0,0,1,1,1], Y = [0,0,1,0,1,1], N = 6); the exact hand derivation +/// is in the test body. +/// Kills: sqrt(N-1) and variance-product mutants in the Z formula. +#[test] +fn z_statistic_matches_hand_computation() { + let x = vec![ + 0, 0, // + 0, 0, // + 0, 1, // + 1, 0, // + 1, 1, // + 1, 1, // + ]; + let res = coef_h(&x, 6, 2).expect("fit"); + // hand: means .5/.5; centered cross products: + // (-.5)(-.5)*2 + (-.5)(.5) + (.5)(-.5) + (.5)(.5)*2 = .5 - .5 + .5 = 0.5 + // S_xy = 0.5/5 = 0.1 ; s_xx = s_yy = (6*.25)/5 = 0.3 + // Smax: sorted-sorted = comonotone = 1.5/5 = 0.3 -> Hij = 1/3 + // Zij = 0.1*sqrt(5)/sqrt(0.09) = 0.1*2.23606.../0.3 = 0.745355... + let expect_z = 0.1 * 5f64.sqrt() / 0.3; + assert!((res.hij[1] - 1.0 / 3.0).abs() < 1e-12, "Hij = {}", res.hij[1]); + assert!((res.zij[1] - expect_z).abs() < 1e-12, "Zij = {}", res.zij[1]); + // total Z for 2 items equals Zij + assert!((res.z - expect_z).abs() < 1e-12); +} + +/// Reads: `aisp` labels on the same fixture: Hij = 1/3 > c = 0.3 but +/// Zij ~ 0.745 < z_crit(0.05) = 1.645, so NO scale may form. +/// Kills: deleting the start-pair Z significance gate (mutant seeds a scale +/// because Hij exceeds c). +#[test] +fn aisp_z_gate_blocks_insignificant_pair() { + let x = vec![ + 0, 0, // + 0, 0, // + 0, 1, // + 1, 0, // + 1, 1, // + 1, 1, // + ]; + let labels = aisp(&x, 6, 2, 0.3, 0.05).expect("aisp"); + assert_eq!(labels, vec![0, 0], "Z-gate must block the scale"); +} + +/// Reads: `aisp` labels plus `MokkenH::hij`/`hi` on a hand-constructed 80x3 +/// contingency design (profile counts: 19x(1,1,1), 11x(1,1,0), 10x(1,0,0), +/// 10x(0,1,1), 11x(0,0,1), 19x(0,0,0); all marginals 0.5). By construction +/// H01 = 0.5 (start pair), H12 = 0.45, H02 = -0.05, and candidate item 2 +/// passes every other gate at c = 0.15: Hi(2) = 0.2 >= c, Zi(2) ~ 2.51 > +/// z_crit, augmented H = 0.3 >= c. Only the negative-Hij (Criterion 1) +/// exclusion keeps it out. +/// Kills: removing the `hij >= 0` candidate filter in the add loop (the +/// mutant then admits item 2, flipping labels to [1,1,1]). +#[test] +fn aisp_excludes_candidate_with_negative_hij() { + let profiles: [([i64; 3], usize); 6] = [ + ([1, 1, 1], 19), + ([1, 1, 0], 11), + ([1, 0, 0], 10), + ([0, 1, 1], 10), + ([0, 0, 1], 11), + ([0, 0, 0], 19), + ]; + let mut x = Vec::with_capacity(80 * 3); + for (row, count) in profiles { + for _ in 0..count { + x.extend_from_slice(&row); + } + } + let res = coef_h(&x, 80, 3).expect("fit"); + // verify the construction via crate values: item 2 negative with item 0 + // yet passes the Hi gate at c = 0.15 + assert!(res.hij[2 * 3] < 0.0, "Hij(2,0) = {}", res.hij[2 * 3]); + assert!((res.hij[2 * 3] - (-0.05)).abs() < 1e-12); + assert!((res.hi[2] - 0.2).abs() < 1e-12, "Hi(2) = {}", res.hi[2]); + assert!((res.hij[1] - 0.5).abs() < 1e-12, "Hij(0,1) = {}", res.hij[1]); + let labels = aisp(&x, 80, 3, 0.15, 0.05).expect("aisp"); + assert_eq!(labels, vec![1, 1, 0], "Criterion 1 must exclude item 2"); +} + +/// Reads: `aisp` labels on a two-cluster simulation (two independent Rasch +/// dimensions). AISP at c = 0.3 must recover the two clusters exactly. +/// Kills: selection-logic mutants (wrong argmax, wrong exclusion of previous +/// scales, missing multi-scale restart). +#[test] +fn aisp_recovers_two_clusters() { + let mut rng = Rng::new(2013); + let n = 1500; + let bs = [-0.8, -0.3, 0.3, 0.8]; + // cluster A: items 0..=3 driven by theta1; cluster B: items 4..=7 by theta2 + let j = 8; + let mut x = vec![0i64; n * j]; + for p in 0..n { + let t1 = rng.next_normal() * 1.6; + let t2 = rng.next_normal() * 1.6; + for (i, &b) in bs.iter().enumerate() { + let pr1 = 1.0 / (1.0 + (-(t1 - b)).exp()); + let pr2 = 1.0 / (1.0 + (-(t2 - b)).exp()); + x[p * j + i] = if rng.next_f64() < pr1 { 1 } else { 0 }; + x[p * j + 4 + i] = if rng.next_f64() < pr2 { 1 } else { 0 }; + } + } + let labels = aisp(&x, n, j, 0.3, 0.05).expect("aisp"); + let first = labels[0]; + let second = labels[4]; + assert!(first > 0 && second > 0 && first != second, "labels = {labels:?}"); + assert!(labels[..4].iter().all(|&l| l == first), "{labels:?}"); + assert!(labels[4..].iter().all(|&l| l == second), "{labels:?}"); +} + +/// Reads: `MokkenH` fields for score-translation invariance: adding a +/// constant to every score of an item must not change any coefficient +/// (covariances are translation-invariant). +/// Kills: accidental use of raw (uncentered) moments. +#[test] +fn coefficients_invariant_to_score_translation() { + let mut rng = Rng::new(99); + let n = 200; + let x = simulate_rasch(&mut rng, n, &[-0.5, 0.0, 0.5], 1.3); + let mut shifted = x.clone(); + for p in 0..n { + shifted[p * 3 + 1] += 3; // item 1 scored 3..4 instead of 0..1 + } + let a = coef_h(&x, n, 3).expect("fit"); + let b = coef_h(&shifted, n, 3).expect("fit"); + assert!((a.h - b.h).abs() < 1e-12); + for i in 0..3 { + assert!((a.hi[i] - b.hi[i]).abs() < 1e-12); + assert!((a.zi[i] - b.zi[i]).abs() < 1e-12); + } +} + +/// Reads: error `Result`s from both entry points. +/// Kills: deletion of the validation guards. +#[test] +fn rejects_bad_inputs() { + let ok = vec![0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1]; + assert!(coef_h(&ok, 2, 2).is_err(), "n_persons < 3"); + assert!(coef_h(&ok[..4], 4, 2).is_err(), "length mismatch"); + assert!(coef_h(&[0, -1, 1, 0, 1, 1], 3, 2).is_err(), "negative score"); + assert!(coef_h(&[1, 0, 1, 1, 1, 0], 3, 2).is_err(), "zero variance item 0"); + assert!(coef_h(&ok, 6, 1).is_err(), "single item"); + assert!(aisp(&ok, 6, 2, 1.2, 0.05).is_err(), "c out of range"); + assert!(aisp(&ok, 6, 2, 0.3, 0.0).is_err(), "alpha out of range"); +} + +/// Reads: `aisp` labels on an exact-tie design: X0 == X3 and X1 == X2 +/// (identical columns, Hij = 1 for both pairs) with the two blocks exactly +/// uncorrelated (balanced half-split vs alternating pattern gives sample +/// cov = 0). mokken's eps tie-break (search.normal.R: penalty row*1e-10 on +/// the LARGER member index) must pick pair {1,2} first, so labels are +/// [2, 1, 1, 2]. +/// Kills: reverting to first-encountered lexicographic tie-breaking, which +/// would start with pair {0,3} and yield [1, 2, 2, 1]. +#[test] +fn aisp_tie_break_matches_mokken_eps_rule() { + let n = 40; + let j = 4; + let mut x = vec![0i64; n * j]; + for p in 0..n { + let a = if p < 20 { 1 } else { 0 }; // half-split + let b = (p % 2) as i64; // alternating; sample cov(a, b) = 0 exactly + x[p * j] = a; + x[p * j + 1] = b; + x[p * j + 2] = b; + x[p * j + 3] = a; + } + let labels = aisp(&x, n, j, 0.3, 0.05).expect("aisp"); + assert_eq!(labels, vec![2, 1, 1, 2], "eps tie-break must favor pair {{1,2}}"); +} + +/// Reads: `aisp` labels; independent items (no common trait) must all remain +/// unscalable at c = 0.3. Smoke check of overall gating (not attributed to a +/// single mutant; the Z-gate kill lives in `aisp_z_gate_blocks_insignificant_pair`). +#[test] +fn aisp_leaves_independent_items_unscaled() { + let mut rng = Rng::new(5); + let n = 500; + let j = 5; + let mut x = vec![0i64; n * j]; + for v in x.iter_mut() { + *v = if rng.next_f64() < 0.5 { 1 } else { 0 }; + } + let labels = aisp(&x, n, j, 0.3, 0.05).expect("aisp"); + assert_eq!(labels, vec![0; j], "labels = {labels:?}"); +} + +/// Monte Carlo: >= 500 replications of a unidimensional Rasch scale +/// (normal and skew-positive traits). Reads crate `h` and `aisp` labels. +/// Asserts distributional behavior: mean H within a plausible band and +/// one-scale full recovery in >= 95% of replications. +/// Limitations stated: this cannot pin exact constants; the algebra anchors +/// live in `coefficients_match_brute_force_oracle` and +/// `z_statistic_matches_hand_computation`. +#[test] +#[ignore] +fn monte_carlo_unidimensional_recovery() { + let bs = [-1.0, -0.5, 0.0, 0.5, 1.0]; + let n = 500; + for (label, skew) in [("normal", false), ("skew", true)] { + let mut full = 0usize; + let mut h_sum = 0.0; + let reps = 500; + for rep in 0..reps { + let mut rng = Rng::new(1000 + rep as u64); + let j = bs.len(); + let mut x = vec![0i64; n * j]; + for p in 0..n { + let mut th = rng.next_normal(); + if skew { + // half-normal shifted: skewed positive trait + th = th.abs() * 1.2 - 0.9; + } + th *= 1.5; + for (i, &b) in bs.iter().enumerate() { + let pr = 1.0 / (1.0 + (-(th - b)).exp()); + x[p * j + i] = if rng.next_f64() < pr { 1 } else { 0 }; + } + } + let res = coef_h(&x, n, bs.len()).expect("fit"); + h_sum += res.h; + let labels = aisp(&x, n, bs.len(), 0.3, 0.05).expect("aisp"); + if labels.iter().all(|&l| l == 1) { + full += 1; + } + } + let mean_h = h_sum / reps as f64; + assert!( + mean_h > 0.35 && mean_h < 0.75, + "{label}: mean H = {mean_h}" + ); + assert!( + full as f64 / reps as f64 >= 0.95, + "{label}: full-recovery rate = {}", + full as f64 / reps as f64 + ); + } +} From a7ef4db49b855cdbc03ac6dd8d56435e1b564bb4 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Sat, 25 Jul 2026 18:21:11 +0900 Subject: [PATCH 2/2] chore: remove temp patch file --- mokken.patch | Bin 93802 -> 0 bytes 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 mokken.patch diff --git a/mokken.patch b/mokken.patch deleted file mode 100644 index e7c2883cbd1ef8a2c951123eac3d7b5be4915ce8..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 93802 zcmeI5`RD59ncz?Kc03k~X=p$K@zBr+fsp*C z@9zD?x4wNNA~UP9XvEGQ`*7&0%8ZP-ao_jGjmZD=zi%z>F7Cv?t;J6K`!;?(iSIWS zCl?yh zM)1Skzy`ecJgzg3+l4Fd7MQ^AlhWH>@Z{aay_h*r>sNl?i)&yW9@V&aW1NJp#=zF% zn`i|kZp43R#NYk+6h73;4{`T#xx?K@(KmPk&f76heFwHjakcH;a}ySRD$jg(2v@>% zqu>gx0f%v~{(kz}kE9%lx;s_@Khom5^q*@?`V{mY*$u!;gvQy~Dd)zY}9Z^_(YM&s?Dv z90it`1+#(&!QW2I?NQ+D$7uO5V3Mr-Sg0q-db@(l@^P}jE$M(~nXUMlUylM7;Sy;9 zUZ@K;8Ksv8=(!VDQ`h~mzdhA;ms` zGhj;10zDqZKP*H-$0&E=s$^ce3E5>7q-;NW5MAzubS0G^7A#vgBj-THXkhAg!5jF$ zRpyILW(2SdM~Hgc1!m-1+zF>b-Gojw1n&DG*4;R>eH^GN3mXL)n_pmSpkjD9sD@SmmqS#U=e&uzL~RlK^&N{aUZIKBXDvM zSZ55{LGxY$pqL~LZj$zf4}T0ez6Z7GD$vMi#$TJ zTT96BD&(hNbM6PEKNQZQtt=ihk##agWoe){5-L9bsLYdJn(xEF=ILe633uW#VMAXF z?{|yc+c~rr&=;N%uUX=mF}O6pr&PdE%#fLEE&Jr#cFYLgOwF}Z+5o97tmm3=JqfBm zidlGc%?)2p9)qM^{LQSNA96B$CXU_h6Bm>=tPkT_j>7U|x`EW=VI{6@IC_fWVeyjM?m*U>TfOaSB2N;v( z-i`0*XZ)T^F`l>_++B)xyYUyRDc#6yv1Ql<(cAtf5cN8Ay$=r!@a(Y{g6ANx^bcju zTQRCV_|=|!|1@?H*`QDGjGW`UzZ3t_m!HISVh~ws`77rOzObyq9#|Iydoilgn>#`J0Pt(F&a;o`ln`2VI5pTg9u<_x<>OZ+ZPreE+_1!{2WCNn;qek;~d* zViiOcpqh?9nXP!~AfN;*XqKIjfsC7W0=nmM_aNHt$ET=vB7PxR#HmN^1sItH5EIXy z3t#YJd_NoiZU;}ESe%bB-iiBAoWgG5R)hi&CAZ!La3l+*QdttFIMy7&6}Gz@hs>4vvLtFOZ^QsZ2MFpFlyG7m?p z+k!#017~0aAM$$9k3E*!=jyhbU18heFKHPw)QGH%ch$yc^U zeYY_WI!jst?S{Rz_i&|$mq)Y>Q?UziLPQTNZw@Aw~zLZ1|9@D=fJ z_9DVP9{-iYx*X$QUiM2g3-{oO$fw&2nY~+Dwlvfz%nfSLmwiCxuxxu33o+(*r5C>a z9RE+mlUHIK;}mOcACqs+K~T6B-=haYzdc=d0}2e*dGSIuzAk|6M7mb@0M))u%F?afB=1bH#h~)Q`7}M z{E$`|)1z z4UWFfu!+65=ZSS|VYBsEy;1?x_lq6GCi4vWa^F251C_SchN~?cqM_@5VxZRQi^WXdfJ)to#>F!mDy6}Rz8L-4|91O+8RAUyS>m| zd!;?POfwvZtoP9Hb;+9GAKi@aK!<;^9aqq{%p~zf&h=(o!G^N30o;29*FOY~Zxk4e zU-CT<$_$`I4~=Z3vQK1woEK$8<>Bn(*}A?y>&vY6<3Bb*Gvx}i&@A)Lsfgj%*=qB0 zZ>Rs7ZLclJe&O6Z#p+uUeJ}9=ceia;+3r6Dh4;gT{}>!}X2Lp8S>nXy(XFtk+)@^q zKZ*?&m$zD#Yg@$^RosL(O|^tP^qsKkmkR#SQX&X#NH~> z&+O0{tmpOd*ZgRU>?HQWvq{}{EP6IBjzuf>#`tJgbjhi`N~ECl(=JE1X7#OJxtXQ#rJoQ^A-taoe*$Mj0a zraf6hw-ta?1LhhBYutFnvR8b5S2UY-+3>zC`|1|9@6)ZA&F#b3d^k>KF{b*!ODB$? zfJe~U+#X~ODsANQc z7xym()x9RvuCWq{;n}kbGj2!wJS0?Ai>~?7 z3b^8IB0c2--3ooISe$RJK}D0YiblrNYjC*=fd8#5e_!{MKC4z!pvQRa6Pn-AbX{wi z6zijTuan?mQIqTX*Lwc-cf4%=F)!+Y+Iycx53JP_2fI>H@kd*&ID87L&r_-dydEAT z{rs+6yBk@jzQ0xNh{s$P*^eaeq8;gO^!|SQzqfLnzSZ<7!caYO{auw&c>M4Pm706; z+1BT@w1Y;&dN5^ow1Onwfp(ceAp=M3aXE70tWChBWNKPi7{8v0sGt$E9&~&0=kU(& zMlW!ft3-A6UR-^1WDe?;+Cy?cl8=5CEiA|at3ir%mFFU!QgssXlX+WpMj1~X#A*e) zxQ&#>^#K}C zp5SSqaP-~dj*3M49HzPO-NKVp&Y5|fF}|&f-!Ntm&ZK-5D_LYm{;kxkeG;RHPSW^UZ>g$MVrx@M-auN(VTQQH zowy6{3ErU}Wpr3uj6t-K*0spm$^xw=E0J<>;U?P+j7y5L^X1iSHVC%&Uyuv+(vl&KCYF$`qz zBb-leN_LR?Dl$}-&z=M2WhtL^tkYVO0wh{vDc_8CfkTw#_6`!{)jLSwt)9bv-?JMYFUN1uH`-0}qS=M!pUj zqTg6aAM4L>pv_x*Mc-%)a$w?|OJP?&jCN>T*&Q?u)Yq>THuXlOD?D6tFxuC-Cq^}0 zto6v6)qU09Td(!+i-di12orh(=v;F^ji~%>c$u|fbUqkiO&Xes%iuY2%1&J8oA9A( zGgf=_P3sd#MqSO&!e<$&eRSbJ_nHXq)Ck`cTZ@kH3igSREZOa2QIUPDXo~k%U`753 zkPkh$Z(>eE5AJxu+Fyc#PePi~qad?<5b~t!tZ?6rf2xP!*S(?}+Bq|Vy*(@`oVv>2 z1A}wNpXzHaWa%B&kOA+|KYSUu+Meb?wa%Ng)H;$u->qKj1aj@@fcyEch(z{hlig-$ zM`9mIIb6`v+M2`~%^r~Tw)Z&U^PqEqgBo5xht#<`re66Y(}ndyM!C-%uz<{cZ}Ce= z0bKk`{Qg6UAifC>+b?p_!k0UyJR*J>yCw7-&tNkor`Q(k12T@D!kXC&Ah!*)mT2}| zv9lp@s(aFJ+e?oeeaBPiVRyaPuT01kner&b$Sq_N$s|US)@Ftnk2H^31SJX!fKqEA{X_uukEjg-x7= zCNW=mmaKUoJ&ZxjmT&;GeRQwti@v$554zupr}578ylJmBrLRlBEi7YFrF};9Mn7Ok zc5m*kJ5@T=&;z^oZP`ECE7hP@~QF5gf)Oy%D&I`R#75VKmJNxIPvS_?$ts9;Qb z!rSojlhQ)ZwNUJpUZFXCoCu!7_v0D%iSEWsq~-cukF2n{-k`5|{avSB{uFRIhG5pV z0YIvHOe*ov?v@O&3755bBKg`}V<)D2QKaN*5IDH62dO1oBeXcU6i+ygQk1~lkV-tw zemwBOeR1rVcLwx*?!kI@OfNDY&%K4L*0Zjx;hHmgmgrvpp$VxCfimcC@Xvhs3nh^? z##|r6@VrK8A-fXzC$cjSN;APT;3xMM*AaL|vrs;xp5%ktM}}(n&(=etGgQ|eYt7~R zG80#!XU5W)Nk0Y40grRXUU*_gx~z_Ki6N1wm)p-hRH+Sm2H=-b_po2);t<9G0+A=L`LTBQe+52$4YO8cyZ zK8w}+IarIM@@ zFWIK?Bs%?hd?E#$+5m*sJLIcn*Wk&dk1X5exYEMKi3N;H|LoVRpK)dYU$NGCshyBW zIK0-mz{*}=>Qj`%tDpTT?lQ-oTeYT@+>ds&Jsn1|+2e4RCE8OJ)~h1j-N2vr#Y0gz z7^pa_fqd3htm=SK@;$eL=Qv@Y#oJEY8`nK>?H5Ihd|3WwHf4-V`pn98RT=NuuK*Wk4=Jo@;!$Ga7z_N znM-ljX4%s{o_AwF`5F&EDDD>hYi$h8ju%;!ZUfUq zZnleKXRSyxt-|f`(RqIya?W`G_+XrmfVUw(K-%hMPTyf}hPCE)bV~Xfv`Sy2rHJID zwnVgsyZIcmNnHQDK>lw__?5}ULJ6Or2j}DWC`*V9vR9@&B%1u;Qc7HHhrH}}l!RLx zHQBeXS88jTfpg?kr+}=+>t>ird=VF2S?apEYqm+BnkQw^(4DOot>qmY!7+&eVdS6j^CY$EBX=RWDuj`J?m;UpSRhd^mqY=Q*i&;v*s{YLivX zJly9*)C(=YD%h(t#+K@yx{4gwD<|BKjlIdpUk4_g@$i_?+CFF>(J}t2(HXm~-8O`t z>&w%w51ca}}<*~Pbd+>sfbV73_*SlPmI3A^M{;R!gI zd_39)O*dTWKuh5Vs}}E|?DcqAE%0SI0YMr=__ChRFWbSdV(#@ix-!3HynqIsq?H>o zV9ba7hbn*e${dO4IXz=9uG{n8D}BOOdqtar3+w@&8EuriQihppcu0GJA$*+{I`koH z4vgEL_2<<(U7gcA<8ht>5^hMPb^D(C-s;r>tpyCNy`mt^LVo{;AxZaRtjy=*okDve zUGfRgmYkN>s%^(*KfTt+)A+fs-+m6g#-S%U&&FDhKa9DexrcfLilBSsv-MK!=e_3N zYT>@Md|7zT=R(3OSX)<T=sdjmx_(Ay@74T=#aaD?v(97z#ci#o@;r2EnY>%h zZo;E6w*6^e9iIi9={lCq{KRT8SwGfs&}Qf_s+@X!^rO~u`fW#+HIg^LU+B-y(-Q&b zv?tQTtKFJCUWAqHWq(7CLI$D14@Im0P;~mU@I6_Hqf!G1@B+acbx$?_@vL|OSnxW~ zqdi`2MH61}m>}DaK3W^2DiSD|k(g{hI_2SH|HdK~NjR+gmZu{Kr*UdnkB1({UxuSB z6I1iDz3{XC%=uCi2(RaiwH}S#ScNh2;^xBZ#hPGwZxqdZvry_rv{2TNTnUi`+)gHO zJK95^+9#lvB-h(+$FDJ8!0?N*z(AzTX_(v5_cVPCsj|}f3@f%h?s}XVSE>1hDQ!&W1^P5()?+TFr!IpWP2vF%G|No6eA_u}8RCfW|3 z!vlF5+(YfpkMSS+@p~?8Q<`Mn^Jtso$DV@?vp5oieUxZN{jW)lgOGjX zueM-dBBkT$Kns|X{-h7&vF-zD+N=+0lGV__G3)rtp3-WVR`t~DT+dT+_87)A+xLq` zktSVlImfd-(r@5G@%FLM8bm7(0-G=5>aX#s8nvxbSMcI%mGyqs72o9xKT&B$YVX3Qu( z!?F2H!|lGL<-rU%7SGaW&)II@f1>x$PT-W&=z8Ai+gqFV*Zms}M`Q^k*d9S9Tbmlk zu_v~B-H5E9WRJGHNADLpX20tomCa-xnpJubo6Twpe@L}_Ny&=vC|~u)B#w~e;ZuTA z`iLHwDG{glfDC0Vsi!(#(GM7=s3U1N{-DuROKSnQ1)9EXj}_&x$S0^qmkK_tFrzvP z>kOf6xA%sTpQYYN)=%~`Z72SRYJHif^O~OhcFFppYotT8D?0TW*ygTwA>reI-7|q( zWyyQIbGoTwm9}L|i`9P0M*%T$zxE$uxjzdi=D+c1xLOMbuc|XZD$okv#f1)SeU!bx zPrs5hzRu73->^nNK8^GIz{7aZwf0_ErO>!u(U)&RTW8&0cUHKr+4!vR!j-@!z833~ zW9zGMd%q4W4(l=KtM;JDuMAhu(q8rcuJnebqRyMU=XQm{U+Rh_U@>$krfBj7Q1?!)&{*o`iXj-q_#1Fv!P>aLU z=rLwbW+$WTu}D&P(SnWy#(U`Xy&>*Dh4=A3#DCLJ0JZrJh=^H+J+*c`t$5Cs?CKm_ zY1n()HFdpjXIVi%@OFfZRWcg#>WB0EtZ&T8&AackXx|lWWefK6p$9vR^xlZsqZ{~? zx5Y2*C4Zx6v-{=#&;|_s3dQEF*2Fg@D*(UIGvD;dUvYrH>65oHthSP9VMX(-va(hx ze*yUKN3WdpAnmsmUJ&p4lFtjC@Q>t8BBSUK>QG#LWh+QT!87#FiZHFPZLD7sN3*Ky zy2jR8Qr}Ho1~sz}0y3=c^Z3OGd9r~r14M>Xb?oS{U#fA*p0Ywu4*hgg!?Ne%-GWJc zj^QpKX_gi?@W7bflZ55R)7Xib(F5E@wQN@Ukdyd1{xgPrGx8EdTGY?6?-*WX=RGkJ z{F}V?L-eNcWyxzAQK@1YT_4bh^Qw^$>2-Qm&69o)qaFE=Nl&12{^X;mi~KZx3ugA{ zvl79&nB-2oTfZuEhRU2#ko40I#c3*Oe%1k&!@oo?SWm5Yc*I^S;%Ta0lecu}TuWB} zsbrxCu<3l=8zG|Y3u<0PsJ0rOjNl~WYjxw{YiV=98*jhZ*!Bg1=f&c&5;(F51vTzcS)hzWHc7^@^#a;MRe81w7tcp4ldPr}NOUlR^w_k4^G zWx+6dLnrH@KR97;T_5M72ihX<%mjJkKi?()>v{kN#)PiCAq3u%+`}i#WU7bN)ctrs zUFv$J>!);T(iXh~Zy_hhi{=EUt)?wHSu)DGgzS6adHR7Pn45Z-=6CRc+rT?@?z)4679t*^@KHUv9OvpmV*H z&T`U;89uFA{s| zZnrl@@P44hc5oy-b1J?QXZA<-==hdFR)dKYXVcq^2B{ z&@=WREiZHF2nXQpB3nb5GA2xYyk^FPLZV17g?_@6CDy5FOApYNt2M^z>ovye>ovye z>oslb>!Ph8ZEXtPotz2YX2p8 zY;$?tIh0qZ``4GdPLp>4Q9#F)K+ryx8M*V@5L;g z8)018;}*iYqWylfvCo0#)=qdl3gIkI?cp3@Rxvo2h#gz@ch1K*VaB#!^T9(QW1)A? zB{kc%5ue%L?|d2Cd!|`fBM!59&^0rWy@%7T1FUX)^!Lr|2Cm3BX>Z$~;#d8y8Mx|o zufxbLlbMm93Qt=%C_}4kA+$V*R&NBvcqCIA9sg}?H{BEu??+!$DBuBmUyb+%TLwnR zz#V1xS(Q!rYyYXOa*!` z$8c7AjjZ$7y2{J01ep@`ircGz1b_9nTn%cqr)S}sAtaFsI_=+tqli5BBa@4C>K$3w z7k0b%K6T%hts54Ep7Qw~hDYHYkjI|VF?CO>@Hjmf-S)sx_>4rjAoVUbj+27N*Qc<1 zXliVh_7vVNwca&WsJFRRW~4b`(TLGm3tPYEM7qha)%AJGho4!21<#OZp9+)jCvjk3 z&QJJLKHb=fd#oz#$3Ln|uvS=Xc{Skg```)Ys9F`INKr5t;B-CF_xIry02%bw8bwVv zWjt6F==r9FPGc*tp_~U*3aTk3S3>W@v-DK?Vfb{~57DUWaP)pCBn#@oYrIjAUg30n zGOw~F%;~EbLVudazCWD1XWmk@u9lpdwfnLJf@?vBeH4!o5oyI}T`xZuN?9+~eL+X{ zfhL@iw0?#<2T=XhR!(uK=8wkO4jzFIh}!H;D56!wt{9|-XSn9$vHiKpDx2W!V;lGM z>l&U*-b@P$zO;(GwX9r%s~?0EXr0Xz&nMs_@9|^B^PAw?;W@EH179CTE~2pSo}Aevk|Cz3QBf2i3^EvOdKubXF*`ru)k4 zt%A^A{}&Ot?Z)XqpTu7*kf}O~ zX=578vfhY?mbhFvbGi=J^*K}}t8Rgs2>eEQ-{$9brxg43HmgfEjy&EQLH%`O90YG` z*A}ojU$H9c@F@P2U6n89x@+w>B$_9Vp*OsI(`=mvbkF`&9vF|Qo4P-B-Q3ker^x+gXhT6r94+cW z8~iN1nOq%s5x>%prA}&UKC*v~Yh{Z~_0K{Y??q46HRQ^$2d&-Fi2(Rid7cm41de@D z+c-79E7~0n+K>P1sG0Mbn#($z+9uR_RpnfPH6yVbOFV<`VI1bUv{4JsT#H`lDSg~2 zNZV&M&fI53F?-zPfr#$L<;{23W7m8)^!?j&7}92r#dzaT7-RlC5ghe0k>F8#@9Ey# zVm$Esq0Kn$=3wbt+D^!(a}~K0cPqn}I|@C&sn1O!S@$8mH%1kp=PeX?ariZya{+AY z}jvC@5h?Jx|5T7 z)URSVRfvBQe)jcPC%GJTf&BeEB8z{Y)THaj$dw&1WII2lBix_0FPUQD)c#&xsj)nt z1B_K;L+?Vc9?G>Ro$xU43=2c;PxaZJDWO@crH&w5>!bJAf6MYg3J3fBx29 zc*YTPW*}*wtIuBV9N@0xKzq;f#H39&g?p%r+A-GY3D4i3*eNZ$S-pbtGjltGSVeF{8JmJ`@j9-}9a8 zWY5TEP$fd%+9w6xk1Kj31N9(KipT=_?^}7TVC*h_1KwEnx?A`r&;R_i^a>wf!-s9b z31cUe-^O>~9f!_$Tc~suXgK}-L3o|!K1Qvtjg9w1NYttkeZL=5*EOoX_v6-Aha)2i z)UhbvJwD&IO$$MN4=6sZ`fDK>-Vra>c(3KZ$9{c}oiy58m2bEZT&r6F;icdeuNH_? zT&trVrtS}MhfbQNneVpv(pBKtEl;;Ns;>>lgTwV!)*d%^3a@c$o?E0{hWlC=##;$D ze*V+IlHPswq~qjEhhq%ew0(RYT!m-wRgBN73i(W~fCIyny{&3J#09bblNH#{`|aF) z9Z0`e_R-Ryy~@beajN9EmgCpO(0T7fPile1P!@u3@-?m&t7vPEEv#wmGllE>re!Tz zW9^3W;I%?N{-KpFFDlueHbKwtx6E z+(VT=r&E9{Vqx;HcrSeGt$ghxukV6} zA5*Gst5sI{;?n02qo=xk&c=0`$Qy68ZRc?*D9BsXzl!zbFCriKb$nim|Nk>ql&(iT z{D*Oc%;A@zr$3KxSDH~Y;y(q(iEru|>Q2sh+RyO&zpO`BNzLXGupRi#b#hsw=OxBA zAE$i5T~LTz13cpuea*%2?M1z`^5}d6pU~3$%r#3D*E#h(f9nb54UiP>>fAc7S@}eh zgsPo$3(K?yV2?WQjTyePg!Rm=fKsy%M23D?4L>zP_EGoHLRV*$GRGiEAL=btqVJX5 zCnx{xh4$tzxzPHT#zUH+mm$wwj5KGq|D>q0r(+SP({S3(!*7=*J-ytTyAlvcnD@0f2P;cO{1tgU6NfQ*mc?ZrM!D(|E-^laOMBo*ACA&Ixq zRkj4fwsTw3PVWPOi`srvSu(FD+RkmcYRu%I)xFg?f}Z4tYMKmv;(nI}a@W31WNO4d zFDG7rhd2*h6{5Gw>=ToU5?G5Qf70R=?X9;px0uLhIA_66kw?|}8JA+Lu?Tk)%LAYH z{u;F+=tm?0Z-F-QZtT^d)9TUTTfY^R@3P(fJMh|B&y#?fuGTCcHL%rI>eYCe=&H|R z#?scd2wdS#{^nOZr>uOGUZAeA#nkk>;)FMOd00&{Wms5gs(O}bcSRTMQ9DWtRmu>( z^e~Qj(P`V!g1N?MtvSBG&H7uG7;S4kY1<{N^w6!>x%T_l?UzH6zKwD5K^Tjyjr@?` zQcokb$up8#yiC>QP$hgNAmRNjx<=c!zLYChN*`dprgiyahn3beTr(gxP?8K`jMFSOQsNAqo^_7U5ee;%_W zL-0l5bSIuWD1Ue2t&%?mcaSYeeY8_>MBmJuH5>Q~&Y~TeM|SR^ORbIM&u|}Xkb~f9 zP7$+i&E6g9nv>xa&FNOqLKIkE^vnJHa+lCu&u5^NA}2dQjLHl zG83N&0mLoNlI}#>Xj*1RgeNZ-oXLtfAL|Up=i!}x9d}t*28&b}B`(G(1f6R+QKZjR z>nrvgBTx38;u7Z&TVM(C`dF(4 z-zUQ>p+~GRcp|@2w~=NctI5e8mVI-GvKiJY%8EdHo%E6VYisehTP;kzZkS`?fB;E@ zazjiy69G>F-}=}=5Vdga7XH9fg}bmA;A#Ck%q~ z;~rL@YAW8v(ntTs2=_;m1{JTv$Xz1E)bV0@TP^DqG%${DSmz4WkJzPkD5==k4bb8OFg@Elv2z}tg5 zqP95RexP?IXLXNuI`R%2GS;`_Nj;Z2Hck}uNG&#>hgKxBg53Lj?K86^Ws2hPa+@t2 zVgjG|OvExD$LDG<7%F+@DIuQ2u+kYl(4OHUOfpZc4KiMf?J>TBj_jqFk~dbA=4q6p z;0X=Z#tlG3ytT?R=X{iLYMY<*R%E!$Z-6FKE`E`=)BI9_@$+T^Y>L0BoY2ekvlNmQS1zHa*RxM}&* z-0sgiLEkal$#!d{9BWM*t(Pe(dhrr-_%(VU9u}tQQ>QI>CEYOjv*$ZvS~MBBMUPoW zrDTD7=_hJj(H*KEz^bw1yzfzM)T^yObjystd=FcG-lc(7FqMZoLs{ASzI|}sRnd=M zXE>PI6bw^mRvA-AgNhlYexJ7geY>uuPOM_O%9FYOhmJ1$R*JrBb9@J)~PK|?Qbu`odI(@<({aXp-BHy_eh ziAlzu2jP79rZD+)j8t=eV)EzUr&X9t$(@2VTa0s~bcD2351-@*PinmW9u12Aqu$@& zc~=g#Qd%Ls6Q5{WorOnyp0)RWV$-MLe4r+$iZ13viW{MCZ-)Q%Zj8YB(KMdp#)FXkw?@PU-i!$6LeX?OIg)c& zWy#2Qoho`N+l#YVU%<|c!Hv&yv-F3~lb!w`C=2(>lW`UFJ< zsFQJDYY+HsnWt$n3MTB^=>2Q$efu5x{ugB-TBv>sOu9<*ac~8>K-TW0i_k>Son0Hi zO`W-XNJIE2)*tABs4?SQuY+<9gSZy&8O?*NJqfxJJ8&8lnGjZHTBt;2ajyC$GlhqU zCk9uGNBBghSz%?3wB{GqA@GFpPtc}bnNKKL3u2~3mfUw6mKszFjh%_w8dE;Q&El0tXk~%GO`_dq1L*6FRcy2D&yfvMznTG z55PVpr>lEpU7jKV=5gxT4 ztURxV#vFjXK+{vZLVtpeN}{WAjd4X6L)ucJ#aZ2FW+D!XuHx+vx{oi9pI47XKEUv{ zBOo8P3ii0Mrt-b9VPyS*LnmytRtwKSlHxvUd(w{dIQldv2~sRE6DFbhg$~juG2>U;BYuo}BnuB#iJ}wt zoK;q1Q}B^16<5k%%eNTdzY2@!#vkHo-cBG)Aa&qc9EXP{%d8y#c9D5@owrzg7SFs6 z2*K>jL8ur*Fs{>s!p(Pq1?)dHsTKwIk3l)&r36R;1gk2?S1P6 zTgDKy?LFRl(PF*TF7OB$V`}ro8&Gn1jdp&c>G)WyMzQp1onj2m_d<5B#TYsRPv_kk zp7l?ki}*F4z8Rc02gcvWco_-lTpP|{I|wWAApVh4R3rtYb8_|~4237%e=YD*1 zZO>dh)b_y`QC`e8f_V5=aL!iDZp=zabJcmsv{KlY77}aF)t@aSCnAFGp~<)RsDJ(4 zJho<}by*B-lcLe-IQrj;R`Yu+i6rsieHxG=Dk5Y0vIk0C09gRI#H%*k=bCEObcbJx z%ZTUe)@RCkplNQ|X)!MUb*EVTodP9V7;kq!;^@q7!$rBOKeh%(90~`@8n+URHpbH> zKY_P_#q(Nz8#%Yo;kOZ0I*xP{FYB-RJ@=DRLY-(GEoLt}(UT%a$Ed@o_^&-yT&?3< z{H=>a|2rO(rul1#2Zj@*Y6?=F&S%@@ngO9ffL|jqvfmVkVah7SSs|70;P(c(ht!M{F2fAwp`Xn*bySLKN!C9p8DC%Sxg)-l6z56xi8uDS?jFuUdt&GJ2zAg4;33yg$xA~C<^5UBE&3BBz|GL^c)Ot!ygqz~(LSyj2}J)}F<&CwaHF zpFJDx^yY8g(awK8%h`u54F6hihRlEytk<^z@yVEna;oATe0$#h0mm?3^CfGqWCY>C z2XR%a*;p*?*R9unITJ=27^=7ZvNndCwXh-UdV@7QC}~RGzaG}{dZ7TEO~pN42hdQ* zAgX9KR0FZb3oXdSv+4{jMO7+&XiYXB{fP|G(+@Et)kUPP;8an|k@Jst3O-t1VyCP| z6Gy6+MNteo4j6h~QSC+Y3Ngy?yxCeC@rFop-I<|ZEf?N*ARe=KOeIzFnBH8^yW2Te zjyJ(8>#fx&AW|lswWoLDp3jYP{PZf|I~Aid5)ckK`BX>(?|{D;<7*{_to-YP5A4mu zNqb%`dJ--4ac~tp#HtE>hqgd-?Z>@mg_dX!cowd=9YivR^%RZRhH})zUM7Qs;4tGUB!&flw`bE%BJaTuE)1xT?oHIr;4~QuccjyJ5+HE zdD(p#3ipx~;96ahR&Z-$)xQ%G=1yf8bK$M1vPft(K_H%DRX}Z~AyOXIc|M&dg-&)) zXjkLPaQGAd>pzMH=iB>Bc{w|?{U6Gz#l@Me-TQR=8e-X%XOn$XrOgTn7-H8a_(Gl& zU+^z0XQHuZhjr@PuBMq~;Yv_U{)98OwKmJ@sXvBXJPUnrDMn_k1?e$I8}B_I)IGf# zZul`~E*vISv1l#dp&L{^1*W`4Kh%_l+r0MU;Y#W1-lzDke@{Ema<)Hf>C#xdvy^Hp zSr%4ok=!RmZr@wV1b!HvmwXz$g7?BN|1eIb(;6>OGDF(von7OUwv@IO{Ys<7@BnuU zhshi8ssoH?&$4F9EMctN&g-X0+55 z(naz&xr3*wHE*tAWz_!a5?{)LCM3L>)iDY!rD;<0{3#&j{nB@fmO|2c&1!yN#g2J1 zX8n$ap8jrcZHEl=oMHjJ6;`9M#)0n!<;UI*Yg_X}+?SW9*pX*?D0}V&0~+^nu^3bi zSTb5k1Q$Rsj=h;<_;xzXB|BWPT@!y5%!6~2l#F}u>hK56~{d)kKj+^lT(B= zqv^T<*9y2o2beUrDhc>CH6xGrEPP5ZE~;q1liy#<{IHxuoZ8BP13WT8`y#%fqkfJ5 z*gP}|w9`6x#zh|6pO)6C-_Mp(<#w%9|1sdBRtjAwI_b0@{8FM%Fp4a=Ruuo6m|8I& z(Po_i!sd|4z}nIS{f?i)(?98zVOtlk6$mzZw7bJSiB&k7RY!f%!y z8A_AkW=FiLTlSY)iWz91Jo{Ox0s#Z+m6f3}J(0ZuIqf;;2z;rx*6vzz8JSZ%6a`2V z>f{zI`TI-!|6w@?_z^y<=W8i)kN#vXueRr*ZuFeN%+te}uA-uFjRw6I?THw{7Fa@` za+m!c^p*$!Y@Y|1LYyh9Vnboq4 zAA6&tbIk=n{yMloGX8Wa-9YcvnQ_b6XC&L4FcM29fRi4fq|Xla_r190YAxAJYEQ9T z!?UEYJ9t6XOR8AM@+wbB&KQ1gEjwAjb(URv3{FZrt6ZII7#^Ck5L$Dw&4U}UueMUQ z+HFRXKFFfOvD6qTFX_mfyJSh^1@09NC$op0$M?}`#{FB|z}t4=b+G5WA6}|5GQfo_ zSrZyE@`0KYx}$zex9lx+v2=|2EcPCNR<+?|klDkN^;n`OquU=aex`GeXiHXghP%1s zClFns6--NVB#eNUcD;PFmnO?LWsj|$R$Zrsq1I_*<*ot7H+{OhF+(lM>1?vW4j;2K!s8b1JqgvgZnTn_}9oe z_n6Z>6#O%=X-}1>@KLjXVQdyt8%=hC@-9q5;wBW}FyB-cDyUsG=XO zL1a&KZC~8)e>payJ>0(~S{VVqbsfe7t3Hygu25#u=yVUujczzve%Dh zIFVm)GCq>2**|Sq=RIwB@i^zg>yafF@^pcdKIZbbasb47lsSa=4OEql$Cw!$fVz8TUJn7r{ch9T7?oVTV%;yu;(3-Zab;RBCG`7#X zA%+&v2 zGUMu4`O9u!28dqxU>6+iK1UIbyH;0$zLOsOBYK*U^U@hH#(=`~M(xRSgE!TdI z7-QXYl2T+7JA_=#l~F@#t>;naK-)vz_b~_ZDH)u)-pfYI^}H__UL#B3nV}=$^Th2s zUgwO%&q5EGi?7A+ADoNVYG@2a8`bTG_ER+vnZLbehw5kH$rE+% z#T{ZDvZndw9mkvSkaD7%>_c`W<9Mo#^dx8a==j~T*oyD)4E}aY@DK}PGjt-PEtb|2 zR!3j%S>v=h=N3OiN??delFY&4b)e0${_Ld2-jQopy{9aT&ip3JiT0AL5f`y$r8t}H zuWCG?DK%;2SmjF**0W+sZL07y58H59=k{?_rdu`MRAyy?Y(J={vBExAmU(`EBbb*7fUM z@brGPw$}!>!&s386-U55KJ6qua^CqueGVKl-y>hI2c2PE?Quv83qJe>d;`@fvReAx zh(@kuvi5V~p(%UlXL{dw)idUo^o_m0bb1^=WBXKVJF~Lma>uqfd%Qw-NSZ+dw7#tM`?^-rz|*aCf1Z_?*Vk05qiArQf~I<& z`5K&A^JYHr7apYZ;s^1a-C`%pU#d;WFyk5FtH@KK&Ya32+O~X1WsP{;x<(^WW1&50 z`R#bAl|t|3we4jq=|lg@`~7wGsh;rV_$>{pZ}hHj+$95ZeiS~ur)R~>%62D9U;d4k z2sZF{p2a`uB7QwTd~Q2(s>frDGx6k`fk9yL^JgcJ_Y-Hz?5T67!dLa`Ckx)OamK~P z$m>#Xcs~BAP0KUZ4K+VgDV{w*4f0CoL#vU;X63Z5i-Ws?X;`tj-YV-XLE=m2sOe2L z%=4nR#P^t2Biqj?Lkm6)oljp#=z8qTfvu+U5Vm@ZA#9RlMX|&tde5A_E6y>oB}Xsl z1b^D`OXIdwUNu+A@Mbeii0An+5{H{yz z2uGO~{~6!mRY`>^!LYq(SuomG{7Z^CqX!wnUq(}~Pwva|-F1svDI?Oh-M7W}_^tOr z?nH8j=fFx^&L1Cu3@=`XYNz~;!QtKvRgLLn3UIa^xMS|h0syV5%)v7mh^(}1eIw2m zWG>!iWK6>ojDaPhkNuD$&urS4F|Cij-{il&8T2B2x4AcyyWF50oNT&eL~4JA5v0NW zjW?8z#i{P2&Y9W!sCCTs;F-THSf!}pp!a%hSfeRQI2JOx6S9CW3opaT_UOG&sg@JV zUJc1nbrJ6j>g#>R8we~``?B6N!QZ=}%3+2hZ6F7c$=n+^m;Awc3g?bHk zzoALc5mZvMy9~*e@8^BDTJ_UgUsPeu`&yqxJ2*xERWrxFg!(!=vcWsMVSQck9oJjE z;5Hf2`C0X;=NQk|1%t0M7JeL->3)0;>j~^9!fRWLKbDybBhcRXQ?*D7+h+kko{9Y& z{4rJMOMl`mOMeQs&!RPU)9tO%>BrY~_M6_$oADr=MojIKgy2qFU)~7P!o{fC|BBpY zMiO3;mM71dd=9@u&4lc9$e!&T=id>V8fqZ*`yQbpcFk?DZTN)jMWa;@erT-z!#hnk1pZn1M{p|=+wE$u0BjHn04xwgx^FZ&{ybF1&8 zSM9(Omcbh~&34DnwicCWiG1zT7_GIa&c%B+r;6N;(Nwaf*o(1|JI|HuKbd_r5j!2> zHQj5m+}dBYpdX&6R!V!dYwm#--Yv?zmIQsRHGu-JL#KjDw()0;kvWz)F)h7<$J2B5 z^-=`}%&Zr;v?cG1_FHX{TOB36bTyv28qcZ%2H#%#9zCmlBGM5T%1RS_ezl9O4M{^TjO?;)J`l%Mrj(yQyn z>LDD%#&p(09Gg4#R2j~imTUjTE9f=70i0S{Y?-1$^4a}*s^b(^6CXuCU=Xb?dxbue zboKl)o$qyz*AJgXN2$(Sy-q)C>O70~$N>_LKBuY=_#VW!gOF_2nOIeMV+r+eUhc_$ z#?)`e)a%y-z$4h%bX5Q?C|?TwdZ*Y;>#8aHNqf$Qg%(*mK(uL>P801=f-rwXNzfGMA>i z8tt;TAusf>v@7SF?^{2w7@ORhAX>dA8c$`@5gCtJlS^R?WhqzHpm<&0{sX#q92@F; zQ;5JpW>A1hF~uC1ur$yucQcMd<}rq_Sq586${ZcTWE_-IC|yrRPVs(HOS+@3qH{lW z1Y(U)bH|)m$o8VPHOD5rs5&-k?Rx8%)^J&4NuKOO^{Y*u(OzJw#apg0c~2TUEEJ*P zvNJ_fy-DE5`?SXkATO>ZV_fe#z&Ra%1A@?`d(7i zbH@v2Un3Tin1b=Zev7lr?`UlwkCRx$UL}~(JIOhB%Mr9!ft`8U%(HB9OpadI#cQZJ zM0=lyR@!oEKV)^Se7<|W_c`1T-ZoN>1ap?tR68m{)N+&azx&=~>DrmXp~P>tq>O_e zIT3KdpSOaBikN^^J3qvkP=H!xEOd=KBICYC^6>OXwRGQyV^V7)~a{n{Jd4mu*c zuV-(yaZ|69knhRbe*BB}B=*gBSDO>iGxy72TFw}ymhcslaLekxWWvod^^bKrY@Tl-aHm+R44eKtlt@|94Wg{--vP8 zZ+$=h>lFNiSGmk%!B2IsFAIxo$Y!wAn3@ZXxwiFIdhfH{%efb`XWq!LR-deOa~I7{ za2=5>z9IOqA0}_Juj9djDvxJ>sNcPrk+WCOlzFSOxH-CgGTNrKul-L|k7CWWI-oP( zX-|fq(T2MaTO4_>L|c7}N5K^`(c>7P((7UPiTNCF-B|6_3g+K+0ZqyAi-JWWRsp@F2W--sAmIfFVL<-X!+X9}sl7-}NJ>iOn5 zY0CGbbIFR~Z4#G^@3ptwk@bcq_GbC^qRaf)PPHTDc`f+ASW@K$v^#>j5!L0^QiKgp zS;In56{1~=K!L?kRla#sK6d|pF6$mcW}J`*WB%s9vE6fnUO6e1J?=!;d%>I7F6EZ-iM8^FKgvq4-Xz1B1K^nz zq1AdBeJsxM8d>{J4rOq?I_Wlfns)X-t%cF0^*dN}s=suw&M8X;E8Z*yDH zUE|l)((4&&v)%{nB^XF}R&SJ59#p#^MKf#z_m2hk*=5Mp+$F);P-jPXDD@mTHtyA+ zrTjZaIJblYJ{*Stx{?nf{>OIGida-9pVk??DLLq&@Cp>WzJ)>l3a7pzt0Dh?HGcZ+ ze`@5HgRQ2pp{b#q3)%E}=(Pczix z26~)__G4_Vk=qi=C&t@+6|ew79V<@b1CQ~C;_J?XXpWxU3u$29{n_2jHkUX+<526p{ zTUpa`j*b3oQ_}*tCWTyfy_VKzrt7nr>$#DYZef!jjvYdWD+-c#r#Gl@MV>nIBo5(h zS@OSBJCUPgT@Y+?N7WDP76*sQR0A`+)6kP(+V4x|8P#jjli#iFTsheBiBzn%oru=N z20oi!=d!udLlon!YRrP2gv1!cnnWBlzB|;u;dQ&LKe^$=T^jj)&|4VU)+mhGe?TG_w?Im z+#Lle%SwcY|3geWCuixXZA`0(zwNd7)3bJjL|x<7;_u3s$QeEsvgc|ct*?x0CZfh# z|41W04ZV?PLx@7mChbG&P($ULWcX8=Gt5uQpp;`_6OL!^gEDTHx4HO2>3mt*3c+(9_UJ=<-@qr>+%IsM8TU=33{ zEo4n(clwNK>mf_mq3@rRoR95rTlxKc^sIT~pP&z=6VSJNF}gVQ=(`K*ciPOcC*T&I z0v<)i_y;YuxsIpwptPx}1D0MzR+1W;bven^E8n(B=O_Dlbv6kP5Uy38QMh>=*wYB( z_^`G*9=H|<&s_tA)2RdBbhXX0?~(ui(UWAJwGJ=cSnDL=Zrw_e=icg-OGPum$y-r{ zgeOOa7>`@k=i)GFCdCoxAUrNK8uBWgf`5S&5I^;X^YevcOYP^$u@vQp~})Iw8fLvLpuFfm`j%qDN_Wa&wx$K3iOYe?0fZC2=VJ ztSyNmGoU77(paJ@Z&UiaxaPGUw43T;Yuff$_#*mHZF}9yF#tTI_ne^Z!4KF6a<7i5 z>b>^RXdyoGUbchAUVzW~Y*TvE88O3mM%jjVrhBFD;kp|-Ju`3fPXA1c+~v^Q?5Q4n_U?li z)!y51eX)NUs_AWG(^&H}#-5%rrm&rQr>gpyTJv!)rj%#CdM#D&UL+1<_0JVX`<=b$#&Q%hvVv zu;*y$n$e(N?JS3~A`qn8b*7-E*AXn?pc68N3vZhC@tZK=& z1vy|Sdofx}?T|FR+FHb4N1_^J;zf90i2*Wi;}8#X8y;t~87F0HY)aFOW7?4|^Ck{@ nnR$EEokGnv&dPXG_|OS_<5=*V-CxTSt7b9Y-w>