Repository navigation
Fix/nmf kl floor - #314
Merged
Merged
Fix/nmf kl floor#314
Conversation
The KL updates in nmf_cpu.py and nmf_gpu.py had no lower bound on W, H or W @ H. A multiplicative update cannot move an entry away from exactly zero, so zeros from the NNDSVD initialization (about half of W and H) or from underflow stayed zero for the whole run. The fit then stalls, and the convergence test can pass early at a worse solution. On pipeline-style input (values raised to at least 1e-4, as in pnmf()), the previous code with nmf_init="nndsvd" stopped at the first allowed check (20,001 iterations) with a KL loss 5-8% higher than with the floor. With the default random start the loss was 0.4-0.6% higher. The original CPU engine (nimfa) had an equivalent floor; it was lost when the CPU loop was replaced by a copy of the GPU loop in 441159a. Changes: - Clamp W @ H, both denominators, and W and H after each update to at least 1e-16 (CPU and GPU, beta = 1). - Compute the KL loss with torch.xlogy. V * log(V / WH) gave NaN wherever V == 0. pnmf() raises V to at least 1e-4, so this affected only direct use of the NMF classes. - Raise FloatingPointError in decipher_signatures() if a replicate's W or H is not finite, instead of passing it to clustering. If a cell of W @ H reached zero, the old update produced inf/NaN with no error (on raw counts with zeros, at iteration 208 in the test below). - Add tests/test_nmf_floor.py: 4 tests, all fail on the previous code and pass now; the GPU variant runs only when CUDA is available.
…ssue 1) parallel_runs() split NMF_replicates into batches of batch_size for both backends, but the CPU pnmf() runs a single replicate per task. With NMF_replicates=10 and batch_size=5, only 2 CPU replicates ran, and the remaining slots of Wall/Hall stayed zero and entered clustering. - Add replicate_batches(): one replicate per task on CPU, batches of batch_size on GPU (unchanged). - Raise RuntimeError if the number of returned replicates differs from NMF_replicates. - Add tests/test_replicate_batches.py. Default settings (batch_size=1) are unaffected.
…ssue 2) sigProfilerExtractor() accepted nnls_add_penalty, nnls_remove_penalty, initial_remove_penalty and collapse_to_SBS96 and recorded them in JOB_METADATA.txt, but never passed them to decomp.spa_analyze(), so the final assignment always used SigProfilerAssignment's defaults. - Pass all four to spa_analyze(). - Remove the unused add_penalty/remove_penalty locals. - Add tests/test_assignment_forwarding.py. SigProfilerAssignment's defaults equal Extractor's (0.05, 0.01, 0.05, True), so runs with default settings are unaffected.
startProcess and endProcess were capped with min(current value, number of samples) inside the mutation-context loop, reusing the value left by the previous context. A sparse early context (for example DBS with 3 samples) therefore limited every later context to ranks 1-3. - Recompute the range for each context from minimum_signatures and maximum_signatures. - Record the rank range evaluated for each context in JOB_METADATA.txt. Only multi-context (VCF) runs were affected. Not run end to end here, because that needs a VCF input and an installed reference genome.
get_normalization_cutoff() fitted sklearn's GaussianMixture without a random_state, so the cutoff depended on NumPy's global random state rather than on the seed recorded in Seeds.txt. - Add a random_state argument to get_normalization_cutoff() (default None, so other callers are unchanged). - In sigProfilerExtractor(), derive it from the root seed with SeedSequence(seed).generate_state(1), which also maps large seeds into GaussianMixture's accepted range. - Add tests/test_normalization_cutoff.py. On the test data the unseeded fit gave the same cutoff for 40 different global seeds, so this makes reproducibility guaranteed rather than fixing an observed difference. The cutoff is already written to JOB_METADATA.txt.
cluster_converge_outerloop() wrapped silhouette_samples() in a bare except and set every coefficient to 1.0 on any error, so a failed or undefined clustering looked perfectly stable. An unknown distance name also ended there (NameError) and became 1.0. Add silhouette_stability(): - rank 1 (one cluster): 1.0 for every signature, as before, since the silhouette is undefined there; - one replicate per signature: raise ValueError asking for at least 2 NMF replicates; - any other error: raise RuntimeError with the rank and the original error. Add tests/test_silhouette_stability.py.
The seed was kept as a NumPy array (np.array(seeds["Seed"])) and later
converted with int(). A multi-row seeds file failed only deep inside the
run, and with NumPy >= 2.5 even a one-row file fails, because int() of
a one-element array is now a TypeError ("only 0-dimensional arrays can
be converted to Python scalars").
- Add read_seed_file(): read the "Seed" column, require exactly one
seed, return it as an int; raise ValueError otherwise.
- Store the random seed as an int too.
- Document the one-seed contract in the docstring.
- Add tests/test_seed_file.py.
…ssue 6) decipher_signatures() filled finalgenomeErrors and finalgenomesReconstructed (mutation types x samples x replicates, float64) for every rank, but nothing reads them: export_information() never uses those positions of loopResults. At 100 replicates this is about 1.2 GB per array for SBS1536 x 1,000 samples and 3.7 GB for SBS4608 x 1,000 samples. Remove both arrays, drop them from the decipher_signatures() return value and from loopResults, and update the converge_information index in export_information(). Checked with a small end-to-end matrix run (ranks 2-4, 5 replicates, fixed seed) before and after: all 30 text/CSV outputs are identical.
The Python API defaulted to maximum_signatures=25 while the CLI parser defaulted to 10, and the API docstring said 10. The README and the input documentation state 25, which the API has used since 2020. - Add sigpro.DEFAULT_MAXIMUM_SIGNATURES = 25 and use it for both the API default and the CLI default and help text, as is already done for DEFAULT_COSMIC_VERSION. - Fix the docstring (10 -> 25). - Add tests/test_default_parity.py. This changes the CLI default from 10 to 25: CLI runs that do not pass --maximum_signatures now search up to 25 signatures.
…iew issue 8)
The documented option "custom" had no branch in normalize_samples(); it
fell into a try/int(norm)/except-pass block, so it silently meant no
normalization. Only a numeric value worked as a manual cutoff, which was
not documented. The docstring also gave the wrong default ("log2") and a
non-existent option ("no_normalization").
- Add check_matrix_normalization(): accept "gmm", "100X", "log2",
"none", or a positive whole number (int or numeric string) used as a
manual cutoff; raise ValueError for anything else, including "custom".
- Call it at the start of sigProfilerExtractor() and in
normalize_samples(); remove the bare except.
- Document the numeric cutoff in the docstring, CLI help, README and
docs, and remove "custom".
- Add tests/test_matrix_normalization.py.
Runs that passed "custom" (or another unrecognised string) now stop with
an error instead of running unnormalized.
Input preprocessing called data.dropna(axis=1), which removed any sample with even one missing value, and then removed all-zero samples, both without telling the user. Add check_input_samples(): - a sample with some missing values raises ValueError naming the sample and the mutation types that are missing; - columns that are entirely empty (for example from a trailing tab) and samples with no mutations are still removed, but are now listed on screen and in JOB_METADATA.txt. Add tests/test_input_samples.py. Inputs with partly missing samples, which used to run with those samples dropped, now stop with an error.
nnmf_cpu() and nnmf_gpu() converted the input with torch.from_numpy(genomes).float(), so with precision="double" the matrix was rounded to float32 before the NMF class cast it back to float64. Add genomes_tensor(), which builds the tensor directly in the requested precision (float64 for "double", float32 otherwise), and use it on both paths. Single precision is unchanged. Add tests/test_precision.py.
`mtypes[0] == "96" or "288" or "384" or "1536" or "4608"` is always true, so a matrix with an unknown number of channels (for example 100) was labelled "SBS100" instead of "CH100". - Use set membership for the SBS channel counts. - Keep "CH..." labels as they are when building mutation_type, so the corrected label does not become "CHCH100". - Add tests/test_context_labels.py, which rejects string constants used as `or` operands in sigpro.py. Checked with a 100-channel matrix run before and after: the only difference is JOB_METADATA.txt, which now records context_types CH100 instead of SBS100. The loop already stripped "SBS" later, so the output folder was CH100 in both cases and all other outputs are identical.
Unavailable or empty VCF contexts were skipped with `continue`, and the run then always wrote "Analysis of mutational signatures completed successfully" to JOB_METADATA.txt, even if every context was skipped. - Track analyzed and skipped contexts. - List skipped contexts, with the reason, on screen and in JOB_METADATA.txt. - If no context was analyzed, write the failure to JOB_METADATA.txt and raise RuntimeError instead of reporting success. Matrix inputs are unaffected (a small end-to-end matrix run gives identical outputs). The VCF path was not run here, because it needs VCF input and an installed reference genome.
Rerunning into an existing output folder wrote over files one by one, so rank folders from an earlier, wider run (for example ranks 6-10 after a rerun with ranks 1-5) stayed next to the new results and could be mistaken for them. Before a context's results are written, move_previous_results() moves an existing, non-empty <output>/<context> folder to <output>/<context>_previous_<timestamp>. Nothing is deleted. The move is reported on screen and in JOB_METADATA.txt. Rerunning into the same folder still works as before; the new results are complete and identical, and the old ones are kept beside them. Folders of contexts that the new run does not analyze are left as they are. Add tests/test_output_directory.py.
export_information() wrote the *_NMF_Activities.txt files from exposureAvg.astype(int), which truncates towards zero. In a normal run this had no effect: the exposures are refitted with SigProfilerAssignment's fit_signatures_pool(), which already returns np.round()ed values, so astype(int) received whole numbers. (Checked on the example matrix: all activity files are identical before and after.) It would only matter if that refit were skipped, when the unrounded NMF averages would be truncated. Add integer_activities(), which rounds to the nearest integer, and use it for the export. The file format is unchanged. Writing full-precision activities instead would change a file that users and SigProfilerPlotting read, so that is left as a separate decision. Add tests/test_activity_export.py.
…ue 16) For BEDPE input the SV matrices were written to os.path.join(os.path.split(input_data)[0], "SV_Matrices"). For "/data/sv" that is /data/SV_Matrices, but for "/data/sv/" it is /data/sv/SV_Matrices, so the same folder written two ways gave two output locations. Add sv_matrix_directory(), which normalizes the path first and always returns SV_Matrices next to the input folder (/data/SV_Matrices in both cases, as without a trailing slash). Users who passed a trailing slash will now find SV_Matrices next to the input folder instead of inside it. Add tests/test_sv_output_path.py.
Importing SigProfilerExtractor.subroutines called
warnings.filterwarnings("ignore") and
multiprocessing.set_start_method("spawn", force=True), which changed
the warning filters and the start method of the whole host program
(notebooks, pipelines that embed the package).
- Remove the module-level filterwarnings("ignore").
- Add SPAWN = multiprocessing.get_context("spawn") and create every
SigProfilerExtractor pool from it (3 places), instead of setting the
global start method.
- Add tests/test_import_side_effects.py (fails on the previous code).
Checked with a full matrix run including SigProfilerAssignment
(cosmic_version=3.5): all result files are identical before and after,
and no warnings are printed in either case. Only the two
*_Signature_Assignment_log.txt files differ, in the order of their
per-sample blocks; they also differ between two runs of the previous
code, because parallel workers write them.
Note: sigProfilerPlotting still calls warnings.filterwarnings("ignore")
when imported (as do SigProfilerMatrixGenerator's SV and CNV scripts),
so importing SigProfilerExtractor still hides warnings through that
dependency. That needs the same fix in those packages.
… issue 18) The per-replicate fit diagnostics (mean L1, L2, KL and related metrics, plus the iteration count) were computed differently on each path: - CPU (nnmf_cpu): the bootstrapped, normalized matrix that NMF fitted against W @ H as returned by NMF; - GPU (a separate branch in decipher_signatures): the original input matrix against W @ H after W was normalized and H denormalized. So CPU and GPU values were not comparable. - Add replicate_diagnostics(), with the CPU definition, and use it in nnmf_cpu() and in the GPU branch of pnmf(), on the bootstrapped, normalized matrix and the raw W and H of each replicate. - Remove the GPU-only recomputation in decipher_signatures(); both paths now return (W, H, diagnostics) from parallel_runs(). - Add tests/test_replicate_diagnostics.py. CPU results are unchanged (a small end-to-end matrix run gives identical outputs, including the NMF convergence files). GPU diagnostics now follow the CPU definition; checked on CPU with a stand-in for nnmf_gpu that returns its exact output structure, since no GPU was available. With batch_size > 1 the GPU path still reports max_iterations as the iteration count, because the batched NMF does not track convergence per replicate.
…s CPU With batch_size > 1 the GPU path did not follow the CPU logic: - all replicates of a batch drew their bootstrap resamples from one shared generator (the other replicates' generators were never used); - they drew their random W/H start from one shared generator, in a different order from nmf_cpu; - H of every replicate was denormalized with the sample totals of the last replicate of the batch; - a batch never stopped at convergence and reported max_iterations. With batch_size = 1 (the default) CPU and GPU already matched. Now CPU and GPU share one pnmf() routine: - prepare_replicate(): bootstrap with the replicate's own generator, 1e-4 floor and normalization, returning the replicate's own totals; - split_generators(): each task gets the generator pairs of exactly its own replicates, so replicate k uses pair k whatever the batch size; - nmf_gpu.NMF accepts one generator per replicate and draws each replicate's W and then H from it, as nmf_cpu does; the NNDSVD variants use the replicate's own generator too; - in a batch, each replicate stops by the same rule as nmf_cpu (every test_conv iterations, relative loss change below the tolerance after min_iterations) and is then frozen; replicate_conv reports its own iteration count, which feeds its diagnostics. Checks: with batch_size = 1, GPU outputs of a small end-to-end run are identical to the previous commit, and CPU outputs are identical to before. On a GPU, each replicate's bootstrap and starting W and H are bit-identical on CPU, GPU alone and GPU in a batch of 4. Final W and H differ only by rounding (largest relative difference 5e-5 GPU vs CPU, 2e-4 batch vs alone), because CPU, single and batched GPU matrix products round differently. Near the tolerance this can move a replicate's stop by one convergence check; with the default tolerance (1e-15) it moves more often, but a batch-4 run gave the same stability table and signatures (lowest cosine to batch 1: 0.99999916).
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Main fix: floor in the NMF updates
The KL multiplicative updates had no lower bound, so values that reached zero (about half of W and H after an NNDSVD start) stayed zero. The fit then stalled and stopped early at a worse solution, and a zero in W @ H could give NaN with no error. W, H and W @ H are now kept at or above 1e-16, the KL loss uses
torch.xlogy, and non-finite replicates raise an error. The CPU code had such a floor until 2020 (nimfa).Other fixes
batch_size > 1collapse_to_SBS96are passed to SigProfilerAssignmentmaximum_signaturesdefault (25) for API and CLImatrix_normalizationvalidatedprecision="double"stays float64CH<n>Also: with
batch_size > 1, each GPU replicate now has its own generators, sample totals and convergence stop, as on CPU.Behaviour changes
--maximum_signatures: 10 -> 25.batch_size > 1: each replicate in a batch now behaves exactly as it would alone (as withbatch_size = 1)