diff --git a/README.md b/README.md index 3759e3d..88471a7 100644 --- a/README.md +++ b/README.md @@ -27,6 +27,8 @@ The engine diagonalizes the GRM once, rotates the common design and every unique - Variable numbers of extra fixed-effect covariates per task through a CSR index. - Satterthwaite t/F fixed-effect inference by default, with optional Kenward-Roger F tests and a joint test of task-specific covariates. +- Named arbitrary linear contrasts over the common fixed-effect design, using + the Satterthwaite or Kenward-Roger method selected for the run. - Atomic block output and safe resume/force semantics. - Generic Python CLI for manifest creation, validation, execution, provenance signatures, status, and result export. @@ -114,11 +116,14 @@ python python/spectra_reml.py make-manifest \ --tasks data/tasks.tsv \ --extra-offsets data/extra_offsets.i64.bin \ --extra-indices data/extra_indices.i32.bin \ + --contrast-matrix data/contrasts.f64.bin \ + --contrast-metadata data/contrasts.tsv \ --output-dir work/blocks \ --n-samples 3523 \ --n-base-covariates 8 \ --n-phenotype-rows 10000 \ - --n-extra-covariate-rows 2000 + --n-extra-covariate-rows 2000 \ + --n-contrasts 5 ``` Validate, run, and export results: @@ -141,6 +146,14 @@ python python/spectra_reml.py finalize \ `results.tsv.gz` retains the complete per-task summary and stores the fixed-effect vector and row-wise packed lower covariance as JSON arrays. +Optional named contrasts are supplied by `--contrast-matrix`, +`--contrast-metadata`, and `--n-contrasts`. Contrast rows act on `base_x`; +task-specific covariate columns receive zero weights automatically. Finalized +results contain contrast estimates, adjusted standard errors, statistics, +numerator and denominator degrees of freedom, and p-values as JSON arrays. +This release uses manifest v2, run-signature v2, block v3, and finalized-output +v2 contracts. It intentionally does not read older contracts. + Fixed-effect inference defaults to Satterthwaite. Select Kenward-Roger or turn inference off with `--fixed-effect-test kenward-roger` or `--fixed-effect-test none`. Coefficient-wise standard errors, statistics, diff --git a/docs/FORMAT.md b/docs/FORMAT.md index 92d1786..2607db7 100644 --- a/docs/FORMAT.md +++ b/docs/FORMAT.md @@ -84,6 +84,22 @@ indices[offsets[i]:offsets[i+1]] An index must be in `[0, n_extra_covariate_rows)`, and a task cannot reference the same row twice. +## Named shared fixed-effect contrasts + +Optional arbitrary hypotheses use two files. `--contrast-matrix` is a +row-major `float64` matrix with shape `n_contrasts x n_base_covariates`; its +columns follow `base_x` exactly. Task-specific extra fixed-effect columns are +assigned zero weights automatically. `--contrast-metadata` is a UTF-8 TSV: + +```text +contrast_index contrast_id rhs +0 group_a_b 0 +``` + +Indices are consecutive and zero based, IDs are unique, and every matrix row +must contain a nonzero weight. Each row tests `L beta = rhs` with the method +selected by `--fixed-effect-test` (Satterthwaite or Kenward-Roger). + ## Block output For block number `KKKKKK`: @@ -96,6 +112,12 @@ block_KKKKKK.fixed_se.f64.bin block_KKKKKK.fixed_stat.f64.bin block_KKKKKK.fixed_ddf.f64.bin block_KKKKKK.fixed_p.f64.bin +block_KKKKKK.contrast_estimate.f64.bin +block_KKKKKK.contrast_se.f64.bin +block_KKKKKK.contrast_stat.f64.bin +block_KKKKKK.contrast_numdf.f64.bin +block_KKKKKK.contrast_ddf.f64.bin +block_KKKKKK.contrast_p.f64.bin block_KKKKKK.complete ``` @@ -119,6 +141,7 @@ grad_inf fixed_test_method fixed_test_status fixed_test_offset +contrast_test_offset extra_joint_num_df extra_joint_den_df extra_joint_f @@ -142,12 +165,13 @@ The covariance array uses the row-wise packed lower triangle: The `.complete` marker is written last and contains tab-separated key/value rows: ```text -format spectra-reml-block-v2 +format spectra-reml-block-v3 block 0 tasks 256 beta_elements 4096 cov_elements 34816 fixed_test_elements 4096 +contrast_test_elements 1280 ``` Consumers must ignore blocks without `.complete`. @@ -180,6 +204,13 @@ one numerator degree of freedom. For every task with extra covariates, the summary also contains an F test of the joint null that all task-specific fixed effects are zero. +The six `contrast_*` arrays share `contrast_test_offset` and follow the +contrast metadata order. A one-row Satterthwaite contrast reports a signed t +statistic; Kenward-Roger reports an F statistic with one numerator degree of +freedom. Both methods retain the signed `L beta` estimate and their adjusted +standard error. Negative `contrast_test_offset` means that no valid contrast +tests were emitted for that task. + `fixed_test_status` is one of: ```text @@ -218,9 +249,13 @@ fixed_effect_se_json fixed_effect_statistic_json fixed_effect_denominator_df_json fixed_effect_p_value_json +contrast_estimate_json +contrast_standard_error_json +contrast_statistic_json +contrast_numerator_df_json +contrast_denominator_df_json +contrast_p_value_json ``` -Project-specific software can attach coefficient names. The public C++ -inference API also accepts general linear hypotheses `L beta = rhs`; the batch -format currently emits coefficient-wise tests and the joint extra-covariate -test. +Contrast JSON arrays follow `contrast_index`; IDs and right-hand sides are +stored in the manifest-referenced contrast metadata TSV. diff --git a/include/spectra_reml/batch_io.hpp b/include/spectra_reml/batch_io.hpp index 4f7f468..be0ecd7 100644 --- a/include/spectra_reml/batch_io.hpp +++ b/include/spectra_reml/batch_io.hpp @@ -20,6 +20,8 @@ struct BatchInputPaths { std::filesystem::path tasks_tsv; std::filesystem::path extra_offsets_i64; std::filesystem::path extra_indices_i32; + std::filesystem::path contrast_matrix_f64; + std::filesystem::path contrast_metadata_tsv; std::filesystem::path output_directory; }; @@ -28,6 +30,7 @@ struct BatchDimensions { std::size_t base_covariate_count = 0; std::size_t phenotype_row_count = 0; std::size_t extra_covariate_row_count = 0; + std::size_t contrast_count = 0; }; struct BatchOptions { diff --git a/include/spectra_reml/types.hpp b/include/spectra_reml/types.hpp index 788dbef..7b1af02 100644 --- a/include/spectra_reml/types.hpp +++ b/include/spectra_reml/types.hpp @@ -164,6 +164,10 @@ struct RemlResult { std::vector beta_covariance_packed_lower; FixedEffectInferenceResult fixed_effect_inference; FixedEffectTestResult extra_fixed_effect_joint_test; + // Named shared contrasts supplied by the batch interface, in manifest + // contrast order. Every contrast acts on the common/base fixed effects; + // task-specific fixed-effect columns are assigned zero weight. + std::vector contrast_tests; std::string error; [[nodiscard]] bool has_estimates() const noexcept { diff --git a/python/spectra_reml.py b/python/spectra_reml.py index ccb9cc0..14cbf04 100644 --- a/python/spectra_reml.py +++ b/python/spectra_reml.py @@ -31,10 +31,10 @@ except ImportError as exc: # pragma: no cover raise SystemExit("NumPy is required: {}".format(exc)) -MANIFEST_FORMAT = "spectra-reml-manifest-v1" -BLOCK_FORMAT = "spectra-reml-block-v2" -RUN_SIGNATURE_FORMAT = "spectra-reml-run-signature-v1" -FINALIZE_FORMAT = "spectra-reml-finalize-v1" +MANIFEST_FORMAT = "spectra-reml-manifest-v2" +BLOCK_FORMAT = "spectra-reml-block-v3" +RUN_SIGNATURE_FORMAT = "spectra-reml-run-signature-v2" +FINALIZE_FORMAT = "spectra-reml-finalize-v2" TASK_HEADER = ( "task_index", "task_id", @@ -59,6 +59,7 @@ SUMMARY_HEADER = ( "fixed_test_method", "fixed_test_status", "fixed_test_offset", + "contrast_test_offset", "extra_joint_num_df", "extra_joint_den_df", "extra_joint_f", @@ -75,6 +76,12 @@ BLOCK_PATTERNS = ( "block_*.fixed_stat.f64.bin", "block_*.fixed_ddf.f64.bin", "block_*.fixed_p.f64.bin", + "block_*.contrast_estimate.f64.bin", + "block_*.contrast_se.f64.bin", + "block_*.contrast_stat.f64.bin", + "block_*.contrast_numdf.f64.bin", + "block_*.contrast_ddf.f64.bin", + "block_*.contrast_p.f64.bin", "block_*.complete", ) @@ -180,6 +187,8 @@ def manifest_paths(manifest: Mapping[str, Any]) -> dict[str, Path | None]: "extra_offsets_i64", "extra_indices_i32", "output_directory", + "contrast_matrix_f64", + "contrast_metadata_tsv", ) missing = [name for name in required if name not in raw] if missing: @@ -199,6 +208,30 @@ def read_tasks(path: Path) -> list[dict[str, str]]: return list(reader) +def read_contrast_metadata(path: Path) -> list[dict[str, str]]: + with path.open("r", encoding="utf-8-sig", newline="") as handle: + reader = csv.DictReader(handle, delimiter="\t") + expected = ("contrast_index", "contrast_id", "rhs") + if tuple(reader.fieldnames or ()) != expected: + fail("Contrast metadata header must be exactly: {}".format("\t".join(expected))) + rows = list(reader) + seen: set[str] = set() + for expected_index, row in enumerate(rows): + try: + index = int(row["contrast_index"]) + rhs = float(row["rhs"]) + except ValueError as exc: + fail("Invalid contrast metadata at row {}: {}".format(expected_index, exc)) + if index != expected_index: + fail("Contrast indices must be consecutive and zero based.") + if not row["contrast_id"] or row["contrast_id"] in seen: + fail("Contrast IDs must be nonempty and unique.") + if not np.isfinite(rhs): + fail("Contrast rhs must be finite.") + seen.add(row["contrast_id"]) + return rows + + def validate_manifest(path: Path) -> dict[str, Any]: path = resolved(path) manifest = read_json(path) @@ -214,6 +247,7 @@ def validate_manifest(path: Path) -> dict[str, Any]: "extra_covariate_row_count", "task_count", "extra_index_count", + "contrast_count", ) parsed: dict[str, int] = {} for name in names: @@ -240,6 +274,28 @@ def validate_manifest(path: Path) -> dict[str, Any]: with paths["grm_id"].open("r", encoding="utf-8-sig") as handle: if sum(bool(line.strip()) for line in handle) != n: fail("GRM ID row count differs from sample_count.") + contrast_paths = (paths["contrast_matrix_f64"], paths["contrast_metadata_tsv"]) + if parsed["contrast_count"]: + if any(path is None for path in contrast_paths): + fail("Both contrast_matrix_f64 and contrast_metadata_tsv are required.") + require_size( + paths["contrast_matrix_f64"], + parsed["contrast_count"] * parsed["base_covariate_count"] * 8, + "contrast matrix", + ) + contrast_rows = read_contrast_metadata(paths["contrast_metadata_tsv"]) + if len(contrast_rows) != parsed["contrast_count"]: + fail("Contrast metadata row count differs from contrast_count.") + contrast_matrix = np.fromfile(paths["contrast_matrix_f64"], dtype=" dict[str, Any]: "extra_offsets_i64": paths["extra_offsets_i64"], "extra_indices_i32": paths["extra_indices_i32"], } + if parsed["contrast_count"]: + controls["contrast_matrix_f64"] = paths["contrast_matrix_f64"] + controls["contrast_metadata_tsv"] = paths["contrast_metadata_tsv"] hashes = manifest.get("control_sha256") if not isinstance(hashes, dict) or set(hashes) != set(controls): fail("control_sha256 does not cover all control files.") @@ -314,6 +373,8 @@ def make_manifest(args: argparse.Namespace) -> Path: "tasks_tsv": resolved(args.tasks), "extra_offsets_i64": resolved(args.extra_offsets), "extra_indices_i32": resolved(args.extra_indices), + "contrast_matrix_f64": resolved(args.contrast_matrix) if args.contrast_matrix else None, + "contrast_metadata_tsv": resolved(args.contrast_metadata) if args.contrast_metadata else None, "output_directory": resolved(args.output_dir), } tasks = read_tasks(paths["tasks_tsv"]) # type: ignore[arg-type] @@ -323,7 +384,13 @@ def make_manifest(args: argparse.Namespace) -> Path: for name in ("grm_id", "extra_covariate_f32"): if paths[name] is not None: sources[name] = file_identity(paths[name]) # type: ignore[arg-type] - controls = ("base_x_f64", "tasks_tsv", "extra_offsets_i64", "extra_indices_i32") + controls = ["base_x_f64", "tasks_tsv", "extra_offsets_i64", "extra_indices_i32"] + if args.n_contrasts: + if paths["contrast_matrix_f64"] is None or paths["contrast_metadata_tsv"] is None: + fail("--contrast-matrix and --contrast-metadata are required when --n-contrasts is positive.") + controls.extend(("contrast_matrix_f64", "contrast_metadata_tsv")) + elif paths["contrast_matrix_f64"] is not None or paths["contrast_metadata_tsv"] is not None: + fail("Contrast paths require a positive --n-contrasts.") manifest = { "format": MANIFEST_FORMAT, "created_utc": utc_now(), @@ -334,6 +401,7 @@ def make_manifest(args: argparse.Namespace) -> Path: "extra_covariate_row_count": args.n_extra_covariate_rows, "task_count": len(tasks), "extra_index_count": extra_index_count, + "contrast_count": args.n_contrasts, }, "paths": {name: str(value) if value is not None else None for name, value in paths.items()}, "source_files": sources, @@ -371,6 +439,7 @@ def build_engine_command(engine: Path, manifest: Mapping[str, Any], args: argpar "--n-base-covariates", str(dims["base_covariate_count"]), "--n-phenotype-rows", str(dims["phenotype_row_count"]), "--n-extra-covariate-rows", str(dims["extra_covariate_row_count"]), + "--n-contrasts", str(dims["contrast_count"]), "--block-size", str(args.block_size), "--threads", str(args.threads), "--fixed-effect-test", args.fixed_effect_test, @@ -379,6 +448,9 @@ def build_engine_command(engine: Path, manifest: Mapping[str, Any], args: argpar command.extend(("--grm-id", str(paths["grm_id"]))) if int(dims["extra_index_count"]): command.extend(("--extra-covariates", str(paths["extra_covariate_f32"]))) + if int(dims["contrast_count"]): + command.extend(("--contrast-matrix", str(paths["contrast_matrix_f64"]))) + command.extend(("--contrast-metadata", str(paths["contrast_metadata_tsv"]))) option_map = ( ("max_iterations", "--max-iterations"), ("line_search_max_evals", "--line-search-max-evals"), @@ -571,16 +643,25 @@ def read_completed(manifest: Mapping[str, Any]) -> list[dict[str, Any]]: fixed_stat_path = output / (stem + ".fixed_stat.f64.bin") fixed_ddf_path = output / (stem + ".fixed_ddf.f64.bin") fixed_p_path = output / (stem + ".fixed_p.f64.bin") + contrast_estimate_path = output / (stem + ".contrast_estimate.f64.bin") + contrast_se_path = output / (stem + ".contrast_se.f64.bin") + contrast_stat_path = output / (stem + ".contrast_stat.f64.bin") + contrast_numdf_path = output / (stem + ".contrast_numdf.f64.bin") + contrast_ddf_path = output / (stem + ".contrast_ddf.f64.bin") + contrast_p_path = output / (stem + ".contrast_p.f64.bin") marker_values = read_marker(marker) + block_format = marker_values.get("format") try: declared_tasks = int(marker_values["tasks"]) beta_elements = int(marker_values["beta_elements"]) cov_elements = int(marker_values["cov_elements"]) fixed_test_elements = int(marker_values["fixed_test_elements"]) + contrast_test_elements = int(marker_values["contrast_test_elements"]) except (KeyError, ValueError) as exc: fail("Malformed completion counts in {}: {}".format(marker, exc)) - if marker_values.get("format") != BLOCK_FORMAT or min( - declared_tasks, beta_elements, cov_elements, fixed_test_elements + if block_format != BLOCK_FORMAT or min( + declared_tasks, beta_elements, cov_elements, fixed_test_elements, + contrast_test_elements, ) < 0: fail("Invalid completion marker: {}".format(marker)) require_size(beta_path, beta_elements * 8, "block beta") @@ -589,6 +670,15 @@ def read_completed(manifest: Mapping[str, Any]) -> list[dict[str, Any]]: require_size(fixed_stat_path, fixed_test_elements * 8, "fixed-effect statistic") require_size(fixed_ddf_path, fixed_test_elements * 8, "fixed-effect denominator df") require_size(fixed_p_path, fixed_test_elements * 8, "fixed-effect p-value") + for path, label in ( + (contrast_estimate_path, "contrast estimate"), + (contrast_se_path, "contrast SE"), + (contrast_stat_path, "contrast statistic"), + (contrast_numdf_path, "contrast numerator df"), + (contrast_ddf_path, "contrast denominator df"), + (contrast_p_path, "contrast p-value"), + ): + require_size(path, contrast_test_elements * 8, label) with summary_path.open("r", encoding="utf-8-sig", newline="") as handle: reader = csv.DictReader(handle, delimiter="\t") if tuple(reader.fieldnames or ()) != SUMMARY_HEADER: @@ -602,9 +692,17 @@ def read_completed(manifest: Mapping[str, Any]) -> list[dict[str, Any]]: fixed_stat = np.fromfile(fixed_stat_path, dtype="= task_count: @@ -612,6 +710,7 @@ def read_completed(manifest: Mapping[str, Any]) -> list[dict[str, Any]]: p = int(row["n_fixed"]) beta_offset, cov_offset = int(row["beta_offset"]), int(row["cov_offset"]) fixed_test_offset = int(row["fixed_test_offset"]) + contrast_test_offset = int(row["contrast_test_offset"]) if beta_offset < 0 or cov_offset < 0: if beta_offset != cov_offset: fail("Only one output offset is negative for task {}".format(index)) @@ -641,6 +740,25 @@ def read_completed(manifest: Mapping[str, Any]) -> list[dict[str, Any]]: fixed_ddf_values = fixed_ddf[fixed_test_offset:fixed_test_offset + p].tolist() fixed_p_values = fixed_p[fixed_test_offset:fixed_test_offset + p].tolist() next_fixed_test_offset += p + if contrast_test_offset < 0: + contrast_estimate_values: list[float] = [] + contrast_se_values: list[float] = [] + contrast_stat_values: list[float] = [] + contrast_numdf_values: list[float] = [] + contrast_ddf_values: list[float] = [] + contrast_p_values: list[float] = [] + else: + if (contrast_test_offset != next_contrast_test_offset or + contrast_test_offset + q > contrast_estimate.size): + fail("Invalid contrast test offset for task {}".format(index)) + section = slice(contrast_test_offset, contrast_test_offset + q) + contrast_estimate_values = contrast_estimate[section].tolist() + contrast_se_values = contrast_se[section].tolist() + contrast_stat_values = contrast_stat[section].tolist() + contrast_numdf_values = contrast_numdf[section].tolist() + contrast_ddf_values = contrast_ddf[section].tolist() + contrast_p_values = contrast_p[section].tolist() + next_contrast_test_offset += q result: dict[str, Any] = dict(row) result["beta_json"] = json.dumps(beta_values, separators=(",", ":")) result["covariance_packed_lower_json"] = json.dumps(cov_values, separators=(",", ":")) @@ -648,9 +766,16 @@ def read_completed(manifest: Mapping[str, Any]) -> list[dict[str, Any]]: result["fixed_effect_statistic_json"] = json.dumps(fixed_stat_values, separators=(",", ":")) result["fixed_effect_denominator_df_json"] = json.dumps(fixed_ddf_values, separators=(",", ":")) result["fixed_effect_p_value_json"] = json.dumps(fixed_p_values, separators=(",", ":")) + result["contrast_estimate_json"] = json.dumps(contrast_estimate_values, separators=(",", ":")) + result["contrast_standard_error_json"] = json.dumps(contrast_se_values, separators=(",", ":")) + result["contrast_statistic_json"] = json.dumps(contrast_stat_values, separators=(",", ":")) + result["contrast_numerator_df_json"] = json.dumps(contrast_numdf_values, separators=(",", ":")) + result["contrast_denominator_df_json"] = json.dumps(contrast_ddf_values, separators=(",", ":")) + result["contrast_p_value_json"] = json.dumps(contrast_p_values, separators=(",", ":")) results[index] = result if (next_beta_offset != beta_elements or next_cov_offset != cov_elements or - next_fixed_test_offset != fixed_test_elements): + next_fixed_test_offset != fixed_test_elements or + next_contrast_test_offset != contrast_test_elements): fail("Block binary arrays contain unused elements: {}".format(stem)) missing = sorted(set(range(task_count)).difference(results)) if missing: @@ -669,6 +794,9 @@ def finalize(manifest: Mapping[str, Any], output: Path) -> Path: "beta_json", "covariance_packed_lower_json", "fixed_effect_se_json", "fixed_effect_statistic_json", "fixed_effect_denominator_df_json", "fixed_effect_p_value_json", + "contrast_estimate_json", "contrast_standard_error_json", + "contrast_statistic_json", "contrast_numerator_df_json", + "contrast_denominator_df_json", "contrast_p_value_json", ] try: with opener(temporary, "wt", encoding="utf-8", newline="") as handle: @@ -771,11 +899,14 @@ def build_parser() -> argparse.ArgumentParser: make.add_argument("--tasks", type=Path, required=True) make.add_argument("--extra-offsets", type=Path, required=True) make.add_argument("--extra-indices", type=Path, required=True) + make.add_argument("--contrast-matrix", type=Path) + make.add_argument("--contrast-metadata", type=Path) make.add_argument("--output-dir", type=Path, required=True) make.add_argument("--n-samples", type=positive_integer, required=True) make.add_argument("--n-base-covariates", type=positive_integer, required=True) make.add_argument("--n-phenotype-rows", type=positive_integer, required=True) make.add_argument("--n-extra-covariate-rows", type=int, default=0) + make.add_argument("--n-contrasts", type=int, default=0) for action in ("validate", "status"): command = subparsers.add_parser(action) command.add_argument("--manifest", type=Path, required=True) diff --git a/src/batch_io.cpp b/src/batch_io.cpp index 110c98d..feb1a25 100644 --- a/src/batch_io.cpp +++ b/src/batch_io.cpp @@ -34,7 +34,96 @@ namespace spectra::reml { namespace { -constexpr const char* kOutputFormatVersion = "spectra-reml-block-v2"; +constexpr const char* kOutputFormatVersion = "spectra-reml-block-v3"; + +std::vector split_tab(const std::string& line); +std::uint64_t parse_u64(const std::string& text, const char* field, + std::size_t line_number); + +struct NamedContrast { + std::string id; + double rhs = 0.0; + std::vector weights; +}; + +std::vector read_named_contrasts( + const BatchInputPaths& paths, const BatchDimensions& dimensions) { + if (dimensions.contrast_count == 0) { + return {}; + } + const ColMajorMatrix matrix = read_row_major_f64_matrix( + paths.contrast_matrix_f64, dimensions.contrast_count, + dimensions.base_covariate_count); + std::ifstream input(paths.contrast_metadata_tsv); + if (!input) { + throw std::runtime_error("cannot open contrast metadata TSV " + + paths.contrast_metadata_tsv.string()); + } + std::string line; + if (!std::getline(input, line)) { + throw std::runtime_error("contrast metadata TSV is empty"); + } + if (!line.empty() && line.back() == '\r') line.pop_back(); + if (split_tab(line) != + std::vector{"contrast_index", "contrast_id", "rhs"}) { + throw std::runtime_error( + "contrast metadata header must be contrast_index, contrast_id, rhs"); + } + std::vector contrasts; + contrasts.reserve(dimensions.contrast_count); + std::unordered_set ids; + std::size_t line_number = 1; + while (std::getline(input, line)) { + ++line_number; + if (!line.empty() && line.back() == '\r') line.pop_back(); + if (line.empty()) continue; + const auto fields = split_tab(line); + if (fields.size() != 3) { + throw std::runtime_error("invalid contrast metadata field count at line " + + std::to_string(line_number)); + } + const auto index = parse_u64(fields[0], "contrast_index", line_number); + if (index != contrasts.size()) { + throw std::runtime_error( + "contrast indices must be consecutive and zero based"); + } + if (fields[1].empty() || !ids.insert(fields[1]).second) { + throw std::runtime_error( + "contrast IDs must be nonempty and unique at line " + + std::to_string(line_number)); + } + std::size_t consumed = 0; + double rhs = std::numeric_limits::quiet_NaN(); + try { + rhs = std::stod(fields[2], &consumed); + } catch (const std::exception&) { + throw std::runtime_error("invalid contrast rhs at line " + + std::to_string(line_number)); + } + if (consumed != fields[2].size() || !std::isfinite(rhs)) { + throw std::runtime_error("invalid contrast rhs at line " + + std::to_string(line_number)); + } + NamedContrast contrast; + contrast.id = fields[1]; + contrast.rhs = rhs; + contrast.weights.resize(dimensions.base_covariate_count); + bool any_nonzero = false; + for (std::size_t col = 0; col < dimensions.base_covariate_count; ++col) { + contrast.weights[col] = matrix(index, col); + any_nonzero = any_nonzero || contrast.weights[col] != 0.0; + } + if (!any_nonzero) { + throw std::runtime_error("contrast " + contrast.id + + " has all-zero weights"); + } + contrasts.push_back(std::move(contrast)); + } + if (contrasts.size() != dimensions.contrast_count) { + throw std::runtime_error("contrast metadata row count mismatch"); + } + return contrasts; +} std::size_t checked_product(std::size_t left, std::size_t right, const char* description) { @@ -295,10 +384,24 @@ void write_block_atomic(const std::filesystem::path& output_directory, output_directory / (stem + ".fixed_stat.f64.bin"); const auto fixed_ddf_path = output_directory / (stem + ".fixed_ddf.f64.bin"); const auto fixed_p_path = output_directory / (stem + ".fixed_p.f64.bin"); + const auto contrast_estimate_path = + output_directory / (stem + ".contrast_estimate.f64.bin"); + const auto contrast_se_path = + output_directory / (stem + ".contrast_se.f64.bin"); + const auto contrast_stat_path = + output_directory / (stem + ".contrast_stat.f64.bin"); + const auto contrast_numdf_path = + output_directory / (stem + ".contrast_numdf.f64.bin"); + const auto contrast_ddf_path = + output_directory / (stem + ".contrast_ddf.f64.bin"); + const auto contrast_p_path = + output_directory / (stem + ".contrast_p.f64.bin"); const auto complete_path = output_directory / (stem + ".complete"); - const std::array final_paths = { + const std::array final_paths = { summary_path, beta_path, covariance_path, fixed_se_path, - fixed_stat_path, fixed_ddf_path, fixed_p_path, complete_path}; + fixed_stat_path, fixed_ddf_path, fixed_p_path, contrast_estimate_path, + contrast_se_path, contrast_stat_path, contrast_numdf_path, + contrast_ddf_path, contrast_p_path, complete_path}; for (const auto& path : final_paths) { if (std::filesystem::exists(path)) { if (!overwrite) { @@ -314,9 +417,16 @@ void write_block_atomic(const std::filesystem::path& output_directory, std::vector fixed_stat_values; std::vector fixed_ddf_values; std::vector fixed_p_values; + std::vector contrast_estimate_values; + std::vector contrast_se_values; + std::vector contrast_stat_values; + std::vector contrast_numdf_values; + std::vector contrast_ddf_values; + std::vector contrast_p_values; std::vector beta_offsets(results.size(), -1); std::vector covariance_offsets(results.size(), -1); std::vector fixed_test_offsets(results.size(), -1); + std::vector contrast_test_offsets(results.size(), -1); for (std::size_t index = 0; index < results.size(); ++index) { const auto& fit = results[index].fit; if (!fit.has_estimates()) { @@ -342,6 +452,19 @@ void write_block_atomic(const std::filesystem::path& output_directory, fixed_ddf_values.push_back(test.denominator_df); fixed_p_values.push_back(test.p_value); } + if (!fit.contrast_tests.empty()) { + contrast_test_offsets[index] = static_cast( + contrast_estimate_values.size()); + for (const auto& test : fit.contrast_tests) { + contrast_estimate_values.push_back(test.estimate); + contrast_se_values.push_back(test.standard_error); + contrast_stat_values.push_back(test.statistic); + contrast_numdf_values.push_back( + static_cast(test.numerator_df)); + contrast_ddf_values.push_back(test.denominator_df); + contrast_p_values.push_back(test.p_value); + } + } } } @@ -352,10 +475,18 @@ void write_block_atomic(const std::filesystem::path& output_directory, const auto fixed_stat_temp = temporary_path(fixed_stat_path); const auto fixed_ddf_temp = temporary_path(fixed_ddf_path); const auto fixed_p_temp = temporary_path(fixed_p_path); + const auto contrast_estimate_temp = temporary_path(contrast_estimate_path); + const auto contrast_se_temp = temporary_path(contrast_se_path); + const auto contrast_stat_temp = temporary_path(contrast_stat_path); + const auto contrast_numdf_temp = temporary_path(contrast_numdf_path); + const auto contrast_ddf_temp = temporary_path(contrast_ddf_path); + const auto contrast_p_temp = temporary_path(contrast_p_path); const auto complete_temp = temporary_path(complete_path); - const std::array temporary_paths = { + const std::array temporary_paths = { summary_temp, beta_temp, covariance_temp, fixed_se_temp, - fixed_stat_temp, fixed_ddf_temp, fixed_p_temp, complete_temp}; + fixed_stat_temp, fixed_ddf_temp, fixed_p_temp, contrast_estimate_temp, + contrast_se_temp, contrast_stat_temp, contrast_numdf_temp, + contrast_ddf_temp, contrast_p_temp, complete_temp}; try { write_binary_doubles(beta_temp, beta_values); write_binary_doubles(covariance_temp, covariance_values); @@ -363,6 +494,12 @@ void write_block_atomic(const std::filesystem::path& output_directory, write_binary_doubles(fixed_stat_temp, fixed_stat_values); write_binary_doubles(fixed_ddf_temp, fixed_ddf_values); write_binary_doubles(fixed_p_temp, fixed_p_values); + write_binary_doubles(contrast_estimate_temp, contrast_estimate_values); + write_binary_doubles(contrast_se_temp, contrast_se_values); + write_binary_doubles(contrast_stat_temp, contrast_stat_values); + write_binary_doubles(contrast_numdf_temp, contrast_numdf_values); + write_binary_doubles(contrast_ddf_temp, contrast_ddf_values); + write_binary_doubles(contrast_p_temp, contrast_p_values); { std::ofstream summary(summary_temp, std::ios::trunc); if (!summary) { @@ -374,6 +511,7 @@ void write_block_atomic(const std::filesystem::path& output_directory, "\tcov_offset\tsigma_g2\tsigma_e2\th2\tlogL\titerations" "\tline_search_steps\tgrad_inf\tfixed_test_method" "\tfixed_test_status\tfixed_test_offset" + "\tcontrast_test_offset" "\textra_joint_num_df\textra_joint_den_df" "\textra_joint_f\textra_joint_p\tfixed_test_error\terror\n"; summary << std::setprecision(17); @@ -396,6 +534,7 @@ void write_block_atomic(const std::filesystem::path& output_directory, << to_string(fit.fixed_effect_inference.method) << '\t' << to_string(fit.fixed_effect_inference.status) << '\t' << fixed_test_offsets[index] << '\t' + << contrast_test_offsets[index] << '\t' << fit.extra_fixed_effect_joint_test.numerator_df << '\t' << fit.extra_fixed_effect_joint_test.denominator_df << '\t' << fit.extra_fixed_effect_joint_test.statistic << '\t' @@ -421,6 +560,9 @@ void write_block_atomic(const std::filesystem::path& output_directory, << "beta_elements\t" << beta_values.size() << '\n' << "cov_elements\t" << covariance_values.size() << '\n' << "fixed_test_elements\t" << fixed_se_values.size() + << '\n' + << "contrast_test_elements\t" + << contrast_estimate_values.size() << '\n'; complete.flush(); if (!complete) { @@ -434,6 +576,12 @@ void write_block_atomic(const std::filesystem::path& output_directory, rename_checked(fixed_stat_temp, fixed_stat_path); rename_checked(fixed_ddf_temp, fixed_ddf_path); rename_checked(fixed_p_temp, fixed_p_path); + rename_checked(contrast_estimate_temp, contrast_estimate_path); + rename_checked(contrast_se_temp, contrast_se_path); + rename_checked(contrast_stat_temp, contrast_stat_path); + rename_checked(contrast_numdf_temp, contrast_numdf_path); + rename_checked(contrast_ddf_temp, contrast_ddf_path); + rename_checked(contrast_p_temp, contrast_p_path); rename_checked(summary_temp, summary_path); rename_checked(complete_temp, complete_path); } catch (...) { @@ -624,6 +772,7 @@ void run_task_batch(const BatchInputPaths& paths, if (tasks.empty()) { throw std::runtime_error("tasks TSV contains no tasks"); } + const auto named_contrasts = read_named_contrasts(paths, dimensions); std::cerr << "Reading and diagonalizing the GRM for " << dimensions.sample_count << " samples...\n"; const SpectralGrm spectral = read_and_diagonalize_gcta_grm( @@ -750,6 +899,7 @@ void run_task_batch(const BatchInputPaths& paths, item.fit = fit_ai_reml_spectral( phenotype, design, spectral.eigenvalues, options.reml); std::vector hypotheses; + const bool has_extra_joint = !task.extra_covariate_rows.empty(); if (!task.extra_covariate_rows.empty()) { FixedEffectHypothesis extra_joint; extra_joint.contrast = ColMajorMatrix( @@ -761,14 +911,35 @@ void run_task_batch(const BatchInputPaths& paths, } hypotheses.push_back(std::move(extra_joint)); } + for (const auto& named : named_contrasts) { + FixedEffectHypothesis hypothesis; + hypothesis.contrast = ColMajorMatrix(1, fixed_count); + for (std::size_t col = 0; + col < dimensions.base_covariate_count; ++col) { + hypothesis.contrast(0, col) = named.weights[col]; + } + hypothesis.rhs = {named.rhs}; + hypotheses.push_back(std::move(hypothesis)); + } item.fit.fixed_effect_inference = infer_fixed_effects_spectral( phenotype, design, spectral.eigenvalues, item.fit, options.reml.fixed_effect_test, hypotheses, options.reml.rank_tolerance_relative, options.reml.covariance_floor_relative); - if (!item.fit.fixed_effect_inference.hypothesis_tests.empty()) { + const auto& tests = + item.fit.fixed_effect_inference.hypothesis_tests; + std::size_t contrast_begin = 0; + if (has_extra_joint && !tests.empty()) { item.fit.extra_fixed_effect_joint_test = - item.fit.fixed_effect_inference.hypothesis_tests.front(); + tests.front(); + contrast_begin = 1; + } + if (tests.size() >= contrast_begin + named_contrasts.size()) { + item.fit.contrast_tests.assign( + tests.begin() + static_cast(contrast_begin), + tests.begin() + static_cast( + contrast_begin + + named_contrasts.size())); } } catch (const std::exception& exception) { item.fit.status = FitStatus::numerical_error; diff --git a/src/main.cpp b/src/main.cpp index 988dadf..b59910c 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -29,12 +29,15 @@ void print_usage(std::ostream& output) { << " --tasks PATH TSV with task_index,task_id,phenotype_row,n_extra_covariates\n" << " --extra-offsets PATH int64 CSR offsets, m+1 entries\n" << " --extra-indices PATH int32 zero-based extra-covariate row indices\n" + << " --contrast-matrix PATH optional row-major q x p0 float64 matrix\n" + << " --contrast-metadata PATH optional TSV: contrast_index,contrast_id,rhs\n" << " --out-dir PATH task-set output directory\n\n" << "Required dimensions:\n" << " --n-samples N\n" << " --n-base-covariates P\n" << " --n-phenotype-rows R\n" << " --n-extra-covariate-rows R\n\n" + << " --n-contrasts Q shared named contrast rows (default 0)\n\n" << "Conditionally required:\n" << " --extra-covariates PATH row-major rows x n float32 matrix; required\n" << " when any task uses extra covariates\n\n" @@ -174,11 +177,14 @@ void reject_unknown_options(const Arguments& arguments) { "--tasks", "--extra-offsets", "--extra-indices", + "--contrast-matrix", + "--contrast-metadata", "--out-dir", "--n-samples", "--n-base-covariates", "--n-phenotype-rows", "--n-extra-covariate-rows", + "--n-contrasts", "--block-size", "--threads", "--fixed-effect-test", @@ -230,6 +236,8 @@ int main(int argc, char** argv) { paths.tasks_tsv = require_value(arguments, "--tasks"); paths.extra_offsets_i64 = require_value(arguments, "--extra-offsets"); paths.extra_indices_i32 = require_value(arguments, "--extra-indices"); + paths.contrast_matrix_f64 = optional_value(arguments, "--contrast-matrix"); + paths.contrast_metadata_tsv = optional_value(arguments, "--contrast-metadata"); paths.output_directory = require_value(arguments, "--out-dir"); BatchDimensions dimensions; @@ -244,6 +252,19 @@ int main(int argc, char** argv) { dimensions.extra_covariate_row_count = parse_size( require_value(arguments, "--n-extra-covariate-rows"), "--n-extra-covariate-rows", true); + dimensions.contrast_count = parse_size( + optional_value(arguments, "--n-contrasts", "0"), + "--n-contrasts", true); + const bool any_contrast_path = !paths.contrast_matrix_f64.empty() || + !paths.contrast_metadata_tsv.empty(); + const bool all_contrast_paths = !paths.contrast_matrix_f64.empty() && + !paths.contrast_metadata_tsv.empty(); + if ((dimensions.contrast_count == 0 && any_contrast_path) || + (dimensions.contrast_count > 0 && !all_contrast_paths)) { + throw std::invalid_argument( + "--contrast-matrix and --contrast-metadata are both required " + "when --n-contrasts is positive, and must be omitted otherwise"); + } BatchOptions options; options.resume = arguments.flags.count("--resume") != 0U; diff --git a/tests/python/test_spectra_reml_cli.py b/tests/python/test_spectra_reml_cli.py index e069e8f..45023d1 100644 --- a/tests/python/test_spectra_reml_cli.py +++ b/tests/python/test_spectra_reml_cli.py @@ -21,7 +21,7 @@ SPEC.loader.exec_module(cli) class SpectraRemlCliTests(unittest.TestCase): - def fixture(self, root: Path) -> tuple[Path, dict]: + def fixture(self, root: Path, with_contrasts: bool = False) -> tuple[Path, dict]: n, p0 = 3, 2 np.zeros(n * (n + 1) // 2, dtype=" None: + def write_blocks(root: Path, contrast_count: int = 0) -> None: root.mkdir(parents=True, exist_ok=True) + first_contrast_offset = 0 if contrast_count else -1 + second_contrast_offset = contrast_count if contrast_count else -1 (root / "block_000000.summary.tsv").write_text( "\t".join(cli.SUMMARY_HEADER) + "\n" "0\ttrait_a\tconverged\t2\t0\t0\t0\t1\t1\t0.5\t-1\t4\t8\t1e-8" - "\tsatterthwaite\tok\t0\t0\tnan\tnan\tnan\t\t\n" + "\tsatterthwaite\tok\t0\t{}\t0\tnan\tnan\tnan\t\t\n" "1\ttrait_b\tconverged_boundary\t3\t1\t2\t3\t0\t1\t0\t-2\t3\t6\t1e-9" - "\tsatterthwaite\tboundary_conditional\t2\t1\t8\t2\t0.2\t\t\n", + "\tsatterthwaite\tboundary_conditional\t2\t{}\t1\t8\t2\t0.2\t\t\n".format( + first_contrast_offset, second_contrast_offset + ), encoding="utf-8", ) np.asarray([1, 2, 3, 4, 5], dtype=" None: + with tempfile.TemporaryDirectory() as temporary: + root = Path(temporary) + _, manifest = self.fixture(root, with_contrasts=True) + self.assertEqual(manifest["dimensions"]["contrast_count"], 2) + args = cli.build_parser().parse_args( + [ + "run", "--manifest", str(root / "manifest.json"), + "--engine", "spectra_reml", "--fixed-effect-test", + "kenward-roger", + ] + ) + command = cli.build_engine_command(Path("spectra_reml"), manifest, args) + self.assertEqual(command.count("--contrast-matrix"), 1) + self.assertEqual(command.count("--contrast-metadata"), 1) + self.assertEqual(command[command.index("--n-contrasts") + 1], "2") + self.assertEqual(command[command.index("--fixed-effect-test") + 1], "kenward-roger") + + def test_new_contract_rejects_missing_contrast_fields(self) -> None: + with tempfile.TemporaryDirectory() as temporary: + root = Path(temporary) + manifest_path, _ = self.fixture(root) + manifest = cli.read_json(manifest_path) + del manifest["dimensions"]["contrast_count"] + del manifest["paths"]["contrast_matrix_f64"] + del manifest["paths"]["contrast_metadata_tsv"] + cli.atomic_json(manifest_path, manifest) + with self.assertRaises(RuntimeError): + cli.validate_manifest(manifest_path) + + def test_both_fixed_effect_methods_are_forwarded(self) -> None: + with tempfile.TemporaryDirectory() as temporary: + root = Path(temporary) + _, manifest = self.fixture(root, with_contrasts=True) + for method in ("satterthwaite", "kenward-roger"): + args = cli.build_parser().parse_args( + [ + "run", "--manifest", str(root / "manifest.json"), + "--engine", "spectra_reml", "--fixed-effect-test", method, + ] + ) + command = cli.build_engine_command(Path("spectra_reml"), manifest, args) + self.assertEqual( + command[command.index("--fixed-effect-test") + 1], method + ) + + def test_contrast_finalize_contract(self) -> None: + with tempfile.TemporaryDirectory() as temporary: + root = Path(temporary) + manifest_path, manifest = self.fixture(root, with_contrasts=True) + engine = root / "spectra_reml" + engine.write_bytes(b"synthetic engine") + args = cli.build_parser().parse_args( + [ + "run", "--manifest", str(manifest_path), "--engine", + str(engine), "--resume", "--fixed-effect-test", + "satterthwaite", + ] + ) + with mock.patch.object( + cli.subprocess, "run", return_value=SimpleNamespace(returncode=0) + ): + cli.run_engine(manifest, args) + output = Path(manifest["paths"]["output_directory"]) + self.write_blocks(output, contrast_count=2) + result = root / "contrast_results.tsv.gz" + cli.finalize(manifest, result) + with gzip.open(result, "rt", encoding="utf-8", newline="") as handle: + rows = list(csv.DictReader(handle, delimiter="\t")) + self.assertEqual(json.loads(rows[0]["contrast_estimate_json"]), [1.0, 2.0]) + self.assertEqual(json.loads(rows[1]["contrast_estimate_json"]), [3.0, 4.0]) + self.assertEqual(json.loads(rows[0]["contrast_numerator_df_json"]), [1.0, 1.0]) + def test_signature_recovery_and_finalize(self) -> None: with tempfile.TemporaryDirectory() as temporary: root = Path(temporary)