{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# fortran function linear4dto1dxx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as N\n",
    "from fortran._interp_ import linear4dto1dxx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### note : Use of N.ogrid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nti = 4\n",
    "# inclusive\n",
    "N.ogrid[0:nti-1:nti*1j]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N.ogrid[0:nti+1:nti]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N.ogrid[0:nti-1:nti]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic initializations for the test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %% Bases\n",
    "\n",
    "nex = 3\n",
    "nexz = 2\n",
    "nxi = 7\n",
    "nyi = 6\n",
    "nzi = 5\n",
    "nti = 4\n",
    "mv = 1.e20\n",
    "no = 1\n",
    "\n",
    "vfunc = lambda t, z, y, x: (1*x + 2.35*y + 3.24*z -0.65*t)\n",
    "\n",
    "# Create t,z,y,x axes over a 4D grid\n",
    "#1j way to use N.ogrid (inclusive)\n",
    "tti, zzi, yyi, xxi = N.mgrid[0:nti-1:nti*1j, 0:nzi-1:nzi*1j,\n",
    "    0:nyi-1:nyi*1j, 0:nxi-1:nxi*1j]\n",
    "print(\"Shapes of tti: {}, zzi: {}, yyi: {}, xxi: {}\".format(tti.shape,zzi.shape,yyi.shape,xxi.shape))\n",
    "\n",
    "# add one dimension\n",
    "print(\"zzi.shape after N.mgrid {}\".format(zzi.shape))\n",
    "zzi = zzi[None]\n",
    "print(\"zzi.shape after zzi[None] {}\".format(zzi.shape))\n",
    "zzi = N.repeat(zzi, nexz, axis=0)\n",
    "print(\"zzi.shape after N.repeat {}\".format(zzi.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nxi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xxi[0,0,0,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xxi[1,1,4,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "tti varies along 3th axis, same t values on z-y-x-dimensions  \n",
    "zzi varies along 2nd axis, same z values on t-y-x-dimensions  \n",
    "yyi varies along 1st axis, same y values on t-z-x-dimensions  \n",
    "xxi varies along 0th axis, same x values on t-z-y-dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N.testing.assert_equal(tti[3,...],tti[1,...]*3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N.testing.assert_equal(yyi[...,5,:],yyi[...,1,:]*5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N.testing.assert_array_equal(zzi[0,...],zzi[1,...])  # after repeat, i.e. extension to the left"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpolation from pure 1D axes for t,z,y and x axes to a unique location"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %% Pure 1D axes along each dimension\n",
    "xi = xxi[0, 0, 0:1, :] # (nyix=1,nxi)\n",
    "yi = yyi[0, 0, :, 0:1] # (nyi,nxiy=1)\n",
    "zi = zzi[0:1, 0:1, :, 0:1, 0:1] # (nexz=1,ntiz=1,nzi,nyiz=1,nxiz=1)\n",
    "ti = tti[:, 0, 0, 0] # (nti)\n",
    "print(\"Size of xi: {}, yi: {}, zi: {}, ti: {},\".format(xi.shape,yi.shape,zi.shape,ti.shape))\n",
    "\n",
    "# 4D data field\n",
    "vi = vfunc(tti, zzi, yyi, xxi)\n",
    "print(\"Size of vi: {} and zzi: {}\".format(vi.shape,zzi.shape))\n",
    "# extend the first dimension of vi\n",
    "vi = N.resize(vi, (nex, )+vi.shape[1:]) # (nex,nti,nzi,nyi,nxi)\n",
    "print(\"Size of vi after N.resize: {}\".format(vi.shape))\n",
    "\n",
    "# Get a unique random location\n",
    "N.random.seed(0) #makes the random numbers predictable\n",
    "xyztomin = -0.5\n",
    "# get one number (no=1) between -0.5 and dim-1.5. saved into a numpy.array\n",
    "xo = N.random.uniform(xyztomin, nxi-1.5, no)\n",
    "yo = N.random.uniform(xyztomin, nyi-1.5, no)\n",
    "zo = N.random.uniform(xyztomin, nzi-1.5, no)\n",
    "to = N.random.uniform(xyztomin, nti-1.5, no)\n",
    "\n",
    "# reference\n",
    "vo_truth = N.ma.array(vfunc(to, zo, yo, xo))\n",
    "\n",
    "# interpolation\n",
    "mv = 1e20\n",
    "vo_interp = linear4dto1dxx(xi,yi,zi,ti,vi,xo,yo,zo,to)\n",
    "print(\"vo_interp :\",vo_interp)\n",
    "#vo_interp = N.ma.masked_values(vo_interp), mv)\n",
    "\n",
    "# None if correct\n",
    "N.testing.assert_almost_equal(vo_interp[0], vo_truth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vo_interp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vo_truth"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ".. important:: Pourquoi 3 valeurs ??? Parce que 5ème dimension avant tzyx est égale à 3. Dimension non prise en compte dans interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### apparté"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vi.shape[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(nex, 9)+vi.shape[1:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpolation from a single point in space (x=cst, y=cst), z and t axes are pure 1D axes to a unique location"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %% Single point in space\n",
    "\n",
    "xi = xxi[0, 0, 0:1, :1] # (nyix=1,nxi)\n",
    "yi = yyi[0, 0, :1, 0:1] # (nyi,nxiy=1)\n",
    "zi = zzi[0:1, 0:1, :, 0:1, 0:1] # (nexz=1,ntiz=1,nzi,nyiz=1,nxiz=1)\n",
    "ti = tti[:, 0, 0, 0] # (nti)\n",
    "vi = vfunc(tti, zzi, yyi, xxi)[:, :, :, :1, :1]\n",
    "vi = N.resize(vi, (nex, )+vi.shape[1:]) # (nex,nti,nzi,1,1)\n",
    "print(\"Size of xi: {}, yi: {}, zi: {}, ti: {},\".format(xi.shape,yi.shape,zi.shape,ti.shape))\n",
    "print(\"Size of vi: {} and zzi: {}\".format(vi.shape,zzi.shape))\n",
    "\n",
    "N.random.seed(0)\n",
    "xyztomin = -0.5\n",
    "xo = N.random.uniform(xyztomin, nxi-1.5, no)\n",
    "yo = N.random.uniform(xyztomin, nyi-1.5, no)\n",
    "zo = N.random.uniform(xyztomin, nzi-1.5, no)\n",
    "to = N.random.uniform(xyztomin, nti-1.5, no)\n",
    "vo_truth = N.ma.array(vfunc(to, zo, yi[0], xi[0]))\n",
    "\n",
    "mv = 1e20\n",
    "vo_interp = linear4dto1dxx(xi,yi,zi,ti,vi,xo,yo,zo,to)\n",
    "#vo_interp = N.ma.masked_values(vo_interp, mv)\n",
    "print(\"vo_interp :\",vo_interp)\n",
    "\n",
    "N.testing.assert_almost_equal(vo_interp[0], vo_truth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vo_truth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N.testing.assert_equal(vo_interp[0], vo_truth)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## No interpolation in time but x, y and z axes are pure 1D axes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %% Constant time\n",
    "\n",
    "xi = xxi[0, 0, 0:1, :] # (nyix=1,nxi)\n",
    "yi = yyi[0, 0, :, 0:1] # (nyi,nxiy=1)\n",
    "zi = zzi[0:1, 0:1, :, 0:1, 0:1] # (ntiz=1,nzi,nyiz=1,nxiz=1)\n",
    "ti = tti[:1, 0, 0, 0] # (1)\n",
    "vi = vfunc(tti, zzi, yyi, xxi)[:, :1, :, :, :]\n",
    "vi = N.resize(vi, (nex, )+vi.shape[1:]) # (nex,1,nzi,nyi,nxi)\n",
    "\n",
    "N.random.seed(0)\n",
    "xyztomin = -0.5\n",
    "xo = N.random.uniform(xyztomin, nxi-1.5, no)\n",
    "yo = N.random.uniform(xyztomin, nyi-1.5, no)\n",
    "zo = N.random.uniform(xyztomin, nzi-1.5, no)\n",
    "to = N.random.uniform(xyztomin, nti-1.5, no)\n",
    "vo_truth = N.ma.array(vfunc(ti, zo, yo, xo))\n",
    "\n",
    "mv = 1e20\n",
    "vo_interp = linear4dto1dxx(xi,yi,zi,ti,vi,xo,yo,zo,to)\n",
    "#vo_interp = N.ma.masked_values(vo_interp, mv)\n",
    "\n",
    "N.testing.assert_almost_equal(vo_interp[0], vo_truth)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Size of xi: {}, yi: {}, zi: {}, ti: {}\".format(xi.shape,yi.shape,zi.shape,ti.shape))\n",
    "print(\"Size of vi: {} and zzi: {}\".format(vi.shape,zzi.shape))\n",
    "print(\"vo_truth : {}\".format(vo_truth))\n",
    "print(\"vo_interp :\",vo_interp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ti"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpolation from 1D x, y and t axes, z is a 5D array to a unique location"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %% Variable depth with 1D X/Y + T\n",
    "\n",
    "xi = xxi[0, 0, 0:1, :] # (nyix=1,nxi)\n",
    "yi = yyi[0, 0, :, 0:1] # (nyi,nxiy=1)\n",
    "zi = zzi[:, :, :, :, :] # (nexz,ntiz=nti,nzi,nyiz=nyi,nxiz=nxi)\n",
    "ti = tti[:, 0, 0, 0] # (nti)\n",
    "vi = vfunc(tti, zzi, yyi, xxi)\n",
    "vi = N.resize(vi, (nex, )+vi.shape[1:]) # (nex,nti,nzi,nyi,nxi)\n",
    "\n",
    "N.random.seed(0)\n",
    "xyztomin = -0.5\n",
    "xo = N.random.uniform(xyztomin, nxi-1.5, no)\n",
    "yo = N.random.uniform(xyztomin, nyi-1.5, no)\n",
    "zo = N.random.uniform(xyztomin, nzi-1.5, no)\n",
    "to = N.random.uniform(xyztomin, nti-1.5, no)\n",
    "vo_truth = N.ma.array(vfunc(to, zo, yo, xo))\n",
    "\n",
    "mv = 1e20\n",
    "vo_interp = linear4dto1dxx(xi,yi,zi,ti,vi,xo,yo,zo,to)\n",
    "#vo_interp = N.ma.masked_values(vo_interp, mv)\n",
    "\n",
    "N.testing.assert_almost_equal(vo_interp[0], vo_truth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Size of xi: {}, yi: {}, zi: {}, ti: {}\".format(xi.shape,yi.shape,zi.shape,ti.shape))\n",
    "print(\"Size of vi: {} and zzi: {}\".format(vi.shape,zzi.shape))\n",
    "print(\"vo_truth : {}\".format(vo_truth))\n",
    "print(\"vo_interp :\",vo_interp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xo,yo,zo,to"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ti"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Size of xi: {}, yi: {}, zi: {}, ti: {}\".format(xi.shape,yi.shape,zi.shape,ti.shape))\n",
    "print(\"Size of vi: {} and zzi: {}\".format(vi.shape,zzi.shape))\n",
    "print(\"vo_truth : {}\".format(vo_truth))\n",
    "print(\"vo_interp :\",vo_interp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ".. error:: Pourquoi nan pour dernier élément ???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpolation from 2D X/Y - no other axes (pure curvilinear) - to a unique location"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %% 2D X/Y with no other axes (pure curvilinear)\n",
    "\n",
    "xi = xxi[0, 0, ...] # (nyix=nyi,nxi)\n",
    "yi = yyi[0, 0, ...]# (nyi,nxiy=nxi)\n",
    "zi = zzi[0:1, 0:1, 0:1, 0:1, 0:1] # (nexz=1,ntiz=1,1,nyiz=1,nxiz=1)\n",
    "ti = tti[:1, 0, 0, 0] # (1)\n",
    "\n",
    "vi = vfunc(tti, zzi, yyi, xxi)[:, :1, :1, :, :]\n",
    "vi = N.resize(vi, (nex, )+vi.shape[1:]) # (nex,1,1,nyi,nxi)\n",
    "\n",
    "N.random.seed(0)\n",
    "xyztomin = -0.5\n",
    "xo = N.random.uniform(xyztomin, nxi-1.5, no)\n",
    "yo = N.random.uniform(xyztomin, nyi-1.5, no)\n",
    "zo = N.random.uniform(xyztomin, nzi-1.5, no)\n",
    "to = N.random.uniform(xyztomin, nti-1.5, no)\n",
    "vo_truth = N.ma.array(vfunc(ti, zi.ravel()[0], yo, xo))\n",
    "\n",
    "mv = 1e20\n",
    "# xi and yi are 2D\n",
    "vo_interp = linear4dto1dxx(xi,yi,zi,ti,vi,xo,yo,zo,to)\n",
    "# xi and yi are vectors - rectilinear grid\n",
    "vo_interp_rect = linear4dto1dxx(xi[:1],yi[:, :1],zi,ti,vi,xo,yo,zo,to)\n",
    "\n",
    "N.testing.assert_almost_equal(vo_interp[0], vo_truth)\n",
    "N.testing.assert_almost_equal(vo_interp_rect[0], vo_truth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Size of xi: {}, yi: {}, zi: {}, ti: {}\".format(xi.shape,yi.shape,zi.shape,ti.shape))\n",
    "print(\"Size of vi: {} and zzi: {}\".format(vi.shape,zzi.shape))\n",
    "print(\"vo_truth : {}\".format(vo_truth))\n",
    "print(\"vo_interp :\",vo_interp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpolation from and to the same coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %% Same coordinates\n",
    "\n",
    "xi = xxi[0, 0, 0:1, :] # (nyix=1,nxi)\n",
    "yi = yyi[0, 0, :, 0:1] # (nyi,nxiy=1)\n",
    "zi = zzi[0:1, 0:1, :, 0:1, 0:1] # (nexz=1,ntiz=1,nzi,nyiz=1,nxiz=1)\n",
    "ti = tti[:, 0, 0, 0] # (nti)\n",
    "vi = vfunc(tti, zzi, yyi, xxi)\n",
    "vi = N.resize(vi, (nex, )+vi.shape[1:]) # (nex,nti,nzi,nyi,nxi)\n",
    "\n",
    "xyzo = N.meshgrid(xi, yi, zi, ti, indexing='ij')\n",
    "xo = xyzo[0].ravel()\n",
    "yo = xyzo[1].ravel()\n",
    "zo = xyzo[2].ravel()\n",
    "to = xyzo[3].ravel()\n",
    "vo_truth = N.ma.array(vfunc(to, zo, yo, xo))\n",
    "\n",
    "mv = 1e20\n",
    "vo_interp = linear4dto1dxx(xi,yi,zi,ti,vi,xo,yo,zo,to)\n",
    "#vo_interp = N.ma.masked_values(vo_interp, mv)\n",
    "\n",
    "N.testing.assert_almost_equal(vo_interp[0], vo_truth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Size of xi: {}, yi: {}, zi: {}, ti: {}\".format(xi.shape,yi.shape,zi.shape,ti.shape))\n",
    "print(\"Size of vi: {} and zzi: {}\".format(vi.shape,zzi.shape))\n",
    "print(\"vo_truth : {}\".format(vo_truth.shape))\n",
    "print(\"vo_interp :\",vo_interp.shape)\n",
    "print('6*7*5*4',6*7*5*4)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.8"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}