chapter3-figures-1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8ed69048-ecee-4a9d-9a13-a09712ffcae2",
   "metadata": {},
   "source": [
    "# 2023-05-16 d7-glucose uptake"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00819f9f-65c9-4263-8384-947cd1bd1f41",
   "metadata": {},
   "outputs": [],
   "source": [
    "from dask.distributed import Client\n",
    "\n",
    "client = Client(processes=False)\n",
    "client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "486a8122-f8cd-47ef-8004-0c04bd2d0af2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "import re\n",
    "import os\n",
    "from pathlib import Path\n",
    "from aicsimageio.readers import TiffGlobReader\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_interactions import hyperslicer\n",
    "import napari\n",
    "import microutil as mu\n",
    "import scipy.ndimage as ndi\n",
    "from scipy.stats import linregress\n",
    "from scipy.spatial import cKDTree\n",
    "from skimage.morphology import disk\n",
    "from skimage.draw import disk as draw_disk\n",
    "from skimage.registration import phase_cross_correlation\n",
    "from skimage.filters import sobel\n",
    "from tqdm.autonotebook import tqdm\n",
    "from skimage.registration import phase_cross_correlation\n",
    "from srs_tools import BackgroundEstimator\n",
    "\n",
    "%matplotlib widget\n",
    "plt.style.use('paper')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "624cfc7f-0935-4916-9199-ddca185a71c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
    "\n",
    "labels_cmap = plt.cm.viridis.copy()\n",
    "labels_cmap.set_under(alpha=0)\n",
    "\n",
    "fig_path = \"/Users/johnrussell/Documents/figures/2023-06-figs-paper/\"\n",
    "savefigs = False"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa1e6864-25a7-46df-b023-a6cafdc41172",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Load and preprocess data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd9aedcb-b158-4313-8cc5-c39f67f715c8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# def indexer_check(idxr):\n",
    "#     for col in idxr:\n",
    "#         consecutive = np.arange(idxr[col].nunique())\n",
    "#         found = np.sort(idxr[col].unique())\n",
    "#         assert np.allclose(consecutive, found), f\"Found invalid values in column {col}\"\n",
    "\n",
    "# experiment_paths = [Path(x) for x in glob.glob(\"/Volumes/Seagate Backup Plus Drive/2023-05-16/dh229*/\")]\n",
    "\n",
    "# experiment_paths = [e for e in experiment_paths if not str(e).endswith(\".zarr\")]\n",
    "\n",
    "# full_idxr = pd.DataFrame()\n",
    "# for i,e in enumerate(experiment_paths):\n",
    "#     files = [x for x in glob.glob(str(e/\"*.tif\")) if \"LUT\" not in x]\n",
    "#     indexer = pd.DataFrame([re.findall(\"(\\d+)\", Path(x).name) for x in files], columns=list('SLZC')).astype(int)\n",
    "#     indexer['Q'] = pd.Series([re.findall(\"fluo|srs\", Path(x).name)[0] for x in files],dtype='category').cat.codes\n",
    "#     indexer['filenames'] = pd.Series(files, copy=True)\n",
    "\n",
    "#     scale = indexer.groupby('S')[\"Q\"].nunique().sum()\n",
    "\n",
    "#     indexer['T']=indexer.groupby(list('SQ'))[\"L\"].transform(lambda x: (x-x.min())//scale)\n",
    "#     indexer['I'] = i\n",
    "#     indexer_check(indexer[list('SQTCZ')])\n",
    "#     full_idxr = pd.concat([full_idxr,indexer])\n",
    "\n",
    "# full_idxr['C'] = full_idxr['C']+full_idxr['Q']\n",
    "# full_idxr = full_idxr.drop(columns='Q')\n",
    "\n",
    "# readers = []\n",
    "# for i in range(2):\n",
    "#     idxr = full_idxr.loc[full_idxr['I']==i]\n",
    "#     readers.append(TiffGlobReader(idxr['filenames'], idxr[list('STCZ')]))\n",
    "\n",
    "# ds0 = readers[1].get_xarray_dask_stack(scene_character='S')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64cf5d2f-a50b-47c6-8314-3d32bd9a25c3",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec4c0172-9b78-48e0-b99b-b151f313842e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import xml.etree.ElementTree as ET\n",
    "\n",
    "# meta_path = \"/Volumes/Seagate Backup Plus Drive/2023-05-16/dh229/*_Properties.xml\"\n",
    "\n",
    "# files = sorted(glob.glob(meta_path))\n",
    "\n",
    "# indexer = full_idxr.loc[full_idxr['I']==0].drop(columns='I')\n",
    "# ldm_lookup = indexer.groupby('L')[['S','T']].min()\n",
    "\n",
    "# timestamps = []\n",
    "# dim_data = []\n",
    "# for s,file in enumerate(files):\n",
    "#     parsed = ET.parse(file)\n",
    "#     ldm = int(re.findall(\"(\\d+)_Properties\",Path(file).name)[0])\n",
    "    \n",
    "#     for x in parsed.iter(\"TimeStamp\"):\n",
    "#         d = x.attrib\n",
    "#         d[\"ldm\"] = ldm\n",
    "#         s,t = ldm_lookup.loc[ldm].values\n",
    "#         d['S'] = s\n",
    "#         d['T'] = t\n",
    "#         timestamps.append(d)\n",
    "    \n",
    "#     for x in parsed.iter(\"DimensionDescription\"):\n",
    "#         d = x.attrib\n",
    "#         d[\"ldm\"] = ldm\n",
    "#         dim_data.append(d)\n",
    "        \n",
    "# timestamps = pd.DataFrame(timestamps)\n",
    "# dim_data = pd.DataFrame(dim_data).set_index(['ldm','DimID']).unstack('DimID')\n",
    "\n",
    "# timestamps['datetime'] = timestamps.apply(\n",
    "#     lambda df: pd.to_datetime(df['Date'] + \" \" + df['Time']) + pd.to_timedelta(int(df['MiliSeconds']), unit=\"ms\"),axis=1)\n",
    "\n",
    "# timestamps = timestamps.sort_values('datetime')\n",
    "\n",
    "# timestamps.to_hdf(\"/Users/johnrussell/Data/2023-05-16-dh229_timestamps.h5\", key='df')\n",
    "\n",
    "# dim_data.to_hdf(\"/Users/johnrussell/Data/2023-05-16-dh229_dimdata`.h5\", key='df')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dae0de3b-fd57-4d9a-b554-3e30beaf32cb",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "edaaadb9-bc6e-455d-94d2-e22f6842c590",
   "metadata": {},
   "outputs": [],
   "source": [
    "timestamps = pd.read_hdf(\"/Users/johnrussell/Data/2023-05-16/2023-05-16-timestamps.h5\", key='df')\n",
    "timestamps['RelativeTime'] = timestamps['RelativeTime'].astype('f4')\n",
    "dim_data = pd.read_hdf(\"/Users/johnrussell/Data/2023-05-16/2023-05-16-dimdata.h5\", key='df')\n",
    "\n",
    "times = timestamps.groupby(list(\"ST\"))['RelativeTime'].min()/60\n",
    "dt = timestamps.groupby(['S','T']).datetime.min().groupby('S').diff()\n",
    "dt_min = dt.mean().total_seconds()/60\n",
    "\n",
    "zpos = timestamps[['S','T','ldm']].set_index('ldm').join(dim_data.loc[:,('Origin','Z')].rename(\"Z\")).groupby(['S','T']).last().astype('f4')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a31ebec8-5ac2-4f3d-8b2d-553439054742",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = xr.open_zarr(\"/Users/johnrussell/Data/2023-05-16/dh229.zarr/\").isel(T=slice(36))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa80333a-1fea-4dff-9acc-1227ab0d6e8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_data = dt_min*ds['T'].values #+28.75 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b96267b4-255e-4fb1-9a56-522584923ef0",
   "metadata": {},
   "outputs": [],
   "source": [
    "srs_squash = ds.images.isel(C=2).mean('Z', dtype='f4')#.load()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd16c904-7892-4e88-b7f6-eb48ca84354b",
   "metadata": {},
   "outputs": [],
   "source": [
    "srs_max = ds.images.isel(C=2).max('Z').astype('f4')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0cb4ff61-934c-4ce6-a4c7-9e8c9e7fe60e",
   "metadata": {},
   "outputs": [],
   "source": [
    "fluo = ds.images.isel(C=0).max('Z')#.load()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e1d9cf3-330e-4047-9941-decb05686ecc",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_srs = srs_squash.isel(T=35).load().data\n",
    "test_fluo = fluo.isel(T=35).load().data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d1e190e-521c-4b83-91f9-aa7adc977b94",
   "metadata": {},
   "outputs": [],
   "source": [
    "shift = [phase_cross_correlation(s,f)[0] for s,f in zip(test_srs, test_fluo)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e26878b9-55c2-4aae-adbe-b1c600816427",
   "metadata": {},
   "outputs": [],
   "source": [
    "fluo.load();\n",
    "ds.fmax_labels.load();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd0ca985-d559-4d7b-a942-c69f18c61211",
   "metadata": {},
   "outputs": [],
   "source": [
    "aligned_fluo = xr.concat([f.shift(dict(zip('YX', s.astype(int))), fill_value=0) for f,s in zip(fluo, shift)], dim='S')\n",
    "aligned_mask = xr.concat([f.shift(dict(zip('YX', s.astype(int))), fill_value=0) for f,s in zip(ds.fmax_labels, shift)], dim='S')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96668550-9649-4622-b0c7-b2849e2c5045",
   "metadata": {},
   "outputs": [],
   "source": [
    "if 'labels' not in ds:\n",
    "    tracked = xr.zeros_like(aligned_mask)\n",
    "    for i,s in enumerate(aligned_mask):\n",
    "        tracked.data[i] = mu.btrack.gogogo_btrack(s.data, \"cell_config.json\", 10, \"/Users/johnrussell/tmp/tracks/h5\")\n",
    "    \n",
    "    ds['labels'] = tracked\n",
    "    ds[['labels']].to_zarr(\"/Users/johnrussell/Data/2023-05-16/dh229.zarr/\", mode='a')\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2224940-e5d1-495d-bb89-b77391f92afc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# be = BackgroundEstimator(srs_squash, tracked)\n",
    "be =  BackgroundEstimator(srs_squash, ds.labels)\n",
    "be.make_cv_labels()\n",
    "# be.sigma_scan(n_samples=5)\n",
    "# be.sigma_opt.load();\n",
    "be.sigma_opt = xr.DataArray([2,66], dims=['k']) #optimum for mean projection\n",
    "be.background_estimate;\n",
    "\n",
    "# be_max =  BackgroundEstimator(srs_max, ds.labels)\n",
    "# be_max.make_cv_labels()\n",
    "# be_max.sigma_opt = xr.DataArray([2,66], dims=['k']) #optimum for mean projection\n",
    "# be_max.background_estimate;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e01d188-07a3-48e4-957f-274bc070baaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "srs_bsub = srs_squash-be.background_estimate\n",
    "# max_sub = srs_max-be_max.background_estimate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b00c662-581c-4e3b-a525-71572c266eff",
   "metadata": {},
   "outputs": [],
   "source": [
    "srs_bsub.load();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "274378a1-414c-4e79-893c-780fdf9a4f9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# bkgd_avg = be.background_estimate.mean(list('SYX'))\n",
    "# bkgd_std = be.background_estimate.std(list('SYX'))\n",
    "\n",
    "# plt.figure()\n",
    "# plt.errorbar(bkgd_avg['T'].values, bkgd_avg, yerr = bkgd_std, fmt='o', capsize=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b3259b2-1564-4442-89f0-abe8441b297c",
   "metadata": {},
   "outputs": [],
   "source": [
    "cell = ['S', 'CellID'] #shorthand for use in groupby"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c03b490-7efe-4d76-a840-788ed8f8a14e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.labels.load();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "060c9279-958f-4515-85b7-e17ab754a6e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "srs_avgs = (mu.single_cell.average(ds,srs_bsub)\n",
    "            .to_series()\n",
    "            .unstack('T')\n",
    "            .dropna(how='all')\n",
    "            .stack('T', dropna=False))\n",
    "\n",
    "# max_avgs = (mu.single_cell.average(ds,max_sub)\n",
    "#             .to_series()\n",
    "#             .unstack('T')\n",
    "#             .dropna(how='all')\n",
    "#             .stack('T', dropna=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45d57f6d-ae86-435e-bbf6-c492df90dbdf",
   "metadata": {},
   "outputs": [],
   "source": [
    "coms = (mu.single_cell.center_of_mass(ds)\n",
    "            .to_series()\n",
    "            .unstack(['com','T'])\n",
    "            .dropna(how='all')\n",
    "            .stack('T', dropna=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9498bdd8-1e51-4d76-8e68-2087a48a7c66",
   "metadata": {},
   "outputs": [],
   "source": [
    "areas = (mu.single_cell.area(ds)\n",
    "            .to_series()\n",
    "            .reorder_levels(['S','CellID','T'])\n",
    "            .loc[coms.index])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17e92b56-948c-40f4-a52b-3c3005064e79",
   "metadata": {},
   "outputs": [],
   "source": [
    "disp = np.sqrt((coms.groupby(['S','CellID']).diff()**2).sum(1, skipna=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59db47c4-2ba0-4700-95c9-2a63b412ed02",
   "metadata": {},
   "outputs": [],
   "source": [
    "first_time = coms.dropna().reset_index('T').groupby(cell).first().set_index('T', append=True).index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff106fb3-a718-48c0-b223-82da6135d7cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "disp.loc[first_time]=0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4be81b27-8d24-45ab-b203-aaef219677ec",
   "metadata": {},
   "source": [
    "**Goal**\n",
    "\n",
    "Select cells that have 15 consecutive observations with displacements less than 7\n",
    "\n",
    "- GLobal shift happens between t=0 and t=1 so ignore that in check \n",
    "  - otherwise do some time to time global registration and then leave it in\n",
    "- Some movements as large as 8 or 9 seem to be real\n",
    "- Some of the tricky tracking errors where it the cell label walks away into the daughter are quite small...\n",
    "\n",
    "*BUT* i think  the 15 consecutive observations is a more stringent cut anyway i.e. many of the tracking errors do not last that long/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9d12bd3-ba81-4e4c-8500-3ea2e08967f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_max_disps = disp.loc[pd.IndexSlice[:,:,2:]].groupby(['S','CellID']).max().dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "143b770c-4033-4d04-97f7-f3ddd59ba51f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,2,sharex=True, figsize=(9,4))\n",
    "ax[0].hist(disp.loc[pd.IndexSlice[:,:,2:]].dropna().values, bins=45, range=[0,15], )\n",
    "ax[0].set_title(\"All displacements\")\n",
    "ax[1].hist(cell_max_disps, bins=45, range=[0,15],)\n",
    "ax[1].set_title(\"Cell Maximum displacements\")\n",
    "for a in ax:\n",
    "    a.semilogy();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e855165d-5c4d-4992-a486-9b82ac8b23e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# x = disp.loc[(0, 150)]\n",
    "# df = disp.loc[(0,15)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfdb6b43-2d59-4765-a23d-7e52d0ab4945",
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_consecutive(x, N=15):\n",
    "    m = x.isna()\n",
    "    s = m.cumsum()\n",
    "    final_mask = s.map(s[~m].value_counts()).ge(N) & ~m\n",
    "    return final_mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adc96b2d-214a-4b6d-9c07-8595ca503beb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specific entries that are a part of a 20+ consecutive observation run\n",
    "nconsec = 15\n",
    "consecutive_cells = disp.groupby(['S','CellID'], group_keys=False).apply(check_consecutive, N=nconsec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75f117fc-130e-48f1-b732-2792a5c6f6d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "disp_consecutive = disp.loc[consecutive_cells].unstack('T')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c936ceb6-8ab8-4527-8181-4b3fbf6d68e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,2,sharex=True, figsize=(6.52,3.2))\n",
    "ax[0].hist(cell_max_disps.values, bins=45, range=[0,15], )\n",
    "ax[0].set_title(\"All cells\")\n",
    "ax[1].hist(disp_consecutive.max(1).values, bins=45, range=[0,15],)\n",
    "ax[1].set_title(\"Continuously observed cells\")\n",
    "ax[0].set_xlabel(\"$N_i{cells}$\")\n",
    "fig.suptitle(\"Maximum single-cell displacements\")\n",
    "for a in ax:\n",
    "    a.semilogy()\n",
    "    a.set_xlabel(\"Displacement (pixels)\")\n",
    "if savefigs:\n",
    "    plt.savefig(fig_path+\"max_disps.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18a94998-79b1-4f5b-8e20-3bee4728040f",
   "metadata": {},
   "outputs": [],
   "source": [
    "big_jump = (disp_consecutive.loc[:,2:].groupby(['S','CellID']).max()>8).stack('T')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4be76e9-8d46-4251-9df7-4cb509b112b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_labels_from_multiindex(labels, index):\n",
    "    check_labels = xr.zeros_like(labels)\n",
    "    \n",
    "    for s in range(labels.sizes['S']):\n",
    "        try:\n",
    "            checks = index.get_loc_level(s)[1].values+1\n",
    "        except KeyError:\n",
    "            continue\n",
    "            \n",
    "        mask = xr.DataArray(np.isin(labels.data[s], checks), dims=labels[s].dims)\n",
    "        check_labels[s] = labels[s].where(mask)\n",
    "    return check_labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ec4e6ff-9ee5-4ecc-9d8a-f365d3efc95f",
   "metadata": {},
   "outputs": [],
   "source": [
    "jump_idx = big_jump.loc[big_jump].reset_index('T').index.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7189f348-e727-4865-af15-09ede8660d7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check = check_labels_from_multiindex(ds.labels, jump_idx)\n",
    "\n",
    "# v = napari.Viewer()\n",
    "# il= v.add_image(srs_bsub)\n",
    "# fl = v.add_image(aligned_fluo)\n",
    "# # ll = v.add_labels(ds.labels)\n",
    "# cl = v.add_labels(check)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29ef13de-68d6-47d6-a902-df7640c5b45c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#big_jump_cell_times = big_jump.loc[big_jump].reset_index('T')['T']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac439df9-20b3-470f-8777-662670e81000",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This doesnt work any more\n",
    "\n",
    "# srs_with_jumps = srs_avgs.loc[consecutive_cells].unstack('T').loc[jump_idx]\n",
    "# pts = srs_with_jumps.loc[jump_idx].max(1)\n",
    "\n",
    "# fig, ax = plt.subplots()\n",
    "# ax.plot(srs_with_jumps.values.T, 'k', linewidth=1, alpha=0.2);\n",
    "# ax.plot(jump_idx.get_level_values('T'), pts.values, 'rx')\n",
    "# plt.ylim(-3, 38);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27759c13-09eb-49d2-a06c-376902113401",
   "metadata": {},
   "outputs": [],
   "source": [
    "no_jumps = big_jump.groupby(cell).sum()==0\n",
    "# no_jump_idx = no_jumps.loc[no_jumps].index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "081bef57-3988-46f1-99cf-33bf2ee61494",
   "metadata": {},
   "outputs": [],
   "source": [
    "consecutive_coms = coms.loc[consecutive_cells]\n",
    "central_cells = ((consecutive_coms['X'].between(5,507))&(consecutive_coms['Y'].between(5,507))).groupby(cell).all()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44a9f3a0-58a1-403f-b2e8-95c4c50e1918",
   "metadata": {},
   "outputs": [],
   "source": [
    "srs_clean = srs_avgs.loc[consecutive_cells].unstack('T').loc[no_jumps&central_cells]\n",
    "# max_clean = max_avgs.loc[consecutive_cells].unstack('T').loc[no_jumps&central_cells]\n",
    "# raw_clean = raw_avgs.loc[consecutive_cells].unstack('T').loc[no_jumps&central_cells]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7eb16597-1440-46c6-bab0-0f2d982eb358",
   "metadata": {},
   "outputs": [],
   "source": [
    "consecutive_cells.unstack('T').sum().max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29816078-e358-4b78-b1db-c57a5b677fda",
   "metadata": {},
   "outputs": [],
   "source": [
    "central_cells.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f484ae58-a319-486d-bf9f-824dcda44a7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "srs_clean.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5b016d9-e418-4bb8-92d8-07fcedff235b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# srs_clean.to_csv(\"/Users/johnrussell/Data/2023-05-16/2023-05-16-srs.csv\")\n",
    "# areas_clean.to_csv(\"/Users/johnrussell/Data/2023-05-16/2023-05-16-area.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb4f5364-305c-4a07-b744-1f2c9b074226",
   "metadata": {},
   "outputs": [],
   "source": [
    "areas_clean = areas.loc[consecutive_cells].unstack('T').loc[no_jumps&central_cells]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fc4b287-cd3d-45a8-b170-3339700aedca",
   "metadata": {},
   "outputs": [],
   "source": [
    "# first_times_clean = areas_clean.stack('T').reset_index('T').groupby(cell).first().set_index('T', append=True).index\n",
    "first = srs_clean.stack('T').reset_index('T').groupby(cell).first().set_index('T', append=True).iloc[:,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "917d6888-f559-43fe-a82e-0303ceb9a842",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_first = srs_avgs.reset_index('T').dropna().groupby(cell).first().set_index('T', append=True).iloc[:,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31cd6d16-4d61-406f-bdb4-ecc4e0607566",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(t_data, areas_clean.values.T, 'k', alpha=0.025);\n",
    "# plt.errorbar(t_data, areas_clean.median().values, yerr=(areas_clean.quantile([0.25,0.75])-areas_clean.median()).abs().values, fmt='mo', capsize=3)\n",
    "plt.ylim([0,400])\n",
    "plt.title(\"Single Cell Area Trajectories\")\n",
    "plt.ylabel(\"Area (pixels)\")\n",
    "plt.xlabel(\"Time (Minutes)\")\n",
    "if savefigs:\n",
    "    plt.savefig(fig_path+\"cell_areas.png\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "691e799a-aec1-4f2c-924c-ffe89450fda4",
   "metadata": {},
   "outputs": [],
   "source": [
    "growth_cells = areas_clean.loc[(areas_clean.stack('T').loc[first.index]<100).droplevel(\"T\")]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd990bb0-b87f-4a34-926a-0382b117340b",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(t_data, growth_cells.values.T, 'k', alpha=0.025);\n",
    "# plt.errorbar(t_data, areas_clean.median().values, yerr=(areas_clean.quantile([0.25,0.75])-areas_clean.median()).abs().values, fmt='mo', capsize=3)\n",
    "plt.ylim([0,400])\n",
    "plt.title(\"Single Cell Area Trajectories\")\n",
    "plt.ylabel(\"Area (pixels)\")\n",
    "plt.xlabel(\"Time (Minutes)\")\n",
    "if savefigs:\n",
    "    plt.savefig(fig_path+\"cell_areas.png\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "982667cf-d140-4524-8f3a-d86fbd1cc4d1",
   "metadata": {},
   "source": [
    "**Note** `fit_all_exponentials` should maybe take an arg about an offset and whether to use absolute or relative t values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2959049b-0928-4b41-a9fa-99b2217b8dfb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# srs_jumps = (srs_clean.diff(axis=1).abs()>15).any(1)\n",
    "# check_idx = srs_jumps.loc[srs_jumps].index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94d7def6-7b28-4e7c-91ac-81c549a3e260",
   "metadata": {},
   "outputs": [],
   "source": [
    "# time_count = srs_avgs.groupby(['S','CellID']).count()\n",
    "\n",
    "# present = (~srs_avgs.unstack(['S','CellID']).isna()).astype(int)\n",
    "# appears_once = (present.diff()==1).astype(int).sum()<2\n",
    "# disappears_once = (present.diff()==-1).astype(int).sum()<2\n",
    "\n",
    "# keep = srs_avgs.unstack('T').loc[(appears_once & disappears_once & (time_count>15)), :35]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76dc96b3-f37f-4e4a-97e6-555a97ed3c57",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,5, sharex=True, sharey=True, figsize=(6.52,3))\n",
    "for s in range(5):\n",
    "    traces = ax[s].plot(srs_clean.loc[s].values.T, 'k', linewidth=1, alpha=0.05)\n",
    "    ax[s].plot(5*(zpos.loc[s]-zpos.loc[s].mean())-10, '.-', label='Scaled Z')\n",
    "    \n",
    "ylim = ax[0].get_ylim()\n",
    "ax[0].plot(20,100,'k', label='Cell Traces')\n",
    "ax[0].set_ylim(ylim)\n",
    "ax[0].legend();\n",
    "plt.suptitle(\"\\\"Sawtoothing\\\" corresponds with changes in Z\")\n",
    "if savefigs:\n",
    "    plt.savefig(fig_path+\"sawtooth_abs_z.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bef7452c-ae2e-444c-a537-585249b53b0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig, ax = plt.subplots(1,5, sharex=True, sharey=True, figsize=(10,3))\n",
    "# for s in range(5):\n",
    "#     traces = ax[s].plot(max_clean.loc[s].values.T, 'k', linewidth=1, alpha=0.05)\n",
    "#     ax[s].plot(5*(zpos.loc[s]-zpos.loc[s].mean())-10, '.-', label='Scaled Z postion')\n",
    "    \n",
    "# ylim = ax[0].get_ylim()\n",
    "# ax[0].plot(20,100,'k', label='Cell Traces')\n",
    "# ax[0].set_ylim(ylim)\n",
    "# ax[0].legend()\n",
    "# plt.suptitle(\"Sawtooth seems to correspond to changes in Z (max projections)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7e76768-a45f-488f-9afe-a571673b73a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "savefigs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e143bf1-8f77-4657-97fc-3a3a7a607d2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "dz = zpos.groupby('S').diff()['Z'].loc[(slice(None), slice(srs_clean.columns.max()))]\n",
    "dsrs = srs_clean.diff(axis=1)#/dt_min\n",
    "\n",
    "fig, ax = plt.subplots(2,1, sharex=True, )\n",
    "for s in range(5):\n",
    "    ax[0].plot(t_data, dsrs.loc[s].mean().values, color=colors[s])\n",
    "    ax[1].plot(t_data, dz.loc[s].values, color=colors[s])    \n",
    "ax[0].plot(t_data, dsrs.mean().values, 'k--', label=\"Mean\")\n",
    "ax[0].legend(loc='lower left')\n",
    "# ax[0].plot(xplt, dy, 'k')\n",
    "ax[1].set_xlabel(\"Time (minutes)\")\n",
    "ax[0].set_ylabel(\"$\\Delta$ SRS (a.u.)\")\n",
    "ax[1].set_ylabel(\"$\\Delta$ Z (µm)\")\n",
    "ax[0].set_title(\"Z-oscillations correspond with \\\"Sawtoothing\\\"\")\n",
    "if savefigs: plt.savefig(fig_path+\"sawtooth_deltas.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "157016bc-056a-427d-b339-13276e72ee5f",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Additive correction - obsolete"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3e43c4e-d57d-499b-93bf-ed625530e99e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# diff = srs_clean.diff(axis=1)\n",
    "# test_diff = diff.loc[0]\n",
    "# test_dz = dz.loc[0]\n",
    "\n",
    "# dz_gb = diff.stack('T').unstack('CellID').groupby(dz.fillna(0)).mean()\n",
    "\n",
    "# plt.figure()\n",
    "# plt.errorbar(dz_gb.index.values, dz_gb.mean(1), yerr=dz_gb.std(1), fmt='o', capsize=3)\n",
    "# plt.xlabel(\"$\\Delta Z$\")\n",
    "# plt.ylabel(\"$\\Delta SRS$\")\n",
    "\n",
    "# fig, ax = plt.subplots(3,1)\n",
    "# ax[0].plot(test_diff.std().values)\n",
    "# # ax[0].plot((mcorr*test_diff).std().values, label='mult. corr.')\n",
    "# ax[0].set_ylabel(\"$\\sigma_{\\Delta SRS}$\")\n",
    "\n",
    "# ax[1].plot(test_dz)\n",
    "# ax[1].set_ylabel(\"$\\Delta$ Z (µm)\")\n",
    "\n",
    "# ax[2].plot(test_diff.mean().values)\n",
    "# ax[2].set_ylabel(\"$\\Delta SRS$\")\n",
    "\n",
    "# fig.suptitle(\"Z-changes do not change the variance\")\n",
    "# if savefigs:\n",
    "#     plt.savefig(fig_path+\"dz_var.png\"))\n",
    "\n",
    "# a,k=fit[0]\n",
    "# dy = a*k*np.exp(-k*xplt)\n",
    "# y_pred = exp_approach(t_data, a,k)\n",
    "# dy_pred = a*k*np.exp(-k*t_data)*dt_min\n",
    "\n",
    "# res_dsrs = (srs_clean.diff(axis=1)- dy_pred).groupby('S').mean()\n",
    "\n",
    "# fig, ax = plt.subplots(1,2, sharey=True)\n",
    "# ax[0].plot(t_data, dsrs.groupby('S').mean().values.T)\n",
    "# ax[0].plot(t_data,np.zeros_like(t_data),'k')\n",
    "# ax[1].plot(t_data, res_dsrs.values.T)\n",
    "# ax[1].plot(t_data,np.zeros_like(t_data),'k')\n",
    "\n",
    "# fig, ax = plt.subplots(1,2, sharex=True, sharey=True)\n",
    "# ax[0].plot(test_diff.values.T, 'k', linewidth=1, alpha=0.1)\n",
    "# # ax[0].plot(5*(zpos.loc[s]-zpos.loc[s].mean())-10, '.-',\n",
    "# # ax[1].plot(test_corrected_diff.values.T, 'k', linewidth=1, alpha=0.1);\n",
    "# ax[1].plot((test_diff - test_diff.mean()).values.T, 'k', linewidth=1, alpha=0.1);\n",
    "# # ax[0].axhline(5, color='r')\n",
    "# # ax[0].axhline(-5, color='b')\n",
    "# # ax[1].axhline(5, color='r')\n",
    "# # ax[1].axhline(-5, color='b')\n",
    "# ax[0].set_title(\"Raw\")\n",
    "# ax[1].set_title(\"Corrected\")\n",
    "# ax[0].set_ylabel(\"$\\Delta SRS$\")\n",
    "# ax[0].set_xlabel(\"Time index\")\n",
    "# ax[1].set_xlabel(\"Time index\")\n",
    "# ax[0].plot(5*test_dz.values-10,'-o')\n",
    "# ax[1].plot(5*test_dz.values-10,'-o')\n",
    "# if savefigs:\n",
    "#     plt.savefig(fig_path+\"test_fov_corr.png\")\n",
    "\n",
    "# fig, ax = plt.subplots(5,2, sharex=True, sharey=True)\n",
    "# for s in range(5):\n",
    "#     ax[s,0].plot(diff.loc[s].values.T, 'k', linewidth=1, alpha=0.02)\n",
    "#     ax[s,1].plot((diff.loc[s]-diff.loc[s].mean()).values.T, 'k', linewidth=1, alpha=0.02);\n",
    "#     ax[s,0].axhline(5, color='r')\n",
    "#     ax[s,0].axhline(-5, color='b')\n",
    "#     ax[s,1].axhline(5, color='r')\n",
    "#     ax[s,1].axhline(-5, color='b')\n",
    "#     ax[s,0].set_ylabel(\"$\\Delta SRS$\")\n",
    "#     # ax[s,2].plot(5*dz.loc[s].fillna(0).values,'-o')\n",
    "#     # ax[s,1].plot(5*dz.loc[s].values-10,'-o')\n",
    "# ax[0,0].set_ylim([-10,10])\n",
    "# ax[0,0].set_title(\"Raw\")\n",
    "# ax[0,1].set_title(\"Corrected\")\n",
    "# # ax[0,2].set_title(\"$\\Delta Z$\")\n",
    "# ax[s,0].set_xlabel(\"Time index\")\n",
    "# ax[s,1].set_xlabel(\"Time index\")\n",
    "\n",
    "# diff = (srs_clean.diff(axis=1) - dy_pred)\n",
    "\n",
    "# offset = (diff.groupby('S').transform(lambda df: df-df.mean())+dy_pred).stack('T', dropna=False)\n",
    "# # offset.loc[first_times_clean] = srs_clean.stack('T').loc[first_times_clean]\n",
    "# offset.loc[first.index] = first\n",
    "\n",
    "# srs_recon = offset.unstack('T').cumsum(axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f288f95b-b09a-4654-ac5d-aa8e3f59edd5",
   "metadata": {},
   "source": [
    "### Multiplicative correction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e409dae3-692f-4771-a95c-3d787e2b6f91",
   "metadata": {},
   "outputs": [],
   "source": [
    "k = np.random.randn((ds.sizes['T'])).astype('f4')\n",
    "k0 = np.zeros(ds.sizes['T']).astype('f4')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fcec7b1-733f-4b35-ae27-acf02e45bdbc",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "851a0f76-058f-447e-aac2-6a4f9df99601",
   "metadata": {},
   "outputs": [],
   "source": [
    "def d2loss_with_reg(k, x, alpha=1):\n",
    "    d = np.diff((1+k)*x, n=2, axis=-1)\n",
    "    d = d*d\n",
    "    w = np.convolve((1+k), [1,2,1], mode='valid')\n",
    "    w = w*w\n",
    "    return np.nanmean(d/w) + alpha*np.mean(k*k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad9056b4-2f21-458b-8b92-6a1d1419a303",
   "metadata": {},
   "outputs": [],
   "source": [
    "def d2loss_with_reg_bcs(k, x, alpha=1.):\n",
    "    s = 1+k\n",
    "    y = s*x\n",
    "    d = np.zeros_like(y)\n",
    "    d[:,1:-1] = -2*y[:,1:-1] + y[:,2:] + y[:,:-2]\n",
    "    d[:,0] = y[:,0] -2*y[:,1] + y[:,2]\n",
    "    d[:,-1] = y[:,-1] -2*y[:,-2] + y[:,-2]\n",
    "    d = d*d\n",
    "    w = np.zeros_like(s)\n",
    "    w[1:-1] = 2*s[1:-1] + s[2:] +s[:-2]\n",
    "    w[0] = s[0]+2*s[1]+s[2]\n",
    "    w[-1] = s[-1]+2*s[-2]+s[-3]\n",
    "    w = w*w\n",
    "    return np.nanmean(d/w) + alpha*np.mean(k*k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28ac4c81-8576-42b8-ab33-d68718fe1e9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# def d1loss_with_reg(k, alpha=1):\n",
    "#     d = np.diff((1+k)*test, n=1, axis=1)\n",
    "#     return np.nanmean(d*d) + alpha*np.mean(k*k)\n",
    "\n",
    "# def central_d1loss_with_reg(k, alpha=1):\n",
    "#     y = np.pad((1+k)*test, 1, mode='edge')\n",
    "#     d = np.roll(y,1, axis=-1)-np.roll(y,-1, axis=-1)\n",
    "#     d = d[:,1:-1]\n",
    "#     return np.nanmean(d*d) + alpha*np.mean(k*k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8e3f01d-eef0-451e-9501-9a704796adc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "na = 9\n",
    "nt = srs_clean.shape[1]\n",
    "ns = 5\n",
    "alphas = np.logspace(-2,2, num=na)\n",
    "avgs = np.zeros((ns, nt, na))\n",
    "kopts = np.zeros((ns, nt, na))\n",
    "tloss = np.zeros((ns, na))\n",
    "\n",
    "for s in tqdm(range(ns)):\n",
    "    y = srs_clean.loc[s].values\n",
    "    val_mask = np.zeros(y.shape[0],dtype='bool')\n",
    "    val_mask[np.random.choice(np.arange(y.shape[0]),size=50, replace=False)]=1\n",
    "    val = y[val_mask]\n",
    "    train = y[~val_mask]\n",
    "\n",
    "    for i,a in enumerate(tqdm(alphas, leave=False)):\n",
    "        res = minimize(d2loss_with_reg_bcs, k0, (train, a))\n",
    "        kopts[s,:,i] = res.x\n",
    "        tloss[s,i] = res.fun\n",
    "        avgs[s,:,i] = np.nanmean((1+res.x)*train, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efb04439-946e-44c4-b411-8931122b6e5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "viridis_colors = plt.cm.viridis(np.linspace(0,1,alphas.shape[0]))\n",
    "fig, axs = plt.subplots(1,5, sharex=True, sharey=True)\n",
    "for s, ax in zip(range(ns), axs):\n",
    "    for x, c, a in zip(avgs[s].T, viridis_colors, alphas):\n",
    "        ax.plot(x, color=c, label=f\"{a:0.2e}\")\n",
    "    ax.plot(srs_clean.loc[s].mean().values, 'k', label='Original')\n",
    "fig.suptitle('Second Derivative Correction')\n",
    "# axs[-1].legend(title=\"Regularization parameter\", loc='upper right',bbox_to_anchor=(1.35, 1.05))\n",
    "if savefigs: plt.savefig(fig_path+\"wt_d2_correction.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b3cb164-1b10-4926-8a26-83009ddc6e0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "kopt= pd.DataFrame(kopts[...,4], index=pd.Index(range(kopts.shape[0]), name='S'), columns=pd.Index(range(kopts.shape[1]), name='T'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00004b09-e2e9-4c05-aa84-696d11ae289f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,2,sharex=True, sharey=True,figsize=(6.52,3))\n",
    "for i in range(5):\n",
    "    ax[0].plot(t_data, srs_clean.loc[i].mean().values, color=colors[i])\n",
    "    ax[1].plot(t_data, (1+kopt.loc[i].values)*srs_clean.loc[i].mean().values, color=colors[i])\n",
    "ax[0].set_title(\"Raw SRS Intensity\")\n",
    "ax[1].set_title(\"SRS Intensity with smoothing\")\n",
    "ax[0].set_xlabel(\"Time (minutes)\")\n",
    "ax[1].set_xlabel(\"Time (minutes)\")\n",
    "ax[0].set_ylabel(\"SRS Intenity (a.u.)\")\n",
    "if savefigs: plt.savefig(fig_path+\"sawtooth_corr.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b9190ab-696e-4261-8af8-ab29fca088c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(3,1, sharex=True, )\n",
    "for s in range(5):\n",
    "    ax[0].plot(t_data, dsrs.loc[s].mean().values, color=colors[s])\n",
    "    ax[1].plot(t_data, dz.loc[s].values, color=colors[s])    \n",
    "    ax[2].plot(t_data, 1+kopt.loc[s], color=colors[s])\n",
    "ax[0].plot(t_data, dsrs.mean().values, 'k--', label=\"Mean\")\n",
    "ax[0].legend(loc='lower left')\n",
    "# ax[0].plot(xplt, dy, 'k')\n",
    "ax[1].set_xlabel(\"Time (minutes)\")\n",
    "ax[0].set_ylabel(\"$\\Delta SRS$\")\n",
    "ax[1].set_ylabel(\"$\\Delta$ Z (µm)\")\n",
    "ax[0].set_title(\"Z-oscillations correspond with \\\"Sawtoothing\\\"\")\n",
    "if savefigs:\n",
    "    plt.savefig(fig_path+\"wt_sawtooth_deltas_weights.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a867a886-5014-4858-af4d-3a5e1a3c7eee",
   "metadata": {},
   "outputs": [],
   "source": [
    "srs_recon = (1+kopt)*srs_clean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f90a595-3dc2-4bae-b728-83357dfef71d",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_recon = (1+kopt)*(srs_avgs.unstack('T'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c166e8c-5be1-4072-bfda-083afc751483",
   "metadata": {},
   "outputs": [],
   "source": [
    "working_dir = \"/Users/johnrussell/Data/2023-05-16/\"\n",
    "ename = 'dh229'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7312ace0-45c7-411a-b75b-d4d660003970",
   "metadata": {},
   "source": [
    "### Inspect rescaled data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03b5b2d3-76a7-4c84-a5c3-4e2e0a333854",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,2,sharex=True, sharey=True, figsize=(6.52,3.2))\n",
    "ax[0].plot(t_data, srs_clean.values.T, 'k', alpha=0.01);\n",
    "ax[0].set_title('Raw')\n",
    "\n",
    "ax[1].plot(t_data, srs_recon.values.T, 'k', alpha=0.01);\n",
    "ax[1].set_title('Corrected')\n",
    "ax[0].set_ylabel(\"SRS (a.u.)\")\n",
    "fig.suptitle(\"Second Derivative Smoothing Reduces Sawtoothing\")\n",
    "for a in ax:\n",
    "    a.set_xlabel(\"Time (minutes)\");\n",
    "    a.set_ylim([-10,40])\n",
    "if savefigs:\n",
    "    plt.savefig(fig_path+\"saw_corr_before_after.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08b97b64-7dc2-4385-81ed-d454173678c7",
   "metadata": {},
   "source": [
    "## Curve fitting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "923da462-3667-409e-b1b7-25d174aa8ef3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86e499d0-563b-4bd7-943b-b91b70b1a68d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def exp_approach(x, a, k):\n",
    "    return a*(1-np.exp(-k*x))#+b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1a2ecac-d277-436f-8483-080e2fbd7724",
   "metadata": {},
   "outputs": [],
   "source": [
    "def exp_approach2(x, a, k,t):\n",
    "    return a*(1-np.exp(-k*(x-t)))#+b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6801350f-5178-4f2e-8d46-71f6eb72ea4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def exp_approach3(x, a, k, b):\n",
    "    return a*(1-np.exp(-k*x))+b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2a07588-7d7d-4402-b193-4e704d55cbc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def exp_decay(x, a, b, k):\n",
    "    return a*np.exp(-k*x)#+b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2583c82-442d-4ca9-91b4-d64d19901367",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df = srs_recon\n",
    "# xdata = dt_min*df.columns.values + 30\n",
    "y = df.mean().values\n",
    "y_offset = y[0]\n",
    "y = y-y_offset\n",
    "sigma = df.sem().values\n",
    "\n",
    "fit = curve_fit(exp_approach, t_data, y, sigma=sigma, p0=[20, 0.1], bounds=(0,np.inf))\n",
    "\n",
    "a,k = fit[0]\n",
    "eqn = f\"Exponential fit: \\n{a:0.1f} (1-exp[-{k:0.3g}t])\"\n",
    "xplt = np.linspace(0, 550, 250)\n",
    "plt.figure()\n",
    "plt.plot(xplt, exp_approach(xplt, *fit[0])+y_offset, zorder=0, linewidth=3, alpha=0.75, label=eqn)\n",
    "# plt.plot(t_data, y+y_offset, 'k.', zorder=-1)\n",
    "plt.errorbar(t_data, y+y_offset,yerr=sigma, fmt='k.', capsize=3, zorder=-1)\n",
    "plt.title(f\"Wild Type Net Glucose Uptake\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"Time (minutes)\")\n",
    "plt.ylabel(\"SRS Intensity (a.u.)\")\n",
    "if savefigs:\n",
    "    plt.savefig(fig_path+\"bulk_fit.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25463c9c-d360-416e-81d6-8e2365da1d91",
   "metadata": {},
   "outputs": [],
   "source": [
    "1-(1/np.e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f64f0a4-b3a7-4a07-b042-c2f4a4ce4682",
   "metadata": {},
   "outputs": [],
   "source": [
    "k*60"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72dc05d3-5e28-40e5-9c2e-04b1ca3006f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sqrt(fit[1][1,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c825fd6a-ca67-4cfd-bde8-94fae7a7d380",
   "metadata": {},
   "outputs": [],
   "source": [
    "# srs_total = srs_recon*areas_clean\n",
    "\n",
    "# fig,ax = plt.subplots(1,2, sharex=True,figsize=(6.52, 3.2))\n",
    "# ax[0].plot(t_data, srs_recon.values.T, 'k', linewidth=1, alpha=0.02)\n",
    "# ax[0].set_xlabel(\"Time (minutes)\")\n",
    "# ax[0].set_ylabel(\"SRS Intensity (a.u.)\")\n",
    "# ax[0].set_title(\"Averages\")\n",
    "# ax[0].set_ylim([-10,40])\n",
    "# ax[1].plot(t_data, srs_total.values.T, 'k', linewidth=1, alpha=0.02)\n",
    "# ax[1].set_xlabel(\"Time (minutes)\")\n",
    "# # ax[1].set_ylabel(\"SRS Intensity (a.u.)\")\n",
    "# ax[1].set_title(\"Totals\")\n",
    "# ax[1].set_ylim([-2000, 8000])\n",
    "# fig.suptitle(\"Single Cell Traces\");\n",
    "# if savefigs:\n",
    "#     plt.savefig(fig_path+\"sc_area_total.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e5cef0c-cd3b-4db5-809e-3bf7afa60664",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Single Cell Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d8f5e31-7766-4b5c-80ae-cb62de1ecdd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit_exp(x, y, p0):\n",
    "    try:\n",
    "        p, cov = curve_fit(exp_approach, x, y, p0=p0, bounds=(0,np.inf), absolute_sigma=True)\n",
    "        out = np.array([*p, *np.diag(cov)])\n",
    "    except RuntimeError:\n",
    "        out = np.full(2*len(p0), np.nan)\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbc50cde-6f5b-41a2-8f3b-2b9ec2f6d76f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# slower to use dask delayed by ~50%\n",
    "df = srs_recon\n",
    "\n",
    "all_srs_params = pd.DataFrame(index=df.index, columns=['a','k', 'sig_a', 'sig_k','x0','y0'], dtype='f4')\n",
    "\n",
    "for i, (idx, s) in enumerate(tqdm(df.iterrows(), total=len(df))):\n",
    "    y = s.values\n",
    "    mask = ~np.isnan(y)\n",
    "    y = y[mask]\n",
    "    y_offset = y[0]\n",
    "    y = y-y_offset\n",
    "    x = t_data[mask]\n",
    "    x_offset = x[0]\n",
    "    x = x - x_offset\n",
    "    out = fit_exp(x,y,fit[0])\n",
    "    all_srs_params.loc[idx] = [*out, x_offset, y_offset]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d5f4bb9-24f7-4e89-b2ba-7ee02a765145",
   "metadata": {},
   "outputs": [],
   "source": [
    "unconverged = all_srs_params.loc[all_srs_params.isna().any(axis=1)].index\n",
    "srs_params = all_srs_params.dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7218818f-e0b9-49ab-9254-6c91600fd7bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# outfile = working_dir+f\"{ename}_srs_tables.h5\"\n",
    "# srs_clean.to_hdf(outfile, key='raw')\n",
    "# srs_recon.to_hdf(outfile, key='recon')\n",
    "# areas_clean.to_hdf(outfile, key='area')\n",
    "# srs_params.to_hdf(outfile, key='params')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e506fbb1-2120-4471-a4b7-aed4743f1253",
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig, ax = plt.subplots(1,2,sharex=True, sharey=True, figsize=(9,4))\n",
    "# ax[0].plot(t_data, srs_recon.loc[srs_params.index].values.T, 'k', alpha=0.02);\n",
    "# ax[0].set_title('converged')\n",
    "# ax[0].set_ylim([-10,40])\n",
    "\n",
    "# ax[1].plot(t_data, srs_recon.loc[unconverged].values.T, 'k', alpha=0.025);\n",
    "# ax[1].set_title('unconverged')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5087602a-dadd-4190-8ae8-6c6b318a6fb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "pcts = np.linspace(0,1, 101)\n",
    "qs = srs_params.quantile(pcts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e5b3340-f92e-4449-bc10-c7a0a2cb5837",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "tc = 3\n",
    "plt.plot(pcts, qs.sig_a.values, '-o', markersize=2, label='sig_a')\n",
    "plt.plot(pcts, qs.sig_k.values,'-o', markersize=2,label='sig_k')\n",
    "plt.legend()\n",
    "plt.semilogy();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7e6cfd3-95d8-4162-8e61-ec4069465cef",
   "metadata": {},
   "outputs": [],
   "source": [
    "vqs = (srs_params['a']*srs_params['k']).quantile(pcts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee8e7f13-8e69-4223-99ee-82a20d3c1f67",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "\n",
    "plt.plot(pcts, qs.a.values, '-o', markersize=2, label='a')\n",
    "plt.plot(pcts, qs.k.values,'-o', markersize=2,label='k')\n",
    "plt.plot(pcts, qs.sig_a.values, '-o', markersize=2, label='sig_a')\n",
    "plt.plot(pcts, qs.sig_k.values,'-o', markersize=2,label='sig_k')\n",
    "plt.plot(pcts, vqs.values, '-o', markersize=2, label='v')\n",
    "plt.legend()\n",
    "plt.semilogy();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cbd8330-4023-4138-9f60-7b864e20809d",
   "metadata": {},
   "outputs": [],
   "source": [
    "low_q = (srs_params<qs.iloc[2]).any(axis=1)\n",
    "low_q = low_q.loc[low_q]\n",
    "print(len(low_q))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3f05207-98de-43f0-8f9e-8d42b06a6819",
   "metadata": {},
   "outputs": [],
   "source": [
    "high_q = (srs_params>qs.iloc[-2]).any(axis=1)\n",
    "high_q = high_q.loc[high_q]\n",
    "print(len(high_q))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "812e3b97-304b-4c37-8ffd-edff96c62ab6",
   "metadata": {},
   "outputs": [],
   "source": [
    "mid_q = ((srs_params>qs.iloc[2]).all(axis=1))&((srs_params<qs.iloc[-2]).all(axis=1))\n",
    "mid_q = mid_q.loc[mid_q]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6d29d5e-49ac-4b70-88d3-54a762cf2afc",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,3, sharey=True)\n",
    "ax[0].plot(srs_recon.loc[low_q.index].values.T, 'k', alpha=0.05)\n",
    "ax[1].plot(srs_recon.loc[mid_q.index].values.T, 'k', alpha=0.05)\n",
    "ax[2].plot(srs_recon.loc[high_q.index].values.T, 'k', alpha=0.05);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12898a07-ddc0-4c76-b23b-84f96cf33a12",
   "metadata": {},
   "outputs": [],
   "source": [
    "k_sig_k =np.sqrt(srs_params['sig_k'])/ srs_params['k']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11c05300-d136-484d-8847-f7d084e35a6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "_ = plt.hist(k_sig_k.loc[np.isfinite(k_sig_k)].values, bins=100, range=(0,5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9878e40a-f101-48f8-9ba3-666220025a54",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = srs_recon.loc[srs_params.index]\n",
    "fig, ax = plt.subplots(1,3,sharex=True, sharey=True, figsize=(10,3))\n",
    "ax[0].plot(t_data, df.loc[srs_params.index].loc[srs_params['sig_k']==0].values.T, 'k', alpha=0.2);\n",
    "ax[0].set_title('$\\sigma_k=0$')\n",
    "\n",
    "ax[1].plot(t_data, df.loc[srs_params['sig_k']>qs.loc[0.99, 'sig_k']].values.T, 'k', alpha=0.2);\n",
    "ax[1].set_title('$\\sigma_k -> inf$')\n",
    "\n",
    "ax[2].plot(t_data, df.loc[srs_params['sig_k'].between(*qs.loc[[0.02, 0.99],'sig_k'])].values.T, 'k', alpha=0.02);\n",
    "ax[2].set_title('$\\sigma_k \\in [1\\%, 99\\%]$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be4e219e-92da-4776-bb96-0f9baa15c07f",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = srs_recon.loc[srs_params.index]\n",
    "fig, ax = plt.subplots()\n",
    "# l1 = plt.plot(t_data, df.loc[(srs_params['sig_k']==0)&(srs_params['sig_a']==0)].values.T, color=colors[0]);\n",
    "l2=  plt.plot(t_data, df.loc[(srs_params['sig_k']==0)&(srs_params['sig_a']>0)].values.T, color=colors[1]);\n",
    "l3=  plt.plot(t_data, df.loc[(srs_params['sig_k']>1)].values.T, color=colors[2]);\n",
    "plt.legend([l[0] for l in ( l2, l3)],['$\\sigma_k=0$ $\\sigma_a=0$', '$\\sigma_k=0$ $\\sigma_a>0$', '$\\sigma_k>1$']);\n",
    "plt.title(\"Traces with uncoverged exponential paramters\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05e2b169-e7dd-472f-bac7-f4274d81de25",
   "metadata": {},
   "outputs": [],
   "source": [
    "sig_k_sel = srs_params['sig_k'].between(*qs.loc[[0.02, 0.98],'sig_k'])#.between(1e-16, 1e4)\n",
    "sig_a_sel = srs_params['sig_a'].between(*qs.loc[[0.02, 0.98],'sig_a'])#.between(1e-16,1e5)#.between(1e-16,10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e392b731-e100-402d-8629-44768a394010",
   "metadata": {},
   "outputs": [],
   "source": [
    "hrange =(0, srs_params['k'].quantile([0.98]).item())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e14494d-51f1-43da-bdd8-64c5f41ae28d",
   "metadata": {},
   "outputs": [],
   "source": [
    "w = 1./srs_params.loc[sig_k_sel,'sig_k']\n",
    "# w = w*w\n",
    "w /= w.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94948111-0e81-4e51-b186-ce6b14bb9476",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,2, figsize=(6.52, 2.5))\n",
    "df = srs_params.loc[sig_k_sel&sig_a_sel]\n",
    "hranges = ([0,0.05], [0,50])\n",
    "for i,x in enumerate('ka'):\n",
    "    # w = 1./df[f'sig_{x}']\n",
    "    # w /= w.sum()\n",
    "    \n",
    "    counts, bins, patches = ax[i].hist(df[x].values, bins=100, range=hranges[i])\n",
    "    print(x, bins[np.argmax(counts)], df[x].mean())\n",
    "    ax[i].axvline(df[x].median(), color='k', alpha=0.5, label=f'Median = {df[x].median():0.3g}')\n",
    "    ax[i].legend()\n",
    "ax[0].set_xlabel(r\"Uptake Rate  $\\left(min^{-1}\\right)$\")\n",
    "ax[1].set_xlabel(r\"Amplitude $\\left( a.u. \\right)$\")\n",
    "\n",
    "ax[0].set_ylabel(\"Cell counts\")\n",
    "\n",
    "fig.suptitle(\"Estimated Single-Cell Metabolic Parameters\")\n",
    "if savefigs:\n",
    "    plt.savefig(fig_path+\"sc_params.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f4913cf-ec4f-4884-826f-be1f0ed1a431",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(srs_recon.loc[srs_recon.index.intersection(srs_params.loc[(np.sqrt(srs_params.sig_k)/srs_params.k)>=1].index)].values.T, 'k', alpha=0.1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42a67876-70ef-4005-940f-272ed15b22e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,2, figsize=(6.52, 2.5))\n",
    "df = srs_params.loc[(np.sqrt(srs_params.sig_k)/srs_params.k)<=1]\n",
    "hranges = ([0,0.05], [0,50])\n",
    "for i,x in enumerate('ka'):\n",
    "    # w = 1./df[f'sig_{x}']\n",
    "    # w /= w.sum()\n",
    "    \n",
    "    counts, bins, patches = ax[i].hist(df[x].values, bins=100, range=hranges[i])\n",
    "    print(x, bins[np.argmax(counts)], df[x].mean())\n",
    "    ax[i].axvline(df[x].median(), color='k', alpha=0.5, label=f'Median = {df[x].median():0.3g}')\n",
    "    ax[i].legend()\n",
    "ax[0].set_xlabel(r\"Uptake Rate  $\\left(min^{-1}\\right)$\")\n",
    "ax[1].set_xlabel(r\"Amplitude $\\left( a.u. \\right)$\")\n",
    "\n",
    "ax[0].set_ylabel(\"Cell counts\")\n",
    "\n",
    "fig.suptitle(\"Estimated Single-Cell Metabolic Parameters\")\n",
    "if savefigs:\n",
    "    plt.savefig(fig_path+\"sc_params.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb7a3dee-0b45-4c6e-9cb1-be5f50babfcb",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.drop((0,96))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27450489-a89e-49d0-8b01-c63770a9c798",
   "metadata": {},
   "outputs": [],
   "source": [
    "bigk = df.loc[df.k>0.025]\n",
    "smallk = df.loc[df.k<0.005]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a497c897-d0ed-4dc4-97e8-54d72c11b51f",
   "metadata": {},
   "outputs": [],
   "source": [
    "xy_big= bigk.groupby('x0')['y0']\n",
    "xy_small= smallk.groupby('x0')['y0']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e5eed24-9202-43d1-82f6-3e7f64053892",
   "metadata": {},
   "outputs": [],
   "source": [
    "xy_big.get_group(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72e7ac12-b59b-4009-86e6-e54b1b1f794b",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.errorbar(xy_big.mean().index, xy_big.mean(),yerr=xy_big.sem(), fmt='.', capsize=3)\n",
    "plt.errorbar(xy_small.mean().index, xy_small.mean(),yerr=xy_small.sem(), fmt='.', capsize=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92ada773-2be6-4fd1-aafd-32e0c9d1ff4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.hist(df['x0'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bac30f8-ea40-44bd-9c90-514b73e8340f",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "counts, bin, patches = plt.hist((srs_params['a']*srs_params['k']).values, bins=100, range=(0,0.5))\n",
    "plt.title(\"Net Glucose Uptake Rates $\\\\left( A\\cdot k\\\\right)$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c16e0a7-22e3-4180-8745-b59768ae3ed8",
   "metadata": {},
   "outputs": [],
   "source": [
    "params_by_age = df.join(first.reset_index('T')['T']).groupby('T').mean()\n",
    "params_by_age_sem = df.join(first.reset_index('T')['T']).groupby('T').sem()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e672bb18-9053-4d59-9383-9e8d86a13fdf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig, ax = plt.subplots(1,2)\n",
    "# ax[0].errorbar(t_data[:22], params_by_age['k'], params_by_age_sem['k'], fmt='o', capsize=3)\n",
    "# ax[0].set_ylabel(\"k\")\n",
    "# ax[0].set_xlabel(\"Time of birth\")\n",
    "# plt.semilogy();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72105df2-a103-4322-903f-003a08564582",
   "metadata": {},
   "outputs": [],
   "source": [
    "bkgd_avgs = srs_bsub.where(ds.labels==0).mean(list('YX'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c81daa10-f538-42a0-9e15-035f344e5e98",
   "metadata": {},
   "outputs": [],
   "source": [
    "bkgd_stds = srs_bsub.where(ds.labels==0).std(list('YX'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8a19931-d357-48c3-a782-856ab39593e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(bkgd_avgs.min().item(), bkgd_avgs.max().item())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4012604-b841-4c51-8ba2-332afc3322dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(bkgd_stds.min().item(), bkgd_stds.max().item())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98bdab53-23a3-47c5-bdc1-9d56b95dff6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "bad = srs_recon.loc[(srs_recon<-1*bkgd_stds.mean().item()).any(axis=1)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b896546-4cbc-418d-9a10-1d09f7f4e290",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(bad.values.T, 'k', alpha=0.2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "046db825-5114-4e1e-9a48-23024d4c8e85",
   "metadata": {},
   "outputs": [],
   "source": [
    "bad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b46a83ec-1f5c-4c83-bfe2-bd6e35488ece",
   "metadata": {},
   "outputs": [],
   "source": [
    "srs_params.join(pd.Series(index=bad.index, name='x'), how='inner')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9338827c-8e62-47db-a248-ca700d82a645",
   "metadata": {},
   "outputs": [],
   "source": [
    "srs_recon.loc[srs_recon[0]<-10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33762146-83d6-44e3-89fa-ae247ee138b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = check_labels_from_multiindex(ds.labels, bad.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00c7afa6-c154-415d-aa3a-02d6b4ff0c03",
   "metadata": {},
   "outputs": [],
   "source": [
    "v = napari.Viewer()\n",
    "sl = v.add_image(srs_bsub)\n",
    "fl = v.add_image(aligned_fluo)\n",
    "ml = v.add_labels(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd8f87e0-d103-403c-942a-6d85d8914bde",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = srs_params.loc[(sig_k_sel&sig_a_sel)]\n",
    "traces = srs_recon.loc[params.index]\n",
    "fig, ax = plt.subplots(1,5, sharey=True)\n",
    "q = params['k'].quantile(np.linspace(0,1,6))\n",
    "for i,(r,x) in enumerate(traces.groupby(pd.cut(params['k'], q.values))):\n",
    "    print(i)\n",
    "    ax[i].plot(x.values.T, 'k', alpha=0.05);\n",
    "    ax[i].set_title(f\"{0.2*i:0.1f}-{0.2*(i+1):0.1f}\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26d73215-3170-4dac-bf55-065ad54d9e07",
   "metadata": {},
   "outputs": [],
   "source": [
    "nanbkgd = srs_bsub.where(ds.labels==0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "057c2ed0-3a1a-40fc-b3f6-5c72deaa1eec",
   "metadata": {},
   "outputs": [],
   "source": [
    "be.cv_labels.load();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0df307db-1320-405d-a173-c21046dd44a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "bkgd_stds = mu.single_cell.standard_dev(be.cv_labels.to_dataset(name='labels'), srs_bsub)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "763cb5e6-f990-4f32-b440-9f9c106309df",
   "metadata": {},
   "outputs": [],
   "source": [
    "bkgd_sds = bkgd_stds.to_series().dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad0dd756-1993-4f74-874f-8b1c0457c9e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = srs_params.loc[(sig_k_sel&sig_a_sel)]\n",
    "traces = srs_recon.loc[params.index].drop((0,96))\n",
    "fig, ax = plt.subplots(figsize=(6.52,3))\n",
    "ax.set_aspect(9, anchor='C');\n",
    "pct = np.linspace(0,1,6)\n",
    "q = (params['k']).quantile(pct)\n",
    "inds = []\n",
    "for i,(r,x) in enumerate(traces.groupby(pd.cut(params['k'], q.values))):\n",
    "    print(i)\n",
    "    inds.append(x.index)\n",
    "    ax.errorbar(t_data, x.mean().values, yerr=x.sem(), fmt='-o', capsize=3, markersize=3, label=f\"{100*pct[i]:0.0f}-{100*pct[i+1]:0.0f}%\");\n",
    "    # ax[i].set_title(r);\n",
    "ax.legend(title=\"Percentiles of $k$\")\n",
    "plt.title(\"Net glucose uptake by $k$-Percentiles\")\n",
    "plt.xlabel(\"Time (minutes)\")\n",
    "plt.ylabel(\"SRS (a.u.)\")\n",
    "# plt.savefig(fig_path+\"kquantiles.png\")\n",
    "# m = bkgd_avgs.groupby('T').mean().values#nanbkgd.mean(list(\"SYX\")).data\n",
    "# s = bkgd_avgs.groupby('T').std().values#nand.std(list(\"SYX\")).data\n",
    "# plt.fill_between(t_data, m+s, m-s, color='k', alpha=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2386bc37-8575-4352-86e5-78d77e65df42",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "params = srs_params.loc[(sig_k_sel&sig_a_sel)]\n",
    "df = srs_recon.loc[params.index]\n",
    "# xdata = dt_min*df.columns.values + 30\n",
    "y = df.mean().values\n",
    "y_offset = y[0]\n",
    "y = y-y_offset\n",
    "sigma = df.sem().values\n",
    "\n",
    "fit = curve_fit(exp_approach, t_data, y, sigma=sigma, p0=[20, 0.05], bounds=(0,np.inf))\n",
    "\n",
    "a,k = fit[0]\n",
    "eqn = f\"Exponential fit: \\n{a:0.1f} (1-exp[-{k:0.2g}t])\"\n",
    "xplt = np.linspace(0, 550, 250)\n",
    "plt.figure()\n",
    "plt.plot(xplt, exp_approach(xplt, *fit[0])+y_offset, zorder=0, linewidth=3, alpha=0.75, label=eqn)\n",
    "# plt.plot(t_data, y+y_offset, 'k.', zorder=-1)\n",
    "plt.errorbar(t_data, y+y_offset,yerr=sigma, fmt='k.', capsize=3, zorder=-1)\n",
    "plt.title(f\"Wild Type Net Glucose Uptake\")\n",
    "plt.legend()\n",
    "plt.ylabel(\"SRS Intensity (a.u.)\")\n",
    "plt.xlabel(\"Time (minutes)\")\n",
    "if savefigs:\n",
    "    plt.savefig(fig_path+\"bulk_fit.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "102b4c42-eeff-4cb9-bd23-d8602ff53b37",
   "metadata": {},
   "outputs": [],
   "source": [
    "# m = check_labels_from_multiindex(ds.labels,inds[0])\n",
    "\n",
    "# v = napari.Viewer()\n",
    "# sl = v.add_image(srs_bsub)\n",
    "# ml = v.add_labels(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6090c55c-ce2c-48a2-9bc3-91173e46166c",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"mean\", srs_params['k'].mean())\n",
    "print(\"median\", srs_params['k'].median())\n",
    "print(\"selection mean\", srs_params.loc[sig_k_sel&sig_a_sel,'k'].mean())\n",
    "print(\"inv var mean\", (srs_params.loc[sig_k_sel&sig_a_sel,'k']*w).sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b7568ca-eadb-4541-b74f-1e6357d474b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,2, figsize=(6.52, 2.5))\n",
    "df = srs_params.loc[~(sig_k_sel&sig_a_sel)]\n",
    "hranges = ([0,0.05], [0,40])\n",
    "for i,x in enumerate('ka'):\n",
    "    # w = 1./df[f'sig_{x}']\n",
    "    # w /= w.sum()\n",
    "    \n",
    "    ax[i].hist(df[x].values, bins=100, range=hranges[i])\n",
    "    ax[i].axvline(df[x].median(), color='k', alpha=0.5, label=f'Median = {df[x].median():0.3g}')\n",
    "    ax[i].legend()\n",
    "ax[0].set_xlabel(r\"Uptake Rate  $\\left(min^{-1}\\right)$\")\n",
    "ax[1].set_xlabel(r\"Amplitude $\\left( a.u. \\right)$\")\n",
    "\n",
    "ax[0].set_ylabel(\"Cell counts\")\n",
    "\n",
    "fig.suptitle(\"Estimated Single-Cell Metabolic Parameters\")\n",
    "if savefigs:\n",
    "    plt.savefig(fig_path+\"bad_sc_params.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33cc63de-db3e-463b-a6fc-6560fced6b2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = srs_params.loc[(sig_k_sel&sig_a_sel)]\n",
    "p = df['a']*df['k']\n",
    "plt.figure()\n",
    "# x = p.loc[p>0.0002]\n",
    "plt.hist(p.values, bins=100, range=(0,0.6));\n",
    "print(p.median())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4aa46ea7-8785-4798-a380-dd602d74f2e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "big_k = srs_params.loc[srs_params['k']>0.02].index\n",
    "big_k_traces = srs_recon.loc[big_k]\n",
    "n_big_k = big_k_traces.shape[0] - big_k_traces.isna().sum()\n",
    "small_k = srs_params.loc[srs_params['k']<0.0025].index\n",
    "small_k_traces = srs_recon.loc[small_k]\n",
    "n_small_k = small_k_traces.shape[0] - small_k_traces.isna().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "242a8e0a-fce8-4b68-bc42-1eb6e402a719",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.errorbar(t_data, big_k_traces.mean(), yerr = big_k_traces.std(), capsize=3, fmt='-o')\n",
    "plt.errorbar(t_data, small_k_traces.mean(), yerr = small_k_traces.std(), capsize=3, fmt='-o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "addc6f66-2213-461f-8976-e81969423b14",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(t_data, n_big_k)\n",
    "plt.plot(t_data, n_small_k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fceecc9f-1005-4d68-bf64-c2b9186f0d4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,3,sharex=True, sharey=True)\n",
    "ax[0].plot(big_k_traces.values.T, 'k', alpha=0.05)\n",
    "ax[0].set_title(\"Big k\")\n",
    "ax[1].plot(small_k_traces.values.T,'k', alpha=0.05) \n",
    "ax[1].set_title(\"Small k\")\n",
    "ax[2].plot(srs_clean.loc[srs_params['k'].loc[srs_params['k'].between(0.0025, 0.2)].index].values.T, 'k', alpha=0.05)\n",
    "ax[2].set_title(\"Medium k\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac55cc43-b98a-4627-9009-31f8c0fb101e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from srs_tools.util import check_labels_from_multiindex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f64c8745-2b44-46e2-a80b-d2b2b90e134f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# big_k_labels = check_labels_from_multiindex(ds.labels, big_k)\n",
    "# small_k_labels = check_labels_from_multiindex(ds.labels, small_k)\n",
    "\n",
    "# v = napari.Viewer()\n",
    "# il = v.add_image(aligned_fluo)\n",
    "# big_l = v.add_labels(big_k_labels>0, color={1:'red'}, name=\"big k\", opacity=0.5)\n",
    "# small_l = v.add_labels(small_k_labels>0, color={1:'blue'}, name='small k', opacity=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24442f8a-7b5a-4b5d-bd95-6519a34364ee",
   "metadata": {},
   "source": [
    "## Age dependence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3792d01c-0824-4a45-8ff0-d84233dc5ef4",
   "metadata": {},
   "outputs": [],
   "source": [
    "tser = first.reset_index('T')['T']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d65872a-a1e4-45b4-be19-8da86c1144c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_t0 = all_first.reset_index('T')['T']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "463c7efb-e0a2-4378-acc6-031c0f1fc62a",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "counts, bins, patches = plt.hist(all_recon.loc[all_t0<5,35], bins=80, range=(0,40),density=True, alpha=0.6)\n",
    "plt.hist(all_recon.loc[all_t0.between(20,25),35], bins=bins, density=True, alpha=0.6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2051b7dd-6920-4812-aa20-f538a6d66263",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "xold= srs_recon.loc[tser<5, 35].values\n",
    "xnew = srs_recon.loc[tser.between(15,21), 35].values\n",
    "counts, bins, patches = plt.hist(xold, bins=80, alpha=0.6, density=True, label=\"Oldest\")\n",
    "plt.hist(xnew, bins=bins, alpha=0.6, density=True, label=\"Newest\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36abccb6-f8fb-41ab-ad59-96503fd65d48",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.quantile(xnew[~np.isnan(xnew)], [0.25, 0.75]))\n",
    "print(np.quantile(xold[~np.isnan(xold)], [0.25, 0.75]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b14d7fc4-8167-4132-bdf6-1a6084af1e33",
   "metadata": {},
   "outputs": [],
   "source": [
    "nold = xold.shape[0]-np.isnan(xold).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88acdc1c-8f8b-4fa8-8efc-40df16a4ac9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "nnew = xnew.shape[0]-np.isnan(xnew).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d83ba14b-90da-4dca-b69a-eaf4d3a58a55",
   "metadata": {},
   "outputs": [],
   "source": [
    "for z in [xold , xnew]:\n",
    "    denom = np.sqrt(z.shape[0]-np.isnan(z).sum())\n",
    "    print(np.nanmean(z), np.nanmedian(z), np.nanstd(z)/denom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95bafa82-c0eb-46cc-ab83-4901fd4dcc65",
   "metadata": {},
   "outputs": [],
   "source": [
    "Tmax = 36\n",
    "ngroups = 3\n",
    "width = (Tmax - nconsec)/ngroups\n",
    "bins = width*np.arange(ngroups+1)-0.5\n",
    "cut = pd.cut(tser, bins)\n",
    "gb = srs_recon.groupby(cut.cat.codes)\n",
    "by_age=gb.mean()\n",
    "by_age_sem=gb.sem()\n",
    "by_age_count = gb.count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "244c9a58-2ca9-4d97-80b6-14fe26686e3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.figure()\n",
    "# plt.plot(srs_recon.loc[all_srs_params['sig_k']==0].values.T, 'k', alpha=0.1);\n",
    "\n",
    "# plt.figure()\n",
    "# plt.plot(srs_recon.loc[all_srs_params['sig_k']==0].values.T, 'k', alpha=0.1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7de5d2f-cd29-4afc-9424-d811ccfa5266",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "for t in range(ngroups):\n",
    "    plt.errorbar(t_data, by_age.loc[t].values, yerr=gb.sem().loc[t].values,fmt='.', capsize=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00bc1f9e-19cf-477a-95a0-104b9d470366",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "params = pd.DataFrame(index=range(ngroups),columns=['a','k','t'])\n",
    "for t in range(ngroups):\n",
    "    y = by_age.loc[t].values\n",
    "    mask = ~np.isnan(y)\n",
    "    x = t_data[mask]\n",
    "    y = y[mask]\n",
    "    sigma = by_age_sem.loc[t].values[mask]\n",
    "    fit = curve_fit(exp_approach2, x, y, p0=[20, 0.01,50], bounds=(0,np.inf))\n",
    "    params.loc[t] = fit[0]\n",
    "    plt.plot(x,y, '.', color=colors[t])\n",
    "    plt.plot(xplt ,exp_approach2(xplt, *fit[0]), color=colors[t], label=\" - \".join([f\"{f:0.3f}\" for f in fit[0]]))\n",
    "plt.legend(title=\"a       -     k     -     t\")\n",
    "plt.axhline(0,color='k', linestyle=':')\n",
    "plt.ylim(-5,25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02227c2b-7c94-4b46-ba22-6549a3ee0aef",
   "metadata": {},
   "outputs": [],
   "source": [
    "xplt = np.linspace(t_data.min(), t_data.max(), 101)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5933ba3f-4e7a-46f8-a59a-b703d8b1a93e",
   "metadata": {},
   "outputs": [],
   "source": [
    "dt.mean().total_seconds()/60"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc8894a2-e90d-4c3b-b02b-cd4a9f5bbb19",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_bins = dt.mean().total_seconds()/60*(bins+0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "043555f7-0a45-44cd-b98b-0c7291670060",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "params = pd.DataFrame(index=range(ngroups),columns=['a','k'])\n",
    "fits = []\n",
    "for t in range(ngroups):\n",
    "    y = by_age.loc[t].values\n",
    "    mask = ~np.isnan(y)\n",
    "    x = t_data[mask]\n",
    "    x_offset = x[0]\n",
    "    x = x - x_offset\n",
    "    y = y[mask]\n",
    "    y_offset = y[0]\n",
    "    y = y-y_offset\n",
    "    fit = curve_fit(exp_approach, x, y, sigma=by_age_sem.loc[t][mask], p0=[20, 0.01], bounds=(0,np.inf))\n",
    "    fits.append(fit)\n",
    "    plt.errorbar(x+x_offset,y+y_offset, yerr=by_age_sem.loc[t][mask], fmt='.', capsize=2, color=colors[t])\n",
    "    xx = xplt[xplt>=x_offset]\n",
    "    plt.plot(xx , exp_approach(xx-x_offset, *fit[0])+y_offset, color=colors[t], label=f\"{t_bins[t]:0.1f} ≤\"+\" $T_{Birth}$ \"+f\"< {t_bins[t+1]:0.1f} \\n {fit[0][0]:0.1f}(1-exp[-{fit[0][1]:0.3g}])\", alpha=0.6)\n",
    "# plt.legend(title=\"    A --- k --- $y_0$\")\n",
    "plt.xlabel(\"Time (minutes)\")\n",
    "plt.ylabel(\"SRS Intensity (a.u.)\")\n",
    "plt.title(\"Net glucose uptake in cells grouped by age\")\n",
    "plt.legend()\n",
    "if savefigs: plt.savefig(fig_path+\"uptake_by_age.png\")\n",
    "# plt.axhline(0,color='k', linestyle=':')\n",
    "# plt.ylim(0,25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcd6bf7f-3702-453b-9127-1b7ce747ef98",
   "metadata": {},
   "outputs": [],
   "source": [
    "for f in fits:\n",
    "    print(np.product(f[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "511fcb9b-3dc7-4fd5-91aa-b15b375ad318",
   "metadata": {},
   "outputs": [],
   "source": [
    "k = [f[0][1] for f in fits]\n",
    "s = [np.sqrt(f[1][1,1]) for f in fits]\n",
    "plt.figure()\n",
    "plt.errorbar(t_bins[:-1] +np.mean(t_bins[:2]), k, yerr=s, fmt='-o', capsize=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd10f689-855e-41ef-8633-dfb7766c39c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "codes = cut.cat.codes[cut.cat.codes>-1].rename(\"codes\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b2fbc70-09af-4b60-857c-1b7f6054fd7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "ps = pd.concat([srs_params.loc[sig_k_sel&sig_a_sel,'k'], codes], axis=1, join='inner')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f0c7632-fc01-4c63-93d9-dfc32c4c94d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "for i in range(3):\n",
    "    x = ps.loc[ps['codes']==i, 'k']\n",
    "    plt.hist(x.values, bins=51, range=(0,0.05), density=True, alpha=0.6)\n",
    "    plt.axvline(x.median(), color=colors[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d9aec3b-80b5-4c24-91b1-cc6b136d2ae3",
   "metadata": {},
   "source": [
    "$$ y = a(1-e^{-k(t-t_0)})$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c85cd4f6-31e6-406a-b540-058afd316b35",
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_count = (~srs_clean.isna()).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aef7b5f2-2298-438a-8314-63be66597936",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_count = srs_avgs.groupby('T').count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94d74d6f-30ca-4f27-bd7d-4508abea69c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "lr = linregress(t_data, np.log2(all_count.values))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f721ebe-2489-4409-9f97-4534b954ebc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "lr2 = linregress(t_data[:20], np.log2(all_count.values)[:20])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc42bc3a-8fa0-432a-a693-96a6b6d598b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(1/lr2.slope)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4a8513b-a4c9-4f24-b78f-72134b6e9ff0",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(1/lr.slope)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "accbda4d-335f-4a5f-bbe5-fb6a1b12ba4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(t_data, cell_count.values, '-.')\n",
    "plt.plot(t_data, all_count.values.astype(int), '.')\n",
    "plt.plot(t_data, np.exp(np.log(2)*(lr.slope*t_data + lr.intercept)), 'k')\n",
    "plt.semilogy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84783e54-8f10-4b16-8bd4-858a2da66ebe",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "params = pd.DataFrame(index=range(5),columns=['a','k'])\n",
    "for t in range(3):\n",
    "    y = by_age.loc[t].values\n",
    "    mask = ~np.isnan(y)\n",
    "    x = t_data[mask]\n",
    "    y = y[mask]\n",
    "    sigma = by_age_sem.loc[t].values[mask]\n",
    "    fit = curve_fit(exp_approach, x, y, p0=[20, 0.01], bounds=(0,np.inf))\n",
    "    params.loc[t] = fit[0]\n",
    "    plt.plot(x,y, '.', color=colors[t])\n",
    "    plt.plot(xplt ,exp_approach(xplt, *fit[0]), color=colors[t], label=\" - \".join([f\"{f:0.3f}\" for f in fit[0]]))\n",
    "plt.legend(title=\"    a       -     k   \")\n",
    "plt.axhline(0,color='k', linestyle=':')\n",
    "plt.ylim(-5,25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f714a22-5448-42cf-a988-197053ca552d",
   "metadata": {},
   "outputs": [],
   "source": [
    "q = times.groupby('T').quantile([0.025, 0.975])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72cb2ecc-a804-4499-8c5a-62b27f75b650",
   "metadata": {},
   "outputs": [],
   "source": [
    "errs = q.unstack(level=-1).values.T"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98f49c81-f757-4967-825c-77e2fc2227d0",
   "metadata": {},
   "source": [
    "## Area growth rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a2faf25-dfb8-4500-bb11-2121ec803ce5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# slower to use dask delayed by ~50%\n",
    "# df = growth_cells\n",
    "\n",
    "# all_area_params = pd.DataFrame(index=df.index, columns=['a','k', 't', 'sig_a', 'sig_k', 'sig_t'], dtype='f4')\n",
    "\n",
    "# for i, (idx, s) in enumerate(tqdm(df.iterrows(), total=len(df))):\n",
    "#     y = s.values\n",
    "#     mask = ~np.isnan(y)\n",
    "#     y = y[mask]\n",
    "#     x = t_data[mask]\n",
    "#     out = fit_exp(x,y,[200,0.05,x[0]])\n",
    "#     all_area_params.loc[idx] = out\n",
    "\n",
    "# # unconverged = all_params.loc[all_params.isna().any(1)].index\n",
    "# area_params = all_area_params.dropna()\n",
    "\n",
    "# combo_rates = pd.concat([area_params['k'].rename('k_srs'), srs_params['k'].rename('k_area')],axis=1, join='inner')\n",
    "\n",
    "# combo_rates.corr(method='spearman') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f6a956b-94c7-4f9b-979c-313755524662",
   "metadata": {},
   "outputs": [],
   "source": [
    "# srs_total_mean = srs_bsub.where(ds.labels>0).mean(list('YX'))\n",
    "# srs_scaled_avg = (1+kopt.values)*srs_total_mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3cd938ed-167f-486d-b368-2b6e305c0aaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "srs_scaled_avg = srs_recon.groupby('S').mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd17d19d-ffd7-4287-b132-4b57543f2268",
   "metadata": {},
   "outputs": [],
   "source": [
    "a0 = areas.unstack('T')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68e48fa5-1a7f-46bc-8096-45a2ded6ca29",
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_counts = (a0>0).groupby('S').sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "437ec5d5-7429-434f-9c7c-95061a55fd0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,2)\n",
    "area_slope = []\n",
    "srs_slope = []\n",
    "log = np.log\n",
    "for i in range(5):\n",
    "    # a = vol_unit*vfactor*(total_area.data[i]**(3./2))\n",
    "    a = cell_counts.values[i]\n",
    "    # a = total_area.data[i]\n",
    "    # a = volumes.values[i]\n",
    "    lr_a = linregress(t_data, log(a))\n",
    "    area_slope.append(lr_a.slope)\n",
    "    ax[0].plot(t_data, a, color=colors[i], label=f\"{lr_a.slope:0.3g}±{lr_a.stderr:0.3g}\")\n",
    "    # s = srs_total_mean.data[i]\n",
    "    # s = srs_scaled_avg.data[i]\n",
    "    s = srs_scaled_avg.values[i]\n",
    "    s = s-s.min()\n",
    "    \n",
    "    # s = srs_recon.loc[i].mean().values\n",
    "    # lr_s = linregress(t_data, np.log(s))\n",
    "    fit_s = fit_exp(t_data, s, [20,0.01])\n",
    "    # srs_slope.append(lr_s.slope)\n",
    "    srs_slope.append(fit_s[1])\n",
    "    # ax[1].plot(t_data, s, color=colors[i], label=f\"{lr_s.slope:0.3g}±{lr_s.stderr:0.3g}\")\n",
    "    ax[1].plot(t_data, s, color=colors[i], label=f\"{fit_s[1]:0.3g}±{fit_s[-1]:0.3g}\")\n",
    "for a in ax:\n",
    "    a.legend(title=\"Slope ± std. err. $\\\\left(min^{-1}\\\\right)$\")\n",
    "    a.semilogy()\n",
    "    a.set_xlabel(\"Time (minutes)\")\n",
    "ax[0].set_title(\"Number of cells (per FOV - log scale)\")\n",
    "ax[0].set_ylabel(\"Approximate Volume $\\\\left(\\mu m^3 \\\\right)$\")\n",
    "ax[1].set_title(\"Average Cellular SRS (per FOV)\")\n",
    "ax[1].set_ylabel(\"SRS intenstity (a.u.)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2631fc39-4227-4979-a1ac-f114bfa6e3f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "r=[]\n",
    "for x,y in zip(srs_slope,area_slope):\n",
    "    r.append(x/y)\n",
    "    print(x, y, x/y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd8944aa-5edd-4b59-9d5e-3dc12282cfa6",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdcc689d-d91a-4598-a040-68775dbdf53c",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a7c6d55-7b41-4172-baae-91616c1bbefd",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.std(r)/np.sqrt(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b94fcb9d-945d-4bc0-aab1-521e033eddcd",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"{np.mean(r):0.3g} ± {np.std(r)/np.sqrt(5):0.3g}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81899ea5-13b4-442b-82fa-12db9ac0fc58",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(srs_scaled_avg.values.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ecf04b4-6ae9-4a45-b73c-d723800dad19",
   "metadata": {},
   "outputs": [],
   "source": [
    "outfile = \"/Users/johnrussell/Data/2023-05-16/dh229_srs_tables.h5\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "040dbef3-e749-402f-9623-bc2ca4b29a49",
   "metadata": {},
   "outputs": [],
   "source": [
    "srs_clean.to_hdf(outfile, key='raw')\n",
    "srs_recon.to_hdf(outfile, key='recon')\n",
    "areas_clean.to_hdf(outfile, key='area')\n",
    "srs_params.to_hdf(outfile, key='params')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}