{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"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 os\nimport math\nimport matplotlib.pyplot as plt\nimport plotly.graph_objects as go\nfrom plotly.subplots import make_subplots\nimport pandas as pd","metadata":{"_uuid":"6f234392-cd2a-46ce-a660-abb53ed0e9c9","_cell_guid":"51902ff4-4153-46ef-86f8-707ac9dcd4f3","collapsed":false,"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2023-02-04T22:27:30.039233Z","iopub.execute_input":"2023-02-04T22:27:30.039804Z","iopub.status.idle":"2023-02-04T22:27:30.046285Z","shell.execute_reply.started":"2023-02-04T22:27:30.039755Z","shell.execute_reply":"2023-02-04T22:27:30.045073Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#This seems to help when graphics stop displaying \nfrom plotly.offline import init_notebook_mode, iplot\ninit_notebook_mode(connected=True)","metadata":{"execution":{"iopub.status.busy":"2023-02-04T22:29:54.460510Z","iopub.execute_input":"2023-02-04T22:29:54.460986Z","iopub.status.idle":"2023-02-04T22:29:54.468511Z","shell.execute_reply.started":"2023-02-04T22:29:54.460951Z","shell.execute_reply":"2023-02-04T22:29:54.467224Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We divide the surface of the sphere into equal sized patches and then count arrivals within each patch.  This is the correct way to determine if there is a preferential arrival direction(s).  If you just bin by zenith angle, it will seem like a large number come from +/- pi/2 because that zenith angle covers a larger spherical angle (for a constant azimuth spread).","metadata":{}},{"cell_type":"code","source":"#CONSTANTS\nDATA_DIR='../input/train-meta-parquet'","metadata":{"execution":{"iopub.status.busy":"2023-02-04T22:27:30.049024Z","iopub.execute_input":"2023-02-04T22:27:30.049408Z","iopub.status.idle":"2023-02-04T22:27:30.067936Z","shell.execute_reply.started":"2023-02-04T22:27:30.049371Z","shell.execute_reply":"2023-02-04T22:27:30.066814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Note that this notebook uses the train-meta-parquet dataset which breaks apart the train_meta parquet file so as to speed load time and reduce memory usage","metadata":{}},{"cell_type":"markdown","source":"**CREATE SPHERICAL BINS OF EQUAL AREA**","metadata":{}},{"cell_type":"code","source":"def compute_zenith_steps(rings=[]):\n    \"\"\"rings are number of pieces each spherical segment is cut into, len(rings) is number of\n    segments\"\"\"\n    zens = []   #list of zeniths that make up the upper edge of each ring.  Last one s/b pi/2\n    total_area = 2*math.pi*1  # half of a sphere of radius 1\n    total_segs = sum(rings)\n    seg_area = total_area/total_segs\n    prior_zenith=0.0\n    for p in rings:\n        #area of all previous rings plus this one is:\n        cum_area = (1-math.cos(prior_zenith))*2*math.pi + p * seg_area\n        #now we know the height of this ring (interesting but not needed)\n        h = cum_area/(2*math.pi) - (1-math.cos(prior_zenith))\n        #The z height of this ring, more useful\n        z = cum_area/(2*math.pi)\n        #and the zenith angle should be\n        zen = math.acos(1-z)\n        zens.append(zen)\n        prior_zenith = zen\n    zens_degrees = [ z*180/math.pi for z in zens]\n    #print(f'got zens: {zens_degrees}')\n    return zens\n\ndef polar_to_cartesian(azimuth, zenith):\n    # see: https://stackoverflow.com/a/10868220/4521646\n    # now works with arrays\n    x = np.cos(azimuth) * np.sin(zenith)\n    y = np.sin(azimuth) * np.sin(zenith)\n    z = np.cos(zenith)\n    return x, y, z\n\ndef generate_table(zens, rings):\n    \"\"\"create np array az0 az1 z0 z1 0 for each bin\n    The last zero will be a placeholder for putting counts\"\"\"\n    ret=[]\n    for half in range(2):\n        for p, z0, z1 in zip(rings, [0]+zens[:-1], zens):\n            if half==1:\n                z0 = math.pi - z0\n                z1 = math.pi - z1\n                z0,z1 = z1,z0\n            for i in range(p):\n                az0 = i*2*math.pi/p\n                az1 = (i+1)*2*math.pi/p\n                ret.append([az0,az1,z0,z1,0])\n    return np.array(ret)\n\ndef generate_lines(zens, rings, halves=1):\n    \"\"\"generate the x,y,z coordinates of the region boundaries for list of zeniths and count of rings\n    in each spherical shell.  Halves is 1 or 2 \"\"\"\n    x=[]  # each will be a tuple of x0,x1 for that line\n    y=[]\n    z=[]\n    for idx in range(halves):\n        #drawing the top half?\n        for p, z0, z1 in zip(rings, [0]+zens[:-1], zens):\n            if idx==1:\n                z0 = math.pi - z0\n                z1 = math.pi - z1\n            #draw the z1 ring\n            N=60\n            for i in range(N):\n                p1x,p1y,p1z=polar_to_cartesian(i*2*math.pi/N, z1)\n                p2x,p2y,p2z=polar_to_cartesian((i+1)*2*math.pi/N, z1)\n                x.extend((p1x,p2x,None))\n                y.extend((p1y,p2y,None))\n                z.extend((p1z,p2z,None))\n            #Draw the cross pieces\n            for i in range(p):\n                p1x,p1y,p1z=polar_to_cartesian(i*2*math.pi/p, z0)\n                p2x,p2y,p2z=polar_to_cartesian(i*2*math.pi/p, z1)\n                x.extend((p1x,p2x,None))\n                y.extend((p1y,p2y,None))\n                z.extend((p1z,p2z,None))\n    return np.array(x),np.array(y),np.array(z)","metadata":{"execution":{"iopub.status.busy":"2023-02-04T22:27:30.069819Z","iopub.execute_input":"2023-02-04T22:27:30.070524Z","iopub.status.idle":"2023-02-04T22:27:30.091168Z","shell.execute_reply.started":"2023-02-04T22:27:30.070489Z","shell.execute_reply":"2023-02-04T22:27:30.089906Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def show_3d_lines(x,y,z, op=.5):\n    \"\"\"plot xyz lines in 3d\"\"\"\n    fig = make_subplots(\n        rows=1, specs=[[{'type': 'scene'}]],\n        subplot_titles=['arrivals']\n    )\n    fig.add_trace(go.Scatter3d(\n        x=x, y=y, z=z, \n        mode='lines', line=dict(color='blue',width=3), opacity=op\n    ), row=1, col=1)\n    fig.update_layout(\n        height=800, width=800, showlegend=False,\n        title_text='Segments',\n    )\n    return fig","metadata":{"execution":{"iopub.status.busy":"2023-02-04T22:27:30.092826Z","iopub.execute_input":"2023-02-04T22:27:30.093223Z","iopub.status.idle":"2023-02-04T22:27:30.110407Z","shell.execute_reply.started":"2023-02-04T22:27:30.093187Z","shell.execute_reply":"2023-02-04T22:27:30.109338Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rings = [4,12,24,36,44,48,52,56,60,64,64,64,64,64,64,64,64,64,64]\nzens = compute_zenith_steps(rings)\nprint(f'      delta-zenith   num segments   segment degrees')\nfor i, (z0,z1,p) in enumerate(zip([0]+zens[:-1],zens, rings)):\n    delta = (z1-z0)*180/math.pi\n    print(f'{i:3}     {delta:10.2f}     {p:10}       {360/p:8.2f}')\nx_lines,y_lines,z_lines = generate_lines(zens, rings)\n#table has rows az0,az1,zen0,zen1 ranges for each bin\ntable = generate_table(zens, rings)\nprint(x_lines.shape, y_lines.shape, z_lines.shape)\n","metadata":{"execution":{"iopub.status.busy":"2023-02-04T22:31:52.507009Z","iopub.execute_input":"2023-02-04T22:31:52.507841Z","iopub.status.idle":"2023-02-04T22:31:52.550691Z","shell.execute_reply.started":"2023-02-04T22:31:52.507797Z","shell.execute_reply":"2023-02-04T22:31:52.549467Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here are the bin patches plotted in 3D.  Each has equal spherical area. This just shows the top half; we will use the entire sphere later.","metadata":{}},{"cell_type":"code","source":"fig = show_3d_lines(x_lines,y_lines,z_lines)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-02-04T22:31:56.114571Z","iopub.execute_input":"2023-02-04T22:31:56.115508Z","iopub.status.idle":"2023-02-04T22:31:56.452731Z","shell.execute_reply.started":"2023-02-04T22:31:56.115450Z","shell.execute_reply":"2023-02-04T22:31:56.451274Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bins = np.zeros_like(table[:,0])\nfor batch_id in range(1,11):\n    df = pd.read_parquet( os.path.join(DATA_DIR, f'train_meta_{batch_id}.parquet') )\n    azimuth = df['azimuth'].values\n    zenith = df['zenith'].values\n    x,y,z = polar_to_cartesian(azimuth, zenith)\n    #xyz = np.stack((x,y,z), axis=1)\n    print(f'xyz shape {x.shape}')\n    #count into buckets\n    for i, row in enumerate(table):\n        #print(row)\n        a1 = np.logical_and(row[0] <= azimuth, azimuth < row[1])\n        a2 = np.logical_and(row[2] <= zenith, zenith < row[3])\n        a3 = np.logical_and(a1,a2)\n        bins[i]+=np.sum(a3)","metadata":{"execution":{"iopub.status.busy":"2023-02-04T22:27:30.377308Z","iopub.execute_input":"2023-02-04T22:27:30.378143Z","iopub.status.idle":"2023-02-04T22:27:42.957652Z","shell.execute_reply.started":"2023-02-04T22:27:30.378088Z","shell.execute_reply":"2023-02-04T22:27:42.956607Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def add_3d_markers(fig,x,y,z,count):\n    \"\"\"Add markers at z,y,z with count \"\"\"\n    mn = np.mean(count)\n    sz = 10 + 10*(count-mn)/mn\n    fig.add_trace(go.Scatter3d(\n        x=x,y=y,z=z,\n        mode='markers', \n        marker=dict(size=sz, color=count, colorscale='gray'),\n        opacity=0.8\n    ), row=1, col=1)","metadata":{"execution":{"iopub.status.busy":"2023-02-04T22:35:09.363373Z","iopub.execute_input":"2023-02-04T22:35:09.363873Z","iopub.status.idle":"2023-02-04T22:35:09.372787Z","shell.execute_reply.started":"2023-02-04T22:35:09.363834Z","shell.execute_reply":"2023-02-04T22:35:09.371419Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"x_lines,y_lines,z_lines = generate_lines(zens, rings, halves=2)\naz_center = np.mean(table[:,0:2], axis=1)\nzen_center = np.mean(table[:,2:4], axis=1)\n#az_center=table[:,0]\n#zen_center=table[:,2]\nx_bins,y_bins,z_bins = polar_to_cartesian(az_center,zen_center)\nfig = show_3d_lines(x_lines,y_lines,z_lines,op=.2)\n#bins = bins/np.max(bins)\nadd_3d_markers(fig, x_bins, y_bins, z_bins, bins)\nfig.show()\n#Show the counts as size of marker","metadata":{"execution":{"iopub.status.busy":"2023-02-04T22:35:11.170358Z","iopub.execute_input":"2023-02-04T22:35:11.170881Z","iopub.status.idle":"2023-02-04T22:35:11.787153Z","shell.execute_reply.started":"2023-02-04T22:35:11.170840Z","shell.execute_reply":"2023-02-04T22:35:11.785958Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(az_center.shape)\nbig = np.stack((az_center,zen_center,bins ), axis=1)\nprint(big.shape)\ndf = pd.DataFrame(big, columns=['az','zen','count'])\n#print(df)\ndf['zena']=(df.zen*1000).astype(np.int32)\n#grp = df.groupby('zena')['zena'].sum()\ngrp = df.groupby('zena')['count'].agg(['sum','count'])\ncount_per_bin = grp.loc[:,'sum']/grp.loc[:,'count']\nplt.plot(grp.index/1000/math.pi*180, count_per_bin)\nplt.xlabel('zenith (degrees)')\nplt.ylabel('count')","metadata":{"execution":{"iopub.status.busy":"2023-02-04T22:27:43.407898Z","iopub.execute_input":"2023-02-04T22:27:43.408325Z","iopub.status.idle":"2023-02-04T22:27:43.656868Z","shell.execute_reply.started":"2023-02-04T22:27:43.408289Z","shell.execute_reply":"2023-02-04T22:27:43.655637Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"For more accuracy, we can run this with a larger number of batches, but you can see the larger number of events that come from smaller zenith.  The larger zenith angles (>120?) are blocked by the Earth partially.","metadata":{}}]}