{"cells":[{"metadata":{"_uuid":"96730c42518755cc720631c8f3f7f58dfcfd225b"},"cell_type":"markdown","source":"# Class 6 metrics\n\nThe class 6 corresponds to microlensing events (see https://en.wikipedia.org/wiki/Gravitational_microlensing for an explanation)\n\nWhat is presented here is a method to fit the data with the theoretical model.\nOutput of the fit can then be employed as features for classification.\n\n### The fit parameters are :\n- The baseline of the model (one per band), corresponding to the flux before amplification\n- The amplitude at maximum (one per band)\n- The time of maximum\n- The Einstein time, reprensenting the duration of the event\n- The Wings parameter, representing the curtosis of the curve, and linked to the physical configuraiton of the event\n\n### Two other numbers are given by the fit :\n- The chi2\n- NDF, number of degrees of freedom\n\n### How to best use these numbers for classification :\n- For members of the class, chi2/NDF should be peaked around 1, while it should be wildly above one for non-members. This is a goodness of fit metric, meaning that an outlier reaching \nchi2/NDF~1 is undistinguishable from a class member. A few class member will have chi2/NDF>>1: this happens when the fit fails.\n- The fit parameters will exhibit the population statistics. Their distribution can be different between ddf=0 or ddf=1, or between training and test sets. They can help to distinguish class members from natural false-positives (like members of class 65)\n\nThis notebook will present the computation of these parameter, and save them in a local file."},{"metadata":{"_uuid":"528c5ef9451a152616911076b292bc7eb50010c2"},"cell_type":"markdown","source":"## Read data and prepare some stuff"},{"metadata":{"trusted":false,"_uuid":"dc88e7a976e89ac57baa5319219d5c87e78a11e2"},"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom scipy import optimize","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"f86b68c667f3437ec64df25d4d1da112bb7f379f"},"cell_type":"code","source":"train_series = pd.read_csv('../input/training_set.csv')\ntrain_metadata = pd.read_csv('../input/training_set_metadata.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"db7dc2caef3b1c095b90760f5d66cfc2a2c1f4c2"},"cell_type":"code","source":"groups = train_series.groupby(['object_id', 'passband'])\ntimes = groups.apply(\n    lambda block: block['mjd'].values).reset_index().rename(columns={0: 'seq'})\nflux = groups.apply(\n    lambda block: block['flux'].values\n).reset_index().rename(columns={0: 'seq'})\nerr = groups.apply(\n    lambda block: block['flux_err'].values\n).reset_index().rename(columns={0: 'seq'})\ndet = groups.apply(\n    lambda block: block['detected'].astype(bool).values\n).reset_index().rename(columns={0: 'seq'})\ntimes_list = times.groupby('object_id').apply(lambda x: x['seq'].tolist()).tolist()\nflux_list = flux.groupby('object_id').apply(lambda x: x['seq'].tolist()).tolist()\nerr_list = err.groupby('object_id').apply(lambda x: x['seq'].tolist()).tolist()\ndet_list = det.groupby('object_id').apply(lambda x: x['seq'].tolist()).tolist()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"10629e77feeefdeee0963c3de8759ee97e8efad8"},"cell_type":"code","source":"object_list=times.groupby('object_id').apply(lambda x: x['object_id'].unique()[0]).tolist()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"b3170f2b8eade16ba85892578b1dd45cdb5eed8d"},"cell_type":"code","source":"colors = ['purple', 'blue', 'green', 'orange', 'red', 'black']\ndef plot_one_object(obj_id,split=False,frequency=0):\n    \n    if split:\n        fig, axes = plt.subplots(2, 3, sharex=True, sharey=True, figsize=(12, 8))\n    else:\n        ax=plt.gca()\n    for band in range(len(colors)):\n        if split:\n            ax = axes[band // 3, band % 3]\n        sample = train_series[(train_series['object_id'] == obj_id) & (train_series['passband']==band)]\n        if frequency:\n            phase=(sample['mjd'] * frequency ) % 1\n        else:\n            phase=sample['mjd']\n        ax.errorbar(x=phase,y=sample['flux'],yerr=sample['flux_err'],c = colors[band],fmt='o',alpha=0.7)\n        ","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"55ead2078b888030aa33805aa8e6bd4573377621"},"cell_type":"code","source":"gal_target=train_metadata[train_metadata['hostgal_photoz']==0.0]['target'].unique()\nt_type={}\nddf_type={}\nwdf_type={}\nfor t in gal_target:\n    t_type[t] = train_metadata[(train_metadata['target']==t)].index\n    ddf_type[t] = train_metadata[(train_metadata['target']==t) & (train_metadata['ddf']==1)].index\n    wdf_type[t] = train_metadata[(train_metadata['target']==t) & (train_metadata['ddf']==0)].index\nt_colors={6:'blue',16:'red',65:'orange',92:'green',53:'purple'}\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"55e9550338fdbe1decb5c9ad7cd3f2456c4cf4dc"},"cell_type":"markdown","source":"### One example of the class"},{"metadata":{"trusted":false,"_uuid":"ee88d1798b13fa712ae37213815d4d4b78494e07"},"cell_type":"code","source":"iobj=train_metadata[(train_metadata['target']==6) & (train_metadata['ddf']==0)]['object_id'].index[5]\n    #plt.figure()\nprint(iobj)\nplot_one_object(object_list[iobj])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"974c0fb9cc07da438ada6b9ee5d6b28ee06c5842"},"cell_type":"markdown","source":"### Theoretical microlensing curve\n\nThe amplitude depends on the angular separation $u$ as\n\n$A(u) = \\dfrac{u^2+2}{u \\sqrt{u^2+4}}$\n\nwhere $u$ depends on the time $t$ as\n\n$u(t) = \\sqrt{u^2_{min} + \\dfrac{(t-t_0)^2}{t_E^2}}$\n\nHowever, the data provided by PLAsTiCC do not include the base flux, \nso we need to infer 2 additional numbers : the baseline (flux outside the microlensing event)\nand the amplification. \n\nPutting everything together, there are also pseudo-degeneracies, so I isolated them in \nthe $w$ parameter representing the \"wings\" to help the convergence.\n\n#### The final formula for the fit is then :\n    \n$F(t) = A S(\\delta t) + B$\n\nwhere $A$ is the amplitude, $B$ the baseline and  \n\n$S(\\delta t) = \\dfrac{\\dfrac{1 + \\delta t^2 + w^2}{\\sqrt{(1 + \\delta t^2)(1+\\delta t^2 + 2w^2)}}-1}\n{\\dfrac{1 + w^2}{\\sqrt{1+2w^2}}-1}$\n\nwhere $\\delta t^2 = \\dfrac{(t-t_0)^2}{t_E^2}$"},{"metadata":{"trusted":false,"_uuid":"442615a1159133fe59731352f73a16b30a5ecf0e"},"cell_type":"code","source":"# This is the formula from Wikipedia (not used here)\n\ndef microamp(t,params):\n    (t0,te,umin)=params\n    u = np.sqrt( umin**2 + ((t-t0)/te)**2 )\n    amp = (u**2 + 2) / (u * np.sqrt(u**2+4))\n    return amp","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"e832d8b34e968aaa542bf3a32531fcc1d6878da4"},"cell_type":"code","source":"# This is the Shape function\n\ndef microshape(t,params):\n    (t0,te,wings) = params\n    tvec=((t-t0)/te)**2\n    flux0 = 1.0 / ((1+wings**2) / (np.sqrt(1 + 2*wings**2)) -1)\n    return flux0*((1 + tvec + wings**2) / np.sqrt((1+tvec)*(1+tvec+2*wings**2))-1)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6b2d0c424e6a5037faa226aee39a131ad3f8600e"},"cell_type":"markdown","source":"#### A few examples..."},{"metadata":{"trusted":false,"_uuid":"f81cd682f1c3eafba93abb2ac791eef692619845"},"cell_type":"code","source":"plt.plot(np.linspace(-10,10,300),microshape(np.linspace(-10,10,300),(0.0,1.0,np.sqrt(2.)))*4+2)\nplt.plot(np.linspace(-10,10,300),microshape(np.linspace(-10,10,300),(0.0,1.0,np.sqrt(200)))*4+2)\nplt.plot(np.linspace(-10,10,300),microshape(np.linspace(-10,10,300),(5.0,0.6,np.sqrt(200)))*3+2)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7d314ef8ee501320a68a0b3fa9acb6962ae705c4"},"cell_type":"markdown","source":"#### Coding the fit engine..."},{"metadata":{"trusted":false,"_uuid":"2604f759c18da623e6f4e76358a4bdefdafcec60"},"cell_type":"code","source":"def fit_microlens(time,flux,err,paraminit=False):\n    wmat=np.zeros((2,2))\n    yvec=np.zeros(2)\n    baseamp = np.zeros((6,2))\n    \n    def _chi2_microlens(params):\n        (t0,te,wings) = params\n        chi2 = 0\n        for b in range(6):\n            shape = microshape(time[b],(t0,te,wings)) \n            # now, predict value of base and amp.\n            wmat[0][0] = np.sum(shape**2/err[b]**2)\n            wmat[0][1] = wmat[1][0] = np.sum(shape/err[b]**2)\n            wmat[1][1] = np.sum(1./err[b]**2)\n            yvec[0] = np.sum(flux[b]*shape/err[b]**2)\n            yvec[1] = np.sum(flux[b]/err[b]**2)\n            try:\n                baseamp[b] = np.linalg.inv(wmat).dot(yvec)\n                baseamp[b][0] = np.max((baseamp[b][0],0))\n            except:\n                baseamp[b][0]=0\n                baseamp[b][1]=np.sum(flux[b]/err[b]**2)/np.sum(1./err[b])\n            \n            chi2 += np.sum(( baseamp[b][0] * shape + baseamp[b][1] - flux[b] ) **2 / err[b]**2)\n        \n        return chi2\n    \n    if not paraminit:\n        t0 = np.concatenate(time)[np.where(np.concatenate(flux)/np.concatenate(err)==np.max(np.concatenate(flux)/np.concatenate(err)))]\n        te = 100.\n        wings = 1.0\n    else:\n        (t0,te,wings) = paraminit\n    \n    #_chi2_microlens((t0,te,wings))\n    #return (baseamp, (t0,te,wings))\n    \n    params=optimize.fmin(_chi2_microlens,(t0,te,wings),full_output=1,disp=0)\n    return (baseamp,params[0], params[1],np.sum([len(f) for f in flux]))\n    #return optimize.fmin(_chi2_microlens,(base,amp0,t0,te,wings))\n\ndef plot_microlens(time,flux,err,params):\n    fig, axes = plt.subplots(2, 3, sharex=True, sharey=True, figsize=(12, 8))\n    #plt.title(target_types[train_metadata.loc[iobj,'target']]) \n    for band in range(6):\n        ax = axes[band // 3, band % 3]\n        ax.errorbar(time[band],flux[band],yerr=err[band],fmt='o',color=colors[band])\n        time_grid=(np.arange(59550,60705,5.))\n        ymicro=microshape(time_grid,params[1])*params[0][band][0] + params[0][band][1]\n        ax.plot(time_grid,ymicro,color=colors[band])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5917a6841189be5205ea3220cd9741bfde0ab66b"},"cell_type":"markdown","source":"### How to fit the model on a given time serie"},{"metadata":{"trusted":false,"_uuid":"6e7701e203db5042adf9646dc7fe41ed13b8153d"},"cell_type":"code","source":"iobj=train_metadata[(train_metadata['target']==6) & (train_metadata['ddf']==0)].index[5]\n\nparams = fit_microlens(times_list[iobj],flux_list[iobj],err_list[iobj]) \nplot_microlens(times_list[iobj],flux_list[iobj],err_list[iobj],params[:2])\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5decdd76357b9e3d83b2a3dd2853606d94490804"},"cell_type":"markdown","source":"# Proceeding to the fit on all galactic training data\n... this will take some time\n\nNote that we don't fit for extragalactic objects to save computation time."},{"metadata":{"trusted":false,"_uuid":"cc959f13ba4e30259056a6f2f80e7b0e45cc015f"},"cell_type":"code","source":"def make_microparams():\n    microparams_list=[0]*len(train_metadata)\n    \n    for iobj in train_metadata[(train_metadata['hostgal_photoz']==0)].index:\n        if iobj%100==0: print(iobj)\n        microparams_list[iobj] = fit_microlens(times_list[iobj],flux_list[iobj],err_list[iobj])\n\n    return microparams_list","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"45c220887688adb1c60a7c73a736d5a0c6916a54"},"cell_type":"markdown","source":"A few warnings will get issued : those are cases of failed fits...\n    \nTime for a coffee... or to go runing on a decent multicore machine."},{"metadata":{"trusted":false,"_uuid":"a6969fb71ccaa3978079bc6accd0b63b7c96e5f0"},"cell_type":"code","source":"microparams=make_microparams()\n\n# then put things into array variables\nmicroamp = np.zeros((len(train_metadata),6))\nmicrobase = np.zeros((len(train_metadata),6))\nmicrot0 = np.zeros(len(train_metadata))\nmicrote = np.zeros(len(train_metadata))\nmicrowings = np.zeros(len(train_metadata))\nmicrochi2 = np.zeros(len(train_metadata))\nmicrondf = np.zeros(len(train_metadata))\nfor iobj in train_metadata[(train_metadata['hostgal_photoz']==0)].index:\n    microamp[iobj] = microparams[iobj][0][:,0]\n    microbase[iobj] = microparams[iobj][0][:,1]\n    microt0[iobj] = np.abs(microparams[iobj][1][0])\n    microte[iobj] = np.abs(microparams[iobj][1][1])\n    microwings[iobj] = microparams[iobj][1][2]\n    microchi2[iobj] = microparams[iobj][2]\n    microndf[iobj] = microparams[iobj][3]-15","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"25280503a51f24aa3b234b41df340bc44ebdeb89"},"cell_type":"markdown","source":"## Analysis of the results\n\nHere we present $chi2/ndf$ vs $t_E$ (time width of the event) on a log-log scale.\nDifferent classes will get different colors.\n\nAs expected, class 6 members have $chi2/ndf$ around 1.\n\nGood news is that this plot also helps to separate other galactic classes.\nAlso, the class 65 (orange) natural false-positive on chi2 have very different time width.\n\nI Hope that using the whole 15+1 extracted features will help your classifications !\n"},{"metadata":{"trusted":false,"_uuid":"3eadee1ba036035341e41db2331f7a3ff1f21357"},"cell_type":"code","source":"fig, axes = plt.subplots(2, 3, figsize=(12, 8))\nfor t in gal_target:\n#for t in [6]:\n    plotx=np.log(microchi2 / microndf)\n    ploty=np.log(microte)\n    ax = axes[0,0]\n    ax.scatter(plotx[t_type[t]],ploty[t_type[t]],color=t_colors[t],alpha=0.2)\n    ax = axes[0,1]\n    ax.scatter(plotx[ddf_type[t]],ploty[ddf_type[t]],color=t_colors[t],alpha=0.2)\n    ax = axes[0,2]\n    ax.scatter(plotx[wdf_type[t]],ploty[wdf_type[t]],color=t_colors[t],alpha=0.2)\n    ax = axes[1,0]\n    ax.scatter(plotx[t_type[t]],ploty[t_type[t]],color=t_colors[t],alpha=0.2)\n    ax = axes[1,1]\n    ax.scatter(plotx[ddf_type[t]],ploty[ddf_type[t]],color=t_colors[t],alpha=0.2)\n    ax = axes[1,2]\n    ax.scatter(plotx[wdf_type[t]],ploty[wdf_type[t]],color=t_colors[t],alpha=0.2)\n\n\nax = axes[0,0]\nax.set_title(\"All\")\nax = axes[0,1]\nax.set_title(\"DDF\")\nax = axes[0,2]\nax.set_title(\"WDF\")\nfor i in range(3):\n    ax = axes[1,i]\n    ax.set_xlim((-2,5))\n    ax.set_ylim((-10,7))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4d3a90dc0402bcde6709bd607635e6ebed3990db"},"cell_type":"markdown","source":"## Writing output to cvs file"},{"metadata":{"_uuid":"159d6e527ca71e3ec7dbce3a4bd58707acaa4a26"},"cell_type":"markdown","source":"Metrics are for now\n- maxsigma, frac3sigma, undersigma, fracunder, ls_frequency,ls_power, microt0, microte, microwings, microchi2, microndf\n\nAnd per band:\n- stds, medians,err_median, microbase, microamp "},{"metadata":{"trusted":false,"_uuid":"228cb819ae9a7545b7d7b1085073a7233761e08d"},"cell_type":"code","source":"\ntrain_metadata['micro_t0'] = microt0\ntrain_metadata['micro_te'] = microte\ntrain_metadata['micro_wings'] = microwings\ntrain_metadata['micro_chi2'] = microchi2\ntrain_metadata['micro_ndf'] = microndf\ntrain_metadata['microbase']=train_metadata.apply(lambda x : microbase[x.name],axis=1)\ntrain_metadata['microamp']=train_metadata.apply(lambda x : microamp[x.name],axis=1)\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"2bdd8a5f629bfdb5ca92075cf53118466028e9f6"},"cell_type":"code","source":"train_metadata.to_csv('training_set_metadata_microlensed.csv')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"cc124b5398cfca0f82e97c8716e68b7b3af4e6a4"},"cell_type":"markdown","source":"# This is the end !\n\nI hope that using the whole 15+1 extracted features will help you with your classifications.\nEnjoy !"}],"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.5"}},"nbformat":4,"nbformat_minor":1}