From b1b31bb98b2b0e18e103615d599eedb28520e884 Mon Sep 17 00:00:00 2001 From: Raphael Shirley Date: Wed, 5 Jun 2024 16:15:51 +0200 Subject: [PATCH 01/14] Prior based on weights --- CMakeLists.txt | 1 + docs/notebooks/Prior.ipynb | 432 +++++++++++++++++++++++++++++++++++++ src/lib/Makefile | 5 +- src/lib/_bindings.cc | 17 ++ src/lib/onesource.cpp | 243 +++++++++++---------- src/lib/onesource.h | 8 +- 6 files changed, 589 insertions(+), 117 deletions(-) create mode 100644 docs/notebooks/Prior.ipynb diff --git a/CMakeLists.txt b/CMakeLists.txt index 3ef880db..81fe7b7d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -46,6 +46,7 @@ set(SOURCES "${SOURCE_DIR}/opa.cpp" "${SOURCE_DIR}/PDF.cpp" "${SOURCE_DIR}/photoz_lib.cpp" + "${SOURCE_DIR}/prior.cpp" "${SOURCE_DIR}/SED.cpp" "${SOURCE_DIR}/SEDLib.cpp" ) diff --git a/docs/notebooks/Prior.ipynb b/docs/notebooks/Prior.ipynb new file mode 100644 index 00000000..e43ba77f --- /dev/null +++ b/docs/notebooks/Prior.ipynb @@ -0,0 +1,432 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "59aa16dc-8c12-4c0c-a75d-9251c7efd3b3", + "metadata": {}, + "source": [ + "# Prior\n", + "\n", + "In this notebook I will demonstrate a basic redshift prior and show how more complex priors might be developed." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ff92daab-96a8-40ba-8d3c-7b8cff3612e2", + "metadata": {}, + "outputs": [], + "source": [ + "import lephare as lp\n", + "import numpy as np\n", + "from matplotlib import pylab as plt\n", + "\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fc8e1f19-837d-4fe2-b1e6-874caf4dd74b", + "metadata": {}, + "outputs": [], + "source": [ + "config = lp.default_cosmos_config.copy()\n", + "config[\"Z_STEP\"] = \"0.1,0.,5.\" # To make notebook very fast to run\n", + "keymap = lp.all_types_to_keymap(config)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c3a37641-afef-4211-89b0-f932e8feb68d", + "metadata": {}, + "outputs": [], + "source": [ + "config" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0b6c33e3-8bf7-4352-b73c-48b89abf1a7e", + "metadata": {}, + "outputs": [], + "source": [ + "# lp.prepare(keymap)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8d6e2bd8-08e8-465e-b0b3-ad84d8259b69", + "metadata": {}, + "outputs": [], + "source": [ + "photz = lp.PhotoZ(keymap)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6d68c664-decb-4357-a847-65f46eae4314", + "metadata": {}, + "outputs": [], + "source": [ + "len(photz.gridz)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "844b1a2f-5a7b-4009-8221-04921e216939", + "metadata": {}, + "outputs": [], + "source": [ + "len(photz.fullLib)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fad4da93-858e-4f15-8dd0-f4b5ab98e89a", + "metadata": {}, + "outputs": [], + "source": [ + "def fulllib_to_weights(fullLib):\n", + " \"\"\"Function which takes the full library of models as an input and returns a prior weight for each.\n", + "\n", + " This is a means to apply a fully generalisable prior which can be called independently for each source.\n", + "\n", + " In this case as a basic example I impose a maximum redshift of 2.\n", + "\n", + " Parameters\n", + " ==========\n", + "\n", + " fullLib : list of lephare.QSOSED, lephare.StarSED, lephare.GalSED\n", + " The full library of models to fit\n", + "\n", + " Returns\n", + " =======\n", + "\n", + " weights : list of floats\n", + " the list of prior weights to apply to each model.\n", + " \"\"\"\n", + " weights = np.array([1 for l in fullLib])\n", + " # Get the redshift using the index\n", + " # Seems odd. Is this really the best way to get the redshift?\n", + " red = np.array([photz.gridz[n - photz.fullLib[n].index_z0] for n in np.arange(len(photz.fullLib))])\n", + "\n", + " weights[red > 2] = 1.0e-9\n", + "\n", + " return weights\n", + "\n", + "\n", + "fulllib_to_weights(photz.fullLib)[:5]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f21e1080-3131-4e7f-bca5-686f4fc0c150", + "metadata": {}, + "outputs": [], + "source": [ + "# Lets look at what we have access to determine the prior weight" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a06ffb78-2701-414e-a1a2-69c8e930c60d", + "metadata": {}, + "outputs": [], + "source": [ + "qso_ex = photz.fullLib[0]\n", + "star_ex = photz.fullLib[900]\n", + "gal_ex = photz.fullLib[3000]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6e4b8afa-c70e-46b4-9efd-40235ac45ad0", + "metadata": {}, + "outputs": [], + "source": [ + "mods = [l.nummod for l in photz.fullLib]\n", + "np.unique(mods)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "61ae1aea-32d2-4036-9bd8-755e2459a441", + "metadata": {}, + "outputs": [], + "source": [ + "photz.fullLib[3000].name, photz.fullLib[7001].name" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2f072770-c20f-4be2-b84d-291441ad7813", + "metadata": {}, + "outputs": [], + "source": [ + "star_ex.index_z0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d846af5c-5de0-4b72-a66f-173c8bbd19f8", + "metadata": {}, + "outputs": [], + "source": [ + "gal_ex.index_z0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5720827c-e3ad-4254-bccd-f3425cf6d391", + "metadata": {}, + "outputs": [], + "source": [ + "red = [photz.gridz[n - photz.fullLib[n].index_z0] for n in np.arange(len(photz.fullLib))]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18256cb5-c9a1-4b14-9bb3-8dd87393c95f", + "metadata": {}, + "outputs": [], + "source": [ + "red[:5]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ecdbfabd-b7b0-4acf-9374-88f9cd17b0be", + "metadata": {}, + "outputs": [], + "source": [ + "np.unique(red)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "37bdaa71-7914-437f-a246-cc9e66ca92c9", + "metadata": {}, + "outputs": [], + "source": [ + "photz.fullLib[3000].ebv, photz.fullLib[4001].ebv" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ac90eabe-dbeb-4314-bc1c-b4f70fa984f3", + "metadata": {}, + "outputs": [], + "source": [ + "photz.fullLib[:5]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "71b56c8d-4a8e-475f-9c9e-b33684469cc1", + "metadata": {}, + "outputs": [], + "source": [ + "# Set the prior function via pybind" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5277626b-e15f-434a-b38d-0f912d1e5e9b", + "metadata": {}, + "outputs": [], + "source": [ + "cat = np.loadtxt(f\"{lp.LEPHAREDIR}/examples/COSMOS.in\")\n", + "id = cat[:, 0]\n", + "fluxes = cat[:, 1:60:2]\n", + "efluxes = cat[:, 2:61:2]\n", + "context = cat[:, 61]\n", + "zspec = cat[:, 62]\n", + "print(\"Check format with context and zspec :\", context, zspec)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8d4cb7c7-2f87-4b3f-9c1f-0b8adb087f62", + "metadata": {}, + "outputs": [], + "source": [ + "w = fulllib_to_weights(photz.fullLib)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8bb4650e-bb43-4b37-8e5f-fb7501883ffe", + "metadata": {}, + "outputs": [], + "source": [ + "srclist = []\n", + "zspec_mask = np.logical_and(zspec > 0.01, zspec < 6)\n", + "for i in np.where(zspec_mask)[0]:\n", + " oneObj = lp.onesource(i, photz.gridz)\n", + " oneObj.readsource(str(id[i]), fluxes[i, :], efluxes[i, :], int(context[i]), zspec[i], \" \")\n", + " oneObj.priorObj.apply_weights = 1\n", + " oneObj.priorObj.weights = w\n", + " photz.prep_data(oneObj)\n", + " srclist.append(oneObj)\n", + "print(\"Sources with a spec-z: \", len(srclist))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "80e3c7e7-5bef-41ef-b19c-6db2b735a772", + "metadata": {}, + "outputs": [], + "source": [ + "# Run the photoz" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "46ac8fe0-7b61-4b61-80d8-cacb53ff7c34", + "metadata": {}, + "outputs": [], + "source": [ + "a0, a1 = photz.run_autoadapt(srclist)\n", + "offsets = \",\".join(np.array(a0).astype(str))\n", + "offsets = \"# Offsets from auto-adapt: \" + offsets + \"\\n\"\n", + "print(offsets)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b03c41ac-d336-4ec7-8a16-72804e608203", + "metadata": {}, + "outputs": [], + "source": [ + "photozlist = []\n", + "for i in range(100):\n", + " oneObj = lp.onesource(i, photz.gridz)\n", + " oneObj.readsource(str(id[i]), fluxes[i, :], efluxes[i, :], int(context[i]), zspec[i], \" \")\n", + " oneObj.priorObj.apply_weights = 1\n", + " oneObj.priorObj.weights = w\n", + " photz.prep_data(oneObj)\n", + " photozlist.append(oneObj)\n", + "print(\"Number of sources to be analysed: \", len(srclist))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3063e28b-0b5a-438d-9dad-f1560e4dea56", + "metadata": {}, + "outputs": [], + "source": [ + "photz.run_photoz(photozlist[:100], a0, a1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e356e632-9a6e-445c-9f55-64aefb3112fd", + "metadata": {}, + "outputs": [], + "source": [ + "t = photz.build_output_tables(photozlist[:100], para_out=None, filename=\"testPrior.fits\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2d288b30-030d-4624-837f-edcd422be307", + "metadata": {}, + "outputs": [], + "source": [ + "t[:5]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "729b795e-2be5-4807-87ce-4a227996780d", + "metadata": {}, + "outputs": [], + "source": [ + "plt.hist(t[\"Z_BEST\"], bins=50)\n", + "plt.xlabel(\"z with z<2 prior\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23b7b424-4924-48d3-b564-551b09248f8e", + "metadata": {}, + "outputs": [], + "source": [ + "photozlist = []\n", + "for i in range(100):\n", + " oneObj = lp.onesource(i, photz.gridz)\n", + " oneObj.readsource(str(id[i]), fluxes[i, :], efluxes[i, :], int(context[i]), zspec[i], \" \")\n", + " # oneObj.priorObj.apply_weights=1\n", + " # oneObj.priorObj.weights=w\n", + " photz.prep_data(oneObj)\n", + " photozlist.append(oneObj)\n", + "print(\"Number of sources to be analysed: \", len(srclist))\n", + "\n", + "photz.run_photoz(photozlist[:100], a0, a1)\n", + "\n", + "t2 = photz.build_output_tables(photozlist[:100], para_out=None, filename=\"testPrior.fits\")\n", + "plt.hist(t2[\"Z_BEST\"], bins=50)\n", + "plt.xlabel(\"z without z<2 prior\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "84085877-f290-436c-aa3a-b886a64d9b9f", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "lincc", + "language": "python", + "name": "lincc" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/lib/Makefile b/src/lib/Makefile index 21e3433e..93a3eb3a 100755 --- a/src/lib/Makefile +++ b/src/lib/Makefile @@ -60,7 +60,7 @@ OBJECTS = \ keyword.o\ flt.o\ SED.o\ - SEDLib.o\ + SEDLib.o\ ext.o\ opa.o\ mag.o\ @@ -68,7 +68,8 @@ OBJECTS = \ cosmology.o\ onesource.o\ PDF.o\ - photoz_lib.o\ + photoz_lib.o\ + prior.o\ all: $(FILES) if [ ! -d "../../bin" ]; then mkdir ../../bin; fi diff --git a/src/lib/_bindings.cc b/src/lib/_bindings.cc index fb393d6b..f7b4ad81 100644 --- a/src/lib/_bindings.cc +++ b/src/lib/_bindings.cc @@ -347,6 +347,8 @@ PYBIND11_MODULE(_lephare, mod) { .def_readwrite("spec", &onesource::spec) .def_readwrite("consiz", &onesource::consiz) .def_readwrite("closest_red", &onesource::closest_red) + // The prior for overwriting the weights + .def_readwrite("priorObj", &onesource::priorObj) .def_readonly("pos", &onesource::pos) .def_readonly("cont", &onesource::cont) .def_readonly("pdfmap", &onesource::pdfmap) @@ -408,6 +410,21 @@ PYBIND11_MODULE(_lephare, mod) { .def_readonly("zsecScale", &onesource::zsecMod) .def_readonly("zsecAge", &onesource::zsecAge); + /******** CLASS prior *********/ + py::class_( + mod, "prior", + "The prior class which allows Python overwriting of the weights") + .def(py::init(), + "Constructor for prior class initialised with a onesource instance") + .def("set_weights_function", &prior::set_weights_function, + "Set a new weights function from Python") + .def_readwrite("apply_nz", &prior::apply_nz, + "int 1 if apply nzlim prior.") + .def_readwrite("apply_weights", &prior::apply_weights, + "int 1 to apply weights.") + .def_readwrite("weights", &prior::weights, + "The weights to apply to the fulllib models."); + py::class_(mod, "PDF") .def(py::init(), py::arg("min"), py::arg("max"), py::arg("size")) diff --git a/src/lib/onesource.cpp b/src/lib/onesource.cpp index a51bb769..8a92ca8c 100644 --- a/src/lib/onesource.cpp +++ b/src/lib/onesource.cpp @@ -17,6 +17,7 @@ #include // print standard file #include #include +#include #include #include @@ -39,6 +40,7 @@ void onesource::readsource(const string &identifier, const vector vals, sab = err_vals; zs = z_spec; str_inp = additional_input; + prior priorObj; // Initialise the prior so that weights can be overwritten. } /* @@ -368,6 +370,13 @@ void onesource::fit(vector &fulllib, const vector> &flux, int number_threads = 1, thread_id = 0; size_t imagm = ab.size(); + // Check weights are correct if used + if (priorObj.apply_weights == 1) { + if (fulllib.size() != priorObj.weights.size()) { + throw std::length_error("weights and fullLib are not the same length."); + } + } + // Do a local minimisation per thread (store chi2 and index) // Catch first the number of threads #ifdef _OPENMP @@ -446,31 +455,38 @@ void onesource::fit(vector &fulllib, const vector> &flux, // Model rejection based on prior // Abs mag rejection - double reds; - int libtype; - reds = (fulllib[il])->red; - libtype = (fulllib[il])->nlib; - // Abs Magnitude from model @ z=0 for rejection for galaxies and AGN - if ((mabsGALprior && libtype == 0) || (mabsAGNprior && libtype == 1)) { - // predicted magnitudes within the babs filter renormalized by the - // scaling - double abs_mag; - abs_mag = (fulllib[il])->mag0 - 2.5 * log10(dmloc[il]); - // specific case when the redshift is 0 - if (reds < 1.e-10) abs_mag = abs_mag - funz0; - // Galaxy rejection - if ((abs_mag <= priorLib[0] || abs_mag >= priorLib[1]) && libtype == 0) - chi2loc[il] = 1.e9; - // AGN rejection - if ((abs_mag <= priorLib[2] || abs_mag >= priorLib[3]) && libtype == 1) - chi2loc[il] = 1.e9; - } - // Prior N(z) - if (bp[0] >= 0 && libtype == 0) { - double pweight = - nzprior((fulllib[il])->luv, (fulllib[il])->lnir, reds, bp); - chi2loc[il] = chi2loc[il] - 2. * log(pweight); - } + // double reds; + // int libtype; + // reds = (fulllib[il])->red; + // libtype = (fulllib[il])->nlib; + // // Abs Magnitude from model @ z=0 for rejection for galaxies and AGN + // if ((mabsGALprior && libtype == 0) || (mabsAGNprior && libtype == 1)) { + // // predicted magnitudes within the babs filter renormalized by the + // // scaling + // double abs_mag; + // abs_mag = (fulllib[il])->mag0 - 2.5 * log10(dmloc[il]); + // // specific case when the redshift is 0 + // if (reds < 1.e-10) abs_mag = abs_mag - funz0; + // // Galaxy rejection + // if ((abs_mag <= priorLib[0] || abs_mag >= priorLib[1]) && libtype == + // 0) + // chi2loc[il] = 1.e9; + // // AGN rejection + // if ((abs_mag <= priorLib[2] || abs_mag >= priorLib[3]) && libtype == + // 1) + // chi2loc[il] = 1.e9; + // } + // // Prior N(z) + // if (bp[0] >= 0 && libtype == 0) { + // double pweight = + // nzprior((fulllib[il])->luv, (fulllib[il])->lnir, reds, bp); + // chi2loc[il] = chi2loc[il] - 2. * log(pweight); + // } + + // General prior + chi2loc[il] = + priorObj.update_chi2(this, chi2loc[il], fulllib[il], il, dmloc[il], + funz0, bp, mabsGALprior, mabsAGNprior); // keep the minimum chi2 if (chi2loc[il] < 0.99e9) { @@ -549,94 +565,95 @@ void onesource::compute_best_fit_physical_quantities(vector &fulllib) { /* Use prior info on N(z) as a function of magnitude (Iband) and of the type. */ -double onesource::nzprior(const double luv, const double lnir, - const double reds, const array bp) { - double val = 1.; - int mod; - double zot = 0., kt = 0., alpt0 = 0.; - - // Define the NUV-R rest-frame color corrected for dust-extinction - double nuvk = -2.5 * (luv - lnir); - - // Apparent magnitude in i band, if possible - double mi = -99.; - if (busnorma[bp[0]] == 1) { - mi = flux2mag(ab[bp[0]]); - } else if (busnorma[bp[1]] == 1) { - mi = flux2mag(ab[bp[1]]); - } else - return val; - - // Bright case I<20, no solution above 1, don't compute the prior - if (mi < 20.) { - if (reds > 1) val = 0; - return val; - } - // Bright case I<22, no solution above 2, but compute the prior - if (mi < 22. && reds > 2) { - val = 0; - return val; - } - - // Set up the parameters to define the redshift distribution - if (nuvk > 4.25) { - // Case E/S0 - // Color UV-K of PHOTO_230506/El_cww.sed.resample.new.resample15.inter - mod = 0; - zot = 0.45181; - kt = 0.13677; - alpt0 = 3.33078; - } else if (nuvk > 3.19 && nuvk < 4.25) { - // Case Sbc - // Color UV-K of PHOTO_230506/Sbc_cww.sed.resample.new.resample8.inter - mod = 1; - zot = 0.16560; - kt = 0.12983; - alpt0 = 1.42815; - } else if (nuvk > 1.9 && nuvk < 3.19) { - // Case Scd - // Color UV-K of PHOTO_230506/Scd_cww.sed.resample.new.resample7.inter - // -19.4878 + 21.1501 - mod = 2; - zot = 0.21072; - kt = 0.14008; - alpt0 = 1.58310; - } else { - // Case Irr - mod = 3; - zot = 0.20418; - kt = 0.13773; - alpt0 = 1.34500; - } - - // P(z|T,m0) - double zmax = zot + kt * (mi - 20); - double pz = pow(reds, alpt0) * exp(-pow((reds / zmax), alpt0)); - - // P(T|m0) - double ktf[4] = {0.47165, 0.30663, 0.12715, -0.34437}; - double ft[4] = {0.43199, 0.07995, 0.31162, 0.21220}; - - // Ratio for each type - double rappSum = - ft[0] * exp(-ktf[0] * (mi - 20.0)) + ft[1] * exp(-ktf[1] * (mi - 20.0)) + - ft[2] * exp(-ktf[2] * (mi - 20.0)) + ft[3] * exp(-ktf[3] * (mi - 20.0)); - double rapp = ft[mod] * exp(-ktf[mod] * (mi - 20.0)); - - // Normalisation of the probability function - // pcal=exp(gammln(1.0/alpt0(mod)+1.0)) - double pcal; - if (mod == 0) pcal = 0.89744; - if (mod == 1) pcal = 0.90868; - if (mod == 2) pcal = 0.89747; - if (mod == 3) pcal = 0.91760; - pcal = pow(zmax, alpt0 + 1) / alpt0 * pcal; - - // Final value - val = pz / pcal * rapp / rappSum; - - return val; -} +// double onesource::nzprior(const double luv, const double lnir, +// const double reds, const array bp) { +// double val = 1.; +// int mod; +// double zot = 0., kt = 0., alpt0 = 0.; + +// // Define the NUV-R rest-frame color corrected for dust-extinction +// double nuvk = -2.5 * (luv - lnir); + +// // Apparent magnitude in i band, if possible +// double mi = -99.; +// if (busnorma[bp[0]] == 1) { +// mi = flux2mag(ab[bp[0]]); +// } else if (busnorma[bp[1]] == 1) { +// mi = flux2mag(ab[bp[1]]); +// } else +// return val; + +// // Bright case I<20, no solution above 1, don't compute the prior +// if (mi < 20.) { +// if (reds > 1) val = 0; +// return val; +// } +// // Bright case I<22, no solution above 2, but compute the prior +// if (mi < 22. && reds > 2) { +// val = 0; +// return val; +// } + +// // Set up the parameters to define the redshift distribution +// if (nuvk > 4.25) { +// // Case E/S0 +// // Color UV-K of PHOTO_230506/El_cww.sed.resample.new.resample15.inter +// mod = 0; +// zot = 0.45181; +// kt = 0.13677; +// alpt0 = 3.33078; +// } else if (nuvk > 3.19 && nuvk < 4.25) { +// // Case Sbc +// // Color UV-K of PHOTO_230506/Sbc_cww.sed.resample.new.resample8.inter +// mod = 1; +// zot = 0.16560; +// kt = 0.12983; +// alpt0 = 1.42815; +// } else if (nuvk > 1.9 && nuvk < 3.19) { +// // Case Scd +// // Color UV-K of PHOTO_230506/Scd_cww.sed.resample.new.resample7.inter +// // -19.4878 + 21.1501 +// mod = 2; +// zot = 0.21072; +// kt = 0.14008; +// alpt0 = 1.58310; +// } else { +// // Case Irr +// mod = 3; +// zot = 0.20418; +// kt = 0.13773; +// alpt0 = 1.34500; +// } + +// // P(z|T,m0) +// double zmax = zot + kt * (mi - 20); +// double pz = pow(reds, alpt0) * exp(-pow((reds / zmax), alpt0)); + +// // P(T|m0) +// double ktf[4] = {0.47165, 0.30663, 0.12715, -0.34437}; +// double ft[4] = {0.43199, 0.07995, 0.31162, 0.21220}; + +// // Ratio for each type +// double rappSum = +// ft[0] * exp(-ktf[0] * (mi - 20.0)) + ft[1] * exp(-ktf[1] * (mi - 20.0)) +// + ft[2] * exp(-ktf[2] * (mi - 20.0)) + ft[3] * exp(-ktf[3] * (mi +// - 20.0)); +// double rapp = ft[mod] * exp(-ktf[mod] * (mi - 20.0)); + +// // Normalisation of the probability function +// // pcal=exp(gammln(1.0/alpt0(mod)+1.0)) +// double pcal; +// if (mod == 0) pcal = 0.89744; +// if (mod == 1) pcal = 0.90868; +// if (mod == 2) pcal = 0.89747; +// if (mod == 3) pcal = 0.91760; +// pcal = pow(zmax, alpt0 + 1) / alpt0 * pcal; + +// // Final value +// val = pz / pcal * rapp / rappSum; + +// return val; +// } /* Check for discrepant band diff --git a/src/lib/onesource.h b/src/lib/onesource.h index b66236cf..8e9e6305 100644 --- a/src/lib/onesource.h +++ b/src/lib/onesource.h @@ -18,6 +18,7 @@ #include "flt.h" // filter class #include "globals.h" #include "opa.h" +#include "prior.h" using namespace std; @@ -60,6 +61,9 @@ class onesource { array priorLib; // Prior with the range in abs mag gal, abs mag AGN + // The general prior + prior priorObj; + vector chibay; vector gridzg, gridLdustIR, gridEbv, gridLIR; PDF PDFebv; @@ -176,8 +180,8 @@ class onesource { const array &bp); void fitIR(vector &fulllib, const vector> &flux, const int imagm, const string fit_frsc, cosmo lcdm); - double nzprior(const double luv, const double lnir, const double reds, - const array bp); + // double nzprior(const double luv, const double lnir, const double reds, + // const array bp); void rm_discrepant(vector &fulllib, const vector> &flux, const vector &valid, const double funz0, const array bp, double thresholdChi2, From 1630a53473f3f68220a099182e994bd1b902334a Mon Sep 17 00:00:00 2001 From: Raphael Shirley Date: Wed, 5 Jun 2024 16:20:36 +0200 Subject: [PATCH 02/14] Added the prior class files --- src/lib/prior.cpp | 179 ++++++++++++++++++++++++++++++++++++++++++++++ src/lib/prior.h | 88 +++++++++++++++++++++++ 2 files changed, 267 insertions(+) create mode 100644 src/lib/prior.cpp create mode 100644 src/lib/prior.h diff --git a/src/lib/prior.cpp b/src/lib/prior.cpp new file mode 100644 index 00000000..da59efa3 --- /dev/null +++ b/src/lib/prior.cpp @@ -0,0 +1,179 @@ +/* + + 24/11/2023 + Implementation of the functions of the prior class. + + This should allow the current absolute magnitude and N(z) priors to be + abstracted outside of the onesource class. + + It should be designed in such a way as to allow overwriting via Python or c++ + with a general prior. + + Multiple priors should be possible. In the first instance we implement + absolute mag and N(z) priors. + +*/ + +#include "prior.h" + +#include +#include + +#include "SED.h" +#include "globals.h" +#include "onesource.h" + +// #include +// namespace py = pybind11; + +using namespace std; + +// The update_chi2 method takes the chi2 and modifies it as a function of the +// SED properties. It must not take source properties which do not change during +// fit. +double prior::update_chi2(onesource* source, double chi2, SED* lib, int il, + double dm, double funz0, const array& bp, + bool mabsGALprior, bool mabsAGNprior) { + // Get all required values from the SED + reds = (lib)->red; + libtype = (lib)->nlib; + luv = (lib)->luv; + lnir = (lib)->lnir; + mag0 = (lib)->mag0; + + // Apply all priors + // nz prior + if (bp[0] >= 0 && libtype == 0 && apply_nz) { + chi2 = nz_prior(source, chi2, luv, lnir, reds, bp); + } + // Abs Magnitude from model @ z=0 for rejection for galaxies and AGN + if ((mabsGALprior && libtype == 0) || (mabsAGNprior && libtype == 1)) { + chi2 = absmag_prior(source, chi2, reds, libtype, dm, mag0, funz0); + } + // General weight + if (apply_weights == 1) { + chi2 = chi2 - 2 * log(weights[il]); + } + // Return final updated chi2 value + return chi2; +} + +// The absolute magnitude prior +double prior::absmag_prior(onesource* source, double chi2, double reds, + int libtype, double dm, double mag0, double funz0) { + // predicted magnitudes within the babs filter renormalized by the scaling + double abs_mag; + abs_mag = mag0 - 2.5 * log10(dm); + // specific case when the redshift is 0 + if (reds < 1.e-10) abs_mag = abs_mag - funz0; + // Galaxy rejection + if ((abs_mag <= source->priorLib[0] || abs_mag >= source->priorLib[1]) && + libtype == 0) + chi2 = 1.e9; + // AGN rejection + if ((abs_mag <= source->priorLib[2] || abs_mag >= source->priorLib[3]) && + libtype == 1) + chi2 = 1.e9; + + return chi2; +} + +/* + Use prior info on N(z) as a function of magnitude (Iband) and of the type. + + For now we simply try to make this a method of the prior class which somehow + gets access to all the variables in the onesource class +*/ +double prior::nz_prior(onesource* source, double chi2, const double luv, + const double lnir, const double reds, + const array bp) { + double val = 1.; + // return val; + // + // } + int mod; + double zot = 0., kt = 0., alpt0 = 0.; + + // Define the NUV-R rest-frame color corrected for dust-extinction + double nuvk = -2.5 * (luv - lnir); + + // Apparent magnitude in i band, if possible + // Apparent magnitude in i band, if possible + double mi = -99.; + if (source->busnorma[bp[0]] == 1) { + mi = flux2mag(source->ab[bp[0]]); + } else if (source->busnorma[bp[1]] == 1) { + mi = flux2mag(source->ab[bp[1]]); + } else + return val; + + // Bright case I<20, no solution above 1, don't compute the prior + if (mi < 20.) { + if (reds > 1) val = 0; + return val; + } + // Bright case I<22, no solution above 2, but compute the prior + if (mi < 22. && reds > 2) { + val = 0; + return val; + } + + // Set up the parameters to define the redshift distribution + if (nuvk > 4.25) { + // Case E/S0 + // Color UV-K of PHOTO_230506/El_cww.sed.resample.new.resample15.inter + mod = 0; + zot = 0.45181; + kt = 0.13677; + alpt0 = 3.33078; + } else if (nuvk > 3.19 && nuvk < 4.25) { + // Case Sbc + // Color UV-K of PHOTO_230506/Sbc_cww.sed.resample.new.resample8.inter + mod = 1; + zot = 0.16560; + kt = 0.12983; + alpt0 = 1.42815; + } else if (nuvk > 1.9 && nuvk < 3.19) { + // Case Scd. Color UV-K of + // PHOTO_230506/Scd_cww.sed.resample.new.resample7.inter -19.4878 + 21.1501 + mod = 2; + zot = 0.21072; + kt = 0.14008; + alpt0 = 1.58310; + } else { + // Case Irr + mod = 3; + zot = 0.20418; + kt = 0.13773; + alpt0 = 1.34500; + } + + // P(z|T,m0) + double zmax = zot + kt * (mi - 20); + double pz = pow(reds, alpt0) * exp(-pow((reds / zmax), alpt0)); + + // P(T|m0) + double ktf[4] = {0.47165, 0.30663, 0.12715, -0.34437}; + double ft[4] = {0.43199, 0.07995, 0.31162, 0.21220}; + + // Ratio for each type + double rappSum = + ft[0] * exp(-ktf[0] * (mi - 20.0)) + ft[1] * exp(-ktf[1] * (mi - 20.0)); + rappSum += + ft[2] * exp(-ktf[2] * (mi - 20.0)) + ft[3] * exp(-ktf[3] * (mi - 20.0)); + double rapp = ft[mod] * exp(-ktf[mod] * (mi - 20.0)); + + // Normalisation of the probability function + // pcal=exp(gammln(1.d0/alpt0(mod)+1.d0)) + double pcal; + if (mod == 0) pcal = 0.89744; + if (mod == 1) pcal = 0.90868; + if (mod == 2) pcal = 0.89747; + if (mod == 3) pcal = 0.91760; + pcal = pow(zmax, alpt0 + 1) / alpt0 * pcal; + + // Final value + val = pz / pcal * rapp / rappSum; + + return chi2 - 2 * log(val); +} diff --git a/src/lib/prior.h b/src/lib/prior.h new file mode 100644 index 00000000..f289b6a9 --- /dev/null +++ b/src/lib/prior.h @@ -0,0 +1,88 @@ +/* + + 24/11/2023 + Class to deploy a prior which modifies chi2 as a function of z and model + properties. + +*/ + +// avoid multiple def of the same class +#ifndef PRIOR_H // check that this keyword has been set already +#define PRIOR_H // define the keyword to be checked + +// #include "onesource.h" +// #include "SED.h" +#include +#include + +using namespace std; + +// Use forward declaration and overwrite these empyy declarations +class SED; +class onesource; + +class prior +// The prior class should be initialised before the onesource fit method is +// called. We apply a number of generic prior forms typically based on redshift. +// We must initialise the prior before the fit to enable over writing in Python. +{ + private: + // Variables + int libtype, luv, lnir; + double mag0; + SED* lib; + + public: + // Variables + int apply_nz, apply_weights; + double chi2, reds; + vector weights; + // const onesource& source; + array bp; + // Constructor + prior() { + apply_nz = 1; // Apply n of z if it is specified unless we switch it off + apply_weights = 0; + luv = 0; + lnir = 0; + chi2 = 1.e9; + reds = 0.; + bp = {0, 0}; + } + // Methods + double update_chi2(onesource* source, double chi2, SED* lib, int il, + double dm, double funz0, const array& bp, + bool mabsGALprior, bool mabsAGNprior); + double absmag_prior(onesource* source, double chi2, double reds, int libtype, + double dm, double mag0, double funz0); + double nz_prior(onesource* source, double chi2, const double luv, + const double lnir, const double reds, const array bp); + + // Python overwriteable prior which takes the fullLib of SEDs and the source + // and returns a vector of weights equal in length to the fullLib. + // In general it is safer to simply overwrite the weights directly. + + // General function allowing pybind overwriting + std::function(vector, onesource*)> weights_function = + [](vector fullLib, onesource* source) { + // Default weight behavior all set to 1. + vector onesVector(fullLib.size(), 1.0f); + return onesVector; + }; + // Default behavior of the method + int set_weights(vector fullLib, onesource* source) { + // Call the current chi2_function + weights = weights_function(fullLib, source); + if (fullLib.size() != weights.size()) { + throw std::length_error("weights and fullLib are not the same length."); + } + return 0; + }; + // Allow setting a new weights_function from Python + void set_weights_function( + std::function(vector, onesource*)> func) { + weights_function = func; + }; +}; + +#endif \ No newline at end of file From 4e40ba41a22cb12787309b363b11fc5e382b7e5f Mon Sep 17 00:00:00 2001 From: Raphael Shirley Date: Wed, 5 Jun 2024 16:31:35 +0200 Subject: [PATCH 03/14] Added incldues and changed variable declaration This is building locally but throwing errors on GitHub --- src/lib/onesource.cpp | 1 + src/lib/prior.h | 8 ++++---- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/lib/onesource.cpp b/src/lib/onesource.cpp index 8a92ca8c..7c33bec2 100644 --- a/src/lib/onesource.cpp +++ b/src/lib/onesource.cpp @@ -25,6 +25,7 @@ #include "SED.h" #include "globals.h" #include "opa.h" +#include "prior.h" using namespace std; diff --git a/src/lib/prior.h b/src/lib/prior.h index f289b6a9..b5a76095 100644 --- a/src/lib/prior.h +++ b/src/lib/prior.h @@ -63,14 +63,14 @@ class prior // In general it is safer to simply overwrite the weights directly. // General function allowing pybind overwriting - std::function(vector, onesource*)> weights_function = - [](vector fullLib, onesource* source) { + std::function(vector, onesource*)> weights_function = + [](vector fullLib, onesource* source) { // Default weight behavior all set to 1. vector onesVector(fullLib.size(), 1.0f); return onesVector; }; // Default behavior of the method - int set_weights(vector fullLib, onesource* source) { + int set_weights(vector fullLib, onesource* source) { // Call the current chi2_function weights = weights_function(fullLib, source); if (fullLib.size() != weights.size()) { @@ -80,7 +80,7 @@ class prior }; // Allow setting a new weights_function from Python void set_weights_function( - std::function(vector, onesource*)> func) { + std::function(vector, onesource*)> func) { weights_function = func; }; }; From 41a9a2c1ea91f523a5efe075a61469ecea0b6435 Mon Sep 17 00:00:00 2001 From: Raphael Shirley Date: Wed, 5 Jun 2024 16:37:30 +0200 Subject: [PATCH 04/14] More incldues added to try to fix github build --- src/lib/_bindings.cc | 2 ++ src/lib/prior.cpp | 1 + src/lib/prior.h | 3 +-- 3 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/lib/_bindings.cc b/src/lib/_bindings.cc index f7b4ad81..c7726f27 100644 --- a/src/lib/_bindings.cc +++ b/src/lib/_bindings.cc @@ -1,3 +1,4 @@ +#include #include #include #include @@ -18,6 +19,7 @@ namespace py = pybind11; #include "oneElLambda.h" #include "opa.h" #include "photoz_lib.h" +#include "prior.h" template void applySEDLibTemplate(modT &m, std::string name) { diff --git a/src/lib/prior.cpp b/src/lib/prior.cpp index da59efa3..0fec352f 100644 --- a/src/lib/prior.cpp +++ b/src/lib/prior.cpp @@ -18,6 +18,7 @@ #include #include +#include #include "SED.h" #include "globals.h" diff --git a/src/lib/prior.h b/src/lib/prior.h index b5a76095..4ed0a53b 100644 --- a/src/lib/prior.h +++ b/src/lib/prior.h @@ -10,10 +10,9 @@ #ifndef PRIOR_H // check that this keyword has been set already #define PRIOR_H // define the keyword to be checked -// #include "onesource.h" -// #include "SED.h" #include #include +#include using namespace std; From df44ddc7b76ae4e1da3105a758a27e109ff0d09c Mon Sep 17 00:00:00 2001 From: Raphael Shirley Date: Wed, 5 Jun 2024 16:48:11 +0200 Subject: [PATCH 05/14] Add array to includes --- src/lib/prior.cpp | 1 + src/lib/prior.h | 1 + 2 files changed, 2 insertions(+) diff --git a/src/lib/prior.cpp b/src/lib/prior.cpp index 0fec352f..e850a20e 100644 --- a/src/lib/prior.cpp +++ b/src/lib/prior.cpp @@ -16,6 +16,7 @@ #include "prior.h" +#include #include #include #include diff --git a/src/lib/prior.h b/src/lib/prior.h index 4ed0a53b..15130c4b 100644 --- a/src/lib/prior.h +++ b/src/lib/prior.h @@ -10,6 +10,7 @@ #ifndef PRIOR_H // check that this keyword has been set already #define PRIOR_H // define the keyword to be checked +#include #include #include #include From 79f4ba83a0618cb8b6d7aab521adca040ded0aa3 Mon Sep 17 00:00:00 2001 From: Raphael Shirley Date: Wed, 5 Jun 2024 16:53:30 +0200 Subject: [PATCH 06/14] Add red to SED accessible variables --- src/lib/_bindings.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/src/lib/_bindings.cc b/src/lib/_bindings.cc index c7726f27..1c18c6eb 100644 --- a/src/lib/_bindings.cc +++ b/src/lib/_bindings.cc @@ -154,6 +154,7 @@ PYBIND11_MODULE(_lephare, mod) { .def_readonly("name", &SED::name) .def_readonly("nummod", &SED::nummod) .def_readonly("mag", &SED::mag) + .def_readonly("red", &SED::red) .def_readwrite("index_z0", &SED::index_z0) .def("is_gal", &SED::is_gal) .def("is_star", &SED::is_star) From 7888033ba3b4e9c14b7fd52c0b2364adc2fa1b60 Mon Sep 17 00:00:00 2001 From: Raphael Shirley Date: Thu, 6 Jun 2024 11:39:07 +0200 Subject: [PATCH 07/14] Cleaned up the notebook --- docs/notebooks/Prior.ipynb | 83 +++++++++++++------------------------- 1 file changed, 29 insertions(+), 54 deletions(-) diff --git a/docs/notebooks/Prior.ipynb b/docs/notebooks/Prior.ipynb index e43ba77f..abe90862 100644 --- a/docs/notebooks/Prior.ipynb +++ b/docs/notebooks/Prior.ipynb @@ -5,9 +5,9 @@ "id": "59aa16dc-8c12-4c0c-a75d-9251c7efd3b3", "metadata": {}, "source": [ - "# Prior\n", + "# Prior functionality\n", "\n", - "In this notebook I will demonstrate a basic redshift prior and show how more complex priors might be developed." + "In this notebook I will demonstrate a basic redshift prior and show how more complex priors might be developed. Essentially the user needs to define a set of weights for each model and pass these weights to the prior for a given source." ] }, { @@ -36,16 +36,6 @@ "keymap = lp.all_types_to_keymap(config)" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "c3a37641-afef-4211-89b0-f932e8feb68d", - "metadata": {}, - "outputs": [], - "source": [ - "config" - ] - }, { "cell_type": "code", "execution_count": null, @@ -53,7 +43,7 @@ "metadata": {}, "outputs": [], "source": [ - "# lp.prepare(keymap)" + "lp.prepare(keymap)" ] }, { @@ -115,7 +105,7 @@ " weights = np.array([1 for l in fullLib])\n", " # Get the redshift using the index\n", " # Seems odd. Is this really the best way to get the redshift?\n", - " red = np.array([photz.gridz[n - photz.fullLib[n].index_z0] for n in np.arange(len(photz.fullLib))])\n", + " red = np.array([photz.fullLib[n].red for n in np.arange(len(photz.fullLib))])\n", "\n", " weights[red > 2] = 1.0e-9\n", "\n", @@ -150,72 +140,54 @@ { "cell_type": "code", "execution_count": null, - "id": "6e4b8afa-c70e-46b4-9efd-40235ac45ad0", - "metadata": {}, - "outputs": [], - "source": [ - "mods = [l.nummod for l in photz.fullLib]\n", - "np.unique(mods)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "61ae1aea-32d2-4036-9bd8-755e2459a441", - "metadata": {}, - "outputs": [], - "source": [ - "photz.fullLib[3000].name, photz.fullLib[7001].name" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "2f072770-c20f-4be2-b84d-291441ad7813", + "id": "37c598f8-4d76-452d-ab59-48a7916cbdb9", "metadata": {}, "outputs": [], "source": [ - "star_ex.index_z0" + "# The redshift\n", + "gal_ex.red" ] }, { "cell_type": "code", "execution_count": null, - "id": "d846af5c-5de0-4b72-a66f-173c8bbd19f8", + "id": "5720827c-e3ad-4254-bccd-f3425cf6d391", "metadata": {}, "outputs": [], "source": [ - "gal_ex.index_z0" + "red = [photz.fullLib[n].red for n in np.arange(len(photz.fullLib))]" ] }, { "cell_type": "code", "execution_count": null, - "id": "5720827c-e3ad-4254-bccd-f3425cf6d391", + "id": "18256cb5-c9a1-4b14-9bb3-8dd87393c95f", "metadata": {}, "outputs": [], "source": [ - "red = [photz.gridz[n - photz.fullLib[n].index_z0] for n in np.arange(len(photz.fullLib))]" + "red[:5]" ] }, { "cell_type": "code", "execution_count": null, - "id": "18256cb5-c9a1-4b14-9bb3-8dd87393c95f", + "id": "ecdbfabd-b7b0-4acf-9374-88f9cd17b0be", "metadata": {}, "outputs": [], "source": [ - "red[:5]" + "np.unique(red)" ] }, { "cell_type": "code", "execution_count": null, - "id": "ecdbfabd-b7b0-4acf-9374-88f9cd17b0be", + "id": "6e4b8afa-c70e-46b4-9efd-40235ac45ad0", "metadata": {}, "outputs": [], "source": [ - "np.unique(red)" + "# The model numbers to set priors on a given template\n", + "nummod = [l.nummod for l in photz.fullLib]\n", + "np.unique(nummod)" ] }, { @@ -225,7 +197,8 @@ "metadata": {}, "outputs": [], "source": [ - "photz.fullLib[3000].ebv, photz.fullLib[4001].ebv" + "# The ebv are all zero. I think this is because they are applied to the photometry not the model\n", + "photz.fullLib[3000].ebv" ] }, { @@ -248,6 +221,16 @@ "# Set the prior function via pybind" ] }, + { + "cell_type": "markdown", + "id": "277e3103-71ee-4c5d-97b5-8b1cb0a60c6f", + "metadata": {}, + "source": [ + "## Compare results with and without the simple redshift prior\n", + "\n", + "Here we will run twice with the redshift prior above applied or not to show how it impacts the predicted redshifts." + ] + }, { "cell_type": "code", "execution_count": null, @@ -398,14 +381,6 @@ "plt.hist(t2[\"Z_BEST\"], bins=50)\n", "plt.xlabel(\"z without z<2 prior\")" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "84085877-f290-436c-aa3a-b886a64d9b9f", - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { From 08ab5e41d3e0813fa183e6c0e987b2f4512ffad8 Mon Sep 17 00:00:00 2001 From: Raphael Shirley Date: Fri, 7 Jun 2024 12:33:49 +0200 Subject: [PATCH 08/14] tabs to spaces in the Makefile --- src/lib/Makefile | 56 ++++++++++++++++++++++++------------------------ 1 file changed, 28 insertions(+), 28 deletions(-) diff --git a/src/lib/Makefile b/src/lib/Makefile index 93a3eb3a..f015f032 100755 --- a/src/lib/Makefile +++ b/src/lib/Makefile @@ -60,7 +60,7 @@ OBJECTS = \ keyword.o\ flt.o\ SED.o\ - SEDLib.o\ + SEDLib.o\ ext.o\ opa.o\ mag.o\ @@ -68,26 +68,26 @@ OBJECTS = \ cosmology.o\ onesource.o\ PDF.o\ - photoz_lib.o\ - prior.o\ + photoz_lib.o\ + prior.o\ all: $(FILES) - if [ ! -d "../../bin" ]; then mkdir ../../bin; fi - cp $(FILES) ../../bin/ + if [ ! -d "../../bin" ]; then mkdir ../../bin; fi + cp $(FILES) ../../bin/ %.o: %.cpp *h - $(FC) $(CFLAGS) -c $*.cpp -o $*.o + $(FC) $(CFLAGS) -c $*.cpp -o $*.o %: %.cpp $(LIB_LEPHARE) *h - $(FC) $(CFLAGS) $*.cpp -o $* $(LIB_LEPHARE) $(LIBS) + $(FC) $(CFLAGS) $*.cpp -o $* $(LIB_LEPHARE) $(LIBS) $(FILES): $(LIB_LEPHARE) $(LIB_LEPHARE): $(OBJECTS) - rm -f $(LIB_LEPHARE) - $(AR) rv $(LIB_LEPHARE) $(OBJECTS) - $(RANLIB) $(LIB_LEPHARE) + rm -f $(LIB_LEPHARE) + $(AR) rv $(LIB_LEPHARE) $(OBJECTS) + $(RANLIB) $(LIB_LEPHARE) # # ################################################# @@ -96,16 +96,16 @@ $(LIB_LEPHARE): $(OBJECTS) # # clean: - rm -f \#* - rm -f *.o - rm -f *.a - rm -f *~ + rm -f \#* + rm -f *.o + rm -f *.a + rm -f *~ - rm -f minuit/\#* - rm -f minuit/*.o - rm -f minuit/*.a - rm -f minuit/*~ - rm -f ../../bin/* + rm -f minuit/\#* + rm -f minuit/*.o + rm -f minuit/*.a + rm -f minuit/*~ + rm -f ../../bin/* # # @@ -115,21 +115,21 @@ clean: # # archi: - rm -f $(ARCHIVE) - tar cvf $(ARCHIVE) *.cpp *.h - gzip $(ARCHIVE) - echo "Archive: $(ARCHIVE).gz" + rm -f $(ARCHIVE) + tar cvf $(ARCHIVE) *.cpp *.h + gzip $(ARCHIVE) + echo "Archive: $(ARCHIVE).gz" # # ################################################ ####### Create working directory #### ################################################ work: - mkdir -p $(LEPHAREWORK) - mkdir -p $(LEPHAREWORK)/filt - mkdir -p $(LEPHAREWORK)/lib_bin - mkdir -p $(LEPHAREWORK)/lib_mag + mkdir -p $(LEPHAREWORK) + mkdir -p $(LEPHAREWORK)/filt + mkdir -p $(LEPHAREWORK)/lib_bin + mkdir -p $(LEPHAREWORK)/lib_mag doc: - doxygen ../doc/Doxyfile + doxygen ../doc/Doxyfile From cb36eceee931c77ae148f5b9ce7fa017525557c4 Mon Sep 17 00:00:00 2001 From: Raphael Shirley Date: Fri, 7 Jun 2024 14:28:56 +0200 Subject: [PATCH 09/14] revert to tabs in some cases in Makefile tabs are required in parts of the makefile depending what languag that part is in --- src/lib/Makefile | 50 ++++++++++++++++++++++++------------------------ 1 file changed, 25 insertions(+), 25 deletions(-) diff --git a/src/lib/Makefile b/src/lib/Makefile index f015f032..ef76bb08 100755 --- a/src/lib/Makefile +++ b/src/lib/Makefile @@ -72,22 +72,22 @@ OBJECTS = \ prior.o\ all: $(FILES) - if [ ! -d "../../bin" ]; then mkdir ../../bin; fi - cp $(FILES) ../../bin/ + if [ ! -d "../../bin" ]; then mkdir ../../bin; fi + cp $(FILES) ../../bin/ %.o: %.cpp *h - $(FC) $(CFLAGS) -c $*.cpp -o $*.o + $(FC) $(CFLAGS) -c $*.cpp -o $*.o %: %.cpp $(LIB_LEPHARE) *h - $(FC) $(CFLAGS) $*.cpp -o $* $(LIB_LEPHARE) $(LIBS) + $(FC) $(CFLAGS) $*.cpp -o $* $(LIB_LEPHARE) $(LIBS) $(FILES): $(LIB_LEPHARE) $(LIB_LEPHARE): $(OBJECTS) - rm -f $(LIB_LEPHARE) - $(AR) rv $(LIB_LEPHARE) $(OBJECTS) - $(RANLIB) $(LIB_LEPHARE) + rm -f $(LIB_LEPHARE) + $(AR) rv $(LIB_LEPHARE) $(OBJECTS) + $(RANLIB) $(LIB_LEPHARE) # # ################################################# @@ -96,16 +96,16 @@ $(LIB_LEPHARE): $(OBJECTS) # # clean: - rm -f \#* - rm -f *.o - rm -f *.a - rm -f *~ + rm -f \#* + rm -f *.o + rm -f *.a + rm -f *~ - rm -f minuit/\#* - rm -f minuit/*.o - rm -f minuit/*.a - rm -f minuit/*~ - rm -f ../../bin/* + rm -f minuit/\#* + rm -f minuit/*.o + rm -f minuit/*.a + rm -f minuit/*~ + rm -f ../../bin/* # # @@ -115,21 +115,21 @@ clean: # # archi: - rm -f $(ARCHIVE) - tar cvf $(ARCHIVE) *.cpp *.h - gzip $(ARCHIVE) - echo "Archive: $(ARCHIVE).gz" + rm -f $(ARCHIVE) + tar cvf $(ARCHIVE) *.cpp *.h + gzip $(ARCHIVE) + echo "Archive: $(ARCHIVE).gz" # # ################################################ ####### Create working directory #### ################################################ work: - mkdir -p $(LEPHAREWORK) - mkdir -p $(LEPHAREWORK)/filt - mkdir -p $(LEPHAREWORK)/lib_bin - mkdir -p $(LEPHAREWORK)/lib_mag + mkdir -p $(LEPHAREWORK) + mkdir -p $(LEPHAREWORK)/filt + mkdir -p $(LEPHAREWORK)/lib_bin + mkdir -p $(LEPHAREWORK)/lib_mag doc: - doxygen ../doc/Doxyfile + doxygen ../doc/Doxyfile From 61f936039c59a32c2ab5c6346e7aea66a0f8924c Mon Sep 17 00:00:00 2001 From: Raphael Shirley Date: Tue, 11 Jun 2024 16:24:05 +0200 Subject: [PATCH 10/14] implement idea for sending prior to process function It might be more efficient to calculate all weights before the loop if the fulllib is large --- src/lephare/process.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/lephare/process.py b/src/lephare/process.py index 0e2d1365..bb5e40a2 100644 --- a/src/lephare/process.py +++ b/src/lephare/process.py @@ -10,7 +10,7 @@ object_types = ["STAR", "GAL", "QSO"] -def process(config, input, col_names=None, standard_names=False, filename=None, offsets=None): +def process(config, input, col_names=None, standard_names=False, filename=None, offsets=None, prior=None): """Run all required steps to produce photometric redshift estimates Parameters @@ -28,6 +28,9 @@ def process(config, input, col_names=None, standard_names=False, filename=None, Output file name for the output catalogue. offsets : list If offsets are set autoadapt is not run but the set values are used. + prior : function + Function that converts a lephare.SED into a weight. We loop over all + SEDs to get the list of weights. Returns ======= @@ -77,7 +80,11 @@ def process(config, input, col_names=None, standard_names=False, filename=None, photozlist = [] for i in range(ng): one_obj = lp.onesource(i, photz.gridz) - one_obj.readsource(str(i), flux[i], flux_err[i], 63, zspec[i], " ") + one_obj.readsource(str(i), flux[i], flux_err[i], context[i], zspec[i], " ") + if prior is not None: + weights = [prior(lib, one_obj) for lib in photz.fullLib] + one_obj.priorObj.apply_weights = 1 + one_obj.priorObj.weights = weights photz.prep_data(one_obj) photozlist.append(one_obj) From 864841240c1bb7cf4978614d04f0a09d4a7cb24f Mon Sep 17 00:00:00 2001 From: Raphael Shirley Date: Tue, 11 Jun 2024 16:33:42 +0200 Subject: [PATCH 11/14] Added base test for prior functionality What else should we test and how do we get high coverage of the c side Need to test the nz prior at least --- tests/lephare/test_prior.py | 42 +++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) create mode 100644 tests/lephare/test_prior.py diff --git a/tests/lephare/test_prior.py b/tests/lephare/test_prior.py new file mode 100644 index 00000000..9fc4819c --- /dev/null +++ b/tests/lephare/test_prior.py @@ -0,0 +1,42 @@ +import os + +import lephare as lp +import numpy as np +from astropy.table import Table + + +def test_prior(test_data_dir: str): + test_dir = os.path.abspath(os.path.dirname(__file__)) + os.environ["LEPHAREDIR"] = os.path.join(test_dir, "../data") + os.environ["LEPHAREWORK"] = os.path.join(test_dir, "../tmp") + # Read the config file. + config_file = os.path.join(test_data_dir, "examples/COSMOS.para") + config = lp.read_config(config_file) + # Run preparation tasks. + lp.prepare(config) + # Read the test input catalogue + input_file = os.path.join(test_data_dir, "examples/COSMOS_first100specz.fits") + input = Table.read(input_file) + # Make a reduced column set for the minimal test + reduced_cols = [] + for c in input.colnames: + if not c.startswith("f"): + reduced_cols.append(c) + elif "IB527" in c: + reduced_cols.append(c) + elif "IB679" in c: + reduced_cols.append(c) + + # initialise the prior + def example_prior(lib, obj): + """An example prior which prohibits redshift above 2""" + weight = 1.0 # Be default do not change prior + if lib.red > 2: + weight = 1.0e-9 + return weight + + output, pdfs, zgrid = lp.process(config, input[reduced_cols], prior=example_prior) + # Check one of the outputs (results are terrible with just one filter and sparse z grid) + assert ~np.isclose(output["Z_BEST"][0], 3.58577857627795) # Should not equal the original redshift + assert len(zgrid) == 51 + assert np.sum(output["Z_BEST"] > 2) == 0 From 498d3127e5b753aa014bed081d5f93b92edcf89e Mon Sep 17 00:00:00 2001 From: Raphael Shirley Date: Tue, 11 Jun 2024 17:45:07 +0200 Subject: [PATCH 12/14] Basic test of absmag prior --- src/lib/_bindings.cc | 4 +++- tests/lephare/test_prior.py | 6 ++++++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/src/lib/_bindings.cc b/src/lib/_bindings.cc index 1c18c6eb..b0b6ad86 100644 --- a/src/lib/_bindings.cc +++ b/src/lib/_bindings.cc @@ -413,7 +413,7 @@ PYBIND11_MODULE(_lephare, mod) { .def_readonly("zsecScale", &onesource::zsecMod) .def_readonly("zsecAge", &onesource::zsecAge); - /******** CLASS prior *********/ + /******** CLASS prior, functions in prior.h *********/ py::class_( mod, "prior", "The prior class which allows Python overwriting of the weights") @@ -421,6 +421,8 @@ PYBIND11_MODULE(_lephare, mod) { "Constructor for prior class initialised with a onesource instance") .def("set_weights_function", &prior::set_weights_function, "Set a new weights function from Python") + .def("absmag_prior", &prior::absmag_prior) + .def("nz_prior", &prior::nz_prior) .def_readwrite("apply_nz", &prior::apply_nz, "int 1 if apply nzlim prior.") .def_readwrite("apply_weights", &prior::apply_weights, diff --git a/tests/lephare/test_prior.py b/tests/lephare/test_prior.py index 9fc4819c..7f715174 100644 --- a/tests/lephare/test_prior.py +++ b/tests/lephare/test_prior.py @@ -40,3 +40,9 @@ def example_prior(lib, obj): assert ~np.isclose(output["Z_BEST"][0], 3.58577857627795) # Should not equal the original redshift assert len(zgrid) == 51 assert np.sum(output["Z_BEST"] > 2) == 0 + + +def test_absmag_prior(): + p = lp.prior() + s = lp.onesource(0, [0, 1]) + assert np.isclose(p.absmag_prior(s, 0.0, 0.0, 0, 0.0, 0.0, 0.0), 1000000000.0) From b982f92362752825a335862ae880682992bcb763 Mon Sep 17 00:00:00 2001 From: Raphael Shirley Date: Tue, 11 Jun 2024 18:04:58 +0200 Subject: [PATCH 13/14] Basic tests for absmag and nzprior --- src/lib/_bindings.cc | 3 ++- tests/lephare/test_prior.py | 9 +++++++++ 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/src/lib/_bindings.cc b/src/lib/_bindings.cc index b0b6ad86..f34bc1e4 100644 --- a/src/lib/_bindings.cc +++ b/src/lib/_bindings.cc @@ -347,11 +347,13 @@ PYBIND11_MODULE(_lephare, mod) { .def("writeFullChi", &onesource::writeFullChi) // .def("write_pdz", &onesource::write_pdz) .def("best_spec_vec", &onesource::best_spec_vec) + .def_readwrite("ab", &onesource::ab) .def_readwrite("spec", &onesource::spec) .def_readwrite("consiz", &onesource::consiz) .def_readwrite("closest_red", &onesource::closest_red) // The prior for overwriting the weights .def_readwrite("priorObj", &onesource::priorObj) + .def_readwrite("busnorma", &onesource::busnorma) .def_readonly("pos", &onesource::pos) .def_readonly("cont", &onesource::cont) .def_readonly("pdfmap", &onesource::pdfmap) @@ -359,7 +361,6 @@ PYBIND11_MODULE(_lephare, mod) { .def_readonly("nbul", &onesource::nbul) .def_readonly("dm", &onesource::dm) .def_readonly("zs", &onesource::zs) - .def_readonly("ab", &onesource::ab) .def_readonly("sab", &onesource::sab) .def_readonly("mab", &onesource::mab) .def_readonly("msab", &onesource::msab) diff --git a/tests/lephare/test_prior.py b/tests/lephare/test_prior.py index 7f715174..aaf48fbf 100644 --- a/tests/lephare/test_prior.py +++ b/tests/lephare/test_prior.py @@ -46,3 +46,12 @@ def test_absmag_prior(): p = lp.prior() s = lp.onesource(0, [0, 1]) assert np.isclose(p.absmag_prior(s, 0.0, 0.0, 0, 0.0, 0.0, 0.0), 1000000000.0) + + +def test_nzprior(): + # What is the key behaviour we should be testing? + p = lp.prior() + s = lp.onesource(0, [0, 1]) + s.ab = [1.0e-9, 1.0e-9] + s.busnorma = [1, 1] + assert np.isclose(p.nz_prior(s, 0.0, 0.0, 0.0, 0.0, [0, 0]), 1) From 87bef7a92edbdbb86601ed8bd39f5bfede5a67e6 Mon Sep 17 00:00:00 2001 From: Raphael Shirley Date: Tue, 11 Jun 2024 18:53:10 +0200 Subject: [PATCH 14/14] Very basic test of prior updatechi2 func --- src/lib/_bindings.cc | 6 +++++- tests/lephare/test_prior.py | 10 ++++++++++ 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/src/lib/_bindings.cc b/src/lib/_bindings.cc index f34bc1e4..192e84e5 100644 --- a/src/lib/_bindings.cc +++ b/src/lib/_bindings.cc @@ -154,8 +154,11 @@ PYBIND11_MODULE(_lephare, mod) { .def_readonly("name", &SED::name) .def_readonly("nummod", &SED::nummod) .def_readonly("mag", &SED::mag) - .def_readonly("red", &SED::red) .def_readwrite("index_z0", &SED::index_z0) + .def_readwrite("luv", &SED::luv) + .def_readwrite("lnir", &SED::lnir) + .def_readwrite("mag0", &SED::mag0) + .def_readwrite("red", &SED::red) .def("is_gal", &SED::is_gal) .def("is_star", &SED::is_star) .def("is_qso", &SED::is_qso) @@ -424,6 +427,7 @@ PYBIND11_MODULE(_lephare, mod) { "Set a new weights function from Python") .def("absmag_prior", &prior::absmag_prior) .def("nz_prior", &prior::nz_prior) + .def("update_chi2", &prior::update_chi2) .def_readwrite("apply_nz", &prior::apply_nz, "int 1 if apply nzlim prior.") .def_readwrite("apply_weights", &prior::apply_weights, diff --git a/tests/lephare/test_prior.py b/tests/lephare/test_prior.py index aaf48fbf..80d50521 100644 --- a/tests/lephare/test_prior.py +++ b/tests/lephare/test_prior.py @@ -55,3 +55,13 @@ def test_nzprior(): s.ab = [1.0e-9, 1.0e-9] s.busnorma = [1, 1] assert np.isclose(p.nz_prior(s, 0.0, 0.0, 0.0, 0.0, [0, 0]), 1) + + +def test_update_chi2(): + p = lp.prior() + s = lp.onesource(0, [0, 1]) + star_sed = "o5v.sed.ext" + sed_filename = os.path.join(lp.LEPHAREDIR, "sed/STAR/", star_sed) + sed = lp.StarSED(star_sed, 1) + sed.read(sed_filename) + assert np.isclose(p.update_chi2(s, 0.0, sed, 0, 0.0, 0.0, [0, 0], False, False), 0.0)