alerts_summary.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy.stats as spstat\n",
    "import mysql.connector as mariadb\n",
    "from astropy.time import Time\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.container { width:90% !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.core.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:90% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.set_option('display.max_columns',100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#uname = os.environ[\"mariadb_username\"]\n",
    "#passwd = os.environ[\"mariadb_password\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = mariadb.connect(user='ztf', database='ztf')\n",
    "cur = con.cursor()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1646515\n",
      "CPU times: user 24.9 s, sys: 2.91 s, total: 27.8 s\n",
      "Wall time: 20.7 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "df = pd.read_sql_query('select objectId, count(objectId) as nobs from alerts GROUP BY objectId HAVING nobs > 1', con)\n",
    "obj_lst = df['objectId'].values.tolist()\n",
    "print(len(df))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def retrieve_lc(objectName):\n",
    "\n",
    "    dflc = pd.read_sql_query(\"SELECT objectId, magpsf, sigmapsf, ra, decl, fid, classtar, rb, \\\n",
    "                             distnr, magnr, sigmagnr, chinr, sharpnr, \\\n",
    "                             sgmag1, srmag1, simag1, sgscore1, \\\n",
    "                             distpsnr1 FROM alerts where objectId='{}' \".format(objectName), con=con)\n",
    "\n",
    "    return dflc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ZTFlc():\n",
    "    '''Light curve object for ZTF formatted data'''\n",
    "\n",
    "    def __init__(self, objId):\n",
    "        '''Read in light curve data'''\n",
    "        DFlc = retrieve_lc(objId)\n",
    "        self.DFlc = DFlc\n",
    "        self.objId = objId\n",
    "        self.filename = objId\n",
    "        self.filter_mag()\n",
    "\n",
    "    def plot_multicolor_lc(self):\n",
    "        '''Plot the 3 band light curve'''\n",
    "        fig, ax = plt.subplots()\n",
    "\n",
    "        now = Time.now().jd\n",
    "        tc = self.DFlc['jd'] - now\n",
    "\n",
    "        g = ax.errorbar(tc.loc[self.DFlc['fid'] == 1],\n",
    "                     self.DFlc['magpsf'].loc[self.DFlc['fid'] == 1],\n",
    "                     self.DFlc['sigmapsf'].loc[self.DFlc['fid'] == 1],\n",
    "                     fmt = 'o', ms = 3, color = '#556B2F', label = \"g\")\n",
    "        r = ax.errorbar(tc.loc[self.DFlc['fid'] == 2],\n",
    "                     self.DFlc['magpsf'].loc[self.DFlc['fid'] == 2],\n",
    "                     self.DFlc['sigmapsf'].loc[self.DFlc['fid'] == 2],\n",
    "                     fmt = 'o', ms = 3, color = '#8B0000', label = \"r\")\n",
    "        i = ax.errorbar(tc.loc[self.DFlc['fid'] == 3],\n",
    "                     self.DFlc['magpsf'].loc[self.DFlc['fid'] == 3],\n",
    "                     self.DFlc['sigmapsf'].loc[self.DFlc['fid'] == 3],\n",
    "                     fmt = 'o', ms =3, color = '#8B008B', label = \"i\")\n",
    "\n",
    "        ax.legend(fancybox = True)\n",
    "        plt.gca().invert_yaxis()\n",
    "        ax.set_xlabel(r\"$\\mathrm{Days Ago}$\")\n",
    "        ax.set_ylabel(r\"$\\mathrm{Magnitude}$\")\n",
    "\n",
    "\n",
    "    def filter_mag(self):\n",
    "        '''Store individual passband mages as object attributes \n",
    "        and # of observations for each band.'''\n",
    "\n",
    "        self.gMag = self.DFlc['magpsf'].loc[self.DFlc['fid'] == 1]\n",
    "        self.gDmag = self.DFlc['sigmapsf'].loc[self.DFlc['fid'] == 1]\n",
    "\n",
    "        self.rMag = self.DFlc['magpsf'].loc[self.DFlc['fid'] == 2]\n",
    "        self.rDmag = self.DFlc['sigmapsf'].loc[self.DFlc['fid'] == 2]\n",
    "\n",
    "        self.iMag = self.DFlc['magpsf'].loc[self.DFlc['fid'] == 3]\n",
    "        self.iDmag = self.DFlc['sigmapsf'].loc[self.DFlc['fid'] == 3]\n",
    "\n",
    "        self.gNobs = len(self.gMag)\n",
    "        self.rNobs = len(self.rMag)\n",
    "        self.iNobs = len(self.iMag)\n",
    "        \n",
    "    def weighted_mean_flux(self):\n",
    "        '''monochromatic flux density using f_nu = 3631 * 10**(-0.4 mag) [Jy]'''\n",
    "        \n",
    "        if not hasattr(self, 'gMag'):\n",
    "            self.filter_mag()\n",
    "            \n",
    "        #flux = lambda mag, dmag: np.sum((3631.*10**(-0.4*mag))*(1./(3631.*10**(-0.4*dmag)))**2)/np.sum((1./(3631.*10**(-0.4*dmag)))**2)\n",
    "        def weighted_mean(mag,dmag):\n",
    "            flux_mag = 3631.*10**(-0.4*mag)\n",
    "            flux_dmag = 3631.*10**(-0.4*dmag)\n",
    "        \n",
    "            return np.sum(flux_mag*(1./flux_dmag)**2)/np.sum((1./flux_dmag)**2)\n",
    "        \n",
    "        self.gMeanJy = weighted_mean(self.gMag, self.gDmag) if self.gNobs > 0 else -999\n",
    "        self.rMeanJy = weighted_mean(self.rMag, self.rDmag) if self.rNobs > 0 else -999\n",
    "        self.iMeanJy = weighted_mean(self.iMag, self.iDmag) if self.iNobs > 0 else -999\n",
    "        \n",
    "    def sample_std_flux(self):\n",
    "        '''Measure standard deviation of sample flux in gri'''\n",
    "\n",
    "        if not hasattr(self, 'gMag'):\n",
    "            self.filter_mag()\n",
    "        \n",
    "        flux = lambda mag: 3631.*10**(-0.4*mag)\n",
    "        \n",
    "        self.gFluxes = flux(self.gMag)\n",
    "        self.rFluxes = flux(self.rMag)\n",
    "        self.iFluxes = flux(self.iMag)\n",
    "        \n",
    "        self.gStdJy = np.std(self.gFluxes, ddof = 1) if self.gNobs > 1 else -999\n",
    "        self.rStdJy = np.std(self.rFluxes, ddof = 1) if self.rNobs > 1 else -999\n",
    "        self.iStdJy = np.std(self.iFluxes, ddof = 1) if self.iNobs > 1 else -999\n",
    "        \n",
    "    def real_bogus(self):\n",
    "        '''Store mean real_bogus score.'''\n",
    "\n",
    "        if not hasattr(self, 'gMag'):\n",
    "            self.filter_mag()\n",
    "\n",
    "        dummy = np.mean(self.gMag)\n",
    "\n",
    "        self.gRB = np.mean(self.DFlc['rb'].loc[self.DFlc['fid'] == 1]) if self.gNobs > 0 else -999\n",
    "        self.rRB = np.mean(self.DFlc['rb'].loc[self.DFlc['fid'] == 2]) if self.rNobs > 0 else -999\n",
    "        self.iRB = np.mean(self.DFlc['rb'].loc[self.DFlc['fid'] == 3]) if self.iNobs > 0 else -999\n",
    "\n",
    "    def star_galaxy(self):\n",
    "        '''Store mean star_galaxy score.'''\n",
    "\n",
    "        if not hasattr(self, 'gMag'):\n",
    "            self.filter_mag()\n",
    "\n",
    "        dummy = np.mean(self.gMag)\n",
    "\n",
    "        self.gClasstar = np.mean(self.DFlc['classtar'].loc[self.DFlc['fid'] == 1]) if self.gNobs > 0 else -999\n",
    "        self.rClasstar = np.mean(self.DFlc['classtar'].loc[self.DFlc['fid'] == 2]) if self.rNobs > 0 else -999\n",
    "        self.iClasstar = np.mean(self.DFlc['classtar'].loc[self.DFlc['fid'] == 3]) if self.iNobs > 0 else -999\n",
    "\n",
    "    def ra_dec(self):\n",
    "        '''Store mean RA & Dec values.'''\n",
    "\n",
    "        self.meanRA = np.mean(self.DFlc['ra'])\n",
    "        self.meanDec = np.mean(self.DFlc['decl'])\n",
    "\n",
    "    def weighted_mean_mag(self):\n",
    "        '''Measure weighted mean mag in gri'''\n",
    "\n",
    "        if not hasattr(self, 'gMag'):\n",
    "            self.filter_mag()\n",
    "\n",
    "        weighted_mean = lambda mag, dmag: np.sum(mag*(1./dmag)**2)/np.sum((1./dmag)**2)\n",
    "\n",
    "        self.gMeanAB = weighted_mean(self.gMag, self.gDmag) if self.gNobs > 0 else -999\n",
    "        self.rMeanAB = weighted_mean(self.rMag, self.rDmag) if self.rNobs > 0 else -999\n",
    "        self.iMeanAB = weighted_mean(self.iMag, self.iDmag) if self.iNobs > 0 else -999\n",
    "\n",
    "    def sample_std_mag(self):\n",
    "        '''Measure standard deviation of sample mag in gri'''\n",
    "\n",
    "        if not hasattr(self, 'gMag'):\n",
    "            self.filter_mag()\n",
    "\n",
    "        self.gStdAB = np.std(self.gMag, ddof = 1) if self.gNobs > 1 else -999\n",
    "        self.rStdAB = np.std(self.rMag, ddof = 1) if self.rNobs > 1 else -999\n",
    "        self.iStdAB = np.std(self.iMag, ddof = 1) if self.iNobs > 1 else -999\n",
    "\n",
    "    def normalized_amplitude(self):\n",
    "        '''Measure the normalized amplitude of variations in gri'''\n",
    "\n",
    "        if not hasattr(self, 'gMag'):\n",
    "            self.filter_mag()\n",
    "\n",
    "        if not hasattr(self, 'gMean'):\n",
    "            self.weighted_mean_mag()\n",
    "\n",
    "        normalized_amplitude = lambda mag, wMeanmag: (np.max(mag) - np.min(mag))/wMeanmag\n",
    "\n",
    "        self.gAmp = normalized_amplitude(self.gMag, self.gMean) if self.gNobs else -999\n",
    "        self.rAmp = normalized_amplitude(self.rMag, self.rMean) if self.rNobs else -999\n",
    "        self.iAmp = normalized_amplitude(self.iMag, self.iMean) if self.iNobs else -999\n",
    "\n",
    "    def normalized_MAD(self):\n",
    "        '''Measure normalized Median Absolute Deviation (MAD) in gri'''\n",
    "\n",
    "        if not hasattr(self, 'gMag'):\n",
    "            self.filter_mag()\n",
    "\n",
    "        if not hasattr(self, 'gMean'):\n",
    "            self.weighted_mean_mag()\n",
    "\n",
    "        normalized_MAD = lambda mag, wMeanmag: np.median(np.abs((mag - np.median(mag))/wMeanmag))\n",
    "\n",
    "        self.gMAD = normalized_MAD(self.gMag, self.gMean) if self.gNobs else -999\n",
    "        self.rMAD = normalized_MAD(self.rMag, self.rMean) if self.rNobs else -999\n",
    "        self.iMAD = normalized_MAD(self.iMag, self.iMean) if self.iNobs else -999\n",
    "        \n",
    "    def normalized_beyond_1std(self):\n",
    "        '''Measure fraction of mag measurements beyond 1 std'''\n",
    "\n",
    "        if not hasattr(self, 'gMag'):\n",
    "            self.filter_mag()\n",
    "\n",
    "        if not hasattr(self, 'gMean'):\n",
    "            self.weighted_mean_mag()\n",
    "\n",
    "        beyond_1std = lambda mag, wMeanmag: sum(np.abs(mag - wMeanmag) > np.std(mag, ddof = 1))/len(mag)\n",
    "\n",
    "        self.gBeyond = beyond_1std(self.gMag, self.gMean) if self.gNobs else -999\n",
    "        self.rBeyond = beyond_1std(self.rMag, self.rMean) if self.rNobs else -999\n",
    "        self.iBeyond = beyond_1std(self.iMag, self.iMean) if self.iNobs else -999\n",
    "\n",
    "    def skew(self):\n",
    "        '''Measure the skew of the mag measurements'''\n",
    "        if not hasattr(self, 'gMag'):\n",
    "            self.filter_mag()\n",
    "\n",
    "        skew = lambda mag: spstat.skew(mag)\n",
    "\n",
    "        self.gSkew = skew(self.gMag) if self.gNobs else -999\n",
    "        self.rSkew = skew(self.rMag) if self.rNobs else -999\n",
    "        self.iSkew = skew(self.iMag) if self.iNobs else -999\n",
    "\n",
    "    def mean_colors(self):\n",
    "        '''Measure the mean g-r, g-i, and r-i colors'''\n",
    "\n",
    "        if not hasattr(self, 'gMag'):\n",
    "            self.filter_mag()\n",
    "\n",
    "        if not hasattr(self, 'gMean'):\n",
    "            self.weighted_mean_mag()\n",
    "\n",
    "        self.gMinusR = (self.gMean - self.rMean) if self.gMean> 0 and self.rMean > 0 else -999\n",
    "        self.rMinusI = (self.rMean - self.iMean) if self.rMean> 0 and self.iMean > 0 else -999\n",
    "        self.gMinusI = (self.gMean - self.iMean) if self.gMean> 0 and self.iMean > 0 else -999\n",
    "        \n",
    "    def nearest_source(self):\n",
    "        '''Nearest source Info in reference image PSF-catalog within 30 arcsec'''\n",
    "        \n",
    "        if not hasattr(self, 'gMag'):\n",
    "            self.filter_mag()\n",
    "\n",
    "        dummy = np.mean(self.gMag)\n",
    "        \n",
    "        self.Distnr = np.mean(self.DFlc['distnr'])\n",
    "        \n",
    "        self.gMagnr = np.mean(self.DFlc['magnr'].loc[self.DFlc['fid'] == 1]) if self.gNobs > 0 else -999\n",
    "        self.gDmagnr = np.mean(self.DFlc['sigmagnr'].loc[self.DFlc['fid'] == 1]) if self.gNobs > 0 else -999\n",
    "\n",
    "        self.rMagnr = np.mean(self.DFlc['magnr'].loc[self.DFlc['fid'] == 2]) if self.rNobs > 0 else -999\n",
    "        self.rDmagnr = np.mean(self.DFlc['sigmagnr'].loc[self.DFlc['fid'] == 2]) if self.rNobs > 0 else -999\n",
    "\n",
    "        self.iMagnr = np.mean(self.DFlc['magnr'].loc[self.DFlc['fid'] == 3]) if self.iNobs > 0 else -999\n",
    "        self.iDmagnr = np.mean(self.DFlc['sigmagnr'].loc[self.DFlc['fid'] == 3]) if self.iNobs > 0 else -999\n",
    "        \n",
    "        self.gChinr = np.mean(self.DFlc['chinr'].loc[self.DFlc['fid'] == 1]) if self.gNobs > 0 else -999\n",
    "        self.rChinr = np.mean(self.DFlc['chinr'].loc[self.DFlc['fid'] == 2]) if self.rNobs > 0 else -999\n",
    "        self.iChinr = np.mean(self.DFlc['chinr'].loc[self.DFlc['fid'] == 3]) if self.iNobs > 0 else -999\n",
    "        \n",
    "        self.gSharpnr = np.mean(self.DFlc['sharpnr'].loc[self.DFlc['fid'] == 1]) if self.gNobs > 0 else -999\n",
    "        self.rSharpnr = np.mean(self.DFlc['sharpnr'].loc[self.DFlc['fid'] == 2]) if self.rNobs > 0 else -999\n",
    "        self.iSharpnr = np.mean(self.DFlc['sharpnr'].loc[self.DFlc['fid'] == 3]) if self.iNobs > 0 else -999\n",
    "\n",
    "    def PS1_closest_source(self):\n",
    "        '''closest source from PS1 catalog; if exists within 30 arcsec'''       \n",
    "        \n",
    "        self.Distnrps1 = np.mean(self.DFlc['distpsnr1'])\n",
    "        \n",
    "        self.gMagps1 = self.DFlc['sgmag1'][0]\n",
    "        self.rMagps1 = self.DFlc['srmag1'][0]\n",
    "        self.iMagps1 = self.DFlc['simag1'][0]\n",
    "        \n",
    "        self.Sgscoreps1 = self.DFlc['sgscore1'][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_features(obj):\n",
    "    lc = ZTFlc(obj)\n",
    "    objectId = lc.objId\n",
    "    lc.ra_dec()\n",
    "    lc.filter_mag()\n",
    "    lc.real_bogus()\n",
    "    lc.star_galaxy()\n",
    "    lc.weighted_mean_mag()\n",
    "    lc.sample_std_mag()\n",
    "    lc.weighted_mean_flux()\n",
    "    lc.sample_std_flux()\n",
    "    lc.nearest_source()\n",
    "    lc.PS1_closest_source()\n",
    "    \n",
    "    \n",
    "    feats = (objectId, lc.meanRA, lc.meanDec,\n",
    "             lc.gNobs,lc.rNobs,lc.iNobs,\n",
    "             lc.gMeanAB,lc.rMeanAB,lc.iMeanAB,\n",
    "             lc.gStdAB,lc.rStdAB,lc.iStdAB,\n",
    "             lc.gMeanJy, lc.rMeanJy, lc.iMeanJy, lc.gStdJy, lc.rStdJy, lc.iStdJy,\n",
    "             lc.gRB,lc.rRB,lc.iRB,\n",
    "             lc.gClasstar,lc.rClasstar,lc.iClasstar,\n",
    "             lc.Distnr, lc.gMagnr, lc.gDmagnr, lc.rMagnr, lc.rDmagnr, lc.iMagnr, lc.iDmagnr,\n",
    "             lc.gChinr, lc.rChinr, lc.iChinr, lc.gSharpnr, lc.rSharpnr, lc.iSharpnr,\n",
    "             lc.Distnrps1, lc.gMagps1, lc.rMagps1, lc.iMagps1, lc.Sgscoreps1)\n",
    "    \n",
    "    return feats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'\\nfor obj in obj_lst:\\n    \\n    feats = calculate_features(obj)\\n\\n    try:\\n        cur.execute(\\'INSERT INTO summary2(objectId, meanRA, meanDec,                                           gNobs, rNobs, iNobs,                                           gMeanAB,rMeanAB,iMeanAB, gStdAB, rStdAB, iStdAB,                                           gMeanJy, rMeanJy, iMeanJy, gStdJy, rStdJy, iStdJy,                                           gRB, rRB, iRB, gClasstar, rClasstar, iClasstar,                                           Distnr, gMagnr, gDmagnr, rMagnr, rDmagnr, iMagnr, iDmagnr,                                           gChinr, rChinr, iChinr, gSharpnr, rSharpnr, iSharpnr,                                           Distnrps1, gMagps1, rMagps1, iMagps1, Sgscoreps1) values {}\\'.format(feats))\\n\\n    except mariadb.Error as error:\\n        print(\"Error: {}\".format(error))\\n\\ncon.commit()\\ncon.close()\\n'"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "for obj in obj_lst:\n",
    "    \n",
    "    feats = calculate_features(obj)\n",
    "\n",
    "    try:\n",
    "        cur.execute('INSERT INTO summary2(objectId, meanRA, meanDec, \\\n",
    "                                          gNobs, rNobs, iNobs, \\\n",
    "                                          gMeanAB,rMeanAB,iMeanAB, gStdAB, rStdAB, iStdAB, \\\n",
    "                                          gMeanJy, rMeanJy, iMeanJy, gStdJy, rStdJy, iStdJy, \\\n",
    "                                          gRB, rRB, iRB, gClasstar, rClasstar, iClasstar, \\\n",
    "                                          Distnr, gMagnr, gDmagnr, rMagnr, rDmagnr, iMagnr, iDmagnr, \\\n",
    "                                          gChinr, rChinr, iChinr, gSharpnr, rSharpnr, iSharpnr, \\\n",
    "                                          Distnrps1, gMagps1, rMagps1, iMagps1, Sgscoreps1) values {}'.format(feats))\n",
    "\n",
    "    except mariadb.Error as error:\n",
    "        print(\"Error: {}\".format(error))\n",
    "\n",
    "con.commit()\n",
    "con.close()\n",
    "\"\"\""
   ]
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}