saro_test.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Testing SARO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import saro\n",
    "import xara\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from time import time\n",
    "from scipy.linalg import sqrtm\n",
    "import scipy.sparse as sparse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating the KPO structure from a pupil text file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thekpo = saro.KPO(\"test_data/pupil_model.txt\")\n",
    "a = thekpo.kpi.plot_pupil_and_uv(cmap=\"plasma\", usesize=True)\n",
    "plt.savefig(\"subaru_pupil_model.pdf\")\n",
    "plt.show(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining basic parameters for the analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wl = 1.6e-6\n",
    "m2pix = xara.mas2rad( 16.7 * 128 / wl)#pscale supposed 16.7\n",
    "thekpo.CWAVEL = wl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importing the test case data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calibrator = np.load(\"test_data/calibrator.npy\")\n",
    "tarimage_single = np.load(\"test_data/tarimage_single.npy\")\n",
    "tarimage_series = np.load(\"test_data/tarimage_series.npy\")\n",
    "params_single = np.load(\"test_data/params_single.npy\")\n",
    "params_series = np.load(\"test_data/params_series.npy\")\n",
    "detpas = np.load(\"test_data/detpas.npy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extracting the visibilities and kernel-phases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cviscal = thekpo.extract_cvis_from_img(calibrator, m2pix)\n",
    "kappacal = thekpo.kpi.KPM.dot(np.angle(cviscal))\n",
    "cvistar_single = thekpo.extract_cvis_from_img(tarimage_single, m2pix)\n",
    "kappasingle = thekpo.kpi.KPM.dot(np.angle(cvistar_single))\n",
    "cvistar_series = np.array([thekpo.extract_cvis_from_img(tarimage_series[i], m2pix) for i in range(detpas.shape[0])])\n",
    "kappa_series = thekpo.kpi.KPM.dot(np.angle(cvistar_series).T).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importing the covariance matrix\n",
    "There are several ways to build a covariance matrix for the data. None of them are perfect.\n",
    "\n",
    "Here the covariance was evaluated from a cube of 10k realizations of independent phase screen and shot noise on a calibrator. It should be very close to ideal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Sigma = np.load(\"test_data/covk_10000_250_50_2_apod_False.npy\")\n",
    "Sigmacal = Sigma.copy()\n",
    "assert Sigma.shape[0] == kappacal.shape[0]\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.imshow(Sigma, cmap=\"inferno\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this matrix is very much NOT diagonal, which shows that the kernel-phases, just like any Fourier plane observables, are heavily correlated. This is the reason we need to use the whitening transforms.\n",
    "\n",
    "Here, we need a whitening matrix for the calibrated data, so the matrix has to be adjusted!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "W = sqrtm(np.linalg.inv(Sigma + Sigmacal))\n",
    "Sigmas = np.array([Sigma for i in range(detpas.shape[0])])\n",
    "Ws = sqrtm(np.linalg.inv(sparse.block_diag(Sigmas).toarray()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Colinearity maps: Single observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#resol = 80\n",
    "#step = 10.\n",
    "resol, gstep = thekpo.define_roi(verbose=True)\n",
    "#We just want prettier maps than the minimum\n",
    "step = 0.5*gstep\n",
    "resol = 2 *resol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thekpo.Mp = W\n",
    "start = time()\n",
    "colinmap, norm, rhos, thetas = thekpo.cpu_sadk_binary_match_map(resol,step,np.array([kappasingle-kappacal]),\n",
    "                                           W=None, deltas=np.array([0.]), full_output=True)\n",
    "single_cmap_time = time()-start\n",
    "print(\"Computation of the map took: %.1f seconds\" % (single_cmap_time))\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.imshow(colinmap,cmap=saro.vdg)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#totalrot.append(dtheta)\n",
    "totalrot = detpas[-1] - detpas[0]\n",
    "thekpo.WlLp = W\n",
    "pgrid = saro.dmap.build_grid(resol, step)\n",
    "start = time()\n",
    "detmap_single = thekpo.get_sadk_contrast_cpu(pgrid, 0.01,0.95,1.6e-6, W=None, dthetas=np.array([0.])).reshape(resol,resol)\n",
    "single_detmap_time = time() - start\n",
    "print(\"Computation of the map took: %.1f seconds\" % (single_detmap_time))\n",
    "logcontrast_single = np.log10(detmap_single)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xx, yy = np.meshgrid(np.arange(-resol//2,resol//2), (np.arange(-resol//2,resol//2)))\n",
    "rad = np.sqrt(xx**2 + yy**2)\n",
    "\n",
    "lengths = np.sqrt(thekpo.kpi.UVC[:,0]**2 + thekpo.kpi.UVC[:,1]**2)\n",
    "D = np.max(lengths)\n",
    "b = np.min(lengths)\n",
    "thekpo.rhomin = 0.5*thekpo.CWAVEL / D * 180/np.pi * 3600 * 1000\n",
    "thekpo.rhomax = 0.5*thekpo.CWAVEL / b * 180/np.pi * 3600 * 1000\n",
    "print(\"rho min\",thekpo.rhomin)\n",
    "print(\"rho max\", thekpo.rhomax)\n",
    "detmask = 1*((rad*step)<thekpo.rhomax)\n",
    "\n",
    "a = plt.figure(figsize=(8,6))\n",
    "plt.imshow(logcontrast_single*detmask, cmap=\"inferno\")\n",
    "plt.colorbar()\n",
    "CSScal = plt.contour((logcontrast_single*detmask),levels=[1.5,2,2.5,3], cmap=\"Greens\")\n",
    "radi = plt.contour(rad*step,levels=[41],linestyles=[\"dashed\"], colors=[\"white\"])\n",
    "plt.clabel(CSScal,  inline=1, fmt=\"%.1f\")\n",
    "plt.title(r\"Contrast energy detector (95% / 1%)  in log ($dashed=1\\lambda/D$)\")\n",
    "plt.tight_layout()\n",
    "#plt.savefig(\"_ADK_detection_map_mag.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Colinearity maps: Series of observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thekpo.Mp = Ws\n",
    "start = time()\n",
    "colinmap, norm, rhos, thetas = thekpo.cpu_sadk_binary_match_map(resol,step,kappa_series-kappacal[None,:],\n",
    "                                           W=None, deltas=detpas, full_output=True)\n",
    "series_cmap_time = time()-start\n",
    "print(\"Computation of the map took: %.1f seconds\" % (series_cmap_time))\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.imshow(colinmap,cmap=saro.vdg)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#totalrot.append(dtheta)\n",
    "totalrot = detpas[-1] - detpas[0]\n",
    "thekpo.WlLp = Ws\n",
    "pgrid = saro.dmap.build_grid(resol, step)\n",
    "start = time()\n",
    "detmap_series = thekpo.get_sadk_contrast_cpu(pgrid, 0.01,0.95,1.6e-6, W=None, dthetas=detpas).reshape(resol,resol)\n",
    "series_detmap_time = time() - start\n",
    "print(\"Computation of the map took: %.1f seconds\" % (series_detmap_time))\n",
    "logcontrast_series = np.log10(detmap_series)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xx, yy = np.meshgrid(np.arange(-resol//2,resol//2), (np.arange(-resol//2,resol//2)))\n",
    "rad = np.sqrt(xx**2 + yy**2)\n",
    "\n",
    "lengths = np.sqrt(thekpo.kpi.UVC[:,0]**2 + thekpo.kpi.UVC[:,1]**2)\n",
    "D = np.max(lengths)\n",
    "b = np.min(lengths)\n",
    "thekpo.rhomin = 0.5*thekpo.CWAVEL / D * 180/np.pi * 3600 * 1000\n",
    "thekpo.rhomax = 0.5*thekpo.CWAVEL / b * 180/np.pi * 3600 * 1000\n",
    "print(\"rho min\",thekpo.rhomin)\n",
    "print(\"rho max\", thekpo.rhomax)\n",
    "detmask = 1*((rad*step)<thekpo.rhomax)\n",
    "\n",
    "a = plt.figure(figsize=(8,6))\n",
    "plt.imshow(logcontrast_series*detmask, cmap=\"inferno\")\n",
    "plt.colorbar()\n",
    "CSScal = plt.contour((logcontrast_series*detmask),levels=[1.5,2,2.5,3], cmap=\"Greens\")\n",
    "radi = plt.contour(rad*step,levels=[41],linestyles=[\"dashed\"], colors=[\"white\"])\n",
    "plt.clabel(CSScal,  inline=1, fmt=\"%.1f\")\n",
    "plt.title(r\"Contrast energy detector (95% / 1%)  in log ($dashed=1\\lambda/D$)\")\n",
    "plt.tight_layout()\n",
    "#plt.savefig(\"_ADK_detection_map_mag.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L1p = thekpo.build_L_matrix(detpas.shape[0], crop=True, giveUf=False).toarray()\n",
    "Lcov = L1p.dot(sparse.block_diag(Sigmas).toarray().dot(L1p.T))\n",
    "Wl = sqrtm(np.linalg.inv(Lcov))\n",
    "WlLp = Wl.dot(L1p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thekpo.Mp = WlLp\n",
    "start = time()\n",
    "colinmap, norm, rhos, thetas = thekpo.cpu_sadk_binary_match_map(resol,step,kappa_series,\n",
    "                                           W=None, deltas=detpas, full_output=True)\n",
    "sadk_cmap_time = time()-start\n",
    "print(\"Computation of the map took: %.1f seconds\" % (sadk_cmap_time))\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.imshow(colinmap,cmap=saro.vdg)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sensitivity maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#totalrot.append(dtheta)\n",
    "totalrot = detpas[-1] - detpas[0]\n",
    "thekpo.WlLp = WlLp\n",
    "pgrid = saro.dmap.build_grid(resol, step)\n",
    "start = time()\n",
    "detmap_adk = thekpo.get_sadk_contrast_cpu(pgrid, 0.01,0.95,1.6e-6, W=None, dthetas=detpas).reshape(resol,resol)\n",
    "sadk_detmap_time = time() - start\n",
    "print(\"Computation of the map took: %.1f seconds\" % (sadk_detmap_time))\n",
    "logcontrast_adk = np.log10(detmap_adk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xx, yy = np.meshgrid(np.arange(-resol//2,resol//2), (np.arange(-resol//2,resol//2)))\n",
    "rad = np.sqrt(xx**2 + yy**2)\n",
    "\n",
    "lengths = np.sqrt(thekpo.kpi.UVC[:,0]**2 + thekpo.kpi.UVC[:,1]**2)\n",
    "D = np.max(lengths)\n",
    "b = np.min(lengths)\n",
    "thekpo.rhomin = 0.5*thekpo.CWAVEL / D * 180/np.pi * 3600 * 1000\n",
    "thekpo.rhomax = 0.5*thekpo.CWAVEL / b * 180/np.pi * 3600 * 1000\n",
    "print(\"rho min\",thekpo.rhomin)\n",
    "print(\"rho max\", thekpo.rhomax)\n",
    "detmask = 1*((rad*step)<thekpo.rhomax)\n",
    "\n",
    "a = plt.figure(figsize=(8,6))\n",
    "plt.imshow(logcontrast_adk*detmask, cmap=\"inferno\")\n",
    "plt.colorbar()\n",
    "CSScal = plt.contour((logcontrast_adk*detmask),levels=[1.5,2,2.5,3], cmap=\"Greens\")\n",
    "radi = plt.contour(rad*step,levels=[41],linestyles=[\"dashed\"], colors=[\"white\"])\n",
    "plt.clabel(CSScal,  inline=1, fmt=\"%.1f\")\n",
    "plt.title(r\"Contrast energy detector (95% / 1%)  in log ($dashed=1\\lambda/D$)\")\n",
    "plt.tight_layout()\n",
    "#plt.savefig(\"_ADK_detection_map_mag.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L1p = thekpo.build_L_matrix(detpas.shape[0], crop=True, giveUf=False).toarray()\n",
    "Lcov = L1p.dot(sparse.block_diag(Sigmas).toarray().dot(L1p.T))\n",
    "Wl = sqrtm(np.linalg.inv(Lcov))\n",
    "WlLp = Wl.dot(L1p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thekpo.Mp = WlLp\n",
    "sglry, sglrmc, phat = thekpo.global_GLR(np.array([kappa_series-kappacal[None,:]]), dtheta=detpas,n=10, N=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thekpo.Mp = W\n",
    "sglry, sglrmc, phat = thekpo.global_GLR(np.array([kappasingle-kappacal]), dtheta=np.array([0.]),n=10, N=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "\n",
    "KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKRMMMRKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK\n",
    "\n",
    "KKKKKKKKKKKKKKKKKKKKKKKKKKM\"`    ,    YKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK\n",
    "\n",
    "KKKKKKKKKKKKKKKKKKKKKKKM\"  ygKKKKKKKN  ]KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK\n",
    "\n",
    "KKKKKKKKKKKKKKKKKKDKK\" ,gKKKKKKKKKKKK  ]KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK\n",
    "\n",
    "KKKKKKKKKKRMM**\"**fF  *MRKKKKKKKKKKKM  #KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK\n",
    "\n",
    "KKKKKM\"     ;wgm#^,#KKKKNBKWm&BKKKKM  #KK5KKKKK*KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK\n",
    "\n",
    "KKM`   y#KKKKKKM,KKKKKKKKKKwKKKKKKM  #KKN KKM,#KKM5555KB555KKM9KKKMKM]555KRKKKKK\n",
    "\n",
    "KM  ,#KKKKDKKKQ#KKKKKKKKKKKKKpYKK^ zKKKKN *aKKKKKN#KKKKJKKKMKN#JKKMBN]KKKKUKKKKK\n",
    "\n",
    "Kp  1KKKKKKN9DKKKKKKKKKKKKKKKKN  y#KKKKKN NpYKKKKN#KKKKJKNC0KN]KNXMBN]KKKKUKKKKK\n",
    "\n",
    "KN   `RKKKKKNBKKKKKKKKKKKKKKKM zm\"KKKKKKN KKKwJ9KN2DDDK2KKKN8N#KKKwKN]DDDKwDDD]K\n",
    "\n",
    "KKKNp    \"*RKKNAKKKKKKKKKKKM,#KMDm\"KKKKKK#KKKKKN#KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK\n",
    "\n",
    "KKKKKKKWw,,      \"\"l\"\"\"*\"Jwg#KKKKKW1KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK\n",
    "\n",
    "KKKKKKKKKKKKKKKKKKN/*DWKKKKKKKKKKKK BKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK\n",
    "\n",
    "KKKKKKKKKKKKKKKKKKKKKW[TKKKKKKKKKKM KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK\n",
    "\n",
    "KKKKKKKKKKKKKKKKKKKKKKKKNwp\"fMMMM\";#KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK\n",
    "\n",
    "RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRKKRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GPU processing\n",
    "See "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}