<jupyter_start><jupyter_text>NESM Python Part 4 - Advanced Topics- Deep learning with Tensorflow- Our image analysis pipeline at a glance- Dask for out of memory computing- Classical machine learning with `scikit-learn`<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
from mpl_interactions import hyperslicer
%matplotlib widget<jupyter_output><empty_output><jupyter_text>Dask for out of memory computing<jupyter_code>memory = 8e9 #8GB 
pixels = 1024*1024
bytes_per_pix = 2 #16 bit unsigned ints
memory/(pixels*bytes_per_pix) #images you can have in memory<jupyter_output><empty_output><jupyter_text>That seems like a lot but that corresponds to less than(20 Time points) x (10 Positions) x (4 Channels) x (5 z-slices) = 4000 Images**Enter Dask Array**<jupyter_code>20*10*4*5
import dask.array as da # To use dask you can call da.xxx when you would have called np.xxx
#impossible_arr = np.random.random((10000,1024,1024))
darr = da.random.random((10000, 1024, 1024))
darr<jupyter_output><empty_output><jupyter_text>What is this thing? Dask arrays function numpy arrays but they are *lazy*. This means that rather than storing all the information of the array in memory, dask keeps track of a *graph* that describes all the computation necessary to generate the array.We can continue to do operations on this array and they will all be lazy. Again this means dask just adds computations to the graph.<jupyter_code>darr.mean(0)
(darr - darr.min())/(darr.std())
# launch a client from dask-labextension
# scale to 3 cores
# Or, the oldschool way
# from dask.distributed import Client
# client = Client()
# client.cluster<jupyter_output><empty_output><jupyter_text>When we actually want to evaluate an array, we call `.compute()` on it. You should watch the dask dashboard (by following the link above). Also beware that you're computer is about to start running aggressively so if you have low battery<jupyter_code>out = darr.mean(0).compute()
out<jupyter_output><empty_output><jupyter_text>Loading data with dask<jupyter_code>import os
# I have had a hard time with dask's imread in the past, here it works great, your mileage may vary
# Also it appears to need the absolute path to the data
cho_arr = da.image.imread(os.getcwd()+'/Fluo-N3DH-CHO/01/*')
cho_arr<jupyter_output><empty_output><jupyter_text>Dask + XarrayXarray supports any<jupyter_code>coords = {'T':9.5*np.arange(cho_arr.shape[0]), 
          'Z':1.0*np.arange(cho_arr.shape[1]),
          'Y':0.202*np.arange(cho_arr.shape[2]),
          'X':0.202*np.arange(cho_arr.shape[3])}
x_data = xr.DataArray(cho_arr, dims=coords.keys(), coords=coords)
x_data
plt.figure()
controls = hyperslicer(data)<jupyter_output><empty_output><jupyter_text>**Note about Dask**- One of my favorite things about dask is that I can develop on my laptop and run with 4 cores but then move to Harvard's computing cluster and run with much more computing power. Dask scales seamlessly between these two settings.- Dask maintains several different APIs. I'd recommend [this page from their documentation](https://docs.dask.org/en/latest/user-interfaces.html) to see what would work for you. In brief there are high level interfaces:  - Array - for data that is a high dimensional rectangle - Will take you far for large imaging datasets and is likely the easiest to use.  - Dataframe - for tabular data. Array:Numpy :: Dataframe:Pandas  - Bag - more like a database format that implements Map-Reduce type operations.  - Dask-ML - scikit-learn (more on this below) like interface for scaling machine learning tasks.  - There are also lower level interfaces for custom computation  - Delayed - For custom python computation that does not necessarily fit the array paradigm. *Importantly* you set up all your computation and tell dask when evaluate it.  - Futures - *Dynamic* custom computation. Things start running in real time and dask decides when things run by evaluating which computations depend on other computations. This is likely the most powerful and most confusing interface. Example: Return a list of some length and do operations on all elements.  - Keep an eye on [`dask_image`](http://image.dask.org/en/latest/) which is gradually wrapping scipy.ndimage functionality for dask arrays.   PCA on Hyperspectral SRS imaging data**What is PCA?**Principal component analysis (PCA) finds the basis vectors which explain most of the variance in a dataset. Below is a picture from the [wikipedia page](https://en.wikipedia.org/wiki/Principal_component_analysis) which shows the principal components of some correlated 2D data. **What is SRS?**[Stimulated Raman Scattering](https://en.wikipedia.org/wiki/Stimulated_Raman_spectroscopy) (SRS) is an optical imaging technique that probes the vibrational energy levels of different molecules. I'm using it to study cellular metabolism and composition but its generally good for chemical mapping of materials with different vibrational energy levels.The toy dataset below is a spectral scan of two different species of beads. We will use PCA to "discover" how many different species are in the sample and what their spectra look like.<jupyter_code>import io
import requests
# Get the dataset directly from github (33MB)
# Feel free to just watch if you dont want to download
response = requests.get(
    "https://github.com/jrussell25/data-sharing/raw/master/srs_beads.npy"
)
response.raise_for_status()
beads = np.load(io.BytesIO(response.content))
# Define the coordinates for the xarray as a dict of name:array pairs
# Wavenums is the relevant spectroscopic unit in cm^-1
# X,Y = actual dimensions of the images in microns from microscope metadata
coords = {'wavenums':np.linspace(2798.65, 3064.95, beads.shape[0]),
          'X':np.linspace(0, 386.44, 512),
          'Y':np.linspace(0, 386.44,512)}

x_beads = xr.DataArray(beads, dims=coords.keys(), coords=coords)
plt.figure()
ctrls = hyperslicer(x_beads)<jupyter_output><empty_output><jupyter_text>How to do PCA in python?Of course we google it first.But the answer is [scikit-learn](https://scikit-learn.org/stable/). The `scikit-learn` (or `sklearn` for short) user manual is exceptionally good and always contains an example. Googling `sklearn some new method` is always the first thing I do when I want to try out a machine learning method.<jupyter_code>from sklearn.decomposition import PCA
pca = PCA(n_components=10)
# need to do some annoying reshapeing because sklearn expects (N_data, N_features)
# N_data = N_pixels = 512*512 and N_features = 126 spectral points
pcs = pca.fit_transform(beads.reshape(beads.shape[0], -1).T)
#instead of 126 spectral points, we have 10 features corresponding to the first 10 PCs
pcs.shape
#Plot the principal components with 
plt.figure()
plt.plot(x_beads['wavenums'],pca.components_[:5].T + np.arange(5)[None,:])
plt.show()
plt.figure()
plt.plot(pca.explained_variance_ratio_)
# a fun visualization
# Use the weights corresponding to the first 3 PCs as RGB values
rgb = pcs[...,:3].reshape(512, 512,3)
# for matplitlib, we need to normalize each color channel to be between 0 and 1
rgb = (rgb-rgb.min(0).min(0))
rgb = rgb/rgb.max(0).max(0)
plt.figure()
plt.imshow(rgb)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0877095746154912, 'fraction_numerical': 0.02133850630455868, 'mean_word_length': 5.106598984771574, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 37, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '370188', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2314, 'n_tokens_neox': 2137, 'n_words': 855, 'path': 'part4-advanced-topics.ipynb', 'repo_name': 'Hekstra-Lab/nesm-python', 'size': '7217'}