Follow-ups for #984: name every sorted result file, credit srivarshni7 (waiting on sign-off)

CONTRIBUTORS.md

@@ -20,5 +20,6 @@ handle are listed by handle, in the same sequence.
 - Richard G. Posner
 - psamuelvijay
 - Herbert M. Sauro
+- srivarshni7
 - Ryan Suderman
 - Brandon R. Thomas

docs/config_keys.rst

@@ -2786,9 +2786,13 @@ For Adaptive MCMC
 
   Prefer :ref:`start_point <start_point>` for new configurations. It is read by every ``job_type``, it is
   matched by **name** rather than by position — ``starting_params`` is positional against declaration
-  order, while every result file PyBNF writes is alphabetical, so round-tripping a result row back into
-  it silently permutes the values — and it is validated against the declared bounds. Note also that
-  ``continue_run = 1`` overrides ``starting_params`` entirely.
+  order, while ``samples.txt``, ``sorted_params*.txt`` and ``bootstrapped_parameter_sets.txt`` sort their
+  parameter columns by name, case-sensitively (``kcat``, ``Km``, ``alpha`` are written ``Km``, ``alpha``,
+  ``kcat``), so round-tripping a result row back into it silently permutes the values — and it is
+  validated against the declared bounds. It is not a drop-in replacement for a sampler, though:
+  ``start_point`` pins only the first chain and the others start from independent draws, while
+  ``starting_params`` sets every chain. Note also that ``continue_run = 1`` overrides
+  ``starting_params`` entirely.
 
   Default: None
 

tests/test_result_file_parameter_order.py

