{"cells":[{"metadata":{"nbpresent":{"id":"aaffd2c3-54c1-4f88-840b-02111a6237a8"},"slideshow":{"slide_type":"slide"},"_uuid":"034e73ae1a13904dd04a7a2b099695023f8f27c0"},"cell_type":"markdown","source":"## Plotting light curves of training data set\nInspired by <a href=\"https://www.kaggle.com/michaelapers/the-plasticc-astronomy-starter-kit\">The PLAsTiCC Astronomy \"Starter Kit\"</a> kernel, I tried to plot light curves of training data set. \nI wonder if this kernel is helpful to classify objects, but I expect to get light curve overview of each target ID.\n\nThank you for creating \"Starter Kit\" kernel !"},{"metadata":{"nbpresent":{"id":"c06c5e08-417f-4806-84ca-55c1fc26b576"},"slideshow":{"slide_type":"skip"},"trusted":true,"_uuid":"c9ad93c9583b438947705ff7dbf89006c6eb5eeb"},"cell_type":"code","source":"# You can edit the font size here to make rendered text more comfortable to read\n# It was built on a 13\" retina screen with 18px\nfrom IPython.core.display import display, HTML\ndisplay(HTML(\"<style>.rendered_html { font-size: 18px; }</style>\"))\n\n# we'll also use this package to read tables\n# it's generally useful for astrophysics work, including this challenge\n# so we'd suggest installing it, even if you elect to work with pandas\nfrom astropy.table import Table","execution_count":null,"outputs":[]},{"metadata":{"nbpresent":{"id":"1fe8a6e9-5eba-4ab1-a2ca-7888f3275dde"},"slideshow":{"slide_type":"skip"},"_uuid":"297ce64b583e8769d2fae4abe2c2fb2a5bcad57a"},"cell_type":"markdown","source":"I arranged LightCurve class in the \"Starter Kit\" kernel in order to:\n* use Lomb-Scargle Periodogram\n* display light curve plots horizontally:  x-axis are 'MJD' and 'Phase'\n* display target ID on light curve plots"},{"metadata":{"nbpresent":{"id":"9d5ae1ff-bcc0-42a8-b366-38703a122397"},"slideshow":{"slide_type":"skip"},"trusted":true,"_uuid":"e842ed69bc5e82ff57d5ba10efea457efcb4a01a"},"cell_type":"code","source":"import os\nimport numpy as np\nimport scipy.stats as spstat\nimport matplotlib.pyplot as plt\nfrom collections import OrderedDict\n\n%matplotlib inline\n\nfrom gatspy.periodic import LombScargleMultiband","execution_count":null,"outputs":[]},{"metadata":{"nbpresent":{"id":"233e0dfe-c579-42b2-ab56-3b1e83492c6b"},"slideshow":{"slide_type":"skip"},"trusted":true,"_uuid":"29ba2f55916b27d0641185d60893005b87a516ec"},"cell_type":"code","source":"class LightCurve(object):\n    '''Light curve object for PLAsTiCC formatted data'''\n    \n    _passbands = OrderedDict([(0,'C4'),\\\n                              (1,'C2'),\\\n                              (2,'C3'),\\\n                              (3,'C1'),\\\n                              (4,'k'),\\\n                              (5,'C5')])\n    \n    _pbnames = ['u','g','r','i','z','y']\n    \n    #def __init__(self, filename):\n    def __init__(self, fluxDF):\n        '''Read in light curve data'''\n\n        #self.DFlc     = Table.read(filename, format='ascii.csv')\n        self.DFlc     = fluxDF\n        #self.filename = filename.replace('.csv','')\n        self._finalize()\n     \n    # this is some simple code to demonstrate how to calculate features on these multiband light curves\n    # we're not suggesting using these features specifically\n    # there also might be additional pre-processing you do before computing anything\n    # it's purely for illustration\n    def _finalize(self):\n        '''Store individual passband fluxes as object attributes'''\n        # in this example, we'll use the weighted mean to normalize the features\n        weighted_mean = lambda flux, dflux: np.sum(flux*(flux/dflux)**2)/np.sum((flux/dflux)**2)\n        \n        # define some functions to compute simple descriptive statistics\n        normalized_flux_std = lambda flux, wMeanFlux: np.std(flux/wMeanFlux, ddof = 1)\n        normalized_amplitude = lambda flux, wMeanFlux: (np.max(flux) - np.min(flux))/wMeanFlux\n        normalized_MAD = lambda flux, wMeanFlux: np.median(np.abs((flux - np.median(flux))/wMeanFlux))\n        beyond_1std = lambda flux, wMeanFlux: sum(np.abs(flux - wMeanFlux) > np.std(flux, ddof = 1))/len(flux)\n        \n        for pb in self._passbands:\n            ind = self.DFlc['passband'] == pb\n            pbname = self._pbnames[pb]\n            \n            if len(self.DFlc[ind]) == 0:\n                setattr(self, f'{pbname}Std', np.nan)\n                setattr(self, f'{pbname}Amp', np.nan)\n                setattr(self, f'{pbname}MAD', np.nan)\n                setattr(self, f'{pbname}Beyond', np.nan)\n                setattr(self, f'{pbname}Skew', np.nan)\n                continue\n            \n            f  = self.DFlc['flux'][ind]\n            df = self.DFlc['flux_err'][ind]\n            m  = weighted_mean(f, df)\n            \n            # we'll save the measurements in each passband to simplify access.\n            setattr(self, f'{pbname}Flux', f)\n            setattr(self, f'{pbname}FluxUnc', df)\n            setattr(self, f'{pbname}Mean', m)\n            \n            # compute the features\n            std = normalized_flux_std(f, df)\n            amp = normalized_amplitude(f, m)\n            mad = normalized_MAD(f, m)\n            beyond = beyond_1std(f, m)\n            skew = spstat.skew(f) \n            \n            # and save the features\n            setattr(self, f'{pbname}Std', std)\n            setattr(self, f'{pbname}Amp', amp)\n            setattr(self, f'{pbname}MAD', mad)\n            setattr(self, f'{pbname}Beyond', beyond)\n            setattr(self, f'{pbname}Skew', skew)\n        \n        # we can also construct features between passbands\n        pbs = list(self._passbands.keys())\n        for i, lpb in enumerate(pbs[0:-1]):\n            rpb = pbs[i+1]\n            \n            lpbname = self._pbnames[lpb]\n            rpbname = self._pbnames[rpb]\n            \n            colname = '{}Minus{}'.format(lpbname, rpbname.upper())\n            lMean = getattr(self, f'{lpbname}Mean', np.nan)\n            rMean = getattr(self, f'{rpbname}Mean', np.nan)\n            col = -2.5*np.log10(lMean/rMean) if lMean> 0 and rMean > 0 else -999\n            setattr(self, colname, col)\n    \n    def plot_multicolor_lc(self, target_id=None):\n        '''Plot the multiband light curve'''\n        \n        # Lomb-Scargle\n        model = LombScargleMultiband(fit_period=True)\n        # we'll window the search range by setting minimums and maximums here\n        # but in general, the search range you want to evaluate will depend on the data\n        # and you will not be able to window like this unless you know something about\n        # the class of the object a priori\n        t_min = max(np.median(np.diff(sorted(self.DFlc['mjd']))), 0.1)\n        t_max = min(10., (self.DFlc['mjd'].max() - self.DFlc['mjd'].min())/2.)\n        \n        model.optimizer.set(period_range=(t_min, t_max), first_pass_coverage=5)\n        model.fit(self.DFlc['mjd'], self.DFlc['flux'], dy=self.DFlc['flux_err'], filts=self.DFlc['passband'])\n        period = model.best_period\n        obj_id = self.DFlc['object_id'][0] # object ID\n        print(f'object ID: {obj_id} has a period of {period} days')\n        \n        phase = (self.DFlc['mjd'] /period) % 1\n        \n        #fig, ax = plt.subplots(figsize=(8,6))\n        fig, (ax, ax2) = plt.subplots(ncols=2, figsize=(16,6))\n\n        #if phase is None:\n        #    phase = []\n        #if len(phase) != len(self.DFlc):\n        #    phase = self.DFlc['mjd']\n        #    xlabel = 'MJD'\n        #else:\n        #    xlabel = 'Phase'\n          \n        for i, pb in enumerate(self._passbands):\n            pbname = self._pbnames[pb]\n            ind = self.DFlc['passband'] == pb\n            if len(self.DFlc[ind]) == 0:\n                continue\n            # errorbar: plot y versus x as lines and/or markers with attached errorbars\n            #ax.errorbar(phase[ind], \n            #         self.DFlc['flux'][ind],\n            #         self.DFlc['flux_err'][ind],\n            #         fmt = 'o', color = self._passbands[pb], label = f'{pbname}')\n            ax.errorbar(self.DFlc['mjd'][ind], \n                     self.DFlc['flux'][ind],\n                     self.DFlc['flux_err'][ind],\n                     fmt = 'o', color = self._passbands[pb], label = f'{pbname}')\n            ax2.errorbar(phase[ind], \n                     self.DFlc['flux'][ind],\n                     self.DFlc['flux_err'][ind],\n                     fmt = 'o', color = self._passbands[pb], label = f'{pbname}')\n        ax.legend(ncol = 4, frameon = True)\n        #ax.set_xlabel(f'{xlabel}', fontsize='large')\n        ax.set_xlabel('MJD', fontsize='large')\n        ax.set_ylabel('Flux', fontsize='large')\n        ax2.legend(ncol = 4, frameon = True)\n        ax2.set_xlabel('Phase', fontsize='large')\n        ax2.set_ylabel('Flux', fontsize='large')\n        #fig.suptitle(self.filename, fontsize='x-large')\n        fig.suptitle('object ID: ' + str(self.DFlc['object_id'][0]) + ', target ID: ' + str(target_id), fontsize='x-large') # graph title = object ID\n        fig.tight_layout(rect=[0, 0, 1, 0.97])\n    \n    def get_features(self):\n        '''Return all the features for this object'''\n        variables = ['Std', 'Amp', 'MAD', 'Beyond', 'Skew']\n        feats = []\n        for i, pb in enumerate(self._passbands):\n            pbname = self._pbnames[pb]\n            feats += [getattr(self, f'{pbname}{x}', np.nan) for x in variables]\n        return feats","execution_count":null,"outputs":[]},{"metadata":{"nbpresent":{"id":"3ff5aa2b-f582-44ab-9552-54c6f9c512ba"},"slideshow":{"slide_type":"slide"},"_uuid":"c6a6543c2adf80f64cf4f3fc71ab2f2f5726ac85"},"cell_type":"markdown","source":"So let's read the training data set:"},{"metadata":{"trusted":true,"_uuid":"c0872b82f295622d98991c99d1a7ae9d65b09735","scrolled":true},"cell_type":"code","source":"trainfilename = '../input/PLAsTiCC-2018/training_set.csv'\ntrain = Table.read(trainfilename, format='csv')\ntrain","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5fedf3487bdbbc2bb27ead43ded28c81248c6686"},"cell_type":"markdown","source":"I make a sample plot (without displaying the target ID) :"},{"metadata":{"trusted":true,"_uuid":"c7f40ccf62cb3e5dba2d9ec00a92bad26cc3f509"},"cell_type":"code","source":"# read a sample object data\nobj_id = train['object_id'][0]\nfluxDF = train[train['object_id'] == obj_id]\nfluxDF","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"78329a9950740fd75c1c7f5b021b87f7514cdc4d"},"cell_type":"code","source":"lc = LightCurve(fluxDF)\nlc.plot_multicolor_lc()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c7378e84be0237fc6f99477f8d0c93bb4b3817c0"},"cell_type":"markdown","source":"Then, I plot 10 light curves for each target ID. There are many objects in each target ID, so I limited the number of plots."},{"metadata":{"trusted":true,"_uuid":"382ac66c513a1969314f341b7c505f093c61105a","scrolled":false},"cell_type":"code","source":"trainmetafilename = '../input/PLAsTiCC-2018/training_set_metadata.csv'\ntrain_meta = Table.read(trainmetafilename, format='csv')\n# set: remove duplicated object_id\n# list: change type to list\ntarget_list = list(set(train_meta['target']))\nall_obj_id_list = []\n\n# create a list: object IDs which has same target ID\nfor target_id in target_list:\n    #print(train_meta[train_meta['target'] == target_id]['object_id'])\n    all_obj_id_list.append(train_meta[train_meta['target'] == target_id]['object_id'])\n    \nfor i in range(len(all_obj_id_list)):\n#for i in range(2):\n    print('********* target ID: ' + str(target_list[i]) + ' *********')\n    each_obj_id_list = all_obj_id_list[i]\n    print('********* ' + str(len(each_obj_id_list)) + ' objects in target type ' + str(target_list[i]) + ' *********')\n    # WARNING: this takes very long time...\n    #for j in range(len(each_obj_id_list)):\n    for j in range(10):\n        obj_id = each_obj_id_list[j]\n        print ('****** object ID: ' + str(obj_id) + ' ******')\n        fluxDF = train[train['object_id'] == obj_id]\n        lc = LightCurve(fluxDF)\n        lc.plot_multicolor_lc(target_list[i])\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1a6c7dc771e16b546da2d6c4c08bea52536acd7c"},"cell_type":"markdown","source":"Thank you for reading !"}],"metadata":{"anaconda-cloud":{},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"livereveal":{"scroll":true},"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":true,"toc_position":{},"toc_section_display":true,"toc_window_display":true}},"nbformat":4,"nbformat_minor":1}