{"cells":[{"metadata":{"_uuid":"178274ed2d0dc2f63afe6e396454786a386290a0"},"cell_type":"markdown","source":"Light curve can tell you a lot about the type of variable object, especially when the object is periodic. This kernel is about to examine, how to extract additional features from the lightcurves using periodograms and phase curves.\n\nWe will use scipy.signal.lobscargle which is something like fourier analysis for unevenly distributed data. Also known as [Least-squares spectral analysis](https://en.wikipedia.org/wiki/Least-squares_spectral_analysis) (LSSA)."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport gc\nfrom matplotlib import pyplot as plt\nfrom scipy.signal import lombscargle\nimport math\nfrom tqdm import tqdm","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"86571ec233921382674f96d7f4855b6348c30c26"},"cell_type":"code","source":"# some help functions\n# angular frequency to period\ndef freq2Period(w):\n    return 2 * math.pi / w\n# period to angular frequency\ndef period2Freq(T):\n    return 2 * math.pi / T","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"gc.enable()\ntrain = pd.read_csv('../input/training_set.csv')\nprint(train['object_id'].unique())\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fd987d675d92cf6623ab6f716d9db8abaec5b971"},"cell_type":"markdown","source":"My basic idea is simple: take lightcurves in different bands and normalize them, so that they fit together. Most periodic changes happen in all bands. The difference is, bands are usually shifted from each other (they have different **mean**), and the amplitude of changes can be also different (they have different **standard deviation**). Normalizing bands will make most of the well-behaved variable object fit together.\n\n*Note: this is just some rough approximation, which may not work for exotic (mainly extra-galactic) or extreme objects (like active black holes or whatever) - in general any object, whose light curves are not similar enough in different bands. But the information about such a mismatch is valuable for classifications by itself. For such objects, you could do periodogram for each band and then merge the results. But not know, maybe in future versions...*"},{"metadata":{"trusted":true,"_uuid":"f3094b696bb7709704aa104581d4232f9f7f120a"},"cell_type":"code","source":"# get data and normalize bands\ndef processData(train, object_id):\n    \n    #load data for given object\n    X = train.loc[train['object_id'] == object_id]\n    x = np.array(X['mjd'].values)\n    y = np.array(X['flux'].values)\n    passband = np.array(X['passband'].values)\n    \n    # normalize bands\n    for i in np.unique(passband):\n        yy = y[np.where(passband==i)]\n        mean = np.mean(yy)\n        std = np.std(yy)\n        y[np.where(passband==i)] = (yy - mean)/std\n    \n    return x, y, passband","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"95c59ced0f8a01cbec4a0645fb88b77caffcc892"},"cell_type":"markdown","source":"Let's get light curve for first object in dataset and plot it."},{"metadata":{"trusted":true,"_uuid":"c9a06e7098ffcd2ebd1525eb138481ec7b6dd86b"},"cell_type":"code","source":"x, y, passband = processData(train, 615)\nplt.scatter(x, y, c=passband)\nplt.xlabel('time (MJD)')\nplt.ylabel('Normalized flux')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"fe156ae920512aa8593b289bf12fa6e9436d74cd"},"cell_type":"markdown","source":"Big mess, right? Let's make some sense in it by periodograms.\nLoosely speaking, periodogram shows you something like the probability for each of possible periods (well, not exactly probability, that's why I call it power, but it's enough for basic understanding).\n\nThe computation takes ages and time is our most valuable resource (we have more than 3M objects in test set). Therefore I look for 5 most \"probable\" periods above some \"probability\" threshold and use them as new features."},{"metadata":{"trusted":true,"_uuid":"2534dfde797fb9c21eb2d840a43f4b87175a16c5"},"cell_type":"code","source":"# calculate periodogram\ndef getPeriodogram(x, y, steps = 10000, minPeriod = None, maxPeriod = None):\n    if not minPeriod:\n        minPeriod = 0.1 # for now, let's ignore very short periodic objects\n    if not maxPeriod:\n        maxPeriod = (np.max(x) - np.min(x))/2 # you cannot detect P > half of your observation period\n\n    maxFreq = np.log2(period2Freq(minPeriod))\n    minFreq = np.log2(period2Freq(maxPeriod))\n    f = np.power(2, np.linspace(minFreq,maxFreq, steps))\n    p = lombscargle(x,y,f,normalize=True)\n    return f, p","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"eaf956f841cbb53e9847102443740b6a1b480041"},"cell_type":"code","source":"%%time\nf,p = getPeriodogram(x, y, steps=20000)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e9d03e99e5f0c1688a0c75dde8dc1253df9dfba2"},"cell_type":"code","source":"plt.semilogx(freq2Period(f),p)\nplt.xlabel('Period (days)')\nplt.ylabel('Power')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2ad8069b6ddda4247135de696b45b4a567db1b89"},"cell_type":"markdown","source":"You can see how a typical periodogram of noisy data looks like. There is huge noise and several peak of different size. We use only 10000 steps for first pass, which means we may not always hit the peak exactly. Let's take all peak candidates of certain height (let's say power > 0.3) and examine them further."},{"metadata":{"trusted":true,"_uuid":"803be9b8b091d23acb396c9b1149ab6f4565da43"},"cell_type":"code","source":"def findBestPeaks(x, y, F, P, threshold=0.3, n=5):\n    \n    # find peaks above threshold\n    indexes = np.where(P>threshold)[0]\n    # if nothing found, look at the highest peaks anyway\n    if len(indexes) == 0:\n        q = np.quantile(P, 0.9995)\n        indexes = np.where(P>q)[0]\n    \n    peaks = []\n    start = 0\n    end = 0\n    for i in indexes:\n        if i - end > 10:\n            peaks.append((start, end))\n            start = i\n            end = i\n        else:\n            end = i\n    \n    peaks.append((start, end))\n        \n    \n    # increase accuracy on the found peaks\n    results = []\n    for start, end in peaks:\n        if end > 0:\n            minPeriod = freq2Period(F[min(F.shape[0]-1, end+1)])\n            maxPeriod = freq2Period(F[max(start-1, 0)])\n            steps = int(100 * np.sqrt(end-start+1)) # the bigger the peak width, the more steps we want - but sensible (linear increase leads to long computation)\n            f, p = getPeriodogram(x, y, steps = steps, minPeriod=minPeriod, maxPeriod=maxPeriod)\n            results.append(np.array([freq2Period(f[np.argmax(p)]), np.max(p)]))\n\n    # sort by normalized periodogram score and return first n results\n    if results:\n        results = np.array(results)\n        results = results[np.flip(results[:,1].argsort())]\n    else:\n        results = np.array([freq2Period(F[np.argmax(P)]), np.max(P)]).reshape(1,2)\n    return results[0:n]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9815d9a36206db31051534d380e143874132a700"},"cell_type":"code","source":"%%time\nresults = findBestPeaks(x, y, f, p)\nprint('Period(days) Power')\nprint(results)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8af02a0fd118795310b7ad53a54e92898227f559"},"cell_type":"markdown","source":"We found 4 peaks - one with a very high power, three others with much lower one. Let's check the results visually:"},{"metadata":{"trusted":true,"_uuid":"b9765fd7332cbd3ad8bb12fc2f1f2782424fdadf"},"cell_type":"code","source":"plt.figure(figsize=(20,25))\n\nfor i in range(results.shape[0]):\n    plt.subplot(results.shape[0],2,i+1)\n    phase = x/results[i][0] % 1\n    plt.scatter(phase, y, c = passband, s=4)\n    plt.xlabel('Phase')\n    plt.ylabel('Normalized flux')\n    plt.title('Period: {:.4f}, power: {:.2f}'.format(results[i][0], results[i][1]))\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"48d3b9c4fc785b37def391b4adcf0286fd8510d2"},"cell_type":"markdown","source":"We can see that only first period makes sense and the rest are false positives. This one looks like some short-period variable star with nicely periodic changes.\n\nYou can use the found periods and their respective powers as a new feature.\n\nWith phase curve, you can also try to extract other features, e.g.:\n* shape\n* symmetry / assymetry of the curve\n* humps, double minimas / maximas\n* fit the phase curve with sin function and calculate residuals for each band - in combination with flux errors, it's a measure of how strong the periodicity is (some objects are nicely periodic, like this one, some are semi-periodic with each minimum/maximum slightly different, which makes the phase curve more noisy)"},{"metadata":{"_uuid":"827ad87dcbc7abd460bafcf4bd072feb3c0f3de4"},"cell_type":"markdown","source":"# Multiprocessing\n\nIn this section, I will try to develop speed optimized technique to precompute basic periodogram features."},{"metadata":{"trusted":true,"_uuid":"15b31a4c92d94f4dcb855d0399a2b61ca02834b3"},"cell_type":"code","source":"from multiprocessing import Pool\nimport multiprocessing as mp\n\nCORES = mp.cpu_count() #4\n\ndef getFeatures(object_id):\n    \n    x, y, passband = processData(train, object_id)\n    f,p = getPeriodogram(x, y)\n    peaks = findBestPeaks(x, y, f, p)\n    features = np.zeros((5,2))\n    features[:peaks.shape[0],:peaks.shape[1]] = peaks\n    \n    return np.append(np.array([object_id]), features.reshape(5*2))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3cedd74c9649d144c87e6df21a2383b3a9f435c5"},"cell_type":"code","source":"object_ids = train['object_id'].unique()[0:100]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"82881d48340efed610b2b9b4c03035c61c2ecaf4"},"cell_type":"markdown","source":"First, let's try to calculate without multi-cpu speed-up."},{"metadata":{"trusted":true,"_uuid":"b9ac9e2302eef15913c16f8d6bfa023feb850732"},"cell_type":"code","source":"%%time\nfeatures = []\nfor object_id in object_ids:\n    results = getFeatures(object_id)\n    features.append(results)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"89ffa044c74fbf82bda8f409fcced9db7d6bf821"},"cell_type":"markdown","source":"That's 0.284 second per star. Training set will then take 2200 seconds to calculate. Testing set would take 852000 seconds, ~10 days. Not good.\nLet's try 4 cores available in Kaggle Kernels."},{"metadata":{"trusted":true,"_uuid":"f1bae71b5ed58763a5fae294df75d4563e38689f"},"cell_type":"code","source":"%%time\np = Pool(CORES)\n\nresults = p.map(getFeatures, object_ids)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2643aa282720930b7ed98d54b8948b0404fd78c7"},"cell_type":"markdown","source":"That's 0.100 seconds per star. Training set will take around 780 seconds. Testing set 300 000 seconds (~3,5 days). With 4x more CPU power, we got roughly 2,8x speedup.\nIf this relation holds linearly with CPU power, we could get testing set computed in 10 hours on 32 CPU on google cloud.\nI will definitelly try and will make it public, if the new features prove to be benefitial on train/validation set."},{"metadata":{"_uuid":"b1cee95203e110957d8b7b6fd321f071f310a23d"},"cell_type":"markdown","source":"# Calculating training set"},{"metadata":{"trusted":true,"_uuid":"7c29eeab0fd894aba3df57548c09eef5327eac7f"},"cell_type":"code","source":"object_ids = train['object_id'].unique()\ncolumns = np.array(['id'])\nfor i in range(5):\n    period_str = 'period_'+str(i+1)\n    power_str = 'power_'+str(i+1)\n    columns = np.append(columns, np.array([period_str, power_str]))\n\nresults = p.map(getFeatures, object_ids)\n\noutput = pd.DataFrame(results, columns=columns)\noutput['id'] = output['id'].astype(np.int32)\noutput.to_csv('./train-periods.csv')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"eb1bf372763fb9db75c02df1b5fd12136dee7fca"},"cell_type":"markdown","source":"# Further development / ideas:\n\n* *the speed is a key.  The process to compute periodograms and features for 3M+ dataset cannot take ages. I.e. we need paralelization and smart optimization of the number steps in periodogram search*\n* sort the observations by phase and band and feed the phase curve into RNN\n* feed the phase curve into CNN with channels = number of bands\n* look for the functions, that fit the curves well - some classes of variable objects can be fitted very precisely by a specific function, which can then help to identify 99 class (increase your chance to have it right)."}],"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}