diff --git a/.streamlit/config.toml b/.streamlit/config.toml new file mode 100644 index 0000000..988efa0 --- /dev/null +++ b/.streamlit/config.toml @@ -0,0 +1,14 @@ +[server] +headless = true +enableCORS = false +enableXsrfProtection = true + +[browser] +gatherUsageStats = false + +[theme] +primaryColor = "#045C64" +backgroundColor = "#FFFFFF" +secondaryBackgroundColor = "#F4F6F4" +textColor = "#2D3436" +font = "sans serif" diff --git a/MANIFEST b/MANIFEST index 7bc0d92..47f6b22 100644 --- a/MANIFEST +++ b/MANIFEST @@ -7,6 +7,13 @@ psisloo.py setup.py stan_models.py stan_utility.py -statistics.py +bioce_statistics.py variationalBayesian.py +fullBayesianMultimodal.py +multimodal_io.py +prepareHDX.py +prepareXLMS.py +bioce_pipeline.py +structure_observables.py +app.py vbw_sc.i diff --git a/README.md b/README.md index acbba55..8a642a7 100644 --- a/README.md +++ b/README.md @@ -32,6 +32,27 @@ python variationalBayesian.py –-help ``` If you see no errors but options menu pops up, you are good to go. +### PyStan 3 (full Bayesian inference) + +Complete Bayesian sampling (`fullBayesian.py`, `fullBayesianTR.py`) uses **PyStan 3**: +install the **`pystan`** package from pip (inside the activated conda env the `bioce.yml` `pip` section already requests `pystan>=3.9,<4`). In code this appears as **`import stan`** (not `import pystan`). See the [PyStan upgrading guide](https://pystan.readthedocs.io/en/latest/upgrading.html). + +**Platforms:** PyStan 3 is supported on **Linux and macOS** only (not Windows). + +**Smoke test** (checks that `stan` compiles and samples, and that `stan_run` / `stan_utility` run): + +```bash +pip install pytest +cd /path/to/bioce +python -m pytest tests/test_pystan_smoke.py -q +``` + +The first test only checks **warmup vs sampling split** logic and passes without PyStan. The other two tests **compile and sample** a tiny model; if `pystan` is missing they are **skipped** (pytest exit code 0). With PyStan installed, all three tests should pass. + +Some httpstan builds reject sampler options such as ``parallel_chains`` or ``show_progress``; ``stan_run.build_and_sample`` then retries with only ``num_chains``, ``num_warmup``, and ``num_samples`` (parallelism and progress follow backend defaults). + +`pytest.ini` in the repo root sets `pythonpath = .` so `stan_run` imports work without installing the package. + ## Problems with installation on OSX 10.14 (Mojave) There is a known issue with xcode installation on OSX 10.14 (Mojave) If you see the following error: @@ -98,6 +119,56 @@ which graphically ilustrated distribution of population weights (shown below) an 3. Script also returns text file containing Q vector, experimental intensity, model intensity and experimental error. +### Streamlit web app (SciLifeLab Serve) + +A browser UI runs the full PDB → observables → Stan pipeline: + +```bash +pip install streamlit biopython "pystan>=3.9,<4" +streamlit run app.py +``` + +Upload PDBs (or zip), experimental SAXS / HDX / XL files, and run inference. See `serve/SERVE.md` for Docker deployment on [SciLifeLab Serve](https://serve.scilifelab.se/docs/). + +SAXS from PDB requires **FoXS** or **Pepsi-SAXS** on `PATH`. HDX and XL use **BioPython** (peptide SASA proxy and Cα distances). + +### HDX-MS and XL-MS (multimodal inference) + +Ensemble-averaged HDX uptake and probabilistic XL-MS distance restraints can be combined with SAXS (or used alone) via `fullBayesianMultimodal.py`. Per-conformer observables must be precomputed (e.g. protection factors or Cα distances); Stan multiplies independent likelihoods over a shared `simplex` weight vector. + +**Prepare inputs** + +```bash +# One file per conformer in a directory (same row order: flattened peptide×time) +python prepareHDX.py -d hdx_predictions/ -o SimulatedHDX.txt + +# One distance file per conformer (rows = crosslinks, same order as xl_restraints.dat) +python prepareXLMS.py -d xl_distances/ -o SimulatedXLdistances.txt +``` + +`hdx_exp.dat`: two columns `uptake sigma` per observation. +`xl_restraints.dat`: `res_i res_j z [d_max] [tau]` (`z=1` observed link; default `d_max=30`, `tau=3` Å). + +**Run inference** + +```bash +# SAXS + HDX + XL +python fullBayesianMultimodal.py \ + -p weights.txt -f structures.txt \ + -e simulated.dat -s SimulatedIntensities.txt \ + -H hdx_exp.dat -S SimulatedHDX.txt \ + -X xl_restraints.dat -D SimulatedXLdistances.txt \ + -i 2000 -c 4 -j 4 + +# HDX only, or XL only (same -p priors) +python fullBayesianMultimodal.py -p weights.txt -H hdx_exp.dat -S SimulatedHDX.txt +python fullBayesianMultimodal.py -p weights.txt -X xl_restraints.dat -D SimulatedXLdistances.txt +python fullBayesianMultimodal.py -p weights.txt -X xl_restraints.dat -D SimulatedXLdistances.txt --xl-fp-mixture +``` + +XL satisfaction uses a soft logistic in linker distance: ensemble probability +`p_sat = sum_k w_k * sigmoid((d_max - d_k) / tau)` with `z_l ~ Bernoulli(p_sat)`. + ### Using chemical shift data 1. In order to use chemical shift data, one needs to install SHIFTX2. This can be done by following instructions at: [SHIFTX2](http://www.shiftx2.ca/download.html). This requires running python 2.6 or later, which won't work with diff --git a/app.py b/app.py new file mode 100644 index 0000000..e40e693 --- /dev/null +++ b/app.py @@ -0,0 +1,242 @@ +""" +Bioce Streamlit app: upload PDB ensemble + SAXS / HDX / XL data → Bayesian weights. + +Run locally: + streamlit run app.py +""" +from __future__ import annotations + +import tempfile +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import streamlit as st + +from bioce_pipeline import PipelineConfig, run_full_pipeline +from scilifelab_theme import ( + AQUA, + GRAPE, + LIME, + TEAL, + apply_matplotlib_theme, + inject_streamlit_theme, + render_footer, + render_header, +) +from structure_observables import _foxs_on_path, _HAS_BIOPYTHON, _pepsi_on_path + +st.set_page_config( + page_title="Bioce | SciLifeLab", + layout="wide", + initial_sidebar_state="expanded", +) +inject_streamlit_theme() +apply_matplotlib_theme() +render_header( + "Bioce — Bayesian ensemble inference", + "Upload PDB conformers and SAXS, HDX-MS, and/or XL-MS data to infer posterior ensemble weights.", +) + +with st.sidebar: + st.markdown( + '

' + "SciLifeLab · Bioce

