This now should me mentally engaging.

Author: stela2502

docs/06-flow.md

@@ -95,7 +95,7 @@ for i in range(2):
 Print every second sample name:
 
 ```python
-samples = list(data["expression"].columns)
+samples = list(data["samples"].index)
 
 for i in range(0, len(samples), 2):
     print(samples[i])
@@ -128,7 +128,7 @@ expr = data["expression"]
 
 gene_means = []
 
-for gene in range(len(expr)):
+for gene in range(expr.shape[0]):
     m = expr[gene].mean()
     gene_means.append(m)
 
@@ -201,13 +201,15 @@ while w <= 5:
 This shows how `while` can stop early when a condition is met.
 
 ```python
+cluster1_samples = []
 i = 0
-samples = list(data["samples"].index)
 
-while i < len(samples) and data["samples"].loc[samples[i], "cluster"] != 1:
+while i+1 < len(data["samples"]):
+    if data["samples"]["cluster"].iloc[i] == 1:
+        cluster1_samples.append(data["samples"].index[i])
     i = i + 1
 
-samples[i]
+cluster1_samples
 ```
 
 ---

docs/07-functions.md

@@ -84,7 +84,7 @@ def get_gene(data, gene):
     """
     if gene not in data["genes"].index:
         raise ValueError(
-            f"Unknown gene '{gene}'. Expected one of: {list(data["genes"].index)}"
+            f"Unknown gene '{gene}'. Expected one of: {list(data['genes'].index)}"
         )
 
     g_idx = data["genes"].index.get_loc(gene)
@@ -138,7 +138,7 @@ in the sample table equals `value`.
    samples_in(data, col_name, value ):
       """
       Returns a list of samples wher <column_name> equals <value>
-      """"
+      """
       samples = data["samples"]
 
       if col_name not in samples.columns:

docs/09-performance.md

@@ -31,6 +31,16 @@ hspc_data = pd.read_csv(
 )
 
 hspc_data
+
+## and now force that into our data type
+def from_pd(expr):
+    genes = pd.DataFrame(index=expr.index)
+    samples = pd.DataFrame(index=expr.columns)
+    return {"expression": np.array(expr), "genes": genes, "samples": samples}
+
+hspc_data = from_pd(hspc_data )
+hspc_data
+
 ```
 
 ---
@@ -79,9 +89,9 @@ def euclid_distance_by_hand(v1, v2):
     return math.sqrt(s)
 
 
-def distance_matrix_df(df, func):
-    X = np.asarray(df)
-    names = df.index
+def distance_matrix_df(data, func):
+    X = data['expression']
+    names = data['genes'].index
     n = X.shape[0]
 
     D = np.zeros((n, n))
@@ -94,10 +104,27 @@ def distance_matrix_df(df, func):
     return pd.DataFrame(D, index=names, columns=names)
 ```
 
-Time it:
+Time it - or better time a subset of the real data:
+
+```python
+## actually first get us a subset of our data:
+def subset_genes(data, gene_idx):
+    gene_idx = np.asarray(gene_idx)
+
+    X2 = data["expression"][ gene_idx,:]
+    genes2 = data["genes"].iloc[gene_idx].copy()
+
+    # samples unchanged
+    samples2 = data["samples"].copy()
+
+    return {"expression": X2, "genes": genes2, "samples": samples2}
+hspc_data_tiny = subset_genes(  hspc_data , np.arange(200) )
+check_data_model( hspc_data_tiny )
+hspc_data_tiny.shape
+```
 
 ```python
-(_, t_slow) = timed(distance_matrix_df, hspc_data.iloc[range(200)], euclid_distance_by_hand )# first 200 genes of the whole data only
+(_, t_slow) = timed(distance_matrix_df, hspc_data_tiny, euclid_distance_by_hand )# first 200 genes of the whole data only
 print("slow:", t_slow)
 ```
 
@@ -119,7 +146,7 @@ def euclid_distance_vec(v1, v2):
 Time it:
 
 ```python
-(_, t_faster) = timed(distance_matrix_df, hspc_data.iloc[range(200)], euclid_distance_vec )
+(_, t_faster) = timed(distance_matrix_df, hspc_data_tiny , euclid_distance_vec )
 print("faster:", t_faster)
 ```
 
@@ -132,15 +159,13 @@ But there are also specialized function that operate on the full matrix:
 ```python
 from scipy.spatial.distance import pdist, squareform
 
-def vectorized_dist_mat(mat):
+def vectorized_dist(data):
     """
-    mat: pandas DataFrame or numpy array (rows = observations, cols = features)
-    returns: full (n x n) Euclidean distance matrix as numpy array
+    data: a dist with "expression" - a numpy array (rows = features, cols = observations)
+    returns: full (nrow x nrow) Euclidean distance matrix as numpy array
     """
-    X = np.asarray(mat, dtype=float)
-
     # upper triangle (condensed form)
-    upper = pdist(X, metric="euclidean")
+    upper = pdist(data['expression'], metric="euclidean")
 
     # convert to full symmetric matrix
     D = squareform(upper)
@@ -164,70 +189,25 @@ If not 100% sure you know all about the possible libraries you could use I recom
 
 ---
 
