{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"scrolled":true},"cell_type":"markdown","source":"This notebook shows basic process of salt-2 template feature described [here](https://www.kaggle.com/c/PLAsTiCC-2018/discussion/75222). See also [official docs](https://sncosmo.readthedocs.io/en/v1.6.x/examples/plot_lc_fit.html) if you want to know in detail."},{"metadata":{"trusted":true,"_uuid":"80e0dcd42e1523d854b35e0fff761fd535c52b08"},"cell_type":"code","source":"import pandas as pd\nimport numpy as np\n!pip install iminuit","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"import sncosmo\nimport pandas as pd\nimport time\nfrom astropy.table import Table\nfrom astropy import wcs, units as u\nfrom sncosmo.bandpasses import read_bandpass\nfrom contextlib import contextmanager\n%matplotlib inline\n\n@contextmanager\ndef timer(name):\n    s = time.time()\n    yield\n    \n    print('[{}] {}'.format(time.time() - s, name))\n\nwith timer('load data'):\n    lc = pd.read_csv('../input/training_set.csv', nrows=10000)\n    meta = pd.read_csv('../input/training_set_metadata.csv')\n    meta.set_index('object_id', inplace=True)\n\n# only use data with signal-to-noise ratio (flux / flux_err) greater than this value\nminsnr = 3","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a38e0bdfd23a4f8269bcd75a350cb90fd3f66027"},"cell_type":"markdown","source":"## Specify a type of template\nAt first, you need to specify a source template (which type of light-curve you want to fit). There are a lot of available source templates (you can see a full list in [source](https://github.com/sncosmo/sncosmo/blob/master/sncosmo/builtins.py)). In this notebook, I choose salt2-extended template."},{"metadata":{"trusted":true,"_uuid":"711abcff29d9ec5da9cf251821cc438168f8e01e","scrolled":false},"cell_type":"code","source":"# template to use\nmodel_type = 'salt2-extended'\nmodel = sncosmo.Model(source=model_type)\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5f5975de0a2c8fd35ee50f7246633b84fcb54df9"},"cell_type":"markdown","source":"## Preprocessing\nYou also need to specify passband to tell sncosmo about observation wavelength. sncosmo already has passband for lsst, so it's process looks quite simple :)"},{"metadata":{"trusted":true,"_uuid":"e086933dd20fbd352c2c57a16aa5d895a72003ee"},"cell_type":"code","source":"\npassbands = ['lsstu','lsstg','lsstr','lssti','lsstz','lssty']\nwith timer('prep'):\n    lc['band'] = lc['passband'].apply(lambda x: passbands[x])\n    lc['zpsys'] = 'ab'\n    lc['zp'] = 25.0\n    ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e4666378cb1bf0990836b4f51c0f8eb8f830a933"},"cell_type":"markdown","source":"## FItting the light curve\nOk, let's start fitting process. "},{"metadata":{"trusted":true,"scrolled":false,"_uuid":"59e95e1187f9b7244b8c8fe7f638f83a996b9117"},"cell_type":"code","source":"object_id = 1598\n\ndata = Table.from_pandas(lc[lc.object_id == object_id])\n\nphotoz = meta.loc[object_id, 'hostgal_photoz']\nphotoz_err = meta.loc[object_id, 'hostgal_photoz_err']\n\n# run the fit\nwith timer('fit_lc'):\n    result, fitted_model = sncosmo.fit_lc(\n        data, model,\n        model.param_names,\n        # sometimes constant bound ('z':(0,1.4)) gives better result, so trying both seems better\n        bounds={'z':(max(1e-8,photoz-photoz_err), photoz+photoz_err)},\n        minsnr=minsnr)  # bounds on parameters\n\nsncosmo.plot_lc(data, model=fitted_model, errors=result.errors, xfigsize=10)\n\n\nprint('chisq:{}'.format(result.chisq))\nprint('hostgal_photoz: {}, hostgal_specz: {}, estimated z by the model: {}'.format(meta.loc[object_id,'hostgal_photoz'],\n                                                                                   meta.loc[object_id,'hostgal_specz'],\n                                                                                   result.parameters[0]))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7f37e055264d5a5751d38e3c77e3c05ead758b1d"},"cell_type":"markdown","source":"You can use estimated parameters and chi-sq value as feature. Be careful that it takes a long time (1sec/object) and don't forget to catch exception when you try it yourself. I also attached all template features on my discussion ([link](https://www.kaggle.com/c/PLAsTiCC-2018/discussion/75222)), so just download and use it if you don't want to wait. Enjoy."},{"metadata":{"trusted":true,"_uuid":"d9b6cb224577d88a834768eccef3ea0d16607078"},"cell_type":"code","source":"df = pd.DataFrame(columns=['chisq'] + model.param_names)\ndf.index.name = 'object_id'\ndf.loc[object_id] = [result.chisq] + list(result.parameters)\ndf","execution_count":null,"outputs":[]}],"metadata":{"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"}},"nbformat":4,"nbformat_minor":1}