", + unsafe_allow_html=True, + ) + st.header("Sampling") + iterations = st.slider("Stan iterations (total per chain)", 100, 3000, 400, 100) + chains = st.slider("Chains", 1, 8, 2) + njobs = st.slider("Parallel chains", 1, 8, 2) + st.header("Modalities") + use_saxs = st.checkbox("SAXS", value=True) + use_hdx = st.checkbox("HDX-MS", value=False) + use_xl = st.checkbox("XL-MS", value=False) + xl_fp = st.checkbox("XL false-positive mixture", value=False, disabled=not use_xl) + st.divider() + st.markdown("**Dependencies**") + st.write("FoXS:", "yes" if _foxs_on_path() else "no") + st.write("Pepsi-SAXS:", "yes" if _pepsi_on_path() else "no") + st.write("BioPython:", "yes" if _HAS_BIOPYTHON else "no") + +st.subheader("Structure library") +pdb_uploads = st.file_uploader( + "PDB files or .zip archive", + type=["pdb", "zip"], + accept_multiple_files=True, + help="One file per conformer in the ensemble.", +) + +col1, col2 = st.columns(2) + +with col1: + st.subheader("SAXS") + saxs_file = st.file_uploader( + "Experimental SAXS (.dat)", + type=["dat", "txt", "csv"], + disabled=not use_saxs, + help="Columns: q, I, sigma (whitespace-separated).", + ) + +with col2: + st.subheader("HDX-MS") + hdx_peptides = st.file_uploader( + "Peptide definitions", + type=["txt", "dat", "csv"], + disabled=not use_hdx, + help="One peptide per line: chain start_res end_res (e.g. A 10 25).", + ) + hdx_exp = st.file_uploader( + "Experimental uptake", + type=["txt", "dat", "csv"], + disabled=not use_hdx, + help="Lines: uptake sigma — or pep_idx time uptake sigma for kinetics.", + ) + +st.subheader("XL-MS") +xl_file = st.file_uploader( + "Crosslink restraints", + type=["txt", "dat"], + disabled=not use_xl, + help="res_i res_j z [d_max] [tau] — or chain_i res_i chain_j res_j z [d_max] [tau].", +) + +with st.expander("File format help"): + st.markdown( + """ +- **SAXS**: three columns `q`, `I`, `error` (same as bioce `fullBayesian.py`). +- **HDX peptides**: `chain start end` per line (1-based residue numbers). +- **HDX data**: `uptake sigma` per peptide, or `pep_idx time uptake sigma` for time series. +- **XL**: `z=1` if the crosslink was observed. Distances are computed as Cα–Cα from each PDB. +- **HDX prediction** uses mean peptide SASA (Shrake–Rupley) as a solvent-exposure proxy—not full HDX kinetics. + """ + ) + +run = st.button("Run inference", type="primary", use_container_width=True) + +if run: + if not pdb_uploads: + st.error("Upload at least one PDB file or a .zip archive.") + st.stop() + if use_saxs and not saxs_file: + st.error("SAXS enabled: upload an experimental curve.") + st.stop() + if use_hdx and (not hdx_peptides or not hdx_exp): + st.error("HDX enabled: upload peptide definitions and experimental uptake.") + st.stop() + if use_xl and not xl_file: + st.error("XL enabled: upload crosslink restraints.") + st.stop() + if use_saxs and not (_foxs_on_path() or _pepsi_on_path()): + st.error("Install FoXS or Pepsi-SAXS on PATH for SAXS simulation from PDBs.") + st.stop() + if (use_hdx or use_xl) and not _HAS_BIOPYTHON: + st.error("Install BioPython for HDX/XL from PDB (pip install biopython).") + st.stop() + if not (use_saxs or use_hdx or use_xl): + st.error("Enable at least one modality.") + st.stop() + + with tempfile.TemporaryDirectory(prefix="bioce_") as tmp: + work = Path(tmp) + config = PipelineConfig( + work_dir=work, + use_saxs=use_saxs, + use_hdx=use_hdx, + use_xl=use_xl, + iterations=iterations, + chains=chains, + njobs=njobs, + xl_fp_mixture=xl_fp, + ) + if use_saxs: + saxs_path = work / "experimental_saxs.dat" + saxs_path.write_bytes(saxs_file.getvalue()) + config.saxs_experimental = saxs_path + if use_hdx: + config.hdx_peptide_lines = hdx_peptides.getvalue().decode().splitlines() + config.hdx_exp_lines = hdx_exp.getvalue().decode().splitlines() + if use_xl: + config.xl_restraint_lines = xl_file.getvalue().decode().splitlines() + + with st.status("Running pipeline…", expanded=True) as status: + try: + st.write("Extracting PDBs and computing observables…") + result = run_full_pipeline(pdb_uploads, config) + status.update(label="Done", state="complete") + except Exception as exc: + status.update(label="Failed", state="error") + st.exception(exc) + st.stop() + + m1, m2, m3, m4 = st.columns(4) + if result.saxs_chi2 is not None: + m1.metric("SAXS χ²", f"{result.saxs_chi2:.3f}") + if result.jsd is not None: + m2.metric("JSD", f"{result.jsd:.4f}") + if result.hdx_rmse is not None: + m3.metric("HDX RMSE", f"{result.hdx_rmse:.4f}") + if result.xl_mean_psat is not None: + m4.metric("Mean XL p(sat)", f"{result.xl_mean_psat:.3f}") + + st.subheader("Posterior mean weights") + df = pd.DataFrame(result.tables["weights"]) + st.dataframe(df, use_container_width=True) + fig, ax = plt.subplots(figsize=(8, 3)) + colors = [TEAL if w == df["weight"].max() else LIME for w in df["weight"]] + ax.bar(df["structure"], df["weight"], color=colors, edgecolor=TEAL, linewidth=0.6) + ax.set_ylabel("Weight") + ax.set_facecolor("#FAFCFA") + ax.grid(axis="y", linestyle="--", alpha=0.5) + ax.tick_params(axis="x", rotation=45) + plt.tight_layout() + st.pyplot(fig) + plt.close(fig) + + if "saxs_fit" in result.plots: + st.subheader("SAXS fit") + curve = np.genfromtxt(result.plots["saxs_fit"]) + fig2, ax2 = plt.subplots(figsize=(7, 4)) + ax2.errorbar( + curve[:, 0], curve[:, 1], yerr=curve[:, 3], + fmt="o", ms=3, color=GRAPE, ecolor=AQUA, label="Experiment", alpha=0.85, + ) + ax2.plot( + curve[:, 0], curve[:, 2], "-", color=TEAL, + label="Weighted model", lw=2, + ) + ax2.set_xlabel("q (Å⁻¹)") + ax2.set_ylabel("I(q)") + ax2.legend(frameon=True) + ax2.set_facecolor("#FAFCFA") + ax2.set_yscale("log") + st.pyplot(fig2) + plt.close(fig2) + + for key, path in result.plots.items(): + if path.suffix == ".png": + st.image(str(path), caption=key) + + if "xl_psat" in result.tables: + st.subheader("XL satisfaction probability (posterior mean)") + psat = result.tables["xl_psat"] + fig_xl, ax_xl = plt.subplots(figsize=(7, 2.8)) + ax_xl.plot(psat, color=TEAL, marker="o", markersize=5, markerfacecolor=LIME) + ax_xl.fill_between(range(len(psat)), psat, alpha=0.2, color=AQUA) + ax_xl.set_ylim(0, 1) + ax_xl.set_xlabel("Crosslink index") + ax_xl.set_ylabel("p(satisfied)") + ax_xl.set_facecolor("#FAFCFA") + ax_xl.grid(axis="y", linestyle="--", alpha=0.5) + st.pyplot(fig_xl) + plt.close(fig_xl) + + st.download_button( + "Download weights CSV", + df.to_csv(index=False).encode(), + file_name="bayesian_weights.csv", + mime="text/csv", + ) + +render_footer() diff --git a/bioce.yml b/bioce.yml index 7aa1a41..2809bce 100644 --- a/bioce.yml +++ b/bioce.yml @@ -6,11 +6,14 @@ channels: - intel dependencies: - seaborn +- biopython +- pandas +- streamlit - clang - openmp - swig - gsl - pip: - - pystan + - "pystan>=3.9,<4" - pandas - git+https://github.com/emblsaxs/sasciftools.git diff --git a/bioce_pipeline.py b/bioce_pipeline.py new file mode 100644 index 0000000..908e7a6 --- /dev/null +++ b/bioce_pipeline.py @@ -0,0 +1,262 @@ +""" +End-to-end pipeline: PDBs + experimental data → observables → Stan inference. +""" +from __future__ import annotations + +import os +from dataclasses import dataclass, field +from pathlib import Path +from typing import Any, Dict, List, Optional, Sequence + +import numpy as np + +from bioce_statistics import calculateChiCrysol, JensenShannonDiv, mean_for_weights +from fullBayesian import calculate_stats, combine_curve, execute_stan +from fullBayesianMultimodal import ( + execute_stan_hdx, + execute_stan_hdx_xl, + execute_stan_saxs_hdx, + execute_stan_saxs_hdx_xl, + execute_stan_saxs_xl, + execute_stan_xl, +) +from structure_observables import ( + extract_pdbs, + generate_hdx_matrix, + generate_saxs_matrix, + generate_xl_distance_matrix, + parse_hdx_experimental, + parse_hdx_peptides, + parse_xl_restraints, + write_hdx_experimental, +) + + +@dataclass +class PipelineConfig: + work_dir: Path + use_saxs: bool = True + use_hdx: bool = False + use_xl: bool = False + iterations: int = 400 + chains: int = 2 + njobs: int = 2 + saxs_experimental: Optional[Path] = None + hdx_peptide_lines: Optional[List[str]] = None + hdx_exp_lines: Optional[List[str]] = None + xl_restraint_lines: Optional[List[str]] = None + xl_fp_mixture: bool = False + + +@dataclass +class PipelineResult: + fit: Any + work_dir: Path + structure_names: List[str] + bayesian_weights: np.ndarray + saxs_chi2: Optional[float] = None + jsd: Optional[float] = None + hdx_rmse: Optional[float] = None + xl_mean_psat: Optional[float] = None + plots: Dict[str, Path] = field(default_factory=dict) + tables: Dict[str, Any] = field(default_factory=dict) + + +def _read_experimental_saxs(path: Path) -> np.ndarray: + return np.genfromtxt(path, dtype=np.float64) + + +def _uniform_priors(n: int) -> np.ndarray: + return np.full(n, 1.0 / n) + + +def prepare_from_pdbs( + pdb_uploads, + config: PipelineConfig, +) -> Dict[str, Path]: + """Extract PDBs and write simulated observable matrices into work_dir.""" + config.work_dir.mkdir(parents=True, exist_ok=True) + pdb_dir = config.work_dir / "pdbs" + paths = extract_pdbs(pdb_uploads, pdb_dir) + n = len(paths) + names = [p.stem for p in paths] + artifacts: Dict[str, Path] = {"pdb_dir": pdb_dir} + + priors = _uniform_priors(n) + priors_path = config.work_dir / "weights.txt" + np.savetxt(priors_path, priors, fmt="%.8g") + names_path = config.work_dir / "structures.txt" + with open(names_path, "w") as handle: + handle.write(" ".join(names) + "\n") + artifacts["priors"] = priors_path + artifacts["names"] = names_path + + if config.use_saxs: + if not config.saxs_experimental: + raise ValueError("SAXS experimental file required.") + sim_p, w_p, n_p, _ = generate_saxs_matrix( + paths, config.saxs_experimental, config.work_dir / "saxs" + ) + artifacts["sim_saxs"] = sim_p + artifacts["priors"] = w_p + artifacts["names"] = n_p + shutil_copy_exp = config.work_dir / "saxs" / config.saxs_experimental.name + artifacts["exp_saxs"] = shutil_copy_exp + + if config.use_xl: + if not config.xl_restraint_lines: + raise ValueError("XL restraints required.") + records = parse_xl_restraints(config.xl_restraint_lines) + xl_p, sim_xl, _ = generate_xl_distance_matrix( + paths, records, config.work_dir / "xl" + ) + artifacts["xl_restraints"] = xl_p + artifacts["sim_xl"] = sim_xl + + if config.use_hdx: + if not config.hdx_peptide_lines or not config.hdx_exp_lines: + raise ValueError("HDX peptides and experimental data required.") + peptides = parse_hdx_peptides(config.hdx_peptide_lines) + uptake, sigma, pep_idx = parse_hdx_experimental( + config.hdx_exp_lines, len(peptides) + ) + hdx_exp_path = config.work_dir / "hdx" / "hdx_exp.dat" + hdx_exp_path.parent.mkdir(parents=True, exist_ok=True) + write_hdx_experimental(uptake, sigma, hdx_exp_path) + sim_h, _ = generate_hdx_matrix( + paths, peptides, config.work_dir / "hdx", pep_indices=pep_idx + ) + artifacts["hdx_exp"] = hdx_exp_path + artifacts["sim_hdx"] = sim_h + + return artifacts + + +def run_inference( + config: PipelineConfig, + artifacts: Dict[str, Path], +) -> PipelineResult: + """Run Stan after prepare_from_pdbs.""" + work = config.work_dir + old_cwd = os.getcwd() + os.chdir(work) + try: + priors = np.genfromtxt(artifacts["priors"]) + names = open(artifacts["names"]).read().split() + + has_saxs = config.use_saxs and "sim_saxs" in artifacts + has_hdx = config.use_hdx and "sim_hdx" in artifacts + has_xl = config.use_xl and "sim_xl" in artifacts + + it, ch, nj = config.iterations, config.chains, config.njobs + + if has_saxs and has_hdx and has_xl: + exp = _read_experimental_saxs(artifacts["exp_saxs"]) + sim = np.genfromtxt(artifacts["sim_saxs"]) + fit = execute_stan_saxs_hdx_xl( + exp, sim, + str(artifacts["hdx_exp"]), str(artifacts["sim_hdx"]), + str(artifacts["xl_restraints"]), str(artifacts["sim_xl"]), + priors, it, ch, nj, + ) + prefix = "saxs_hdx_xl" + elif has_saxs and has_hdx: + exp = _read_experimental_saxs(artifacts["exp_saxs"]) + sim = np.genfromtxt(artifacts["sim_saxs"]) + fit = execute_stan_saxs_hdx( + exp, sim, + str(artifacts["hdx_exp"]), str(artifacts["sim_hdx"]), + priors, it, ch, nj, + ) + prefix = "saxs_hdx" + elif has_saxs and has_xl: + exp = _read_experimental_saxs(artifacts["exp_saxs"]) + sim = np.genfromtxt(artifacts["sim_saxs"]) + fit = execute_stan_saxs_xl( + exp, sim, + str(artifacts["xl_restraints"]), str(artifacts["sim_xl"]), + priors, it, ch, nj, + ) + prefix = "saxs_xl" + elif has_hdx and has_xl: + fit = execute_stan_hdx_xl( + str(artifacts["hdx_exp"]), str(artifacts["sim_hdx"]), + str(artifacts["xl_restraints"]), str(artifacts["sim_xl"]), + priors, it, ch, nj, + ) + prefix = "hdx_xl" + elif has_saxs: + exp = _read_experimental_saxs(artifacts["exp_saxs"]) + sim = np.genfromtxt(artifacts["sim_saxs"]) + fit = execute_stan(exp, sim, priors, it, ch, nj) + prefix = "saxs" + elif has_hdx: + fit = execute_stan_hdx( + str(artifacts["hdx_exp"]), str(artifacts["sim_hdx"]), + priors, it, ch, nj, + ) + prefix = "hdx" + elif has_xl: + fit = execute_stan_xl( + str(artifacts["xl_restraints"]), str(artifacts["sim_xl"]), + priors, it, ch, nj, use_fp=config.xl_fp_mixture, + ) + prefix = "xl" + else: + raise ValueError("No modality enabled.") + + weights = np.asarray(mean_for_weights(fit)) + result = PipelineResult( + fit=fit, + work_dir=work, + structure_names=names, + bayesian_weights=weights, + ) + + plots = {} + for name in ( + "{}_weights.png".format(prefix), + "stan_weights.png", + "{}_scale.png".format(prefix), + "stan_scale.png", + ): + p = work / name + if p.exists(): + plots["weights" if "weight" in name else "scale"] = p + + if has_saxs: + exp = _read_experimental_saxs(artifacts["exp_saxs"]) + sim = np.genfromtxt(artifacts["sim_saxs"]) + bw, jsd, chi2 = calculate_stats(fit, exp, sim) + result.bayesian_weights = bw + result.jsd = float(jsd) + result.saxs_chi2 = float(chi2) + curve_path = work / "bayesianEstimateCurve.txt" + if curve_path.exists(): + plots["saxs_fit"] = curve_path + + if has_hdx: + from multimodal_io import read_hdx_experimental + tu, _ = read_hdx_experimental(str(artifacts["hdx_exp"])) + sim_h = np.genfromtxt(artifacts["sim_hdx"]) + pred = sim_h.T @ result.bayesian_weights + result.hdx_rmse = float(np.sqrt(np.mean((tu - pred) ** 2))) + + if has_xl and "p_sat" in fit: + psat = np.asarray(fit["p_sat"]) + result.xl_mean_psat = float(np.mean(psat)) + result.tables["xl_psat"] = np.mean(psat, axis=1) + + result.plots = plots + result.tables["weights"] = { + "structure": names, + "weight": result.bayesian_weights, + } + return result + finally: + os.chdir(old_cwd) + + +def run_full_pipeline(pdb_uploads, config: PipelineConfig) -> PipelineResult: + artifacts = prepare_from_pdbs(pdb_uploads, config) + return run_inference(config, artifacts) diff --git a/statistics.py b/bioce_statistics.py similarity index 87% rename from statistics.py rename to bioce_statistics.py index c6ccfd9..05e2487 100644 --- a/statistics.py +++ b/bioce_statistics.py @@ -64,17 +64,14 @@ def JensenShannonDiv(weights_a, weights_b): def mean_for_weights(fit): """ - Calculation mean weight from stan fitting object + Calculation mean weight from stan fitting object (PyStan 3 ``stan`` fit). :param fit: stan fit model :return: """ - stan_chain = fit.extract() - weights = stan_chain["weights"] - vlen = np.shape(weights) - mean_weights = [] - for ind in range(vlen[1]): - mean_weights.append(np.mean(weights[:,ind])) - return mean_weights + weights = np.asarray(fit["weights"]) + if weights.ndim == 1: + weights = weights.reshape(1, -1) + return list(np.mean(weights, axis=1)) def me_log_lik(fit): """ @@ -82,8 +79,7 @@ def me_log_lik(fit): :param log_lik: :return: """ - stan_chain = fit.extract() - log_lik = stan_chain["loglikes"] + log_lik = np.asarray(fit["loglikes"]).T return np.sum(np.mean(log_lik, axis=1)) diff --git a/data/bayesianEstimateCurve.txt b/data/bayesianEstimateCurve.txt new file mode 100644 index 0000000..f453f9f --- /dev/null +++ b/data/bayesianEstimateCurve.txt @@ -0,0 +1,868 @@ +5.658709999999999955e-03 2.605899000000000063e-02 2.616062206258566031e-02 6.538312999999999590e-03 +6.048959999999999927e-03 2.604476999999999834e-02 2.614642553244841960e-02 2.102154000000000175e-03 +6.439219999999999840e-03 2.602932999999999983e-02 2.613093788388811539e-02 1.494653000000000084e-03 +6.829469999999999812e-03 2.601389000000000132e-02 2.611545728295707608e-02 9.480715999999999710e-04 +7.219729999999999724e-03 2.599862000000000006e-02 2.610009012072817686e-02 8.027344999999999967e-04 +7.609979999999999696e-03 2.598347000000000157e-02 2.608481282367591803e-02 6.931519000000000135e-04 +8.000240000000000476e-03 2.596833000000000127e-02 2.606953600827629711e-02 6.532776999999999811e-04 +8.390490000000000448e-03 2.595048000000000146e-02 2.605169141831172389e-02 6.185780999999999788e-04 +8.780750000000000360e-03 2.593263000000000165e-02 2.603384489006434904e-02 5.481044000000000237e-04 +9.171000000000000332e-03 2.591364999999999988e-02 2.601488655363431990e-02 4.967650000000000321e-04 +9.561260000000000245e-03 2.589324000000000139e-02 2.599452072152295609e-02 4.533153000000000269e-04 +9.951510000000000217e-03 2.587282999999999944e-02 2.597415471944051341e-02 4.240569999999999821e-04 +1.034179999999999995e-02 2.585194000000000034e-02 2.595310307588851589e-02 4.169377000000000044e-04 +1.073200000000000022e-02 2.583097999999999991e-02 2.593195733719368470e-02 3.988217999999999971e-04 +1.112229999999999989e-02 2.580864000000000075e-02 2.590946588873705922e-02 3.705421999999999774e-04 +1.151250000000000016e-02 2.578327000000000049e-02 2.588401921184589616e-02 3.640655000000000128e-04 +1.190279999999999984e-02 2.575790999999999845e-02 2.585857098873837329e-02 3.507869999999999738e-04 +1.229300000000000011e-02 2.573218999999999854e-02 2.583275907384328246e-02 3.442941999999999943e-04 +1.268329999999999978e-02 2.570634999999999934e-02 2.580682189420015898e-02 3.282116999999999869e-04 +1.307350000000000005e-02 2.568010000000000084e-02 2.578047611010549772e-02 3.200821999999999760e-04 +1.346379999999999973e-02 2.565202000000000107e-02 2.575234086343195297e-02 3.101631999999999783e-04 +1.385409999999999940e-02 2.562394999999999951e-02 2.572420620700981500e-02 3.058952999999999781e-04 +1.424429999999999967e-02 2.559552000000000008e-02 2.569570709559651153e-02 2.950391999999999923e-04 +1.463459999999999935e-02 2.556686999999999849e-02 2.566696773545071308e-02 3.078153999999999938e-04 +1.502479999999999961e-02 2.553804000000000143e-02 2.563803874317140477e-02 3.012835000000000166e-04 +1.541509999999999929e-02 2.550625000000000114e-02 2.560611445868011052e-02 2.733028000000000085e-04 +1.580530000000000129e-02 2.547446999999999906e-02 2.557420230719907248e-02 2.734953000000000037e-04 +1.619560000000000097e-02 2.544131999999999991e-02 2.554096425242296034e-02 2.639868000000000184e-04 +1.658580000000000124e-02 2.540685000000000165e-02 2.550644328570561212e-02 2.644279000000000030e-04 +1.697610000000000091e-02 2.537236000000000005e-02 2.547190510867966898e-02 2.621254999999999961e-04 +1.736630000000000118e-02 2.533708999999999961e-02 2.543644993791249753e-02 2.570121000000000189e-04 +1.775660000000000086e-02 2.530175999999999953e-02 2.540092452031132708e-02 2.476986999999999939e-04 +1.814680000000000112e-02 2.526528999999999928e-02 2.536432687947702200e-02 2.443009000000000109e-04 +1.853710000000000080e-02 2.522695000000000146e-02 2.532594727996597198e-02 2.489613000000000112e-04 +1.892730000000000107e-02 2.518861000000000017e-02 2.528757043286513836e-02 2.444742000000000182e-04 +1.931760000000000074e-02 2.514983999999999970e-02 2.524865907802703938e-02 2.372960999999999929e-04 +1.970780000000000101e-02 2.511098000000000149e-02 2.520964089121384422e-02 2.312833000000000087e-04 +2.009810000000000069e-02 2.507155000000000147e-02 2.517006786087300566e-02 2.374785999999999933e-04 +2.048830000000000096e-02 2.503052000000000055e-02 2.512892251980788916e-02 2.293792000000000010e-04 +2.087860000000000063e-02 2.498946999999999974e-02 2.508775786016333512e-02 2.149956999999999921e-04 +2.126880000000000090e-02 2.494715000000000127e-02 2.504529666142804936e-02 2.118533000000000135e-04 +2.165910000000000057e-02 2.490424999999999930e-02 2.500223509887694942e-02 2.192177999999999868e-04 +2.204930000000000084e-02 2.486106999999999900e-02 2.495890729806614361e-02 2.140549999999999961e-04 +2.243960000000000052e-02 2.481596999999999970e-02 2.491366631818287686e-02 2.131601000000000058e-04 +2.282980000000000079e-02 2.477087000000000039e-02 2.486843542063147533e-02 2.061498000000000096e-04 +2.322010000000000046e-02 2.472517000000000117e-02 2.482258766903275987e-02 2.121330999999999900e-04 +2.361030000000000073e-02 2.467901000000000122e-02 2.477627286130649337e-02 2.020128000000000049e-04 +2.400060000000000041e-02 2.463283999999999960e-02 2.472994973412406122e-02 2.001299999999999965e-04 +2.439080000000000067e-02 2.458416000000000073e-02 2.468109980321592234e-02 1.995239999999999958e-04 +2.478110000000000035e-02 2.453547999999999840e-02 2.463223546107055212e-02 1.976458000000000110e-04 +2.517130000000000062e-02 2.448625999999999928e-02 2.458284903446922298e-02 2.016204000000000119e-04 +2.556160000000000029e-02 2.443632999999999916e-02 2.453275775761378646e-02 2.024715000000000011e-04 +2.595180000000000056e-02 2.438641000000000073e-02 2.448268596933279515e-02 1.909114999999999966e-04 +2.634210000000000024e-02 2.433557000000000081e-02 2.443166877408394985e-02 2.043170000000000023e-04 +2.673230000000000051e-02 2.428461000000000161e-02 2.438053628808508830e-02 1.936863999999999893e-04 +2.712260000000000018e-02 2.423283999999999855e-02 2.432858389157794773e-02 1.903317000000000127e-04 +2.751280000000000045e-02 2.417932000000000067e-02 2.427488591406776197e-02 2.043404999999999928e-04 +2.790310000000000012e-02 2.412577999999999945e-02 2.422117591591303049e-02 1.948816000000000079e-04 +2.829330000000000039e-02 2.407114999999999949e-02 2.416634125602292946e-02 1.970847999999999871e-04 +2.868360000000000007e-02 2.401613999999999832e-02 2.411112750156875831e-02 1.975338000000000094e-04 +2.907380000000000034e-02 2.396099999999999966e-02 2.405578088241197565e-02 2.035050000000000108e-04 +2.946400000000000061e-02 2.390524000000000121e-02 2.399986446249081223e-02 2.001492000000000115e-04 +2.985430000000000028e-02 2.384947000000000109e-02 2.394393501132088473e-02 1.964019000000000008e-04 +3.024450000000000055e-02 2.379231000000000054e-02 2.388658683189068960e-02 1.985501000000000084e-04 +3.063480000000000023e-02 2.373429000000000164e-02 2.382837053411693967e-02 1.955279999999999887e-04 +3.102500000000000049e-02 2.367622000000000129e-02 2.377009792278205608e-02 1.928091999999999903e-04 +3.141530000000000017e-02 2.361707999999999863e-02 2.371078538673483324e-02 1.953644999999999893e-04 +3.180550000000000044e-02 2.355795000000000111e-02 2.365148633163891012e-02 1.973626999999999933e-04 +3.219579999999999664e-02 2.349814999999999890e-02 2.359147076443030705e-02 1.887171000000000096e-04 +3.258599999999999691e-02 2.343773000000000037e-02 2.353078441378800637e-02 1.907273999999999927e-04 +3.297619999999999718e-02 2.337728999999999849e-02 2.347009045720021767e-02 1.884219000000000128e-04 +3.336650000000000033e-02 2.331555999999999906e-02 2.340816778604896606e-02 1.905653999999999898e-04 +3.375670000000000059e-02 2.325375999999999832e-02 2.334618798816043284e-02 1.929118000000000120e-04 +3.414699999999999680e-02 2.319146999999999875e-02 2.328369695884240087e-02 1.909141999999999957e-04 +3.453719999999999707e-02 2.312840999999999855e-02 2.322042092153627293e-02 1.899299000000000126e-04 +3.492750000000000021e-02 2.306534000000000015e-02 2.315712794439268515e-02 1.871816000000000105e-04 +3.531770000000000048e-02 2.300048999999999844e-02 2.309206454526258404e-02 1.875427999999999893e-04 +3.570790000000000075e-02 2.293525000000000078e-02 2.302660681167787229e-02 1.775640000000000096e-04 +3.609819999999999696e-02 2.286989000000000036e-02 2.296102696796075565e-02 1.800720000000000066e-04 +3.648839999999999723e-02 2.280425000000000160e-02 2.289514687198180784e-02 1.738188000000000103e-04 +3.687870000000000037e-02 2.273860000000000117e-02 2.282925837235896122e-02 1.691143999999999948e-04 +3.726890000000000064e-02 2.267243999999999857e-02 2.276287633017694267e-02 1.711170000000000085e-04 +3.765910000000000091e-02 2.260602999999999918e-02 2.269626427694436052e-02 1.703656000000000007e-04 +3.804939999999999711e-02 2.253940000000000110e-02 2.262940981175838931e-02 1.651758999999999985e-04 +3.843959999999999738e-02 2.247125999999999985e-02 2.256096766360838474e-02 1.651951000000000135e-04 +3.882979999999999765e-02 2.240311999999999859e-02 2.249253149108529359e-02 1.570115000000000128e-04 +3.922010000000000080e-02 2.233434000000000114e-02 2.242349985648887239e-02 1.585086000000000090e-04 +3.961030000000000106e-02 2.226507999999999959e-02 2.235403099202725591e-02 1.607230000000000128e-04 +4.000059999999999727e-02 2.219580000000000164e-02 2.228453954574925733e-02 1.557550999999999955e-04 +4.039079999999999754e-02 2.212518999999999944e-02 2.221370317630664420e-02 1.543247000000000032e-04 +4.078099999999999781e-02 2.205458999999999892e-02 2.214287531483003640e-02 1.524874000000000064e-04 +4.117130000000000095e-02 2.198353999999999933e-02 2.207156502516190513e-02 1.512198999999999987e-04 +4.156150000000000122e-02 2.191194000000000128e-02 2.199968055860630420e-02 1.568401000000000029e-04 +4.195170000000000149e-02 2.184033999999999975e-02 2.192779741440599103e-02 1.512809000000000105e-04 +4.234199999999999769e-02 2.176865999999999871e-02 2.185582541119675953e-02 1.485494999999999980e-04 +4.273219999999999796e-02 2.169698000000000113e-02 2.178386728115323800e-02 1.494238999999999908e-04 +4.312239999999999823e-02 2.162477999999999970e-02 2.171138286684865318e-02 1.502061000000000119e-04 +4.351270000000000138e-02 2.155142000000000169e-02 2.163775730809103237e-02 1.504898999999999972e-04 +4.390290000000000165e-02 2.147806000000000021e-02 2.156414080090841498e-02 1.444441000000000084e-04 +4.429310000000000191e-02 2.140384999999999857e-02 2.148966756723826727e-02 1.433128999999999946e-04 +4.468339999999999812e-02 2.132933999999999872e-02 2.141489893531567726e-02 1.426393999999999882e-04 +4.507359999999999839e-02 2.125476999999999922e-02 2.134006930622955184e-02 1.384957000000000028e-04 +4.546379999999999866e-02 2.117990000000000150e-02 2.126494343211247223e-02 1.401301999999999886e-04 +4.585410000000000180e-02 2.110501999999999864e-02 2.118979691702054316e-02 1.368936000000000068e-04 +4.624430000000000207e-02 2.102901000000000076e-02 2.111348932754061880e-02 1.376049000000000011e-04 +4.663450000000000234e-02 2.095233000000000165e-02 2.103648214514772935e-02 1.410387999999999887e-04 +4.702479999999999855e-02 2.087566000000000074e-02 2.095949130355157300e-02 1.375681000000000018e-04 +4.741499999999999881e-02 2.079947999999999866e-02 2.088303758512468494e-02 1.362973000000000071e-04 +4.780519999999999908e-02 2.072330000000000005e-02 2.080657668109291494e-02 1.386550000000000065e-04 +4.819539999999999935e-02 2.064625999999999961e-02 2.072925708084319071e-02 1.386187999999999949e-04 +4.858570000000000250e-02 2.056830999999999937e-02 2.065105059226504322e-02 1.403200000000000118e-04 +4.897590000000000277e-02 2.049039000000000069e-02 2.057286657319130560e-02 1.400907999999999971e-04 +4.936610000000000303e-02 2.041286000000000142e-02 2.049503381128415722e-02 1.407901999999999992e-04 +4.975639999999999924e-02 2.033532999999999868e-02 2.041720223102738013e-02 1.419589999999999871e-04 +5.014659999999999951e-02 2.025754999999999917e-02 2.033911230906939455e-02 1.429395000000000025e-04 +5.053679999999999978e-02 2.017930999999999891e-02 2.026056093967545588e-02 1.444802000000000130e-04 +5.092700000000000005e-02 2.010106999999999866e-02 2.018200928790250467e-02 1.396145999999999965e-04 +5.131730000000000319e-02 2.002090999999999940e-02 2.010155013096629265e-02 1.403297000000000127e-04 +5.170750000000000346e-02 1.994032999999999917e-02 2.002067715298700562e-02 1.399915999999999964e-04 +5.209769999999999679e-02 1.985984999999999834e-02 1.993989845014924878e-02 1.370638999999999939e-04 +5.248789999999999706e-02 1.977968000000000087e-02 1.985942781678529112e-02 1.392661999999999915e-04 +5.287809999999999733e-02 1.969950999999999994e-02 1.977894999785135763e-02 1.382088999999999974e-04 +5.326840000000000047e-02 1.961970000000000033e-02 1.969884494753471743e-02 1.374512000000000097e-04 +5.365860000000000074e-02 1.954010999999999942e-02 1.961894403765851447e-02 1.363897999999999998e-04 +5.404880000000000101e-02 1.946029000000000161e-02 1.953881890439435748e-02 1.397796999999999993e-04 +5.443900000000000128e-02 1.937890999999999919e-02 1.945712248120418097e-02 1.365279000000000113e-04 +5.482929999999999748e-02 1.929749999999999868e-02 1.937539567089395162e-02 1.324386999999999895e-04 +5.521949999999999775e-02 1.921609999999999985e-02 1.929367191118193728e-02 1.323242000000000027e-04 +5.560969999999999802e-02 1.913466999999999946e-02 1.921191688852882593e-02 1.303377000000000040e-04 +5.599989999999999829e-02 1.905325000000000074e-02 1.913016916381843346e-02 1.282578000000000038e-04 +5.639009999999999856e-02 1.897110999999999936e-02 1.904770657616991839e-02 1.292793999999999942e-04 +5.678029999999999883e-02 1.888897000000000145e-02 1.896524398853013937e-02 1.298956000000000036e-04 +5.717060000000000197e-02 1.880679999999999852e-02 1.888275375382969523e-02 1.257393000000000041e-04 +5.756080000000000224e-02 1.872464000000000073e-02 1.880027243959095373e-02 1.250178000000000010e-04 +5.795100000000000251e-02 1.864247999999999947e-02 1.871779171569387321e-02 1.245061999999999905e-04 +5.834120000000000278e-02 1.855971999999999830e-02 1.863470501784624325e-02 1.245713999999999979e-04 +5.873140000000000305e-02 1.847688999999999929e-02 1.855153903662197859e-02 1.268389000000000029e-04 +5.912160000000000332e-02 1.839404999999999860e-02 1.846836326163447603e-02 1.221787999999999965e-04 +5.951189999999999952e-02 1.831115999999999994e-02 1.838514979396218066e-02 1.202135999999999971e-04 +5.990209999999999979e-02 1.822829000000000116e-02 1.830195129957841102e-02 1.205508000000000046e-04 +6.029230000000000006e-02 1.814539000000000082e-02 1.821869290184508633e-02 1.199057999999999998e-04 +6.068250000000000033e-02 1.806249000000000049e-02 1.813541556231786436e-02 1.196926000000000066e-04 +6.107270000000000060e-02 1.797939999999999955e-02 1.805196406641292745e-02 1.215309000000000056e-04 +6.146290000000000087e-02 1.789556000000000133e-02 1.796779688773924438e-02 1.190594999999999940e-04 +6.185310000000000114e-02 1.781171000000000143e-02 1.788362838675769395e-02 1.213999000000000030e-04 +6.224330000000000140e-02 1.772804000000000046e-02 1.779963720853651801e-02 1.203586999999999967e-04 +6.263349999999999473e-02 1.764446999999999890e-02 1.771575186057049861e-02 1.177072999999999970e-04 +6.302380000000000482e-02 1.756090000000000081e-02 1.763187033451096147e-02 1.180958999999999953e-04 +6.341399999999999815e-02 1.747786999999999952e-02 1.754845157169272404e-02 1.175629999999999991e-04 +6.380420000000000536e-02 1.739484000000000169e-02 1.746504173716099256e-02 1.154626999999999948e-04 +6.419439999999999868e-02 1.731098999999999832e-02 1.738082843969512153e-02 1.158778999999999972e-04 +6.458460000000000589e-02 1.722631000000000162e-02 1.729581794735354089e-02 1.142516000000000023e-04 +6.497479999999999922e-02 1.714165000000000133e-02 1.721081475300145422e-02 1.179614000000000053e-04 +6.536500000000000643e-02 1.705813000000000121e-02 1.712698946708405071e-02 1.134273000000000041e-04 +6.575519999999999976e-02 1.697469999999999882e-02 1.704324910729296955e-02 1.138225000000000033e-04 +6.614539999999999309e-02 1.689097000000000168e-02 1.695918753140932336e-02 1.141267999999999997e-04 +6.653560000000000030e-02 1.680672000000000069e-02 1.687458106879960104e-02 1.139221999999999983e-04 +6.692579999999999363e-02 1.672246000000000149e-02 1.678996713824609799e-02 1.118457999999999990e-04 +6.731600000000000084e-02 1.663859000000000171e-02 1.670571967990813755e-02 1.111126999999999975e-04 +6.770619999999999417e-02 1.655480999999999966e-02 1.662155061615788049e-02 1.102739999999999949e-04 +6.809640000000000137e-02 1.647097000000000144e-02 1.653734311087814296e-02 1.108059999999999959e-04 +6.848659999999999470e-02 1.638695999999999903e-02 1.645301678229664152e-02 1.078567000000000002e-04 +6.887680000000000191e-02 1.630293000000000020e-02 1.636868222549060087e-02 1.079058999999999997e-04 +6.926699999999999524e-02 1.621883000000000005e-02 1.628425341034724777e-02 1.088324000000000051e-04 +6.965720000000000245e-02 1.613469000000000014e-02 1.619977826552451011e-02 1.054888000000000055e-04 +7.004739999999999578e-02 1.605060999999999988e-02 1.611537158983979146e-02 1.065766000000000055e-04 +7.043760000000000299e-02 1.596699999999999856e-02 1.603144840729229387e-02 1.057296000000000050e-04 +7.082779999999999632e-02 1.588339999999999891e-02 1.594752536272911123e-02 1.054971000000000034e-04 +7.121800000000000352e-02 1.579985999999999891e-02 1.586363310298559079e-02 1.064822999999999952e-04 +7.160819999999999685e-02 1.571635999999999866e-02 1.577975762019730860e-02 1.057505999999999968e-04 +7.199840000000000406e-02 1.563287000000000010e-02 1.569588376776566083e-02 1.060603000000000050e-04 +7.238859999999999739e-02 1.554901000000000026e-02 1.561168883189871688e-02 1.041179000000000014e-04 +7.277880000000000460e-02 1.546516000000000036e-02 1.552749659743643887e-02 1.048078999999999965e-04 +7.316899999999999793e-02 1.538139999999999993e-02 1.544340075760752481e-02 1.042733999999999968e-04 +7.355920000000000514e-02 1.529776000000000052e-02 1.535942896224578094e-02 1.006476000000000018e-04 +7.394939999999999847e-02 1.521410999999999944e-02 1.527544981134015649e-02 9.874159999999999665e-05 +7.433950000000000280e-02 1.513095000000000065e-02 1.519191771079473653e-02 9.938579999999999973e-05 +7.472969999999999613e-02 1.504783999999999983e-02 1.510843537956833502e-02 1.014812000000000000e-04 +7.511990000000000334e-02 1.496485999999999997e-02 1.502507671237522484e-02 1.014849999999999948e-04 +7.551009999999999667e-02 1.488215000000000024e-02 1.494201267182827461e-02 9.770799999999999353e-05 +7.590030000000000387e-02 1.479944000000000051e-02 1.485894279364298638e-02 9.800800000000000083e-05 +7.629049999999999720e-02 1.471698999999999924e-02 1.477614466245610186e-02 9.677200000000000114e-05 +7.668070000000000441e-02 1.463461999999999923e-02 1.469343688698711639e-02 9.775429999999999341e-05 +7.707089999999999774e-02 1.455221999999999939e-02 1.461068916925587283e-02 9.549610000000000206e-05 +7.746100000000000207e-02 1.446968000000000039e-02 1.452781386919502378e-02 9.513919999999999875e-05 +7.785119999999999540e-02 1.438712999999999971e-02 1.444491660580846104e-02 9.770649999999999702e-05 +7.824140000000000261e-02 1.430488000000000072e-02 1.436233045946878721e-02 9.529270000000000329e-05 +7.863159999999999594e-02 1.422283000000000054e-02 1.427993459091195065e-02 9.643389999999999840e-05 +7.902180000000000315e-02 1.414078000000000036e-02 1.419753652017903545e-02 9.496639999999999932e-05 +7.941199999999999648e-02 1.405880999999999971e-02 1.411522168003312601e-02 9.281660000000000171e-05 +7.980210000000000081e-02 1.397687000000000061e-02 1.403292925548272306e-02 9.309979999999999927e-05 +8.019229999999999414e-02 1.389486000000000020e-02 1.395056495690346829e-02 9.322759999999999489e-05 +8.058250000000000135e-02 1.381279000000000014e-02 1.386815161594708544e-02 9.204930000000000653e-05 +8.097269999999999468e-02 1.373072000000000008e-02 1.378573869526804677e-02 9.424679999999999354e-05 +8.136290000000000189e-02 1.364999000000000073e-02 1.370465831832146444e-02 9.418200000000000053e-05 +8.175300000000000622e-02 1.356939000000000062e-02 1.362370536413650204e-02 9.375499999999999936e-05 +8.214319999999999955e-02 1.348890999999999979e-02 1.354288902835608377e-02 9.299580000000000162e-05 +8.253340000000000676e-02 1.340868999999999915e-02 1.346234826006290942e-02 9.085749999999999534e-05 +8.292360000000000009e-02 1.332848000000000019e-02 1.338181495972442722e-02 9.101269999999999684e-05 +8.331370000000000442e-02 1.324830999999999925e-02 1.330131282203845904e-02 9.108979999999999649e-05 +8.370389999999999775e-02 1.316813000000000011e-02 1.322079307304159512e-02 8.960849999999999425e-05 +8.409410000000000496e-02 1.308808000000000020e-02 1.314039956958151616e-02 8.983630000000000586e-05 +8.448419999999999541e-02 1.300845999999999947e-02 1.306043089402079624e-02 8.921549999999999988e-05 +8.487440000000000262e-02 1.292882000000000059e-02 1.298044935338811892e-02 9.003020000000000413e-05 +8.526459999999999595e-02 1.284920999999999980e-02 1.290050074779767499e-02 8.868019999999999842e-05 +8.565480000000000316e-02 1.276963000000000056e-02 1.282057691750264697e-02 8.695680000000000065e-05 +8.604489999999999361e-02 1.269016000000000068e-02 1.274075932068990946e-02 8.962229999999999741e-05 +8.643510000000000082e-02 1.261128000000000075e-02 1.266155427410518668e-02 8.753039999999999801e-05 +8.682529999999999415e-02 1.253241000000000077e-02 1.258234905755320661e-02 8.837510000000000027e-05 +8.721539999999999848e-02 1.245405999999999944e-02 1.250367044019390361e-02 8.655230000000000139e-05 +8.760560000000000569e-02 1.237609999999999925e-02 1.242538581230035431e-02 8.655489999999999354e-05 +8.799579999999999902e-02 1.229815000000000075e-02 1.234710660171956485e-02 8.638530000000000113e-05 +8.838590000000000335e-02 1.222026000000000015e-02 1.226889996120519334e-02 8.725110000000000222e-05 +8.877609999999999668e-02 1.214234999999999967e-02 1.219067611926958825e-02 8.563389999999999338e-05 +8.916620000000000101e-02 1.206469999999999938e-02 1.211270039273322670e-02 8.539790000000000445e-05 +8.955639999999999434e-02 1.198735000000000078e-02 1.203501427124121712e-02 8.543029999999999418e-05 +8.994660000000000155e-02 1.191000000000000045e-02 1.195732803739646584e-02 8.535040000000000194e-05 +9.033670000000000588e-02 1.183300999999999971e-02 1.188002885239964126e-02 8.536150000000000595e-05 +9.072689999999999921e-02 1.175604000000000059e-02 1.180276148483007728e-02 8.610179999999999613e-05 +9.111700000000000355e-02 1.167938999999999956e-02 1.172579702615247216e-02 8.831479999999999680e-05 +9.150719999999999688e-02 1.160340999999999977e-02 1.164949131458491334e-02 8.455170000000000567e-05 +9.189730000000000121e-02 1.152744999999999985e-02 1.157319563153584932e-02 8.551700000000000133e-05 +9.228749999999999454e-02 1.145181999999999971e-02 1.149722739159510232e-02 8.333830000000000133e-05 +9.267770000000000175e-02 1.137631000000000059e-02 1.142138567875185425e-02 8.406990000000000632e-05 +9.306780000000000608e-02 1.130088999999999920e-02 1.134563681146314937e-02 8.385449999999999480e-05 +9.345799999999999941e-02 1.122578999999999938e-02 1.127024298802847997e-02 8.294990000000000349e-05 +9.384810000000000374e-02 1.115070999999999944e-02 1.119487558927098093e-02 8.268419999999999687e-05 +9.423829999999999707e-02 1.107604000000000054e-02 1.111988931460827952e-02 8.292019999999999935e-05 +9.462840000000000140e-02 1.100166999999999985e-02 1.104519137259795357e-02 8.294990000000000349e-05 +9.501859999999999473e-02 1.092732000000000078e-02 1.097051175932093128e-02 8.221800000000000462e-05 +9.540869999999999906e-02 1.085367999999999923e-02 1.089657276738932069e-02 8.236279999999999686e-05 +9.579890000000000627e-02 1.078001999999999953e-02 1.082260749246298157e-02 8.193070000000000485e-05 +9.618899999999999673e-02 1.070659000000000020e-02 1.074888658763396669e-02 8.178810000000000387e-05 +9.657920000000000393e-02 1.063335999999999969e-02 1.067537224943441651e-02 8.181649999999999839e-05 +9.696929999999999439e-02 1.056015000000000079e-02 1.060188360391375699e-02 8.147700000000000612e-05 +9.735939999999999872e-02 1.048799000000000085e-02 1.052939438393979237e-02 8.278839999999999496e-05 +9.774960000000000593e-02 1.041589000000000056e-02 1.045696799431052211e-02 8.139279999999999768e-05 +9.813969999999999638e-02 1.034379000000000028e-02 1.038454969622005369e-02 8.078309999999999571e-05 +9.852990000000000359e-02 1.027162000000000040e-02 1.031209247583799934e-02 8.075579999999999683e-05 +9.891999999999999404e-02 1.019948000000000035e-02 1.023965541846634693e-02 8.166849999999999912e-05 +9.931009999999999838e-02 1.012806999999999943e-02 1.016795739511015123e-02 8.070110000000000564e-05 +9.970030000000000558e-02 1.005682999999999924e-02 1.009642937686532418e-02 8.032129999999999955e-05 +1.000899999999999984e-01 9.985924999999999663e-03 1.002523421813927840e-02 7.978789999999999548e-05 +1.004810000000000009e-01 9.915571999999999234e-03 9.954599718145421341e-03 8.159719999999999864e-05 +1.008710000000000023e-01 9.845411999999999567e-03 9.884159526810371368e-03 7.986729999999999340e-05 +1.012610000000000038e-01 9.775653999999999941e-03 9.814106784750434298e-03 7.986729999999999340e-05 +1.016510000000000052e-01 9.706091000000000024e-03 9.744242462897068080e-03 7.918749999999999357e-05 +1.020410000000000067e-01 9.636580999999999966e-03 9.674432580379145954e-03 7.900750000000000546e-05 +1.024309999999999943e-01 9.567463999999999483e-03 9.605033723073736970e-03 8.018249999999999336e-05 +1.028209999999999957e-01 9.498345000000000052e-03 9.535632953108497123e-03 7.885409999999999436e-05 +1.032109999999999972e-01 9.429552000000000767e-03 9.466559960887240693e-03 7.842439999999999405e-05 +1.036019999999999996e-01 9.360867999999999620e-03 9.397596168922404958e-03 7.870160000000000556e-05 +1.039920000000000011e-01 9.292359000000000036e-03 9.328807703208418303e-03 7.832429999999999818e-05 +1.043820000000000026e-01 9.225139999999999729e-03 9.261315498934339155e-03 7.877779999999999646e-05 +1.047720000000000040e-01 9.157956000000000082e-03 9.193857792025946288e-03 7.797689999999999537e-05 +1.051620000000000055e-01 9.090878999999999627e-03 9.126505445784194315e-03 7.807569999999999517e-05 +1.055520000000000069e-01 9.023957000000000853e-03 9.059303304327512296e-03 7.775609999999999911e-05 +1.059419999999999945e-01 8.957042000000000129e-03 8.992108492832507835e-03 7.758559999999999795e-05 +1.063319999999999960e-01 8.890629000000000587e-03 8.925430038632565719e-03 7.768290000000000123e-05 +1.067229999999999984e-01 8.824138999999999941e-03 8.858676837224304865e-03 7.763419999999999609e-05 +1.071129999999999999e-01 8.758095999999999937e-03 8.792371405876674498e-03 7.712839999999999833e-05 +1.075030000000000013e-01 8.692715999999999638e-03 8.726727138619471899e-03 7.724800000000000308e-05 +1.078930000000000028e-01 8.627343000000000858e-03 8.661089610923416310e-03 7.665580000000000560e-05 +1.082830000000000042e-01 8.562518999999999339e-03 8.596002692194178868e-03 7.663239999999999494e-05 +1.086730000000000057e-01 8.497892999999999419e-03 8.531113124543358156e-03 7.804370000000000632e-05 +1.090629999999999933e-01 8.433404000000000247e-03 8.466362483832982205e-03 7.651570000000000333e-05 +1.094529999999999947e-01 8.369599999999999609e-03 8.402303593678224491e-03 7.670269999999999324e-05 +1.098439999999999972e-01 8.305626999999999455e-03 8.338075339093179231e-03 7.598570000000000671e-05 +1.102339999999999987e-01 8.242275000000000296e-03 8.274459938656674590e-03 7.580390000000000110e-05 +1.106240000000000001e-01 8.179225000000000037e-03 8.211141233379041712e-03 7.644600000000000635e-05 +1.110140000000000016e-01 8.116197000000000342e-03 8.147845584243957814e-03 7.546660000000000011e-05 +1.114040000000000030e-01 8.053737999999999522e-03 8.085144385881561785e-03 7.610000000000000661e-05 +1.117940000000000045e-01 7.991279999999999911e-03 8.022444817826972913e-03 7.703320000000000642e-05 +1.121840000000000059e-01 7.929313000000000333e-03 7.960237384892022983e-03 7.653899999999999343e-05 +1.125739999999999935e-01 7.867884000000000336e-03 7.898568572810989522e-03 7.756130000000000565e-05 +1.129639999999999950e-01 7.806457999999999627e-03 7.836901811329656520e-03 7.686739999999999524e-05 +1.133549999999999974e-01 7.745417000000000100e-03 7.775603703224632027e-03 7.765850000000000194e-05 +1.137449999999999989e-01 7.684589999999999824e-03 7.714518792126701108e-03 7.802629999999999527e-05 +1.141350000000000003e-01 7.624022000000000264e-03 7.653699613847527415e-03 7.729590000000000646e-05 +1.145250000000000018e-01 7.563936999999999744e-03 7.593375507885479631e-03 7.768290000000000123e-05 +1.149150000000000033e-01 7.503852000000000091e-03 7.533051992232098377e-03 7.729590000000000646e-05 +1.153050000000000047e-01 7.444553999999999928e-03 7.473524223322518083e-03 7.751289999999999439e-05 +1.156950000000000062e-01 7.385476999999999667e-03 7.414221451600616630e-03 7.842439999999999405e-05 +1.160849999999999937e-01 7.326397999999999591e-03 7.354915841988299459e-03 7.736799999999999515e-05 +1.164749999999999952e-01 7.267306999999999760e-03 7.295595498027451081e-03 7.770719999999999353e-05 +1.168649999999999967e-01 7.208209000000000144e-03 7.236268106446680157e-03 7.758559999999999795e-05 +1.172549999999999981e-01 7.150253999999999985e-03 7.178079956390027473e-03 7.703649999999999333e-05 +1.176460000000000006e-01 7.092756999999999881e-03 7.120347429028719630e-03 7.822460000000000318e-05 +1.180360000000000020e-01 7.035422000000000377e-03 7.062779703614394979e-03 7.773169999999999982e-05 +1.184260000000000035e-01 6.978250000000000078e-03 7.005383064547699226e-03 7.727630000000000414e-05 +1.188160000000000049e-01 6.921072000000000335e-03 6.947981178150877506e-03 7.815000000000000223e-05 +1.192060000000000064e-01 6.864224999999999841e-03 6.890917802534510077e-03 7.760990000000000380e-05 +1.195959999999999940e-01 6.807675000000000011e-03 6.834157822465209212e-03 7.805100000000000200e-05 +1.199859999999999954e-01 6.751117000000000055e-03 6.777390964746413289e-03 7.734399999999999674e-05 +1.203759999999999969e-01 6.695979000000000340e-03 6.722025499520075005e-03 7.768290000000000123e-05 +1.207659999999999983e-01 6.640886999999999832e-03 6.666705536198270013e-03 7.815000000000000223e-05 +1.211559999999999998e-01 6.586002999999999996e-03 6.611601297508906336e-03 7.724800000000000308e-05 +1.215460000000000013e-01 6.531437000000000180e-03 6.556825657495615745e-03 7.800160000000000210e-05 +1.219370000000000037e-01 6.476736000000000298e-03 6.501915450311195922e-03 7.715230000000000330e-05 +1.223270000000000052e-01 6.422755999999999951e-03 6.447712721281343221e-03 7.758559999999999795e-05 +1.227170000000000066e-01 6.368890999999999789e-03 6.393623928437020161e-03 7.797689999999999537e-05 +1.231069999999999942e-01 6.315285000000000343e-03 6.339801163925852870e-03 7.736799999999999515e-05 +1.234969999999999957e-01 6.262353000000000017e-03 6.286673901904199328e-03 7.753709999999999325e-05 +1.238869999999999971e-01 6.209428000000000343e-03 6.233553799906901317e-03 7.729590000000000646e-05 +1.242769999999999986e-01 6.156970000000000325e-03 6.180889476861710884e-03 7.729590000000000646e-05 +1.246670000000000000e-01 6.104708999999999830e-03 6.128417887546812673e-03 7.763419999999999609e-05 +1.250570000000000015e-01 6.052553999999999920e-03 6.076055362815024029e-03 7.689040000000000501e-05 +1.254470000000000030e-01 6.000977999999999694e-03 6.024284341728588545e-03 7.717819999999999910e-05 +1.258370000000000044e-01 5.949409000000000121e-03 5.972520926627121647e-03 7.727189999999999450e-05 +1.262270000000000059e-01 5.898373999999999631e-03 5.921288942881298457e-03 7.682040000000000060e-05 +1.266170000000000073e-01 5.847736000000000080e-03 5.870453452929540530e-03 7.729590000000000646e-05 +1.270070000000000088e-01 5.797102999999999631e-03 5.819622487404882936e-03 7.667920000000000270e-05 +1.273980000000000112e-01 5.747145999999999852e-03 5.769479612444817183e-03 7.751289999999999439e-05 +1.277880000000000127e-01 5.697312000000000383e-03 5.719460098137257864e-03 7.756130000000000565e-05 +1.281779999999999864e-01 5.647738000000000237e-03 5.669694945821524695e-03 7.642279999999999614e-05 +1.285679999999999878e-01 5.598477999999999614e-03 5.620236180117848396e-03 7.672409999999999950e-05 +1.289579999999999893e-01 5.549219999999999674e-03 5.570779635037102426e-03 7.677369999999999984e-05 +1.293479999999999908e-01 5.500625000000000306e-03 5.521996757653743933e-03 7.691469999999999731e-05 +1.297379999999999922e-01 5.452114000000000092e-03 5.473299762841953714e-03 7.660759999999999478e-05 +1.301279999999999937e-01 5.403907999999999802e-03 5.424907309546476639e-03 7.627970000000000085e-05 +1.305179999999999951e-01 5.356323999999999981e-03 5.377137328982770428e-03 7.783710000000000054e-05 +1.309079999999999966e-01 5.308740000000000160e-03 5.329367768746692824e-03 7.771550000000000495e-05 +1.312979999999999980e-01 5.261789000000000223e-03 5.282231884939614722e-03 7.905880000000000275e-05 +1.316879999999999995e-01 5.215028000000000025e-03 5.235286909652696503e-03 7.852990000000000176e-05 +1.320780000000000010e-01 5.168378999999999855e-03 5.188456452459017446e-03 7.853420000000000442e-05 +1.324680000000000024e-01 5.122149999999999793e-03 5.142055266766147777e-03 7.963070000000000315e-05 +1.328580000000000039e-01 5.075912000000000132e-03 5.095645714863236524e-03 8.081219999999999853e-05 +1.332480000000000053e-01 5.030459000000000146e-03 5.050013739781279536e-03 7.934110000000000511e-05 +1.336380000000000068e-01 4.985449999999999778e-03 5.004820076744203976e-03 7.866659999999999658e-05 +1.340280000000000082e-01 4.940476000000000069e-03 4.959663593969062742e-03 7.825859999999999642e-05 +1.344180000000000097e-01 4.895980999999999771e-03 4.915012916732218635e-03 7.892710000000000535e-05 +1.348080000000000112e-01 4.851486000000000340e-03 4.870362377487781960e-03 7.965670000000000595e-05 +1.351990000000000136e-01 4.807284999999999683e-03 4.825995781660621237e-03 7.846529999999999564e-05 +1.355889999999999873e-01 4.763584999999999868e-03 4.782119847129753971e-03 7.989869999999999449e-05 +1.359789999999999888e-01 4.719883999999999712e-03 4.738244220564536216e-03 8.006769999999999914e-05 +1.363689999999999902e-01 4.676684000000000396e-03 4.694881402238172165e-03 7.941140000000000340e-05 +1.367589999999999917e-01 4.633507000000000250e-03 4.651541923787324770e-03 7.987600000000000570e-05 +1.371489999999999931e-01 4.590636999999999843e-03 4.608508992499429739e-03 7.922749999999999996e-05 +1.375389999999999946e-01 4.548243999999999898e-03 4.565950688357889829e-03 7.928880000000000563e-05 +1.379289999999999961e-01 4.505850999999999953e-03 4.523392384216710742e-03 7.876470000000000162e-05 +1.383189999999999975e-01 4.464336000000000172e-03 4.481712575352084227e-03 7.825370000000000600e-05 +1.387089999999999990e-01 4.423024999999999700e-03 4.440235406367055299e-03 7.813860000000000434e-05 +1.390990000000000004e-01 4.381847000000000346e-03 4.398893938955567018e-03 7.822839999999999796e-05 +1.394890000000000019e-01 4.341051999999999758e-03 4.357940383301032658e-03 7.697190000000000076e-05 +1.398790000000000033e-01 4.300258000000000379e-03 4.316988287989037033e-03 7.687109999999999657e-05 +1.402690000000000048e-01 4.259889000000000210e-03 4.276468154694128450e-03 7.571339999999999917e-05 +1.406590000000000062e-01 4.219719000000000247e-03 4.236148722935041168e-03 7.414069999999999981e-05 +1.410490000000000077e-01 4.179676999999999698e-03 4.195957146193621604e-03 7.480040000000000300e-05 +1.414390000000000092e-01 4.140493999999999772e-03 4.156619114688175042e-03 7.588220000000000338e-05 +1.418290000000000106e-01 4.101312000000000187e-03 4.117283021458397786e-03 7.531849999999999385e-05 +1.422190000000000121e-01 4.062477999999999576e-03 4.078298015633048841e-03 7.607619999999999509e-05 +1.426090000000000135e-01 4.023918999999999933e-03 4.039590265035316631e-03 7.624829999999999976e-05 +1.429989999999999872e-01 3.985362000000000106e-03 4.000884735052255732e-03 7.414020000000000549e-05 +1.433889999999999887e-01 3.947797000000000285e-03 3.963170753358553726e-03 7.331370000000000248e-05 +1.437789999999999901e-01 3.910231999999999597e-03 3.925457022029048108e-03 7.485140000000000641e-05 +1.441689999999999916e-01 3.872865999999999982e-03 3.887943785065704327e-03 7.537839999999999645e-05 +1.445589999999999931e-01 3.835751000000000091e-03 3.850681680258943255e-03 7.395579999999999418e-05 +1.449489999999999945e-01 3.798642000000000078e-03 3.813427181385629482e-03 7.325290000000000469e-05 +1.453389999999999960e-01 3.762179000000000044e-03 3.776822312248529542e-03 7.403250000000000650e-05 +1.457289999999999974e-01 3.725818000000000099e-03 3.740321278818892030e-03 7.324929999999999679e-05 +1.461189999999999989e-01 3.689692999999999966e-03 3.704052713725200681e-03 7.281560000000000127e-05 +1.465090000000000003e-01 3.654090000000000029e-03 3.668303033334851301e-03 7.394060000000000151e-05 +1.468990000000000018e-01 3.618497000000000033e-03 3.632561889120675819e-03 7.294230000000000126e-05 +1.472890000000000033e-01 3.583296000000000137e-03 3.597226939049489453e-03 7.582549999999999425e-05 +1.476790000000000047e-01 3.548233000000000029e-03 3.562034062721425656e-03 7.260239999999999457e-05 +1.480690000000000062e-01 3.513329999999999839e-03 3.526999471655489250e-03 7.098020000000000186e-05 +1.484590000000000076e-01 3.479165000000000209e-03 3.492691353754984532e-03 7.115269999999999386e-05 +1.488490000000000091e-01 3.444998999999999804e-03 3.458383713795992569e-03 7.203430000000000250e-05 +1.492390000000000105e-01 3.411118999999999783e-03 3.424373427789681370e-03 7.151860000000000335e-05 +1.496290000000000120e-01 3.377404999999999990e-03 3.390536762628489856e-03 6.975539999999999963e-05 +1.500190000000000135e-01 3.343750999999999841e-03 3.356759908741979194e-03 7.207830000000000411e-05 +1.504089999999999872e-01 3.311071999999999974e-03 3.323938945049851115e-03 7.000219999999999869e-05 +1.507989999999999886e-01 3.278390999999999857e-03 3.291116828980062293e-03 6.909410000000000649e-05 +1.511889999999999901e-01 3.245917999999999980e-03 3.258518306957129598e-03 7.074379999999999850e-05 +1.515789999999999915e-01 3.213656999999999972e-03 3.226147669106086465e-03 6.995620000000000625e-05 +1.519689999999999930e-01 3.181401999999999841e-03 3.193782586540288605e-03 6.967729999999999779e-05 +1.523589999999999944e-01 3.149945999999999822e-03 3.162195504152960288e-03 7.165079999999999507e-05 +1.527489999999999959e-01 3.118567999999999948e-03 3.130683740703821619e-03 7.121440000000000040e-05 +1.531389999999999973e-01 3.087377000000000090e-03 3.099362059268761377e-03 7.058950000000000576e-05 +1.535289999999999988e-01 3.056539000000000044e-03 3.068397982216640105e-03 7.117140000000000098e-05 +1.539190000000000003e-01 3.025705999999999968e-03 3.037440050754521004e-03 6.978990000000000074e-05 +1.543090000000000017e-01 2.995299999999999786e-03 3.006917707629095521e-03 7.064020000000000173e-05 +1.546990000000000032e-01 2.965016000000000006e-03 2.976521192391699128e-03 7.175539999999999404e-05 +1.550890000000000046e-01 2.934939000000000124e-03 2.946330042971181343e-03 7.090020000000000263e-05 +1.554790000000000061e-01 2.905603999999999999e-03 2.916875599319826602e-03 6.824010000000000416e-05 +1.558690000000000075e-01 2.876259999999999841e-03 2.887412789456781422e-03 7.015070000000000582e-05 +1.562590000000000090e-01 2.847224000000000196e-03 2.858261355699302873e-03 7.053629999999999753e-05 +1.566490000000000105e-01 2.818324999999999997e-03 2.829248270715712353e-03 6.910960000000000659e-05 +1.570390000000000119e-01 2.789516000000000131e-03 2.800324278336249317e-03 7.000039999999999474e-05 +1.574290000000000134e-01 2.761415000000000137e-03 2.772100292214679969e-03 6.826470000000000389e-05 +1.578189999999999871e-01 2.733315000000000050e-03 2.743876476065622814e-03 6.868950000000000023e-05 +1.582089999999999885e-01 2.705531000000000168e-03 2.715979443758927483e-03 6.784339999999999490e-05 +1.585979999999999890e-01 2.678091000000000065e-03 2.688434681094967280e-03 6.768259999999999468e-05 +1.589879999999999904e-01 2.650581999999999851e-03 2.660821708235847852e-03 7.048529999999999412e-05 +1.593779999999999919e-01 2.623696999999999922e-03 2.633825199887842184e-03 6.678670000000000212e-05 +1.597679999999999934e-01 2.596838000000000184e-03 2.606854667905693888e-03 6.806260000000000120e-05 +1.601579999999999948e-01 2.570183000000000190e-03 2.580092981264940372e-03 6.871039999999999863e-05 +1.605479999999999963e-01 2.543807999999999833e-03 2.553619222140252835e-03 6.722320000000000378e-05 +1.609379999999999977e-01 2.517438999999999788e-03 2.527151890940786474e-03 6.860290000000000008e-05 +1.613279999999999992e-01 2.491541999999999802e-03 2.501157969061430005e-03 6.642819999999999530e-05 +1.617180000000000006e-01 2.465730000000000179e-03 2.475251014726331193e-03 6.829539999999999667e-05 +1.621080000000000021e-01 2.440117999999999801e-03 2.449537867086657578e-03 6.741309999999999328e-05 +1.624980000000000036e-01 2.414995000000000076e-03 2.424297177188450281e-03 6.718890000000000311e-05 +1.628880000000000050e-01 2.389871999999999917e-03 2.399057527302398806e-03 6.741690000000000161e-05 +1.632780000000000065e-01 2.365219999999999909e-03 2.374304457557258833e-03 6.624600000000000237e-05 +1.636680000000000079e-01 2.340750999999999856e-03 2.349741823667695923e-03 6.692020000000000347e-05 +1.640580000000000094e-01 2.316367000000000165e-03 2.325261849698087840e-03 6.688960000000000413e-05 +1.644480000000000108e-01 2.292428000000000000e-03 2.301221692539575595e-03 6.762290000000000607e-05 +1.648380000000000123e-01 2.268488999999999835e-03 2.277181057447325192e-03 6.835449999999999751e-05 +1.652280000000000137e-01 2.244897999999999946e-03 2.253496645638895877e-03 6.869079999999999631e-05 +1.656169999999999864e-01 2.221619000000000198e-03 2.230129570371323054e-03 6.746150000000000454e-05 +1.660069999999999879e-01 2.198294999999999850e-03 2.206717817169590208e-03 6.907780000000000463e-05 +1.663969999999999894e-01 2.175713000000000126e-03 2.184036898316521740e-03 6.889120000000000204e-05 +1.667869999999999908e-01 2.153129000000000154e-03 2.161355251162463010e-03 6.772520000000000677e-05 +1.671769999999999923e-01 2.130687000000000032e-03 2.138817610188813428e-03 7.153609999999999429e-05 +1.675669999999999937e-01 2.108416000000000110e-03 2.116453762942960889e-03 6.937949999999999532e-05 +1.679569999999999952e-01 2.086145000000000187e-03 2.094090393634430881e-03 6.916460000000000522e-05 +1.683469999999999966e-01 2.064514999999999840e-03 2.072371495021373220e-03 6.855229999999999755e-05 +1.687369999999999981e-01 2.042971000000000197e-03 2.050739599191951038e-03 6.994760000000000094e-05 +1.691269999999999996e-01 2.021560999999999844e-03 2.029243498964222823e-03 6.755549999999999381e-05 +1.695170000000000010e-01 2.000431999999999905e-03 2.008034416458182297e-03 6.868370000000000107e-05 +1.699070000000000025e-01 1.979303999999999873e-03 1.986825924251870573e-03 6.795630000000000529e-05 +1.702970000000000039e-01 1.958590999999999978e-03 1.966031270102142937e-03 6.902040000000000074e-05 +1.706860000000000044e-01 1.938052000000000004e-03 1.945411555716657342e-03 6.863630000000000555e-05 +1.710760000000000058e-01 1.917571999999999983e-03 1.924849603339735555e-03 6.846359999999999957e-05 +1.714660000000000073e-01 1.897524000000000042e-03 1.904718635712856394e-03 6.896460000000000036e-05 +1.718560000000000088e-01 1.877467000000000068e-03 1.884578877799051580e-03 6.800080000000000122e-05 +1.722460000000000102e-01 1.857825000000000015e-03 1.864854490435936644e-03 6.760270000000000244e-05 +1.726360000000000117e-01 1.838407000000000019e-03 1.845355983491966381e-03 6.726500000000000057e-05 +1.730260000000000131e-01 1.819022000000000000e-03 1.825890163223910885e-03 6.831339999999999548e-05 +1.734159999999999868e-01 1.799957999999999940e-03 1.806752033628942307e-03 6.840200000000000002e-05 +1.738059999999999883e-01 1.780894000000000097e-03 1.787613313724997332e-03 6.856929999999999417e-05 +1.741959999999999897e-01 1.762116999999999893e-03 1.768766170387312507e-03 6.971699999999999675e-05 +1.745849999999999902e-01 1.743671000000000023e-03 1.750251885518847462e-03 6.820779999999999432e-05 +1.749749999999999917e-01 1.725175000000000016e-03 1.731688090524433543e-03 6.736290000000000517e-05 +1.753649999999999931e-01 1.707289000000000072e-03 1.713723014906770809e-03 6.912970000000000323e-05 +1.757549999999999946e-01 1.689441999999999984e-03 1.695796694041889997e-03 6.874500000000000673e-05 +1.761449999999999960e-01 1.671713999999999926e-03 1.677992205595691842e-03 6.779379999999999456e-05 +1.765349999999999975e-01 1.654202999999999923e-03 1.660411151334141057e-03 6.832499999999999381e-05 +1.769249999999999989e-01 1.636686000000000044e-03 1.642822937097313321e-03 6.812130000000000117e-05 +1.773150000000000004e-01 1.619659999999999980e-03 1.625720564719550188e-03 6.800459999999999600e-05 +1.777050000000000018e-01 1.602738999999999943e-03 1.608721818498840559e-03 6.811849999999999503e-05 +1.780940000000000023e-01 1.585943999999999939e-03 1.591854310222239359e-03 6.869170000000000506e-05 +1.784840000000000038e-01 1.569340000000000014e-03 1.575195460186004041e-03 6.833049999999999909e-05 +1.788740000000000052e-01 1.552743999999999999e-03 1.558543965727642294e-03 6.938700000000000499e-05 +1.792640000000000067e-01 1.536533000000000083e-03 1.542272193801553278e-03 6.898549999999999875e-05 +1.796540000000000081e-01 1.520524000000000094e-03 1.526199178039503064e-03 6.900180000000000061e-05 +1.800440000000000096e-01 1.504544000000000020e-03 1.510156012833899653e-03 6.884539999999999649e-05 +1.804340000000000110e-01 1.488862000000000085e-03 1.494410031161222960e-03 6.918419999999999399e-05 +1.808240000000000125e-01 1.473188999999999966e-03 1.478673147738473375e-03 6.781660000000000390e-05 +1.812130000000000130e-01 1.457792000000000072e-03 1.463214792436034768e-03 6.849660000000000416e-05 +1.816029999999999867e-01 1.442556999999999911e-03 1.447921733548443793e-03 6.983220000000000540e-05 +1.819929999999999881e-01 1.427321000000000058e-03 1.432627606424102243e-03 6.843779999999999720e-05 +1.823829999999999896e-01 1.412658999999999920e-03 1.417905559205587989e-03 6.880020000000000580e-05 +1.827729999999999910e-01 1.398004999999999907e-03 1.403192696987993300e-03 6.995859999999999796e-05 +1.831629999999999925e-01 1.383488999999999899e-03 1.388617841926898617e-03 6.783569999999999834e-05 +1.835529999999999939e-01 1.369161999999999940e-03 1.374233844136508618e-03 6.942659999999999695e-05 +1.839419999999999944e-01 1.354876999999999991e-03 1.359891272603341384e-03 6.821979999999999353e-05 +1.843319999999999959e-01 1.340961000000000106e-03 1.345919825199748052e-03 6.796620000000000667e-05 +1.847219999999999973e-01 1.327122999999999931e-03 1.332026637582098489e-03 6.796980000000000101e-05 +1.851119999999999988e-01 1.313382000000000091e-03 1.318232867819654807e-03 6.891140000000000568e-05 +1.855020000000000002e-01 1.299879000000000095e-03 1.304681860970273184e-03 6.910809999999999653e-05 +1.858920000000000017e-01 1.286381999999999977e-03 1.291138014099025167e-03 6.863720000000000075e-05 +1.862820000000000031e-01 1.273152000000000051e-03 1.277861459579325261e-03 6.880290000000000494e-05 +1.866710000000000036e-01 1.260058000000000098e-03 1.264720936661419860e-03 6.981570000000000310e-05 +1.870610000000000051e-01 1.246998000000000030e-03 1.251613997147216738e-03 6.799179999999999504e-05 +1.874510000000000065e-01 1.234293999999999929e-03 1.238858164487802302e-03 6.890579999999999340e-05 +1.878410000000000080e-01 1.221597999999999955e-03 1.226111462067930357e-03 6.919630000000000019e-05 +1.882310000000000094e-01 1.209161999999999954e-03 1.213626342147920874e-03 6.872659999999999349e-05 +1.886210000000000109e-01 1.196896999999999935e-03 1.201313942123454848e-03 6.908550000000000118e-05 +1.890100000000000113e-01 1.184680000000000021e-03 1.189048940457348991e-03 6.897539999999999694e-05 +1.894000000000000128e-01 1.172826999999999984e-03 1.177152247689267637e-03 6.886089999999999659e-05 +1.897899999999999865e-01 1.160982999999999980e-03 1.165264483203993368e-03 7.001879999999999443e-05 +1.901799999999999879e-01 1.149257000000000100e-03 1.153495245985217696e-03 6.856119999999999673e-05 +1.905699999999999894e-01 1.137664999999999944e-03 1.141860735426192134e-03 7.087970000000000511e-05 +1.909599999999999909e-01 1.126073000000000005e-03 1.130226224863439734e-03 6.940489999999999680e-05 +1.913489999999999913e-01 1.114895000000000079e-03 1.119007555244638666e-03 6.862689999999999849e-05 +1.917389999999999928e-01 1.103739000000000066e-03 1.107810294446727288e-03 6.882030000000000244e-05 +1.921289999999999942e-01 1.092667000000000074e-03 1.096697523321008266e-03 6.936189999999999739e-05 +1.925189999999999957e-01 1.081785999999999946e-03 1.085779477513103969e-03 6.912439999999999839e-05 +1.929089999999999971e-01 1.070898999999999940e-03 1.074854246105282291e-03 6.976779999999999972e-05 +1.932979999999999976e-01 1.060366999999999994e-03 1.064284449820336338e-03 7.027020000000000358e-05 +1.936879999999999991e-01 1.049908999999999912e-03 1.053787282946540369e-03 7.000329999999999433e-05 +1.940780000000000005e-01 1.039502999999999998e-03 1.043343662972487200e-03 6.885150000000000309e-05 +1.944680000000000020e-01 1.029313000000000016e-03 1.033115555796461837e-03 6.894200000000000501e-05 +1.948580000000000034e-01 1.019123999999999941e-03 1.022888658603138181e-03 7.096749999999999435e-05 +1.952470000000000039e-01 1.009179000000000022e-03 1.012909300069499835e-03 6.915410000000000253e-05 +1.956370000000000053e-01 9.993420000000000688e-04 1.003038293878567248e-03 6.918809999999999576e-05 +1.960270000000000068e-01 9.895248999999999625e-04 9.931879096075624247e-04 6.942610000000000263e-05 +1.964170000000000083e-01 9.800364000000000211e-04 9.836718300761993806e-04 6.848079999999999662e-05 +1.968070000000000097e-01 9.705494999999999747e-04 9.741573808560532012e-04 6.902649999999999379e-05 +1.971960000000000102e-01 9.612031000000000264e-04 9.647777774550493363e-04 6.887369999999999755e-05 +1.975860000000000116e-01 9.519647999999999695e-04 9.555009113983180041e-04 6.977089999999999974e-05 +1.979760000000000131e-01 9.427205999999999607e-04 9.462180377396460160e-04 7.010900000000000248e-05 +1.983659999999999868e-01 9.338194999999999656e-04 9.372868192542283244e-04 7.041399999999999363e-05 +1.987559999999999882e-01 9.249433000000000047e-04 9.283810851779537073e-04 6.937639999999999530e-05 +1.991449999999999887e-01 9.161593999999999684e-04 9.195695964686032052e-04 6.894009999999999407e-05 +1.995349999999999902e-01 9.074693000000000429e-04 9.108550750942090825e-04 6.926260000000000327e-05 +1.999249999999999916e-01 8.987793999999999688e-04 9.021406917109429241e-04 6.999849999999999735e-05 +2.003149999999999931e-01 8.904117000000000455e-04 8.937485997501862364e-04 6.912399999999999751e-05 +2.007049999999999945e-01 8.821229000000000309e-04 8.854354138798793249e-04 7.009270000000000062e-05 +2.010939999999999950e-01 8.739231000000000082e-04 8.772095543280502552e-04 6.913019999999999755e-05 +2.014839999999999964e-01 8.659477999999999967e-04 8.692034264066068883e-04 6.921679999999999771e-05 +2.018739999999999979e-01 8.579660999999999711e-04 8.611907288178628989e-04 6.942480000000000656e-05 +2.022639999999999993e-01 8.501022000000000006e-04 8.533009988635987999e-04 7.077370000000000308e-05 +2.026529999999999998e-01 8.423226999999999711e-04 8.454982493515854694e-04 6.875199999999999497e-05 +2.030430000000000013e-01 8.345451000000000475e-04 8.376976063981077117e-04 6.847060000000000136e-05 +2.034330000000000027e-01 8.269878000000000038e-04 8.301202155110144710e-04 6.866130000000000616e-05 +2.038230000000000042e-01 8.194364000000000205e-04 8.225486622492074983e-04 6.901299999999999807e-05 +2.042120000000000046e-01 8.120967999999999460e-04 8.151883666480172876e-04 6.848079999999999662e-05 +2.046020000000000061e-01 8.048926999999999991e-04 8.079628113859400581e-04 6.878680000000000352e-05 +2.049920000000000075e-01 7.976956000000000541e-04 8.007444161019963268e-04 6.940600000000000599e-05 +2.053820000000000090e-01 7.906798000000000255e-04 7.937056473873297432e-04 6.942829999999999390e-05 +2.057710000000000095e-01 7.836861000000000521e-04 7.866890986928393040e-04 6.880680000000000671e-05 +2.061610000000000109e-01 7.767628000000000163e-04 7.797459928857612211e-04 6.865040000000000259e-05 +2.065510000000000124e-01 7.699818999999999674e-04 7.729502003901422011e-04 6.917379999999999829e-05 +2.069410000000000138e-01 7.631924000000000216e-04 7.661458717190005631e-04 6.926610000000000416e-05 +2.073299999999999865e-01 7.566879000000000096e-04 7.596194446442036240e-04 6.823499999999999976e-05 +2.077199999999999880e-01 7.502236000000000451e-04 7.531319456828068941e-04 6.853699999999999788e-05 +2.081099999999999894e-01 7.437897000000000388e-04 7.466770908252821437e-04 6.796980000000000101e-05 +2.084999999999999909e-01 7.374404999999999676e-04 7.403128179979843668e-04 6.834640000000000008e-05 +2.088889999999999914e-01 7.311061999999999899e-04 7.339634608725748239e-04 6.839109999999999645e-05 +2.092789999999999928e-01 7.249501999999999899e-04 7.277898952841397456e-04 6.781660000000000390e-05 +2.096689999999999943e-01 7.188729999999999730e-04 7.216941500352032747e-04 6.851600000000000604e-05 +2.100589999999999957e-01 7.128269000000000445e-04 7.156302823263786155e-04 6.842769999999999539e-05 +2.104479999999999962e-01 7.069889999999999898e-04 7.097783818959492105e-04 6.902180000000000381e-05 +2.108379999999999976e-01 7.011352999999999683e-04 7.039106355198940625e-04 6.835980000000000235e-05 +2.112279999999999991e-01 6.953802999999999937e-04 6.981409677519885516e-04 6.792310000000000025e-05 +2.116180000000000005e-01 6.896969999999999663e-04 6.924425967012601615e-04 6.844879999999999422e-05 +2.120070000000000010e-01 6.840299000000000422e-04 6.867603795611481547e-04 6.762110000000000213e-05 +2.123970000000000025e-01 6.785852999999999894e-04 6.813038585791437269e-04 6.755829999999999995e-05 +2.127870000000000039e-01 6.731420999999999802e-04 6.758486855729307257e-04 6.856259999999999980e-05 +2.131760000000000044e-01 6.677984999999999997e-04 6.704913894465928687e-04 6.912399999999999751e-05 +2.135660000000000058e-01 6.625446000000000330e-04 6.652216110399162760e-04 7.012630000000000653e-05 +2.139560000000000073e-01 6.572852999999999595e-04 6.599464729458228178e-04 6.977480000000000151e-05 +2.143460000000000087e-01 6.521788999999999841e-04 6.548279369688987054e-04 7.041510000000000282e-05 +2.147350000000000092e-01 6.471002000000000221e-04 6.497376769277467701e-04 7.062169999999999505e-05 +2.151250000000000107e-01 6.420698000000000507e-04 6.446946792448021322e-04 7.086940000000000286e-05 +2.155150000000000121e-01 6.371779000000000375e-04 6.397881268960462636e-04 7.109459999999999521e-05 +2.159040000000000126e-01 6.322927999999999578e-04 6.348883222406692716e-04 7.168090000000000009e-05 +2.162939999999999863e-01 6.275971999999999888e-04 6.301800963007824939e-04 7.212059999999999521e-05 +2.166839999999999877e-01 6.229798000000000150e-04 6.255508677629553500e-04 7.188699999999999799e-05 +2.170739999999999892e-01 6.183621999999999730e-04 6.209211039055799681e-04 7.209579999999999505e-05 +2.174629999999999896e-01 6.137747000000000038e-04 6.163210494368342616e-04 7.233290000000000672e-05 +2.178529999999999911e-01 6.091799999999999646e-04 6.117138551753275276e-04 7.376390000000000030e-05 +2.182429999999999926e-01 6.047643000000000333e-04 6.072841543603626352e-04 7.447619999999999685e-05 +2.186319999999999930e-01 6.004547999999999960e-04 6.029600849150960461e-04 7.445059999999999493e-05 +2.190219999999999945e-01 5.961450999999999990e-04 5.986359846397048710e-04 7.423320000000000612e-05 +2.194119999999999959e-01 5.919172999999999579e-04 5.944001196198046497e-04 7.694409999999999401e-05 +2.198009999999999964e-01 5.876931999999999859e-04 5.901680213972076462e-04 7.600549999999999592e-05 +2.201909999999999978e-01 5.835492000000000066e-04 5.860143891473464442e-04 7.501700000000000360e-05 +2.205809999999999993e-01 5.794997000000000516e-04 5.819536888610457506e-04 7.525920000000000613e-05 +2.209699999999999998e-01 5.754609999999999847e-04 5.779037453747489558e-04 7.535879999999999413e-05 +2.213600000000000012e-01 5.716140000000000197e-04 5.740447495814852150e-04 7.507610000000000444e-05 +2.217500000000000027e-01 5.677743000000000503e-04 5.701929119472412110e-04 7.616279999999999524e-05 +2.221390000000000031e-01 5.639681000000000121e-04 5.663751303869838001e-04 7.717400000000000345e-05 +2.225290000000000046e-01 5.601752999999999897e-04 5.625715545116477447e-04 7.489909999999999580e-05 +2.229190000000000060e-01 5.563828999999999953e-04 5.587683989680501803e-04 7.528420000000000673e-05 +2.233080000000000065e-01 5.527140999999999738e-04 5.550917580136095243e-04 7.629890000000000229e-05 +2.236980000000000079e-01 5.490538999999999576e-04 5.514245115115100262e-04 7.543390000000000295e-05 +2.240880000000000094e-01 5.454408999999999907e-04 5.478033392683911854e-04 7.638869999999999591e-05 +2.244770000000000099e-01 5.419897000000000126e-04 5.443404896498550954e-04 7.586849999999999367e-05 +2.248670000000000113e-01 5.385267000000000222e-04 5.408656732480038613e-04 7.697740000000000604e-05 +2.252570000000000128e-01 5.351541999999999931e-04 5.374823751808695954e-04 7.633019999999999638e-05 +2.256460000000000132e-01 5.318344000000000456e-04 5.341521670767570293e-04 7.570559999999999562e-05 +2.260359999999999869e-01 5.285119000000000448e-04 5.308194525007712021e-04 7.635630000000000618e-05 +2.264259999999999884e-01 5.252755000000000228e-04 5.275739897428699489e-04 7.607249999999999375e-05 +2.268149999999999888e-01 5.220592999999999501e-04 5.243488599752575653e-04 7.625800000000000070e-05 +2.272049999999999903e-01 5.188963000000000212e-04 5.211747521796385695e-04 7.635129999999999521e-05 +2.275949999999999918e-01 5.157939999999999747e-04 5.180595612221928180e-04 7.641360000000000307e-05 +2.279839999999999922e-01 5.127041999999999624e-04 5.149569273508962765e-04 7.657019999999999408e-05 +2.283739999999999937e-01 5.096855000000000179e-04 5.119318132104029627e-04 7.675929999999999537e-05 +2.287629999999999941e-01 5.066858000000000473e-04 5.089259154116908095e-04 7.713689999999999664e-05 +2.291529999999999956e-01 5.037321000000000420e-04 5.059626831677892581e-04 7.749009999999999861e-05 +2.295429999999999970e-01 5.008699000000000137e-04 5.030857209643141946e-04 8.049690000000000512e-05 +2.299319999999999975e-01 4.980200999999999855e-04 5.002211458746868081e-04 7.742350000000000165e-05 +2.303219999999999990e-01 4.952927999999999888e-04 4.974829497830962925e-04 7.756269999999999517e-05 +2.307120000000000004e-01 4.925853999999999476e-04 4.947654876549569805e-04 7.768140000000000472e-05 +2.311010000000000009e-01 4.898656000000000147e-04 4.920368189838526893e-04 7.824030000000000372e-05 +2.314910000000000023e-01 4.870828000000000115e-04 4.892485885268587189e-04 7.785819999999999937e-05 +2.318810000000000038e-01 4.842986999999999986e-04 4.864590357192099916e-04 7.792439999999999545e-05 +2.322700000000000042e-01 4.817013000000000089e-04 4.838507201814716660e-04 7.842570000000000368e-05 +2.326600000000000057e-01 4.791815000000000267e-04 4.813176280713186267e-04 7.779229999999999717e-05 +2.330490000000000061e-01 4.766713999999999912e-04 4.787944541501742242e-04 7.904830000000000005e-05 +2.334390000000000076e-01 4.742080999999999770e-04 4.763197785544048151e-04 7.839200000000000432e-05 +2.338290000000000091e-01 4.717464000000000206e-04 4.738467332739635112e-04 7.914910000000000424e-05 +2.342180000000000095e-01 4.693370999999999893e-04 4.714269911357107717e-04 7.886809999999999795e-05 +2.346080000000000110e-01 4.669499000000000132e-04 4.690299703306956064e-04 7.852709999999999562e-05 +2.349970000000000114e-01 4.645712999999999883e-04 4.666416556780947247e-04 7.977469999999999364e-05 +2.353870000000000129e-01 4.622499000000000074e-04 4.643079802055170261e-04 7.883379999999999728e-05 +2.357769999999999866e-01 4.599267999999999889e-04 4.619726744175325093e-04 8.025940000000000613e-05 +2.361659999999999870e-01 4.576649000000000233e-04 4.596989765219442283e-04 7.984699999999999632e-05 +2.365559999999999885e-01 4.554652000000000180e-04 4.574881665785909628e-04 7.910940000000000528e-05 +2.369449999999999890e-01 4.532679000000000178e-04 4.552797472259941902e-04 7.998979999999999773e-05 +2.373349999999999904e-01 4.511566999999999964e-04 4.531568947683767648e-04 7.928319999999999335e-05 +2.377249999999999919e-01 4.490598000000000267e-04 4.510483630473043049e-04 7.963420000000000405e-05 +2.381139999999999923e-01 4.469772999999999802e-04 4.489540424166306985e-04 8.045939999999999744e-05 +2.385039999999999938e-01 4.449292999999999890e-04 4.468936883618616169e-04 8.020559999999999658e-05 +2.388929999999999942e-01 4.428782000000000248e-04 4.448303652101485354e-04 8.114180000000000296e-05 +2.392829999999999957e-01 4.408912000000000182e-04 4.428314666446142556e-04 7.988270000000000006e-05 +2.396719999999999962e-01 4.389266000000000066e-04 4.408551351836859479e-04 8.049590000000000293e-05 +2.400619999999999976e-01 4.369667000000000256e-04 4.388833855576689544e-04 8.121580000000000259e-05 +2.404519999999999991e-01 4.350612000000000012e-04 4.369652799777722367e-04 8.075249999999999637e-05 +2.408409999999999995e-01 4.331595000000000258e-04 4.350510253312271799e-04 8.170340000000000111e-05 +2.412310000000000010e-01 4.313143999999999984e-04 4.331928183481683595e-04 8.053239999999999487e-05 +2.416200000000000014e-01 4.295174000000000019e-04 4.313821914650272206e-04 8.138619999999999676e-05 +2.420100000000000029e-01 4.277177000000000062e-04 4.295688660260522126e-04 8.119829999999999810e-05 +2.423990000000000034e-01 4.259333000000000238e-04 4.277723940753346845e-04 8.108620000000000302e-05 +2.427890000000000048e-01 4.241458999999999941e-04 4.259728032374355427e-04 8.161339999999999350e-05 +2.431790000000000063e-01 4.224106000000000041e-04 4.242241106834040748e-04 8.191940000000000040e-05 +2.435680000000000067e-01 4.207327000000000180e-04 4.225314046490879160e-04 8.168959999999999796e-05 +2.439580000000000082e-01 4.190514000000000109e-04 4.208352680127957048e-04 8.168959999999999796e-05 +2.443470000000000086e-01 4.173979999999999769e-04 4.191706055024867596e-04 8.222880000000000120e-05 +2.447370000000000101e-01 4.157496000000000217e-04 4.175113447655430131e-04 8.191940000000000040e-05 +2.451260000000000105e-01 4.141207999999999738e-04 4.158710927291010937e-04 8.226769999999999840e-05 +2.455160000000000120e-01 4.125228000000000206e-04 4.142605567231768984e-04 8.273919999999999550e-05 +2.459050000000000125e-01 4.109351000000000019e-04 4.126601872126887342e-04 8.188090000000000407e-05 +2.462949999999999862e-01 4.093753999999999904e-04 4.110863739720261581e-04 8.285830000000000593e-05 +2.466839999999999866e-01 4.078436999999999861e-04 4.095400122767430570e-04 8.365229999999999866e-05 +2.470739999999999881e-01 4.063197999999999853e-04 4.080012848605372978e-04 8.262060000000000649e-05 +2.474629999999999885e-01 4.048394999999999988e-04 4.065053693058205556e-04 8.242399999999999553e-05 +2.478529999999999900e-01 4.033589999999999984e-04 4.050091457871405633e-04 8.258109999999999442e-05 +2.482429999999999914e-01 4.019059000000000173e-04 4.035418672309337576e-04 8.293799999999999773e-05 +2.486319999999999919e-01 4.004659000000000040e-04 4.020884662091811132e-04 8.199640000000000661e-05 +2.490219999999999934e-01 3.990258999999999907e-04 4.006350075816027603e-04 8.281849999999999997e-05 +2.494109999999999938e-01 3.976412000000000241e-04 3.992355094572715710e-04 8.285799999999999849e-05 +2.498009999999999953e-01 3.962619000000000015e-04 3.978413710208203100e-04 8.437110000000000269e-05 +2.501900000000000235e-01 3.949418000000000004e-04 3.965056445224996741e-04 8.325919999999999729e-05 +2.505800000000000249e-01 3.936823999999999901e-04 3.952296165820911363e-04 8.188090000000000407e-05 +2.509689999999999976e-01 3.924250999999999911e-04 3.939556968932646447e-04 8.313829999999999647e-05 +2.513589999999999991e-01 3.910785000000000064e-04 3.925955111759000886e-04 8.250240000000000481e-05 +2.517480000000000273e-01 3.897311000000000200e-04 3.912347717919792997e-04 8.285830000000000593e-05 +2.521379999999999733e-01 3.884205999999999857e-04 3.899102961176904121e-04 8.301799999999999696e-05 +2.525270000000000015e-01 3.871979000000000219e-04 3.886728415363535226e-04 8.230670000000000260e-05 +2.529170000000000029e-01 3.859741000000000083e-04 3.874341769715277956e-04 8.350249999999999546e-05 +2.533059999999999756e-01 3.847324000000000165e-04 3.861768757179676612e-04 8.207370000000000670e-05 +2.536959999999999771e-01 3.834845000000000246e-04 3.849129516077297937e-04 8.289809999999999833e-05 +2.540850000000000053e-01 3.822664999999999830e-04 3.836790206950234914e-04 8.329960000000000456e-05 +2.544750000000000068e-01 3.811107000000000101e-04 3.825071998295300683e-04 8.269960000000000353e-05 +2.548639999999999795e-01 3.799572999999999883e-04 3.813377695555121810e-04 8.317859999999999674e-05 +2.552539999999999809e-01 3.788126000000000058e-04 3.801785214937243476e-04 8.258109999999999442e-05 +2.556430000000000091e-01 3.776833999999999964e-04 3.790356802368470535e-04 8.254179999999999634e-05 +2.560330000000000106e-01 3.765493999999999764e-04 3.778877002759607902e-04 8.269960000000000353e-05 +2.564219999999999833e-01 3.754230000000000002e-04 3.767455282664382418e-04 8.273919999999999550e-05 +2.568110000000000115e-01 3.742891000000000143e-04 3.755957183447737992e-04 8.285830000000000593e-05 +2.572010000000000129e-01 3.732103999999999869e-04 3.745010042164901242e-04 8.184259999999999463e-05 +2.575899999999999856e-01 3.721775999999999991e-04 3.734519514324489520e-04 8.281849999999999997e-05 +2.579799999999999871e-01 3.711429000000000139e-04 3.724009859942455591e-04 8.207370000000000670e-05 +2.583690000000000153e-01 3.700840999999999961e-04 3.713283543449538574e-04 8.254179999999999634e-05 +2.587590000000000168e-01 3.690182999999999766e-04 3.702488488100828725e-04 8.238480000000000445e-05 +2.591479999999999895e-01 3.679803000000000045e-04 3.691954650533356896e-04 8.218990000000000399e-05 +2.595379999999999909e-01 3.669831000000000135e-04 3.681800772745679991e-04 8.258109999999999442e-05 +2.599270000000000191e-01 3.659876000000000058e-04 3.671664897833429937e-04 8.203499999999999638e-05 +2.603170000000000206e-01 3.650066999999999895e-04 3.661697568317291880e-04 8.262060000000000649e-05 +2.607059999999999933e-01 3.640278000000000047e-04 3.651756072583566904e-04 8.218990000000000399e-05 +2.610959999999999948e-01 3.630583000000000269e-04 3.641912245035230683e-04 8.250240000000000481e-05 +2.614850000000000230e-01 3.620726000000000001e-04 3.631919706771840435e-04 8.269960000000000353e-05 +2.618739999999999957e-01 3.610870000000000074e-04 3.621928868228600513e-04 8.191940000000000040e-05 +2.622639999999999971e-01 3.601376999999999790e-04 3.612282867518446048e-04 8.289809999999999833e-05 +2.626530000000000253e-01 3.592089999999999823e-04 3.602834142411698824e-04 8.262060000000000649e-05 +2.630430000000000268e-01 3.582801000000000258e-04 3.593384622506609334e-04 8.234569999999999325e-05 +2.634319999999999995e-01 3.573768999999999969e-04 3.584197742395971425e-04 8.285830000000000593e-05 +2.638220000000000010e-01 3.564659000000000187e-04 3.574932783443720880e-04 8.211240000000000346e-05 +2.642109999999999737e-01 3.555758999999999916e-04 3.565886486139654022e-04 8.297800000000000412e-05 +2.646000000000000019e-01 3.546822999999999838e-04 3.556811108155387099e-04 8.234569999999999325e-05 +2.649900000000000033e-01 3.537950999999999900e-04 3.547800009455959929e-04 8.277880000000000101e-05 +2.653789999999999760e-01 3.529100999999999874e-04 3.538787492590604643e-04 8.289860000000000620e-05 +2.657689999999999775e-01 3.520199999999999805e-04 3.529722182722326530e-04 8.273919999999999550e-05 +2.661580000000000057e-01 3.511594999999999842e-04 3.520969858921923590e-04 8.289809999999999833e-05 +2.665480000000000071e-01 3.503359999999999743e-04 3.512608985775068850e-04 8.234569999999999325e-05 +2.669369999999999798e-01 3.495048999999999747e-04 3.504172053304884546e-04 8.277880000000000101e-05 +2.673260000000000081e-01 3.486849000000000198e-04 3.495811798816702274e-04 8.258109999999999442e-05 +2.677160000000000095e-01 3.478703999999999889e-04 3.487500496461259621e-04 8.305800000000000335e-05 +2.681049999999999822e-01 3.470608999999999825e-04 3.479246821251523675e-04 8.297800000000000412e-05 +2.684949999999999837e-01 3.462585000000000120e-04 3.471087628822637060e-04 8.242399999999999553e-05 +2.688840000000000119e-01 3.454654000000000145e-04 3.463021118659109494e-04 8.265999999999999802e-05 +2.692729999999999846e-01 3.446712000000000213e-04 3.454942881872301113e-04 8.246320000000000017e-05 +2.696629999999999860e-01 3.438798000000000072e-04 3.446891478714741454e-04 8.285830000000000593e-05 +2.700520000000000143e-01 3.430975999999999861e-04 3.438934777345286492e-04 8.254179999999999634e-05 +2.704420000000000157e-01 3.422990999999999903e-04 3.430824890824102922e-04 8.218990000000000399e-05 +2.708309999999999884e-01 3.415086000000000120e-04 3.422795586736045235e-04 8.285830000000000593e-05 +2.712200000000000166e-01 3.407489999999999998e-04 3.415058869621151135e-04 8.184259999999999463e-05 +2.716100000000000181e-01 3.400189000000000184e-04 3.407605521124414962e-04 8.254179999999999634e-05 +2.719989999999999908e-01 3.392892000000000108e-04 3.400154996013372865e-04 8.180420000000000530e-05 +2.723880000000000190e-01 3.385379999999999900e-04 3.392519909690060675e-04 8.215110000000000023e-05 +2.727780000000000205e-01 3.377866999999999892e-04 3.384883123643084800e-04 8.246320000000000017e-05 +2.731669999999999932e-01 3.370491999999999781e-04 3.377384225828707932e-04 8.157540000000000505e-05 +2.735569999999999946e-01 3.363228000000000116e-04 3.369993744033151443e-04 8.222880000000000120e-05 +2.739460000000000228e-01 3.355961999999999769e-04 3.362602138563092983e-04 8.161339999999999350e-05 +2.743349999999999955e-01 3.348962999999999940e-04 3.355484832650906172e-04 8.222880000000000120e-05 +2.747249999999999970e-01 3.342042000000000146e-04 3.348447535514761497e-04 8.195789999999999673e-05 +2.751140000000000252e-01 3.335179999999999872e-04 3.341463513027778383e-04 8.165150000000000251e-05 +2.755029999999999979e-01 3.328184999999999781e-04 3.334330144753725107e-04 8.218990000000000399e-05 +2.758929999999999993e-01 3.321258999999999908e-04 3.327265552837520160e-04 8.184259999999999463e-05 +2.762820000000000276e-01 3.314357000000000089e-04 3.320250677655894249e-04 8.203499999999999638e-05 +2.766710000000000003e-01 3.307432000000000015e-04 3.313223192860144536e-04 8.203499999999999638e-05 +2.770610000000000017e-01 3.300472999999999732e-04 3.306159145585010656e-04 8.207370000000000670e-05 +2.774499999999999744e-01 3.293694000000000114e-04 3.299267388617606453e-04 8.246320000000000017e-05 +2.778390000000000026e-01 3.286914999999999954e-04 3.292375631648275080e-04 8.168959999999999796e-05 +2.782290000000000041e-01 3.280377000000000119e-04 3.285727736867610988e-04 8.234569999999999325e-05 +2.786179999999999768e-01 3.274161000000000042e-04 3.279403822881075091e-04 8.165150000000000251e-05 +2.790079999999999782e-01 3.267888999999999842e-04 3.273023232365097337e-04 8.211240000000000346e-05 +2.793970000000000065e-01 3.261428000000000244e-04 3.266454499426300570e-04 8.289809999999999833e-05 +2.797859999999999792e-01 3.254962000000000025e-04 3.259879863446824922e-04 8.153749999999999648e-05 +2.801750000000000074e-01 3.248607000000000252e-04 3.253424293529359659e-04 8.246320000000000017e-05 +2.805650000000000088e-01 3.242284999999999807e-04 3.247012024111746746e-04 8.188090000000000407e-05 +2.809539999999999815e-01 3.236021999999999965e-04 3.240658130936892832e-04 8.270060000000000573e-05 +2.813430000000000097e-01 3.229795999999999731e-04 3.234336121555275809e-04 8.218990000000000399e-05 +2.817330000000000112e-01 3.223511999999999775e-04 3.227955159873011331e-04 8.222880000000000120e-05 +2.821219999999999839e-01 3.217324000000000030e-04 3.221673814313388731e-04 8.250240000000000481e-05 +2.825110000000000121e-01 3.211372999999999788e-04 3.215636167837490991e-04 8.222880000000000120e-05 +2.829010000000000136e-01 3.205344999999999852e-04 3.209522142245878839e-04 8.226769999999999840e-05 +2.832899999999999863e-01 3.199211999999999889e-04 3.203307399967297636e-04 8.199640000000000661e-05 +2.836790000000000145e-01 3.193045000000000258e-04 3.197059901843678654e-04 8.234569999999999325e-05 +2.840690000000000159e-01 3.186921999999999910e-04 3.190856144799374583e-04 8.269960000000000353e-05 +2.844579999999999886e-01 3.181235000000000247e-04 3.185084241350567814e-04 8.188090000000000407e-05 +2.848470000000000169e-01 3.175531000000000209e-04 3.179296034750427220e-04 8.281849999999999997e-05 +2.852359999999999896e-01 3.169784999999999942e-04 3.173476001360786862e-04 8.180420000000000530e-05 +2.856259999999999910e-01 3.163948999999999885e-04 3.167573903861529592e-04 8.246320000000000017e-05 +2.860150000000000192e-01 3.158189000000000265e-04 3.161747609431665089e-04 8.265999999999999802e-05 +2.864039999999999919e-01 3.152512000000000218e-04 3.156000446020477230e-04 8.238480000000000445e-05 +2.867939999999999934e-01 3.146771999999999829e-04 3.150190383251282684e-04 8.242399999999999553e-05 +2.871830000000000216e-01 3.141202999999999747e-04 3.144554989098492489e-04 8.218990000000000399e-05 +2.875719999999999943e-01 3.135693000000000269e-04 3.138983131810690322e-04 8.281849999999999997e-05 +2.879610000000000225e-01 3.130173000000000092e-04 3.133400874412470651e-04 8.188090000000000407e-05 +2.883510000000000240e-01 3.124709999999999836e-04 3.127877833527109575e-04 8.254179999999999634e-05 +2.887399999999999967e-01 3.119258999999999877e-04 3.122365987546242474e-04 8.281849999999999997e-05 +2.891290000000000249e-01 3.113915999999999885e-04 3.116958610556516554e-04 8.234569999999999325e-05 +2.895190000000000263e-01 3.108664999999999822e-04 3.111639123794826631e-04 8.305800000000000335e-05 +2.899079999999999990e-01 3.103345999999999882e-04 3.106250860680211741e-04 8.289809999999999833e-05 +2.902970000000000272e-01 3.098102000000000039e-04 3.100968226842441269e-04 8.329960000000000456e-05 +2.906859999999999999e-01 3.092856000000000056e-04 3.095692622435483187e-04 8.407869999999999851e-05 +2.910760000000000014e-01 3.087618999999999889e-04 3.090422500213310989e-04 8.313829999999999647e-05 +2.914649999999999741e-01 3.082444000000000264e-04 3.085199056372787030e-04 8.496570000000000544e-05 +2.918540000000000023e-01 3.077266999999999958e-04 3.079974232607622263e-04 8.358409999999999820e-05 +2.922429999999999750e-01 3.072207999999999775e-04 3.074874499239215841e-04 8.470940000000000587e-05 +2.926329999999999765e-01 3.067187000000000082e-04 3.069816375423161476e-04 8.445540000000000458e-05 +2.930220000000000047e-01 3.062264000000000197e-04 3.064856836911707930e-04 8.479460000000000296e-05 +2.934109999999999774e-01 3.057250000000000180e-04 3.059810751879292164e-04 8.548539999999999980e-05 +2.938000000000000056e-01 3.052234999999999822e-04 3.054763286922513687e-04 8.479460000000000296e-05 +2.941889999999999783e-01 3.047333999999999850e-04 3.049831701519195625e-04 8.592580000000000324e-05 +2.945789999999999798e-01 3.042586000000000010e-04 3.045056688852436333e-04 8.535459999999999760e-05 +2.949680000000000080e-01 3.037785999999999785e-04 3.040229202980747640e-04 8.592580000000000324e-05 +2.953569999999999807e-01 3.032835999999999909e-04 3.035251375254263413e-04 8.610389999999999396e-05 +2.957460000000000089e-01 3.027831000000000251e-04 3.030218250930092130e-04 8.916239999999999864e-05 +2.961360000000000103e-01 3.023052999999999938e-04 3.025416491063127834e-04 9.207629999999999797e-05 +2.965249999999999830e-01 3.018715000000000048e-04 3.021063378256179681e-04 9.246090000000000103e-05 +2.969140000000000112e-01 3.014297999999999782e-04 3.016631382737857336e-04 9.416600000000000610e-05 +2.973029999999999839e-01 3.009938999999999779e-04 3.012255466088361732e-04 9.615489999999999649e-05 +2.976920000000000122e-01 3.005495999999999862e-04 3.007794956393054645e-04 9.703769999999999421e-05 +2.980820000000000136e-01 3.001017999999999936e-04 3.003300250798394836e-04 9.915039999999999856e-05 +2.984709999999999863e-01 2.996515000000000159e-04 2.998785606043726068e-04 1.008286000000000056e-04 +2.988600000000000145e-01 2.992000999999999884e-04 2.994261978608201650e-04 1.023691999999999956e-04 +2.992489999999999872e-01 2.987414000000000193e-04 2.989671735280575635e-04 1.031280000000000060e-04 +2.996380000000000154e-01 2.982851000000000014e-04 2.985110828290101458e-04 1.052644000000000013e-04 +3.000280000000000169e-01 2.978276000000000079e-04 2.980536954049686359e-04 1.071930000000000018e-04 +3.004169999999999896e-01 2.974362999999999833e-04 2.976622144129323659e-04 1.076264000000000035e-04 +3.008060000000000178e-01 2.970434999999999894e-04 2.972692730780028265e-04 1.099684000000000024e-04 +3.011949999999999905e-01 2.966461999999999788e-04 2.968716871976620399e-04 1.111651999999999974e-04 +3.015840000000000187e-01 2.962471999999999848e-04 2.964723252996577518e-04 1.130720999999999978e-04 +3.019729999999999914e-01 2.958465000000000075e-04 2.960711631146576502e-04 1.151507999999999953e-04 +3.023629999999999929e-01 2.954343000000000117e-04 2.956596889000617074e-04 1.155804999999999956e-04 +3.027520000000000211e-01 2.950240000000000132e-04 2.952502947077579741e-04 1.179139999999999962e-04 +3.031409999999999938e-01 2.946258000000000210e-04 2.948528261236499728e-04 1.201248999999999991e-04 +3.035300000000000220e-01 2.942433000000000158e-04 2.944709295258800047e-04 1.223413000000000072e-04 +3.039189999999999947e-01 2.938530999999999869e-04 2.940812570235265587e-04 1.231830999999999963e-04 +3.043080000000000229e-01 2.934821999999999800e-04 2.937116032257194777e-04 1.248210000000000031e-04 +3.046980000000000244e-01 2.931133000000000046e-04 2.933442057771573738e-04 1.263055999999999923e-04 +3.050869999999999971e-01 2.927405000000000003e-04 2.929730133024328062e-04 1.301193000000000130e-04 +3.054760000000000253e-01 2.923455999999999950e-04 2.925798728506316757e-04 1.311964000000000099e-04 +3.058649999999999980e-01 2.919530999999999949e-04 2.921891229900896742e-04 1.326278999999999978e-04 +3.062540000000000262e-01 2.915967999999999794e-04 2.918346506990721371e-04 1.327878999999999963e-04 +3.066429999999999989e-01 2.912583000000000164e-04 2.914980668315784724e-04 1.359770999999999961e-04 +3.070330000000000004e-01 2.909181999999999957e-04 2.911598013983691426e-04 1.359777000000000110e-04 +3.074219999999999731e-01 2.905863999999999865e-04 2.908298558052878455e-04 1.388615000000000053e-04 +3.078110000000000013e-01 2.902456999999999780e-04 2.904909819297539974e-04 1.386927999999999946e-04 +3.081999999999999740e-01 2.899120000000000256e-04 2.901592602158933323e-04 1.372109000000000046e-04 +3.085890000000000022e-01 2.895684999999999839e-04 2.898179864775426041e-04 1.372097000000000019e-04 +3.089779999999999749e-01 2.892239000000000008e-04 2.894755603607760937e-04 1.373872000000000049e-04 +3.093670000000000031e-01 2.889093999999999820e-04 2.891639687079120588e-04 1.404807000000000050e-04 +3.097559999999999758e-01 2.885952999999999912e-04 2.888528549917964154e-04 1.388170999999999893e-04 +3.101459999999999773e-01 2.882841999999999935e-04 2.885446645662028027e-04 1.393598000000000070e-04 +3.105350000000000055e-01 2.879916000000000160e-04 2.882546986399607655e-04 1.387952000000000023e-04 +3.109239999999999782e-01 2.876926000000000244e-04 2.879584747578470185e-04 1.381104999999999985e-04 +3.113130000000000064e-01 2.873911999999999734e-04 2.876596672921420661e-04 1.395496000000000030e-04 +3.117019999999999791e-01 2.870918000000000081e-04 2.873629424524224091e-04 1.421229999999999943e-04 +3.120910000000000073e-01 2.867985999999999886e-04 2.870725305695521383e-04 1.431030000000000019e-04 +3.124799999999999800e-01 2.865276999999999841e-04 2.868047998748426714e-04 1.397399999999999868e-04 +3.128690000000000082e-01 2.862494999999999840e-04 2.865297136077542078e-04 1.424918999999999969e-04 +3.132579999999999809e-01 2.859730000000000215e-04 2.862566201767633291e-04 1.447575000000000044e-04 +3.136470000000000091e-01 2.856897999999999968e-04 2.859768865536277253e-04 1.437219999999999904e-04 +3.140370000000000106e-01 2.854147000000000238e-04 2.857053785418446016e-04 1.503511999999999980e-04 +3.144259999999999833e-01 2.851482000000000018e-04 2.854426330389806194e-04 1.491061999999999921e-04 +3.148150000000000115e-01 2.848902999999999852e-04 2.851884237164071427e-04 1.485733999999999894e-04 +3.152039999999999842e-01 2.846341000000000062e-04 2.849362166330986776e-04 1.487511000000000063e-04 +3.155930000000000124e-01 2.843874000000000141e-04 2.846936177211068300e-04 1.462811000000000113e-04 +3.159819999999999851e-01 2.841335999999999862e-04 2.844440031813924488e-04 1.491451000000000029e-04 +3.163710000000000133e-01 2.839121000000000223e-04 2.842258852343302125e-04 1.484608000000000000e-04 +3.167599999999999860e-01 2.836825000000000067e-04 2.839995390720561480e-04 1.529481000000000070e-04 +3.171490000000000142e-01 2.834544999999999947e-04 2.837750277808021807e-04 1.570753000000000036e-04 +3.175379999999999869e-01 2.832267999999999765e-04 2.835510327617293607e-04 1.478384999999999975e-04 +3.179270000000000151e-01 2.830067000000000021e-04 2.833346756543718556e-04 1.539524000000000068e-04 +3.183159999999999878e-01 2.828026999999999885e-04 2.831343504477663907e-04 1.523301999999999870e-04 +3.187050000000000161e-01 2.825972000000000054e-04 2.829325136488639015e-04 1.544651999999999928e-04 +3.190939999999999888e-01 2.823982999999999962e-04 2.827374075339390228e-04 1.538359999999999955e-04 +3.194830000000000170e-01 2.821968000000000219e-04 2.825397433566017506e-04 1.548649999999999886e-04 +3.198719999999999897e-01 2.819898000000000153e-04 2.823365495193574828e-04 1.572579000000000110e-04 +3.202619999999999911e-01 2.817887000000000148e-04 2.821397133289695589e-04 1.556504999999999965e-04 +3.206510000000000193e-01 2.815925999999999846e-04 2.819480184465777799e-04 1.583218000000000060e-04 +3.210399999999999920e-01 2.813914000000000042e-04 2.817512142358586606e-04 1.589162000000000083e-04 +3.214290000000000203e-01 2.812264999999999882e-04 2.815899062007329464e-04 1.612678999999999946e-04 +3.218179999999999930e-01 2.810627000000000221e-04 2.814298081491948948e-04 1.659742999999999874e-04 +3.222070000000000212e-01 2.808965999999999763e-04 2.812677553906732736e-04 1.615507999999999983e-04 +3.225959999999999939e-01 2.807290999999999953e-04 2.811046755693940571e-04 1.621755000000000061e-04 +3.229850000000000221e-01 2.805689999999999898e-04 2.809490636878003776e-04 1.639381000000000085e-04 +3.233739999999999948e-01 2.804433999999999855e-04 2.808265596311325298e-04 1.633284999999999999e-04 +3.237630000000000230e-01 2.803101999999999916e-04 2.806964176623623613e-04 1.670692000000000097e-04 +3.241519999999999957e-01 2.801767999999999837e-04 2.805663875138520043e-04 1.623965999999999962e-04 +3.245410000000000239e-01 2.800222999999999905e-04 2.804158037509697017e-04 1.691728999999999943e-04 +3.249299999999999966e-01 2.798744000000000253e-04 2.802717859079632877e-04 1.690384999999999978e-04 +3.253190000000000248e-01 2.797472999999999974e-04 2.801487714205305648e-04 1.698040999999999960e-04 +3.257079999999999975e-01 2.796283000000000210e-04 2.800338434058708662e-04 1.683768000000000038e-04 +3.260970000000000257e-01 2.795105000000000202e-04 2.799199516517741492e-04 1.717026999999999985e-04 +3.264859999999999984e-01 2.794000999999999950e-04 2.798130599727718243e-04 1.714243000000000114e-04 +3.268750000000000266e-01 2.792884000000000143e-04 2.797048459433649223e-04 1.809872999999999875e-04 +3.272639999999999993e-01 2.791758999999999777e-04 2.795964360080406097e-04 1.710141999999999999e-04 +3.276530000000000276e-01 2.790673000000000242e-04 2.794922665078223754e-04 1.745460000000000056e-04 +3.280420000000000003e-01 2.789588999999999762e-04 2.793880676323133670e-04 1.790577000000000119e-04 +3.284309999999999730e-01 2.788690999999999826e-04 2.793018478638439406e-04 1.755163999999999921e-04 +3.288200000000000012e-01 2.787795000000000030e-04 2.792157980678862199e-04 1.770495000000000131e-04 +3.292089999999999739e-01 2.787062999999999780e-04 2.791455373694078264e-04 1.819953000000000022e-04 +3.295980000000000021e-01 2.786501999999999838e-04 2.790919715385649586e-04 1.798648000000000130e-04 +3.299869999999999748e-01 2.785921000000000122e-04 2.790362974550876095e-04 1.806505999999999878e-04 +3.303760000000000030e-01 2.785096999999999942e-04 2.789575760155024189e-04 1.929718000000000080e-04 +3.307649999999999757e-01 2.784263999999999945e-04 2.788778987026790815e-04 1.865740000000000063e-04 +3.311529999999999752e-01 2.783627999999999905e-04 2.788176080128882978e-04 1.816882000000000132e-04 +3.315420000000000034e-01 2.783212999999999876e-04 2.787786984474626599e-04 1.850993000000000050e-04 +3.319309999999999761e-01 2.782808000000000005e-04 2.787408288929876994e-04 1.885498999999999953e-04 +3.323200000000000043e-01 2.782226999999999747e-04 2.786856020376570105e-04 1.943861000000000124e-04 +3.327089999999999770e-01 2.781713999999999910e-04 2.786373945610160570e-04 1.922694000000000128e-04 +3.330980000000000052e-01 2.781282999999999846e-04 2.785972991827907405e-04 1.936216000000000098e-04 +3.334869999999999779e-01 2.781045000000000001e-04 2.785762619442790902e-04 1.895111999999999958e-04 +3.338760000000000061e-01 2.780904999999999965e-04 2.785650512567331653e-04 1.995795000000000023e-04 +3.342649999999999788e-01 2.780640999999999929e-04 2.785415119668670550e-04 1.914451000000000011e-04 +3.346540000000000070e-01 2.780376000000000093e-04 2.785177770799753894e-04 2.080398999999999866e-04 +3.350429999999999797e-01 2.780125000000000153e-04 2.784953901687556711e-04 1.967505999999999997e-04 +3.354320000000000079e-01 2.779993999999999933e-04 2.784851218586190352e-04 2.038550999999999992e-04 +3.358209999999999806e-01 2.779861000000000116e-04 2.784745455842722222e-04 2.015398999999999912e-04 +3.362100000000000088e-01 2.779875000000000011e-04 2.784781080324947060e-04 2.103934999999999974e-04 +3.365989999999999815e-01 2.779834999999999923e-04 2.784757101450725000e-04 2.026726000000000016e-04 +3.369869999999999810e-01 2.779880000000000090e-04 2.784819582003916641e-04 2.038153000000000068e-04 +3.373760000000000092e-01 2.780122000000000214e-04 2.785084795041505349e-04 2.283696000000000098e-04 +3.377649999999999819e-01 2.780378000000000233e-04 2.785364931303875081e-04 2.369937999999999873e-04 +3.381540000000000101e-01 2.780599000000000243e-04 2.785606129473946254e-04 2.330691000000000078e-04 +3.385429999999999828e-01 2.780775000000000087e-04 2.785794535527421202e-04 2.513144999999999941e-04 +3.389320000000000110e-01 2.780867000000000018e-04 2.785899279491779322e-04 2.499806000000000033e-04 +3.393209999999999837e-01 2.781112999999999879e-04 2.786163928818992993e-04 2.681560000000000259e-04 +3.397100000000000120e-01 2.781419999999999943e-04 2.786490929886539488e-04 2.847969999999999907e-04 +3.400989999999999847e-01 2.781673000000000023e-04 2.786761153378254988e-04 2.799604000000000242e-04 +3.404869999999999841e-01 2.782223000000000009e-04 2.787325158340676944e-04 2.986880000000000242e-04 +3.408760000000000123e-01 2.782697000000000100e-04 2.787813103991318257e-04 3.000426000000000265e-04 +3.412649999999999850e-01 2.783198000000000182e-04 2.788325749176380356e-04 3.488532999999999825e-04 +3.416540000000000132e-01 2.783662999999999914e-04 2.788801017760343367e-04 3.978211000000000052e-04 +3.420429999999999859e-01 2.784075000000000005e-04 2.789222689472517546e-04 4.270850999999999879e-04 +3.424320000000000142e-01 2.784590999999999781e-04 2.789749536591799753e-04 5.026348000000000144e-04 +3.428209999999999869e-01 2.785028000000000265e-04 2.790197500992105343e-04 6.213787999999999603e-04 +3.432100000000000151e-01 2.785644999999999788e-04 2.790825270358666043e-04 7.163046999999999860e-04 +3.435980000000000145e-01 2.786366999999999880e-04 2.791557882985907439e-04 9.452773000000000035e-04 diff --git a/data/stan_scale.png b/data/stan_scale.png new file mode 100644 index 0000000..e8d4ac6 Binary files /dev/null and b/data/stan_scale.png differ diff --git a/data/stan_weights.png b/data/stan_weights.png index 04c5b5e..cd15adc 100644 Binary files a/data/stan_weights.png and b/data/stan_weights.png differ diff --git a/docs/METHODS_HDX_XL.md b/docs/METHODS_HDX_XL.md new file mode 100644 index 0000000..fe0eed8 --- /dev/null +++ b/docs/METHODS_HDX_XL.md @@ -0,0 +1,137 @@ +# Materials and Methods: HDX-MS and XL-MS extensions in bioce + +This section describes the hydrogen–deuterium exchange mass spectrometry (HDX-MS) and cross-linking mass spectrometry (XL-MS) likelihoods added to the bioce Bayesian ensemble framework. Small-angle X-ray scattering (SAXS) inference in bioce is described in Potrzebowski et al. (PLOS Computational Biology, 2018); here we assume the same conformational library and population-weight formalism, and document only the new modalities and their integration. + +**Workflow diagram:** [WORKFLOW.md](WORKFLOW.md) + +--- + +## Overview + +We extend bioce to combine HDX-MS and XL-MS data with a shared discrete ensemble of *K* structural conformers {*S*₁,…,*S*ₖ}. Each conformer has a non-negative population weight *w*ₖ with Σₖ *w*ₖ = 1. Experimental observables are predicted as ensemble averages over conformers, and independent likelihoods for each modality are multiplied in a joint Stan model (`fullBayesianMultimodal.py`, `stan_models.py`). When SAXS is included in the same run, it shares the same weight vector; the SAXS scale factor applies only to the scattering term and is not part of the HDX or XL blocks. + +--- + +## Structural library and per-conformer observables + +The input structural library is the same as in the original bioce workflow: one three-dimensional model per ensemble member (PDB format). For each modality we precompute a matrix of predicted observables with rows indexed by experimental observations (HDX peptides or time points; XL restraints) and columns indexed by conformers. + +### HDX-MS predicted uptake + +For each peptide *i*, defined by chain identifier, and inclusive residue range [*start*ᵢ, *end*ᵢ], we compute a structure-based solvent-exposure proxy from conformer *S*ₖ. Residue solvent-accessible surface area (SASA) is calculated with the Shrake–Rupley algorithm (BioPython `ShrakeRupley`, residue-level). The mean SASA over standard amino-acid residues in the peptide range is divided by 200 Ų and capped at 1: + +$$D_{\mathrm{pred},i}(S_k) = \min\!\left(1,\; \frac{1}{|P_i|}\sum_{r \in P_i} \frac{\mathrm{SASA}_r(S_k)}{200}\right)$$ + +where *P*ᵢ is the set of residues in peptide *i* and |*P*ᵢ| is the number of residues contributing. Higher values correspond to greater apparent solvent exposure (a simple proxy for deuterium uptake). This implementation does not fit exchange kinetics (e.g. logistic or Weibull curves); users may supply externally computed per-conformer uptake matrices via `prepareHDX.py` instead of the built-in SASA proxy (`structure_observables.py`). + +### XL-MS predicted distances + +For each crosslink restraint *l* between residues (*chain*ᵢ, *res*ᵢ) and (*chain*ⱼ, *res*ⱼ), we compute the Euclidean distance between preferred backbone atoms in conformer *S*ₖ, using Cα when present, otherwise Cβ or N. The distance *d*ₗ(*S*ₖ) is stored in Å. If atoms or residues are missing, a penalty distance of 999 Å is assigned so the restraint is effectively unsatisfied in that conformer. + +--- + +## HDX-MS likelihood + +Experimental HDX data are reduced to peptide-level deuterium uptake values *y*ᵢ with associated uncertainties *σ*ᵢ. Observations may represent a single time point per peptide (one row per peptide) or multiple time points encoded as flattened rows with peptide index and time (`multimodal_io.py`). + +The ensemble-averaged predicted uptake is + +$$\hat{D}_i = \sum_{k=1}^{K} w_k\, D_{\mathrm{pred},i}(S_k).$$ + +Each observation is modeled with a Gaussian likelihood: + +$$y_i \sim \mathcal{N}(\hat{D}_i,\, \sigma_i^2).$$ + +In Stan (`stan_code_HDX`, and jointly in `stan_code_SAXS_HDX`, `stan_code_HDX_XL`, `stan_code_SAXS_HDX_XL`), the predicted vector is formed as **sim_hdx** × **w**, where **sim_hdx** is the *n*ₒᵦₛ × *K* matrix of *D*ₚᵣₑ𝒹,ᵢ(*S*ₖ). + +--- + +## XL-MS likelihood + +Crosslinks are treated as probabilistic distance restraints rather than hard cutoffs. For restraint *l*, let *z*ₗ ∈ {0,1} indicate whether the crosslink was reported (*z*ₗ = 1) or listed as absent/decoy (*z*ₗ = 0). User-specified or default parameters are the effective maximum linker distance *d*ₘₐₓ,ₗ (default 30 Å) and softness *τ*ₗ (default 3 Å), which absorb linker length, flexibility, and assignment uncertainty. + +Per conformer, a soft satisfaction score is + +$$s_{l,k} = \mathrm{logit}^{-1}\!\left(\frac{d_{\max,l} - d_l(S_k)}{\tau_l}\right) = \frac{1}{1 + \exp\!\left(\frac{d_l(S_k) - d_{\max,l}}{\tau_l}\right)}.$$ + +The ensemble probability that the restraint is satisfied is + +$$p_{\mathrm{sat},l} = \sum_{k=1}^{K} w_k\, s_{l,k}.$$ + +Observed crosslinks are modeled with a Bernoulli likelihood: + +$$z_l \sim \mathrm{Bernoulli}(p_{\mathrm{sat},l}).$$ + +An optional false-positive mixture (`stan_code_XL_FP`, CLI flag `--xl-fp-mixture`) introduces a global false-positive rate *π*ₑₚ with prior *π*ₑₚ ~ Beta(1, 9): + +$$P(z_l = 1) = (1 - \pi_{\mathrm{fp}})\, p_{\mathrm{sat},l} + \pi_{\mathrm{fp}}\, (1 - p_{\mathrm{sat},l}).$$ + +**Non-observed crosslinks** (links not reported in the restraint list) are not scored as negatives, avoiding penalty for conformers when detectability, digestion, or sampling limits crosslink identification. + +Posterior satisfaction probabilities *p*ₛₐₜ,ₗ are computed in generated quantities (`xl_posterior_predict`) for each MCMC draw and summarized (e.g. posterior mean per restraint). + +--- + +## Priors and Bayesian inference + +Population weights are constrained to the simplex with a Dirichlet prior, + +$$\mathbf{w} \sim \mathrm{Dirichlet}(\boldsymbol{\alpha}),$$ + +where **α** is read from a user file (e.g. uniform αₖ = 1, or energy-informed priors when combined with the SAXS+energy models in the original bioce release). HDX and XL modules use the same **w** as SAXS when run jointly. + +Inference uses Hamiltonian Monte Carlo with the No-U-Turn Sampler (NUTS) via PyStan 3 (`stan_run.build_and_sample`). Total iterations per chain are split equally into warmup and sampling (PyStan 2–compatible default: half warmup, half samples). Typical CLI defaults are 2000 iterations and 4 chains (`fullBayesianMultimodal.py`); the Streamlit app uses lower defaults for interactive use. + +Software: Python 3, NumPy, Stan/PyStan 3 (httpstan backend), BioPython for PDB parsing and SASA when using the automated pipeline (`bioce_pipeline.py`, `app.py`). + +--- + +## Data preparation and software + +**Command-line workflow** + +| Step | Script | Output | +|------|--------|--------| +| HDX simulated matrix | `prepareHDX.py` | `SimulatedHDX.txt` (rows = observations, columns = conformers) | +| XL distance matrix | `prepareXLMS.py` | `SimulatedXLdistances.txt` | +| Inference | `fullBayesianMultimodal.py` | Posterior samples, weight plots | + +**Experimental file formats** + +- `hdx_exp.dat`: columns `uptake` `sigma`, or `peptide_index` `time` `uptake` `sigma` for time series. +- `hdx_peptides.txt` (pipeline/app): `chain` `start` `end` per line. +- `xl_restraints.dat`: `res_i` `res_j` `z` [`d_max`] [`tau`], or `chain_i` `res_i` `chain_j` `res_j` `z` [`d_max`] [`tau`]. + +**Web interface** + +The Streamlit application (`app.py`) runs the same Stan models: users upload PDB conformers and experimental files; the app builds simulated matrices (SASA-based HDX proxy, Cα distances for XL) and calls `bioce_pipeline.run_full_pipeline`. + +--- + +## Post hoc validation + +- **HDX:** root-mean-square error between observed uptake and posterior mean ensemble prediction *D̂* (`bioce_pipeline.py`). +- **XL:** posterior mean of *p*ₛₐₜ,ₗ per restraint from generated quantities. + +These metrics are diagnostic summaries; formal model comparison (e.g. Bayes factors, LOO) follows the original bioce SAXS workflow where implemented. + +--- + +## Limitations + +1. **HDX:** The default structure-to-uptake mapping uses mean peptide SASA, not hydrogen-exchange kinetics, LC-MS peak integration, or condition-specific Weibull/logistic models (cf. Crook et al., J. Proteome Res. 2023). Users requiring kinetic HDX models should supply precomputed `SimulatedHDX.txt` from external tools. +2. **XL:** Distances use Cα (or fallback) pairs, not full linker chemistry or atom-specific crosslinker geometry. There is no explicit model for undetected crosslinks beyond omitting them from the restraint list. +3. **Ensemble selection:** Variational Bayesian model selection (`variationalBayesian.py`) has not been extended to multimodal HDX/XL likelihoods; users may pre-filter conformers with SAXS-only VB before multimodal inference. + +--- + +## Key references + +1. Potrzebowski W., Trewhella J., Andre I. Bayesian inference of protein conformational ensembles from limited structural data. *PLOS Computational Biology* **14**, e1006641 (2018). +2. Crook O. M., Gittens N., Chung C.-w., Deane C. M. A Functional Bayesian Model for Hydrogen-Deuterium Exchange Mass Spectrometry. *Journal of Proteome Research* **22**, 2959–2972 (2023). *(Related methodology; kinetic HDX layer not implemented in bioce.)* +3. Stan Development Team. Stan modeling language users guide and reference manual (versions cited with PyStan 3 install). +4. BioPython project: PDB parsing and Shrake–Rupley SASA (Cock et al., *Bioinformatics* 2009). + +--- + +*Document version: bioce multimodal extension (HDX/XL). Implementation: `stan_models.py`, `structure_observables.py`, `fullBayesianMultimodal.py`, `multimodal_io.py`.* diff --git a/docs/WORKFLOW.md b/docs/WORKFLOW.md new file mode 100644 index 0000000..b8653cf --- /dev/null +++ b/docs/WORKFLOW.md @@ -0,0 +1,174 @@ +# Bioce multimodal workflow + +End-to-end flow for Bayesian ensemble inference with SAXS, HDX-MS, and/or XL-MS. The same logic runs via **CLI** (`fullBayesianMultimodal.py`) or the **Streamlit app** (`app.py` → `bioce_pipeline.py`). + +See also: [METHODS_HDX_XL.md](METHODS_HDX_XL.md) for equations and [README.md](../README.md) for usage commands. + +--- + +## Overview + +```mermaid +flowchart TB + subgraph inputs [User inputs] + PDB["PDB ensemble\n(zip or files)"] + SAXSexp["SAXS curve\nq, I, sigma"] + HDXpep["HDX peptides\nchain start end"] + HDXexp["HDX uptake\nuptake sigma"] + XLrest["XL restraints\nres_i res_j z"] + end + + subgraph prep [Per-conformer observables] + FoXS["FoXS / Pepsi-SAXS"] + SASA["Peptide SASA\nShrake-Rupley"] + CAdist["C-alpha distances"] + SimSAXS["SimulatedIntensities.txt"] + SimHDX["SimulatedHDX.txt"] + SimXL["SimulatedXLdistances.txt"] + end + + subgraph stan [Stan NUTS inference] + W["Population weights w\nDirichlet prior"] + LikSAXS["SAXS: Normal likelihood"] + LikHDX["HDX: Normal likelihood"] + LikXL["XL: Bernoulli likelihood"] + Joint["Joint posterior\nshared w"] + end + + subgraph outputs [Outputs] + Weights["Posterior mean weights"] + Chi2["SAXS chi-squared"] + RMSE["HDX RMSE"] + Psat["XL p_sat per link"] + Plots["Fit plots and CSV"] + end + + PDB --> FoXS + PDB --> SASA + PDB --> CAdist + SAXSexp --> FoXS + FoXS --> SimSAXS + HDXpep --> SASA + HDXexp --> SimHDX + SASA --> SimHDX + XLrest --> CAdist + CAdist --> SimXL + + SimSAXS --> LikSAXS + SimHDX --> LikHDX + SimXL --> LikXL + SAXSexp --> LikSAXS + HDXexp --> LikHDX + XLrest --> LikXL + + W --> LikSAXS + W --> LikHDX + W --> LikXL + LikSAXS --> Joint + LikHDX --> Joint + LikXL --> Joint + + Joint --> Weights + Joint --> Chi2 + Joint --> RMSE + Joint --> Psat + Joint --> Plots +``` + +--- + +## Modality selection + +Any non-empty subset of SAXS, HDX, and XL is supported. Stan model choice is automatic: + +```mermaid +flowchart LR + subgraph mods [Enabled modalities] + S[SAXS] + H[HDX] + X[XL] + end + + S --> M1[stan_code] + H --> M2[stan_code_HDX] + X --> M3[stan_code_XL] + + S --> SH[stan_code_SAXS_HDX] + H --> SH + S --> SX[stan_code_SAXS_XL] + X --> SX + H --> HX[stan_code_HDX_XL] + X --> HX + S --> SHX[stan_code_SAXS_HDX_XL] + H --> SHX + X --> SHX +``` + +| SAXS | HDX | XL | Script / model | +|:----:|:---:|:--:|----------------| +| | ✓ | | `stan_code_HDX` | +| | | ✓ | `stan_code_XL` | +| ✓ | | | `stan_code` (SAXS-only) | +| ✓ | ✓ | | `stan_code_SAXS_HDX` | +| ✓ | | ✓ | `stan_code_SAXS_XL` | +| | ✓ | ✓ | `stan_code_HDX_XL` | +| ✓ | ✓ | ✓ | `stan_code_SAXS_HDX_XL` | + +--- + +## CLI vs web app + +```mermaid +flowchart TB + subgraph cli [Command line] + P1[prepareBayesian.py / prepareHDX.py / prepareXLMS.py] + P2[fullBayesianMultimodal.py] + P1 --> P2 + end + + subgraph web [Streamlit app.py] + U[Upload PDBs and data files] + BP[bioce_pipeline.run_full_pipeline] + U --> BP + end + + subgraph core [Shared core] + SO[structure_observables.py] + SM[stan_models.py] + SR[stan_run.build_and_sample] + SO --> SM --> SR + end + + P2 --> core + BP --> SO + BP --> SM +``` + +--- + +## HDX and XL at a glance + +```mermaid +flowchart LR + subgraph hdx [HDX-MS] + Dpred["D_pred,i(S_k)\nmean peptide SASA / 200"] + Dhat["D_hat_i = sum w_k D_pred"] + Y["y_i ~ Normal(D_hat, sigma)"] + Dpred --> Dhat --> Y + end + + subgraph xl [XL-MS] + dlk["d_l(S_k)\nCA distance"] + sat["s_l,k = sigmoid((d_max - d)/tau)"] + psat["p_sat = sum w_k s_l,k"] + Z["z_l ~ Bernoulli(p_sat)"] + dlk --> sat --> psat --> Z + end + + w["weights w"] --> Dhat + w --> psat +``` + +--- + +*Render these diagrams in GitHub, VS Code (Mermaid extension), or [mermaid.live](https://mermaid.live).* diff --git a/fullBayesian.py b/fullBayesian.py index 24d0934..57574f4 100755 --- a/fullBayesian.py +++ b/fullBayesian.py @@ -11,12 +11,11 @@ import optparse import numpy as np -import pystan import psisloo -import matplotlib.pyplot as plt import stan_utility +import stan_run -from statistics import calculateChiCrysol, calculateChemShiftsChi, JensenShannonDiv, waic, mean_for_weights, me_log_lik +from bioce_statistics import calculateChiCrysol, calculateChemShiftsChi, JensenShannonDiv, waic, mean_for_weights, me_log_lik from stan_models import stan_code, stan_code_CS, stan_code_EP, stan_code_EP_CS, \ psisloo_quanities, posterior_predict_quanities @@ -39,23 +38,20 @@ def execute_stan(experimental, simulated, priors, iterations, chains, njobs): "n_structures" : np.shape(simulated)[1], "priors":priors} - #sm = pystan.StanModel(model_code=stan_code+psisloo_quanities) - sm = pystan.StanModel(model_code=stan_code) - fit = sm.sampling(data=stan_dat, iter=iterations, chains=chains, - n_jobs=njobs, sample_file="saved_samples.txt") + fit = stan_run.build_and_sample( + stan_code, + stan_dat, + iterations, + chains, + njobs=njobs, + ) - #initial_values = [{"weight[0]":0.05, "weight[1]":0.1, "weight[2]":0.15, + # initial_values = [{"weight[0]":0.05, "weight[1]":0.1, "weight[2]":0.15, # "weight[3]":0.3, "weight[4]":0.4, "scale":1}] - #fit = sm.optimizing(data=stan_dat, init=initial_values, algorithm="BFGS") + # PyStan 3 does not expose optimizing(); use CmdStanPy if needed. - fig = fit.plot(pars="weights") - #ax.set_color_cycle(['red', 'black', 'yellow', 'green', 'blue']) - fig.subplots_adjust(wspace=0.8) - fig.savefig("stan_weights.png", dpi=300) - - fig = fit.plot(pars="scale") - fig.subplots_adjust(wspace=0.8) - fig.savefig("stan_scale.png", dpi=300) + stan_run.plot_parameter_marginals(fit, "weights", "stan_weights.png", dpi=300) + stan_run.plot_parameter_marginals(fit, "scale", "stan_scale.png", dpi=300) #np.savetxt("target_curve_full.csv", fit.summary()['summary'][-869:-1][:,:2]) return fit @@ -79,20 +75,15 @@ def execute_stan_EP(experimental, simulated, priors, iterations, chains, njobs): "n_measures" : np.shape(experimental)[0], "n_structures" : np.shape(simulated)[1], "energy_priors":priors} - sm = pystan.StanModel(model_code=stan_code_EP) - fit = sm.sampling(data=stan_dat, iter=iterations, chains=chains, n_jobs=njobs, sample_file="saved_samples.txt") - - fig = fit.plot(pars="weights") - fig.subplots_adjust(wspace=0.8) - fig.savefig("stan_weights_EP.png", dpi=300) - - fig = fit.plot(pars="boltzmann_shift") - fig.subplots_adjust(wspace=0.8) - fig.savefig("stan_boltzmann_EP.png", dpi=300) + fit = stan_run.build_and_sample( + stan_code_EP, stan_dat, iterations, chains, njobs=njobs + ) - fig = fit.plot(pars="scale") - fig.subplots_adjust(wspace=0.8) - fig.savefig("stan_scale_EP.png", dpi=300) + stan_run.plot_parameter_marginals(fit, "weights", "stan_weights_EP.png", dpi=300) + stan_run.plot_parameter_marginals( + fit, "boltzmann_shift", "stan_boltzmann_EP.png", dpi=300 + ) + stan_run.plot_parameter_marginals(fit, "scale", "stan_scale_EP.png", dpi=300) return fit @@ -126,20 +117,15 @@ def execute_stan_EP_CS(experimental, simulated, priors, "n_structures" : np.shape(simulated)[1], "energy_priors":priors} - sm = pystan.StanModel(model_code=stan_code_EP_CS) - fit = sm.sampling(data=stan_dat, iter=iterations, chains=chains, n_jobs=njobs, sample_file="saved_samples.txt") - - fig = fit.plot(pars="weights") - fig.subplots_adjust(wspace=0.8) - fig.savefig("stan_weights_EP_CS.png", dpi=300) - - fig = fit.plot(pars="boltzmann_shift") - fig.subplots_adjust(wspace=0.8) - fig.savefig("stan_boltzmann_EP_CS.png", dpi=300) + fit = stan_run.build_and_sample( + stan_code_EP_CS, stan_dat, iterations, chains, njobs=njobs + ) - fig = fit.plot(pars="scale") - fig.subplots_adjust(wspace=0.8) - fig.savefig("stan_scale_EP_CS.png", dpi=300) + stan_run.plot_parameter_marginals(fit, "weights", "stan_weights_EP_CS.png", dpi=300) + stan_run.plot_parameter_marginals( + fit, "boltzmann_shift", "stan_boltzmann_EP_CS.png", dpi=300 + ) + stan_run.plot_parameter_marginals(fit, "scale", "stan_scale_EP_CS.png", dpi=300) return fit @@ -172,16 +158,12 @@ def execute_stan_CS(experimental, simulated, priors, "n_structures" : np.shape(simulated)[1], "priors":priors} - sm = pystan.StanModel(model_code=stan_code_CS) - fit = sm.sampling(data=stan_dat, iter=iterations, chains=chains, n_jobs=njobs, sample_file="saved_samples.txt") + fit = stan_run.build_and_sample( + stan_code_CS, stan_dat, iterations, chains, njobs=njobs + ) - fig = fit.plot(pars="weights") - fig.subplots_adjust(wspace=0.8) - fig.savefig("stan_weights_CS.png", dpi=300) - - fig = fit.plot(pars="scale") - fig.subplots_adjust(wspace=0.8) - fig.savefig("stan_scale.png", dpi=300) + stan_run.plot_parameter_marginals(fit, "weights", "stan_weights_CS.png", dpi=300) + stan_run.plot_parameter_marginals(fit, "scale", "stan_scale.png", dpi=300) return fit def execute_bws(experimental, simulated, priors, file_names, threshold, @@ -215,7 +197,7 @@ def execute_bws(experimental, simulated, priors, file_names, threshold, last_loo = None last_waic = None model_comp_diff = 1 - sm = pystan.StanModel(model_code=stan_code+psisloo_quanities) + model_code_bws = stan_code + psisloo_quanities iteration = 0 repeat_iteration = 0 @@ -234,15 +216,16 @@ def execute_bws(experimental, simulated, priors, file_names, threshold, "n_structures" : n_structures, "priors" : alphas} - fit = sm.sampling(data=stan_dat, iter=iterations, chains=chains, - n_jobs=njobs) + fit = stan_run.build_and_sample( + model_code_bws, stan_dat, iterations, chains, njobs=njobs + ) #Calculating psis loo - stan_chain=fit.extract() + loglikes = np.asarray(fit["loglikes"]).T #If less than 100 models start psisloo otherwise waic if n_structures < structure_cutoff: - current_loo = psisloo.psisloo(stan_chain['loglikes']) + current_loo = psisloo.psisloo(loglikes) log_file.write("greater than 0.5 ") log_file.write(str(current_loo.print_summary()[0])+"\n") log_file.write("greater than 1 ") @@ -254,7 +237,7 @@ def execute_bws(experimental, simulated, priors, file_names, threshold, log_file.write(str(model_comp['diff'])+"\n") log_file.write(str(model_comp['se_diff'])+"\n") else: - current_waic = waic(stan_chain['loglikes']) + current_waic = waic(loglikes) if last_waic: model_comp_diff = current_waic - last_waic log_file.write("WAIC model comparison: \n") @@ -279,15 +262,15 @@ def execute_bws(experimental, simulated, priors, file_names, threshold, else: last_waic = current_waic repeat_iteration = 0 - current_weights = fit.summary()['summary'][:,0][:n_structures] + current_weights = np.asarray(fit["weights"]).mean(axis=1) sim_curves = sim_curves[:,current_weights>threshold] alphas = alphas[current_weights>threshold] n_structures = np.shape(sim_curves)[1] file_names = file_names[current_weights>threshold] #np.savetxt("fit.txt",fit.summary()['summary'][:,0],delimiter=" ") - fig = fit.plot() - fig.subplots_adjust(wspace=0.8) - fig.savefig("stan_fit_"+str(iteration)+".png", dpi=300) + stan_run.plot_parameter_marginals( + fit, "weights", "stan_fit_" + str(iteration) + ".png", dpi=300 + ) bayesian_weights, jsd, crysol_chi2 = calculate_stats(fit, experimental, simulated) log_file.write("JSD: "+str(jsd)) @@ -329,32 +312,23 @@ def execute_bws(experimental, simulated, priors, file_names, threshold, def calculate_stats(fit, experimental, simulated, cs_simulated=None, cs_rms=None, cs_experimental=None): """ - Calculates statistics based on stan model + Calculates statistics based on stan model (PyStan 3 ``stan`` fit). :param fit: :return: """ - #la = fit.extract(permuted=True) # return a dictionary of arrays - #mu = la['weights'] - - ## return an array of three dimensions: iterations, chains, parameters - results_array = fit.extract(permuted=False, inc_warmup=False) + weights_draws = np.asarray(fit["weights"]) + if weights_draws.ndim == 1: + weights_draws = weights_draws.reshape(1, -1) + n_struct = np.shape(simulated)[1] + weights_draws = weights_draws[:n_struct, :] + n_draws = weights_draws.shape[1] + bayesian_weights = weights_draws.mean(axis=1) - nsamples = 0 jsd_sum = 0.0 - bayesian_weights = np.zeros(np.shape(simulated)[1]) - for iteration in results_array: - for parameters in iteration: - current_weights = parameters[:np.shape(simulated)[1]] - bayesian_weights+=current_weights - nsamples+=1 - bayesian_weights=bayesian_weights/nsamples - - for iteration in results_array: - for parameters in iteration: - current_weights = parameters[:np.shape(simulated)[1]] - jsd_sum+=JensenShannonDiv(current_weights, bayesian_weights) - jsd = (np.sqrt(jsd_sum/nsamples)) - + for j in range(n_draws): + current_weights = weights_draws[:, j] + jsd_sum += JensenShannonDiv(current_weights, bayesian_weights) + jsd = np.sqrt(jsd_sum / n_draws) crysol_chi2 = calculateChiCrysol(np.dot(bayesian_weights, np.transpose(simulated)), experimental[:,1], @@ -365,8 +339,8 @@ def calculate_stats(fit, experimental, simulated, cs_simulated=None, except: chemical_shifts_on = False - scale = fit.summary(pars='scale')['summary'][0][0] - print("Scale from summary", scale) + scale = float(np.mean(np.asarray(fit["scale"]))) + print("Scale from posterior mean", scale) combine_curve(experimental, simulated, bayesian_weights, scale) if chemical_shifts_on: diff --git a/fullBayesianMultimodal.py b/fullBayesianMultimodal.py new file mode 100644 index 0000000..1c773c3 --- /dev/null +++ b/fullBayesianMultimodal.py @@ -0,0 +1,342 @@ +""" +Bayesian ensemble inference with HDX-MS and/or XL-MS likelihoods (PyStan 3). + +Combines with existing SAXS inputs using the same file layouts as fullBayesian.py. +""" +from __future__ import print_function + +import optparse +import os + +import numpy as np +import stan_run +import stan_utility +from bioce_statistics import JensenShannonDiv, mean_for_weights +from multimodal_io import ( + check_sim_matrix, + read_hdx_experimental, + read_matrix, + read_xl_restraints, +) +from stan_models import ( + stan_code_HDX, + stan_code_HDX_XL, + stan_code_SAXS_HDX, + stan_code_SAXS_HDX_XL, + stan_code_SAXS_XL, + stan_code_XL, + stan_code_XL_FP, + xl_posterior_predict, +) + + +def _plot_weights(fit, prefix="stan"): + stan_run.plot_parameter_marginals(fit, "weights", "{}_weights.png".format(prefix), dpi=300) + + +def execute_stan_hdx(hdx_exp, sim_hdx, priors, iterations, chains, njobs): + target_hdx, target_hdxerr = read_hdx_experimental(hdx_exp) + sim = read_matrix(sim_hdx) + n_structures = check_sim_matrix(sim, target_hdx.size, "HDX") + stan_dat = { + "n_obs": int(target_hdx.size), + "n_structures": int(n_structures), + "target_hdx": target_hdx, + "target_hdxerr": target_hdxerr, + "sim_hdx": sim, + "priors": priors, + } + print("Starting HDX-only inference ({} obs, {} structures)".format( + stan_dat["n_obs"], n_structures + )) + fit = stan_run.build_and_sample( + stan_code_HDX, stan_dat, iterations, chains, njobs=njobs + ) + _plot_weights(fit, "hdx") + return fit + + +def execute_stan_xl(xl_restraints, sim_xl, priors, iterations, chains, njobs, use_fp=False): + xl_obs, d_max, tau = read_xl_restraints(xl_restraints) + sim = read_matrix(sim_xl) + n_structures = check_sim_matrix(sim, xl_obs.size, "XL") + stan_dat = { + "n_xl": int(xl_obs.size), + "n_structures": int(n_structures), + "xl_obs": xl_obs, + "d_max": d_max, + "tau": tau, + "xl_distances": sim, + "priors": priors, + } + code = stan_code_XL_FP if use_fp else stan_code_XL + xl_posterior_predict + print("Starting XL-only inference ({} links, {} structures)".format( + stan_dat["n_xl"], n_structures + )) + fit = stan_run.build_and_sample(code, stan_dat, iterations, chains, njobs=njobs) + _plot_weights(fit, "xl") + return fit + + +def execute_stan_saxs_hdx( + experimental, simulated, hdx_exp, sim_hdx, priors, iterations, chains, njobs +): + target_hdx, target_hdxerr = read_hdx_experimental(hdx_exp) + sim_h = read_matrix(sim_hdx) + check_sim_matrix(sim_h, target_hdx.size, "HDX") + stan_dat = { + "n_measures": int(experimental.shape[0]), + "n_hdx": int(target_hdx.size), + "n_structures": int(simulated.shape[1]), + "target_saxs": experimental[:, 1], + "target_saxserr": experimental[:, 2], + "sim_saxs": simulated, + "target_hdx": target_hdx, + "target_hdxerr": target_hdxerr, + "sim_hdx": sim_h, + "priors": priors, + } + if sim_h.shape[1] != stan_dat["n_structures"]: + raise ValueError("SAXS and HDX simulated matrices have different structure counts") + print("Starting SAXS+HDX inference") + fit = stan_run.build_and_sample( + stan_code_SAXS_HDX, stan_dat, iterations, chains, njobs=njobs + ) + _plot_weights(fit, "saxs_hdx") + stan_run.plot_parameter_marginals(fit, "scale", "saxs_hdx_scale.png", dpi=300) + return fit + + +def execute_stan_saxs_xl( + experimental, + simulated, + xl_restraints, + sim_xl, + priors, + iterations, + chains, + njobs, +): + xl_obs, d_max, tau = read_xl_restraints(xl_restraints) + sim_x = read_matrix(sim_xl) + check_sim_matrix(sim_x, xl_obs.size, "XL") + stan_dat = { + "n_measures": int(experimental.shape[0]), + "n_xl": int(xl_obs.size), + "n_structures": int(simulated.shape[1]), + "target_saxs": experimental[:, 1], + "target_saxserr": experimental[:, 2], + "sim_saxs": simulated, + "xl_obs": xl_obs, + "d_max": d_max, + "tau": tau, + "xl_distances": sim_x, + "priors": priors, + } + if sim_x.shape[1] != stan_dat["n_structures"]: + raise ValueError("SAXS and XL simulated matrices have different structure counts") + print("Starting SAXS+XL inference") + code = stan_code_SAXS_XL + xl_posterior_predict + fit = stan_run.build_and_sample(code, stan_dat, iterations, chains, njobs=njobs) + _plot_weights(fit, "saxs_xl") + stan_run.plot_parameter_marginals(fit, "scale", "saxs_xl_scale.png", dpi=300) + return fit + + +def execute_stan_saxs_hdx_xl( + experimental, + simulated, + hdx_exp, + sim_hdx, + xl_restraints, + sim_xl, + priors, + iterations, + chains, + njobs, +): + target_hdx, target_hdxerr = read_hdx_experimental(hdx_exp) + xl_obs, d_max, tau = read_xl_restraints(xl_restraints) + sim_h = read_matrix(sim_hdx) + sim_x = read_matrix(sim_xl) + n_structures = int(simulated.shape[1]) + check_sim_matrix(sim_h, target_hdx.size, "HDX") + check_sim_matrix(sim_x, xl_obs.size, "XL") + if sim_h.shape[1] != n_structures or sim_x.shape[1] != n_structures: + raise ValueError("All simulated matrices must have the same number of columns") + stan_dat = { + "n_measures": int(experimental.shape[0]), + "n_hdx": int(target_hdx.size), + "n_xl": int(xl_obs.size), + "n_structures": n_structures, + "target_saxs": experimental[:, 1], + "target_saxserr": experimental[:, 2], + "sim_saxs": simulated, + "target_hdx": target_hdx, + "target_hdxerr": target_hdxerr, + "sim_hdx": sim_h, + "xl_obs": xl_obs, + "d_max": d_max, + "tau": tau, + "xl_distances": sim_x, + "priors": priors, + } + print("Starting SAXS+HDX+XL inference") + code = stan_code_SAXS_HDX_XL + xl_posterior_predict + fit = stan_run.build_and_sample(code, stan_dat, iterations, chains, njobs=njobs) + _plot_weights(fit, "saxs_hdx_xl") + stan_run.plot_parameter_marginals(fit, "scale", "saxs_hdx_xl_scale.png", dpi=300) + return fit + + +def execute_stan_hdx_xl(hdx_exp, sim_hdx, xl_restraints, sim_xl, priors, iterations, chains, njobs): + target_hdx, target_hdxerr = read_hdx_experimental(hdx_exp) + xl_obs, d_max, tau = read_xl_restraints(xl_restraints) + sim_h = read_matrix(sim_hdx) + sim_x = read_matrix(sim_xl) + n_structures = check_sim_matrix(sim_h, target_hdx.size, "HDX") + if sim_x.shape[1] != n_structures: + raise ValueError("HDX and XL simulated matrices have different structure counts") + check_sim_matrix(sim_x, xl_obs.size, "XL") + stan_dat = { + "n_hdx": int(target_hdx.size), + "n_xl": int(xl_obs.size), + "n_structures": int(n_structures), + "target_hdx": target_hdx, + "target_hdxerr": target_hdxerr, + "sim_hdx": sim_h, + "xl_obs": xl_obs, + "d_max": d_max, + "tau": tau, + "xl_distances": sim_x, + "priors": priors, + } + print("Starting HDX+XL inference") + code = stan_code_HDX_XL + xl_posterior_predict + fit = stan_run.build_and_sample(code, stan_dat, iterations, chains, njobs=njobs) + _plot_weights(fit, "hdx_xl") + return fit + + +def print_weight_summary(fit, file_names=None): + weights = mean_for_weights(fit) + if file_names is not None: + for index, fname in enumerate(file_names): + print(fname, weights[index]) + else: + print("Posterior mean weights:", weights) + return weights + + +def read_file_safe(filename, dtype="float64"): + try: + return np.genfromtxt(filename, dtype=dtype) + except IOError as err: + print(os.strerror(err.errno)) + raise + + +if __name__ == "__main__": + doc = """ + Bayesian inference with HDX-MS and/or XL-MS (optionally with SAXS). + Usage: python fullBayesianMultimodal.py --help + """ + print(doc) + usage = "usage: %prog [options]" + parser = optparse.OptionParser(usage=usage, version="0.1") + + parser.add_option("-p", "--priors", dest="priors_file", default=None, + help="Dirichlet prior weights (one per structure)") + parser.add_option("-f", "--file_names", dest="names_file", default=None, + help="Structure names (optional, for printed output)") + parser.add_option("-i", "--iterations", dest="iterations", default=2000, type="int") + parser.add_option("-j", "--jobs", dest="njobs", default=4, type="int") + parser.add_option("-c", "--chains", dest="chains", default=4, type="int") + + parser.add_option("-e", "--experimental", dest="experimental_file", default=None, + help="Experimental SAXS (q, I, sigma) — required for SAXS modes") + parser.add_option("-s", "--simulated", dest="simulated_file", default=None, + help="Simulated SAXS matrix") + + parser.add_option("-H", "--hdx_experimental", dest="hdx_exp", default=None, + help="HDX uptake file (uptake sigma per line)") + parser.add_option("-S", "--hdx_simulated", dest="hdx_sim", default=None, + help="Simulated HDX matrix") + + parser.add_option("-X", "--xl_restraints", dest="xl_restraints", default=None, + help="Crosslink restraint list (res_i res_j z [d_max] [tau])") + parser.add_option("-D", "--xl_distances", dest="xl_sim", default=None, + help="Simulated XL distance matrix") + parser.add_option("--xl-fp-mixture", dest="xl_fp", action="store_true", default=False, + help="XL-only: use false-positive mixture model") + + options, _ = parser.parse_args() + + if not options.priors_file: + parser.error("-p/--priors is required") + + priors = read_file_safe(options.priors_file) + file_names = None + if options.names_file: + file_names = read_file_safe(options.names_file, dtype="unicode") + + has_saxs = options.experimental_file and options.simulated_file + has_hdx = options.hdx_exp and options.hdx_sim + has_xl = options.xl_restraints and options.xl_sim + + if not (has_hdx or has_xl): + parser.error("Provide HDX (-H/-S) and/or XL (-X/-D) inputs") + + iterations = options.iterations + chains = options.chains + njobs = options.njobs + + if has_saxs and has_hdx and has_xl: + experimental = read_file_safe(options.experimental_file) + simulated = read_file_safe(options.simulated_file) + fit = execute_stan_saxs_hdx_xl( + experimental, simulated, + options.hdx_exp, options.hdx_sim, + options.xl_restraints, options.xl_sim, + priors, iterations, chains, njobs, + ) + elif has_saxs and has_hdx: + experimental = read_file_safe(options.experimental_file) + simulated = read_file_safe(options.simulated_file) + fit = execute_stan_saxs_hdx( + experimental, simulated, + options.hdx_exp, options.hdx_sim, + priors, iterations, chains, njobs, + ) + elif has_saxs and has_xl: + experimental = read_file_safe(options.experimental_file) + simulated = read_file_safe(options.simulated_file) + fit = execute_stan_saxs_xl( + experimental, simulated, + options.xl_restraints, options.xl_sim, + priors, iterations, chains, njobs, + ) + elif has_hdx and has_xl: + fit = execute_stan_hdx_xl( + options.hdx_exp, options.hdx_sim, + options.xl_restraints, options.xl_sim, + priors, iterations, chains, njobs, + ) + elif has_hdx: + fit = execute_stan_hdx( + options.hdx_exp, options.hdx_sim, + priors, iterations, chains, njobs, + ) + elif has_xl: + fit = execute_stan_xl( + options.xl_restraints, options.xl_sim, + priors, iterations, chains, njobs, + use_fp=options.xl_fp, + ) + else: + parser.error("Unrecognized modality combination") + + print_weight_summary(fit, file_names) + stan_utility.check_treedepth(fit) + stan_utility.check_energy(fit) + stan_utility.check_div(fit) diff --git a/fullBayesianTR.py b/fullBayesianTR.py index a1f12c4..082e20d 100644 --- a/fullBayesianTR.py +++ b/fullBayesianTR.py @@ -9,15 +9,13 @@ import os import optparse -import pickle import numpy as np -import pystan import psisloo -import matplotlib.pyplot as plt import stan_utility +import stan_run -from statistics import ( +from bioce_statistics import ( calculateChiCrysol, calculateChemShiftsChi, JensenShannonDiv, @@ -164,14 +162,8 @@ def execute_stan(experimental, simulated, priors, iterations, chains, njobs): "priors": priors, } - # sm = pystan.StanModel(model_code=stan_code+psisloo_quanities) - sm = pystan.StanModel(model_code=stan_code) - fit = sm.sampling( - data=stan_dat, - iter=iterations, - chains=chains, - n_jobs=njobs, - sample_file="saved_samples.txt", + fit = stan_run.build_and_sample( + stan_code, stan_dat, iterations, chains, njobs=njobs ) # initial_values = [{"weight[0]":0.05, "weight[1]":0.1, "weight[2]":0.15, @@ -202,24 +194,19 @@ def calculate_stats( # la = fit.extract(permuted=True) # return a dictionary of arrays # mu = la['weights'] - ## return an array of three dimensions: iterations, chains, parameters - results_array = fit.extract(permuted=False, inc_warmup=False) + weights_draws = np.asarray(fit["weights"]) + if weights_draws.ndim == 1: + weights_draws = weights_draws.reshape(1, -1) + n_struct = np.shape(simulated)[1] + weights_draws = weights_draws[:n_struct, :] + n_draws = weights_draws.shape[1] + bayesian_weights = weights_draws.mean(axis=1) - nsamples = 0 jsd_sum = 0.0 - bayesian_weights = np.zeros(np.shape(simulated)[1]) - for iteration in results_array: - for parameters in iteration: - current_weights = parameters[: np.shape(simulated)[1]] - bayesian_weights += current_weights - nsamples += 1 - bayesian_weights = bayesian_weights / nsamples - - for iteration in results_array: - for parameters in iteration: - current_weights = parameters[: np.shape(simulated)[1]] - jsd_sum += JensenShannonDiv(current_weights, bayesian_weights) - jsd = np.sqrt(jsd_sum / nsamples) + for j in range(n_draws): + current_weights = weights_draws[:, j] + jsd_sum += JensenShannonDiv(current_weights, bayesian_weights) + jsd = np.sqrt(jsd_sum / n_draws) crysol_chi2 = calculateChiCrysol( np.dot(bayesian_weights, np.transpose(simulated)), @@ -388,12 +375,6 @@ def calculatePostChi(): # Even distribution of priors priors = np.ones(number_of_states + 1) * 1.0 / (number_of_states + 1) file_names = np.append(file_names, "diff") - # Pickling for teh first time compilation - sm = pystan.StanModel(model_code=stan_code) - # with open('model.pkl', 'wb') as f: - # pickle.dump(sm, f) - # sm = pickle.load(open('model.pkl', 'rb')) - iteration_folder = "iteration_" + str(options.iteration) # previous_iteration_folder = 'iteration_' + str(options.iteration-1) if os.path.exists(iteration_folder) == False: @@ -404,8 +385,6 @@ def calculatePostChi(): # fout.write('jsd, chi2\n') # List of experiemntal files to scan over - stan_initialized = False - # TODO: For each frame load diff file and combine it with simulated # What to do with CTD simulated? # make if statement for non-diff cases @@ -452,8 +431,6 @@ def calculatePostChi(): # TODO: Priors should be defined as a matrix and defined for specific matrix stan_dat = { "sim_curves": full_simulated, - # "ctd_curves" : ctd_simulated, - "qvector": qvector, "target_curve": experimental[:, 1], "target_errors": target_errors, "n_measures": np.shape(experimental)[0], @@ -461,8 +438,13 @@ def calculatePostChi(): "priors": priors_i, } - fit = sm.sampling( - data=stan_dat, iter=iterations, chains=chains, n_jobs=njobs, refresh=-1 + fit = stan_run.build_and_sample( + stan_code, + stan_dat, + iterations, + chains, + njobs=njobs, + show_progress=False, ) print(fit) diff --git a/hdx_weights.png b/hdx_weights.png new file mode 100644 index 0000000..2e0a245 Binary files /dev/null and b/hdx_weights.png differ diff --git a/multimodal_io.py b/multimodal_io.py new file mode 100644 index 0000000..5d7e941 --- /dev/null +++ b/multimodal_io.py @@ -0,0 +1,87 @@ +""" +I/O helpers for HDX-MS and XL-MS inputs to multimodal Stan models. + +File formats (same column layout as SAXS simulated matrices: rows = observations, +columns = conformers): + + HDX experimental (hdx_exp.dat): uptake sigma [per flattened peptide×time obs] + HDX simulated (SimulatedHDX.txt): one row per observation, space-separated columns + + XL restraints (xl_restraints.dat): + res_i res_j z [d_max] [tau] + z = 1 observed crosslink, 0 absent / decoy; d_max and tau optional (defaults 30, 3 Å). + + XL distances (SimulatedXLdistances.txt): rows = crosslinks, cols = conformers (Cα or + linker-compatible distance in Å). +""" +from __future__ import annotations + +import numpy as np + +DEFAULT_XL_DMAX = 30.0 +DEFAULT_XL_TAU = 3.0 + + +def read_matrix(path: str, dtype="float64") -> np.ndarray: + return np.genfromtxt(path, dtype=dtype) + + +def read_hdx_experimental(path: str) -> tuple[np.ndarray, np.ndarray]: + """Return (uptake, sigma) vectors.""" + data = read_matrix(path) + if data.ndim == 1: + if data.size % 2 != 0: + raise ValueError("HDX experimental file must have uptake and sigma columns") + data = data.reshape(-1, 2) + if data.shape[1] < 2: + raise ValueError("HDX experimental file needs at least two columns: uptake sigma") + return data[:, 0], data[:, 1] + + +def read_xl_restraints(path: str) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """ + Parse xl_restraints.dat. + + Returns (xl_obs, d_max, tau) as 1-D float/int arrays of length n_xl. + """ + rows = [] + with open(path) as handle: + for line in handle: + line = line.strip() + if not line or line.startswith("#"): + continue + rows.append(line.split()) + if not rows: + raise ValueError("No crosslink restraints found in {}".format(path)) + + xl_obs = [] + d_max = [] + tau = [] + for parts in rows: + if len(parts) < 3: + raise ValueError( + "Each XL restraint line needs: res_i res_j z [d_max] [tau]: {}".format( + parts + ) + ) + xl_obs.append(int(float(parts[2]))) + d_max.append(float(parts[3]) if len(parts) > 3 else DEFAULT_XL_DMAX) + tau.append(float(parts[4]) if len(parts) > 4 else DEFAULT_XL_TAU) + return ( + np.asarray(xl_obs, dtype=np.int32), + np.asarray(d_max, dtype=np.float64), + np.asarray(tau, dtype=np.float64), + ) + + +def check_sim_matrix(sim: np.ndarray, n_obs: int, label: str) -> int: + """Ensure simulated matrix row count matches observations; return n_structures.""" + if sim.ndim != 2: + raise ValueError("{} simulated matrix must be 2-D".format(label)) + if sim.shape[0] != n_obs: + raise ValueError( + "{} simulated rows ({}) != {} observations ({})".format( + label, sim.shape[0], label, n_obs + ) + ) + return sim.shape[1] diff --git a/prepareHDX.py b/prepareHDX.py new file mode 100644 index 0000000..48b8e17 --- /dev/null +++ b/prepareHDX.py @@ -0,0 +1,90 @@ +#! /usr/bin/env python +""" +Build HDX simulated matrix for multimodal Bayesian inference. + +Expected per-structure uptake files in a directory (one file per conformer), each with +the same number of lines (flattened peptide×time observations). Lines are either +``uptake`` or ``time uptake``. + +Usage: python prepareHDX.py --help +""" +from __future__ import print_function + +import os +import optparse + +import numpy as np + + +def _read_uptake_column(path): + data = np.genfromtxt(path, dtype=np.float64) + if data.ndim == 1: + return data + if data.shape[1] == 1: + return data[:, 0] + if data.shape[1] >= 2: + return data[:, 1] + raise ValueError("Cannot parse uptake from {}".format(path)) + + +def build_hdx_matrix(dir_name, suffix=".hdx", out_sim="SimulatedHDX.txt"): + files = sorted( + f + for f in os.listdir(dir_name) + if f.endswith(suffix) or (suffix == ".hdx" and f.endswith(".txt")) + ) + if not files: + raise RuntimeError( + "No *{} files in {}; use --suffix to match your naming".format( + suffix, dir_name + ) + ) + + columns = [] + n_obs = None + for fname in files: + col = _read_uptake_column(os.path.join(dir_name, fname)) + if n_obs is None: + n_obs = col.size + elif col.size != n_obs: + raise ValueError( + "Inconsistent observation count in {} ({} vs {})".format( + fname, col.size, n_obs + ) + ) + columns.append(col) + + sim = np.column_stack(columns) + np.savetxt(out_sim, sim, fmt="%.8g") + print("Wrote {} ({} observations × {} structures)".format( + out_sim, sim.shape[0], sim.shape[1] + )) + return out_sim + + +if __name__ == "__main__": + usage = "usage: %prog [options]" + parser = optparse.OptionParser(usage=usage, version="0.1") + parser.add_option( + "-d", + "--hdx_dir", + dest="dir_name", + help="Directory with per-structure predicted HDX uptake files", + ) + parser.add_option( + "-o", + "--output", + dest="output", + default="SimulatedHDX.txt", + help="Output simulated matrix [default: SimulatedHDX.txt]", + ) + parser.add_option( + "--suffix", + dest="suffix", + default=".hdx", + help="Filename suffix for structure files [default: .hdx]", + ) + options, _ = parser.parse_args() + if not options.dir_name: + parser.error("-d/--hdx_dir is required") + build_hdx_matrix(options.dir_name, suffix=options.suffix, out_sim=options.output) diff --git a/prepareXLMS.py b/prepareXLMS.py new file mode 100644 index 0000000..ab8de96 --- /dev/null +++ b/prepareXLMS.py @@ -0,0 +1,87 @@ +#! /usr/bin/env python +""" +Build XL-MS distance matrix for multimodal Bayesian inference. + +Each per-structure file lists one distance per crosslink (same order as +xl_restraints.dat). Lines: ``distance`` or ``res_i res_j distance``. + +Usage: python prepareXLMS.py --help +""" +from __future__ import print_function + +import os +import optparse + +import numpy as np + + +def _read_distance_column(path): + data = np.genfromtxt(path, dtype=np.float64) + if data.ndim == 0: + return np.asarray([float(data)]) + if data.ndim == 1: + return data + return data[:, -1] + + +def build_xl_matrix(dir_name, suffix=".xldist", out_sim="SimulatedXLdistances.txt"): + files = sorted( + f + for f in os.listdir(dir_name) + if f.endswith(suffix) or (suffix == ".xldist" and f.endswith(".txt")) + ) + if not files: + raise RuntimeError( + "No *{} files in {}; use --suffix to match your naming".format( + suffix, dir_name + ) + ) + + columns = [] + n_xl = None + for fname in files: + col = _read_distance_column(os.path.join(dir_name, fname)) + if n_xl is None: + n_xl = col.size + elif col.size != n_xl: + raise ValueError( + "Inconsistent crosslink count in {} ({} vs {})".format( + fname, col.size, n_xl + ) + ) + columns.append(col) + + sim = np.column_stack(columns) + np.savetxt(out_sim, sim, fmt="%.8g") + print("Wrote {} ({} crosslinks × {} structures)".format( + out_sim, sim.shape[0], sim.shape[1] + )) + return out_sim + + +if __name__ == "__main__": + usage = "usage: %prog [options]" + parser = optparse.OptionParser(usage=usage, version="0.1") + parser.add_option( + "-d", + "--xl_dir", + dest="dir_name", + help="Directory with per-structure crosslink distance files", + ) + parser.add_option( + "-o", + "--output", + dest="output", + default="SimulatedXLdistances.txt", + help="Output distance matrix [default: SimulatedXLdistances.txt]", + ) + parser.add_option( + "--suffix", + dest="suffix", + default=".xldist", + help="Filename suffix [default: .xldist]", + ) + options, _ = parser.parse_args() + if not options.dir_name: + parser.error("-d/--xl_dir is required") + build_xl_matrix(options.dir_name, suffix=options.suffix, out_sim=options.output) diff --git a/psisloo.py b/psisloo.py index 8bbe0fc..6a8d001 100644 --- a/psisloo.py +++ b/psisloo.py @@ -80,11 +80,11 @@ def psisloo(log_likelihood): and approximate Leave-One-Out cross-validation (LOO). Takes as input an ndarray of posterior log likelihood terms [ p( y_i | theta^s ) ] - per observation unit. + per observation unit (rows = posterior draws, cols = observations). - e.x. if using pystan: + Example with PyStan 3 (``import stan``):: - loosummary = stanity.psisloo(stan_fit.extract()['log_lik']) + loosummary = psisloo.psisloo(stan_fit["loglikes"].T) Returns a Psisloo object. Useful methods such as print_summary() & plot(). diff --git a/pytest.ini b/pytest.ini new file mode 100644 index 0000000..3082da9 --- /dev/null +++ b/pytest.ini @@ -0,0 +1,6 @@ +[pytest] +pythonpath = . +testpaths = tests +filterwarnings = + ignore::DeprecationWarning:httpstan + ignore::DeprecationWarning:importlib.resources._legacy diff --git a/scilifelab_theme.py b/scilifelab_theme.py new file mode 100644 index 0000000..1095f5a --- /dev/null +++ b/scilifelab_theme.py @@ -0,0 +1,149 @@ +""" +SciLifeLab visual identity colors and UI helpers. + +Brand palette (https://www.scilifelab.se/visual-identity): + Lime #A7C947 — main + Teal #045C64 + Aqua #4C979F + Grape #491F53 — accent (use sparingly) +""" +from __future__ import annotations + +import matplotlib as mpl + +LIME = "#A7C947" +TEAL = "#045C64" +AQUA = "#4C979F" +GRAPE = "#491F53" +WHITE = "#FFFFFF" +GRAY_LIGHT = "#F4F6F4" +GRAY_MID = "#6B7280" +GRAY_DARK = "#2D3436" + +PALETTE = [TEAL, LIME, AQUA, GRAPE, "#7BA428", "#2E7D84"] + +STREAMLIT_CSS = """ + +""" + + +def _build_streamlit_css() -> str: + return ( + STREAMLIT_CSS.replace("__TEAL__", TEAL) + .replace("__AQUA__", AQUA) + .replace("__LIME__", LIME) + .replace("__WHITE__", WHITE) + .replace("__GRAY__", GRAY_MID) + .replace("__GRAY_BG__", GRAY_LIGHT) + ) + + +def inject_streamlit_theme() -> None: + import streamlit as st + + st.markdown(_build_streamlit_css(), unsafe_allow_html=True) + + +def render_header(title: str, subtitle: str) -> None: + import streamlit as st + + st.markdown( + '

{}

{}

'.format(title, subtitle), + unsafe_allow_html=True, + ) + + +def render_footer() -> None: + import streamlit as st + + st.markdown( + '
Bioce · Bayesian ensemble inference · ' + 'SciLifeLab
'.format( + teal=TEAL + ), + unsafe_allow_html=True, + ) + + +def apply_matplotlib_theme() -> None: + mpl.rcParams.update({ + "axes.prop_cycle": mpl.cycler(color=PALETTE), + "axes.edgecolor": GRAY_MID, + "axes.labelcolor": GRAY_DARK, + "axes.titlecolor": TEAL, + "axes.titlesize": 12, + "axes.titleweight": "bold", + "figure.facecolor": WHITE, + "axes.facecolor": WHITE, + "grid.color": "#E5EBE5", + "grid.alpha": 0.7, + "font.family": "sans-serif", + "font.sans-serif": ["Lato", "DejaVu Sans", "Arial"], + "legend.framealpha": 0.9, + }) diff --git a/serve/Dockerfile b/serve/Dockerfile new file mode 100644 index 0000000..fd7c79f --- /dev/null +++ b/serve/Dockerfile @@ -0,0 +1,25 @@ +# Bioce Streamlit app for SciLifeLab Serve (or local Docker) +FROM condaforge/mambaforge:24.7.1-0 + +WORKDIR /app +COPY . /app + +RUN mamba env create -f serve/environment.yml && mamba clean -ay +SHELL ["mamba", "run", "-n", "bioce-serve", "/bin/bash", "-c"] + +# Optional: install ATSAS FoXS for SAXS (uncomment if your image may fetch ATSAS) +# RUN apt-get update && apt-get install -y wget && \ +# wget -q https://www.bioisis.net/download/ATSAS/atsas-3.2.1-1_amd64.deb && \ +# apt-get install -y ./atsas-3.2.1-1_amd64.deb || true + +ENV PATH="/opt/conda/envs/bioce-serve/bin:$PATH" +ENV PYTHONPATH=/app + +EXPOSE 8501 +HEALTHCHECK CMD curl --fail http://localhost:8501/_stcore/health || exit 1 + +CMD ["streamlit", "run", "app.py", \ + "--server.port=8501", \ + "--server.address=0.0.0.0", \ + "--server.headless=true", \ + "--browser.gatherUsageStats=false"] diff --git a/serve/SERVE.md b/serve/SERVE.md new file mode 100644 index 0000000..ed0a2a1 --- /dev/null +++ b/serve/SERVE.md @@ -0,0 +1,26 @@ +# Deploying Bioce on SciLifeLab Serve + +1. Register at [SciLifeLab Serve](https://serve.scilifelab.se/) and create a **Project**. +2. Build and push the Docker image (or let Serve build from this repo): + +```bash +docker build -f serve/Dockerfile -t bioce-app . +docker run -p 8501:8501 bioce-app +``` + +3. In Serve, create an app of type **Streamlit**, point to this image, port **8501**. + +## FoXS / SAXS + +SAXS simulation from PDBs requires **FoXS** (ATSAS) or **Pepsi-SAXS** on `PATH`. The default Dockerfile does not install ATSAS; add it to the image or run **HDX + XL only** until FoXS is available. + +## Example uploads + +See `serve/examples/` for HDX peptide list, uptake, and XL restraint templates. Use your own PDB ensemble and SAXS `.dat` file (same format as bioce `data/16_extrap_5concs-10p.dat`). + +## Resource hints + +- Limit conformers (e.g. ≤ 30) and use moderate Stan settings in the sidebar. +- First run compiles Stan models (several minutes). + +Documentation: [Serve user guide](https://serve.scilifelab.se/docs/) diff --git a/serve/environment.yml b/serve/environment.yml new file mode 100644 index 0000000..e92202a --- /dev/null +++ b/serve/environment.yml @@ -0,0 +1,16 @@ +name: bioce-serve +channels: + - conda-forge +dependencies: + - python=3.11 + - numpy + - matplotlib + - pandas + - biopython + - gsl + - clang + - cxx-compiler + - pip + - pip: + - "pystan>=3.9,<4" + - streamlit>=1.28 diff --git a/serve/examples/hdx_exp.txt b/serve/examples/hdx_exp.txt new file mode 100644 index 0000000..714e4c3 --- /dev/null +++ b/serve/examples/hdx_exp.txt @@ -0,0 +1,2 @@ +0.35 0.05 +0.62 0.06 diff --git a/serve/examples/hdx_peptides.txt b/serve/examples/hdx_peptides.txt new file mode 100644 index 0000000..ac06384 --- /dev/null +++ b/serve/examples/hdx_peptides.txt @@ -0,0 +1,3 @@ +# chain start end +A 2 10 +A 15 25 diff --git a/serve/examples/xl_restraints.txt b/serve/examples/xl_restraints.txt new file mode 100644 index 0000000..edb3167 --- /dev/null +++ b/serve/examples/xl_restraints.txt @@ -0,0 +1,3 @@ +# res_i res_j z d_max tau +12 48 1 30 3 +20 55 1 30 3 diff --git a/setup.py b/setup.py index 2c97c02..186ffd8 100644 --- a/setup.py +++ b/setup.py @@ -42,9 +42,12 @@ url='https://www.andrelab.org', ext_package='vbwSC', ext_modules=extensions, - py_modules=['psis','psisloo','stan_models','stan_utility', - 'statistics','variationalBayesian','fullBayesian', - 'prepareBayesian', 'prepareChemicalShifts'], + py_modules=['psis','psisloo','stan_models','stan_utility','stan_run', + 'bioce_statistics','variationalBayesian','fullBayesian', + 'fullBayesianMultimodal','multimodal_io','bioce_pipeline', + 'structure_observables', 'scilifelab_theme', + 'prepareBayesian', 'prepareChemicalShifts', + 'prepareHDX', 'prepareXLMS'], package_data={'bioce': ['bioce.yml']}, include_package_data=True, ) diff --git a/stan_models.py b/stan_models.py index f8c8919..ff29461 100644 --- a/stan_models.py +++ b/stan_models.py @@ -145,4 +145,225 @@ loglikes[i] = normal_lpdf(target_curve[i] | pred_curve[i], target_errors[i]); } } +""" + +# --- HDX-MS ensemble-averaged uptake (precomputed per conformer) --- +stan_code_HDX = """ +data { + int n_obs; + int n_structures; + vector[n_obs] target_hdx; + vector[n_obs] target_hdxerr; + matrix[n_obs, n_structures] sim_hdx; + vector[n_structures] priors; +} + +parameters { + simplex[n_structures] weights; +} + +model { + vector[n_obs] pred_hdx; + weights ~ dirichlet(priors); + pred_hdx = sim_hdx * weights; + target_hdx ~ normal(pred_hdx, target_hdxerr); +} +""" + +# --- XL-MS probabilistic distance restraints (soft satisfaction) --- +stan_code_XL = """ +data { + int n_xl; + int n_structures; + array[n_xl] int xl_obs; + vector[n_xl] d_max; + vector[n_xl] tau; + matrix[n_xl, n_structures] xl_distances; + vector[n_structures] priors; +} + +parameters { + simplex[n_structures] weights; +} + +model { + weights ~ dirichlet(priors); + for (l in 1:n_xl) { + xl_obs[l] ~ bernoulli( + dot_product(weights, inv_logit((d_max[l] - xl_distances[l]) / tau[l])) + ); + } +} +""" + +# False-positive mixture: P(z=1) = (1-pi_fp)*p_sat + pi_fp*(1-p_sat) +stan_code_XL_FP = """ +data { + int n_xl; + int n_structures; + array[n_xl] int xl_obs; + vector[n_xl] d_max; + vector[n_xl] tau; + matrix[n_xl, n_structures] xl_distances; + vector[n_structures] priors; +} + +parameters { + simplex[n_structures] weights; + real pi_fp; +} + +model { + real p_sat; + pi_fp ~ beta(1, 9); + weights ~ dirichlet(priors); + for (l in 1:n_xl) { + p_sat = dot_product(weights, inv_logit((d_max[l] - xl_distances[l]) / tau[l])); + xl_obs[l] ~ bernoulli((1 - pi_fp) * p_sat + pi_fp * (1 - p_sat)); + } +} +""" + +stan_code_SAXS_HDX = """ +data { + int n_measures; + int n_hdx; + int n_structures; + vector[n_measures] target_saxs; + vector[n_measures] target_saxserr; + matrix[n_measures, n_structures] sim_saxs; + vector[n_hdx] target_hdx; + vector[n_hdx] target_hdxerr; + matrix[n_hdx, n_structures] sim_hdx; + vector[n_structures] priors; +} + +parameters { + simplex[n_structures] weights; + real scale; +} + +model { + vector[n_measures] pred_saxs; + vector[n_hdx] pred_hdx; + weights ~ dirichlet(priors); + pred_saxs = sim_saxs * weights * scale; + pred_hdx = sim_hdx * weights; + target_saxs ~ normal(pred_saxs, target_saxserr); + target_hdx ~ normal(pred_hdx, target_hdxerr); +} +""" + +stan_code_SAXS_XL = """ +data { + int n_measures; + int n_xl; + int n_structures; + vector[n_measures] target_saxs; + vector[n_measures] target_saxserr; + matrix[n_measures, n_structures] sim_saxs; + array[n_xl] int xl_obs; + vector[n_xl] d_max; + vector[n_xl] tau; + matrix[n_xl, n_structures] xl_distances; + vector[n_structures] priors; +} + +parameters { + simplex[n_structures] weights; + real scale; +} + +model { + vector[n_measures] pred_saxs; + real p_sat; + weights ~ dirichlet(priors); + pred_saxs = sim_saxs * weights * scale; + target_saxs ~ normal(pred_saxs, target_saxserr); + for (l in 1:n_xl) { + p_sat = dot_product(weights, inv_logit((d_max[l] - xl_distances[l]) / tau[l])); + xl_obs[l] ~ bernoulli(p_sat); + } +} +""" + +stan_code_SAXS_HDX_XL = """ +data { + int n_measures; + int n_hdx; + int n_xl; + int n_structures; + vector[n_measures] target_saxs; + vector[n_measures] target_saxserr; + matrix[n_measures, n_structures] sim_saxs; + vector[n_hdx] target_hdx; + vector[n_hdx] target_hdxerr; + matrix[n_hdx, n_structures] sim_hdx; + array[n_xl] int xl_obs; + vector[n_xl] d_max; + vector[n_xl] tau; + matrix[n_xl, n_structures] xl_distances; + vector[n_structures] priors; +} + +parameters { + simplex[n_structures] weights; + real scale; +} + +model { + vector[n_measures] pred_saxs; + vector[n_hdx] pred_hdx; + real p_sat; + weights ~ dirichlet(priors); + pred_saxs = sim_saxs * weights * scale; + pred_hdx = sim_hdx * weights; + target_saxs ~ normal(pred_saxs, target_saxserr); + target_hdx ~ normal(pred_hdx, target_hdxerr); + for (l in 1:n_xl) { + p_sat = dot_product(weights, inv_logit((d_max[l] - xl_distances[l]) / tau[l])); + xl_obs[l] ~ bernoulli(p_sat); + } +} +""" + +stan_code_HDX_XL = """ +data { + int n_hdx; + int n_xl; + int n_structures; + vector[n_hdx] target_hdx; + vector[n_hdx] target_hdxerr; + matrix[n_hdx, n_structures] sim_hdx; + array[n_xl] int xl_obs; + vector[n_xl] d_max; + vector[n_xl] tau; + matrix[n_xl, n_structures] xl_distances; + vector[n_structures] priors; +} + +parameters { + simplex[n_structures] weights; +} + +model { + vector[n_hdx] pred_hdx; + real p_sat; + weights ~ dirichlet(priors); + pred_hdx = sim_hdx * weights; + target_hdx ~ normal(pred_hdx, target_hdxerr); + for (l in 1:n_xl) { + p_sat = dot_product(weights, inv_logit((d_max[l] - xl_distances[l]) / tau[l])); + xl_obs[l] ~ bernoulli(p_sat); + } +} +""" + +xl_posterior_predict = """ +generated quantities { + vector[n_xl] p_sat; + for (l in 1:n_xl) { + p_sat[l] = dot_product(weights, inv_logit((d_max[l] - xl_distances[l]) / tau[l])); + } +} """ \ No newline at end of file diff --git a/stan_run.py b/stan_run.py new file mode 100644 index 0000000..521716e --- /dev/null +++ b/stan_run.py @@ -0,0 +1,102 @@ +""" +PyStan 3 (import stan) helpers mirroring PyStan 2 sampling defaults. + +PyStan 2 used ``iter`` total iterations per chain with the first half as warmup +by default. CmdStan / PyStan 3 split ``num_warmup`` and ``num_samples`` explicitly. +""" +from __future__ import annotations + +from typing import Any, Optional, Tuple + + +def iterations_to_warmup_and_samples(iterations: int) -> Tuple[int, int]: + """Match PyStan 2 default split: ``warmup = iter // 2``, sampling = remainder.""" + iterations = int(iterations) + num_warmup = max(1, iterations // 2) + num_samples = max(1, iterations - num_warmup) + return num_warmup, num_samples + + +def build_and_sample( + model_code: str, + data: dict, + iterations: int, + chains: int, + njobs: Optional[int] = None, + random_seed: Optional[int] = None, + show_progress: bool = True, +) -> Any: + """ + Compile with ``stan.build`` and draw posterior samples. + + Parameters + ---------- + njobs + When the backend supports it, may map to parallel chain execution; ignored if + the installed httpstan rejects ``parallel_chains`` (falls back to defaults). + show_progress + When supported, passed to the sampler; ignored if the backend has no such field. + """ + import stan + + num_warmup, num_samples = iterations_to_warmup_and_samples(iterations) + posterior = stan.build( + model_code, data=data, random_seed=random_seed + ) + base_kwargs = { + "num_chains": int(chains), + "num_warmup": num_warmup, + "num_samples": num_samples, + } + optional = {} + if njobs is not None and int(njobs) > 0: + optional["parallel_chains"] = min(int(chains), int(njobs)) + if not show_progress: + optional["show_progress"] = False + + kwargs = {**base_kwargs, **optional} + try: + return posterior.sample(**kwargs) + except (TypeError, ValueError) as exc: + err = str(exc).lower() + # Newer httpstan / PyStan backends reject some CmdStan kwargs (ValueError) or + # raise TypeError for unexpected keywords. + if optional and ( + "unknown field" in err + or "unexpected keyword" in err + or "got an unexpected keyword" in err + ): + return posterior.sample(**base_kwargs) + raise + + +def plot_parameter_marginals(fit, par: str, filename: str, dpi: int = 300) -> None: + """Histogram of each scalar component of ``par`` (similar spirit to PyStan 2 ``fit.plot``).""" + import matplotlib.pyplot as plt + import numpy as np + + arr = np.asarray(fit[par]) + if arr.ndim == 1: + arr = arr.reshape(1, -1) + n_par = arr.shape[0] + ncols = min(5, max(1, n_par)) + nrows = (n_par + ncols - 1) // ncols + fig, axes = plt.subplots( + nrows, + ncols, + figsize=(3.2 * ncols, 2.4 * nrows), + squeeze=False, + ) + idx = 0 + for r in range(nrows): + for c in range(ncols): + ax = axes[r][c] + if idx < n_par: + ax.hist(arr[idx, :].ravel(), bins=40, density=True, color="steelblue") + ax.set_title("{}[{}]".format(par, idx + 1) if n_par > 1 else par) + else: + ax.axis("off") + idx += 1 + fig.subplots_adjust(wspace=0.35, hspace=0.45) + fig.savefig(filename, dpi=dpi) + plt.close(fig) diff --git a/stan_utility.py b/stan_utility.py index d728e4b..a6ac5b4 100644 --- a/stan_utility.py +++ b/stan_utility.py @@ -1,142 +1,128 @@ -import pystan -import pickle +""" +Diagnostics and helpers for PyStan 3 (``import stan``) fits. + +Sampler diagnostics use CmdStan-style names ending in ``__``, indexed like +``fit["divergent__"]``. See https://pystan.readthedocs.io/en/latest/upgrading.html +""" +from __future__ import print_function + import numpy +import stan + + +def _sampler_as_draws_by_chain(fit, name): + """Return array shape (num_saved_draws, num_chains) for a scalar sampler diagnostic.""" + x = numpy.asarray(fit[name]) + x = x.reshape(-1) + num_saved = fit._draws.shape[1] + num_chains = fit.num_chains + return x.reshape(num_saved, num_chains, order="F") + + def check_div(fit): - """Check transitions that ended with a divergence""" - sampler_params = fit.get_sampler_params(inc_warmup=False) - divergent = [x for y in sampler_params for x in y['divergent__']] - n = sum(divergent) - N = len(divergent) - print('{} of {} iterations ended with a divergence ({}%)'.format(n, N, - 100 * n / N)) + """Check transitions that ended with a divergence.""" + d = numpy.asarray(fit["divergent__"]).ravel() + n = int(d.sum()) + N = len(d) + print( + "{} of {} iterations ended with a divergence ({}%)".format( + n, N, 100 * n / N if N else 0 + ) + ) if n > 0: - print(' Try running with larger adapt_delta to remove the divergences') + print(" Try running with larger adapt_delta to remove the divergences") + -def check_treedepth(fit, max_depth = 10): - """Check transitions that ended prematurely due to maximum tree depth limit""" - sampler_params = fit.get_sampler_params(inc_warmup=False) - depths = [x for y in sampler_params for x in y['treedepth__']] +def check_treedepth(fit, max_depth=10): + """Check transitions that hit the maximum tree depth.""" + mat = _sampler_as_draws_by_chain(fit, "treedepth__") + depths = mat.ravel(order="C") n = sum(1 for x in depths if x == max_depth) N = len(depths) - print(('{} of {} iterations saturated the maximum tree depth of {}' - + ' ({}%)').format(n, N, max_depth, 100 * n / N)) + print( + ("{} of {} iterations saturated the maximum tree depth of {} ({}%)").format( + n, N, max_depth, 100 * n / N if N else 0 + ) + ) if n > 0: - print(' Run again with max_depth set to a larger value to avoid saturation') + print( + " Run again with max_depth passed through to sample() " + "to avoid saturation (CmdStan: max_depth keyword)." + ) + def check_energy(fit): - """Checks the energy Bayesian fraction of missing information (E-BFMI)""" - sampler_params = fit.get_sampler_params(inc_warmup=False) + """E-BFMI from energy__ within each chain.""" + mat = _sampler_as_draws_by_chain(fit, "energy__") no_warning = True - for chain_num, s in enumerate(sampler_params): - energies = s['energy__'] - numer = sum((energies[i] - energies[i - 1])**2 for i in range(1, len(energies))) / len(energies) + for chain_num in range(mat.shape[1]): + energies = mat[:, chain_num] + if len(energies) < 2: + continue + numer = ( + sum((energies[i] - energies[i - 1]) ** 2 for i in range(1, len(energies))) + / len(energies) + ) denom = numpy.var(energies) - if numer / denom < 0.2: - print('Chain {}: E-BFMI = {}'.format(chain_num, numer / denom)) + if denom == 0: + continue + ratio = numer / denom + if ratio < 0.2: + print("Chain {}: E-BFMI = {}".format(chain_num, ratio)) no_warning = False if no_warning: - print('E-BFMI indicated no pathological behavior') + print("E-BFMI indicated no pathological behavior") else: - print(' E-BFMI below 0.2 indicates you may need to reparameterize your model') + print( + " E-BFMI below 0.2 indicates you may need to reparameterize your model" + ) + def check_n_eff(fit): - """Checks the effective sample size per iteration""" - fit_summary = fit.summary(probs=[0.5]) - n_effs = [x[4] for x in fit_summary['summary']] - names = fit_summary['summary_rownames'] - n_iter = len(fit.extract()['lp__']) + """Not supported for PyStan 3 fits (no ``fit.summary``); use ArviZ on draws.""" + print( + "check_n_eff: skipped under PyStan 3 (no built-in summary table). " + "Use ArviZ or manual ESS from ``fit.to_frame()``." + ) - no_warning = True - for n_eff, name in zip(n_effs, names): - ratio = n_eff / n_iter - if (ratio < 0.001): - print('n_eff / iter for parameter {} is {}!'.format(name, ratio)) - print('E-BFMI below 0.2 indicates you may need to reparameterize your model') - no_warning = False - if no_warning: - print('n_eff / iter looks reasonable for all parameters') - else: - print(' n_eff / iter below 0.001 indicates that the effective sample size has likely been overestimated') def check_rhat(fit): - """Checks the potential scale reduction factors""" - from math import isnan - from math import isinf + """Not supported for PyStan 3 fits; use ArviZ on draws.""" + print( + "check_rhat: skipped under PyStan 3 (no built-in summary table). " + "Use ArviZ or manual R-hat from chain-split draws." + ) - fit_summary = fit.summary(probs=[0.5]) - rhats = [x[5] for x in fit_summary['summary']] - names = fit_summary['summary_rownames'] - - no_warning = True - for rhat, name in zip(rhats, names): - if (rhat > 1.1 or isnan(rhat) or isinf(rhat)): - print('Rhat for parameter {} is {}!'.format(name, rhat)) - no_warning = False - if no_warning: - print('Rhat looks reasonable for all parameters') - else: - print(' Rhat above 1.1 indicates that the chains very likely have not mixed') def check_all_diagnostics(fit): - """Checks all MCMC diagnostics""" + """Checks all MCMC diagnostics implemented for PyStan 3.""" check_n_eff(fit) check_rhat(fit) check_div(fit) check_treedepth(fit) check_energy(fit) -def _by_chain(unpermuted_extraction): - num_chains = len(unpermuted_extraction[0]) - result = [[] for _ in range(num_chains)] - for c in range(num_chains): - for i in range(len(unpermuted_extraction)): - result[c].append(unpermuted_extraction[i][c]) - return numpy.array(result) - -def _shaped_ordered_params(fit): - ef = fit.extract(permuted=False, inc_warmup=False) # flattened, unpermuted, by (iteration, chain) - ef = _by_chain(ef) - ef = ef.reshape(-1, len(ef[0][0])) - ef = ef[:, 0:len(fit.flatnames)] # drop lp__ - shaped = {} - idx = 0 - for dim, param_name in zip(fit.par_dims, fit.extract().keys()): - length = int(numpy.prod(dim)) - shaped[param_name] = ef[:,idx:idx + length] - shaped[param_name].reshape(*([-1] + dim)) - idx += length - return shaped def partition_div(fit): - """ Returns parameter arrays separated into divergent and non-divergent transitions""" - sampler_params = fit.get_sampler_params(inc_warmup=False) - div = numpy.concatenate([x['divergent__'] for x in sampler_params]).astype('int') - params = _shaped_ordered_params(fit) - nondiv_params = dict((key, params[key][div == 0]) for key in params) - div_params = dict((key, params[key][div == 1]) for key in params) - return nondiv_params, div_params + """Not ported for PyStan 3 (requires ordered unconstrained draws per transition).""" + raise NotImplementedError( + "partition_div is not implemented for PyStan 3 fits. " + "Use ArviZ or CmdStanPy if you need divergent/non-divergent draw splits." + ) + -def compile_model(filename, model_name=None, **kwargs): - """This will automatically cache models - great if you're just running a - script on the command line. +def compile_model(filename, model_name=None, data=None, random_seed=None, **kwargs): + """Compile a Stan program with ``stan.build`` (PyStan 3). - See http://pystan.readthedocs.io/en/latest/avoiding_recompilation.html""" - from hashlib import md5 + Unlike PyStan 2, a program with a nonempty ``data`` block must receive matching + ``data`` here. Optional disk caching of compiled models is not supported. + See https://pystan.readthedocs.io/en/latest/upgrading.html + """ with open(filename) as f: model_code = f.read() - code_hash = md5(model_code.encode('ascii')).hexdigest() - if model_name is None: - cache_fn = 'cached-model-{}.pkl'.format(code_hash) - else: - cache_fn = 'cached-{}-{}.pkl'.format(model_name, code_hash) - try: - sm = pickle.load(open(cache_fn, 'rb')) - except: - sm = pystan.StanModel(model_code=model_code) - with open(cache_fn, 'wb') as f: - pickle.dump(sm, f) - else: - print("Using cached StanModel") - return sm + payload = {} if data is None else data + return stan.build( + model_code, data=payload, random_seed=random_seed + ) diff --git a/structure_observables.py b/structure_observables.py new file mode 100644 index 0000000..cd47ec4 --- /dev/null +++ b/structure_observables.py @@ -0,0 +1,351 @@ +""" +Compute SAXS, HDX, and XL-MS observables from PDB conformers. + +SAXS: FoXS (preferred) or Pepsi-SAXS if available on PATH. +HDX: peptide mean relative SASA (Shrake–Rupley) as a protection proxy. +XL: Cα–Cα distances for residue pairs in the restraint list. +""" +from __future__ import annotations + +import os +import shlex +import shutil +import zipfile +from pathlib import Path +from subprocess import PIPE, Popen +from typing import List, Optional, Sequence, Tuple + +import numpy as np + +try: + from Bio.PDB import PDBParser, ShrakeRupley + from Bio.PDB.Polypeptide import is_aa + _HAS_BIOPYTHON = True +except ImportError: + _HAS_BIOPYTHON = False + +PDBParser_instance = PDBParser(QUIET=True) if _HAS_BIOPYTHON else None + + +def _foxs_on_path() -> bool: + return shutil.which("foxs") is not None + + +def _pepsi_on_path() -> bool: + return shutil.which("Pepsi-SAXS") is not None + + +def extract_pdbs(uploaded, dest_dir: Path) -> List[Path]: + """Save uploaded Streamlit files or a zip archive into dest_dir; return .pdb paths.""" + dest_dir.mkdir(parents=True, exist_ok=True) + paths: List[Path] = [] + for item in uploaded: + name = getattr(item, "name", str(item)) + if name.lower().endswith(".zip"): + with zipfile.ZipFile(item) as zf: + zf.extractall(dest_dir) + else: + data = item.getvalue() if hasattr(item, "getvalue") else Path(item).read_bytes() + out = dest_dir / Path(name).name + out.write_bytes(data) + for pdb in sorted(dest_dir.rglob("*.pdb")): + paths.append(pdb) + for pdb in sorted(dest_dir.rglob("*.PDB")): + if pdb not in paths: + paths.append(pdb) + if not paths: + raise ValueError("No PDB files found in upload.") + return paths + + +def _run_foxs(pdb_path: Path, exp_path: Path, work_dir: Path) -> List[float]: + cmd = "foxs {} {}".format(pdb_path, exp_path) + proc = Popen(shlex.split(cmd), stdout=PIPE, cwd=str(work_dir), shell=False) + proc.wait() + out_name = pdb_path.stem + "_" + exp_path.name + foxs_file = work_dir / out_name + if not foxs_file.exists(): + foxs_file = pdb_path.parent / out_name + if not foxs_file.exists(): + raise RuntimeError("FoXS did not produce {}".format(out_name)) + lines = foxs_file.read_text().splitlines()[1:] + norm_c = 1e8 + scaling_c = norm_c * float(lines[0].split(",")[1][12:]) + intensities = [] + for line in lines[2:]: + intensities.append(float(line.split()[2]) / scaling_c) + return intensities + + +def _run_pepsi(pdb_path: Path, exp_path: Path) -> np.ndarray: + cmd = ( + "Pepsi-SAXS --dp_min 0 --dp_max 0 --dp_N 1 --r0_min_factor 1 " + "--r0_max_factor 1 --r0_N 1 {} {}".format(pdb_path, exp_path) + ) + proc = Popen(shlex.split(cmd), stdout=PIPE, shell=False) + proc.wait() + fit_file = Path(str(pdb_path)[:-4] + "-" + str(exp_path)[:-3] + "fit") + if not fit_file.exists(): + fit_file = pdb_path.parent / fit_file.name + data = np.genfromtxt(fit_file, skip_header=6) + intensity = data[:, 3] + log_path = Path(str(pdb_path)[:-4] + "-" + str(exp_path)[:-3] + "log") + log_line = log_path.read_text().splitlines()[62] + if not log_line.startswith("Scaling"): + raise RuntimeError("Pepsi-SAXS scaling line not found in log") + scale = float(log_line.split(":")[1]) + return intensity / scale + + +def generate_saxs_matrix( + pdb_paths: Sequence[Path], + experimental_path: Path, + out_dir: Path, +) -> Tuple[Path, Path, Path, np.ndarray]: + """ + Write SimulatedIntensities.txt, weights.txt, structures.txt; return paths + matrix. + """ + out_dir.mkdir(parents=True, exist_ok=True) + exp_copy = out_dir / experimental_path.name + shutil.copy(experimental_path, exp_copy) + + if _foxs_on_path(): + backend = "FoXS" + columns = [] + for pdb in pdb_paths: + columns.append(_run_foxs(pdb, exp_copy, pdb.parent)) + elif _pepsi_on_path(): + backend = "Pepsi-SAXS" + columns = [_run_pepsi(pdb, exp_copy) for pdb in pdb_paths] + else: + raise RuntimeError( + "Neither FoXS nor Pepsi-SAXS found on PATH. " + "Install ATSAS FoXS or Pepsi-SAXS for SAXS simulation." + ) + + sim = np.column_stack(columns) + names = [p.stem for p in pdb_paths] + n = len(names) + priors = np.full(n, 1.0 / n) + + sim_path = out_dir / "SimulatedIntensities.txt" + weights_path = out_dir / "weights.txt" + names_path = out_dir / "structures.txt" + np.savetxt(sim_path, sim, fmt="%.8g") + np.savetxt(weights_path, priors, fmt="%.8g") + with open(names_path, "w") as handle: + handle.write(" ".join(names) + "\n") + + (out_dir / "saxs_backend.txt").write_text(backend) + return sim_path, weights_path, names_path, sim + + +def _get_ca_coord(structure, chain_id: str, resseq: int): + model = structure[0] + chain = model[chain_id] if chain_id in [c.id for c in model] else model.child_list[0] + res = chain[resseq] + for atom_name in ("CA", "CB", "N"): + if atom_name in res: + return res[atom_name].coord + raise KeyError("No CA/CB/N for {}:{}".format(chain_id, resseq)) + + +def ca_distance(pdb_path: Path, chain_i: str, res_i: int, chain_j: str, res_j: int) -> float: + if not _HAS_BIOPYTHON: + raise RuntimeError("BioPython is required for XL/HDX from PDB (pip install biopython).") + structure = PDBParser_instance.get_structure("bioce", str(pdb_path)) + ci = _get_ca_coord(structure, chain_i, res_i) + cj = _get_ca_coord(structure, chain_j, res_j) + return float(np.linalg.norm(ci - cj)) + + +def parse_xl_restraints(lines: List[str]) -> List[dict]: + """ + Parse crosslink lines. Supported formats per line: + res_i res_j z [d_max] [tau] + chain_i res_i chain_j res_j z [d_max] [tau] + """ + records = [] + for raw in lines: + line = raw.strip() + if not line or line.startswith("#"): + continue + parts = line.split() + if len(parts) < 3: + continue + if len(parts) in (3, 4, 5): + res_i, res_j, z = int(parts[0]), int(parts[1]), int(float(parts[2])) + d_max = float(parts[3]) if len(parts) > 3 else 30.0 + tau = float(parts[4]) if len(parts) > 4 else 3.0 + records.append({ + "chain_i": "A", "res_i": res_i, "chain_j": "A", "res_j": res_j, + "z": z, "d_max": d_max, "tau": tau, + }) + else: + chain_i, res_i = parts[0], int(parts[1]) + chain_j, res_j = parts[2], int(parts[3]) + z = int(float(parts[4])) + d_max = float(parts[5]) if len(parts) > 5 else 30.0 + tau = float(parts[6]) if len(parts) > 6 else 3.0 + records.append({ + "chain_i": chain_i, "res_i": res_i, "chain_j": chain_j, "res_j": res_j, + "z": z, "d_max": d_max, "tau": tau, + }) + if not records: + raise ValueError("No XL restraints parsed.") + return records + + +def write_xl_restraints_std(records: List[dict], path: Path) -> None: + with open(path, "w") as handle: + for r in records: + handle.write( + "{} {} {} {} {} {} {}\n".format( + r["chain_i"], r["res_i"], r["chain_j"], r["res_j"], + r["z"], r["d_max"], r["tau"], + ) + ) + + +def generate_xl_distance_matrix( + pdb_paths: Sequence[Path], + records: List[dict], + out_dir: Path, +) -> Tuple[Path, Path, np.ndarray]: + out_dir.mkdir(parents=True, exist_ok=True) + n_xl = len(records) + n_struct = len(pdb_paths) + dist = np.zeros((n_xl, n_struct)) + for j, pdb in enumerate(pdb_paths): + for i, rec in enumerate(records): + try: + dist[i, j] = ca_distance( + pdb, rec["chain_i"], rec["res_i"], rec["chain_j"], rec["res_j"] + ) + except (KeyError, IndexError) as exc: + dist[i, j] = 999.0 + xl_path = out_dir / "xl_restraints.dat" + write_xl_restraints_for_stan(records, xl_path) + write_xl_restraints_std(records, out_dir / "xl_restraints_full.dat") + sim_path = out_dir / "SimulatedXLdistances.txt" + np.savetxt(sim_path, dist, fmt="%.4f") + return xl_path, sim_path, dist + + +def parse_hdx_peptides(lines: List[str]) -> List[dict]: + """Each line: chain start_res end_res (1-based inclusive).""" + peptides = [] + for raw in lines: + line = raw.strip() + if not line or line.startswith("#"): + continue + parts = line.split() + if len(parts) < 3: + continue + peptides.append({ + "chain": parts[0], + "start": int(parts[1]), + "end": int(parts[2]), + }) + if not peptides: + raise ValueError("No HDX peptides defined.") + return peptides + + +def parse_hdx_experimental( + lines: List[str], n_peptides: int +) -> Tuple[np.ndarray, np.ndarray, Optional[np.ndarray]]: + """ + Formats: + uptake sigma (n_peptides lines, single time point) + pep_idx time uptake sigma (multiple time points; pep_idx 0-based) + Returns (uptake, sigma, pep_indices) where pep_indices is None for single-time. + """ + rows = [] + for raw in lines: + line = raw.strip() + if not line or line.startswith("#"): + continue + rows.append([float(x) for x in line.split()]) + if not rows: + raise ValueError("No HDX experimental data.") + arr = np.asarray(rows) + if arr.shape[1] == 2: + if arr.shape[0] != n_peptides: + raise ValueError( + "Expected {} peptide uptake rows, got {}".format(n_peptides, arr.shape[0]) + ) + return arr[:, 0], arr[:, 1], None + if arr.shape[1] >= 4: + order = np.lexsort((arr[:, 1], arr[:, 0])) + arr = arr[order] + pep_idx = arr[:, 0].astype(int) + return arr[:, 2], arr[:, 3], pep_idx + raise ValueError("HDX experimental file: use 'uptake sigma' or 'pep_idx time uptake sigma'.") + + +def peptide_mean_relative_sasa(pdb_path: Path, chain_id: str, start: int, end: int) -> float: + if not _HAS_BIOPYTHON: + raise RuntimeError("BioPython is required for HDX from PDB.") + structure = PDBParser_instance.get_structure("bioce", str(pdb_path)) + sr = ShrakeRupley() + sr.compute(structure, level="R") + model = structure[0] + chain = model[chain_id] if chain_id in [c.id for c in model] else model.child_list[0] + values = [] + for res in chain: + if not is_aa(res, standard=True): + continue + if res.id[1] < start or res.id[1] > end: + continue + sasa = res.sasa if hasattr(res, "sasa") else 0.0 + values.append(sasa) + if not values: + return 0.5 + mean_sasa = float(np.mean(values)) + return min(1.0, mean_sasa / 200.0) + + +def generate_hdx_matrix( + pdb_paths: Sequence[Path], + peptides: List[dict], + out_dir: Path, + pep_indices: Optional[np.ndarray] = None, +) -> Tuple[Path, np.ndarray]: + """Build simulated HDX matrix; rows follow pep_indices or one row per peptide.""" + out_dir.mkdir(parents=True, exist_ok=True) + if pep_indices is None: + obs_pep = list(range(len(peptides))) + else: + obs_pep = [int(i) for i in pep_indices] + n_obs = len(obs_pep) + n_struct = len(pdb_paths) + sim = np.zeros((n_obs, n_struct)) + cache = {} + for j, pdb in enumerate(pdb_paths): + for row_i, pidx in enumerate(obs_pep): + key = (pidx, j) + if key not in cache: + pep = peptides[pidx] + cache[key] = peptide_mean_relative_sasa( + pdb, pep["chain"], pep["start"], pep["end"] + ) + sim[row_i, j] = cache[key] + sim_path = out_dir / "SimulatedHDX.txt" + np.savetxt(sim_path, sim, fmt="%.6f") + return sim_path, sim + + +def write_xl_restraints_for_stan(records: List[dict], path: Path) -> None: + """Stan reader format: res_i res_j z [d_max] [tau].""" + with open(path, "w") as handle: + for r in records: + handle.write( + "{} {} {} {} {}\n".format( + r["res_i"], r["res_j"], r["z"], r["d_max"], r["tau"] + ) + ) + + +def write_hdx_experimental(uptake: np.ndarray, sigma: np.ndarray, path: Path) -> None: + np.savetxt(path, np.column_stack([uptake, sigma]), fmt="%.6f") diff --git a/tests/fixtures/two_chain_mini_A.pdb b/tests/fixtures/two_chain_mini_A.pdb new file mode 100644 index 0000000..205353c --- /dev/null +++ b/tests/fixtures/two_chain_mini_A.pdb @@ -0,0 +1,24 @@ +ATOM 1 N ALA A 2 0.000 0.000 0.000 1.00 0.00 N +ATOM 2 CA ALA A 2 1.458 0.000 0.000 1.00 0.00 C +ATOM 3 C ALA A 2 2.009 1.420 0.000 1.00 0.00 C +ATOM 4 O ALA A 2 1.251 2.390 0.000 1.00 0.00 O +ATOM 5 CB ALA A 2 1.958 -0.770 1.210 1.00 0.00 C +ATOM 6 N GLY A 3 3.300 1.520 0.000 1.00 0.00 N +ATOM 7 CA GLY A 3 4.000 2.800 0.000 1.00 0.00 C +ATOM 8 C GLY A 3 5.500 2.600 0.000 1.00 0.00 C +ATOM 9 O GLY A 3 6.100 1.550 0.000 1.00 0.00 O +ATOM 10 N ALA A 10 10.000 0.000 0.000 1.00 0.00 N +ATOM 11 CA ALA A 10 11.458 0.000 0.000 1.00 0.00 C +ATOM 12 C ALA A 10 12.009 1.420 0.000 1.00 0.00 C +ATOM 13 O ALA A 10 11.251 2.390 0.000 1.00 0.00 O +ATOM 14 CB ALA A 10 11.958 -0.770 1.210 1.00 0.00 C +ATOM 15 N GLY A 15 20.000 0.000 0.000 1.00 0.00 N +ATOM 16 CA GLY A 15 21.458 0.000 0.000 1.00 0.00 C +ATOM 17 C GLY A 15 22.009 1.420 0.000 1.00 0.00 C +ATOM 18 O GLY A 15 21.251 2.390 0.000 1.00 0.00 O +ATOM 19 N ALA A 25 40.000 0.000 0.000 1.00 0.00 N +ATOM 20 CA ALA A 25 41.458 0.000 0.000 1.00 0.00 C +ATOM 21 C ALA A 25 42.009 1.420 0.000 1.00 0.00 C +ATOM 22 O ALA A 25 41.251 2.390 0.000 1.00 0.00 O +ATOM 23 CB ALA A 25 41.958 -0.770 1.210 1.00 0.00 C +END diff --git a/tests/fixtures/two_chain_mini_B.pdb b/tests/fixtures/two_chain_mini_B.pdb new file mode 100644 index 0000000..c95ffac --- /dev/null +++ b/tests/fixtures/two_chain_mini_B.pdb @@ -0,0 +1,24 @@ +ATOM 1 N ALA A 2 0.000 0.000 0.000 1.00 0.00 N +ATOM 2 CA ALA A 2 1.458 0.000 0.000 1.00 0.00 C +ATOM 3 C ALA A 2 2.009 1.420 0.000 1.00 0.00 C +ATOM 4 O ALA A 2 1.251 2.390 0.000 1.00 0.00 O +ATOM 5 CB ALA A 2 1.958 -0.770 1.210 1.00 0.00 C +ATOM 6 N GLY A 3 3.300 1.520 0.000 1.00 0.00 N +ATOM 7 CA GLY A 3 4.000 2.800 0.000 1.00 0.00 C +ATOM 8 C GLY A 3 5.500 2.600 0.000 1.00 0.00 C +ATOM 9 O GLY A 3 6.100 1.550 0.000 1.00 0.00 O +ATOM 10 N ALA A 10 12.000 0.000 0.000 1.00 0.00 N +ATOM 11 CA ALA A 10 13.458 0.000 0.000 1.00 0.00 C +ATOM 12 C ALA A 10 14.009 1.420 0.000 1.00 0.00 C +ATOM 13 O ALA A 10 13.251 2.390 0.000 1.00 0.00 O +ATOM 14 CB ALA A 10 13.958 -0.770 1.210 1.00 0.00 C +ATOM 15 N GLY A 15 18.000 0.000 0.000 1.00 0.00 N +ATOM 16 CA GLY A 15 19.458 0.000 0.000 1.00 0.00 C +ATOM 17 C GLY A 15 20.009 1.420 0.000 1.00 0.00 C +ATOM 18 O GLY A 15 19.251 2.390 0.000 1.00 0.00 O +ATOM 19 N ALA A 25 35.000 0.000 0.000 1.00 0.00 N +ATOM 20 CA ALA A 25 36.458 0.000 0.000 1.00 0.00 C +ATOM 21 C ALA A 25 37.009 1.420 0.000 1.00 0.00 C +ATOM 22 O ALA A 25 36.251 2.390 0.000 1.00 0.00 O +ATOM 23 CB ALA A 25 36.958 -0.770 1.210 1.00 0.00 C +END diff --git a/tests/test_multimodal_stan.py b/tests/test_multimodal_stan.py new file mode 100644 index 0000000..2fca66e --- /dev/null +++ b/tests/test_multimodal_stan.py @@ -0,0 +1,140 @@ +""" +Compile/sample tests for HDX and XL-MS Stan models (skip if pystan missing). +""" +from __future__ import annotations + +import tempfile +from pathlib import Path + +import numpy as np +import pytest + + +@pytest.fixture +def priors(): + return np.array([1.0, 1.0, 1.0]) + + +@pytest.fixture +def tmp_dir(): + with tempfile.TemporaryDirectory() as d: + yield Path(d) + + +def test_read_xl_restraints_defaults(): + from multimodal_io import read_xl_restraints, DEFAULT_XL_DMAX, DEFAULT_XL_TAU + + with tempfile.NamedTemporaryFile(mode="w", suffix=".dat", delete=False) as f: + f.write("10 20 1\n") + f.write("30 40 0\n") + path = f.name + obs, dmax, tau = read_xl_restraints(path) + assert list(obs) == [1, 0] + assert dmax[0] == DEFAULT_XL_DMAX + assert tau[1] == DEFAULT_XL_TAU + + +def test_stan_hdx_only(priors): + pytest.importorskip("stan") + from stan_run import build_and_sample + from stan_models import stan_code_HDX + + n_obs, n_structures = 6, 3 + rng = np.random.default_rng(0) + sim = rng.uniform(0, 1, (n_obs, n_structures)) + w_true = np.array([0.5, 0.3, 0.2]) + target = sim @ w_true + stan_dat = { + "n_obs": n_obs, + "n_structures": n_structures, + "target_hdx": target, + "target_hdxerr": np.full(n_obs, 0.05), + "sim_hdx": sim, + "priors": priors, + } + fit = build_and_sample( + stan_code_HDX, stan_dat, iterations=40, chains=2, njobs=2, show_progress=False + ) + assert "weights" in fit + + +def test_stan_xl_only(priors): + pytest.importorskip("stan") + from stan_run import build_and_sample + from stan_models import stan_code_XL + + n_xl, n_structures = 4, 3 + distances = np.array([ + [10.0, 25.0, 28.0], + [12.0, 11.0, 35.0], + [30.0, 32.0, 31.0], + [15.0, 16.0, 14.0], + ]) + stan_dat = { + "n_xl": n_xl, + "n_structures": n_structures, + "xl_obs": np.array([1, 1, 0, 1], dtype=np.int32), + "d_max": np.full(n_xl, 30.0), + "tau": np.full(n_xl, 3.0), + "xl_distances": distances, + "priors": priors, + } + fit = build_and_sample( + stan_code_XL, stan_dat, iterations=40, chains=2, njobs=2, show_progress=False + ) + assert "weights" in fit + + +def test_stan_saxs_hdx_xl(priors): + pytest.importorskip("stan") + from stan_run import build_and_sample + from stan_models import stan_code_SAXS_HDX_XL, xl_posterior_predict + + n_measures, n_hdx, n_xl, n_structures = 5, 4, 3, 3 + rng = np.random.default_rng(1) + sim_saxs = rng.uniform(0.01, 0.05, (n_measures, n_structures)) + sim_hdx = rng.uniform(0, 1, (n_hdx, n_structures)) + xl_distances = rng.uniform(8, 35, (n_xl, n_structures)) + stan_dat = { + "n_measures": n_measures, + "n_hdx": n_hdx, + "n_xl": n_xl, + "n_structures": n_structures, + "target_saxs": sim_saxs.mean(axis=1), + "target_saxserr": np.full(n_measures, 0.001), + "sim_saxs": sim_saxs, + "target_hdx": sim_hdx.mean(axis=1), + "target_hdxerr": np.full(n_hdx, 0.05), + "sim_hdx": sim_hdx, + "xl_obs": np.array([1, 1, 0], dtype=np.int32), + "d_max": np.full(n_xl, 30.0), + "tau": np.full(n_xl, 3.0), + "xl_distances": xl_distances, + "priors": priors, + } + fit = build_and_sample( + stan_code_SAXS_HDX_XL + xl_posterior_predict, + stan_dat, + iterations=40, + chains=2, + njobs=2, + show_progress=False, + ) + assert "weights" in fit + assert "p_sat" in fit + + +def test_fullBayesianMultimodal_hdx_cli(tmp_dir, priors): + pytest.importorskip("stan") + from fullBayesianMultimodal import execute_stan_hdx + + n_obs, n_structures = 5, 3 + sim = np.ones((n_obs, n_structures)) * np.linspace(0.2, 0.8, n_structures) + exp_path = tmp_dir / "hdx_exp.dat" + sim_path = tmp_dir / "SimulatedHDX.txt" + np.savetxt(exp_path, np.column_stack([sim.mean(axis=1), np.full(n_obs, 0.1)])) + np.savetxt(sim_path, sim) + fit = execute_stan_hdx( + str(exp_path), str(sim_path), priors, iterations=32, chains=2, njobs=2 + ) + assert "weights" in fit diff --git a/tests/test_pipeline.py b/tests/test_pipeline.py new file mode 100644 index 0000000..554869d --- /dev/null +++ b/tests/test_pipeline.py @@ -0,0 +1,71 @@ +"""Pipeline tests (HDX+XL from PDB; SAXS if FoXS/Pepsi available).""" +from __future__ import annotations + +from pathlib import Path + +import numpy as np +import pytest + +FIXTURES = Path(__file__).parent / "fixtures" + + +class _FakeUpload: + def __init__(self, path: Path): + self.name = path.name + self._path = path + + def getvalue(self): + return self._path.read_bytes() + + +@pytest.fixture +def mini_pdbs(): + return [ + _FakeUpload(FIXTURES / "two_chain_mini_A.pdb"), + _FakeUpload(FIXTURES / "two_chain_mini_B.pdb"), + ] + + +def test_observables_xl_hdx(mini_pdbs, tmp_path): + pytest.importorskip("Bio.PDB") + from structure_observables import ( + extract_pdbs, + generate_hdx_matrix, + generate_xl_distance_matrix, + parse_hdx_peptides, + parse_xl_restraints, + ) + + pdb_dir = tmp_path / "pdbs" + paths = extract_pdbs(mini_pdbs, pdb_dir) + xl_lines = ["2 25 1 30 3", "10 15 0 30 3"] + records = parse_xl_restraints(xl_lines) + _, sim_xl, dist = generate_xl_distance_matrix(paths, records, tmp_path / "xl") + assert dist.shape == (2, 2) + assert dist[0, 0] < dist[0, 1] + + peptides = parse_hdx_peptides(["A 2 10", "A 15 25"]) + _, sim_h = generate_hdx_matrix(paths, peptides, tmp_path / "hdx") + assert sim_h.shape == (2, 2) + + +def test_pipeline_hdx_xl_inference(mini_pdbs, tmp_path): + pytest.importorskip("stan") + pytest.importorskip("Bio.PDB") + from bioce_pipeline import PipelineConfig, run_full_pipeline + + config = PipelineConfig( + work_dir=tmp_path / "run", + use_saxs=False, + use_hdx=True, + use_xl=True, + iterations=60, + chains=2, + njobs=2, + hdx_peptide_lines=["A 2 10", "A 15 25"], + hdx_exp_lines=["0.3 0.05", "0.55 0.05"], + xl_restraint_lines=["2 25 1", "10 15 0"], + ) + result = run_full_pipeline(mini_pdbs, config) + assert result.bayesian_weights.shape[0] == 2 + assert np.isclose(result.bayesian_weights.sum(), 1.0, atol=0.02) diff --git a/tests/test_pystan_smoke.py b/tests/test_pystan_smoke.py new file mode 100644 index 0000000..d23ff04 --- /dev/null +++ b/tests/test_pystan_smoke.py @@ -0,0 +1,78 @@ +""" +Smoke tests for PyStan 3 (``import stan``) and the bioce ``stan_run`` helpers. + +Requires: pip install pystan pytest +Stan-dependent tests skip automatically if ``stan`` is not installed. +""" +from __future__ import annotations + +import pytest + + +MINIMAL_STAN = """ +data { + int n; + array[n] real y; +} +parameters { + real mu; +} +model { + y ~ normal(mu, 1); +} +""" + + +def test_iterations_to_warmup_and_samples_matches_pystan2_style(): + from stan_run import iterations_to_warmup_and_samples + + w, s = iterations_to_warmup_and_samples(2000) + assert w == 1000 + assert s == 1000 + assert w + s == 2000 + + w2, s2 = iterations_to_warmup_and_samples(3) + assert w2 >= 1 and s2 >= 1 + assert w2 + s2 == 3 + + +def test_build_and_sample_minimal_model(): + """Compile and draw a few samples (exercises httpstan/CmdStan bridge).""" + pytest.importorskip("stan", reason="pystan not installed (pip install pystan)") + from stan_run import build_and_sample + + data = {"n": 4, "y": [0.0, 0.1, -0.05, 0.02]} + fit = build_and_sample( + MINIMAL_STAN, + data, + iterations=32, + chains=2, + njobs=2, + random_seed=12345, + show_progress=False, + ) + assert "mu" in fit + mu = fit["mu"] + assert mu.size >= 1 + # PyStan 3: scalar params often shaped (1, num_draws) + assert mu.ravel().shape[0] >= 4 + + +def test_stan_utility_diagnostics_do_not_crash(): + pytest.importorskip("stan", reason="pystan not installed (pip install pystan)") + from stan_run import build_and_sample + import stan_utility + + data = {"n": 3, "y": [0.0, 0.0, 0.0]} + fit = build_and_sample( + MINIMAL_STAN, + data, + iterations=24, + chains=2, + njobs=2, + random_seed=99, + show_progress=False, + ) + stan_utility.check_div(fit) + stan_utility.check_treedepth(fit, max_depth=10) + stan_utility.check_energy(fit)