{"cells":[{"metadata":{"trusted":true,"_uuid":"d01baac3757b22b92747599ed87524c1ddc0ad1d"},"cell_type":"markdown","source":"# Chapter 1: Introduction\n\nThe kernel demonstrates the clustering and expending we used to get to 7# place. Running this kernel on training event 1000 will score ~0.635 after the clustering stage, and 0.735 after expending stage.\n\nEvery stage takes about 8-10 min on Kaggle, and about half the time on my laptop. In the clustering part of the kernel the algorithm uses 5.500 pairs of z0, 1/2R (more on it below). By increasing the number to about 100.000 pairs the score will plateau at about 0.765 (after expending).\n\nHow does it work:\nIn each clustering loop the algorithm try to find all tracks originating from (0,0,z0) and with a radius of 1/(2*kt).\nIf a hit (x,y,z) is on a track the helix can be fully defined by the following features (1), (2)\n\nrr=(x^2+y^2)^0.5\ntheta_=arctan(y/x)\ndtheta = arcsin(kt*rr)\n(1)\tTheta=theta_+dtheta\n(2)\t(z-z0)*kt/dtheta\n\nTo solve the +pi,-pi problem we use sin, cos for theta.\nTo make (2) more uniform, we use arctan((z-z0)/(3.3*dtheta/kt))\n\nAfter calculating the features, the algorithm tries to cluster all the hits with the same features. This is done by sparse binning – using np.unique.\n\nThe disadvantage of sparse binning over dbscan is it’s sensitivity, the advantages are its speed and its sensitivity (almost no outliners).\n\nAfter clustering every hit choose if his cluster is good according to the clusters length.\nEvery 500 loops all hits belonging to tracks which are long enough are removed from the dataset\nIf two hits from the same detector are on the same track, the one which is closest to the track’s center of mass is chosen.\nThe z0, kt pairs are chosen randomly. While running, the algorithm changes the bin width and the length of the minimum track to be extracted from the dataset.\n\nExpending is done by selecting the un-clustered hits which are close to the center of mass of the track."},{"metadata":{"_uuid":"510cc36749be9fcd4cd768d3ab112bc7617e5e71"},"cell_type":"markdown","source":"# Chapter 2: Clustering"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport datetime\nimport os\n#print(os.listdir(\"../input\"))\nfrom ipywidgets import FloatProgress,FloatText\nfrom IPython.display import display\n\nimport time\nimport pdb\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.utils import shuffle\nfrom trackml.dataset import load_event\nfrom trackml.randomize import shuffle_hits\nfrom trackml.score import score_event\nfrom trackml.dataset import load_dataset\nfrom sklearn.cluster import DBSCAN\nfrom sklearn.preprocessing import StandardScaler\nfrom mpl_toolkits.mplot3d import Axes3D\nfrom matplotlib import cm\nfrom itertools import product\nimport gc\nimport cProfile\nfrom tqdm import tqdm\n\n%matplotlib inline\n#make wider graphs\nsns.set(rc={'figure.figsize':(12,5)})\nplt.figure(figsize=(12,5))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"1d60502a63403cc916c67b1137eeb36ca591c129"},"cell_type":"code","source":"#path = '../input/train_1/'\npath='../input/unzippedcernsample/'\nlabel_shift_M=1000000","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"06cbcc6c3210db2daf69b2ae244f8d0e07512bbc"},"cell_type":"code","source":"def create_one_event_submission(event_id, hits, labels):\n    sub_data = np.column_stack(([event_id]*len(hits), hits.hit_id.values, labels))\n    submission = pd.DataFrame(data=sub_data, columns=[\"event_id\", \"hit_id\", \"track_id\"]).astype(int)\n    return submission","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"225efbd270270d28321b3a0a5dd0fb0c6fed49b7"},"cell_type":"code","source":"def hit_score(res,truth):\n    tt=res.merge(truth[['hit_id','particle_id','weight']],on='hit_id',how='left')\n    un,inv,count = np.unique(tt['track_id'],return_inverse=True, return_counts=True)\n    tt['track_len']=count[inv]\n    un,inv,count = np.unique(tt['particle_id'],return_inverse=True, return_counts=True)\n    tt['real_track_len']=count[inv]\n    gp=tt.groupby('track_id')\n    gp=gp['particle_id'].value_counts().rename('par_freq').reset_index()\n    tt=tt.merge(gp,on=['track_id','particle_id'],how='left')\n    gp=gp.groupby('track_id').head(1)\n    gp=gp.rename(index=str, columns={'particle_id': 'common_particle_id'})\n    tt = tt.merge(gp.drop(['par_freq'],axis=1),on='track_id',how='left')\n    tt['to_score']=(2*tt['par_freq']>tt['track_len']) & (2*tt['par_freq']>tt['real_track_len'])\n    tt['score']=tt['weight']*tt['to_score']\n    return tt","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"fd911596627b039c2af3d2c06912ee8a5ee7246b"},"cell_type":"code","source":"def calc_features(hits,hipos,phik,double_sided=False):\n    \n    if not 'rr' in list(hits.columns):\n        hits['theta_']=np.arctan2(hits.y,hits.x)\n        hits['rr']=np.sqrt(np.square(hits.x)+np.square(hits.y))\n    ktrr=hits.rr*hipos.kt\n    hits['dtheta']=np.where((np.abs(ktrr)<1),np.arcsin(ktrr,where=(np.abs(ktrr)<1) ),ktrr)\n    hits['theta'] = hits.theta_+hits.dtheta\n    hits['phi'] = np.arctan2((hits.z-hipos.z0) ,phik*hits.dtheta/hipos.kt)*2.0/np.pi\n    hits['sint']=np.sin(hits['theta'])\n    hits['cost']=np.cos(hits['theta'])\n    hits['fault']=(np.abs(ktrr)>1).astype('int')\n    if double_sided:\n        hits['phi2'] = np.arctan2((hits.z-hipos.z0) ,phik*(np.pi-hits.dtheta)/hipos.kt)*2.0/np.pi\n        hits['theta2'] = hits.theta_+np.pi-hits.dtheta\n        hits['sint2']=np.sin(hits['theta2'])\n        hits['cost2']=np.cos(hits['theta2'])\n    return hits","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2283196746e1d9f8e4889505a5954b5b4c4971f3"},"cell_type":"code","source":"def tag_bins(cat):\n    un,inv,count = np.unique(cat,return_inverse=True, return_counts=True)\n    bin_tag=inv\n    bin_count=count[inv]\n    return bin_tag,bin_count","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2eabe9ad6b83c925244a344b13682180c49aa502"},"cell_type":"code","source":"def sparse_bin(features,bin_num,randomize=True,fault=None):\n    err=np.random.rand(features.shape[1])*randomize\n    cat=np.zeros(features.shape[0]).astype('int64')\n    factore=1\n    for i,feature in enumerate(features.columns):\n        cat=cat+(features[feature]*bin_num._asdict()[feature]+err[i]).astype('int64')*factore\n        factore=factore*(2*bin_num._asdict()[feature]+1)\n    if not fault is None:\n        cat=cat+(factore*features.index*fault).astype('int64')\n    return tag_bins(cat)\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"3c276e6101f8de47741f3e97873cc061bcb0a3f0"},"cell_type":"code","source":"def clustering(hits,stds,filters,phik=1.0,nu=500,weights=None,res=None,truth=None,history=None,pre_test_points=None):\n    start = time.time()\n    rest = hits.copy()\n    if weights is None:\n        weights={'phi':1, 'theta':0.15}\n    calc_score = not truth is None\n    if not history is None:\n        hist_list=[]\n    \n    if calc_score:\n        rest = rest.merge(truth[['hit_id','particle_id','weight']],on='hit_id',how='left')\n        dum,rest['particle_track_len']=tag_bins(rest['particle_id'])\n        score = 0 \n        hit_num=0\n        total_num=0\n        frs=FloatText(value=0, description=\"full score:\")\n        display(frs)\n        fs=FloatText(value=0, description=\"score:\")\n        display(fs)\n        fss=FloatText(value=0, description=\"s rate:\")\n        display(fss)\n        fsd=FloatText(value=0, description=\"add score:\")\n        display(fsd)\n    \n    ft = FloatText(value=rest.shape[0], description=\"Rest size:\")\n    display(ft)\n    fg = FloatText(value=rest.shape[0], description=\"Group size:\")\n    display(fg)\n    fgss = FloatText(description=\"filter:\")\n    display(fgss)\n    \n    if res is None:\n        rest['track_len']=1\n        rest['track_id']=-rest.index\n        rest['kt']=1e-6\n        rest['z0']=0\n    else:\n        rest=rest.merge(res[['hit_id','track_id','kt','z0']],on='hit_id',how='left')\n        dum,rest['track_len']=tag_bins(rest['track_id'])\n\n    res_list=[]\n    rest['sensor']=rest.volume_id+rest.layer_id*100+100000*rest.module_id\n    rest['layers']=rest.volume_id+rest.layer_id*100\n    if pre_test_points is None:\n        maxprog= filters.npoints.sum()\n    else:\n        maxprog = filters.shape[0]*pre_test_points.shape[0]\n    pbar = tqdm(total=maxprog,mininterval=5.0)\n    rest['pre_track_id']=rest['track_id']\n    p=-1\n    feature_cols=['theta','sint','cost','phi','rr','theta_','dtheta','fault']\n    for filt in filters.itertuples():\n        if pre_test_points is None:\n            test_points=pd.DataFrame()\n            for col in stds:\n                test_points[col] = np.random.normal(scale=stds[col],size=filt.npoints)\n        else:\n            test_points=pre_test_points.sample(frac=filt.npoints).reset_index(drop=True)\n        \n        for row in test_points.itertuples():\n            p=p+1\n            pbar.update()\n            calc_features(rest,row,phik)\n            rest['new_track_id'],rest['new_track_len']=sparse_bin(rest[['phi','sint','cost']],filt,fault=rest.fault)\n            rest['new_track_id']=rest['new_track_id']+(p+1)*label_shift_M\n            better = (rest.new_track_len>rest.track_len) & (rest.new_track_len<19)\n            rest['new_track_id']=rest['new_track_id'].where(better,rest.track_id)\n            dum,rest['new_track_len']=tag_bins(rest['new_track_id'])\n            better = (rest.new_track_len>rest.track_len) & (rest.new_track_len<19)\n            rest['track_id']=rest['track_id'].where(~better,rest['new_track_id']) \n            rest['track_len']=rest['track_len'].where(~better,rest['new_track_len'])\n            rest['kt']=rest['kt'].where(~better,row.kt)\n            rest['z0']=rest['z0'].where(~better,row.z0)\n            \n            if (((row.Index+1)%nu == 0) or (row.Index + 1 == test_points.shape[0])):\n                dum,rest['track_len']=tag_bins(rest['track_id'])\n                calc_features(rest,rest[['kt','z0']],phik)\n                gp = rest.groupby(['track_id']).agg({'phi': np.mean , \n                    'sint':np.mean, 'cost':np.mean}).rename(columns={ 'phi': 'mean_phi', \n                                'sint':'mean_sint', 'cost':'mean_cost'}).reset_index()\n                cols_to_drop = rest.columns.intersection(gp.columns).drop('track_id')\n                rest = rest.drop(cols_to_drop,axis=1).reset_index().merge(gp,on=['track_id'],how = 'left').set_index('index')\n                rest['dist'] = weights['theta']*np.square(rest.sint-rest.mean_sint)+ weights['theta']*np.square(rest.cost-rest.mean_cost)+ weights['phi']*np.square(rest.phi-rest.mean_phi)\n                rest=rest.sort_values('dist')\n                rest['closest']=rest.groupby(['track_id','sensor'])['dist'].cumcount()\n                rest['closest2']=rest.groupby(['track_id','layers'])['dist'].cumcount()\n                select = (rest['closest']!=0) | (rest['closest2']>2)  \n                rest['track_id']=rest['track_id'].where(~select,rest['pre_track_id'])\n                dum,rest['track_len']=tag_bins(rest['track_id'])\n                fgss.value=filt.phi\n                fg.value=filt.min_group\n                ft.value = rest[rest.track_len<=filt.min_group].shape[0]\n\n                select = (rest['track_len']>filt.min_group)\n                \n                #The next lines are just for printing\n                if calc_score:\n                    tm=rest[select]                   \n                    gp = tm.groupby(['track_id','particle_id'])['hit_id'].count().rename('par_count').reset_index()\n                    tm=tm.merge(gp,on=['track_id','particle_id'],how='left')\n                    gp = rest.groupby(['track_id','particle_id'])['hit_id'].count().rename('par_count').reset_index()\n                    rs=rest.merge(gp,on=['track_id','particle_id'],how='left')\n                    to_full_score=(rs.weight*((rs.par_count*2>rs.track_len) & (rs.par_count*2>rs.particle_track_len)))\n                    frs.value=to_full_score.sum()+fs.value\n                    to_score=(tm.weight*((tm.par_count*2>tm.track_len) & (tm.par_count*2>tm.particle_track_len)))\n                    hit_num=hit_num+(to_score>0).sum()\n                    total_num=total_num+tm.weight.sum()\n                    fs.value=fs.value+to_score.sum()\n                    fss.value=fs.value/total_num\n                    fsd.value=to_score.sum()\n                    gp = rest.groupby(['track_id','particle_id'])['hit_id'].count().rename('par_count').reset_index()\n                    rs=rest.merge(gp,on=['track_id','particle_id'],how='left')\n                    to_full_score=(rs.weight*((rs.par_count*2>rs.track_len) & (rs.par_count*2>rs.particle_track_len)))\n                    frs.value=to_full_score.sum()+fs.value-to_score.sum()\n                    if not history is None:\n                        hist_list.append(pd.DataFrame({'P':p,'ftheta':filt.phi,'added_score':to_score.sum(),'min_group':filt.min_group,\n                                                    'full_score':frs.value,'score':fsd.value,'correct':fss.value,\n                                                    'clustered':tm.shape[0],'left':rest.shape[0]-tm.shape[0]}, index=[0]))\n\n                #end of printing part \n                \n                tm=rest[select][['hit_id','track_id','kt','z0']]\n                res_list.append(tm)\n                rest = rest[~select]\n                dum,rest['track_len']=tag_bins(rest['track_id'])\n                rest['pre_track_id']=rest['track_id']\n\n    ft.value = rest.shape[0]\n    res_list.append(rest[['hit_id','track_id','kt','z0']].copy())\n    res = pd.concat(res_list, ignore_index=True)\n    pbar.close()\n    rest['track_id'],dum=tag_bins(rest['track_id'])\n     \n    if not history is None:\n        history.append(pd.concat(hist_list,ignore_index=False))\n    print ('took {:.5f} sec'.format(time.time()-start))\n    return res \n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"028e58874978f0689caa5d9e95198154ef003436"},"cell_type":"code","source":"def refine_hipos(res,hits,stds,nhipos,phik=3.3,weights=None): \n    cols=list(res.columns)\n    if weights is None:\n        weights={'theta':0.15, 'phi':1.0}\n\n    groups = res.merge(hits,on='hit_id',how='left')\n    if not groups.columns.contains('kt'):\n        groups['kt']=0\n        groups['z0']=0\n        print(\"No kt's, calculating\")\n    calc_features(groups,groups[['kt','z0']],phik)\n\n    gp=groups.groupby('track_id').agg({'phi': np.std , 'sint' : np.std,\n            'cost' : np.std}).rename(columns={ 'phi': 'phi_std', \n            'sint' : 'sint_std', 'cost':'cost_std'}).reset_index()\n    groups=groups.merge(gp,on='track_id',how='left')\n    groups['theta_std']=np.sqrt(weights['theta']*np.square(groups.sint_std)+weights['theta']*np.square(groups.cost_std))\n    hipos=pd.DataFrame()\n    for col in stds:\n        hipos[col]=np.random.normal(scale=stds[col],size=nhipos)\n\n    for hipo in tqdm(hipos.itertuples(),total=nhipos):\n\n        groups['kt_new']=groups['kt']+hipo.kt\n        groups['z0_new']=groups['z0']+hipo.z0\n        calc_features(groups,groups[['kt_new','z0_new']].rename(columns={\"kt_new\": \"kt\", \"z0_new\": \"z0\"}),phik)\n        gp=groups.groupby('track_id').agg({'phi': np.std , 'sint' : np.std,\n            'cost' : np.std}).rename(columns={ 'phi': 'new_phi_std', \n            'sint' : 'new_sint_std', 'cost':'new_cost_std'}).reset_index()\n        groups=groups.merge(gp,on='track_id',how='left')\n        groups['new_theta_std']=np.sqrt(weights['theta']*np.square(groups.new_sint_std)+weights['theta']*np.square(groups.new_cost_std))\n\n        old_std=np.sqrt(np.square(groups.theta_std)+weights['phi']*np.square(groups.phi_std))\n        new_std=np.sqrt(np.square(groups.new_theta_std)+np.square(groups.new_phi_std))\n        cond=(old_std<=new_std) \n        groups['kt']=groups['kt'].where(cond,groups.kt_new)\n        groups['z0']=groups['z0'].where(cond,groups.z0_new)\n        groups['theta_std']=groups['theta_std'].where(cond,groups.new_theta_std)\n        groups['sint_std']=groups['sint_std'].where(cond,groups.new_sint_std)\n        groups['cost_std']=groups['cost_std'].where(cond,groups.new_cost_std)\n        groups['phi_std']=groups['phi_std'].where(cond,groups.new_phi_std)\n        groups=groups.drop(['new_theta_std','new_phi_std','new_sint_std','new_cost_std'],axis=1)\n\n        #pdb.set_trace()\n    to_return=groups[cols+['theta_std','phi_std','sint_std','cost_std']]\n    return to_return\n    \n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c1681033621feb237981cd38709db1fe91ab1620"},"cell_type":"code","source":"def expand_tracks(res,hits,min_track_len,max_track_len,max_expand,to_track_len,mstd=1.0,dstd=0.0,phik=3.3,max_dtheta=10,mstd_size=None,mstd_vol=None,drop=0,nhipo=1000,weights=None):\n    if weights is None:\n        weights={'theta':0.25, 'phi':1.0}\n\n    if mstd_size is None:\n        mstd_size=[0 for i in range(20)]\n    if mstd_vol is None:\n        mstd_vol={7:0,8:0,9:0,12:0,13:0,14:0,16:0,17:0,18:0}\n    gp=res.groupby('track_id').first().reset_index()\n    orig_hipo=gp[['track_id','kt','z0']]\n    eres=res.copy()\n    res_list=[]\n    stds={'kt':7e-5,'z0':0.8}\n    eres=refine_hipos(eres,hits,stds,nhipo,phik=phik,weights=weights)\n    dum,eres['track_len']=tag_bins(eres['track_id'])\n    eres['max_track_len']=np.clip(eres.track_len+max_expand,0,max_track_len) \n    eres['max_track_len']=2*(  eres['max_track_len']/2).astype('int')+1\n    eres=eres.sort_values('track_len')\n    eres = eres.merge(hits,on='hit_id',how='left')\n    eres['sensor']=eres.volume_id+eres.layer_id*100+100000*eres.module_id\n    group_sensors=eres.groupby('track_id').sensor.unique()\n    groups=eres[eres.track_len>min_track_len].groupby('track_id').first().reset_index().copy()\n    groups['order']=-groups.track_len \n    groups=groups.sort_values('order').reset_index(drop=True)\n    groups=groups.head(int((1.0-drop)*groups.shape[0])).copy()\n    select=eres.track_len<to_track_len\n    grouped=eres[~select]\n    regrouped=eres[select].copy()\n    regrouped['min_dist']=100\n    regrouped['new_track_len']=0\n    regrouped['new_track_id']=regrouped['track_id']\n    regrouped['new_kt']=regrouped['kt']\n    regrouped['new_z0']=regrouped['z0']\n    regrouped['new_max_size'] = max_track_len\n\n    f = FloatProgress(min=0, max=groups.shape[0], description='calculating:') # instantiate the bar\n    display(f) # display the bar\n\n    for group_tul in tqdm(groups.itertuples(),total=groups.shape[0]):\n        if group_tul.Index%20 ==0: f.value=group_tul.Index\n        if group_tul.track_len>=max_track_len: continue\n        group=eres[eres.track_id==group_tul.track_id].copy()\n        calc_features(group,group[['kt','z0']],phik)\n        group['abs_z']=np.abs(group.z)\n        group['abs_theta']=np.abs(group.theta)\n        phi_mean=group.phi.mean()\n        sint_mean=group.sint.mean()            \n        cost_mean=group.cost.mean()\n        max_z=group.abs_z.max()\n        max_theta=group.abs_theta.max()\n        regrouped['abs_z']=np.abs(regrouped.z)\n        calc_features(regrouped,group_tul,phik,double_sided=True)\n        regrouped['dist'] =np.sqrt(weights['theta']*np.square(regrouped.sint-sint_mean)+weights['theta']*np.square(regrouped.cost-cost_mean)+weights['phi']*np.square(regrouped.phi-phi_mean))\n        regrouped['dist2'] =np.sqrt(weights['theta']*np.square(regrouped.sint2-sint_mean)+weights['theta']*np.square(regrouped.cost2-cost_mean)+weights['phi']*np.square(regrouped.phi2-phi_mean))\n        select = (regrouped.abs_z>max_z)  & (max_dtheta >max_dtheta) & (regrouped.dist2<regrouped.dist)\n        regrouped['dist']=regrouped['dist'].where(~select,regrouped['dist2'])    \n        cmstd=regrouped.volume_id.map(mstd_vol)+mstd_size[group_tul.track_len]+mstd\n        if (dstd==0.0):\n            sdstd==group.dstd\n        else:\n            sdstd=dstd\n        better =( regrouped.dist<cmstd*sdstd) & ( regrouped.dist<regrouped.min_dist) & (~regrouped.sensor.isin(group_sensors.loc[group_tul.track_id]))\n        regrouped['min_dist']=np.where(better,regrouped.dist,regrouped.min_dist)\n        regrouped['new_track_id']=np.where(better,group_tul.track_id,regrouped.new_track_id)\n        regrouped['new_z0']=np.where(better,group_tul.z0,regrouped.new_z0)\n        regrouped['new_kt']=np.where(better,group_tul.kt,regrouped.new_kt)\n        regrouped['new_track_len']=np.where(better,group_tul.track_len,regrouped.new_track_len)\n        regrouped['new_max_size']=np.where(better,group_tul.max_track_len,regrouped.new_max_size)\n    f.value=group_tul.Index\n    regrouped=regrouped.sort_values('min_dist')\n    regrouped['closest']=regrouped.groupby('new_track_id')['min_dist'].cumcount()\n    better=regrouped.closest+regrouped.new_track_len>=regrouped.new_max_size\n    regrouped['track_id']=regrouped['track_id'].where(better,regrouped['new_track_id'])\n    res_list.append(grouped[['hit_id','track_id']])\n    res_list.append(regrouped[['hit_id','track_id']])\n    to_return=pd.concat(res_list)\n    to_return=to_return.merge(orig_hipo,on='track_id',how='left')\n    return to_return\n        \n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"af56879ed0a129fa102b13e2123f537662b547f1"},"cell_type":"code","source":"# the following 2 functions are taken from outrunner's kernel: https://www.kaggle.com/outrunner/trackml-2-solution-example\n\ndef get_event(event):\n    hits= pd.read_csv(path+'%s-hits.csv'%event)\n    cells= pd.read_csv(path+'%s-cells.csv'%event)\n    truth= pd.read_csv(path+'%s-truth.csv'%event)\n    particles= pd.read_csv(path+'%s-particles.csv'%event)\n    return hits, cells, particles, truth\n\ndef score_event_fast(truth, submission):\n    truth = truth[[\"hit_id\", \"particle_id\", \"weight\"]].merge(submission, how='left', on='hit_id')\n    df = truth.groupby(['track_id', 'particle_id']).hit_id.count().to_frame('count_both').reset_index()\n    truth = truth.merge(df, how='left', on=['track_id', 'particle_id'])\n    \n    df1 = df.groupby(['particle_id']).count_both.sum().to_frame('count_particle').reset_index()\n    truth = truth.merge(df1, how='left', on='particle_id')\n    df1 = df.groupby(['track_id']).count_both.sum().to_frame('count_track').reset_index()\n    truth = truth.merge(df1, how='left', on='track_id')\n    truth.count_both *= 2\n    score = truth[(truth.count_both > truth.count_particle) & (truth.count_both > truth.count_track)].weight.sum()\n    particles = truth[(truth.count_both > truth.count_particle) & (truth.count_both > truth.count_track)].particle_id.unique()\n\n    return score","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"00cff1de64d7a7f24df2e9430344f5ae90dc953c"},"cell_type":"code","source":"history=[]\nevent_num = 0\nevent_prefix = 'event00000100{}'.format(event_num)\n#hits, cells, particles, truth = load_event(os.path.join(path, event_prefix))\nhits, cells, particles, truth = get_event(event_prefix)\nweights={'pi':1,'theta':0.15}\nstds={'z0':7.5, 'kt':7.5e-4}\nd =    {'sint':[225,110,110,110,110,110],\n        'cost':[225,110,110,110,110,110],\n          'phi':[550,260,260,260,260,260],\n        'min_group':[11,11,10,9,8,7],\n        #'npoints':[50,50,50,50,50,50]}  # half minute run for testing\n        'npoints':[500,2000,1000,1000,500,500]}\n\n\nfilters=pd.DataFrame(d)\nnu=500\nresa=clustering(hits,stds,filters,phik=3.3,nu=nu,truth=truth,history=history)\nresa[\"event_id\"]=event_num\nscore = score_event_fast(truth, resa.rename(index=str, columns={\"label\": \"track_id\"}))\nprint(\"Your score: \", score)\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6e40624261df49241f12dbefdc111f361d2bda33"},"cell_type":"markdown","source":"# Chapter 3: Expand tracks"},{"metadata":{"trusted":true,"_uuid":"f17bf5518006da909b788b04fd6646317b65dc9b"},"cell_type":"code","source":"mstd_vol={7:0,8:0,9:0,12:2,13:1,14:2,16:3,17:2,18:3}\nmstd_size=[4,4,4,4,3,3,3,2,2,2,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]\nweights={'theta':0.1, 'phi':1}\nnresa=expand_tracks(resa,hits,5,16,5,7,mstd=8,dstd=0.00085,phik=3.3,max_dtheta=0.9*np.pi/2,mstd_vol=mstd_vol,mstd_size=mstd_size,weights=weights,nhipo=100)\nnresa['event_id']=0\nscore = score_event_fast(truth, nresa)\nprint(\"Your score: \", score)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"317879ec4824598fad62b7faba8f981eed065681"},"cell_type":"markdown","source":"# Chapter 4: Employ Machine Learning\n## 4.1 General strategy\n- Produce different submission candidates sub_1, sub_2, ..., sub_N\n- Create a machine learning model, which gives probabilities between 0 and 1 for each track candidate\n- Merge two submission candidates by assigning to each hit the track, which has higher probability\n- Merge the submission candidates successively (sub_1 and sub_2 to sub_12, sub_12 and sub_3 to sub_123, etc.) to get the final submission\n\n[Note: For our final solution, the methods described in this chapter gave around +0.01 to the LB score. Certain benefits from these methods have already been captured by the function, which expands the tracks.]\n\n## 4.2 Create machine learning model (LightGBM)\n- Training data: Use the truth file of the 13 train events 3-15 to get true tracks (target=1). To get wrong tracks (target=0) use generated clustering submissions on those latter events. In particular, we consider each track candidate, for which not all hits belong to the same particle_id, as a wrong track. Also, choose only tracks which have at least length 4, to slightly optimize compute time at negligible cost.\n- Use 13 features per track: \n - variance of x,y,z (these are the most important)\n - minimum of x,y,z\n - maximum of x,y,z\n - mean of z\n - volume_id of first hit \n - number of clusters per track (i.e. are there many hits, which are close together)\n - number of hits / number of clusters   \n- Validation data: Same as with training data, but use the 3 training events 0,1,2\n\n[Note: We also tried many more features (compare below dataframe df_train), but they gave only small additional gains, so we just kept those 13 features for simplicity. In fact, many other interesting features, such as number of hits, number of different volumes crossed etc are closely related to the above used ones.]\n\nWe have prepared the training and test data and load it directly from a pkl-file:"},{"metadata":{"trusted":true,"_uuid":"5800822ec8d5798c6bc9ecba2df5f20964fe13cc"},"cell_type":"code","source":"import pickle\ndef load_obj(filename):\n    with open(filename, 'rb') as f:\n        return pickle.load(f)\n    \ndf_train=load_obj('../input/trackml-validation-data-for-ml-pandas-df/df_train_v1.pkl')\ndf_test=load_obj('../input/trackml-validation-data-for-ml-pandas-df/df_test_v1.pkl')\ny_train=df_train.target.values\ny_test=df_test.target.values\nprint(\"The dataframe with all features:\")\ndisplay(df_train.head())\nprint(\"Features for each track:\",df_train.columns.values)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ce829f5fff7e54686b77149957379b04955511f1"},"cell_type":"markdown","source":"In the competition, we used roughly 250 events for training, but the additional improvement to just using 13 events is not too big. We now create the LightGBM model, using the mentioned training data and features:"},{"metadata":{"trusted":true,"_uuid":"a22e57031668f11d7069e34f781fce92de5cf135"},"cell_type":"code","source":"import lightgbm\ns=time.time()\n# choose which features of the tracks we want to use:\ncolumns=['svolume','nclusters', 'nhitspercluster', 'xmax','ymax','zmax', 'xmin','ymin','zmin', 'zmean', 'xvar','yvar','zvar']\nrounds=1000\nround_early_stop=50\nparameters = { 'subsample_for_bin':800, 'max_bin': 512, 'num_threads':8, \n               'application': 'binary','objective': 'binary','metric': 'auc','boosting': 'gbdt',\n               'num_leaves': 128,'feature_fraction': 0.7,'learning_rate': 0.05,'verbose': 0}\ntrain_data = lightgbm.Dataset(df_train[columns].values, label=y_train)\ntest_data = lightgbm.Dataset(df_test[columns].values, label=y_test)\nmodel = lightgbm.train(parameters,train_data,valid_sets=test_data,num_boost_round=rounds,early_stopping_rounds=round_early_stop,verbose_eval=50)\nprint('took',time.time()-s,'seconds')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"850aa36925328de14ac021516953e6be3180aa38"},"cell_type":"markdown","source":"## 4.3 Judge machine learning model\n\nWe doublecheck the model's performance, by calculating its precision, recall and accuracy on the validation set:\n\n[Note: For ML to be helpful in our situation, it needs to distinguish correct from wrong tracks with very high precision=true_positives/(false_positives+true_positives)). Also, it needs to do so for various sets of track candidates (especially such, which are generated if one tries to find tracks which originate far away from the origin; in those situations, often a lot of bad candidates are produced).]"},{"metadata":{"trusted":true,"_uuid":"8767bd6204e4233d21196288fa54cc62d06560ef"},"cell_type":"code","source":"def precision_and_recall(y_true, y_pred,threshold=0.5):\n    tp,fp,fn,tn=0,0,0,0\n\n    for i in range(0,len(y_true)):\n        if y_pred[i]>=threshold:\n            if y_true[i]>0:\n                tp+=1\n            else:\n                fp+=1\n        elif y_true[i]==0:\n            tn+=1\n        else:\n            fn+=1\n    precision=tp/(tp+fp) if (tp+fp != 0) else 0\n    recall=tp/(tp+fn) if (tp+fn != 0) else 0\n    accuracy=(tp+tn)/(tp+tn+fp+fn)\n    print('Threshold',threshold,' --- Precision: {:5.4f}, Recall: {:5.4f}, Accuracy: {:5.4f}'.format(precision,recall,accuracy))\n    return precision, recall, accuracy\n\ny_test_pred=model.predict(df_test[columns].values)\nprecision, recall, accuracy=precision_and_recall(y_test, y_test_pred,threshold=0.1)\nprecision, recall, accuracy=precision_and_recall(y_test, y_test_pred,threshold=0.5)\nprecision, recall, accuracy=precision_and_recall(y_test, y_test_pred,threshold=0.9)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2d682c11edcd597fc5f77a87fd4b872d76cd5043"},"cell_type":"markdown","source":"## 4.4 Use machine learning model\n\nLet us first load the event data of training event 0, as well as two generated submissions we have prepared:"},{"metadata":{"trusted":true,"_uuid":"2b7ce25ccbb90fc2da2fc560d9a48f5be46b1744"},"cell_type":"code","source":"### create one more submission\nhistory=[]\nresa2=clustering(hits,stds,filters,phik=3.3,nu=nu,truth=truth,history=history)\nresa2[\"event_id\"]=event_num\nscore = score_event_fast(truth, resa2.rename(index=str, columns={\"label\": \"track_id\"}))\nprint(\"Your score: \", score)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2ee989d56424c16256ca4646d9f1702316ec3c2c"},"cell_type":"markdown","source":"Calculate the probabilities for the tracks in those two submissions (small optimization: take also length of track into account, and after a couple of merged submissions, ask the probability of the track from the new submission to be at least C higher than the current probability; this latter option is not used in this kernel, but was used when merging >= 4 submissions)"},{"metadata":{"trusted":true,"_uuid":"1eb29b93df5606540e538cd94253a761ef5339af"},"cell_type":"code","source":"from tqdm import tqdm_notebook\nfrom sklearn.cluster.dbscan_ import dbscan\n\ndef get_features(track_hits,cluster_size=10):    \n    \"\"\"\n    Input: dataframe with hits of 1 track\n    Output: array with features of track\n    \"\"\"\n    nhits = len(track_hits)\n    svolume=track_hits['volume_id'].values.min()\n    X=np.column_stack([track_hits.x.values, track_hits.y.values, track_hits.z.values])\n    _, labels = dbscan(X, eps=cluster_size, min_samples=1, algorithm='ball_tree', metric='euclidean')\n    uniques = np.unique(labels)\n    nclusters = len(uniques)\n    nhitspercluster = nhits/nclusters\n    xmax=track_hits['x'].values.max()\n    xmin=track_hits['x'].values.min()\n    xvar=track_hits['x'].values.var()\n    ymax=track_hits['y'].values.max()\n    ymin=track_hits['y'].values.min()\n    yvar=track_hits['y'].values.var()\n    zmax=track_hits['z'].values.max()\n    zmin=track_hits['z'].values.min()\n    zvar=track_hits['z'].values.var()\n    zmean=track_hits['z'].values.mean()\n    features=np.array([svolume,nclusters,nhitspercluster,xmax,ymax,zmax,xmin,ymin,zmin,zmean,xvar,yvar,zvar])\n    return features\n\ndef get_predictions(sub,hits,model,min_length=4):  \n    \"\"\"\n    Input: dataframe sub with track id for each hit, \n           dataframe hits with hit information, \n           model=ML model to get prediction\n    Output: dataframe with predicted probability for each track\n    \"\"\"\n    preds=pd.DataFrame()\n    sub=sub.merge(hits,on='hit_id',how='left')\n    trackids_long=[]\n    trackids_short=[]\n    features=[]\n    \n    trackids=np.unique(sub['track_id']).astype(\"int64\")\n    for track_id in tqdm_notebook(trackids):        \n        track_hits=sub[sub['track_id']==track_id]\n        if len(track_hits) < min_length:\n            trackids_short.append(track_id)\n        else:\n            features.append(get_features(track_hits))\n            trackids_long.append(track_id)\n\n    probabilities_long=model.predict(np.array(features))\n    probabilities_short=np.array([0]*len(trackids_short))\n    \n    preds['quality']=np.concatenate((probabilities_long,probabilities_short))\n    preds['track_id']=np.concatenate((trackids_long,trackids_short))\n    preds['quality']=preds['quality'].fillna(1)  # assume it is a good track, if no probability can be calculated\n    return preds\n\npreds={}\npreds[0]=get_predictions(resa,hits,model)\npreds[1]=get_predictions(resa2,hits,model)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5f20e9e5d400a5041d6f1493edd135c6f426b8b7"},"cell_type":"markdown","source":"Merge both submissions, based on the probabilities of its track candidates:"},{"metadata":{"trusted":true,"_uuid":"777aad2ee971e421608ef5b2cd0a0ab4b68d611a"},"cell_type":"code","source":"def merge_with_probabilities(sub1,sub2,preds1,preds2,truth=None,length_factor=0,at_least_more=0):\n    \"\"\"\n    Input:  sub1 and sub2 are two dataframes, which assign track_ids to hits\n            preds1 and preds2 are dataframes, which assign a quality (probability to be correct) to each track_id\n            truth: if given, then calculate score of the merge\n            length_factor: merge not only by quality, but also by length\n            at_least_more: ask new quality to be at_least_more than old, to overwrite\n    \n    Output: new dataframe with updated track_ids\n    \"\"\"\n    un,inv,count = np.unique(sub1['track_id'],return_inverse=True, return_counts=True)\n    sub1['group_size']=count[inv]\n    un,inv,count = np.unique(sub2['track_id'],return_inverse=True, return_counts=True)\n    sub2['group_size']=count[inv]\n    sub1=sub1.merge(preds1,on='track_id',how='left')\n    sub2=sub2.merge(preds2,on='track_id',how='left')\n    \n    sub1['quality']=sub1['quality']+length_factor*sub1['group_size']\n    sub2['quality']=sub2['quality']+length_factor*sub2['group_size']\n    \n    sub=sub1.merge(sub2,on='hit_id',suffixes=('','_new'))\n    mm=sub.track_id.max()+1\n    sub['track_id_new']=sub['track_id_new']+mm\n    \n    sub['quality']=sub['quality']+at_least_more\n    cond=(sub['quality']>=sub['quality_new'])\n    for col in ['track_id','z0','kt']:\n        sub[col]=sub[col].where(cond,sub[col+'_new'])\n    \n    sub=sub[['hit_id','track_id','event_id','kt','z0']]\n    if not truth is None:\n        print('Score',score_event(truth,sub))\n    \n    # calculate track_ids again to make them smaller\n    un,inv,count = np.unique(sub['track_id'],return_inverse=True, return_counts=True)\n    sub['track_id']=inv\n    return sub\n\nprint('Merge submission 0 and 1 into sub01:')\nsub01=merge_with_probabilities(resa,resa2,preds[0],preds[1],truth,length_factor=0.5)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d3bb7b2f850b5bbc622c883ccbb9ca11de150999"},"cell_type":"markdown","source":"To continue, take additional submission sub2 and merge it onto sub01. Therefore, need to first calculcate probabilities of the tracks in sub01, as well as sub2. Similarly for more submissions."},{"metadata":{"trusted":true,"_uuid":"2cee39460e80a73ca8c4d86e0dcd4b8b7823dc16"},"cell_type":"code","source":"# Expanding the merged submission of the two clustering solutions gives:\nmstd_vol={7:0,8:0,9:0,12:2,13:1,14:2,16:3,17:2,18:3}\nmstd_size=[4,4,4,4,3,3,3,2,2,2,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]\nweights={'theta':0.1, 'phi':1}\nnresa2=expand_tracks(sub01,hits,5,16,5,7,mstd=8,dstd=0.00085,phik=3.3,max_dtheta=0.9*np.pi/2,mstd_vol=mstd_vol,mstd_size=mstd_size,weights=weights,nhipo=100)\nnresa2['event_id']=0\nscore = score_event_fast(truth, nresa2)\nprint(\"Your score: \", score)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f80c71d8356e923724b33362651a9705307af396"},"cell_type":"markdown","source":"\nWe achieved 0.757, which is a +0.02 improvement in the case of this notebook."},{"metadata":{"trusted":true,"_uuid":"ce764c46430500a6e5087ce2503b483bf14c957f"},"cell_type":"code","source":"","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}