{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":59093,"databundleVersionId":7469972,"sourceType":"competition"}],"dockerImageVersionId":30664,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy\nimport pandas\nimport plotly.express as px\nimport matplotlib.pyplot as plt\nimport numba","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-04-01T13:14:34.029039Z","iopub.execute_input":"2024-04-01T13:14:34.030337Z","iopub.status.idle":"2024-04-01T13:14:34.036829Z","shell.execute_reply.started":"2024-04-01T13:14:34.030287Z","shell.execute_reply":"2024-04-01T13:14:34.035015Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def read_eeg_file(file_id):\n    path = '/kaggle/input/hms-harmful-brain-activity-classification/train_eegs/' + str(file_id) + '.parquet'\n    spec = pandas.read_parquet(path)\n    spec = spec.fillna(0).values[:, 1:].T # fill NaN values with 0, transpose for (Time, Freq) -> (Freq, Time)\n    spec = spec.astype(\"float32\")\n    return spec","metadata":{"execution":{"iopub.status.busy":"2024-04-01T13:14:34.039293Z","iopub.execute_input":"2024-04-01T13:14:34.039726Z","iopub.status.idle":"2024-04-01T13:14:34.064559Z","shell.execute_reply.started":"2024-04-01T13:14:34.039683Z","shell.execute_reply":"2024-04-01T13:14:34.063345Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<center>$$ c_{(i, j)}  = \\frac{\\sum_{\\tau=t-T}^t(x_i(\\tau) - \\overline x_i)(x_j(\\tau) - \\overline x_j)}{\\sqrt{\\sum_{\\tau=t-T}^t(x_i(\\tau) - \\overline x_i^2)} \\sqrt{\\sum_{\\tau\\prime=t-T}^t(x_j(\\tau') - \\overline x_j^2)}}$$</center>","metadata":{}},{"cell_type":"code","source":"#@numba.njit\ndef computePearsonCoefficient(x, y, current_time, time_horizon, delta_t): # Python recognises two identically-named functions as different if they take different arguments. Weird, right?\n    if len(x) == 0 or len(y) == 0: # If either ticker doesn't have valid stock data associated:\n        return 0\n    x = numpy.array([(x[current_time + t + delta_t] + x[current_time + t]) / x[current_time + t] for t in range(time_horizon)])\n    y = numpy.array([(y[current_time + t + delta_t] + y[current_time + t]) / y[current_time + t] for t in range(time_horizon)])\n    covariance = numpy.sum(x - x.mean()) * numpy.sum(y - y.mean())\n    pearsons_coefficient = covariance / (numpy.var(y) * numpy.var(x))\n    return pearsons_coefficient","metadata":{"execution":{"iopub.status.busy":"2024-04-01T13:44:47.367395Z","iopub.execute_input":"2024-04-01T13:44:47.36812Z","iopub.status.idle":"2024-04-01T13:44:47.37518Z","shell.execute_reply.started":"2024-04-01T13:44:47.368078Z","shell.execute_reply":"2024-04-01T13:44:47.374153Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#computePearsonCoefficient(numpy.random.rand(10), numpy.random.rand(10), 5, 5, 1)\n# Should produce a coefficient very close to zero (two random arrays should not be well correlated).","metadata":{"execution":{"iopub.status.busy":"2024-04-01T13:45:41.779683Z","iopub.execute_input":"2024-04-01T13:45:41.7801Z","iopub.status.idle":"2024-04-01T13:45:41.84785Z","shell.execute_reply.started":"2024-04-01T13:45:41.780068Z","shell.execute_reply":"2024-04-01T13:45:41.84651Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"t = 50 # Current time\nT = 15 # Time horizon (how many seconds into the future? In this case, 15).\ndelta_t = 10\n\ntraining_sample = '1000317312.parquet'\n\nspec = pandas.read_parquet('/kaggle/input/hms-harmful-brain-activity-classification/train_spectrograms/' + training_sample)\nspec = spec.fillna(0).values[:, 1:].T # fill NaN values with 0, transpose for (Time, Freq) -> (Freq, Time)\nspec = spec.astype(\"float32\")\nnumrows, numcols = spec.shape\ncorrelationMatrix = numpy.zeros((numrows,numcols))\nfor row in range(numrows):\n    for column in range(numcols):\n        correlationMatrix[row, column] = computePearsonCoefficient(\n            spec[row],\n            spec[column],\n            t,\n            T,\n            delta_t\n    )\nplt.imshow(correlationMatrix)","metadata":{"execution":{"iopub.status.busy":"2024-04-01T13:45:01.212541Z","iopub.execute_input":"2024-04-01T13:45:01.213062Z","iopub.status.idle":"2024-04-01T13:45:15.331231Z","shell.execute_reply.started":"2024-04-01T13:45:01.213024Z","shell.execute_reply":"2024-04-01T13:45:15.330192Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.animation as animation\n\ndef computePearsonCoefficients(spec, t):\n    numrows, numcols = spec.shape\n    correlationMatrix = numpy.zeros((numrows,numcols))\n    for row in range(numrows):\n        for column in range(numcols):\n            correlationMatrix[row, column] = computePearsonCoefficient(\n                spec[row],\n                spec[column],\n                1,\n                t,\n                10\n        )\n    return correlationMatrix\n\ndef updatePearsonCoefficientPlot(frameNumber):\n    current_time = frameNumber\n    next_frame = computePearsonCoefficients(spec, current_time)\n    return next_frame","metadata":{"execution":{"iopub.status.busy":"2024-04-01T13:48:00.556422Z","iopub.execute_input":"2024-04-01T13:48:00.556807Z","iopub.status.idle":"2024-04-01T13:48:00.564167Z","shell.execute_reply.started":"2024-04-01T13:48:00.55678Z","shell.execute_reply":"2024-04-01T13:48:00.56289Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"updatePearsonCoefficientPlot(10)","metadata":{"execution":{"iopub.status.busy":"2024-04-01T13:48:04.722601Z","iopub.execute_input":"2024-04-01T13:48:04.722966Z","iopub.status.idle":"2024-04-01T13:48:17.742657Z","shell.execute_reply.started":"2024-04-01T13:48:04.722938Z","shell.execute_reply":"2024-04-01T13:48:17.741508Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots()\ncoefficient_plot = plt.imshow(updatePearsonCoefficientPlot(10))\n\ndef plot_next_timestep(iteration):\n  coefficient_plot.set_array(updatePearsonCoefficientPlot(iteration + 10))\n  return coefficient_plot\n\nani = animation.FuncAnimation(fig=fig, func=plot_next_timestep, frames=389)\nani.save(filename=\"coefficient_tests.mp4\", writer=\"ffmpeg\", fps='60')","metadata":{"execution":{"iopub.status.busy":"2024-04-01T14:27:40.300612Z","iopub.execute_input":"2024-04-01T14:27:40.301064Z","iopub.status.idle":"2024-04-01T14:29:26.108864Z","shell.execute_reply.started":"2024-04-01T14:27:40.301032Z","shell.execute_reply":"2024-04-01T14:29:26.107623Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"px.line(spec)","metadata":{"execution":{"iopub.status.busy":"2024-04-01T13:14:35.620072Z","iopub.execute_input":"2024-04-01T13:14:35.620487Z","iopub.status.idle":"2024-04-01T13:14:38.908502Z","shell.execute_reply.started":"2024-04-01T13:14:35.62045Z","shell.execute_reply":"2024-04-01T13:14:38.907114Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import plotly.graph_objects as go\n\nspectrogram = pandas.read_parquet('/kaggle/input/hms-harmful-brain-activity-classification/train_spectrograms/1000655456.parquet')\nspectro_numpy = spectrogram.to_numpy()\nfig = go.Figure()\n\nmean_spectrogram = sum(spectrogram.mean(axis=1)) / len(spectrogram.mean(axis=1))\nmean_311 = spectrogram[311]\n\nconversion_ratio = mean_311 / mean_spectrogram\n\nfig.add_trace(go.Scatter(y=reference[:250], mode='lines', name='Mean death value of<br>0-dimensional homologies'))\nfig.add_trace(go.Scatter(y=spectrogram.mean(axis=1)[:250]*conversion_ratio, mode='lines', name='Mean Spectrogram Reading'))\nfig.update_layout(\n    xaxis_title=\"Time (s)\",\n    showlegend=True\n)","metadata":{"execution":{"iopub.status.busy":"2024-04-01T13:22:24.053489Z","iopub.execute_input":"2024-04-01T13:22:24.054256Z","iopub.status.idle":"2024-04-01T13:22:25.274432Z","shell.execute_reply.started":"2024-04-01T13:22:24.054215Z","shell.execute_reply":"2024-04-01T13:22:25.272859Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#del spectrogram['time']\nspectrogram","metadata":{"execution":{"iopub.status.busy":"2024-04-01T13:14:40.393259Z","iopub.status.idle":"2024-04-01T13:14:40.393765Z","shell.execute_reply.started":"2024-04-01T13:14:40.39351Z","shell.execute_reply":"2024-04-01T13:14:40.39353Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''t = 220\nT = 20 # t > T.\ndelta_t = 10'''\n\ncorrelationMatrix = numpy.zeros((400,400))\n\nspectro_numpy = spectrogram.to_numpy()\nspectro_numpy.shape","metadata":{"execution":{"iopub.status.busy":"2024-04-01T13:14:40.395137Z","iopub.status.idle":"2024-04-01T13:14:40.395624Z","shell.execute_reply.started":"2024-04-01T13:14:40.395372Z","shell.execute_reply":"2024-04-01T13:14:40.395392Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"spectro_numpy[:, 0]","metadata":{"execution":{"iopub.status.busy":"2024-04-01T13:14:40.397013Z","iopub.status.idle":"2024-04-01T13:14:40.397486Z","shell.execute_reply.started":"2024-04-01T13:14:40.397245Z","shell.execute_reply":"2024-04-01T13:14:40.397265Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for row in range(400):\n        for column in range(400):\n            correlationMatrix[row, column] = computePearsonCoefficient(\n                spectro_numpy[:, row],\n                spectro_numpy[:, column],\n                10,\n                0,\n                15\n        )\n            \ncorrelationMatrix","metadata":{"execution":{"iopub.status.busy":"2024-04-01T13:14:40.399625Z","iopub.status.idle":"2024-04-01T13:14:40.400142Z","shell.execute_reply.started":"2024-04-01T13:14:40.399867Z","shell.execute_reply":"2024-04-01T13:14:40.399889Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"distanceMatrix = numpy.clip(correlationMatrix, -1, 1)\ndistanceMatrix = (2*(1-distanceMatrix))**0.5","metadata":{"execution":{"iopub.status.busy":"2024-04-01T13:14:40.401409Z","iopub.status.idle":"2024-04-01T13:14:40.401922Z","shell.execute_reply.started":"2024-04-01T13:14:40.401669Z","shell.execute_reply":"2024-04-01T13:14:40.40169Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from ripser import ripser\nfrom persim import plot_diagrams\n\ndata = ripser(distanceMatrix, distance_matrix=True)\nresult = data['dgms']\nh0_results = result[0]\nh0_results = h0_results[ h0_results < 2]\nplot_diagrams(result, show=True)","metadata":{"execution":{"iopub.status.busy":"2024-04-01T13:14:40.403224Z","iopub.status.idle":"2024-04-01T13:14:40.403747Z","shell.execute_reply.started":"2024-04-01T13:14:40.403461Z","shell.execute_reply":"2024-04-01T13:14:40.403481Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from ripser import ripser\nimport numpy\n\ndef computeH0Average(t):\n    correlationMatrix = numpy.zeros((400,400))\n    for row in range(400):\n        for column in range(400):\n            correlationMatrix[row, column] = computePearsonCoefficient(\n                spectro_numpy[:, row],\n                spectro_numpy[:, column],\n                10,\n                t,\n                20\n            )\n    distanceMatrix = pearsonCorrelationDataframe.to_numpy()\n    distanceMatrix = numpy.clip(distanceMatrix, -1, 1)\n    distanceMatrix = (2*(1-distanceMatrix))**0.5\n    data = ripser(distanceMatrix, distance_matrix=True)\n    result = data['dgms']\n    h0_results = result[0]\n    death_values = [value[1] for value in h0_results if value[1] <= 2]\n    return sum(death_values) / len(death_values)","metadata":{"execution":{"iopub.status.busy":"2024-04-01T13:14:40.404883Z","iopub.status.idle":"2024-04-01T13:14:40.405381Z","shell.execute_reply.started":"2024-04-01T13:14:40.405134Z","shell.execute_reply":"2024-04-01T13:14:40.405154Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"@numba.njit\ndef computePearsonCoefficient(x, y, delta_t, t, T): # Python recognises two identically-named functions as different if they take different arguments. Weird, right?\n    if len(x) == 0 or len(y) == 0: # If either ticker doesn't have valid stock data associated:\n        return 0\n    x = x[t-T:t] # Crop the array to the time horizon we're interested in. Doing this at the start is more efficient.\n    y = y[t-T:t]\n    p\n    \n    x = [ (x[value + delta_t]/x[value]) - 1 for value in range(len(x) - delta_t)] # We compute the adjusted stock price for x.\n    y = [ (y[value + delta_t]/y[value]) - 1 for value in range(len(y) - delta_t)] # And now for y.\n    mean_x = sum(x) / len(x)\n    mean_y = sum(y) / len(y)\n    topline = sum([ (x[value] - mean_x) * (y[value] - mean_y) for value in range(len(x)) ])\n    variance_i = sum([ (x[value] - mean_x)**2 for value in range(len(x)) ]) ** 0.5\n    variance_j = sum([ (y[value] - mean_y)**2 for value in range(len(y)) ]) ** 0.5\n    pearsons_coefficient = topline / (variance_i * variance_j)\n    return pearsons_coefficient\n\nmean_h0_value = []\nfor t in range(399 - 15 - 1):\n    mean_h0_value.append(computeH0Average(t))","metadata":{"execution":{"iopub.status.busy":"2024-04-01T13:14:40.406764Z","iopub.status.idle":"2024-04-01T13:14:40.407273Z","shell.execute_reply.started":"2024-04-01T13:14:40.407013Z","shell.execute_reply":"2024-04-01T13:14:40.407034Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mean_h0_value","metadata":{"execution":{"iopub.status.busy":"2024-04-01T13:14:40.408693Z","iopub.status.idle":"2024-04-01T13:14:40.409201Z","shell.execute_reply.started":"2024-04-01T13:14:40.408928Z","shell.execute_reply":"2024-04-01T13:14:40.408949Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"For how many positive integers $m$ does the equation \\[\\vert \\vert x-1 \\vert -2 \\vert=\\frac{m}{100}\\] have $4$ integer solutions?","metadata":{"execution":{"iopub.status.busy":"2024-04-02T14:42:14.729116Z","iopub.execute_input":"2024-04-02T14:42:14.729665Z","iopub.status.idle":"2024-04-02T14:42:14.759354Z","shell.execute_reply.started":"2024-04-02T14:42:14.729634Z","shell.execute_reply":"2024-04-02T14:42:14.758645Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"For how many positive integers $m$ does the equation \\[$\\vert \\vert x-1 \\vert -2 \\vert=\\frac{m}{100}$\\] have $4$ integer solutions?","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}