Add reconstruction error calculation with filtered spectra

.github/workflows/python-tests.yml

@@ -19,10 +19,35 @@ jobs:
           python-version: "3.10"
           cache: pip
 
+      - name: Read synthetic regression source
+        id: synthetic-regression
+        run: |
+          python - <<'PY'
+          import json
+          import os
+          from pathlib import Path
+
+          spec = json.loads(
+              Path("tests/fixtures/synthetic_cnmf_k3/run_spec.json").read_text()
+          )
+          with open(os.environ["GITHUB_OUTPUT"], "a", encoding="utf-8") as output:
+              print(f"repository={spec['upstream']['repository']}", file=output)
+              print(f"commit={spec['upstream']['commit']}", file=output)
+          PY
+
+      - name: Check out pinned upstream cNMF
+        uses: actions/checkout@v4
+        with:
+          repository: ${{ steps.synthetic-regression.outputs.repository }}
+          ref: ${{ steps.synthetic-regression.outputs.commit }}
+          path: .upstream/dylkot-cnmf
+
       - name: Install dependencies
         run: |
           python -m pip install --upgrade pip
           python -m pip install -e ".[gpu,test]"
 
       - name: Run tests
+        env:
+          CNMF_UPSTREAM_SOURCE: ${{ github.workspace }}/.upstream/dylkot-cnmf
         run: bash run_tests.sh

src/cnmf/cnmf.py

@@ -935,6 +935,8 @@ def consensus(self, k, density_threshold=0.5, local_neighborhood_size=0.30, show
                     columns = ['stats'])
 
             return(consensus_stats)           
+
+        density_filtered_refit_usages = rf_usages
 
         # Re-order usage by total contribution
         norm_usages = rf_usages.div(rf_usages.sum(axis=1), axis=0)      
@@ -1095,6 +1097,8 @@ def consensus(self, k, density_threshold=0.5, local_neighborhood_size=0.30, show
             l2_spectra=l2_spectra,
             local_density=local_density,
             kmeans_cluster_labels=kmeans_cluster_labels,
+            normalized_counts=norm_counts.X,
+            refit_usages=density_filtered_refit_usages.values,
         )
 
 

src/cnmf/qc.py

@@ -138,6 +138,7 @@ def get_r2(X, U, H, log2p1=False, scale=False):
 
     return r2, sse, tss
 
+
 def plot_heatmap(df, outfile, scale_rows=False):
     # df: rows=sources, cols=conditions
     if scale_rows:
@@ -210,11 +211,20 @@ def run_qc(
     l2_spectra,
     local_density,
     kmeans_cluster_labels,
+    normalized_counts,
+    refit_usages,
 ):
     """Calculate and save QC using state already produced by ``consensus``."""
     # local_density contains all spectra, while l2_spectra has already been density-filtered.
     nruns = len(local_density)
     nruns_per_gep = nruns/k
+    if sp.issparse(normalized_counts):
+        normalized_counts = np.asarray(normalized_counts.todense())
+    run_r2, run_sse, run_tss = get_r2(
+        normalized_counts,
+        refit_usages,
+        median_spectra.values,
+    )
 
     #-----------------------------------------------------------------------
     # Calculate how many iterations feed into each GEP    
@@ -258,6 +268,9 @@ def run_qc(
     annotation['run_davies_bouldin'] = davies_bouldin_score(l2_spectra.values, kmeans_cluster_labels)
     annotation['run_median_density'] = np.median(local_density.iloc[:, 0])
     annotation['run_mean_density'] = np.mean(local_density.iloc[:, 0])
+    annotation['run_r2'] = run_r2
+    annotation['run_sse'] = run_sse
+    annotation['run_tss'] = run_tss
 
     annotation.index = range(1, len(annotation)+1)
     annotation = annotation.loc[reorder.index,:]