-# What about pandas apply?
-
-In R, `apply()` can be a speed trick.  
-In Python, `DataFrame.apply()` usually still runs a Python function once per row/column.
-
-That means it often behaves like a loop (and can be slow). But eucledian distance is probably the worst example for apply as it works on one row/column only.
-So for this we switch to a simple mean calculation.
-
-```python
-def mean_vec(vec1):
-    return vec1.mean() # very fast here
-```
-
-```python
-(_, t_vec) = timed( hspc_data.apply, mean_vec , axis=1)
-print("apply time:", t_vec)
-```
-
-This is rather fast, but using the inbuild pandas mean function is even faster:
-
-
-```python
-(_, t_vec) = timed( hspc_data.mean,  axis=1)
-print("pandas time:", t_vec)
-```
-
-If we would have our data as a numpy array instead we could get this even faster using numpy's vectorization:
-
-
-```python
-arr = np.array(hspc_data)
-(c, t_fast) = timed(arr.mean, axis=1 )
-print("numpy:", t_fast)
-```
-
-**Take Home** If numpy has a function for your problem use that!
-
----
 
 # Interpreting the result
 
 Usually:
 
 - vectorized version is very fast
-- apply is much slower
-
-Why? Because `apply` calls Python code repeatedly.
-
----
-
-# When is apply OK?
-
-`apply` can be reasonable when:
+- pure python as well as manual for loops are much slower
 
-- there is no clean vectorized solution
-- the dataset is small
-- readability matters more than speed
+Why? Because even a for loop calls Python code repeatedly whereas the vectorized function (c or fortran) works on one big memory block.
 
-But for large biological matrices (genes x samples), prefer vectorization.
 
 ---
 
 # Exercise
 
-Just keep the ``timed`` function and apply it later on whenever you like ;-)
+Take the function ``zscore_rows`` and convert it from using a numpy ndarray to using our own data structure.
+While doing that change tha action to modifying the data in place. 
+
+**Note:** Mutable objects (like lists, dictionaries, and arrays) can be changed inside a function, while immutable objects (like numbers and strings) cannot. Think of it like "Small objects like numbers or strings can be copied, but putatively large ones like matrices or dictionaries should not be copied". 
+
 
 In the next section, we will apply these ideas to real expression data: selecting variable genes and scaling (z-scores).

docs/10-variable-genes.md

@@ -52,25 +52,30 @@ So we often focus on the **top N most variable genes**.
 # Select the top variable genes
 
 ```python
-def get_top_variable_genes(mat, top_n=500):
-    variances = mat.var(axis=1, ddof=1)          # variance per row (gene)
-    top_genes = variances.nlargest(top_n).index  # gene names of the largest variances
-    return mat.loc[top_genes]
+def get_top_variable_genes(data, top_n=500):
+    mat = data['expression']
+    var = mat.var(axis=1, ddof=1)          # variance per row (gene)
+    
+    # 2) get sorted gene indices (small → large variance)
+    sorted_ids = np.argsort(var)
+
+    # return a subset (bottom top_n)
+    return subset_genes(data, sorted_ids[-top_n:])
 ```
 
 Use it:
 
 ```python
 hspc_var = get_top_variable_genes(hspc_data, top_n=500)
-hspc_var.shape
+hspc_var['expression'].shape
 ```
 
 ---
 
 # Check the value range
 
 ```python
-hspc_var.values.min(), hspc_var.values.max()
+hspc_var['expression'].min(), hspc_var['expression'].max()
 ```
 
 ---
@@ -100,16 +105,17 @@ This makes genes comparable even if they have very different expression ranges.
 
 ```python
 def zscore_rows(mat):
-    m = mat.mean(axis=1)
-    s = mat.std(axis=1, ddof=1)
-    return mat.sub(m, axis=0).div(s, axis=0)
+    m = np.mean(mat, axis=1, keepdims=True)
+    s = np.std(mat, axis=1, ddof=1, keepdims=True)
+    return (mat - m) / s
 ```
 
 Run it:
 
 ```python
-hspc_zs = zscore_rows(hspc_var)
-hspc_zs.shape
+hspc_zs = hspc_data_tiny.copy()
+hspc_zs['expression'] = zscore_rows(hspc_data_tiny['expression'])
+hspc_zs['expression'].shape
 ```
 
 ---
@@ -133,18 +139,30 @@ print(hspc_zs.std(axis=1, ddof=1).head())
 
 # Visual check: boxplots before and after scaling
 
+There is a boxplot function in pandas, but our dict strores the data as numpy array.
+The easiest here is to define one more function that actually converts out dict to a pandas DatFrame - the reverse of our from_df function earlier:
+
+```python
+def expression_df(data):
+    return pd.DataFrame(
+        data["expression"],
+        columns=data["samples"].index,
+        index=data["genes"].index
+    )
+```
+
 Before:
 
 ```python
-hspc_var.boxplot(rot=90)
+expression_df(hspc_var).T.boxplot(rot=90)
 plt.title("Before scaling (raw values)")
 plt.show()
 ```
 
 After:
 
 ```python
-hspc_zs.boxplot(rot=90)
+expression_df(hspc_zs).T.boxplot(rot=90)
 plt.title("After scaling (z-scores)")
 plt.show()
 ```
@@ -153,9 +171,8 @@ plt.show()
 
 # Exercise
 
-1. Select the top 200 variable genes into `hspc_var200`.
-2. Z-score them into `hspc_zs200`.
-3. Verify that row means are ~0 and row standard deviations are ~1.
+We used ``expression_df(hspc_var).T`` there, but we likely also need a transform for our own data.
+Implement a ``def transform()`` that returns a transformed data structure. 
 
 ---