{"cells":[{"metadata":{"_uuid":"674a818e56bbd5fb04d01097ffdebb36c8ad813f"},"cell_type":"markdown","source":"Here I will merge [my Something Different kernel ](https://www.kaggle.com/jimpsull/something-different) with [Alexander Firsov's](https://www.kaggle.com/alexfir) [fast-test-set-reading kernel](https://www.kaggle.com/alexfir/fast-test-set-reading)\n- There is an approach here that is different from much of what is in the publicly available kernels\n- I am not an astronomer.  A different perspective can sometimes be a good thing.  You may find that a few of these features can be combined with your features"},{"metadata":{"_uuid":"e0cd47528247c4227089d320912778d37bee12c5"},"cell_type":"markdown","source":"**Begin Alexander Firsov's code**"},{"metadata":{"trusted":true,"_uuid":"3e3bdc8a9b07fc6953b739c8108f9512a74a7221"},"cell_type":"markdown","source":"The kernel is an attempt to overcome kernels limitation (RAM, disk space) and allow fast test_set reading from kernels. It also can be helpful for owners of low memory computers.\n\nThe idea is to store each column as numpy array, then access them using numpy.memmap and combine column into pandas when needed. Reading is available by object id or in chunks. A chunk may contain multiple object ids, object ids are not spread between chunks.\n\nData preparation is done in two kernels because of max disk space limitation. See [kernel 1](https://www.kaggle.com/alexfir/test-set-columns-part-1) and [kernel 2](https://www.kaggle.com/alexfir/test-set-columns-part-2). Output of the kernels is added as input to this kernel."},{"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 os.path\nimport time\n\nCOLUMN_TO_TYPE = {\n    'object_id': np.int32,\n    'mjd': np.float32,\n    'passband': np.int8,\n    'flux': np.float32,\n    'flux_err': np.float32,\n    'detected': np.int8\n}\n\npart1_directory = r'../input/test-set-columns-part-1'\npart2_directory = r'../input/test-set-columns-part-2'\n\nCOLUMN_TO_FOLDER = {\n    'object_id': part2_directory,\n    'mjd': part2_directory,\n    'passband': part2_directory,\n    'flux': part1_directory,\n    'flux_err': part1_directory,\n    'detected': part1_directory\n}\n\n\ndef init_reading():\n    info = {}\n    object_range_file_path = os.path.join(COLUMN_TO_FOLDER['object_id'], 'object_id_range.h5')\n    print('reading {}'.format(object_range_file_path))\n    object_id_to_range = pd.read_hdf(object_range_file_path, 'data')\n    info['object_id_to_range'] = object_id_to_range\n    id_to_range = object_id_to_range.set_index('object_id')\n    info['object_id_start'] = id_to_range['start'].to_dict()\n    info['object_id_end'] = id_to_range['end'].to_dict()\n\n    records_number = object_id_to_range['end'].max()\n\n    mmaps = {}\n    for column, dtype in COLUMN_TO_TYPE.items():\n        directory = COLUMN_TO_FOLDER[column]\n        file_path = os.path.join(directory, 'test_set_{}.bin'.format(column))\n        mmap = np.memmap(file_path, dtype=COLUMN_TO_TYPE[column], mode='r', shape=(records_number,))\n        mmaps[column] = mmap\n\n    info['mmaps'] = mmaps\n\n    return info\n\n\ndef read_object_info(info, object_id, as_pandas=True, columns=None):\n    start = info['object_id_start'][object_id]\n    end = info['object_id_end'][object_id]\n\n    data = read_object_by_index_range(info, start, end, as_pandas, columns)\n    return data\n\n\ndef read_object_by_index_range(info, start, end, as_pandas=True, columns=None):\n    data = {}\n    for column, mmap in info['mmaps'].items():\n        if columns is None or column in columns:\n            data[column] = mmap[start: end]\n\n    if as_pandas:\n        data = pd.DataFrame(data)\n\n    return data\n\n\ndef get_chunks(info, chunk_size=1000):\n    object_id_to_range = info['object_id_to_range']\n    end_of_file_offset = object_id_to_range['end'].max()\n    start_offsets = object_id_to_range['start'].values[::chunk_size]\n    end_offsets = object_id_to_range['end'].values[(chunk_size - 1)::chunk_size]\n\n    end_offsets = list(end_offsets) + [end_of_file_offset]\n\n    chunks = pd.DataFrame({'start': start_offsets, 'end': end_offsets})\n    chunks = chunks.values.tolist()\n\n    return chunks","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1ffba8dd401bcdc68420cd963ecbaf80fc5039e3"},"cell_type":"markdown","source":"Before reading data, call init function and get info object that will be used later."},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"info = init_reading()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e29d3dc61dcdd7d73dec31fb8c2c258f1757561c"},"cell_type":"markdown","source":"Reading can be done just by object id, like in code below."},{"metadata":{"trusted":true,"_uuid":"06baac27f130843fd6a9131ca6dd151aa60b35eb"},"cell_type":"code","source":"# single object read as pandas object, first object\nobject_info13 = read_object_info(info, 13)\nobject_info13.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1fcf2f1a66ec9fa1f3821408716288993a386dba"},"cell_type":"code","source":"# last object from test_set\nobject_info104853812 = read_object_info(info, 104853812)\nobject_info104853812.tail()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"13eb5429762fd51fbb0c83724a183e4d71dfed71"},"cell_type":"markdown","source":"Data can be returned as pandas object or as dict of numpy arrays (numpy memmap). "},{"metadata":{"trusted":true,"_uuid":"c5449f7efacf30381ea25f1bddda5cd4db06d653"},"cell_type":"code","source":"object_info104853812 = read_object_info(info, 104853812, as_pandas=False, columns=['flux', 'flux_err'])\nobject_info104853812['flux'][-5:]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c95d418e637874e14e9dd5f53b8495cc4c07a25c"},"cell_type":"markdown","source":"Let's read all test set and perform some operation (we need doing something as data is not read until we use them) and measure required time."},{"metadata":{"trusted":true,"_uuid":"6fe270f4f00176c8037c1414dbd3a86bd15861e4"},"cell_type":"code","source":"    object_ids = info['object_id_to_range']['object_id'].values.tolist()\n    start_time = time.time()\n    records_read = 0\n    for object_id in object_ids:\n        object_info = read_object_info(info, object_id, columns=['flux'], as_pandas=False)\n        flux = object_info['flux']\n        records_read += flux.shape[0]\n        max = flux.max()\n\n    print(\"Single field reading took {:6.4f} secs, records = {}\".format((time.time() - start_time), records_read))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b336399748691048fa9c696bc7fbe5cf019bc680"},"cell_type":"markdown","source":"Not bad, but still only reading took a lot of time, even so single field was read and no pandas object was created.\nNow let's see what is performance of chunk reading. Each chunk contains 10000 objects and pandas object is created."},{"metadata":{"trusted":true,"_uuid":"c1df626563a8002a0d32e0df4d32780255106bfe"},"cell_type":"code","source":"    start_time = time.time()\n    records_read = 0\n    chunks = get_chunks(info, 10_000)\n    for index_start, index_end in chunks:\n        data = read_object_by_index_range(info, index_start, index_end)\n        flux = data['flux']\n        max = flux.max()\n        records_read += data.shape[0]\n\n    print(\"Chunks reading took {:6.4f} secs, records = {}\".format((time.time() - start_time), records_read))\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"baf8f798ab0aac4d5e9a40f37c762477c5104723"},"cell_type":"markdown","source":"**Begin Jim Sullivan's code**"},{"metadata":{"trusted":true,"_uuid":"c98983db0abdfb3a9ad7380b92b2c4add7a4fdbf"},"cell_type":"code","source":"import warnings\nfrom matplotlib import pyplot as plt\nimport copy\nimport scipy.stats as ss\n\nprint(os.getcwd())\nprint(os.listdir('../input/PLAsTiCC-2018'))\nprint('Getting Base (meta) dataFrame')\nbdf=pd.read_csv('../input/PLAsTiCC-2018/test_set_metadata.csv')\nprint(bdf.shape)\n\n#this is what Alex's code avoids having to do\n#print('Getting raw (lightcurve) dataFrame')\n#rdf=pd.read_csv('../input/training_set.csv')\n#print(rdf.shape)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"09beeab576d0a3c21a52246c0ffad86f86272065"},"cell_type":"markdown","source":"**I'd like to replace getLCDFArchive with getLCDF based on Alex's code**"},{"metadata":{"trusted":true,"_uuid":"c532b0df8533fd54d7f94e2648b2ac9e2018d725"},"cell_type":"code","source":"\ndef getLCDF(objid, info, show=False, lepb=2, hepb=5):\n    # single object read as pandas object, first object\n    # Alex's code for getting a light curve data frame\n    #original unfiltered light curve data frame\n    ulcdf = read_object_info(info, objid)\n    #unfiltered lcdf\n    #this may not be necessary and may be causing my slowdown as size increases\n    #ulcdf=copy.deepcopy(oulcdf)\n    filterLow = ulcdf.loc[:,'passband']==lepb\n    filterHigh = ulcdf.loc[:,'passband']==hepb\n    filterPb = filterLow | filterHigh\n    lcdf=ulcdf.loc[filterPb,:]\n    if show:\n        plt.plot(lcdf.loc[:,'mjd'],lcdf.loc[:,'flux'])\n        plt.show()\n        \n    return lcdf\n\nobjA=bdf.loc[0,'object_id']\nobjB=bdf.loc[3000000, 'object_id']\n\nalcdf=getLCDF(objA, info, show=True)\nblcdf=getLCDF(objB, info, show=True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e3f916df7457f541b24cbc01d022c6c6c30ff8a5"},"cell_type":"markdown","source":"**I hope the hardest part is over**\n- I've now integrated my code with Alex's code.\n- I can extract lightcurves and extract features without loading the whole test set into one dataFrame\n- More comments and markdown on my original Feature Extractor Kernel"},{"metadata":{"trusted":true,"_uuid":"9783c168541c495de12fd5ca301ee05219f30898"},"cell_type":"code","source":"def divideLcdf(elcdf, ddf, lep=2, hep=5):\n    #lcdf=copy.deepcopy(elcdf)\n    #this simple date cutting works on the ddf objects\n    if ddf:\n        minDate=np.min(elcdf.loc[:,'mjd'])\n        maxDate=np.max(elcdf.loc[:,'mjd'])\n\n        halfPoint=np.average([minDate, maxDate])\n        firstCut=np.average([minDate, halfPoint])\n        secondCut=np.average([halfPoint, maxDate])\n        minDate=np.min(elcdf.loc[:,'mjd'])\n        maxDate=np.max(elcdf.loc[:,'mjd'])\n\n        halfPoint=np.average([minDate, maxDate])\n        firstCut=np.average([minDate, halfPoint])\n        secondCut=np.average([halfPoint, maxDate])\n\n        #early\n        efilter=elcdf.loc[:,'mjd']<=firstCut\n        #late\n        lfilter=elcdf.loc[:,'mjd']>=secondCut\n        #mid\n        mfilter=(efilter | lfilter)==False\n    \n        edf=elcdf.loc[efilter]\n        mdf=elcdf.loc[mfilter]\n        ldf=elcdf.loc[lfilter]\n        \n        ledf = edf[edf['passband']==lep]\n        hedf = edf[edf['passband']==hep]\n        lmdf = mdf[mdf['passband']==lep]\n        hmdf = mdf[mdf['passband']==hep]\n        lldf = ldf[ldf['passband']==lep]\n        hldf = ldf[ldf['passband']==hep]\n    \n    #using the datecutting method often leads to zero population sizes with non-ddf objects\n    else:\n        \n        lowdf=elcdf[elcdf['passband']==lep]\n        highdf=elcdf[elcdf['passband']==hep]\n        lenLow=lowdf.shape[0]\n        lenHigh=highdf.shape[0]\n        \n        minSizeLow = int(lenLow / 3)\n        minSizeHigh = int(lenHigh / 3)\n        \n        lldf=lowdf.nlargest(minSizeLow, 'mjd')\n        hldf=highdf.nlargest(minSizeHigh, 'mjd')\n        ledf=lowdf.nsmallest(minSizeLow, 'mjd')\n        hedf=highdf.nsmallest(minSizeHigh, 'mjd')\n        lmdf=lowdf.nlargest(lenLow-minSizeLow, 'mjd').nsmallest(lenLow-2*minSizeLow, 'mjd')\n        hmdf=highdf.nlargest(lenHigh-minSizeHigh, 'mjd').nsmallest(lenHigh-2*minSizeHigh, 'mjd')\n    \n    return ledf, hedf, lmdf, hmdf, lldf, hldf\n\naalcdf, balcdf, calcdf, dalcdf, ealcdf, falcdf=divideLcdf(alcdf, True)\nprint(aalcdf.shape)\nprint(balcdf.shape)\nprint(calcdf.shape)\nprint(dalcdf.shape)\nprint(ealcdf.shape)\nprint(falcdf.shape)\n        ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"dbe36f7a65777e3f996d9c17a4d3183937f8c230"},"cell_type":"code","source":"def getSubPopFeats(pbdf, outSig=3.0):\n    \n    average=np.average(pbdf.loc[:,'flux'])\n    median=np.median(pbdf.loc[:,'flux'])\n    stdev=np.std(pbdf.loc[:,'flux'])\n    maxflux=np.max(pbdf.loc[:,'flux'])\n    minflux=np.min(pbdf.loc[:,'flux'])\n    stdflerr=np.std(pbdf.loc[:,'flux_err'])\n    medflerr=np.median(pbdf.loc[:,'flux_err'])\n    \n    #We want a means to extract the rate of decay or rise from minima or maxima\n    #This is grabbing within the population\n    #We also will look between populations\n    maxmjd=np.max(pbdf[pbdf['flux']==maxflux].loc[:,'mjd'])\n    minmjd=np.max(pbdf[pbdf['flux']==minflux].loc[:,'mjd'])\n    \n    #at what date does the max occur?\n    aftmaxdf=pbdf[pbdf['mjd']>maxmjd]\n    \n    #if there are data points after the max, what is the value and date of the lowest?\n    if aftmaxdf.shape[0]>0:\n        minaft=np.min(aftmaxdf.loc[:,'flux'])\n        aftminmjd=np.min(aftmaxdf[aftmaxdf['flux']==minaft].loc[:,'mjd'])\n        #(val at t0 - val at t1) / (t0 - t1) sb neg\n        decaySlope=(maxflux-minaft)/(maxmjd-aftminmjd)\n    \n    else:\n        decaySlope=0\n        \n    aftmindf=pbdf[pbdf['mjd']<minmjd]\n    if aftmindf.shape[0]>0:\n        maxaft=np.max(aftmindf.loc[:,'flux'])\n        aftmaxmjd=np.max(aftmindf[aftmindf['flux']==maxaft].loc[:,'mjd'])\n        #(val at t0 - val at t1) / (t0 - t1) sb pos\n        riseSlope=(minflux - maxaft)/(aftmaxmjd-minmjd)\n    \n    else:\n        riseSlope=0\n        \n    return average, stdev, median, medflerr, stdflerr, maxflux, \\\n            maxmjd, decaySlope, minflux, minmjd, riseSlope\n\na,b,c,d,e,f,g, h,i,j,k=getSubPopFeats(balcdf)\nprint(a)\nprint(b)\nprint(c)\nprint(d)\nprint(e)\nprint(f)\nprint(g)\nprint(h)\nprint(i)\nprint(j)\nprint(k)\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"61b3f065d8caa21a1855edd63a06069d56ff164c"},"cell_type":"markdown","source":"### The processing isn't done, but we will do as much as possible array style later\n- TBD how much slowdown this change causes.  I am now returning 12 instead of 7\n- Want to get decay rates within subsection in case peak at last subsection"},{"metadata":{"trusted":true,"_uuid":"c776983d5f9cab669f0a757a51d98164b009ebc2"},"cell_type":"code","source":"def processLc(objid, elcdf, ddf, lep=2, hep=5):\n    \n    lcdf=copy.deepcopy(elcdf)\n    \n    #feature borrowed from Grzegorz Sionkowski (../sionek)\n    #dt[detected==1, mjd_diff:=max(mjd)-min(mjd), by=object_id]\n    #detectMjds=elcdf[elcdf['detected']==1].loc[:,'mjd']\n    #deltaDetect=np.max(detectMjds) - np.min(detectMjds)\n    \n    #divide the incoming light curve to 6 subpopulations\n    ledf, hedf, lmdf, hmdf, lldf, hldf=divideLcdf(lcdf, ddf,lep=lep, hep=hep)\n    #return average, stdev, median, medflerr, stdflerr, maxflux, \\\n    #        maxmjd, decayslope, minflux, minmjd, riseSlope\n    \n    leavg, lestd, lemed, lemfl, lesfl, lemax, lemxd, ledsl, lemin, lemnd, lersl=getSubPopFeats(ledf)\n    heavg, hestd, hemed, hemfl, hesfl, hemax, hemxd, hedsl, hemin, hemnd, hersl=getSubPopFeats(hedf)\n    lmavg, lmstd, lmmed, lmmfl, lmsfl, lmmax, lmmxd, lmdsl, lmmin, lmmnd, lmrsl=getSubPopFeats(lmdf)\n    hmavg, hmstd, hmmed, hmmfl, hmsfl, hmmax, hmmxd, hmdsl, hmmin, hmmnd, hmrsl=getSubPopFeats(hmdf)\n    llavg, llstd, llmed, llmfl, llsfl, llmax, llmxd, lldsl, llmin, llmnd, llrsl=getSubPopFeats(lldf)\n    hlavg, hlstd, hlmed, hlmfl, hlsfl, hlmax, hlmxd, hldsl, hlmin, hlmnd, hlrsl=getSubPopFeats(hldf)\n    \n    \n    feats= [objid, leavg, lestd, lemed, lemfl, lesfl, lemax, \n            lemxd, ledsl, lemin, lemnd, lersl,\n            heavg, hestd, hemed, hemfl, hesfl, hemax, hemxd,\n            hedsl, hemin, hemnd, hersl,\n            lmavg, lmstd, lmmed, lmmfl, lmsfl, lmmax, lmmxd,\n            lmdsl, lmmin, lmmnd, lmrsl,\n            hmavg, hmstd, hmmed, hmmfl, hmsfl, hmmax, hmmxd, \n            hmdsl, hmmin, hmmnd, hmrsl,\n            llavg, llstd, llmed, llmfl, llsfl, llmax, llmxd, \n            lldsl, llmin, llmnd, llrsl,\n            hlavg, hlstd, hlmed, hlmfl, hlsfl, hlmax, hlmxd, \n            hldsl, hlmin, hlmnd, hlrsl]\n    \n    return feats\n\nfeats=processLc(objB, blcdf, bdf.loc[3000000,'ddf'])\n\nprint(feats)\n    ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"729a1e465c3649caceb2461b4c3da66913fb2751"},"cell_type":"markdown","source":"https://stackoverflow.com/questions/41888080/python-efficient-way-to-add-rows-to-dataframe#comment70958764_41888080\n\n\n\nfor row in iterable_object:\n    csv_writer.writerow(row)\n\nreturn df"},{"metadata":{"trusted":true,"_uuid":"d30684b381fe8f66674d7c21dc40702cfbeb9c71"},"cell_type":"code","source":"from io import StringIO\nfrom csv import writer \nimport time\n\ndef writeAChunk(firstRecord, lastRecord, info, statusFreq=500):\n    output = StringIO()\n    csv_writer = writer(output)\n\n    fdf=pd.DataFrame(columns=['object_id', 'leavg', 'lestd', 'lemed', 'lemfl', 'lesfl', 'lemax', \n                'lemxd', 'ledsl', 'lemin', 'lemnd', 'lersl',\n                'heavg', 'hestd', 'hemed', 'hemfl', 'hesfl', 'hemax', 'hemxd',\n                'hedsl', 'hemin', 'hemnd', 'hersl',\n                'lmavg', 'lmstd', 'lmmed', 'lmmfl', 'lmsfl', 'lmmax', 'lmmxd',\n                'lmdsl', 'lmmin', 'lmmnd', 'lmrsl',\n                'hmavg', 'hmstd', 'hmmed', 'hmmfl', 'hmsfl', 'hmmax', 'hmmxd', \n                'hmdsl', 'hmmin', 'hmmnd', 'hmrsl',\n                'llavg', 'llstd', 'llmed', 'llmfl', 'llsfl', 'llmax', 'llmxd', \n                'lldsl', 'llmin', 'llmnd', 'llrsl',\n                'hlavg', 'hlstd', 'hlmed', 'hlmfl', 'hlsfl', 'hlmax', 'hlmxd', \n                'hldsl', 'hlmin', 'hlmnd', 'hlrsl'])\n\n    theColumns=fdf.columns\n    \n    csv_writer.writerow(theColumns)\n    started=time.time()\n    for rindex in range(firstRecord, lastRecord):\n        #if you want to monitor progress\n        #ddf 18 sec per 100 on my macAir\n        #non ddf 25 sec per 100 on my macAir\n        if rindex%statusFreq==(statusFreq-1):\n            print(rindex)\n            print(\"Processing took {:6.4f} secs, records = {}\".format((time.time() - started), statusFreq))\n            started=time.time()\n            #fdf=pd.merge(fdf, tdf, on='key')\n        objid = bdf.loc[rindex,'object_id']\n        ddf=bdf.loc[rindex,'ddf']==1\n        #ig=bdf.loc[rindex,'hostgal_specz']==0\n        lcdf=getLCDF(objid, info)\n        feats=processLc(objid, lcdf, ddf)\n        #fdf.loc[rindex,:]=feats\n        csv_writer.writerow(feats)\n\n    output.seek(0) # we need to get back to the start of the BytesIO\n    chdf = pd.read_csv(output)\n    chdf.columns=theColumns\n    \n    return chdf\n\ntheColumns=['object_id', 'leavg', 'lestd', 'lemed', 'lemfl', 'lesfl', 'lemax', \n                'lemxd', 'ledsl', 'lemin', 'lemnd', 'lersl',\n                'heavg', 'hestd', 'hemed', 'hemfl', 'hesfl', 'hemax', 'hemxd',\n                'hedsl', 'hemin', 'hemnd', 'hersl',\n                'lmavg', 'lmstd', 'lmmed', 'lmmfl', 'lmsfl', 'lmmax', 'lmmxd',\n                'lmdsl', 'lmmin', 'lmmnd', 'lmrsl',\n                'hmavg', 'hmstd', 'hmmed', 'hmmfl', 'hmsfl', 'hmmax', 'hmmxd', \n                'hmdsl', 'hmmin', 'hmmnd', 'hmrsl',\n                'llavg', 'llstd', 'llmed', 'llmfl', 'llsfl', 'llmax', 'llmxd', \n                'lldsl', 'llmin', 'llmnd', 'llrsl',\n                'hlavg', 'hlstd', 'hlmed', 'hlmfl', 'hlsfl', 'hlmax', 'hlmxd', \n                'hldsl', 'hlmin', 'hlmnd', 'hlrsl']\n\nfdf=pd.DataFrame(columns=theColumns)\nchunksize=700\nfirstLoop=0\nlastLoop=100\nloops=lastLoop-firstLoop\nveryFirstRow=firstLoop*chunksize\nveryLastRow=lastLoop*chunksize-1\nfor i in range(firstLoop, lastLoop):\n    startRow=i*chunksize\n    stopRow=(i+1)*chunksize\n    chdf=writeAChunk(startRow, stopRow, info, statusFreq=int(chunksize/2))\n    fdf= pd.concat([fdf, chdf])\n    print(fdf.shape)\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c8a6acde0444cd20fcf96ae6a15605b2679ca583"},"cell_type":"markdown","source":"**Merge features with header information**"},{"metadata":{"trusted":true,"_uuid":"3f000a29f7c240b1b69ba7d448e68aed092bc5a3"},"cell_type":"code","source":"print(bdf.shape)\nprint(fdf.shape)\n\n#.merge complained without doing this\nfdf.index=range(veryFirstRow,veryLastRow+1)\nfdf.loc[:,'object_id']=fdf.loc[:,'object_id'].astype(str)\nbdf.loc[:,'object_id']=bdf.loc[:,'object_id'].astype(str)\n\nmdf=fdf.merge(bdf, on='object_id', how='left')\nprint(mdf.shape)\nmdf.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"4956be0e29553d28406c71af3bc69108dc3d9940"},"cell_type":"markdown","source":"**Psuedo box-plot**\n- checking to see if populations overlap one another\n- transients can slope up, down, or peak in middle\n- but at least one of early, middle, late will look different from others\n- especially if repeated in both passbands"},{"metadata":{"trusted":true,"_uuid":"d7b57d6abbb3c6d54634f510a53b5f0deb419cdc"},"cell_type":"code","source":"def testForOutlier(bdf, energy='high', sigmas=1.0):\n    \n    if energy=='high':\n        valCols=['heavg', 'hmavg', 'hlavg']\n        sigCols=['hestd', 'hmstd', 'hlstd']\n    else:\n        valCols=['leavg', 'lmavg', 'llavg']\n        sigCols=['lestd', 'lmstd', 'llstd']\n    \n    fdf=copy.deepcopy(bdf)\n    \n\n    \n    fdf.loc[:,energy + 'Energy_transitory_' + str(round(sigmas,1)) + '_TF']=False\n    for i in range(len(valCols)):\n        fdf.loc[:,'min' + str(valCols[i])] = fdf.loc[:,valCols[i]] - sigmas*fdf.loc[:,sigCols[i]]\n        fdf.loc[:,'max' + str(valCols[i])] = fdf.loc[:,valCols[i]] + sigmas*fdf.loc[:,sigCols[i]]\n    \n    for i in range(len(valCols)):\n        #fdf.loc[:,'earlySet']=range(fdf.loc[:,'minX100' + str(valCols[0])],fdf.loc[:, 'maxX100' + str(valCols[0])])\n        #earlyMaxLessThanMedMin\n        for j in range(len(valCols)):\n            if j!=i:\n                maxFailsOverlap=fdf.loc[:,'max' + str(valCols[i])]<fdf.loc[:,'min' + str(valCols[j])]\n                minFailsOverlap=fdf.loc[:,'min' + str(valCols[i])]>fdf.loc[:,'max' + str(valCols[j])]\n                fdf.loc[:,energy + 'Energy_transitory_' + str(round(sigmas,1)) + '_TF']= \\\n                (fdf.loc[:,energy + 'Energy_transitory_' + str(round(sigmas,1)) + '_TF'] | minFailsOverlap | maxFailsOverlap)\n    \n    for i in range(len(valCols)):\n        fdf=fdf.drop('min' + str(valCols[i]), axis=1)\n        fdf=fdf.drop('max' + str(valCols[i]), axis=1)\n        \n    return fdf\n\n\nenergy='high'\nsigmas=1.0\nfdf=testForOutlier(mdf)\nfdf.shape\nprint(fdf.loc[:,energy + 'Energy_transitory_' + str(round(sigmas,1)) + '_TF'].sum())\n\nsigmas=1.5\nfdf=testForOutlier(fdf, energy=energy, sigmas=sigmas)\nprint(fdf.loc[:,energy + 'Energy_transitory_' + str(round(sigmas,1)) + '_TF'].sum())\n\nenergy='low'\nsigmas=1.0\nfdf=testForOutlier(fdf, energy=energy, sigmas=sigmas)\nprint(fdf.loc[:,energy + 'Energy_transitory_' + str(round(sigmas,1)) + '_TF'].sum())\n\nsigmas=1.5\nfdf=testForOutlier(fdf, energy=energy, sigmas=sigmas)\nprint(fdf.loc[:,energy + 'Energy_transitory_' + str(round(sigmas,1)) + '_TF'].sum())\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"eb17d88fea668933a2e2fb5cde16ef28e729be4f"},"cell_type":"markdown","source":"**Consolidate outlier information**"},{"metadata":{"trusted":true,"_uuid":"464994bdf482c12b6cdd229f64164561a190a152"},"cell_type":"code","source":"\n\nfdf.loc[:,'outlierString']=fdf.loc[:,'highEnergy_transitory_1.5_TF'].astype(str) + \\\n                             fdf.loc[:,'highEnergy_transitory_1.0_TF'].astype(str) + \\\n                             fdf.loc[:,'lowEnergy_transitory_1.5_TF'].astype(str) + \\\n                             fdf.loc[:,'lowEnergy_transitory_1.0_TF'].astype(str)\n\n\ndef getOutlierScore(row):\n    tdict={'TrueTrueTrueTrue':8, 'FalseTrueTrueTrue':7, 'TrueTrueFalseTrue':7,\n       'FalseTrueFalseTrue':6, 'FalseFalseTrueTrue':3, 'TrueTrueFalseFalse':3,\n       'FalseTrueFalseFalse':3, 'FalseFalseFalseTrue':3, 'FalseFalseFalseFalse':0}\n    \n    #tdict could be a parameter\n    #both very significant - 8\n    #one fairly, one very - 7\n    #both fairly, 6\n\n    #the reason for the dropoff from both fairly to one very is that a deviation in one band could be measurement error\n    #the reason it doesn't drop lower is that the lack of deviation in one band could ALSO be measurement error\n    #the reason one fairly = one very is that measurement errors are likely to have high magnitude\n    #one very 3\n    #one fairly 3\n    #none 0\n\n    return tdict[row['outlierString']]\n\nfdf['outlierScore']=fdf.apply(getOutlierScore, axis=1)\n    \nfdf=fdf.drop('outlierString', axis=1)\n\n#fdf.to_csv('fastestFeatureTableWithTransitoryFlags.csv')\nprint(fdf.shape)\nprint(fdf.columns)\nprint(np.average(fdf.loc[:,'outlierScore']))\nprint(np.min(fdf.loc[:,'outlierScore']))\nprint(np.max(fdf.loc[:,'outlierScore']))\nprint(np.median(fdf.loc[:,'outlierScore']))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"af18b033fcc2ca65ff664b3e7c9455bf4a1bc5d8"},"cell_type":"markdown","source":"## Rate of decay\n- If outlierScore is zero we only care about the intraPop slopes and its amplitude related\n- If it is non-zero we want to see "},{"metadata":{"trusted":true,"_uuid":"b6dc273d70ec1a2c057ab9d8b6615448315717a1"},"cell_type":"code","source":"fdf['hipd']=0\nfdf['hipr']=0\nfdf['htpd']=0\nfdf['htpr']=0\n\nfdf['lipd']=0\nfdf['lipr']=0\nfdf['ltpd']=0\nfdf['ltpr']=0\n\noutlierFilter=(fdf['outlierScore']>0)\nprint(outlierFilter.sum())\n\nhipdFilter = (fdf['hmmax']>fdf['hemax']) & (fdf['hmmax']>fdf['hlmax']) & outlierFilter\nhtpdFilter = (fdf['hemax']>fdf['hmmax']) & (fdf['hemax']>fdf['hlmax']) & outlierFilter\nlipdFilter = (fdf['lmmax']>fdf['lemax']) & (fdf['lmmax']>fdf['llmax']) & outlierFilter\nltpdFilter = (fdf['lemax']>fdf['lmmax']) & (fdf['lemax']>fdf['llmax']) & outlierFilter\n\nprint(hipdFilter.sum())\nprint(htpdFilter.sum())\nprint(lipdFilter.sum())\nprint(ltpdFilter.sum())\n\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f3e3ebd0c6133dbe2b08a7c73970eb3e2f2da6a2"},"cell_type":"markdown","source":"## Use the filters to set decay values without ifs and loops"},{"metadata":{"trusted":true,"_uuid":"47c81bccaac4343efea0294306eb35fcc77760a0"},"cell_type":"code","source":"#peak to peak\n#these are light curves where the peak was in the middle\nfdf.loc[hipdFilter,'hipd']=(fdf.loc[hipdFilter,'hmmax']-fdf.loc[hipdFilter,'hlmax']) / \\\n     (fdf.loc[hipdFilter,'hmmxd']-fdf.loc[hipdFilter,'hlmxd'])\nfdf.loc[lipdFilter,'lipd']=(fdf.loc[lipdFilter,'lmmax']-fdf.loc[lipdFilter,'llmax']) / \\\n     (fdf.loc[lipdFilter,'lmmxd']-fdf.loc[lipdFilter,'llmxd'])\n\n#these are light curves where the peak was in the beginning\nfdf.loc[htpdFilter,'hipd']=(fdf.loc[htpdFilter,'hemax']-fdf.loc[htpdFilter,'hmmax']) / \\\n     (fdf.loc[htpdFilter,'hemxd']-fdf.loc[htpdFilter,'hmmxd'])\nfdf.loc[ltpdFilter,'lipd']=(fdf.loc[ltpdFilter,'lemax']-fdf.loc[ltpdFilter,'lmmax']) / \\\n     (fdf.loc[ltpdFilter,'lemxd']-fdf.loc[ltpdFilter,'lmmxd'])\nfdf.loc[htpdFilter,'htpd']=(fdf.loc[htpdFilter,'hmmax']-fdf.loc[htpdFilter,'hlmax']) / \\\n     (fdf.loc[htpdFilter,'hmmxd']-fdf.loc[htpdFilter,'hlmxd'])\nfdf.loc[ltpdFilter,'ltpd']=(fdf.loc[ltpdFilter,'lmmax']-fdf.loc[ltpdFilter,'llmax']) / \\\n     (fdf.loc[ltpdFilter,'lmmxd']-fdf.loc[ltpdFilter,'llmxd'])\n\n#print(fdf.loc[lipdFilter,'lipd'])\nfdf[outlierFilter].head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ad9091ef337bb0ad520751f6e0d67aaabf8a6e75"},"cell_type":"markdown","source":"## Repeat decay stats for rise\n- Not sure if this will have significance but not ruling it out\n- Many of the rises will have negative slopes - which is counterintuitive at first but is because most outliers are flux spikes so even min to min is a decay\n- If an outlier is actually a trough then the 'decay' could have a positive slope"},{"metadata":{"trusted":true,"_uuid":"136d359b67d2f9261a196b5d20f79db6c2e226fd"},"cell_type":"code","source":"hiprFilter = (fdf['hmmin']<fdf['hemin']) & (fdf['hmmin']<fdf['hlmin']) & outlierFilter\nhtprFilter = (fdf['hemin']<fdf['hmmin']) & (fdf['hemin']<fdf['hlmin']) & outlierFilter\nliprFilter = (fdf['lmmin']<fdf['lemin']) & (fdf['lmmin']<fdf['llmin']) & outlierFilter\nltprFilter = (fdf['lemin']<fdf['lmmin']) & (fdf['lemin']<fdf['llmin']) & outlierFilter\n\n#these are light curves where the peak was in the middle\nfdf.loc[hipdFilter,'hipr']=(fdf.loc[hipdFilter,'hmmin']-fdf.loc[hipdFilter,'hlmin']) / \\\n     (fdf.loc[hipdFilter,'hmmnd']-fdf.loc[hipdFilter,'hlmnd'])\nfdf.loc[lipdFilter,'lipr']=(fdf.loc[lipdFilter,'lmmin']-fdf.loc[lipdFilter,'llmin']) / \\\n     (fdf.loc[lipdFilter,'lmmnd']-fdf.loc[lipdFilter,'llmnd'])\n\n#these are light curves where the peak was in the beginning\nfdf.loc[htpdFilter,'hipr']=(fdf.loc[htpdFilter,'hemin']-fdf.loc[htpdFilter,'hmmin']) / \\\n     (fdf.loc[htpdFilter,'hemnd']-fdf.loc[htpdFilter,'hmmnd'])\nfdf.loc[ltpdFilter,'lipr']=(fdf.loc[ltpdFilter,'lemin']-fdf.loc[ltpdFilter,'lmmin']) / \\\n     (fdf.loc[ltpdFilter,'lemnd']-fdf.loc[ltpdFilter,'lmmnd'])\nfdf.loc[htpdFilter,'htpr']=(fdf.loc[htpdFilter,'hmmin']-fdf.loc[htpdFilter,'hlmin']) / \\\n     (fdf.loc[htpdFilter,'hmmnd']-fdf.loc[htpdFilter,'hlmnd'])\nfdf.loc[ltpdFilter,'ltpr']=(fdf.loc[ltpdFilter,'lmmin']-fdf.loc[ltpdFilter,'llmin']) / \\\n     (fdf.loc[ltpdFilter,'lmmnd']-fdf.loc[ltpdFilter,'llmnd'])\n\nfdf[outlierFilter].head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"68b943bc67f420f69cd8c9e9222c23a9ec6faef6"},"cell_type":"markdown","source":"save output"},{"metadata":{"trusted":true,"_uuid":"90669c0b82b44ae24d7a8b727865d951b8df47cf"},"cell_type":"code","source":"fdf.to_csv('testFeaturesFrom' + str(veryFirstRow) + 'TO' + str(veryLastRow) + '.csv')\nprint('saved file: testFeaturesFrom_' + str(veryFirstRow) + '_to_' + str(veryLastRow) + '.csv')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ef11d1f71dcccd2c4e819c9045663ed498508253"},"cell_type":"markdown","source":""}],"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}