diff --git a/examples/rubin_dp2/CheckDescShearMatch.ipynb b/examples/rubin_dp2/CheckDescShearMatch.ipynb index 08e0355..b0e0c68 100644 --- a/examples/rubin_dp2/CheckDescShearMatch.ipynb +++ b/examples/rubin_dp2/CheckDescShearMatch.ipynb @@ -14,13 +14,23 @@ "id": "1e088a65-bfa6-476d-8c77-9068925d8ba5", "metadata": {}, "outputs": [], - "source": "# Standard imports\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport tables_io\nimport hpmcm" + "source": [ + "# Standard imports\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import tables_io\n", + "import hpmcm" + ] }, { "cell_type": "markdown", "id": "8b53c63c", - "source": "## Configuration\n\nSet the tract number and define input file paths.", - "metadata": {} + "metadata": {}, + "source": [ + "## Configuration\n", + "\n", + "Set the tract number and define input file paths." + ] }, { "cell_type": "code", @@ -28,13 +38,31 @@ "id": "156c3f8c-11d6-465d-8373-4e0bb33f5e74", "metadata": {}, "outputs": [], - "source": "tract = 8006\n\n# Load shear catalogs and match output tables\nshear_desc_in = tables_io.read(f'{tract}/mdet-run-{tract}.parq')\nshear_desc_ns = tables_io.read(f'{tract}/mdet-run-reduced_uncleaned_{tract}_ns.parq')\nshear_desc_1m = tables_io.read(f'{tract}/mdet-run-reduced_uncleaned_{tract}_1m.parq')\nshear_desc_1p = tables_io.read(f'{tract}/mdet-run-reduced_uncleaned_{tract}_1p.parq')\nshear_desc_obj_stats = tables_io.read(f'{tract}/desc_md_match_object_stats.pq')\nshear_desc_obj_assoc = tables_io.read(f'{tract}/desc_md_match_object_assoc.pq')\ninput_files = [\n f'{tract}/mdet-run-reduced_uncleaned_{tract}_ns.parq',\n f'{tract}/mdet-run-reduced_uncleaned_{tract}_1m.parq',\n f'{tract}/mdet-run-reduced_uncleaned_{tract}_1p.parq',\n] \ncatalog_ids = [0, 1, 2]" + "source": [ + "tract = 8006\n", + "\n", + "# Load shear catalogs and match output tables\n", + "shear_desc_in = tables_io.read(f'{tract}/mdet-run-{tract}.parq')\n", + "shear_desc_ns = tables_io.read(f'{tract}/mdet-run-reduced_uncleaned_{tract}_ns.parq')\n", + "shear_desc_1m = tables_io.read(f'{tract}/mdet-run-reduced_uncleaned_{tract}_1m.parq')\n", + "shear_desc_1p = tables_io.read(f'{tract}/mdet-run-reduced_uncleaned_{tract}_1p.parq')\n", + "shear_desc_obj_stats = tables_io.read(f'{tract}/desc_md_match_object_stats.pq')\n", + "shear_desc_obj_assoc = tables_io.read(f'{tract}/desc_md_match_object_assoc.pq')\n", + "input_files = [\n", + " f'{tract}/mdet-run-reduced_uncleaned_{tract}_ns.parq',\n", + " f'{tract}/mdet-run-reduced_uncleaned_{tract}_1m.parq',\n", + " f'{tract}/mdet-run-reduced_uncleaned_{tract}_1p.parq',\n", + "] \n", + "catalog_ids = [0, 1, 2]" + ] }, { "cell_type": "markdown", "id": "a6268c28", - "source": "## Validation plots — positional checks", - "metadata": {} + "metadata": {}, + "source": [ + "## Validation plots — positional checks" + ] }, { "cell_type": "code", @@ -62,7 +90,15 @@ "id": "c8db1973-1088-4a83-8fd0-3a85d99d981f", "metadata": {}, "outputs": [], - "source": "# Join source catalog data onto the object stats table\njoined_shear_obj = hpmcm.output_tables.buildJoinedObjectTable(\n shear_desc_obj_stats,\n shear_desc_obj_assoc,\n input_files,\n catalog_ids,\n)" + "source": [ + "# Join source catalog data onto the object stats table\n", + "joined_shear_obj = hpmcm.output_tables.buildJoinedObjectTable(\n", + " shear_desc_obj_stats,\n", + " shear_desc_obj_assoc,\n", + " input_files,\n", + " catalog_ids,\n", + ")" + ] }, { "cell_type": "code", @@ -126,4 +162,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} \ No newline at end of file +} diff --git a/examples/rubin_dp2/CheckJoins.ipynb b/examples/rubin_dp2/CheckJoins.ipynb index b3962f8..34b1cf6 100644 --- a/examples/rubin_dp2/CheckJoins.ipynb +++ b/examples/rubin_dp2/CheckJoins.ipynb @@ -3,8 +3,12 @@ { "cell_type": "markdown", "id": "e2cba6bd", - "source": "# Check joined tables\n\nVerifies that buildJoinedObjectTable correctly merges object statistics with per-catalog source columns.", - "metadata": {} + "metadata": {}, + "source": [ + "# Check joined tables\n", + "\n", + "Verifies that buildJoinedObjectTable correctly merges object statistics with per-catalog source columns." + ] }, { "cell_type": "code", @@ -12,7 +16,13 @@ "id": "9cb3a6e4-5813-4b95-8b20-0449bc89bfb2", "metadata": {}, "outputs": [], - "source": "# Standard imports\nimport tables_io\nimport hpmcm\nimport matplotlib.pyplot as plt\nimport numpy as np" + "source": [ + "# Standard imports\n", + "import tables_io\n", + "import hpmcm\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np" + ] }, { "cell_type": "code", @@ -20,7 +30,18 @@ "id": "d56eba26-ca76-4cd1-b545-b0c3ef8cc75f", "metadata": {}, "outputs": [], - "source": "tract = 8006\n\n# Input catalogs and their IDs for the meta-match\nall_files = [\n f'{tract}/object_reduced_{tract}.parq',\n f'{tract}/anacal_reduced_{tract}.parq',\n f'{tract}/desc_md_match_object_stats.pq',\n f'{tract}/rubin_md_match_object_stats.pq',\n]\nall_ids = list(range(len(all_files)))" + "source": [ + "tract = 8006\n", + "\n", + "# Input catalogs and their IDs for the meta-match\n", + "all_files = [\n", + " f'{tract}/object_reduced_{tract}.parq',\n", + " f'{tract}/anacal_reduced_{tract}.parq',\n", + " f'{tract}/desc_md_match_object_joined.pq',\n", + " f'{tract}/rubin_md_match_object_joined.pq',\n", + "]\n", + "all_ids = list(range(len(all_files)))" + ] }, { "cell_type": "code", @@ -28,7 +49,11 @@ "id": "9b34104d-db8c-4efe-a2d3-7b4c5bfc0bff", "metadata": {}, "outputs": [], - "source": "# Load association and stats tables\nobject_assoc = hpmcm.output_tables.ObjectAssocTable.read(f'{tract}/all_match_object_assoc.pq', [])\nobject_stats = hpmcm.output_tables.ObjectStatsTable.read(f'{tract}/all_match_object_stats.pq', [])" + "source": [ + "# Load association and stats tables\n", + "object_assoc = hpmcm.output_tables.ObjectAssocTable.read(f'{tract}/all_match_object_assoc.pq', [])\n", + "object_stats = hpmcm.output_tables.ObjectStatsTable.read(f'{tract}/all_match_object_stats.pq', [])" + ] }, { "cell_type": "code", @@ -36,7 +61,15 @@ "id": "ed073cb6-005c-4e0e-95aa-05ef3f67161c", "metadata": {}, "outputs": [], - "source": "# Build wide joined table — one row per object, source cols suffixed by catalog_id\njoined = hpmcm.output_tables.buildJoinedObjectTable(\n object_assoc,\n all_files,\n all_ids,\n)" + "source": [ + "# Build wide joined table — one row per object, source cols suffixed by catalog_id\n", + "joined = hpmcm.output_tables.buildJoinedObjectTable(\n", + " object_stats,\n", + " object_assoc,\n", + " all_files,\n", + " all_ids,\n", + ")" + ] }, { "cell_type": "code", @@ -44,7 +77,10 @@ "id": "5688ec6c-8981-432f-99da-b0e5929d5cfd", "metadata": {}, "outputs": [], - "source": "# Count objects detected in catalog 0 (bit 0 set)\n(np.bitwise_and(joined['catalog_mask'],1) == 1).sum()" + "source": [ + "# Count objects detected in catalog 0 (bit 0 set)\n", + "(np.bitwise_and(joined['catalog_mask'],1) == 1).sum()" + ] }, { "cell_type": "code", @@ -132,7 +168,13 @@ "id": "3341a247-cf9b-4a26-8d71-2ebcdc3ed687", "metadata": {}, "outputs": [], - "source": "# Histogram of fluxes split by which catalogs detected each object\nfor i in range(1,16):\n _ = plt.hist(joined['g_gaapPsfFlux_0'][joined['catalog_mask']==i], bins=np.logspace(0, 6, 61))\n_ = plt.xscale('log')\n_ = plt.yscale('log')" + "source": [ + "# Histogram of fluxes split by which catalogs detected each object\n", + "for i in range(1,16):\n", + " _ = plt.hist(joined['g_gaapPsfFlux_0'][joined['catalog_mask']==i], bins=np.logspace(0, 6, 61))\n", + "_ = plt.xscale('log')\n", + "_ = plt.yscale('log')" + ] }, { "cell_type": "code", @@ -164,4 +206,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} \ No newline at end of file +} diff --git a/examples/rubin_dp2/CheckObjectMatch.ipynb b/examples/rubin_dp2/CheckObjectMatch.ipynb index 2147607..07c65ab 100644 --- a/examples/rubin_dp2/CheckObjectMatch.ipynb +++ b/examples/rubin_dp2/CheckObjectMatch.ipynb @@ -3,8 +3,12 @@ { "cell_type": "markdown", "id": "4b91779a", - "source": "# Check object matching results\n\nExamines per-catalog coverage, SNR distributions, and spatial layouts of the multi-catalog object match.", - "metadata": {} + "metadata": {}, + "source": [ + "# Check object matching results\n", + "\n", + "Examines per-catalog coverage, SNR distributions, and spatial layouts of the multi-catalog object match." + ] }, { "cell_type": "code", @@ -12,7 +16,13 @@ "id": "b6d5608f-ba28-4a27-9c17-20da93f778cc", "metadata": {}, "outputs": [], - "source": "# Standard imports\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport tables_io\nimport hpmcm" + "source": [ + "# Standard imports\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import tables_io\n", + "import hpmcm" + ] }, { "cell_type": "code", @@ -22,7 +32,25 @@ "scrolled": true }, "outputs": [], - "source": "tract = 8006\n\n# Input catalog files and IDs for the object meta-match\nall_obj_stats = tables_io.read(f'{tract}/all_match_object_stats.pq')\nall_obj_assoc = tables_io.read(f'{tract}/all_match_object_assoc.pq')\nobject_in = tables_io.read(f'{tract}/object_reduced_{tract}.parq')\nanacal_in = tables_io.read(f'{tract}/anacal_reduced_{tract}.parq')\ndesc_md_in = tables_io.read(f'{tract}/desc_md_match_object_stats.pq')\nrubin_md_in = tables_io.read(f'{tract}/rubin_md_match_object_stats.pq')\n\nall_files = [\n f'{tract}/object_reduced_{tract}.parq',\n f'{tract}/anacal_reduced_{tract}.parq',\n f'{tract}/desc_md_match_object_stats.pq',\n f'{tract}/rubin_md_match_object_stats.pq',\n]\nall_ids = list(range(len(all_files)))" + "source": [ + "tract = 8006\n", + "\n", + "# Input catalog files and IDs for the object meta-match\n", + "all_obj_stats = tables_io.read(f'{tract}/all_match_object_stats.pq')\n", + "all_obj_assoc = tables_io.read(f'{tract}/all_match_object_assoc.pq')\n", + "object_in = tables_io.read(f'{tract}/object_reduced_{tract}.parq')\n", + "anacal_in = tables_io.read(f'{tract}/anacal_reduced_{tract}.parq')\n", + "desc_md_in = tables_io.read(f'{tract}/desc_md_match_object_stats.pq')\n", + "rubin_md_in = tables_io.read(f'{tract}/rubin_md_match_object_stats.pq')\n", + "\n", + "all_files = [\n", + " f'{tract}/object_reduced_{tract}.parq',\n", + " f'{tract}/anacal_reduced_{tract}.parq',\n", + " f'{tract}/desc_md_match_object_stats.pq',\n", + " f'{tract}/rubin_md_match_object_stats.pq',\n", + "]\n", + "all_ids = list(range(len(all_files)))" + ] }, { "cell_type": "code", @@ -40,7 +68,15 @@ "id": "f37611d3-aa7c-42c6-b925-63e7854f1b80", "metadata": {}, "outputs": [], - "source": "# Build wide joined table\njoined_obj = hpmcm.output_tables.buildJoinedObjectTable(\n all_obj_stats,\n all_obj_assoc,\n all_files,\n all_ids,\n)" + "source": [ + "# Build wide joined table\n", + "joined_obj = hpmcm.output_tables.buildJoinedObjectTable(\n", + " all_obj_stats,\n", + " all_obj_assoc,\n", + " all_files,\n", + " all_ids,\n", + ")" + ] }, { "cell_type": "code", @@ -48,7 +84,16 @@ "id": "3d658e93-ee0a-45e0-bd7e-e73925c22cf6", "metadata": {}, "outputs": [], - "source": "# SNR distributions split by which catalogs detected each object\nbins = np.logspace(0, 5, 51)\nfor i in range(4):\n mask_bit = 1 << i\n mask = np.bitwise_and(joined_obj.catalog_mask, mask_bit) == mask_bit\n _ = plt.hist(joined_obj[mask].snr, bins=bins, alpha=0.1)\n_ = plt.xscale('log')\n_ = plt.yscale('log')" + "source": [ + "# SNR distributions split by which catalogs detected each object\n", + "bins = np.logspace(0, 5, 51)\n", + "for i in range(4):\n", + " mask_bit = 1 << i\n", + " mask = np.bitwise_and(joined_obj.catalog_mask, mask_bit) == mask_bit\n", + " _ = plt.hist(joined_obj[mask].snr, bins=bins, alpha=0.1)\n", + "_ = plt.xscale('log')\n", + "_ = plt.yscale('log')" + ] }, { "cell_type": "code", @@ -56,7 +101,11 @@ "id": "88a6fdeb-fb5f-426d-9835-7b5b75fe63fd", "metadata": {}, "outputs": [], - "source": "# 2D histogram: which catalogs detected each object vs SNR\n_ = plt.hist2d(joined_obj.catalog_mask, joined_obj.snr, bins=(np.linspace(0.5, 15.5, 16), np.logspace(0, 5, 51)), norm='log')\n_ = plt.yscale('log')" + "source": [ + "# 2D histogram: which catalogs detected each object vs SNR\n", + "_ = plt.hist2d(joined_obj.catalog_mask, joined_obj.snr, bins=(np.linspace(0.5, 15.5, 16), np.logspace(0, 5, 51)), norm='log')\n", + "_ = plt.yscale('log')" + ] }, { "cell_type": "code", @@ -99,4 +148,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} \ No newline at end of file +} diff --git a/examples/rubin_dp2/CheckReducedInputs.ipynb b/examples/rubin_dp2/CheckReducedInputs.ipynb index 0471991..99848d9 100644 --- a/examples/rubin_dp2/CheckReducedInputs.ipynb +++ b/examples/rubin_dp2/CheckReducedInputs.ipynb @@ -3,8 +3,12 @@ { "cell_type": "markdown", "id": "1171dc87", - "source": "# Check reduced inputs\n\nValidates that reducing input tables preserves coordinates by comparing RA/Dec before and after reduction for Anacal, Object, and Shear catalogs.", - "metadata": {} + "metadata": {}, + "source": [ + "# Check reduced inputs\n", + "\n", + "Validates that reducing input tables preserves coordinates by comparing RA/Dec before and after reduction for Anacal, Object, and Shear catalogs." + ] }, { "cell_type": "code", @@ -12,21 +16,30 @@ "id": "d9a5fd3c-f16a-4d43-badc-00917f2e47df", "metadata": {}, "outputs": [], - "source": "# Standard imports\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport tables_io" + "source": [ + "# Standard imports\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import tables_io" + ] }, { "cell_type": "markdown", "id": "fca79fbf", - "source": "## Configuration", - "metadata": {} + "metadata": {}, + "source": [ + "## Configuration" + ] }, { "cell_type": "code", + "execution_count": null, "id": "16490dd7", - "source": "tract = 8006", "metadata": {}, - "execution_count": null, - "outputs": [] + "outputs": [], + "source": [ + "tract = 8006" + ] }, { "cell_type": "code", @@ -34,7 +47,13 @@ "id": "2013f078-201c-4cb6-a2e8-bbed32f728bf", "metadata": {}, "outputs": [], - "source": "# Load raw and reduced Anacal tables\nanacal_input = tables_io.read(f'{tract}/anacal_merged_{tract}.parq')\nanacal_reduced = tables_io.read(f'{tract}/anacal_reduced_{tract}.parq')\nanacal_input.index = anacal_input.object_id\nanacal_reduced.index = anacal_reduced.id" + "source": [ + "# Load raw and reduced Anacal tables\n", + "anacal_input = tables_io.read(f'{tract}/anacal_merged_{tract}.parq')\n", + "anacal_reduced = tables_io.read(f'{tract}/anacal_reduced_{tract}.parq')\n", + "anacal_input.index = anacal_input.object_id\n", + "anacal_reduced.index = anacal_reduced.id" + ] }, { "cell_type": "code", @@ -42,7 +61,10 @@ "id": "718ca408-0abe-4a83-aadc-0e2ab727d252", "metadata": {}, "outputs": [], - "source": "# Match on id to compare raw vs reduced\nanacal_joined = anacal_input.join(anacal_reduced, how='left', lsuffix='_in', rsuffix='_red')" + "source": [ + "# Match on id to compare raw vs reduced\n", + "anacal_joined = anacal_input.join(anacal_reduced, how='left', lsuffix='_in', rsuffix='_red')" + ] }, { "cell_type": "code", @@ -50,7 +72,13 @@ "id": "fc418fab-21a6-4858-a4de-419f86545ec7", "metadata": {}, "outputs": [], - "source": "# Check RA/Dec round-trip error for Anacal\nra_diff = anacal_joined.ra_in - anacal_joined.ra_red\ndec_diff = anacal_joined.dec_in - anacal_joined.dec_red\nprint(f\"RA {ra_diff.min()} - {ra_diff.max()}\")\nprint(f\"DEC {dec_diff.min()} - {dec_diff.max()}\")" + "source": [ + "# Check RA/Dec round-trip error for Anacal\n", + "ra_diff = anacal_joined.ra_in - anacal_joined.ra_red\n", + "dec_diff = anacal_joined.dec_in - anacal_joined.dec_red\n", + "print(f\"RA {ra_diff.min()} - {ra_diff.max()}\")\n", + "print(f\"DEC {dec_diff.min()} - {dec_diff.max()}\")" + ] }, { "cell_type": "code", @@ -58,7 +86,13 @@ "id": "67571a52-077f-4e43-86b9-b102ffe056c0", "metadata": {}, "outputs": [], - "source": "# Load raw and reduced Object catalog tables\nobject_in = tables_io.read(f'{tract}/object_{tract}.parq')\nobject_red = tables_io.read(f'{tract}/object_reduced_{tract}.parq')\nobject_in.index = object_in.objectId\nobject_red.index = object_red.id" + "source": [ + "# Load raw and reduced Object catalog tables\n", + "object_in = tables_io.read(f'{tract}/object_{tract}.parq')\n", + "object_red = tables_io.read(f'{tract}/object_reduced_{tract}.parq')\n", + "object_in.index = object_in.objectId\n", + "object_red.index = object_red.id" + ] }, { "cell_type": "code", @@ -89,7 +123,15 @@ "id": "ec596e7b-72d9-40d0-bb6c-41be9c418cc0", "metadata": {}, "outputs": [], - "source": "# Load raw and reduced Rubin shear tables\nshear_rubin_in = tables_io.read(f'{tract}/shear_object_{tract}.parq')\nshear_rubin_ns = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_ns.parq')\nshear_rubin_1m = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_1m.parq')\nshear_rubin_1p = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_1p.parq')\nshear_rubin_2m = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_2m.parq')\nshear_rubin_2p = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_2p.parq') \n" + "source": [ + "# Load raw and reduced Rubin shear tables\n", + "shear_rubin_in = tables_io.read(f'{tract}/shear_object_{tract}.parq')\n", + "shear_rubin_ns = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_ns.parq')\n", + "shear_rubin_1m = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_1m.parq')\n", + "shear_rubin_1p = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_1p.parq')\n", + "shear_rubin_2m = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_2m.parq')\n", + "shear_rubin_2p = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_2p.parq') \n" + ] }, { "cell_type": "code", @@ -201,7 +243,17 @@ "id": "9f820c45-e862-4449-8d5f-7453b32ade47", "metadata": {}, "outputs": [], - "source": "# Load raw and reduced DESC shear tables\nshear_desc_in = tables_io.read(f'{tract}/mdet-run-{tract}.parq')\nshear_desc_ns = tables_io.read(f'{tract}/mdet-run-{tract}_uncleaned_{tract}_ns.parq')\nshear_desc_1m = tables_io.read(f'{tract}/mdet-run-{tract}_uncleaned_{tract}_1m.parq')\nshear_desc_1p = tables_io.read(f'{tract}/mdet-run-{tract}_uncleaned_{tract}_1p.parq')\nshear_desc_in.index = shear_desc_in.id\nshear_desc_ns.index = shear_desc_ns.id\nshear_desc_1m.index = shear_desc_1m.id\nshear_desc_1p.index = shear_desc_1p.id" + "source": [ + "# Load raw and reduced DESC shear tables\n", + "shear_desc_in = tables_io.read(f'{tract}/mdet-run-{tract}.parq')\n", + "shear_desc_ns = tables_io.read(f'{tract}/mdet-run-{tract}_uncleaned_{tract}_ns.parq')\n", + "shear_desc_1m = tables_io.read(f'{tract}/mdet-run-{tract}_uncleaned_{tract}_1m.parq')\n", + "shear_desc_1p = tables_io.read(f'{tract}/mdet-run-{tract}_uncleaned_{tract}_1p.parq')\n", + "shear_desc_in.index = shear_desc_in.id\n", + "shear_desc_ns.index = shear_desc_ns.id\n", + "shear_desc_1m.index = shear_desc_1m.id\n", + "shear_desc_1p.index = shear_desc_1p.id" + ] }, { "cell_type": "code", @@ -304,4 +356,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} \ No newline at end of file +} diff --git a/examples/rubin_dp2/CheckRubinShearMatch.ipynb b/examples/rubin_dp2/CheckRubinShearMatch.ipynb index 949746e..eb83104 100644 --- a/examples/rubin_dp2/CheckRubinShearMatch.ipynb +++ b/examples/rubin_dp2/CheckRubinShearMatch.ipynb @@ -14,7 +14,13 @@ "id": "3446af0f-9c99-44c4-8ea1-fdd0df73fba9", "metadata": {}, "outputs": [], - "source": "# Standard imports\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport tables_io\nimport hpmcm" + "source": [ + "# Standard imports\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import tables_io\n", + "import hpmcm" + ] }, { "cell_type": "code", @@ -33,13 +39,33 @@ "id": "c702bf68-25e1-453e-b191-668deb4d3197", "metadata": {}, "outputs": [], - "source": "# Input shear catalog files and matching output tables\nshear_rubin_in = tables_io.read('{tract}/shear_object_{tract}.parq')\nshear_rubin_ns = tables_io.read('{tract}/shear_object_reduced_uncleaned_{tract}_ns.parq')\nshear_rubin_1m = tables_io.read('{tract}/shear_object_reduced_uncleaned_{tract}_1m.parq')\nshear_rubin_1p = tables_io.read('{tract}/shear_object_reduced_uncleaned_{tract}_1p.parq')\nshear_rubin_2m = tables_io.read('{tract}/shear_object_reduced_uncleaned_{tract}2m.parq')\nshear_rubin_2p = tables_io.read('{tract}/shear_object_reduced_uncleaned_{tract}_2p.parq') \nshear_rubin_obj_stats = tables_io.read(f'{tract}/rubin_md_match_object_stats.pq')\nshear_rubin_obj_assoc = tables_io.read(f'{tract}/rubin_md_match_object_assoc.pq')\ninput_files = [\n '{tract}/shear_object_reduced_uncleaned_{tract}_ns.parq',\n '{tract}/shear_object_reduced_uncleaned_{tract}_1m.parq',\n '{tract}/shear_object_reduced_uncleaned_{tract}_1p.parq',\n '{tract}/shear_object_reduced_uncleaned_{tract}_2m.parq',\n '{tract}/shear_object_reduced_uncleaned_{tract}_2p.parq',\n] \ncatalog_ids = [0, 1, 2, 3, 4]" + "source": [ + "# Input shear catalog files and matching output tables\n", + "shear_rubin_in = tables_io.read(f'{tract}/shear_object_{tract}.parq')\n", + "shear_rubin_ns = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_ns.parq')\n", + "shear_rubin_1m = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_1m.parq')\n", + "shear_rubin_1p = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_1p.parq')\n", + "shear_rubin_2m = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_2m.parq')\n", + "shear_rubin_2p = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_2p.parq') \n", + "shear_rubin_obj_stats = tables_io.read(f'{tract}/rubin_md_match_object_stats.pq')\n", + "shear_rubin_obj_assoc = tables_io.read(f'{tract}/rubin_md_match_object_assoc.pq')\n", + "input_files = [\n", + " f'{tract}/shear_object_reduced_uncleaned_{tract}_ns.parq',\n", + " f'{tract}/shear_object_reduced_uncleaned_{tract}_1m.parq',\n", + " f'{tract}/shear_object_reduced_uncleaned_{tract}_1p.parq',\n", + " f'{tract}/shear_object_reduced_uncleaned_{tract}_2m.parq',\n", + " f'{tract}/shear_object_reduced_uncleaned_{tract}_2p.parq',\n", + "] \n", + "catalog_ids = [0, 1, 2, 3, 4]" + ] }, { "cell_type": "markdown", "id": "8b62b38e", - "source": "## Validation plots", - "metadata": {} + "metadata": {}, + "source": [ + "## Validation plots" + ] }, { "cell_type": "code", @@ -138,4 +164,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} \ No newline at end of file +} diff --git a/examples/rubin_dp2/RunMatching.ipynb b/examples/rubin_dp2/RunMatching.ipynb index bdcfb1b..30c1b8a 100644 --- a/examples/rubin_dp2/RunMatching.ipynb +++ b/examples/rubin_dp2/RunMatching.ipynb @@ -34,7 +34,7 @@ "# Set the tract\n", "tract = 8006\n", "# This will need to change with tract\n", - "ref_dir=(340.5000000000, -9.6694214876),\n", + "ref_dir=(340.5000000000, -9.6694214876)\n", "our_geom = hpmcm.shear_utils.ShearCellGeometry(ref_dir=ref_dir)" ] }, @@ -78,8 +78,8 @@ "all_files = [\n", " f'{tract}/object_reduced_{tract}.parq',\n", " f'{tract}/anacal_reduced_{tract}.parq',\n", - " f'{tract}/desc_md_match_object_stats.pq',\n", - " f'{tract}/rubin_md_match_object_stats.pq',\n", + " f'{tract}/desc_md_match_object_joined.pq',\n", + " f'{tract}/rubin_md_match_object_joined.pq',\n", "]\n", "all_ids = list(range(len(all_files)))" ] @@ -92,8 +92,8 @@ "outputs": [], "source": [ "# Allow for re-running various matchers\n", - "rematch_rubin_shear = True\n", - "rematch_desc_shear = True\n", + "rematch_rubin_shear = False\n", + "rematch_desc_shear = False\n", "rematch_objects = True\n", "\n", "if rematch_rubin_shear:\n", @@ -111,7 +111,7 @@ " rubin_stats.update(**rubin_shear_match.extractShearStats())\n", " tables_io.write(rubin_stats, f'{tract}/rubin_md_match_')\n", " rubin_shear_match.save(f'{tract}/rubin_md_match')\n", - "\n", + " \n", "\n", "if rematch_desc_shear:\n", " desc_shear_match = hpmcm.shear_match.ShearMatch.createShearMatch(\n", @@ -128,29 +128,79 @@ " desc_stats.update(**desc_shear_match.extractShearStats())\n", " tables_io.write(desc_stats, f'{tract}/desc_md_match_') \n", " desc_shear_match.save(f'{tract}/desc_md_match')\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "608b0387-feb2-42e1-aa47-bf93ceb39030", + "metadata": {}, + "outputs": [], + "source": [ + "shear_rubin_obj_stats = tables_io.read(f'{tract}/rubin_md_match_object_stats.pq')\n", + "shear_rubin_obj_assoc = tables_io.read(f'{tract}/rubin_md_match_object_assoc.pq')\n", + "shear_rubin_obj_shear = tables_io.read(f'{tract}/rubin_md_match_object_shear.pq')\n", "\n", + "joined_rubin_shear_obj = hpmcm.output_tables.buildJoinedObjectTable(\n", + " shear_rubin_obj_stats,\n", + " shear_rubin_obj_assoc,\n", + " rubin_shear_match_files,\n", + " rubin_shear_match_ids,\n", + " object_shear=f'{tract}/rubin_md_match_object_shear.pq',\n", + ")\n", + "stats_cols = [f\"{band}_gaussFlux\" for band in 'griz']\n", + "stats_cols += [f\"{band}_gaussFluxErr\" for band in 'griz']\n", + "drop_cols = []\n", + "for col_prefix in [\n", + " 'tract', 'ra', 'dec', 'gauss_snr', 'patch_x', 'patch_y', \n", + " 'cell_idx_x', 'cell_idx_y', 'shear', 'meta_step',\n", + "]:\n", + " drop_cols += [f\"{col_prefix}_{id}\" for id in rubin_shear_match_ids]\n", "\n", - "if rematch_objects:\n", - " pixel_scale = 0.2\n", - " n_pix_region = 30000\n", - " region_size_deg = n_pix_region*pixel_scale/3600\n", - " wcs_match = hpmcm.WcsMatch.create(\n", - " ref_dir=ref_dir,\n", - " region_size=(region_size_deg, region_size_deg),\n", - " pixel_size=pixel_scale/3600,\n", - " cell_size=150,\n", - " cell_buffer=5,\n", - " cell_max_object=10000,\n", - " )\n", - " \n", - " wcs_match.reduceData(all_files, all_ids)\n", - " if test_only:\n", - " wcs_match.analysisLoop(range(100, 101), range(100, 101))\n", - " else:\n", - " wcs_match.analysisLoop()\n", - " stats = wcs_match.extractStats()\n", - " tables_io.write(stats, f'{tract}/all_match_') \n", - " wcs_match.save(f'{tract}/all_match')" + "rd_joined_rubin_shear_obj = hpmcm.reduceJoinedTable(\n", + " joined_rubin_shear_obj,\n", + " catalog_ids=rubin_shear_match_ids,\n", + " drop_cols=drop_cols,\n", + " stats_cols=stats_cols,\n", + ")\n", + "tables_io.write(rd_joined_rubin_shear_obj, f'{tract}/rubin_md_match_object_joined.pq')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0ccc0274-d134-434b-9865-31e2a9416003", + "metadata": {}, + "outputs": [], + "source": [ + "shear_desc_obj_stats = tables_io.read(f'{tract}/desc_md_match_object_stats.pq')\n", + "shear_desc_obj_assoc = tables_io.read(f'{tract}/desc_md_match_object_assoc.pq')\n", + "shear_desc_obj_shear = tables_io.read(f'{tract}/desc_md_match_object_shear.pq')\n", + "\n", + "joined_desc_shear_obj = hpmcm.output_tables.buildJoinedObjectTable(\n", + " shear_desc_obj_stats,\n", + " shear_desc_obj_assoc,\n", + " desc_shear_match_files,\n", + " desc_shear_match_ids,\n", + " object_shear=f'{tract}/desc_md_match_object_shear.pq',\n", + ")\n", + "stats_cols = [f\"flux_{band}\" for band in 'riz']\n", + "stats_cols += [f\"flux_err_{band}\" for band in 'riz']\n", + "drop_cols = []\n", + "for col_prefix in [\n", + " 'ra', 'dec', 'gauss_snr', 'patch_x', 'patch_y', \n", + " 'cell_idx_x', 'cell_idx_y', 'shear', 'meta_step',\n", + "]:\n", + " drop_cols += [f\"{col_prefix}_{id}\" for id in rubin_shear_match_ids]\n", + "\n", + "rd_joined_desc_shear_obj = hpmcm.reduceJoinedTable(\n", + " joined_desc_shear_obj,\n", + " catalog_ids=desc_shear_match_ids,\n", + " drop_cols=drop_cols,\n", + " stats_cols=stats_cols,\n", + ")\n", + "tables_io.write(rd_joined_desc_shear_obj, f'{tract}/desc_md_match_object_joined.pq')" ] }, { @@ -191,6 +241,36 @@ ")" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "5bd02341-0983-4dc5-85ab-368230a0ac41", + "metadata": {}, + "outputs": [], + "source": [ + "if rematch_objects:\n", + " pixel_scale = 0.2\n", + " n_pix_region = 30000\n", + " region_size_deg = n_pix_region*pixel_scale/3600\n", + " wcs_match = hpmcm.WcsMatch.create(\n", + " ref_dir=ref_dir,\n", + " region_size=(region_size_deg, region_size_deg),\n", + " pixel_size=pixel_scale/3600,\n", + " cell_size=150,\n", + " cell_buffer=5,\n", + " cell_max_object=10000,\n", + " )\n", + " \n", + " wcs_match.reduceData(all_files, all_ids)\n", + " if test_only:\n", + " wcs_match.analysisLoop(range(100, 101), range(100, 101))\n", + " else:\n", + " wcs_match.analysisLoop()\n", + " stats = wcs_match.extractStats()\n", + " tables_io.write(stats, f'{tract}/all_match_') \n", + " wcs_match.save(f'{tract}/all_match')" + ] + }, { "cell_type": "code", "execution_count": null, @@ -209,6 +289,14 @@ " cell_data,\n", ")" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f588b1fc-64ab-4d4d-8d6c-d4b8e6d63e3e", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/examples/rubin_dp2/StudyDescShearMatch.ipynb b/examples/rubin_dp2/StudyDescShearMatch.ipynb index 27a9242..d2bb9bb 100644 --- a/examples/rubin_dp2/StudyDescShearMatch.ipynb +++ b/examples/rubin_dp2/StudyDescShearMatch.ipynb @@ -3,8 +3,12 @@ { "cell_type": "markdown", "id": "4fe425c9", - "source": "# Study DESC shear matching\n\nAnalyzes the DESC metadetect shear match: SNR distributions, catalog coverage, and spatial distributions of good vs bad matches.", - "metadata": {} + "metadata": {}, + "source": [ + "# Study DESC shear matching\n", + "\n", + "Analyzes the DESC metadetect shear match: SNR distributions, catalog coverage, and spatial distributions of good vs bad matches." + ] }, { "cell_type": "code", @@ -12,21 +16,31 @@ "id": "18d7b7d3-9d40-47e4-82a0-326b0987ba78", "metadata": {}, "outputs": [], - "source": "# Standard imports\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport tables_io\nimport hpmcm" + "source": [ + "# Standard imports\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import tables_io\n", + "import hpmcm" + ] }, { "cell_type": "markdown", "id": "f8dc62ad", - "source": "## Configuration", - "metadata": {} + "metadata": {}, + "source": [ + "## Configuration" + ] }, { "cell_type": "code", + "execution_count": null, "id": "5d281432", - "source": "tract = 8006", "metadata": {}, - "execution_count": null, - "outputs": [] + "outputs": [], + "source": [ + "tract = 8006" + ] }, { "cell_type": "code", @@ -34,7 +48,21 @@ "id": "11517fff-fb91-45bc-9544-fe6be9cdaf91", "metadata": {}, "outputs": [], - "source": "# Load DESC shear match output tables and input catalog files\nshear_desc_in = tables_io.read(f'{tract}/mdet-run-{tract}.parq')\nshear_desc_ns = tables_io.read(f'{tract}/mdet-run-reduced_uncleaned_{tract}_ns.parq')\nshear_desc_1m = tables_io.read(f'{tract}/mdet-run-reduced_uncleaned_{tract}_1m.parq')\nshear_desc_1p = tables_io.read(f'{tract}/mdet-run-reduced_uncleaned_{tract}_1p.parq')\nshear_desc_obj_stats = tables_io.read(f'{tract}/desc_md_match_object_stats.pq')\nshear_desc_obj_assoc = tables_io.read(f'{tract}/desc_md_match_object_assoc.pq')\ninput_files = [\n f'{tract}/mdet-run-reduced_uncleaned_{tract}_ns.parq',\n f'{tract}/mdet-run-reduced_uncleaned_{tract}_1m.parq',\n f'{tract}/mdet-run-reduced_uncleaned_{tract}_1p.parq',\n] \ncatalog_ids = [0, 1, 2]" + "source": [ + "# Load DESC shear match output tables and input catalog files\n", + "shear_desc_in = tables_io.read(f'{tract}/mdet-run-{tract}.parq')\n", + "shear_desc_ns = tables_io.read(f'{tract}/mdet-run-reduced_uncleaned_{tract}_ns.parq')\n", + "shear_desc_1m = tables_io.read(f'{tract}/mdet-run-reduced_uncleaned_{tract}_1m.parq')\n", + "shear_desc_1p = tables_io.read(f'{tract}/mdet-run-reduced_uncleaned_{tract}_1p.parq')\n", + "shear_desc_obj_stats = tables_io.read(f'{tract}/desc_md_match_object_stats.pq')\n", + "shear_desc_obj_assoc = tables_io.read(f'{tract}/desc_md_match_object_assoc.pq')\n", + "input_files = [\n", + " f'{tract}/mdet-run-reduced_uncleaned_{tract}_ns.parq',\n", + " f'{tract}/mdet-run-reduced_uncleaned_{tract}_1m.parq',\n", + " f'{tract}/mdet-run-reduced_uncleaned_{tract}_1p.parq',\n", + "] \n", + "catalog_ids = [0, 1, 2]" + ] }, { "cell_type": "code", @@ -42,7 +70,15 @@ "id": "12fe8065-a3c0-4619-afb5-bea47f3c2aa5", "metadata": {}, "outputs": [], - "source": "# Build wide joined table — one row per object with all catalog columns\njoined_shear_obj = hpmcm.output_tables.buildJoinedObjectTable(\n shear_desc_obj_stats,\n shear_desc_obj_assoc,\n input_files,\n catalog_ids,\n)" + "source": [ + "# Build wide joined table — one row per object with all catalog columns\n", + "joined_shear_obj = hpmcm.output_tables.buildJoinedObjectTable(\n", + " shear_desc_obj_stats,\n", + " shear_desc_obj_assoc,\n", + " input_files,\n", + " catalog_ids,\n", + ")" + ] }, { "cell_type": "code", @@ -50,7 +86,16 @@ "id": "287813a8-e429-49be-8de1-3ac30427caa7", "metadata": {}, "outputs": [], - "source": "# SNR histograms split by catalog detection bits\nbins = np.logspace(0, 5, 51)\nfor i in range(3):\n mask_bit = 1 << i\n mask = np.bitwise_and(joined_shear_obj.catalog_mask, mask_bit) == mask_bit\n _ = plt.hist(joined_shear_obj[mask].snr, bins=bins, alpha=0.1)\n_ = plt.xscale('log')\n_ = plt.yscale('log')" + "source": [ + "# SNR histograms split by catalog detection bits\n", + "bins = np.logspace(0, 5, 51)\n", + "for i in range(3):\n", + " mask_bit = 1 << i\n", + " mask = np.bitwise_and(joined_shear_obj.catalog_mask, mask_bit) == mask_bit\n", + " _ = plt.hist(joined_shear_obj[mask].snr, bins=bins, alpha=0.1)\n", + "_ = plt.xscale('log')\n", + "_ = plt.yscale('log')" + ] }, { "cell_type": "code", @@ -69,7 +114,13 @@ "id": "aea91454-7a5e-40c2-ad48-4ed5f1e2f403", "metadata": {}, "outputs": [], - "source": "# Objects detected in all 3 DESC catalogs (ns, 1m, 1p) have catalog_mask == 7\ngood_mask = joined_shear_obj.catalog_mask == 7\nbad_mask = joined_shear_obj.catalog_mask != 7\ngood_sub = joined_shear_obj[good_mask]\nbad_sub = joined_shear_obj[bad_mask]" + "source": [ + "# Objects detected in all 3 DESC catalogs (ns, 1m, 1p) have catalog_mask == 7\n", + "good_mask = joined_shear_obj.catalog_mask == 7\n", + "bad_mask = joined_shear_obj.catalog_mask != 7\n", + "good_sub = joined_shear_obj[good_mask]\n", + "bad_sub = joined_shear_obj[bad_mask]" + ] }, { "cell_type": "code", @@ -77,7 +128,11 @@ "id": "aa364bca-9b82-4a9b-b330-b9a42a21a632", "metadata": {}, "outputs": [], - "source": "# Spatial distribution of good/bad matches in cell 1\n_ = plt.scatter(good_sub.x_cell_coadd_1[::10], good_sub.y_cell_coadd_1[::10], s=1)\n_ = plt.scatter(bad_sub.x_cell_coadd_1[::5], bad_sub.y_cell_coadd_1[::5], s=1)" + "source": [ + "# Spatial distribution of good/bad matches in cell 1\n", + "_ = plt.scatter(good_sub.x_cell_coadd_1[::10], good_sub.y_cell_coadd_1[::10], s=1)\n", + "_ = plt.scatter(bad_sub.x_cell_coadd_1[::5], bad_sub.y_cell_coadd_1[::5], s=1)" + ] }, { "cell_type": "code", @@ -131,4 +186,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} \ No newline at end of file +} diff --git a/examples/rubin_dp2/StudyObjectMatches.ipynb b/examples/rubin_dp2/StudyObjectMatches.ipynb index 1c4e8c3..efef92b 100644 --- a/examples/rubin_dp2/StudyObjectMatches.ipynb +++ b/examples/rubin_dp2/StudyObjectMatches.ipynb @@ -3,8 +3,12 @@ { "cell_type": "markdown", "id": "b6097577", - "source": "# Study object matching\n\nExamines the all-catalog object match: SNR by catalog coverage and spatial distribution of fully-matched vs partially-matched objects.", - "metadata": {} + "metadata": {}, + "source": [ + "# Study object matching\n", + "\n", + "Examines the all-catalog object match: SNR by catalog coverage and spatial distribution of fully-matched vs partially-matched objects." + ] }, { "cell_type": "code", @@ -12,21 +16,31 @@ "id": "02118d68-3f5d-4a38-a1f9-e8724bb05b3b", "metadata": {}, "outputs": [], - "source": "# Standard imports\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport tables_io\nimport hpmcm" + "source": [ + "# Standard imports\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import tables_io\n", + "import hpmcm" + ] }, { "cell_type": "markdown", "id": "6baed9c7", - "source": "## Configuration", - "metadata": {} + "metadata": {}, + "source": [ + "## Configuration" + ] }, { "cell_type": "code", + "execution_count": null, "id": "2a49b375", - "source": "tract = 8006", "metadata": {}, - "execution_count": null, - "outputs": [] + "outputs": [], + "source": [ + "tract = 8006" + ] }, { "cell_type": "code", @@ -34,7 +48,23 @@ "id": "3b59e0d0-ab90-4f0e-8212-0b7fcc757675", "metadata": {}, "outputs": [], - "source": "# Input catalog files and IDs for the object meta-match\nall_object = tables_io.read(f'{tract}/object_reduced_{tract}.parq')\nall_anacal = tables_io.read(f'{tract}/anacal_reduced_{tract}.parq')\nall_desc_md = tables_io.read(f'{tract}/desc_md_match_object_stats.pq')\nall_rubin_md = tables_io.read(f'{tract}/rubin_md_match_object_stats.pq')\nobj_stats = tables_io.read(f'{tract}/all_match_object_stats.pq')\nobj_assoc = tables_io.read(f'{tract}/all_match_object_assoc.pq')\n\nall_files = [\n f'{tract}/object_reduced_{tract}.parq',\n f'{tract}/anacal_reduced_{tract}.parq',\n f'{tract}/desc_md_match_object_stats.pq',\n f'{tract}/rubin_md_match_object_stats.pq',\n]\nall_ids = list(range(len(all_files)))\n" + "source": [ + "# Input catalog files and IDs for the object meta-match\n", + "all_object = tables_io.read(f'{tract}/object_reduced_{tract}.parq')\n", + "all_anacal = tables_io.read(f'{tract}/anacal_reduced_{tract}.parq')\n", + "all_desc_md = tables_io.read(f'{tract}/desc_md_match_object_stats.pq')\n", + "all_rubin_md = tables_io.read(f'{tract}/rubin_md_match_object_stats.pq')\n", + "obj_stats = tables_io.read(f'{tract}/all_match_object_stats.pq')\n", + "obj_assoc = tables_io.read(f'{tract}/all_match_object_assoc.pq')\n", + "\n", + "all_files = [\n", + " f'{tract}/object_reduced_{tract}.parq',\n", + " f'{tract}/anacal_reduced_{tract}.parq',\n", + " f'{tract}/desc_md_match_object_stats.pq',\n", + " f'{tract}/rubin_md_match_object_stats.pq',\n", + "]\n", + "all_ids = list(range(len(all_files)))\n" + ] }, { "cell_type": "code", @@ -42,7 +72,15 @@ "id": "4f283920-b521-414b-9f3c-f502ebc21df1", "metadata": {}, "outputs": [], - "source": "# Build wide joined table\njoined_shear_obj = hpmcm.output_tables.buildJoinedObjectTable(\n obj_stats,\n obj_assoc,\n all_files,\n all_ids,\n)" + "source": [ + "# Build wide joined table\n", + "joined_shear_obj = hpmcm.output_tables.buildJoinedObjectTable(\n", + " obj_stats,\n", + " obj_assoc,\n", + " all_files,\n", + " all_ids,\n", + ")" + ] }, { "cell_type": "code", @@ -50,7 +88,16 @@ "id": "e353574d-9e59-4504-b2af-d60656e3e53a", "metadata": {}, "outputs": [], - "source": "# SNR distributions split by catalog detection bits\nbins = np.logspace(0, 5, 51)\nfor i in range(4):\n mask_bit = 1 << i\n mask = np.bitwise_and(joined_shear_obj.catalog_mask, mask_bit) == mask_bit\n _ = plt.hist(joined_shear_obj[mask].snr, bins=bins, alpha=0.1)\n_ = plt.xscale('log')\n_ = plt.yscale('log')" + "source": [ + "# SNR distributions split by catalog detection bits\n", + "bins = np.logspace(0, 5, 51)\n", + "for i in range(4):\n", + " mask_bit = 1 << i\n", + " mask = np.bitwise_and(joined_shear_obj.catalog_mask, mask_bit) == mask_bit\n", + " _ = plt.hist(joined_shear_obj[mask].snr, bins=bins, alpha=0.1)\n", + "_ = plt.xscale('log')\n", + "_ = plt.yscale('log')" + ] }, { "cell_type": "code", @@ -69,7 +116,13 @@ "id": "d88f0341-3a9c-454b-bb43-34bb91439742", "metadata": {}, "outputs": [], - "source": "# Objects detected in all 4 catalogs have catalog_mask == 15\ngood_mask = joined_shear_obj.catalog_mask == 15\nbad_mask = joined_shear_obj.catalog_mask != 15\ngood_sub = joined_shear_obj[good_mask]\nbad_sub = joined_shear_obj[bad_mask]" + "source": [ + "# Objects detected in all 4 catalogs have catalog_mask == 15\n", + "good_mask = joined_shear_obj.catalog_mask == 15\n", + "bad_mask = joined_shear_obj.catalog_mask != 15\n", + "good_sub = joined_shear_obj[good_mask]\n", + "bad_sub = joined_shear_obj[bad_mask]" + ] }, { "cell_type": "code", @@ -130,4 +183,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} \ No newline at end of file +} diff --git a/examples/rubin_dp2/StudyRubinShearMatch.ipynb b/examples/rubin_dp2/StudyRubinShearMatch.ipynb index c54695d..7fc07df 100644 --- a/examples/rubin_dp2/StudyRubinShearMatch.ipynb +++ b/examples/rubin_dp2/StudyRubinShearMatch.ipynb @@ -3,8 +3,12 @@ { "cell_type": "markdown", "id": "d0be3384", - "source": "# Study Rubin shear matching\n\nAnalyzes the Rubin metadetect shear match: SNR distributions, catalog coverage, and spatial distributions of good vs bad matches.", - "metadata": {} + "metadata": {}, + "source": [ + "# Study Rubin shear matching\n", + "\n", + "Analyzes the Rubin metadetect shear match: SNR distributions, catalog coverage, and spatial distributions of good vs bad matches." + ] }, { "cell_type": "code", @@ -12,21 +16,31 @@ "id": "3d64fa81-9bea-47c3-940b-820a464699d0", "metadata": {}, "outputs": [], - "source": "# Standard imports\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport tables_io\nimport hpmcm" + "source": [ + "# Standard imports\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import tables_io\n", + "import hpmcm" + ] }, { "cell_type": "markdown", "id": "ff0976d8", - "source": "## Configuration", - "metadata": {} + "metadata": {}, + "source": [ + "## Configuration" + ] }, { "cell_type": "code", + "execution_count": null, "id": "0384f8ef", - "source": "tract = 8006", "metadata": {}, - "execution_count": null, - "outputs": [] + "outputs": [], + "source": [ + "tract = 8006" + ] }, { "cell_type": "code", @@ -34,7 +48,25 @@ "id": "6d951102-a03d-4065-b7f2-7f51bc68158c", "metadata": {}, "outputs": [], - "source": "# Load Rubin shear match output tables and input catalog files\nshear_rubin_in = tables_io.read(f'{tract}/shear_object_{tract}.parq')\nshear_rubin_ns = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_ns.parq')\nshear_rubin_1m = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_1m.parq')\nshear_rubin_1p = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_1p.parq')\nshear_rubin_2m = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_2m.parq')\nshear_rubin_2p = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_2p.parq') \nshear_rubin_obj_stats = tables_io.read(f'{tract}/rubin_md_match_object_stats.pq')\nshear_rubin_obj_assoc = tables_io.read(f'{tract}/rubin_md_match_object_assoc.pq')\ninput_files = [\n f'{tract}/shear_object_reduced_uncleaned_{tract}_ns.parq',\n f'{tract}/shear_object_reduced_uncleaned_{tract}_1m.parq',\n f'{tract}/shear_object_reduced_uncleaned_{tract}_1p.parq',\n f'{tract}/shear_object_reduced_uncleaned_{tract}_2m.parq',\n f'{tract}/shear_object_reduced_uncleaned_{tract}_2p.parq',\n] \ncatalog_ids = [0, 1, 2, 3, 4]" + "source": [ + "# Load Rubin shear match output tables and input catalog files\n", + "shear_rubin_in = tables_io.read(f'{tract}/shear_object_{tract}.parq')\n", + "shear_rubin_ns = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_ns.parq')\n", + "shear_rubin_1m = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_1m.parq')\n", + "shear_rubin_1p = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_1p.parq')\n", + "shear_rubin_2m = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_2m.parq')\n", + "shear_rubin_2p = tables_io.read(f'{tract}/shear_object_reduced_uncleaned_{tract}_2p.parq') \n", + "shear_rubin_obj_stats = tables_io.read(f'{tract}/rubin_md_match_object_stats.pq')\n", + "shear_rubin_obj_assoc = tables_io.read(f'{tract}/rubin_md_match_object_assoc.pq')\n", + "input_files = [\n", + " f'{tract}/shear_object_reduced_uncleaned_{tract}_ns.parq',\n", + " f'{tract}/shear_object_reduced_uncleaned_{tract}_1m.parq',\n", + " f'{tract}/shear_object_reduced_uncleaned_{tract}_1p.parq',\n", + " f'{tract}/shear_object_reduced_uncleaned_{tract}_2m.parq',\n", + " f'{tract}/shear_object_reduced_uncleaned_{tract}_2p.parq',\n", + "] \n", + "catalog_ids = [0, 1, 2, 3, 4]" + ] }, { "cell_type": "code", @@ -42,7 +74,15 @@ "id": "fa3223ca-d99d-4bf9-b164-ab5aa403ab1d", "metadata": {}, "outputs": [], - "source": "# Build wide joined table\njoined_shear_obj = hpmcm.output_tables.buildJoinedObjectTable(\n shear_rubin_obj_stats,\n shear_rubin_obj_assoc,\n input_files,\n catalog_ids,\n)" + "source": [ + "# Build wide joined table\n", + "joined_shear_obj = hpmcm.output_tables.buildJoinedObjectTable(\n", + " shear_rubin_obj_stats,\n", + " shear_rubin_obj_assoc,\n", + " input_files,\n", + " catalog_ids,\n", + ")" + ] }, { "cell_type": "code", @@ -50,7 +90,16 @@ "id": "029dccce-6b1e-4df1-98ae-4cd93663271d", "metadata": {}, "outputs": [], - "source": "# SNR histograms by catalog detection bits\nbins = np.logspace(0, 5, 51)\nfor i in range(5):\n mask_bit = 1 << i\n mask = np.bitwise_and(joined_shear_obj.catalog_mask, mask_bit) == mask_bit\n _ = plt.hist(joined_shear_obj[mask].snr, bins=bins, alpha=0.1)\n_ = plt.xscale('log')\n_ = plt.yscale('log')" + "source": [ + "# SNR histograms by catalog detection bits\n", + "bins = np.logspace(0, 5, 51)\n", + "for i in range(5):\n", + " mask_bit = 1 << i\n", + " mask = np.bitwise_and(joined_shear_obj.catalog_mask, mask_bit) == mask_bit\n", + " _ = plt.hist(joined_shear_obj[mask].snr, bins=bins, alpha=0.1)\n", + "_ = plt.xscale('log')\n", + "_ = plt.yscale('log')" + ] }, { "cell_type": "code", @@ -69,7 +118,13 @@ "id": "b6f63755-5889-49d2-9df6-be45422f07d3", "metadata": {}, "outputs": [], - "source": "# Objects detected in all 5 Rubin catalogs have catalog_mask == 31\ngood_mask = joined_shear_obj.catalog_mask == 31\nbad_mask = joined_shear_obj.catalog_mask != 31\ngood_sub = joined_shear_obj[good_mask]\nbad_sub = joined_shear_obj[bad_mask]\n" + "source": [ + "# Objects detected in all 5 Rubin catalogs have catalog_mask == 31\n", + "good_mask = joined_shear_obj.catalog_mask == 31\n", + "bad_mask = joined_shear_obj.catalog_mask != 31\n", + "good_sub = joined_shear_obj[good_mask]\n", + "bad_sub = joined_shear_obj[bad_mask]\n" + ] }, { "cell_type": "code", @@ -145,4 +200,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} \ No newline at end of file +} diff --git a/src/hpmcm/__init__.py b/src/hpmcm/__init__.py index a08b5bd..b56bd33 100644 --- a/src/hpmcm/__init__.py +++ b/src/hpmcm/__init__.py @@ -67,12 +67,16 @@ from .object import ObjectData, ShearObjectData from .output_tables import ( ClusterAssocTable, + ClusterShearTable, ClusterStatsTable, ObjectAssocTable, + ObjectShearTable, ObjectStatsTable, ShearTable, buildJoinedClusterTable, buildJoinedObjectTable, + computeColumnStats, + reduceJoinedTable, ) from .shear_data import ShearData from .shear_match import ShearMatch @@ -83,14 +87,18 @@ "CellData", "ClusterAssocTable", "ClusterData", + "ClusterShearTable", "ClusterStatsTable", "buildJoinedClusterTable", "buildJoinedObjectTable", + "computeColumnStats", + "reduceJoinedTable", "Footprint", "FootprintSet", "Match", "ObjectAssocTable", "ObjectData", + "ObjectShearTable", "ObjectStatsTable", "ShearCellData", "ShearClusterData", diff --git a/src/hpmcm/output_tables.py b/src/hpmcm/output_tables.py index 4fd6d8a..0b9d6a3 100644 --- a/src/hpmcm/output_tables.py +++ b/src/hpmcm/output_tables.py @@ -2,6 +2,7 @@ from __future__ import annotations +from pathlib import Path from typing import TYPE_CHECKING import numpy as np @@ -106,6 +107,8 @@ class ObjectStatsTable(TableInterface): cell_idx=TableColumnInfo(int, "Index of associated cell"), has_ref_cat=TableColumnInfo(bool, "Has source from the reference catalog"), catalog_mask=TableColumnInfo(int, "Mask of which catalogs are in object"), + n_sib_objects=TableColumnInfo(int, "Number of objects in parent cluster"), + n_sib_srcs=TableColumnInfo(int, "Total sources across all objects in parent cluster"), ) @staticmethod @@ -134,6 +137,8 @@ def buildFromCellData(cell_data: CellData) -> ObjectStatsTable: snr_rms = np.zeros((n_obj), dtype=float) has_ref_cat = np.zeros((n_obj), dtype=bool) catalog_mask = np.zeros((n_obj), dtype=int) + n_sib_objects = np.zeros((n_obj), dtype=int) + n_sib_srcs = np.zeros((n_obj), dtype=int) for idx, obj in enumerate(cell_data.object_dict.values()): cluster_ids[idx] = obj.parent_cluster.i_cluster @@ -147,6 +152,9 @@ def buildFromCellData(cell_data: CellData) -> ObjectStatsTable: snr_rms[idx] = obj.snr_rms has_ref_cat[idx] = obj.hasRefCatalog() catalog_mask[idx] = obj.catalogMask(catalog_id_map) + sibs = obj.parent_cluster.objects + n_sib_objects[idx] = len(sibs) + n_sib_srcs[idx] = sum(o.n_src for o in sibs) ra, dec = cell_data.getRaDec(x_cents, y_cents) dist_rms *= cell_data.matcher.pixToArcsec() @@ -170,6 +178,8 @@ def buildFromCellData(cell_data: CellData) -> ObjectStatsTable: cell_idx=np.repeat(cell_data.idx, len(dist_rms)).astype(int), has_ref_cat=has_ref_cat, catalog_mask=catalog_mask, + n_sib_objects=n_sib_objects, + n_sib_srcs=n_sib_srcs, ) @@ -339,7 +349,7 @@ def buildFromCellData(cell_data: CellData) -> ClusterStatsTable: class ShearTable(TableInterface): - """Interface of table with shear information""" + """Base interface of table with shear information (no id columns).""" _schema = TableInterface._schema.copy() _schema["good"] = TableColumnInfo(bool, "Has unique match") @@ -358,7 +368,7 @@ class ShearTable(TableInterface): ) @classmethod - def buildObjectShearStats(cls, cell_data: CellData) -> ShearTable: + def buildObjectShearStats(cls, cell_data: CellData) -> ObjectShearTable: """Create shear stats table for objects in a cell Parameters @@ -368,19 +378,20 @@ def buildObjectShearStats(cls, cell_data: CellData) -> ShearTable: Returns ------- - Shear stats table + Object shear stats table (includes ``object_id`` and ``cluster_id``) """ n_obj = cell_data.n_objects - out_dict = ShearTable.emtpyNumpyDict(n_obj) + out_dict = ObjectShearTable.emtpyNumpyDict(n_obj) for idx, obj in enumerate(cell_data.object_dict.values()): assert isinstance(obj, ShearObjectData) - obj_stats = obj.shearStats() - for key, val in obj_stats.items(): + out_dict["object_id"][idx] = obj.object_id + out_dict["cluster_id"][idx] = obj.parent_cluster.i_cluster + for key, val in obj.shearStats().items(): out_dict[key][idx] = val - return ShearTable(**out_dict) + return ObjectShearTable(**out_dict) @classmethod - def buildClusterShearStats(cls, cell_data: CellData) -> ShearTable: + def buildClusterShearStats(cls, cell_data: CellData) -> ClusterShearTable: """Create shear stats table for clusters in a cell Parameters @@ -390,16 +401,35 @@ def buildClusterShearStats(cls, cell_data: CellData) -> ShearTable: Returns ------- - Shear stats table + Cluster shear stats table (includes ``cluster_id``) """ n_clusters = cell_data.n_clusters - out_dict = ShearTable.emtpyNumpyDict(n_clusters) + out_dict = ClusterShearTable.emtpyNumpyDict(n_clusters) for idx, clus in enumerate(cell_data.cluster_dict.values()): assert isinstance(clus, ShearClusterData) - clus_stats = clus.shearStats() - for key, val in clus_stats.items(): + out_dict["cluster_id"][idx] = clus.i_cluster + for key, val in clus.shearStats().items(): out_dict[key][idx] = val - return ShearTable(**out_dict) + return ClusterShearTable(**out_dict) + + +class ObjectShearTable(ShearTable): + """Shear stats table for objects — adds ``object_id`` and ``cluster_id``.""" + + _schema = { + "object_id": TableColumnInfo(int, "Unique Object ID"), + "cluster_id": TableColumnInfo(int, "Parent Cluster Unique ID"), + **ShearTable._schema, + } + + +class ClusterShearTable(ShearTable): + """Shear stats table for clusters — adds ``cluster_id``.""" + + _schema = { + "cluster_id": TableColumnInfo(int, "Unique Cluster ID"), + **ShearTable._schema, + } SourceColsType = list[str] | dict[int, list[str]] | None @@ -429,12 +459,14 @@ def buildJoinedObjectTable( input_files: list[str], catalog_ids: list[int], source_cols: SourceColsType = None, + object_shear: str | Path | None = None, ) -> pandas.DataFrame: """Build a wide joined table from an ObjectStatsTable and source catalogs. Produces one row per object: all stats columns from ``object_stats`` plus source-level data from each catalog joined in as additional columns, - suffixed by ``_{catalog_id}``. + suffixed by ``_{catalog_id}``. Optionally joins shear statistics from an + ``ObjectShearTable`` parquet file. Rows are matched using ``object_assoc.object_id == object_stats.object_id`` and ``object_assoc.source_id == input_catalog.id``. @@ -459,6 +491,11 @@ def buildJoinedObjectTable( - ``dict[int, list[str]]``: per-catalog column lists, keyed by ``catalog_id``. Catalogs absent from the dict get no columns (only ``distance`` is added for them). + object_shear: + Optional path to an ``ObjectShearTable`` parquet file. When provided, + shear statistics are joined on ``object_id``. The ``cluster_id`` + column is dropped before joining because it is already present in + ``object_stats``. Returns ------- @@ -469,6 +506,10 @@ def buildJoinedObjectTable( catalog_file_map = dict(zip(catalog_ids, input_files)) base = object_stats.set_index("object_id") + if object_shear is not None: + shear_df = pandas.read_parquet(object_shear).drop(columns=["cluster_id"], errors="ignore") + base = base.join(shear_df.set_index("object_id"), how="left") + for cat_id, f_name in catalog_file_map.items(): mask = object_assoc["catalog_id"] == cat_id if not mask.any(): @@ -499,17 +540,20 @@ def buildJoinedObjectTable( def buildJoinedClusterTable( + cluster_stats: pandas.DataFrame, cluster_assoc: pandas.DataFrame, input_files: list[str], catalog_ids: list[int], source_cols: SourceColsType = None, + cluster_shear: str | Path | None = None, ) -> pandas.DataFrame: """Build a wide joined table from a ClusterStatsTable and source catalogs. Produces one row per cluster: all stats columns from ``cluster_stats`` plus source-level data from each catalog joined in as additional columns, - suffixed by ``_{catalog_id}``. + suffixed by ``_{catalog_id}``. Optionally joins shear statistics from a + ``ClusterShearTable`` parquet file. Parameters ---------- @@ -531,6 +575,9 @@ def buildJoinedClusterTable( - ``dict[int, list[str]]``: per-catalog column lists, keyed by ``catalog_id``. Catalogs absent from the dict get no columns (only ``distance`` is added for them). + cluster_shear: + Optional path to a ``ClusterShearTable`` parquet file. When provided, + shear statistics are joined on ``cluster_id``. Returns ------- @@ -541,6 +588,10 @@ def buildJoinedClusterTable( catalog_file_map = dict(zip(catalog_ids, input_files)) base = cluster_stats.set_index("cluster_id") + if cluster_shear is not None: + shear_df = pandas.read_parquet(cluster_shear) + base = base.join(shear_df.set_index("cluster_id"), how="left") + for cat_id, f_name in catalog_file_map.items(): mask = cluster_assoc["catalog_id"] == cat_id if not mask.any(): @@ -558,3 +609,142 @@ def buildJoinedClusterTable( base = base.join(rows, how="left") return base.reset_index() + + +def computeColumnStats( + df: pandas.DataFrame, + col_prefix: str, + catalog_ids: list[int], +) -> pandas.DataFrame: + """Compute per-row mean and std across per-catalog columns. + + For each row in ``df``, collects the values of all columns named + ``{col_prefix}_{cat_id}`` for each ``cat_id`` in ``catalog_ids`` that + exists in ``df``, then computes the mean and sample standard deviation + while ignoring NaN values (missing matches). + + Parameters + ---------- + df: + DataFrame produced by :func:`buildJoinedObjectTable` or + :func:`buildJoinedClusterTable`. + col_prefix: + Prefix shared by the per-catalog columns to aggregate, e.g. ``"flux"`` + to operate on ``flux_0``, ``flux_1``, ... + catalog_ids: + Catalog IDs whose columns to include. Any catalog ID whose column + ``{col_prefix}_{cat_id}`` is absent from ``df`` is silently skipped. + + Returns + ------- + DataFrame with the same index as ``df`` and two columns: + + * ``{col_prefix}_mean`` -- row-wise mean of non-NaN values. + * ``{col_prefix}_std`` -- row-wise sample std (ddof=1); NaN when fewer + than two non-NaN values are available. + """ + cols = [ + f"{col_prefix}_{cat_id}" + for cat_id in catalog_ids + if f"{col_prefix}_{cat_id}" in df.columns + ] + if not cols: + return pandas.DataFrame( + { + f"{col_prefix}_mean": np.full(len(df), np.nan), + f"{col_prefix}_std": np.full(len(df), np.nan), + }, + index=df.index, + ) + + arr = df[cols].to_numpy(dtype=float) + n = np.sum(~np.isnan(arr), axis=1) + + mean = np.full(len(df), np.nan) + std = np.full(len(df), np.nan) + + has_any = n >= 1 + has_two = n >= 2 + if has_any.any(): + mean[has_any] = np.nanmean(arr[has_any], axis=1) + if has_two.any(): + std[has_two] = np.nanstd(arr[has_two], axis=1, ddof=1) + + return pandas.DataFrame( + { + f"{col_prefix}_mean": mean, + f"{col_prefix}_std": std, + }, + index=df.index, + ) + + +def reduceJoinedTable( + df: pandas.DataFrame, + catalog_ids: list[int], + keep_cols: list[str] | None = None, + drop_cols: list[str] | None = None, + stats_cols: list[str] | None = None, +) -> pandas.DataFrame: + """Reduce a joined table by selecting columns and aggregating per-catalog ones. + + Applies three independent operations in sequence: + + 1. **Column selection** — keep only the columns in ``keep_cols``; or, if + ``keep_cols`` is ``None``, start with all columns and remove any in + ``drop_cols``. (``drop_cols`` is ignored when ``keep_cols`` is given.) + 2. **Per-catalog column removal** — for each prefix in ``stats_cols``, + the individual per-catalog columns ``{prefix}_{cat_id}`` are removed + from the selected set. + 3. **Stats aggregation** — for each prefix in ``stats_cols``, + :func:`computeColumnStats` is called on the *original* ``df`` and the + three summary columns (``{prefix}_mean``, ``{prefix}_std``, + ``{prefix}_n``) are appended to the result. + + Parameters + ---------- + df: + DataFrame produced by :func:`buildJoinedObjectTable` or + :func:`buildJoinedClusterTable`. + catalog_ids: + Catalog IDs used to identify per-catalog columns and to pass to + :func:`computeColumnStats`. + keep_cols: + Columns to retain as-is. ``None`` keeps all columns (subject to + ``drop_cols``). + drop_cols: + Columns to remove. Only used when ``keep_cols`` is ``None``; silently + ignored otherwise. + stats_cols: + Column-name prefixes whose per-catalog variants should be replaced by + aggregated statistics. Each prefix ``p`` causes ``p_0``, ``p_1``, … + **and any column named exactly** ``p`` to be dropped, and ``p_mean``, + ``p_std``, ``p_n`` to be added. + + Returns + ------- + Reduced DataFrame with the same row order as ``df``. + """ + # 1. Build the base column list + if keep_cols is not None: + base_cols = [c for c in keep_cols if c in df.columns] + else: + base_cols = list(df.columns) + if drop_cols: + excluded = set(drop_cols) + base_cols = [c for c in base_cols if c not in excluded] + + # 2. Remove columns that will be replaced by stats: + # both the bare prefix names and all per-catalog variants {prefix}_{cat_id} + if stats_cols: + to_drop = set(stats_cols) | {f"{p}_{cid}" for p in stats_cols for cid in catalog_ids} + base_cols = [c for c in base_cols if c not in to_drop] + + result = df[base_cols].copy() + + # 3. Append aggregated stats columns + if stats_cols: + stats_frames = [computeColumnStats(df, p, catalog_ids) for p in stats_cols] + result = pandas.concat([result, *stats_frames], axis=1) + + return result diff --git a/tests/test_output_tables_unit.py b/tests/test_output_tables_unit.py index c701c2b..509576a 100644 --- a/tests/test_output_tables_unit.py +++ b/tests/test_output_tables_unit.py @@ -12,6 +12,8 @@ _resolve_cols, buildJoinedClusterTable, buildJoinedObjectTable, + computeColumnStats, + reduceJoinedTable, ) @@ -37,6 +39,7 @@ def object_tables(): object_stats = pd.DataFrame( { "object_id": [100, 101, 102], + "cluster_id": [200, 200, 201], "ra": [0.1, 0.2, 0.3], "dec": [0.1, 0.2, 0.3], "catalog_mask": [3, 1, 2], @@ -228,3 +231,244 @@ def test_buildJoinedClusterTable_preserves_all_stats_rows(source_parquets, clust result = buildJoinedClusterTable(cluster_stats, cluster_assoc, input_files, catalog_ids) assert set(result.cluster_id) == set(cluster_stats.cluster_id) + + +# --------------------------------------------------------------------------- +# object_shear / cluster_shear optional join +# --------------------------------------------------------------------------- + + +@pytest.fixture +def shear_parquets(tmp_path): + """Write minimal shear parquets matching the object/cluster fixtures.""" + object_shear = pd.DataFrame( + { + "object_id": [100, 101, 102], + "cluster_id": [200, 200, 201], # should be dropped before join + "good": [True, True, False], + "g_1_ns": [0.01, 0.02, 0.03], + } + ) + cluster_shear = pd.DataFrame( + { + "cluster_id": [200, 201], + "good": [True, False], + "g_1_ns": [0.015, 0.025], + } + ) + obj_f = str(tmp_path / "object_shear.parquet") + clust_f = str(tmp_path / "cluster_shear.parquet") + object_shear.to_parquet(obj_f) + cluster_shear.to_parquet(clust_f) + return obj_f, clust_f, object_shear, cluster_shear + + +def test_buildJoinedObjectTable_with_shear(source_parquets, object_tables, shear_parquets): + input_files, catalog_ids, _, _ = source_parquets + object_stats, object_assoc = object_tables + obj_shear_f, _, object_shear, _ = shear_parquets + + result = buildJoinedObjectTable( + object_stats, object_assoc, input_files, catalog_ids, + object_shear=obj_shear_f, + ) + + assert "g_1_ns" in result.columns + assert "good" in result.columns + # cluster_id from shear table is dropped; only one copy (from object_stats) remains + assert result.columns.tolist().count("cluster_id") == 1 + assert result.loc[result.object_id == 100, "cluster_id"].iloc[0] == 200 + # values match the shear fixture + row100 = result.loc[result.object_id == 100].iloc[0] + assert row100["g_1_ns"] == pytest.approx(0.01) + row102 = result.loc[result.object_id == 102].iloc[0] + assert row102["good"] is False or row102["good"] == False # noqa: E712 + + +def test_buildJoinedClusterTable_with_shear(source_parquets, cluster_tables, shear_parquets): + input_files, catalog_ids, _, _ = source_parquets + cluster_stats, cluster_assoc = cluster_tables + _, clust_shear_f, _, cluster_shear = shear_parquets + + result = buildJoinedClusterTable( + cluster_stats, cluster_assoc, input_files, catalog_ids, + cluster_shear=clust_shear_f, + ) + + assert "g_1_ns" in result.columns + assert "good" in result.columns + row200 = result.loc[result.cluster_id == 200].iloc[0] + assert row200["g_1_ns"] == pytest.approx(0.015) + row201 = result.loc[result.cluster_id == 201].iloc[0] + assert row201["good"] is False or row201["good"] == False # noqa: E712 + + +# --------------------------------------------------------------------------- +# computeColumnStats +# --------------------------------------------------------------------------- + + +@pytest.fixture +def joined_df(): + """Minimal joined DataFrame with flux columns from three catalogs.""" + return pd.DataFrame( + { + "object_id": [1, 2, 3, 4], + "flux_0": [1.0, 2.0, np.nan, np.nan], + "flux_1": [3.0, np.nan, 5.0, np.nan], + "flux_2": [5.0, 4.0, 7.0, np.nan], + } + ) + + +def test_computeColumnStats_basic(joined_df): + result = computeColumnStats(joined_df, "flux", [0, 1, 2]) + + assert list(result.columns) == ["flux_mean", "flux_std"] + assert len(result) == len(joined_df) + + # row 0: values [1, 3, 5] → mean=3, std=2 + assert result["flux_mean"].iloc[0] == pytest.approx(3.0) + assert result["flux_std"].iloc[0] == pytest.approx(2.0) + + # row 1: values [2, NaN, 4] → mean=3, std≈1.414 + assert result["flux_mean"].iloc[1] == pytest.approx(3.0) + assert result["flux_std"].iloc[1] == pytest.approx(np.sqrt(2)) + + # row 2: values [NaN, 5, 7] → mean=6 + assert result["flux_mean"].iloc[2] == pytest.approx(6.0) + + +def test_computeColumnStats_single_value_has_nan_std(joined_df): + """With only one non-NaN value, sample std (ddof=1) is NaN.""" + result = computeColumnStats(joined_df, "flux", [1]) + assert np.isnan(result["flux_std"].iloc[2]) + + +def test_computeColumnStats_all_nan_row(joined_df): + """Row with all-NaN values gets NaN mean and std.""" + result = computeColumnStats(joined_df, "flux", [0, 1, 2]) + assert np.isnan(result["flux_mean"].iloc[3]) + assert np.isnan(result["flux_std"].iloc[3]) + + +def test_computeColumnStats_missing_catalog_skipped(joined_df): + """Catalog IDs with no matching column are silently skipped.""" + result = computeColumnStats(joined_df, "flux", [0, 99]) # 99 doesn't exist + assert result["flux_mean"].iloc[0] == pytest.approx(1.0) + + +def test_computeColumnStats_no_matching_columns(joined_df): + """When no column matches the prefix, returns NaN columns of correct length.""" + result = computeColumnStats(joined_df, "nonexistent", [0, 1, 2]) + assert len(result) == len(joined_df) + assert np.all(np.isnan(result["nonexistent_mean"])) + + +# --------------------------------------------------------------------------- +# reduceJoinedTable +# --------------------------------------------------------------------------- + + +@pytest.fixture +def wide_df(): + """Joined DataFrame with several per-catalog columns and metadata columns.""" + return pd.DataFrame( + { + "object_id": [1, 2, 3], + "ra": [0.1, 0.2, 0.3], + "dec": [0.1, 0.2, 0.3], + "cluster_id": [10, 10, 20], + "flux_0": [1.0, 2.0, np.nan], + "flux_1": [3.0, np.nan, 5.0], + "flux_2": [5.0, 4.0, 7.0], + "size_0": [0.5, 0.6, 0.7], + "size_1": [0.4, np.nan, 0.8], + "size_2": [0.6, 0.5, np.nan], + } + ) + + +def test_reduceJoinedTable_keep_cols(wide_df): + """keep_cols selects only the named columns.""" + result = reduceJoinedTable(wide_df, [0, 1, 2], keep_cols=["object_id", "ra"]) + assert list(result.columns) == ["object_id", "ra"] + assert len(result) == len(wide_df) + + +def test_reduceJoinedTable_drop_cols(wide_df): + """drop_cols removes columns when keep_cols is None.""" + result = reduceJoinedTable(wide_df, [0, 1, 2], drop_cols=["cluster_id", "dec"]) + assert "cluster_id" not in result.columns + assert "dec" not in result.columns + assert "object_id" in result.columns + + +def test_reduceJoinedTable_drop_cols_ignored_with_keep(wide_df): + """drop_cols is ignored when keep_cols is provided.""" + result = reduceJoinedTable( + wide_df, [0, 1, 2], + keep_cols=["object_id", "cluster_id"], + drop_cols=["cluster_id"], + ) + assert "cluster_id" in result.columns + + +def test_reduceJoinedTable_stats_cols(wide_df): + """stats_cols replaces per-catalog columns with mean and std (no _n).""" + result = reduceJoinedTable(wide_df, [0, 1, 2], stats_cols=["flux"]) + assert "flux_mean" in result.columns + assert "flux_std" in result.columns + assert "flux_n" not in result.columns + # raw per-catalog columns are gone + for cid in [0, 1, 2]: + assert f"flux_{cid}" not in result.columns + # other columns untouched + assert "size_0" in result.columns + assert "object_id" in result.columns + + +def test_reduceJoinedTable_stats_values(wide_df): + """Stats in the reduced table match direct computeColumnStats output.""" + result = reduceJoinedTable(wide_df, [0, 1, 2], stats_cols=["flux"]) + expected = computeColumnStats(wide_df, "flux", [0, 1, 2]) + np.testing.assert_array_almost_equal(result["flux_mean"], expected["flux_mean"]) + np.testing.assert_array_almost_equal(result["flux_std"], expected["flux_std"]) + assert "flux_n" not in result.columns + + +def test_reduceJoinedTable_multiple_stats_cols(wide_df): + """Multiple stats_cols are all aggregated.""" + result = reduceJoinedTable(wide_df, [0, 1, 2], stats_cols=["flux", "size"]) + for prefix in ["flux", "size"]: + assert f"{prefix}_mean" in result.columns + for cid in [0, 1, 2]: + assert f"{prefix}_{cid}" not in result.columns + + +def test_reduceJoinedTable_keep_and_stats(wide_df): + """keep_cols combined with stats_cols: raw per-catalog cols are still removed.""" + result = reduceJoinedTable( + wide_df, [0, 1, 2], + keep_cols=["object_id", "flux_0", "flux_1", "flux_2"], + stats_cols=["flux"], + ) + assert "flux_mean" in result.columns + for cid in [0, 1, 2]: + assert f"flux_{cid}" not in result.columns + assert "object_id" in result.columns + + +def test_reduceJoinedTable_bare_prefix_col_dropped(wide_df): + """A column named exactly after a stats prefix is also dropped.""" + df_with_bare = wide_df.copy() + df_with_bare["flux"] = 99.0 # bare column with the same name as the prefix + result = reduceJoinedTable(df_with_bare, [0, 1, 2], stats_cols=["flux"]) + assert "flux" not in result.columns + assert "flux_mean" in result.columns + + +def test_reduceJoinedTable_no_reduction(wide_df): + """With no stats_cols and no drops, the DataFrame is returned unchanged.""" + result = reduceJoinedTable(wide_df, [0, 1, 2]) + assert list(result.columns) == list(wide_df.columns)