{"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":"!pip install --quiet tables","metadata":{"execution":{"iopub.status.busy":"2022-08-25T08:11:48.773832Z","iopub.execute_input":"2022-08-25T08:11:48.774278Z","iopub.status.idle":"2022-08-25T08:12:03.726787Z","shell.execute_reply.started":"2022-08-25T08:11:48.774242Z","shell.execute_reply":"2022-08-25T08:12:03.725660Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nimport scipy.sparse as sps\nfrom tqdm import tqdm as tqdm\nimport gc","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-08-25T08:12:03.729226Z","iopub.execute_input":"2022-08-25T08:12:03.729596Z","iopub.status.idle":"2022-08-25T08:12:03.736801Z","shell.execute_reply.started":"2022-08-25T08:12:03.729557Z","shell.execute_reply":"2022-08-25T08:12:03.735593Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DATA_DIR = \"/kaggle/input/open-problems-multimodal/\"\n\nSUBMISSON = os.path.join(DATA_DIR,\"sample_submission.csv\")\n\nEVALUATION_IDS = os.path.join(DATA_DIR,\"evaluation_ids.csv\")\n\nFP_CELL_METADATA = os.path.join(DATA_DIR,\"metadata.csv\")\n\nFP_CITE_TRAIN_INPUTS = os.path.join(DATA_DIR,\"train_cite_inputs.h5\")\nFP_CITE_TRAIN_TARGETS = os.path.join(DATA_DIR,\"train_cite_targets.h5\")\nFP_CITE_TEST_INPUTS = os.path.join(DATA_DIR,\"test_cite_inputs.h5\")\n\nFP_MULTIOME_TRAIN_INPUTS = os.path.join(DATA_DIR,\"train_multi_inputs.h5\")\nFP_MULTIOME_TRAIN_TARGETS = os.path.join(DATA_DIR,\"train_multi_targets.h5\")\nFP_MULTIOME_TEST_INPUTS = os.path.join(DATA_DIR,\"test_multi_inputs.h5\")\n\nFP_SUBMISSION = os.path.join(DATA_DIR,\"sample_submission.csv\")\nFP_EVALUATION_IDS = os.path.join(DATA_DIR,\"evaluation_ids.csv\")","metadata":{"execution":{"iopub.status.busy":"2022-08-25T08:18:19.612711Z","iopub.execute_input":"2022-08-25T08:18:19.613155Z","iopub.status.idle":"2022-08-25T08:18:19.621773Z","shell.execute_reply.started":"2022-08-25T08:18:19.613122Z","shell.execute_reply":"2022-08-25T08:18:19.620735Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Multiome Dataset","metadata":{}},{"cell_type":"markdown","source":"According to https://www.kaggle.com/code/ambrosm/msci-multiome-quickstart, Multiome dataset is way to large to fit into the 16GB memory available on Kaggle. In fact:\n- train inputs: 105942 * 228942 float32 values (97 GByte)\n- train targets: 105942 * 23418 float32 values (10 GByte)\n- test inputs: 55935 * 228942 float32 values (13 GByte)","metadata":{}},{"cell_type":"markdown","source":"## Problem\nAs we can see from the competition datasets, Multiome data are instrinsically sparse. To prove this statement, we can measure the sparsity rate of the Train-Multi-Inputs dataset. As described above, the entire dataset cannot be load in memory, thus we limit our study to the first 5000 rows","metadata":{}},{"cell_type":"code","source":"df = pd.read_hdf(FP_MULTIOME_TRAIN_INPUTS, start=0, stop=5000)","metadata":{"execution":{"iopub.status.busy":"2022-08-25T08:18:35.442881Z","iopub.execute_input":"2022-08-25T08:18:35.443628Z","iopub.status.idle":"2022-08-25T08:19:01.739900Z","shell.execute_reply.started":"2022-08-25T08:18:35.443585Z","shell.execute_reply":"2022-08-25T08:19:01.738630Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.info(memory_usage='deep')","metadata":{"execution":{"iopub.status.busy":"2022-08-25T08:19:32.586432Z","iopub.execute_input":"2022-08-25T08:19:32.586875Z","iopub.status.idle":"2022-08-25T08:19:49.001501Z","shell.execute_reply.started":"2022-08-25T08:19:32.586839Z","shell.execute_reply":"2022-08-25T08:19:49.000044Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Count Non-Zero Values in Each column","metadata":{"execution":{"iopub.status.busy":"2022-08-25T08:21:15.848095Z","iopub.execute_input":"2022-08-25T08:21:15.848489Z","iopub.status.idle":"2022-08-25T08:21:15.856363Z","shell.execute_reply.started":"2022-08-25T08:21:15.848459Z","shell.execute_reply":"2022-08-25T08:21:15.854241Z"}}},{"cell_type":"code","source":"nnz = df.astype(bool).sum()\nnnz.sort_values()","metadata":{"execution":{"iopub.status.busy":"2022-08-25T08:28:41.066916Z","iopub.execute_input":"2022-08-25T08:28:41.067350Z","iopub.status.idle":"2022-08-25T08:28:43.833581Z","shell.execute_reply.started":"2022-08-25T08:28:41.067316Z","shell.execute_reply":"2022-08-25T08:28:43.832359Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To measure the total sparsity of the DataFrame, we can extract the fraction of NNZ values over the total number of values","metadata":{}},{"cell_type":"code","source":"total_nnz = nnz.sum()\ntotal_values = df.shape[0] * df.shape[1]\ntotal_nnz / total_values","metadata":{"execution":{"iopub.status.busy":"2022-08-25T08:25:22.016400Z","iopub.execute_input":"2022-08-25T08:25:22.016777Z","iopub.status.idle":"2022-08-25T08:25:22.026014Z","shell.execute_reply.started":"2022-08-25T08:25:22.016746Z","shell.execute_reply":"2022-08-25T08:25:22.024993Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see, the dataset is extreamly sparse, since the Number of Non-Zero values correspond to just `2%` of the entire dataset loaded. It is reasonable to state that the same behaviour holds in the rest of the dataset. \nWe are able to tackle this waste of memory by adopting a different data structure","metadata":{}},{"cell_type":"code","source":"del df, nnz, total_nnz, total_values","metadata":{"execution":{"iopub.status.busy":"2022-08-25T08:33:49.297234Z","iopub.execute_input":"2022-08-25T08:33:49.298444Z","iopub.status.idle":"2022-08-25T08:33:49.308171Z","shell.execute_reply.started":"2022-08-25T08:33:49.298389Z","shell.execute_reply":"2022-08-25T08:33:49.307085Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-08-25T08:33:54.817122Z","iopub.execute_input":"2022-08-25T08:33:54.817520Z","iopub.status.idle":"2022-08-25T08:33:56.962226Z","shell.execute_reply.started":"2022-08-25T08:33:54.817489Z","shell.execute_reply":"2022-08-25T08:33:56.960795Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Memory Optimization with Sparse Matrices\nGiven the intrinsic sparse nature of the data in Multiome datasets, we can leverage on Sparse Matrices to optimize the space required to load data in memory. In particular, we can use [Compressed Sparse Row](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csr_matrix.html) matrices to reduce considerably the memory used.\n\nCSR Matrix are built upon three different one-dimensional arrays:\n- Data Array: Shape: (Number Non-Zero values). It contains non-zero values that corresponds to our data.\n- Indices Array: Shape: (Number Non-Zero values). It contains the column indices\n- Indptr Array: Shape: (Number of Rows + 1). It represents the extent of each row with respect to the other two (data/indices) arrays. To access data of a particular row *i* in the matrix, we can slice the Data Array with Indptr Array as follows: `data[indptr[i]:indptr[i+1]]`. Same for Indices Array","metadata":{}},{"cell_type":"markdown","source":"Since we are not able to load the entire Train-Multi-Inputs dataset in memory, we are going to manually build the three arrays by loading chunk of data at a time.","metadata":{}},{"cell_type":"markdown","source":"### Utility functions","metadata":{}},{"cell_type":"markdown","source":"To speed up the computation, we compute the indptr array by exploiting Cython. In this way, we can halve the time required to compress the huge array of row indices to extract the indptr array","metadata":{}},{"cell_type":"code","source":"%load_ext Cython","metadata":{"execution":{"iopub.status.busy":"2022-08-25T08:49:45.197199Z","iopub.execute_input":"2022-08-25T08:49:45.197599Z","iopub.status.idle":"2022-08-25T08:49:45.859437Z","shell.execute_reply.started":"2022-08-25T08:49:45.197569Z","shell.execute_reply":"2022-08-25T08:49:45.858288Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%cython\n\nimport cython\ncimport cython\ncimport numpy as np\nimport numpy as np\nfrom tqdm import tqdm, trange\n\nctypedef np.int64_t INT64_t\n\n@cython.boundscheck(False)\n@cython.wraparound(False)\ncpdef np.ndarray[INT64_t, ndim=1] create_indptr(INT64_t[:] row_indices, int start_pos, int nrows):\n    cdef int shape = row_indices.shape[0]\n    res = np.zeros(nrows, dtype=np.int64)\n    cdef INT64_t[:] res_view = res\n    \n    cdef int i\n    cdef int curr_row = 0\n    cdef int prev = row_indices[0]\n    \n    for i in range(shape):\n        if row_indices[i] != prev:\n            curr_row += 1\n            res_view[curr_row] = i\n            prev = row_indices[i]\n    # res_view[curr_row + 1] = shape\n    return res + start_pos\n","metadata":{"execution":{"iopub.status.busy":"2022-08-25T08:49:45.861858Z","iopub.execute_input":"2022-08-25T08:49:45.862360Z","iopub.status.idle":"2022-08-25T08:49:53.831971Z","shell.execute_reply.started":"2022-08-25T08:49:45.862314Z","shell.execute_reply":"2022-08-25T08:49:53.830592Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def create_csr_arrays(h5_file_path):\n    def check_size(xs, ys, datas):\n        return (xs.nbytes + ys.nbytes + datas.nbytes) * 1e-9\n\n    print(f\"\\n\\nProcessing File {h5_file_path}\")\n    pbar = tqdm()\n\n    # Initialize Variables\n    chunksize = 1000 # Keep it low\n    loaded_rows = chunksize\n    start = 0\n    start_pos = 0\n    file_pointer = 0\n\n    # Initialize CSR arrays\n    indptr = np.array([], dtype=np.int64)\n    indices = np.array([], dtype=np.int32)\n    data_s = np.array([], dtype=np.float32)\n    \n    prefix_filename = h5_file_path.split('/')[-1].replace('.h5', '')\n\n    while chunksize == loaded_rows:\n\n        # Check current size: if the total sum of sizes are > 7GB, then save three arrays and re-initialize them\n        size_gb = check_size(indptr, indices, data_s)\n        if size_gb > 7.0:\n            pbar.set_description(f\"Total size is {size_gb}. Saving ..\")\n            np.save(f\"{prefix_filename}_indptr_{file_pointer}.npy\", indptr)\n            np.save(f\"{prefix_filename}_indices_{file_pointer}.npy\", indices)\n            np.save(f\"{prefix_filename}_data_{file_pointer}.npy\", data_s)\n            # Re-initialize\n            indptr = np.array([], dtype=np.int64)\n            indices = np.array([], dtype=np.int32)\n            data_s = np.array([], dtype=np.float32)\n            # Increment pointer\n            file_pointer += 1\n\n        pbar.set_description(\"Reading .h5 chunk\")\n        df = pd.read_hdf(h5_file_path, start=start, stop=start+chunksize)\n        pbar.set_description(\"Extracting non-zero values\")\n        x_coords, y_coords = df.values.nonzero()\n        tmp_data = df.values[df.values != 0.0]\n\n        loaded_rows = df.shape[0]\n\n        # Convert types\n        y_coords = y_coords.astype(np.int32, copy=False)\n        tmp_data = tmp_data.astype(np.float32, copy=False)\n\n        # Compress x_coords\n        pbar.set_description(\"Compressing rows values\")\n        x_coords = create_indptr(x_coords, start_pos=start_pos, nrows=loaded_rows)\n\n        gc.collect()\n\n        # Update variables\n        pbar.set_description(\"Update variables\")\n        start_pos += y_coords.shape[0]\n        start += chunksize\n        # Append data at the end of each array\n        indptr = np.hstack((indptr, x_coords))\n        indices = np.hstack((indices, y_coords))\n        data_s = np.hstack((data_s, tmp_data))\n\n        pbar.update(loaded_rows)\n\n    print('Done. Save last files')\n    np.save(f\"{prefix_filename}_indptr_{file_pointer}.npy\", indptr)\n    np.save(f\"{prefix_filename}_indices_{file_pointer}.npy\", indices)\n    np.save(f\"{prefix_filename}_data_{file_pointer}.npy\", data_s)\n    \n    del indptr, indices, data_s\n","metadata":{"execution":{"iopub.status.busy":"2022-08-25T09:00:11.116488Z","iopub.execute_input":"2022-08-25T09:00:11.117398Z","iopub.status.idle":"2022-08-25T09:00:11.134176Z","shell.execute_reply.started":"2022-08-25T09:00:11.117353Z","shell.execute_reply":"2022-08-25T09:00:11.133011Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# create_csr_arrays(FP_MULTIOME_TRAIN_INPUTS) # This will create three different arrays","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The previous command will create and save three different array in .npy format:\n- train_multi_inputs_indptr_0.npy\n- train_multi_inputs_indices_0.npy\n- train_multi_inputs_data_0.npy","metadata":{}},{"cell_type":"code","source":"# indptr = np.load('train_multi_inputs_indptr_0.npy')\n# indices = np.load('train_multi_inputs_indices_0.npy')\n# data = np.load('train_multi_inputs_data_0.npy')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Since indptr array has shape (Number of Rows) instead of (Number of Rows + 1), we can add the last element to the array, which corresponds to the length of indices or data arrays. ","metadata":{}},{"cell_type":"code","source":"# indptr = np.append(indptr, indptr[-1] + indices[indptr[-1]:].shape)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Eventually, we can build out csr_matrix as follows:","metadata":{}},{"cell_type":"code","source":"N_ROWS = 105942\nN_COLS = 228942\n# csr_matrix = sps.csr_matrix((data, indices, indptr), shape=(N_ROWS, N_COLS))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# sps.save_npz('train_multiome_input_sparse.npz', csr_matrix)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# del csr_matrix, indices, indptr, data","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can repeat the same process for the other Multiome Datasets, namely `train_multi_targets.h5` and `test_multi_inputs.h5` to obtain the corresponding Compressed Sparse Row matrices.\nI wrapped up these CSR matrices in the following Kaggle Dataset: https://www.kaggle.com/datasets/sbunzini/open-problems-msci-multiome-sparse-matrices","metadata":{}},{"cell_type":"markdown","source":"# Compression Rate","metadata":{}},{"cell_type":"code","source":"train_input = sps.load_npz('../input/open-problems-msci-multiome-sparse-matrices/train_multiome_input_sparse.npz')","metadata":{"execution":{"iopub.status.busy":"2022-08-25T07:33:32.601286Z","iopub.execute_input":"2022-08-25T07:33:32.601659Z","iopub.status.idle":"2022-08-25T07:34:41.510065Z","shell.execute_reply.started":"2022-08-25T07:33:32.601628Z","shell.execute_reply":"2022-08-25T07:34:41.509089Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_size(sparse_m):\n    size_gb = (sparse_m.indices.nbytes + sparse_m.indptr.nbytes + sparse_m.data.nbytes) * 1e-9\n    return f\"Size: {size_gb} GB\"","metadata":{"execution":{"iopub.status.busy":"2022-08-25T07:34:41.511972Z","iopub.execute_input":"2022-08-25T07:34:41.512385Z","iopub.status.idle":"2022-08-25T07:34:41.518702Z","shell.execute_reply.started":"2022-08-25T07:34:41.512353Z","shell.execute_reply":"2022-08-25T07:34:41.517618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"get_size(train_input)","metadata":{"execution":{"iopub.status.busy":"2022-08-25T08:10:14.005071Z","iopub.execute_input":"2022-08-25T08:10:14.005531Z","iopub.status.idle":"2022-08-25T08:10:14.012902Z","shell.execute_reply.started":"2022-08-25T08:10:14.005495Z","shell.execute_reply":"2022-08-25T08:10:14.012137Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Memory Usage: `4.85883614 GB`","metadata":{}},{"cell_type":"code","source":"# Percentage of Reduction\n(1.0 - (4.85883614 / 97)) * 100","metadata":{"execution":{"iopub.status.busy":"2022-08-24T19:36:53.033261Z","iopub.execute_input":"2022-08-24T19:36:53.033641Z","iopub.status.idle":"2022-08-24T19:36:53.042925Z","shell.execute_reply.started":"2022-08-24T19:36:53.033614Z","shell.execute_reply":"2022-08-24T19:36:53.040937Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Reduced Memory Usage: `94.99%`\n","metadata":{}},{"cell_type":"markdown","source":"Same memory usage reduction can be applied to the other Multiome files (train_targets and test_inputs). Lots of state-of-the-art models can accept a sparse matrix as input for training, thus avoiding painful and slow iterators and speeding up the computation","metadata":{}},{"cell_type":"markdown","source":"# !! Update !!\nThe memory usage can be further shrinked by using float16 to represent data values and int16 to represent indices of columns. A new version of the dataset will be available with this kind of optimization which will allow to achieve a **97%** of compression","metadata":{}}]}