@@ -0,0 +1,244 @@
+"""Which order each result file lists the free parameters in (lanl/PyBNF#932).
+
+The ``starting_params`` entry in docs/config_keys.rst warns that ``starting_params`` is
+positional against declaration order, while ``samples.txt``, ``sorted_params*.txt`` and
+``bootstrapped_parameter_sets.txt`` sort their parameter columns by name, case-sensitively, so a row
+copied out of one of those files into ``starting_params`` lands its values on the wrong
+parameters. docs/troubleshooting.rst says ``start_point.txt`` is in declaration order. The
+same entry says ``continue_run = 1`` overrides ``starting_params``. These tests pin each of
+those statements against what PyBNF really writes, so the documentation cannot go stale
+silently when a writer changes.
+
+The parameters are declared ``kcat``, ``Km``, ``alpha``. That order is neither sorted nor
+case-insensitive alphabetical (``alpha``, ``kcat``, ``Km``), and the files hold
+``Km``, ``alpha``, ``kcat``: the sort is Python's case-sensitive string order, which puts
+every uppercase letter before every lowercase one. The expected orders below are written
+out by hand rather than computed with ``sorted()``, so they are an oracle independent of
+the writers.
+
+Each parameter's box is disjoint from the other two, so a value written under the wrong
+header lands outside the box its header names. That is how the row checks below tell a
+correctly labelled column from a permuted one without knowing what the fit converged to.
+"""
+import os
+import types
+
+import numpy as np
+
+from . import integration_harness as H
+from .context import algorithms, config
+
+DECLARED = ['kcat', 'Km', 'alpha']
+SORTED_IN_FILES = ['Km', 'alpha', 'kcat']   # case-sensitive: uppercase first
+BOXES = {'kcat': (1.0, 9.0), 'Km': (50.0, 90.0), 'alpha': (0.1, 0.9)}
+TRUE = {'kcat': 7.0, 'Km': 70.0, 'alpha': 0.3}
+
+
+def test_the_fixture_distinguishes_the_three_orders():
+    """Declaration order, case-insensitive alphabetical order and the order the files
+    use must all differ, or the tests below could not tell them apart."""
+    assert sorted(DECLARED, key=str.casefold) == ['alpha', 'kcat', 'Km']
+    assert len({tuple(DECLARED), tuple(SORTED_IN_FILES), ('alpha', 'kcat', 'Km')}) == 3
+    for name, (lo, hi) in BOXES.items():
+        for other, (lo2, hi2) in BOXES.items():
+            if other != name:
+                assert hi < lo2 or hi2 < lo, (name, other)
+        assert lo < TRUE[name] < hi
+
+
+def _vars_conf(declare):
+    return ''.join('uniform_var = %s %g %g\n' % (n, *BOXES[n]) for n in declare)
+
+
+def _read_table(path, n_lead):
+    """Header names after the ``n_lead`` leading non-parameter columns, and the rows as
+    ``{name: value}`` dicts keyed by those header names."""
+    with open(path) as f:
+        lines = [ln for ln in f.read().splitlines() if ln.strip()]
+    header = lines[0].lstrip('#').split()
+    names = header[n_lead:]
+    rows = [dict(zip(names, map(float, ln.split()[n_lead:]))) for ln in lines[1:]]
+    return names, rows
+
+
+def _assert_each_value_in_the_box_its_header_names(rows, where):
+    assert rows, where
+    for row in rows:
+        for name, value in row.items():
+            lo, hi = BOXES[name]
+            assert lo <= value <= hi, (where, name, value)
+
+
+def _declared_rows(path):
+    """``(name, first value)`` for each row of a one-row-per-parameter file
+    (start_point.txt, credible*.txt), skipping comments and the ``starts`` lines."""
+    out = []
+    with open(path) as f:
+        for ln in f.read().splitlines():
+            parts = ln.split('\t')
+            if ln.startswith('#') or not ln.strip() or parts[0] in ('starts', 'starts_pinned'):
+                continue
+            out.append((parts[0], float(parts[1])))
+    return out
+
+
+# --------------------------------------------------------------------------- #
+# A sampler: samples.txt and sorted_params are sorted, the per-parameter files are
+# in declaration order, and starting_params is applied by declaration order.
+# --------------------------------------------------------------------------- #
+
+_SAMPLER_CONF = (
+    'edition = 2\n'
+    'objective = expression\n'
+    'expression = 0.5*((kcat - 7)/0.1)^2 + 0.5*((Km - 70)/1)^2 + 0.5*((alpha - 0.3)/0.01)^2\n'
+    + _vars_conf(DECLARED) +
+    'starting_params = 7 70 0.3\n'
+    'job_type = am\n'
+    'max_iterations = 40\n'
+    'burn_in = 10\n'
+    'adaptive = 10\n'
+    'sample_every = 1\n'
+    'population_size = 2\n'
+    'output_hist_every = 1000000000\n'
+    'credible_intervals = 68\n'
+    'random_seed = 1\n'
+    'output_dir = out\n'
+    'verbosity = 0\n'
+    'backup_every = 1000000000\n'
+    'output_every = 1000000000\n'
+    'diagnostics_every = 1000000000\n')
+
+
+def _run_conf(tmp_path, monkeypatch, conf_text):
+    from pybnf.parse import load_config
+    from pybnf.pybnf import _create_algorithm
+    H.install(monkeypatch)
+    monkeypatch.chdir(tmp_path)
+    (tmp_path / 'order.conf').write_text(conf_text)
+    conf = load_config(str(tmp_path / 'order.conf'))
+    alg = _create_algorithm(conf)
+    H.drive(alg)
+    return conf, alg, tmp_path / 'out' / 'Results'
+
+
+def test_sampler_files_and_starting_params_order(tmp_path, monkeypatch):
+    _, alg, results = _run_conf(tmp_path, monkeypatch, _SAMPLER_CONF)
+    assert [v.name for v in alg.variables] == DECLARED
+
+    # starting_params = 7 70 0.3 is read by declaration order: kcat 7, Km 70, alpha 0.3.
+    # start_point.txt lists the rows in declaration order.
+    assert _declared_rows(results / 'start_point.txt') == [
+        ('kcat', 7.0), ('Km', 70.0), ('alpha', 0.3)]
+
+    # samples.txt: Name, Ln_probability, then the parameters in sorted order.
+    names, rows = _read_table(results / 'samples.txt', 2)
+    assert names == SORTED_IN_FILES
+    _assert_each_value_in_the_box_its_header_names(rows, 'samples.txt')
+
+    # sorted_params_final.txt: Simulation, Obj, then the parameters in sorted order.
+    names, rows = _read_table(results / 'sorted_params_final.txt', 2)
+    assert names == SORTED_IN_FILES
+    _assert_each_value_in_the_box_its_header_names(rows, 'sorted_params_final.txt')
+
+    # The credible-interval rows are in declaration order (#932).
+    credible = sorted(results.glob('credible68*_final.txt'))
+    assert len(credible) == 1
+    assert [n for n, _ in _declared_rows(credible[0])] == DECLARED
+
+    # am's own posterior files are in declaration order too.
+    with open(results / 'A_MCMC' / 'Runs' / 'params_0.txt') as f:
+        assert f.readline().split() == DECLARED
+
+
+def test_continue_run_overrides_starting_params(tmp_path, monkeypatch):
+    """With ``continue_run = 1`` every chain starts from the previous run's
+    ``adaptive_files/MLE_params.txt`` (one value per parameter, declaration order), and
+    ``starting_params`` is ignored."""
+    from pybnf.parse import load_config
+    from pybnf.pybnf import _create_algorithm
+    H.install(monkeypatch)
+    monkeypatch.chdir(tmp_path)
+    adaptive = tmp_path / 'out' / 'adaptive_files'
+    adaptive.mkdir(parents=True)
+    previous_mle = [3.0, 60.0, 0.5]                       # kcat, Km, alpha
+    np.savetxt(adaptive / 'MLE_params.txt', previous_mle)
+    np.savetxt(adaptive / 'diff.txt', [0.1])
+    np.savetxt(adaptive / 'diffMatrix.txt', np.eye(3))
+    (tmp_path / 'order.conf').write_text(_SAMPLER_CONF + 'continue_run = 1\n')
+    alg = _create_algorithm(load_config(str(tmp_path / 'order.conf')))
+    os.makedirs(alg.res_dir, exist_ok=True)
+
+    first = alg.start_run()
+
+    assert len(first) == 2
+    for pset in first:
+        assert [(p.name, p.value) for p in pset] == list(zip(DECLARED, previous_mle))
+
+
+# --------------------------------------------------------------------------- #
+# An optimizer with bootstrapping: sorted_params.txt (delete_old_files = 2), each
+# replicate's sorted_params.txt, and bootstrapped_parameter_sets.txt are sorted.
+# --------------------------------------------------------------------------- #
+
+class _NamedLinearModel(H.LinearGaussianModel):
+    """``f(x) = A x`` with ``x = (kcat, Km, alpha)`` read by name, scored by chi_sq
+    against real data, so the bootstrap has data rows to resample."""
+
+    def execute(self, folder, filename, timeout):
+        from pybnf.data import Data
+        x = np.array([self._pset[n] for n in DECLARED])
+        y = self.A @ x
+        d = Data(arr=np.column_stack([np.arange(len(y), dtype=float), y]))
+        d.cols = {'index': 0, 'y': 1}
+        d.headers = {0: 'index', 1: 'y'}
+        return {'target': d}
+
+
+def _bootstrap_config(tmp_path):
+    A = np.vstack([np.eye(3), np.eye(3)])                  # each parameter measured twice
+    truth = np.array([TRUE[n] for n in DECLARED])
+    sigma = np.tile([0.1, 1.0, 0.01], 2)
+    tgt, _ = H.write_target(tmp_path, H.gaussian_spec([0.0] * 3, [1.0] * 3))
+    exp = tmp_path / 'target.exp'
+    exp.write_text('# index\ty\ty_SD\n' + ''.join(
+        '%d\t%.17g\t%.17g\n' % (i, v, s) for i, (v, s) in enumerate(zip(A @ truth, sigma))))
+    base = {
+        'output_dir': str(tmp_path / 'out'),
+        'models': {tgt}, tgt: [str(exp)], 'exp_data': {str(exp)},
+        'objfunc': 'chi_sq', 'fit_type': 'de', 'initialization': 'lh',
+        'population_size': 6, 'max_iterations': 3, 'bootstrap': 2,
+        'delete_old_files': 2, 'verbosity': 0, 'wall_time_sim': 0, 'random_seed': 7,
+    }
+    base.update({('uniform_var', n): list(BOXES[n]) for n in DECLARED})
+    conf = config.Configuration(base)
+    conf.models['g'] = _NamedLinearModel(A, name='g')
+    return conf
+
+
+def test_optimizer_and_bootstrap_files_are_sorted(tmp_path, monkeypatch):
+    from pybnf.pybnf import _run_bootstrapping
+    H.install(monkeypatch)
+    conf = _bootstrap_config(tmp_path)
+    alg = algorithms.DifferentialEvolution(conf)
+    assert [v.name for v in alg.variables] == DECLARED
+    H.drive(alg)
+    out = tmp_path / 'out'
+    assert [n for n, _ in _declared_rows(out / 'Results' / 'start_point.txt')] == DECLARED
+
+    # delete_old_files = 2 collapses every sorted_params_<name>.txt into sorted_params.txt.
+    names, rows = _read_table(out / 'Results' / 'sorted_params.txt', 2)
+    assert names == SORTED_IN_FILES
+    _assert_each_value_in_the_box_its_header_names(rows, 'sorted_params.txt')
+
+    _run_bootstrapping(conf, alg, types.SimpleNamespace(client=H.FakeClient()), False)
+
+    names, rows = _read_table(out / 'Results' / 'bootstrapped_parameter_sets.txt', 2)
+    assert names == SORTED_IN_FILES
+    assert len(rows) == 2
+    _assert_each_value_in_the_box_its_header_names(rows, 'bootstrapped_parameter_sets.txt')
+    for n in range(2):
+        names, rows = _read_table(out / ('Results-boot%d' % n) / 'sorted_params.txt', 2)
+        assert names == SORTED_IN_FILES
+        _assert_each_value_in_the_box_its_header_names(rows, 'Results-boot%d' % n)
+        assert [r for r, _ in _declared_rows(out / ('Results-boot%d' % n) / 'start_point.txt')] \
+            == DECLARED