{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.14","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":84969,"databundleVersionId":10033515,"sourceType":"competition"},{"sourceId":9867543,"sourceType":"datasetVersion","datasetId":6040935},{"sourceId":9876108,"sourceType":"datasetVersion","datasetId":6040928}],"dockerImageVersionId":30787,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"#!pip download connected-components-3d\n#!pip download zarr\n#!pip download numcodecs\n\ntry:\n    import zarr\nexcept: \n    !cp -r '/kaggle/input/hengck-czii-cryo-et-01/wheel_file' '/kaggle/working/'\n    !pip install /kaggle/working/wheel_file/asciitree-0.3.3/asciitree-0.3.3\n    !pip install --no-index --find-links=/kaggle/working/wheel_file zarr\n    !pip install --no-index --find-links=/kaggle/working/wheel_file connected-components-3d\n\nprint('PIP INSTALL OK!!!')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-11T18:53:01.594697Z","iopub.execute_input":"2024-11-11T18:53:01.595576Z","iopub.status.idle":"2024-11-11T18:53:01.606673Z","shell.execute_reply.started":"2024-11-11T18:53:01.595536Z","shell.execute_reply":"2024-11-11T18:53:01.605794Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\nfrom datetime import datetime\nimport pytz\nprint('LOGGING TIME OF START:',  datetime.strftime(datetime.now(pytz.timezone('Asia/Singapore')), \"%Y-%m-%d %H:%M:%S\"))\n\nimport sys\nsys.path.append('/kaggle/input/hengck-czii-cryo-et-01')\n\nfrom czii_helper import *\nfrom dataset import *\nfrom model2 import *\nimport numpy as np\nfrom scipy.optimize import linear_sum_assignment\nimport glob\nimport cc3d\nimport cv2\n\n\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d import Axes3D\n\n\nprint('IMPORT OK!!!')","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2024-11-11T18:53:01.608299Z","iopub.execute_input":"2024-11-11T18:53:01.608596Z","iopub.status.idle":"2024-11-11T18:53:01.624112Z","shell.execute_reply.started":"2024-11-11T18:53:01.608566Z","shell.execute_reply":"2024-11-11T18:53:01.623292Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"DATA_KAGGLE_DIR = '/kaggle/input/czii-cryo-et-object-identification'\n\nMODE='submit'\n\nif MODE=='local':\n    valid_dir =f'{DATA_KAGGLE_DIR}/train'\n    valid_id = ['TS_73_6', ]\n    \nif MODE=='submit':\n    valid_dir =f'{DATA_KAGGLE_DIR}/test' \n    valid_id = glob.glob(f'{valid_dir}/static/ExperimentRuns/*')\n    valid_id = [f.split('/')[-1] for f in valid_id]\n    \nprint('valid_id:',len(valid_id), valid_id)\n\ncfg = dotdict(\n    arch ='resnet34d',\n    checkpoint=\\\n    '/kaggle/input/hengck-czii-cryo-et-weights-01/resnet34d-00010164-extern.pth',\n   # '/kaggle/input/hengck-czii-cryo-et-weights-01/resnet34d-aug-noise-00003956.pth',\n       #  '/kaggle/input/hengck-czii-cryo-et-weights-01/resnet34d-00002300.pth',\n       # '/kaggle/input/hengck-czii-cryo-et-weights-01/00003531.pth',\n    threshold={ \n        'apo-ferritin': 0.05,\n        'beta-amylase': 0.05,\n        'beta-galactosidase': 0.05,\n        'ribosome': 0.05,\n        'thyroglobulin': 0.05,\n        'virus-like-particle': 0.05,\n    },\n)\n\nprint('MODE:', MODE)\nprint('SETTING OK!!!')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-11T18:53:01.625171Z","iopub.execute_input":"2024-11-11T18:53:01.625565Z","iopub.status.idle":"2024-11-11T18:53:01.637146Z","shell.execute_reply.started":"2024-11-11T18:53:01.625501Z","shell.execute_reply":"2024-11-11T18:53:01.636332Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"net = Net(pretrained=False, cfg=cfg)\nstate_dict = torch.load(cfg.checkpoint, map_location=lambda storage, loc: storage)['state_dict']\nprint(net.load_state_dict(state_dict, strict=False))\n\nprint(net.arch)\nprint('MODEL OK!!!')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-11T18:53:01.639240Z","iopub.execute_input":"2024-11-11T18:53:01.639871Z","iopub.status.idle":"2024-11-11T18:53:02.216306Z","shell.execute_reply.started":"2024-11-11T18:53:01.639825Z","shell.execute_reply":"2024-11-11T18:53:02.215128Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def probability_to_location(probability,cfg):\n    _,D,H,W = probability.shape\n\n    location={}\n    for p in PARTICLE:\n        p = dotdict(p)\n        l = p.label\n\n        cc, P = cc3d.connected_components(probability[l]>cfg.threshold[p.name], return_N=True)\n        stats = cc3d.statistics(cc)\n        zyx=stats['centroids'][1:]*10.012444\n        xyz = np.ascontiguousarray(zyx[:,::-1]) \n        location[p.name]=xyz\n        '''\n            j=1\n            z,y,x = np.where(cc==j)\n            z=z.mean()\n            y=y.mean()\n            x=x.mean()\n            print([x,y,z])\n        '''\n    return location\n\ndef location_to_df(location):\n    location_df = []\n    for p in PARTICLE:\n        p = dotdict(p)\n        xyz = location[p.name]\n        if len(xyz)>0:\n            df = pd.DataFrame(data=xyz, columns=['x','y','z'])\n            #df.loc[:,'particle_type']= p.name\n            df.insert(loc=0, column='particle_type', value=p.name)\n            location_df.append(df)\n    location_df = pd.concat(location_df)\n    return location_df\n\n\n#start here !!!! ====================================================================================\ndef run_submit():\n    \n    net.cuda()\n    net.eval()\n    net.output_type = ['infer']\n\n    submit_df = []\n    start_timer = timer()\n    for i,id in enumerate(valid_id):\n        print(i, id, '---------------')\n        volume = read_one_data(id, static_dir=f'{valid_dir}/static/ExperimentRuns')\n        D, H, W = volume.shape\n        print(D, H, W)\n\n        probability = np.zeros((7, D, H, W), dtype=np.float32)\n        count = np.zeros((7, D, H, W), dtype=np.float32)\n        pad_volume = np.pad(volume, [[0, 0], [0, 640 - H], [0, 640 - W]], mode='constant', constant_values=0)\n        \n        num_slice=64\n        zz = list(range(0, D - num_slice, num_slice//2)) + [D - num_slice]\n        for z in zz:\n            print('\\r',f'z:{z}', end='',flush=True)\n            image = pad_volume[z:z + num_slice]\n            batch = dotdict(\n                image=torch.from_numpy(image).unsqueeze(0),\n            )\n            with torch.amp.autocast('cuda', enabled=True):\n                with torch.no_grad():\n                    output = net(batch)\n            prob = output['particle'][0].cpu().numpy()\n            probability[:, z:z + num_slice] += prob[:, :, :H, :W]\n            count[:, z:z + num_slice] += 1\n        probability = probability / (count + 0.0001)\n        location = probability_to_location(probability, cfg)\n        df = location_to_df(location)\n        df.insert(loc=0, column='experiment', value=id)\n        submit_df.append(df)\n        print('')\n        print(time_to_str(timer() - start_timer, 'sec'))\n\n    print('\\ndone!')\n    total_time = timer() - start_timer\n    num_volume = len(valid_id)\n    print(f'Total time for {num_volume} volumes:', time_to_str(total_time, 'min'))\n    print(f'Total time for 500 volumes:', time_to_str(total_time/num_volume*500, 'min'))\n    print('')\n    submit_df = pd.concat(submit_df)\n    submit_df.insert(loc=0, column='id', value=np.arange(len(submit_df)))\n    return submit_df\n\n\n\nif 1:\n    submit_df = run_submit()\n    print('submit_df', submit_df.shape)\n    print(submit_df)\n    submit_df.to_csv('submission.csv', index=False)\n\n#!ls\nprint('SUBMIT OK!!!')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-11T18:53:02.217663Z","iopub.execute_input":"2024-11-11T18:53:02.217985Z","iopub.status.idle":"2024-11-11T18:53:28.027722Z","shell.execute_reply.started":"2024-11-11T18:53:02.217951Z","shell.execute_reply":"2024-11-11T18:53:28.026830Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def do_one_eval(truth, predict, threshold):\n    P=len(predict)\n    T=len(truth)\n\n    if P==0:\n        hit=[[],[]]\n        miss=np.arange(T).tolist()\n        fp=[]\n        metric = [P,T,len(hit[0]),len(miss),len(fp)]\n        return hit, fp, miss, metric\n\n    if T==0:\n        hit=[[],[]]\n        fp=np.arange(P).tolist()\n        miss=[]\n        metric = [P,T,len(hit[0]),len(miss),len(fp)]\n        return hit, fp, miss, metric\n\n    #---\n    distance = predict.reshape(P,1,3)-truth.reshape(1,T,3)\n    distance = distance**2\n    distance = distance.sum(axis=2)\n    distance = np.sqrt(distance)\n    p_index, t_index = linear_sum_assignment(distance)\n\n    valid = distance[p_index, t_index] <= threshold\n    p_index = p_index[valid]\n    t_index = t_index[valid]\n    hit = [p_index.tolist(), t_index.tolist()]\n    miss = np.arange(T)\n    miss = miss[~np.isin(miss,t_index)].tolist()\n    fp = np.arange(P)\n    fp = fp[~np.isin(fp,p_index)].tolist()\n\n    metric = [P,T,len(hit[0]),len(miss),len(fp)] #for lb metric F-beta copmutation\n    return hit, fp, miss, metric\n\n\ndef compute_lb(submit_df, overlay_dir):\n    valid_id = list(submit_df['experiment'].unique())\n    print(valid_id)\n\n    eval_df = []\n    for id in valid_id:\n        truth = read_one_truth(id, overlay_dir) #=f'{valid_dir}/overlay/ExperimentRuns')\n        id_df = submit_df[submit_df['experiment'] == id]\n        for p in PARTICLE:\n            p = dotdict(p)\n            print('\\r', id, p.name, end='', flush=True)\n            xyz_truth = truth[p.name]\n            xyz_predict = id_df[id_df['particle_type'] == p.name][['x', 'y', 'z']].values\n            hit, fp, miss, metric = do_one_eval(xyz_truth, xyz_predict, p.radius* 0.5)\n            eval_df.append(dotdict(\n                id=id, particle_type=p.name,\n                P=metric[0], T=metric[1], hit=metric[2], miss=metric[3], fp=metric[4],\n            ))\n    print('')\n    eval_df = pd.DataFrame(eval_df)\n    gb = eval_df.groupby('particle_type').agg('sum').drop(columns=['id'])\n    gb.loc[:, 'precision'] = gb['hit'] / gb['P']\n    gb.loc[:, 'precision'] = gb['precision'].fillna(0)\n    gb.loc[:, 'recall'] = gb['hit'] / gb['T']\n    gb.loc[:, 'recall'] = gb['recall'].fillna(0)\n    gb.loc[:, 'f-beta4'] = 17 * gb['precision'] * gb['recall'] / (16 * gb['precision'] + gb['recall'])\n    gb.loc[:, 'f-beta4'] = gb['f-beta4'].fillna(0)\n\n    gb = gb.sort_values('particle_type').reset_index(drop=False)\n    # https://www.kaggle.com/competitions/czii-cryo-et-object-identification/discussion/544895\n    gb.loc[:, 'weight'] = [1, 0, 2, 1, 2, 1]\n    lb_score = (gb['f-beta4'] * gb['weight']).sum() / gb['weight'].sum()\n    return gb, lb_score\n\n\n#debug\nif 1:\n    if MODE=='local':\n    #if 1:\n        submit_df=pd.read_csv(\n           'submission.csv'\n            # '/kaggle/input/hengck-czii-cryo-et-weights-01/submission.csv'\n        )\n        gb, lb_score = compute_lb(submit_df, f'{valid_dir}/overlay/ExperimentRuns')\n        print(gb)\n        print('lb_score:',lb_score)\n        print('')\n\n\n        #show one ----------------------------------\n        fig = plt.figure(figsize=(18, 8))\n\n        id = valid_id[0]\n        truth = read_one_truth(id,overlay_dir=f'{valid_dir}/overlay/ExperimentRuns')\n\n        submit_df = submit_df[submit_df['experiment']==id]\n        for p in PARTICLE:\n            p = dotdict(p)\n            xyz_truth = truth[p.name]\n            xyz_predict = submit_df[submit_df['particle_type']==p.name][['x','y','z']].values\n            hit, fp, miss, _ = do_one_eval(xyz_truth, xyz_predict, p.radius)\n            print(id, p.name)\n            print('\\t num truth   :',len(xyz_truth) )\n            print('\\t num predict :',len(xyz_predict) )\n            print('\\t num hit  :',len(hit[0]) )\n            print('\\t num fp   :',len(fp) )\n            print('\\t num miss :',len(miss) )\n\n            ax = fig.add_subplot(2, 3, p.label, projection='3d')\n            if hit[0]:\n                pt = xyz_predict[hit[0]]\n                ax.scatter(pt[:, 0], pt[:, 1], pt[:, 2], alpha=0.5, color='r')\n                pt = xyz_truth[hit[1]]\n                ax.scatter(pt[:,0], pt[:,1], pt[:,2], s=80, facecolors='none', edgecolors='r')\n            if fp:\n                pt = xyz_predict[fp]\n                ax.scatter(pt[:, 0], pt[:, 1], pt[:, 2], alpha=1, color='k')\n            if miss:\n                pt = xyz_truth[miss]\n                ax.scatter(pt[:, 0], pt[:, 1], pt[:, 2], s=160, alpha=1, facecolors='none', edgecolors='k')\n\n            ax.set_title(f'{p.name} ({p.difficulty})')\n\n        plt.tight_layout()\n        plt.show()\n        \n        #--- \n        zz=0","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-11T18:53:28.029817Z","iopub.execute_input":"2024-11-11T18:53:28.030139Z","iopub.status.idle":"2024-11-11T18:53:29.430246Z","shell.execute_reply.started":"2024-11-11T18:53:28.030107Z","shell.execute_reply":"2024-11-11T18:53:29.429309Z"}},"outputs":[],"execution_count":null}]}