magic fitting workflow.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Magic ratio AFM force curve analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ipywidgets\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.style\n",
    "matplotlib.style.use(['seaborn-v0_8-notebook','seaborn-v0_8-whitegrid','seaborn-v0_8-colorblind'])\n",
    "# from scipy.optimize import curve_fit, root_scalar\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matplotlib.rcParams['lines.markersize']=3\n",
    "matplotlib.rcParams['figure.dpi']=96\n",
    "matplotlib.rcParams['lines.linewidth']=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from magic_afm.calculation import secant, mylinspace, curve_fit, brentq, secant"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load and clean up data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from magic_afm.data_readers import mmap_path_read_only, ARDFFile\n",
    "\n",
    "# filename=r\"X:\\Data\\AFM\\Cypher\\2019-10-25\\PDMSW0006.ARDF\"\n",
    "# filename=r\"X:\\Data\\AFM\\Cypher\\2019-10-15\\TMPMP0006.ARDF\"\n",
    "# filename=r\"X:\\Data\\AFM\\Cypher\\2019-07-25\\magic0001.ARDF\"\n",
    "filename = r\"example.data\\Image0002.ARDF\"\n",
    "fvfile = ARDFFile.parse(mmap_path_read_only(filename))\n",
    "fs = 1/fvfile.t_step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fvfile.volumes[0].get_all_curves().shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xz, Xd = np.moveaxis(fvfile.volumes[0].get_all_curves(),2,0)\n",
    "\n",
    "# # Reshape to \"dataset\" (samples,features) and generate views of deflection and displacement\n",
    "im_r, im_c, *_ = Xz.shape\n",
    "Xz=np.reshape(Xz,(im_r*im_c,-1))\n",
    "Xd=np.reshape(Xd,(im_r*im_c,-1))\n",
    "X_cols = Xz.shape[-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Resample dataset with linear or Fourier interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Asylum default filter settings result in a cutoff frequency at 250xA_rate in the deflection channel and 25xA_rate in the z channel. Nyquist sampling would therefore be 500 points per indent, but 512 is better for FFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from magic_afm.calculation import resample_dset\n",
    "npts=512"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualize filter effect in power density spectrum\n",
    "i=np.random.randint(len(Xd),size=1)\n",
    "# j=np.random.randint(im_c,size=1)\n",
    "z=Xz[i]#,j].reshape(-1)\n",
    "d=Xd[i]#,j].reshape(-1)\n",
    "zsmp=resample_dset(z,npts,fourier=False)\n",
    "# zsmp=np.interp(np.linspace(0,1,npts,False),np.linspace(0,1,z.shape[-1],False),z.squeeze())[None,:]\n",
    "dsmp=resample_dset(d,npts,fourier=True)\n",
    "\n",
    "from scipy.signal import periodogram\n",
    "fig,(a1,a2)=plt.subplots(1,2)\n",
    "opts = dict(window=('kaiser',6),detrend='linear',scaling='density')\n",
    "a1.loglog(*periodogram(z.squeeze(),fs,**opts),label='z')\n",
    "a1.loglog(*periodogram(zsmp.squeeze(),fs*npts/X_cols,**opts),label='z resampled')\n",
    "a1.legend()\n",
    "# print(zsmp[:,-20:].mean())\n",
    "a2.loglog(*periodogram(d.squeeze(),fs,**opts),label='d')\n",
    "fd, pdd_den=periodogram(dsmp.squeeze(),fs*npts/X_cols,**opts)\n",
    "rms_noise = (pdd_den[-20:].mean()*fs*npts/X_cols/2)**.5\n",
    "print(rms_noise)\n",
    "a2.loglog(fd, pdd_den,label='d resampled')\n",
    "a2.legend()\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# resample and regenerate views\n",
    "Xz=resample_dset(Xz,npts,fourier=True)\n",
    "Xd=resample_dset(Xd,npts,fourier=True)\n",
    "#X=np.concatenate((Xz,Xd),axis=1)\n",
    "#X_cols = X.shape[1]//2\n",
    "#Xz=X[:,:X_cols]\n",
    "#Xd=X[:,X_cols:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspect dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "im_r, im_c, Xz.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,6))\n",
    "n=5\n",
    "indices = np.random.randint(len(Xz),size=n)\n",
    "plt.plot(Xz[indices,:].T,Xd[indices,:].T,alpha=1/n,color='black')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remove nuisance parameters if necessary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #%%time\n",
    "\n",
    "# # for piezo nuisance parameters\n",
    "# t = np.linspace(0,2*np.pi,X_cols//2,endpoint=False)\n",
    "# z = -np.cos(t)\n",
    "# zprime = np.sin(t)\n",
    "# zoffset= np.ones_like(t)\n",
    "# zramp = np.linspace(0,1,X_cols//2,endpoint=False)\n",
    "# # zoffsetdelta = np.concatenate((np.zeros(X_cols//4),np.ones(X_cols//4)))\n",
    "# zrampdelta = np.concatenate((np.zeros(X_cols//4),np.linspace(0,.5,X_cols//4,endpoint=False)))\n",
    "\n",
    "# # remove nuisance from z channel\n",
    "# A=np.stack([z,zoffset,zramp],axis=-1)\n",
    "# fit, res, rank, s = np.linalg.lstsq(A,Xz.T,rcond=None)\n",
    "# znuisance=(A[:,1:]@fit[1:]).T\n",
    "# Xz-=znuisance\n",
    "\n",
    "# # remove nuisance from deflection channel\n",
    "# A=np.stack([zoffset,],axis=-1)\n",
    "# keep = np.ones(X_cols//2,dtype='bool')\n",
    "# keep[X_cols//6:-X_cols//6] = False\n",
    "# fit, res, rank, s = np.linalg.lstsq(A[keep,:],Xd.T[keep,:],rcond=None)\n",
    "# dnuisance=(A@fit).T\n",
    "# Xd-=dnuisance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model fitting: Schwarz contact "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from magic_afm.calculation import schwarz_red"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "for red_fc in np.linspace(-2,-1.5,4).tolist():\n",
    "    x0 = np.linspace((7*red_fc+8)/3,red_fc,1000,endpoint=False)\n",
    "    x1 = np.linspace(red_fc,2,1000)\n",
    "    y0 =schwarz_red(x0,red_fc,-1,0)\n",
    "    y1 =schwarz_red(x1,red_fc,1,0)\n",
    "    line=plt.plot(np.concatenate((y0,y1)),\n",
    "                  np.concatenate((x0,x1)),\n",
    "                  '--',\n",
    "                 )[0]\n",
    "    \n",
    "    # k=0 instability point (df/ddelta=inf)\n",
    "    ai=((red_fc+2)/3)**(1/3)\n",
    "    di=ai**2-4*(ai*(red_fc+2)/3)**.5\n",
    "    plt.plot(di,(7*red_fc+8)/3,marker='P',ls='',c=line._color)\n",
    "    # critical force (df/ddelta=0 trivial solution)\n",
    "#     plt.axhline(red_fc,ls=':',c=line._color) #\n",
    "    dc = schwarz_red(red_fc,red_fc,1,0)\n",
    "    plt.plot(dc,red_fc,marker='X',ls='',c=line._color) \n",
    "    # imaginary maximum unstable force (df/ddelta=0 other solution)\n",
    "#     plt.axhline(4*red_fc+6,ls=':',c=line._color)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpolate to locate spring snap off point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "w=np.linspace((di-dc)*.95+dc,dc,1003)\n",
    "plt.plot(w,np.interp(w,y0,x0))\n",
    "plt.plot(w,np.interp(w,y0,np.gradient(x0,y0)))\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fc = -300\n",
    "# ref_force = np.pi*gamma*radius\n",
    "# red_fc = fc/ref_force\n",
    "ref_force = fc/red_fc\n",
    "radius= 70\n",
    "K=.25\n",
    "ref_radius = (ref_force*radius/K)**(1/3)\n",
    "ref_depth = ref_radius*ref_radius/radius\n",
    "red_k=-fvfile.k/(ref_force)*(ref_depth)\n",
    "print(red_k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import root_scalar\n",
    "root_scalar(lambda x: np.interp(x,y0,np.gradient(x0,y0))-red_k, bracket=(di,dc), method='brentq')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## model snap-in with L-J type potential\n",
    "as Lin, D. C., Dimitriadis, E. K. and Horkay, F. (2007) ‘Robust Strategies for Automated AFM Force Curve Analysis—II: Adhesion-Influenced Indentation of Soft, Elastic Materials’, Journal of Biomechanical Engineering, 129(6), p. 904 but with a 9-3 LJ potential and finite k snap-in/out as KL, J. and JA, G. (1997) ‘An Adhesion Map for the Contact of Elastic Spheres’, Journal Of Colloid And Interface Science, 192(2), pp. 326–333. \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from magic_afm.calculation import lj_force, lj_gradient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta = np.linspace(-2,0,1000)\n",
    "force = lj_force(delta,-.6,red_fc,0,0)\n",
    "gradient = lj_gradient(delta[delta<-.5],-.6,red_fc,0,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(delta,force)\n",
    "plt.plot(delta[delta<-.5],gradient)\n",
    "plt.plot(delta,np.gradient(force,delta),'--')\n",
    "plt.ylim([-2,2])\n",
    "_=plt.plot(w,np.interp(w,y0,x0),'--')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from magic_afm.calculation import lj_limit_factor\n",
    "args = (-.6,red_k,0,red_k)\n",
    "print('reduced k:',red_k)\n",
    "print('min gradient:', lj_gradient(lj_limit_factor*args[0],*args))\n",
    "print('bracket:',(2*args[1]-args[-2],lj_limit_factor*args[0]-args[-2]))\n",
    "root_scalar(lj_gradient,args=args,bracket=(2*args[0]-args[-2],lj_limit_factor*args[0]-args[-2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate dimensionless extend and retract curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from magic_afm.calculation import red_extend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "red_extend(np.linspace(-2,2,10),-1.752,red_k*.01,-.6,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x=np.linspace(-5,2,10000)\n",
    "y=red_extend(x,-1.52,red_k*1,-.6,)\n",
    "# y=red_extend(x,-1.525,7.530771527978149,-10,10)\n",
    "plt.figure()\n",
    "plt.plot(x,y)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from magic_afm.calculation import red_retract"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "red_retract(np.linspace(-2,2,10),-1.752,red_k,-.33333,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y=red_retract(x,-1.52,red_k,-.6,)#4,0.1)\n",
    "# plt.figure()\n",
    "plt.plot(x,y)\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scale model to real-space dimensions and fit to the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from magic_afm.calculation import force_curve,delta_curve,schwarz_wrap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i=np.random.randint(len(Xz))\n",
    "# i=1889\n",
    "# i=(63-63)*64+1\n",
    "fudgefactork=1.0\n",
    "# z=np.mean(Xz, axis=0)\n",
    "# d=np.mean(Xd, axis=0)/np.sqrt(fudgefactork)\n",
    "z=Xz[i]\n",
    "d=Xd[i]/np.sqrt(fudgefactork)\n",
    "f=d*fvfile.k*fudgefactork\n",
    "delta = z-d\n",
    "split = npts//2\n",
    "extsl=slice(split)\n",
    "retsl=slice(split,None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(d)\n",
    "plt.axvline(split,ls=':')\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from magic_afm.calculation import rapid_forcecurve_estimate\n",
    "p0 = rapid_forcecurve_estimate(delta[retsl], f[retsl], radius)\n",
    "print(p0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def partial_force_curve(delta, K, fc,  delta_shift, force_shift, \n",
    "                        lj_delta_scale,):\n",
    "    return force_curve(red_retract, delta, fvfile.k*fudgefactork, radius, K, -fc, tau, delta_shift, force_shift, \n",
    "                       lj_delta_scale,)\n",
    "        \n",
    "tau = .5\n",
    "x,y = delta[retsl], f[retsl]\n",
    "beta,cov,info=curve_fit(partial_force_curve, x,y, p0=p0,# xtol=1e-9,ftol=1e-8,\n",
    "                   bounds=np.transpose((\n",
    "                        (0.0, np.inf),  # M\n",
    "                        (0.0, np.inf),  # fc\n",
    "                        # (0, 1),           # tau\n",
    "                        (np.min(delta), np.max(delta)),  # delta_shift\n",
    "                        (np.min(f), np.max(f)),  # force_shift\n",
    "                        (-6.0, 6.0),  # lj_delta_scale\n",
    "                   )),\n",
    "                  )#method='trf', verbose=2, jac='2-point')\n",
    "\n",
    "with np.printoptions(suppress=True,precision=3,linewidth=100):\n",
    "    print(('K', 'fc', 'delta_shift', 'force_shift', 'lj_delta_scale', ))\n",
    "    print(repr(np.stack((beta,np.sqrt(np.diag(cov)),100*np.sqrt(np.diag(cov))/beta))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# with matplotlib.style.context('seaborn-ticks'):\n",
    "#     plt.figure(figsize=(4,4))\n",
    "#     _=plt.imshow(cov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# confidence interval of predictions\n",
    "from magic_afm.calculation import RT_EPS\n",
    "xsort = np.linspace(x.min(),x.max(),1000)\n",
    "yprimes = np.empty_like(y,shape=(len(beta),len(xsort)))\n",
    "ypred = partial_force_curve(xsort, *beta,)\n",
    "for j in range(len(yprimes)):\n",
    "    thisbeta=beta.copy()\n",
    "    dbeta = RT_EPS*beta[j]\n",
    "    thisbeta[j]+=dbeta\n",
    "    thisy = partial_force_curve(xsort, *thisbeta,)\n",
    "    yprimes[j,:]= (thisy-ypred)/dbeta\n",
    "    \n",
    "yci = np.sqrt(np.diag(yprimes.T@cov@yprimes))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(4,3))\n",
    "plt.plot(x-beta[2],(y-beta[3]-beta[1]),'.')\n",
    "\n",
    "# plt.plot(delta[extsl],f[extsl],'.')\n",
    "# f2=partial_force_curve(xsort-beta[2], *p0,)\n",
    "# f3=partial_force_curve(x, *beta,) #ypred\n",
    "# plt.plot(xsort-beta[2], f2-beta[3],':')\n",
    "plt.plot(xsort-beta[2], (ypred-beta[3]-beta[1]),'-')\n",
    "# plt.fill_between(xsort-beta[2],ypred-beta[3]+3*yci,ypred-beta[3]-3*yci,color='black',alpha=0.3)\n",
    "# plt.axvline(beta[2])\n",
    "# plt.xlim(beta[2]-2.5,beta[2]+2.5)\n",
    "# plt.ylim(-110,-90)\n",
    "\n",
    "plt.xlabel(r'Indentation $\\delta$ (nm)')\n",
    "plt.ylabel(r'Force $F$ (nN)')\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(x,y,'.')\n",
    "\n",
    "plt.plot(delta[extsl],f[extsl],'.')\n",
    "f2=partial_force_curve(xsort, *p0,)\n",
    "# f3=partial_force_curve(x, *beta,) #ypred\n",
    "plt.plot(xsort, f2,':')\n",
    "plt.plot(xsort, ypred,'k:')\n",
    "plt.fill_between(xsort,ypred+3*yci,ypred-3*yci,color='black',alpha=0.3)\n",
    "plt.axvline(beta[2],ls=':')\n",
    "plt.xlabel(r'Indentation $\\delta$ (nm)')\n",
    "plt.ylabel(r'Force $F$ (nN)')\n",
    "# plt.xlim(beta[2]-2.5,beta[2]+2.5)\n",
    "# plt.ylim(-110,-90)\n",
    "#plt.savefig('test.svg')\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.xlabel(r'Indentation depth $\\delta$ (nm)')\n",
    "# plt.ylabel(r'Indentation force $F$ (nN)')\n",
    "# plt.title(' '.join(('F =', str(int(force_setpoint)),'nN')))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.gcf().set_size_inches(4,3)\n",
    "# plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(cov[2,2]**.5,np.interp(0,xsort,yci))\n",
    "z_true_surface = float(beta[2] + (partial_force_curve(beta[2], *beta,))/fvfile.k)\n",
    "z_true_surface_std=np.sqrt(cov[2,2]+np.interp(0,xsort,yci)**2)\n",
    "print(f'z_true_surface {z_true_surface:.3f} ± {z_true_surface_std:.3f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from magic_afm.calculation import calc_def_ind_ztru, FitMode\n",
    "print(beta,radius,fvfile.k,tau)\n",
    "print(calc_def_ind_ztru(d[retsl],beta,radius,fvfile.k,tau,FitMode.RETRACT))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## seems OK so try to fit a bunch of data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from magic_afm.calculation import fitfun\n",
    "from tqdm.notebook import tqdm "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fudgefactork=1.05\n",
    "radius = 30\n",
    "tau=0\n",
    "n=len(Xz)\n",
    "inds = np.arange(n)\n",
    "# inds = np.random.choice(len(X),size=n,replace=False)\n",
    "fudgefactors=fudgefactork**-0.5\n",
    "\n",
    "betas = np.empty((n,2,10)) # shape: datapoints, #sens + 1, #fitparms*2\n",
    "for i in tqdm(inds, smoothing=1/24/5, desc='Fitting progress'):\n",
    "    z=Xz[i]\n",
    "    d=Xd[i]/np.sqrt(fudgefactork)\n",
    "    d_fudge=d*fudgefactors\n",
    "    f=d*fvfile.k\n",
    "    f_fudge=d*fvfile.k*fudgefactork\n",
    "    delta = z-d\n",
    "    delta_fudge = z-d_fudge\n",
    "    betas[i,0,:5],betas[i,0,5:],*_=fitfun(delta,f,fvfile.k,radius,tau,FitMode.RETRACT)\n",
    "    betas[i,1,:5],betas[i,1,5:],*_=fitfun(delta_fudge,f_fudge,fvfile.k*fudgefactork,radius,tau,FitMode.RETRACT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "abs_sens=(betas[:,1,:]-betas[:,0,:])/(fudgefactork-1)/fvfile.k\n",
    "rel_sens=(betas[:,1,:]-betas[:,0,:])/(fudgefactork-1)/betas[:,0,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import trimboth\n",
    "histdata=trimboth(rel_sens[:,0].T,.05)\n",
    "histdata=histdata[~np.isnan(histdata)]\n",
    "_=(histdata.mean(),histdata.std(),np.median(histdata))\n",
    "print(*_,_[1]/_[0]*100)\n",
    "plt.figure()\n",
    "_=plt.hist(histdata, bins='auto', cumulative=False, density=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## calculate indentation ratio for each force curve from the data and fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from magic_afm import calc_def_ind_ztru\n",
    "# inds = np.random.choice(np.arange(len(X)),size=2,replace=False)\n",
    "# inds = np.arange(n)\n",
    "# radius,tau = 10,0\n",
    "deflections=np.full(inds.shape,np.nan)\n",
    "indentations=np.full(inds.shape,np.nan)\n",
    "z_true_surfaces=np.full(inds.shape,np.nan)\n",
    "for i in range(n):\n",
    "    if np.any(np.isnan(betas[i,0,:])):\n",
    "        continue\n",
    "    K, fc, delta_shift, force_shift, lj_delta_scale,*_ = betas[i,0,:]\n",
    "    z = Xz[i]\n",
    "    z = z[retsl]\n",
    "    d = Xd[i]\n",
    "    d = d[retsl]\n",
    "    delta=z-d\n",
    "    force = force_curve(red_retract, delta, fvfile.k, radius, K, -fc, tau, \n",
    "                               delta_shift, force_shift, lj_delta_scale,) # for plotting in next cells\n",
    "\n",
    "    deflections[i],indentations[i],z_true_surfaces[i],_=calc_def_ind_ztru(force,betas[i,0,:],radius,fvfile.k,tau,FitMode.RETRACT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('deflection,   indentation,   deflection/indentation')\n",
    "print(deflections[i],indentations[i], deflections[i]/indentations[i])\n",
    "delta=z-d\n",
    "with plt.style.context({'axes.titlesize':17,'axes.labelsize':15,'xtick.labelsize':12,'ytick.labelsize':12}):\n",
    "    plt.figure(figsize=(4*1.5,3*1.5))\n",
    "    _=plt.plot(delta-delta_shift,d-force_shift/fvfile.k,'.',alpha=.75)\n",
    "    _=plt.plot(delta-delta_shift,(force-force_shift)/fvfile.k,'k')\n",
    "    snapoff_offset=(delta-delta_shift)[np.argmin(force)]\n",
    "    _=plt.plot(snapoff_offset,fc/fvfile.k,'s')\n",
    "    _=plt.plot(indentations[i]+snapoff_offset,deflections[i]+fc/fvfile.k,'rX',ms=10)\n",
    "#     _=plt.xlim(-25,5)\n",
    "#     _=plt.ylim(bottom=-.5)\n",
    "    _=plt.xlabel(r'Indentation $\\delta$ (nm)')\n",
    "    _=plt.ylabel(r'Deflection $d$ (nm)')\n",
    "    _=plt.title(r'Example $d$ vs. $\\delta$ curve, shifted')\n",
    "_=plt.plot(0,(force_curve(red_retract, delta_shift, fvfile.k, radius, K, -fc, tau, \n",
    "                               delta_shift, force_shift, lj_delta_scale,)-force_shift)/fvfile.k,'rP',ms=10)\n",
    "# )CAN IT BE DONE for F<Fmax BY TRUNCATING DATA?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('deflection,   indentation,   deflection/indentation')\n",
    "# i=np.random.choice(np.arange(len(X)))\n",
    "print(deflections[i],indentations[i], deflections[i]/indentations[i], z_true_surfaces[i])\n",
    "plt.figure()\n",
    "_=plt.plot(z,d,'.')\n",
    "_=plt.plot(z,(force)/fvfile.k,'k')\n",
    "# _=plt.plot(z[i],deflections[i]+fc/k,'x')\n",
    "# _=plt.plot((z)[np.argmin(force)],(fc+force_shift)/k,'P')\n",
    "_=plt.plot(z_true_surfaces[i],(force_curve(red_retract, delta_shift, fvfile.k, radius, K, -fc, tau, \n",
    "                               delta_shift, force_shift, lj_delta_scale,))/fvfile.k,'rP',ms=10)\n",
    "_=plt.xlabel(r'Z displacement $z$ (nm)')\n",
    "_=plt.ylabel(r'Deflection w/offset $d$ (nm)')\n",
    "_=plt.title('Location of true surface visualized')\n",
    "# )CAN IT BE DONE for F<Fmax BY TRUNCATING DATA?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Put it all together in summary plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save above analysis in a dict according to filename\n",
    "saved = {}\n",
    "saved[filename]= deflections, indentations, betas, fvfile.k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sensfig,((sensratax,sensindax,sensdefax),\n",
    "         (modratax,modindax,moddefax))=plt.subplots(2,3,figsize=(9,6),sharex='col',sharey='row')\n",
    "thresh = .03\n",
    "s=10\n",
    "centroid=np.median\n",
    "poisson=0.350\n",
    "\n",
    "leg=[]\n",
    "for fn in saved:\n",
    "    if fn in {\n",
    "#         r\"X:\\Data\\AFM\\Cypher\\2019-07-25\\magic0010.ARDF\",\n",
    "#         r\"X:\\Data\\AFM\\Cypher\\2019-07-25\\magic0009.ARDF\",\n",
    "#         r\"X:\\Data\\AFM\\Cypher\\2019-07-25\\magic0002.ARDF\",\n",
    "#         r\"X:\\Data\\AFM\\Cypher\\2019-07-25\\magic0005.ARDF\",\n",
    "    }:\n",
    "        continue\n",
    "    deflections, indentations, betas, k = saved[fn]\n",
    "#     force_setpoint, z_rate=int(round(force_setpoint)), int(round(z_rate))\n",
    "    # if force_setpoint != 50:\n",
    "    #     pass\n",
    "    \n",
    "    ratios=(deflections)/indentations\n",
    "#     abs_sens=(betas[:,1,:]-betas[:,0,:])/(fudgefactork-1)/k\n",
    "    rel_sens=(betas[:,1,:]-betas[:,0,:])/(fudgefactork-1)/betas[:,0,:]\n",
    "    sens=rel_sens[:,0]\n",
    "    moduli = (1-poisson**2)*betas[:,0,0]/(4/3)\n",
    "\n",
    "    lo,hi = np.nanquantile(sens,(thresh,1-thresh))\n",
    "    keep = (lo < sens) & (sens < hi)\n",
    "    lo,hi = np.nanquantile(ratios,(thresh,1-thresh))\n",
    "    keep &= (lo < ratios) & (ratios < hi) #& (.1 < ratios)\n",
    "    alpha = min(100/sum(keep),1)\n",
    "#     alpha /= (np.pi*sigmas*sigmas)\n",
    "#     alpha /= alpha.max()*5\n",
    "#     c=np.array(to_rgb(line.get_color()))\n",
    "#     c=np.concatenate((np.broadcast_to(c,(len(sigmas),3)),alpha[:,np.newaxis]),axis=-1)\n",
    "\n",
    "    deflections = deflections[keep]\n",
    "    indentations = indentations[keep]\n",
    "    ratios = ratios[keep]\n",
    "    sens = sens[keep]\n",
    "    moduli = moduli[keep]\n",
    "    \n",
    "    # style settings\n",
    "    # could set color here instead of using default cycle\n",
    "    scatter_style = dict(s=s,alpha=alpha, linewidth=0) \n",
    "    centroid_style = dict(marker='x',color='black',zorder=3)\n",
    "    \n",
    "    # sensratax => sensitivity vs ratio axes\n",
    "    p=np.polyfit(ratios,sens,1)\n",
    "    x=np.linspace(lo,hi)\n",
    "    y=np.polyval(p,x)\n",
    "    line,=sensratax.plot(x,y,':')\n",
    "    # sensax.hist2d(ratios,sens,\n",
    "    #            bins=150,cmin=0,density=True,cmap='Greys',\n",
    "    #           range=((0,.7),(0,.5)),\n",
    "    # #           range=((0,.75),(0,.0175)),\n",
    "    #           )\n",
    "#     sigmas=betas[keep,0,5]/betas[keep,0,0]\n",
    "#     s = 300*sigmas\n",
    "    sensratax.scatter(ratios,sens,**scatter_style)\n",
    "    sensratax.scatter(centroid(ratios),centroid(sens),**centroid_style)\n",
    "#     sensax.set_xlabel(r'Indentation ratio $\\frac{d}{\\delta}$')\n",
    "    #     sensax.set_xlim((0,.7))\n",
    "    #     sensax.set_ylim((0,.5))\n",
    "    \n",
    "    # modratax => modulus vs ratio axes\n",
    "    p=np.polyfit(ratios,moduli,1)\n",
    "    x=np.linspace(lo,hi)\n",
    "    y=np.polyval(p,x)\n",
    "    modratax.plot(x,y,':')\n",
    "    modratax.scatter(ratios,moduli,**scatter_style)\n",
    "    modratax.scatter(centroid(ratios),centroid(moduli),**centroid_style)\n",
    "    \n",
    "    # sensindax => sensitivity vs indentation axes\n",
    "    lo,hi = np.nanquantile(indentations,(thresh,1-thresh))\n",
    "    p=np.polyfit(indentations,sens,1)\n",
    "    x=np.linspace(lo,hi)\n",
    "    y=np.polyval(p,x)\n",
    "    sensindax.plot(x,y,':')\n",
    "    sensindax.scatter(indentations,sens,**scatter_style)\n",
    "    sensindax.scatter(centroid(indentations),centroid(sens),**centroid_style)\n",
    "    \n",
    "    # modindax => modulus vs indentation axes\n",
    "    p=np.polyfit(indentations,moduli,1)\n",
    "    x=np.linspace(lo,hi)\n",
    "    y=np.polyval(p,x)\n",
    "    modindax.plot(x,y,':')\n",
    "    modindax.scatter(indentations,moduli,**scatter_style)\n",
    "    modindax.scatter(centroid(indentations),centroid(moduli),**centroid_style)\n",
    "    \n",
    "    # sensindax => sensitivity vs deflections axes\n",
    "    lo,hi = np.nanquantile(deflections,(thresh,1-thresh))\n",
    "    p=np.polyfit(deflections,sens,1)\n",
    "    x=np.linspace(lo,hi)\n",
    "    y=np.polyval(p,x)\n",
    "    sensdefax.plot(x,y,':')\n",
    "    sensdefax.scatter(deflections,sens,**scatter_style)\n",
    "    sensdefax.scatter(centroid(deflections),centroid(sens),**centroid_style)\n",
    "    \n",
    "    # modindax => modulus vs deflections axes\n",
    "    p=np.polyfit(deflections,moduli,1)\n",
    "    x=np.linspace(lo,hi)\n",
    "    y=np.polyval(p,x)\n",
    "    moddefax.plot(x,y,':')\n",
    "    moddefax.scatter(deflections,moduli,**scatter_style)\n",
    "    moddefax.scatter(centroid(deflections),centroid(moduli),**centroid_style)\n",
    "    \n",
    "    # leg.append((force_setpoint,z_rate))\n",
    "    \n",
    "\n",
    "# sensratax.relim()\n",
    "# modratax.relim()\n",
    "# sensindax.relim()\n",
    "# modindax.relim()\n",
    "# sensdefax.relim()\n",
    "# moddefax.relim()\n",
    "# sensratax.set_ylim(top=1.0)\n",
    "# modratax.set_ylim(bottom=-0.025,top=.375)\n",
    "# modratax.set_xlim(left=0.5)\n",
    "# modindax.set_xlim(left=0.0)\n",
    "\n",
    "# sensdefax.legend([f'{(force_setpoint)} nN, {(z_rate)} Hz' for force_setpoint,z_rate in leg],loc='upper right')\n",
    "sensratax.axhline(0,color='black',linestyle=':')\n",
    "sensindax.axhline(0,color='black',linestyle=':')\n",
    "sensdefax.axhline(0,color='black',linestyle=':')\n",
    "sensratax.set_ylabel(r'Rel. modulus sens. $\\frac{dE}{dk} \\times \\frac{k}{E}$')\n",
    "modratax.set_xlabel(r'Indentation ratio $\\frac{d_m}{\\delta_m}$')\n",
    "modratax.set_ylabel(r'Modulus $E$')\n",
    "modindax.set_xlabel(r'Indentation depth $\\delta_m$')\n",
    "moddefax.set_xlabel(r'Deflection $d_m$')\n",
    "sensfig.tight_layout()\n",
    "sensfig.align_labels()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sensfig.savefig('ratio vs sens.png',dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import probplot\n",
    "plt.figure()\n",
    "probplot(ratios,plot=plt)\n",
    "probplot(rel_sens[:,0],plot=plt)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# maxforces = Xd[:,:len(Xd[0])//25].mean(axis=-1)*k/\n",
    "# deflections, indentations, betas, k, force_setpoint, z_rate=saved[r\"X:\\Data\\AFM\\Cypher\\2019-07-25\\magic0003.ARDF\"]\n",
    "# histdata=trimboth(betas[:,0,0],.001)\n",
    "# print(np.nanmean(histdata),np.nanstd(histdata),len(histdata))\n",
    "# plt.figure()\n",
    "# _=plt.hist(histdata, bins='auto', cumulative=False, density=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "saved.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}