diff --git a/skypy/halos/Literature_model_values (SHAM paper).txt b/skypy/halos/Literature_model_values (SHAM paper).txt new file mode 100644 index 000000000..9100b1edb --- /dev/null +++ b/skypy/halos/Literature_model_values (SHAM paper).txt @@ -0,0 +1,180 @@ +Arrays for literature models and the specific values used as featured in Davidson et al. (2025). Please see the paper for more information. + +Variables: +smhm: Ratio of galaxy and halo mass, unitless +halo_bins: Bins in halo mass used to average the Davidson et al. (2025) model, units of solar mass +halo_mids: Halo mass used to plot the smhm values, units of solar mass +error: Error for each smhm value, used in calculation of chi^2 + +### Model values (ratios and masses in log_10) + +# Moster et al. (2010) Central +smhm = [-2.57004016 -2.42705081 -2.28589751 -2.14773997 -2.01439977 -1.88862798 + -1.77432421 -1.67648241 -1.60053465 -1.55094212 -1.52953496 -1.53470726 + -1.56212721 -1.60637663 -1.66240233 -1.72622366 -1.79501688 -1.86690587 + -1.94069371 -2.01563547 -2.09127455 -2.16733355 -2.24364493 -2.32010779] +halo_bins = [[10.568965517241379, 10.706896551724139], [10.706896551724139, 10.844827586206897], [10.844827586206897, 10.982758620689655], [10.982758620689655, 11.120689655172413], [11.120689655172413, 11.258620689655173], [11.258620689655173, 11.39655172413793], [11.39655172413793, 11.53448275862069], [11.53448275862069, 11.672413793103448], [11.672413793103448, 11.810344827586206], [11.810344827586206, 11.948275862068966], [11.948275862068966, 12.086206896551724], [12.086206896551724, 12.224137931034482], [12.224137931034482, 12.362068965517242], [12.362068965517242, 12.5], [12.5, 12.637931034482758], [12.637931034482758, 12.775862068965518], [12.775862068965518, 12.913793103448276], [12.913793103448276, 13.051724137931034], [13.051724137931034, 13.189655172413794], [13.189655172413794, 13.327586206896552], [13.327586206896552, 13.46551724137931], [13.46551724137931, 13.60344827586207], [13.60344827586207, 13.741379310344827], [13.741379310344827, 13.879310344827585]] +halo_mids = [10.63793103 10.77586207 10.9137931 11.05172414 11.18965517 11.32758621 + 11.46551724 11.60344828 11.74137931 11.87931034 12.01724138 12.15517241 + 12.29310345 12.43103448 12.56896552 12.70689655 12.84482759 12.98275862 + 13.12068966 13.25862069 13.39655172 13.53448276 13.67241379 13.81034483] +error = [0.15672331 0.14696246 0.13516398 0.11777702 0.10907865 0.0911158 + 0.07569422 0.06522406 0.05195222 0.04356688 0.04001839 0.03944438 + 0.04274867 0.044372 0.0460492 0.04942211 0.05153041 0.05513094 + 0.05910394 0.06061898 0.06253746 0.06479315 0.06646424 0.06818571] + +# Moster et al. (2010) Satellite +smhm = [-2.75856607 -2.65145891 -2.54479023 -2.43882704 -2.33399508 -2.23096701 + -2.13079049 -2.03506001 -1.94611533 -1.86719836 -1.80241684 -1.75629629 + -1.73279702 -1.73405907 -1.75957841 -1.80640987] +halo_bins = [[10.21896551724138, 10.381034482758622], [10.381034482758622, 10.543103448275861], [10.543103448275861, 10.705172413793104], [10.705172413793104, 10.867241379310347], [10.867241379310347, 11.029310344827586], [11.029310344827586, 11.191379310344828], [11.191379310344828, 11.35344827586207], [11.35344827586207, 11.51551724137931], [11.51551724137931, 11.677586206896553], [11.677586206896553, 11.839655172413792], [11.839655172413792, 12.001724137931035], [12.001724137931035, 12.163793103448278], [12.163793103448278, 12.325862068965517], [12.325862068965517, 12.48793103448276], [12.48793103448276, 12.649999999999999], [12.649999999999999, 12.812068965517241]] +halo_mids = [10.3 10.46206897 10.62413793 10.7862069 10.94827586 11.11034483 + 11.27241379 11.43448276 11.59655172 11.75862069 11.92068966 12.08275862 + 12.24482759 12.40689655 12.56896552 12.73103448] +error = [0.17291564 0.17291564 0.15837369 0.14629925 0.13179159 0.11399535 + 0.09797352 0.07784644 0.06580277 0.05086428 0.04319014 0.0384812 + 0.0423966 0.04385776 0.0460492 0.04961714] + +# Moster et al. (2013) Central +smhm = [-2.28432357 -2.10346477 -1.92996775 -1.76921583 -1.62920783 -1.51981875 + -1.44973361 -1.42202404 -1.43223424 -1.4710998 -1.52895089 -1.59829823 + -1.67422813 -1.75380621 -1.83536956 -1.9180025 -2.0012086 -2.08472094 + -2.16839666 -2.25215945] +halo_bins = [[10.548275862068966, 10.713793103448277], [10.713793103448277, 10.879310344827587], [10.879310344827587, 11.044827586206896], [11.044827586206896, 11.210344827586209], [11.210344827586209, 11.375862068965517], [11.375862068965517, 11.541379310344828], [11.541379310344828, 11.706896551724139], [11.706896551724139, 11.872413793103448], [11.872413793103448, 12.03793103448276], [12.03793103448276, 12.203448275862069], [12.203448275862069, 12.36896551724138], [12.36896551724138, 12.53448275862069], [12.53448275862069, 12.7], [12.7, 12.865517241379312], [12.865517241379312, 13.03103448275862], [13.03103448275862, 13.196551724137931], [13.196551724137931, 13.362068965517242], [13.362068965517242, 13.527586206896551], [13.527586206896551, 13.693103448275863], [13.693103448275863, 13.858620689655172]] +halo_mids = [10.63103448 10.79655172 10.96206897 11.12758621 11.29310345 11.45862069 + 11.62413793 11.78965517 11.95517241 12.12068966 12.2862069 12.45172414 + 12.61724138 12.78275862 12.94827586 13.1137931 13.27931034 13.44482759 + 13.61034483 13.77586207] +error = [0.23522954 0.33757769 0.3419407 0.32097663 0.2984984 0.27033753 + 0.24000522 0.21196993 0.19068384 0.17968647 0.18283211 0.18886339 + 0.19534785 0.20429844 0.21042223 0.22251548 0.22256098 0.22256098 + 0.22256098 0.22256098] + +# Moster et al. (2013) Satellite +smhm = [-2.65647959 -2.4692828 -2.28432357 -2.10346477 -1.92996775 -1.76921583 + -1.62920783 -1.51981875 -1.44973361 -1.42202404 -1.43223424 -1.4710998 + -1.52895089 -1.59829823 -1.67422813 -1.75380621] +halo_bins = [[10.217241379310344, 10.382758620689657], [10.382758620689657, 10.548275862068966], [10.548275862068966, 10.713793103448277], [10.713793103448277, 10.879310344827587], [10.879310344827587, 11.044827586206896], [11.044827586206896, 11.210344827586209], [11.210344827586209, 11.375862068965517], [11.375862068965517, 11.541379310344828], [11.541379310344828, 11.706896551724139], [11.706896551724139, 11.872413793103448], [11.872413793103448, 12.03793103448276], [12.03793103448276, 12.203448275862069], [12.203448275862069, 12.36896551724138], [12.36896551724138, 12.53448275862069], [12.53448275862069, 12.7], [12.7, 12.865517241379312]] +halo_mids = [10.3 10.46551724 10.63103448 10.79655172 10.96206897 11.12758621 + 11.29310345 11.45862069 11.62413793 11.78965517 11.95517241 12.12068966 + 12.2862069 12.45172414 12.61724138 12.78275862] +error = [0.07317073 0.13284553 0.23522954 0.33757769 0.3419407 0.32097663 + 0.2984984 0.27033753 0.24000522 0.21196993 0.19068384 0.17968647 + 0.18283211 0.18886339 0.19534785 0.20429844] + +# Behroozi et al. (2010) Central +smhm = [-2.485207 -2.42245351 -2.35970378 -2.29696055 -2.23422859 -2.17151614 + -2.10883743 -2.04621709 -1.98369758 -1.921352 -1.85930582 -1.79777336 + -1.7371182 -1.67795071 -1.62127997 -1.56873684 -1.5228708 -1.48747894 + -1.4678339 -1.47055649 -1.50283566 -1.5709424 -1.67858543 -1.82616384 + -2.01162857 -2.23252182 -2.48798535 -2.77985791] +halo_bins = [[10.760549266211303, 10.809864741175637], [10.809864741175637, 10.859182094594342], [10.859182094594342, 10.908504579865326], [10.908504579865326, 10.957835952177302], [10.957835952177302, 11.007182710617883], [11.007182710617883, 11.056556096892267], [11.056556096892267, 11.105975540071805], [11.105975540071805, 11.155474577917708], [11.155474577917708, 11.205110995132], [11.205110995132, 11.254984078852235], [11.254984078852235, 11.305263726230866], [11.305263726230866, 11.356238885039033], [11.356238885039033, 11.408396525726285], [11.408396525726285, 11.46254637679387], [11.46254637679387, 11.520008406631039], [11.520008406631039, 11.582872787487615], [11.582872787487615, 11.654312800567126], [11.654312800567126, 11.738863312360706], [11.738863312360706, 11.84247105576249], [11.84247105576249, 11.972040900800486], [11.972040900800486, 12.134302819418874], [12.134302819418874, 12.334246669804703], [12.334246669804703, 12.573926359239696], [12.573926359239696, 12.852516895488126], [12.852516895488126, 13.167764846844236], [13.167764846844236, 13.518012204371459], [13.518012204371459, 13.90374921580872], [13.90374921580872, 14.32823756018442]] +halo_mids = [10.785207 10.83452248 10.88384171 10.93316745 10.98250446 11.03186097 + 11.08125123 11.13069985 11.1802493 11.22997269 11.27999547 11.33053198 + 11.38194579 11.43484726 11.49024549 11.54977132 11.61597425 11.69265135 + 11.78507528 11.89986684 12.04421497 12.22439067 12.44410267 12.70375005 + 13.00128374 13.33424595 13.70177845 14.10571998] +error = [0.24022346 0.24022346 0.24022346 0.24022346 0.24022346 0.24022346 + 0.24110716 0.24345663 0.24266 0.24186056 0.24105631 0.24024381 + 0.24082122 0.24145178 0.24211211 0.24282163 0.24117727 0.24117727 + 0.24117727 0.2406515 0.24253213 0.24301676 0.24246406 0.24024221 + 0.23959823 0.24022346 0.24167363 0.24070613] + +# Behroozi et al. (2010) Satellite +smhm = [-2.485207 -2.42245351 -2.35970378 -2.29696055 -2.23422859 -2.17151614 + -2.10883743 -2.04621709 -1.98369758 -1.921352 -1.85930582 -1.79777336 + -1.7371182 -1.67795071 -1.62127997 -1.56873684 -1.5228708 -1.48747894 + -1.4678339 -1.47055649 -1.50283566 -1.5709424 -1.67858543 -1.82616384] +halo_bins = [[10.760549266211303, 10.809864741175637], [10.809864741175637, 10.859182094594342], [10.859182094594342, 10.908504579865326], [10.908504579865326, 10.957835952177302], [10.957835952177302, 11.007182710617883], [11.007182710617883, 11.056556096892267], [11.056556096892267, 11.105975540071805], [11.105975540071805, 11.155474577917708], [11.155474577917708, 11.205110995132], [11.205110995132, 11.254984078852235], [11.254984078852235, 11.305263726230866], [11.305263726230866, 11.356238885039033], [11.356238885039033, 11.408396525726285], [11.408396525726285, 11.46254637679387], [11.46254637679387, 11.520008406631039], [11.520008406631039, 11.582872787487615], [11.582872787487615, 11.654312800567126], [11.654312800567126, 11.738863312360706], [11.738863312360706, 11.84247105576249], [11.84247105576249, 11.972040900800486], [11.972040900800486, 12.134302819418874], [12.134302819418874, 12.334246669804703], [12.334246669804703, 12.573926359239696], [12.573926359239696, 12.852516895488126]] +halo_mids = [10.785207 10.83452248 10.88384171 10.93316745 10.98250446 11.03186097 + 11.08125123 11.13069985 11.1802493 11.22997269 11.27999547 11.33053198 + 11.38194579 11.43484726 11.49024549 11.54977132 11.61597425 11.69265135 + 11.78507528 11.89986684 12.04421497 12.22439067 12.44410267 12.70375005] +error = [0.24022346 0.24022346 0.24022346 0.24022346 0.24022346 0.24022346 + 0.24110716 0.24345663 0.24266 0.24186056 0.24105631 0.24024381 + 0.24082122 0.24145178 0.24211211 0.24282163 0.24117727 0.24117727 + 0.24117727 0.2406515 0.24253213 0.24301676 0.24246406 0.24024221] + +# Behroozi et al. (2013) +smhm = [-2.48751602 -2.41660749 -2.30213401 -2.1246685 -1.91944642 -1.74299858 + -1.62731965 -1.57598542 -1.57845192 -1.62107266 -1.69200217 -1.78240931 + -1.88616788 -1.9991694 -2.11869479 -2.24294762 -2.37073493 -2.5012567 + -2.63396911] +halo_bins = [[10.60344827586207, 10.775862068965518], [10.775862068965518, 10.948275862068966], [10.948275862068966, 11.120689655172415], [11.120689655172415, 11.293103448275863], [11.293103448275863, 11.465517241379311], [11.465517241379311, 11.637931034482758], [11.637931034482758, 11.810344827586206], [11.810344827586206, 11.982758620689655], [11.982758620689655, 12.155172413793103], [12.155172413793103, 12.327586206896552], [12.327586206896552, 12.5], [12.5, 12.672413793103448], [12.672413793103448, 12.844827586206897], [12.844827586206897, 13.017241379310345], [13.017241379310345, 13.189655172413794], [13.189655172413794, 13.362068965517242], [13.362068965517242, 13.53448275862069], [13.53448275862069, 13.706896551724139], [13.706896551724139, 13.879310344827585]] +halo_mids = [10.68965517 10.86206897 11.03448276 11.20689655 11.37931034 11.55172414 + 11.72413793 11.89655172 12.06896552 12.24137931 12.4137931 12.5862069 + 12.75862069 12.93103448 13.10344828 13.27586207 13.44827586 13.62068966 + 13.79310345] +error = [0.14718669 0.14457671 0.14155105 0.13658943 0.13335122 0.13292164 + 0.13375841 0.13508384 0.1328937 0.13070356 0.13165321 0.13275358 + 0.13259402 0.13170842 0.13219501 0.13429368 0.13531976 0.13423426 + 0.13401214] + +# Rodríguez-Puebla et al. (2015) Red +smhm = [-2.99736031 -2.53354078 -2.13795798 -1.8509961 -1.6693305 -1.56631564 + -1.51615562 -1.50136918 -1.51113442 -1.53863006 -1.57933449 -1.63010548 + -1.68868395 -1.75341011 -1.8230479 -1.89666877 -1.97357062 -2.05321954 + -2.13520696 -2.21921798 -2.30500774] +halo_bins = [[11.032758620689656, 11.167241379310344], [11.167241379310344, 11.301724137931036], [11.301724137931036, 11.436206896551724], [11.436206896551724, 11.570689655172412], [11.570689655172412, 11.705172413793104], [11.705172413793104, 11.839655172413792], [11.839655172413792, 11.974137931034484], [11.974137931034484, 12.108620689655172], [12.108620689655172, 12.24310344827586], [12.24310344827586, 12.377586206896552], [12.377586206896552, 12.51206896551724], [12.51206896551724, 12.64655172413793], [12.64655172413793, 12.78103448275862], [12.78103448275862, 12.915517241379309], [12.915517241379309, 13.05], [13.05, 13.184482758620689], [13.184482758620689, 13.318965517241379], [13.318965517241379, 13.453448275862069], [13.453448275862069, 13.587931034482757], [13.587931034482757, 13.722413793103449], [13.722413793103449, 13.856896551724137]] +halo_mids = [11.1 11.23448276 11.36896552 11.50344828 11.63793103 11.77241379 + 11.90689655 12.04137931 12.17586207 12.31034483 12.44482759 12.57931034 + 12.7137931 12.84827586 12.98275862 13.11724138 13.25172414 13.3862069 + 13.52068966 13.65517241 13.78965517] +error = [0.12912088 0.13120149 0.14406974 0.14594754 0.14282864 0.13777258 + 0.13549995 0.13707432 0.1393947 0.13712814 0.13830438 0.13747584 + 0.13757254 0.1398297 0.14009933 0.13615089 0.13401917 0.13851688 + 0.13442656 0.13342085 0.13678363] + +# Rodríguez-Puebla et al. (2015) Blue +smhm = [-2.21891921 -2.06184743 -1.87505554 -1.68903858 -1.53029606 -1.41215075 + -1.33641199] +halo_bins = [[11.032758620689656, 11.167241379310344], [11.167241379310344, 11.301724137931036], [11.301724137931036, 11.436206896551724], [11.436206896551724, 11.570689655172412], [11.570689655172412, 11.705172413793104], [11.705172413793104, 11.839655172413792], [11.839655172413792, 11.974137931034484]] +halo_mids = [11.1 11.23448276 11.36896552 11.50344828 11.63793103 11.77241379 + 11.90689655] +error = [0.11538462 0.11993201 0.12175015 0.11608564 0.11823144 0.11533601 + 0.11682322 0.11987174 0.11906686 0.11969658 0.11614356 0.12088651] + +# Birrer et al. (2014) Central +smhm = [-3.09976892 -3.00600462 -2.78189533 -2.50945864 -2.12898722 -1.86049243 + -1.60692518 -1.45728534 -1.58330625 -1.7875487 -2.03217974 -2.22871097 + -2.46614663 -2.63304501] +halo_bins = [[10.387887131411066, 10.649586264491044], [10.649586264491044, 10.91146284869845], [10.91146284869845, 11.153380403421384], [11.153380403421384, 11.382635398620078], [11.382635398620078, 11.567985216224141], [11.567985216224141, 11.698340738842486], [11.698340738842486, 11.894441622224594], [11.894441622224594, 12.1505607391627], [12.1505607391627, 12.393550931121602], [12.393550931121602, 12.751677197238596], [12.751677197238596, 13.162859085298011], [13.162859085298011, 13.522827946210597], [13.522827946210597, 13.849897134636286], [13.849897134636286, 14.131299086763534]] +halo_mids = [10.51779347 10.78137906 11.04154664 11.26521417 11.50005663 11.63591381 + 11.76076767 12.02811557 12.2730059 12.51409596 12.98925844 13.33645973 + 13.70919616 13.99059811] +error = [0.07177935 0.09813791 0.1212464 0.12795643 0.1350017 0.13907741 + 0.14282303 0.15168693 0.16638035 0.17957712 0.16532225 0.13471862 + 0.10954019 0.09547009] + +# Birrer et al. (2014) Red +smhm = [-2.22203409 -1.86038535 -1.64416933 -1.49453798 -1.59081952 -1.7875487 + -2.03217974 -2.22871097 -2.46614663 -2.63304501] +halo_bins = [[11.44833528645327, 11.602656154691013], [11.602656154691013, 11.720258224777439], [11.720258224777439, 11.892543500280942], [11.892543500280942, 12.135884337863782], [12.135884337863782, 12.380738746820189], [12.380738746820189, 12.751677197238596], [12.751677197238596, 13.162859085298011], [13.162859085298011, 13.522827946210597], [13.522827946210597, 13.849897134636286], [13.849897134636286, 14.131299086763534]] +halo_mids = [11.52549572 11.67981659 11.76069986 12.02438714 12.24738154 12.51409596 + 12.98925844 13.33645973 13.70919616 13.99059811] +error = [0.18974504 0.18820183 0.187393 0.18451226 0.18005237 0.17471808 + 0.16521483 0.13471862 0.10954019 0.09547009] + +# Birrer et al. (2014) Blue +smhm = [-2.50945021 -2.08050613 -1.87166536 -1.60319063 -1.31950041] +halo_bins = [[11.138725793841601, 11.397313679989932], [11.397313679989932, 11.456034937801823], [11.511104470634775, 11.70016345738368], [11.70016345738368, 11.892740921058529], [11.892740921058529, 12.14935637415585]] +halo_mids = [11.26887271 11.52575465 11.63589372 11.76443319 12.02104865] +error = [0.12537745 0.12845472 0.12184638 0.11413401 0.09894757] + +# Rodríguez-Puebla et al. (2012) +smhm = [-2.36182991 -2.29421066 -2.22521308 -2.17233385 -2.0448248 -1.96986527 + -1.90224601 -1.85913401 -1.75049998 -1.69097017 -1.63554509 -1.59651799 + -1.50561126 -1.46832248 -1.43266169 -1.4105227 -1.37934481 -1.37200081 + -1.37999256 -1.39043151 -1.42462393 -1.46085755 -1.50204824 -1.53737123 + -1.62471421 -1.6800135 ] +halo_bins = [[10.809389202970845, 10.84558252965303], [10.84558252965303, 10.87951377341758], [10.87951377341758, 10.90665876842922], [10.90665876842922, 10.949638343864315], [10.949638343864315, 10.999404168052322], [10.999404168052322, 11.035597494734507], [11.035597494734507, 11.06726665558142], [11.06726665558142, 11.112508313934153], [11.112508313934153, 11.164536221039796], [11.164536221039796, 11.20751579647489], [11.20751579647489, 11.248233288992349], [11.248233288992349, 11.309309527768539], [11.309309527768539, 11.379434098215274], [11.379434098215274, 11.43598617115619], [11.43598617115619, 11.49027616117947], [11.49027616117947, 11.576235312049661], [11.576235312049661, 11.680291126260945], [11.680291126260945, 11.77077444296641], [11.77077444296641, 11.84542317924842], [11.84542317924842, 11.944954827624432], [11.944954827624432, 12.067107305176808], [12.067107305176808, 12.16663895355282], [12.16663895355282, 12.252598104423011], [12.252598104423011, 12.37475058197539], [12.37475058197539, 12.510475557033587], [12.510475557033587, 12.6100072054096]] +halo_mids = [10.82748587 10.86367919 10.89534835 10.91796918 10.9813075 11.01750083 + 11.05369416 11.08083915 11.14417747 11.18489497 11.23013663 11.26632995 + 11.3522891 11.40657909 11.46539325 11.51515907 11.63731155 11.7232707 + 11.81827818 11.87256817 12.01734148 12.11687313 12.21640478 12.28879143 + 12.46070973 12.56024138] +error = [0.24804178 0.24804178 0.24804178 0.24804178 0.24804178 0.24804178 + 0.24804178 0.24804178 0.24804178 0.24804178 0.24804178 0.24804178 + 0.24804178 0.24804178 0.24804178 0.24804178 0.24804178 0.24804178 + 0.24804178 0.24804178 0.24804178 0.24804178 0.24804178 0.24804178 + 0.24804178 0.24804178] diff --git a/skypy/halos/__init__.py b/skypy/halos/__init__.py index e5ded3234..3f40d254b 100644 --- a/skypy/halos/__init__.py +++ b/skypy/halos/__init__.py @@ -14,9 +14,18 @@ __all__ = [ 'colossus_mf', + 'quenching_funct', + 'find_min', + 'run_file', + 'gen_sub_cat', + 'galaxy_cat', + 'assignment', + 'sham_arrays', + 'run_sham', ] from . import abundance_matching # noqa F401,F403 from . import mass # noqa F401,F403 from . import quenching # noqa F401,F403 from ._colossus import colossus_mf # noqa F401,F403 +from . import sham # noqa F401,F403 diff --git a/skypy/halos/bc_file.yaml b/skypy/halos/bc_file.yaml new file mode 100644 index 000000000..8301c72e6 --- /dev/null +++ b/skypy/halos/bc_file.yaml @@ -0,0 +1,20 @@ +m_star: !numpy.power [10, 10.65] +phi_star: !numpy.power [10, -2.98] +alpha_val: -1.48 +m_min: !numpy.power [10, 7.0] +m_max: !numpy.power [10, 14.0] +sky_area: 100.0 deg2 +z_range: !numpy.linspace [0.010004503125093451, 0.09999991945507442, 100] +cosmology: !astropy.cosmology.FlatLambdaCDM + H0: 70.0 + Om0: 0.3 +tables: + galaxy: + z, sm: !skypy.galaxies.schechter_smf + redshift: $z_range + m_star: $m_star + phi_star: $phi_star + alpha: $alpha_val + m_min: $m_min + m_max: $m_max + sky_area: $sky_area diff --git a/skypy/halos/bs_file.yaml b/skypy/halos/bs_file.yaml new file mode 100644 index 000000000..ca04db43b --- /dev/null +++ b/skypy/halos/bs_file.yaml @@ -0,0 +1,20 @@ +m_star: !numpy.power [10, 10.55] +phi_star: !numpy.power [10, -3.96] +alpha_val: -1.53 +m_min: !numpy.power [10, 7.0] +m_max: !numpy.power [10, 14.0] +sky_area: 100.0 deg2 +z_range: !numpy.linspace [0.010004503125093451, 0.09999991945507442, 100] +cosmology: !astropy.cosmology.FlatLambdaCDM + H0: 70.0 + Om0: 0.3 +tables: + galaxy: + z, sm: !skypy.galaxies.schechter_smf + redshift: $z_range + m_star: $m_star + phi_star: $phi_star + alpha: $alpha_val + m_min: $m_min + m_max: $m_max + sky_area: $sky_area diff --git a/skypy/halos/halo.yaml b/skypy/halos/halo.yaml new file mode 100644 index 000000000..8145b28e2 --- /dev/null +++ b/skypy/halos/halo.yaml @@ -0,0 +1,26 @@ +parameters: + model: 'sheth99' + mdef: 'fof' + m_min: 1.e+9 + m_max: 1.e+15 + sky_area: 600. deg2 # Same as galaxies + sigma_8: 0.8 + ns: 0.96 +cosmology: !astropy.cosmology.FlatLambdaCDM + H0: 70 + Om0: 0.3 + name: 'FlatLambdaCDM' # Requires user defined name + Ob0: 0.045 + Tcmb0: 2.7 +z_range: !numpy.linspace [0.01, 0.1, 100] # Same as galaxies +tables: + halo: + z, mass: !skypy.halos.colossus_mf + redshift: $z_range + model: $model + mdef: $mdef + m_min: $m_min + m_max: $m_max + sky_area: $sky_area + sigma8: $sigma_8 + ns: $ns diff --git a/skypy/halos/mass.py b/skypy/halos/mass.py index 919b97e1f..1d783fe6b 100644 --- a/skypy/halos/mass.py +++ b/skypy/halos/mass.py @@ -348,7 +348,7 @@ def number_subhalos(halo_mass, alpha, beta, gamma_M, x, m_min, noise=True): subhalos. m_min : array_like Original mass of the least massive subhalo, in units of solar mass. - Current stripped mass is given by :math:`x m_{min}`. + Current stripped mass is given by :math:`m_{min}/x`. noise : bool, optional Poisson-sample the number of subhalos. Default is `True`. @@ -383,6 +383,7 @@ def number_subhalos(halo_mass, alpha, beta, gamma_M, x, m_min, noise=True): # Random number of subhalos following a Poisson distribution # with mean n_subhalos + # TODO Poisson rounds the value, otherwise doesn't. mass_sampler needs rounded numbers return np.random.poisson(n_subhalos) if noise else n_subhalos @@ -440,9 +441,12 @@ def subhalo_mass_sampler(halo_mass, nsubhalos, alpha, beta, halo_mass = np.atleast_1d(halo_mass) nsubhalos = np.atleast_1d(nsubhalos) subhalo_list = [] + # TODO fix for not using Poisson in above function + # nsubhalos = [int(x) for x in nsubhalos] + for M, n in zip(halo_mass, nsubhalos): x_min = m_min / (x * beta * M) x_max = 0.5 / (x * beta) - subhalo_mass = schechter(alpha, x_min, x_max, resolution, size=n, scale=x * beta * M) + subhalo_mass = schechter(-alpha, x_min, x_max, resolution, size=n, scale=x * beta * M) subhalo_list.append(subhalo_mass) return np.concatenate(subhalo_list) diff --git a/skypy/halos/rc_file.yaml b/skypy/halos/rc_file.yaml new file mode 100644 index 000000000..c59cfb174 --- /dev/null +++ b/skypy/halos/rc_file.yaml @@ -0,0 +1,20 @@ +m_star: !numpy.power [10, 10.58] +phi_star: !numpy.power [10, -2.77] +alpha_val: -0.33 +m_min: !numpy.power [10, 7.0] +m_max: !numpy.power [10, 14.0] +sky_area: 100.0 deg2 +z_range: !numpy.linspace [0.010004503125093451, 0.09999991945507442, 100] +cosmology: !astropy.cosmology.FlatLambdaCDM + H0: 70.0 + Om0: 0.3 +tables: + galaxy: + z, sm: !skypy.galaxies.schechter_smf + redshift: $z_range + m_star: $m_star + phi_star: $phi_star + alpha: $alpha_val + m_min: $m_min + m_max: $m_max + sky_area: $sky_area diff --git a/skypy/halos/red_central.yaml b/skypy/halos/red_central.yaml new file mode 100644 index 000000000..6d41a452f --- /dev/null +++ b/skypy/halos/red_central.yaml @@ -0,0 +1,20 @@ +m_star: !numpy.power [10, 10.75] +phi_star: !numpy.power [10, -2.37] +alpha_val: -0.18 +m_min: !numpy.power [10, 7.0] +m_max: !numpy.power [10, 14.0] +sky_area: 600.0 deg2 +z_range: !numpy.linspace [0.01, 0.1, 100] +cosmology: !astropy.cosmology.FlatLambdaCDM + H0: 70.0 + Om0: 0.3 +tables: + galaxy: + z, sm: !skypy.galaxies.schechter_smf + redshift: $z_range + m_star: $m_star + phi_star: $phi_star + alpha: $alpha_val + m_min: $m_min + m_max: $m_max + sky_area: $sky_area diff --git a/skypy/halos/rs_file.yaml b/skypy/halos/rs_file.yaml new file mode 100644 index 000000000..e008b1c7d --- /dev/null +++ b/skypy/halos/rs_file.yaml @@ -0,0 +1,20 @@ +m_star: !numpy.power [10, 10.64] +phi_star: !numpy.power [10, -4.24] +alpha_val: -1.54 +m_min: !numpy.power [10, 7.0] +m_max: !numpy.power [10, 14.0] +sky_area: 100.0 deg2 +z_range: !numpy.linspace [0.010004503125093451, 0.09999991945507442, 100] +cosmology: !astropy.cosmology.FlatLambdaCDM + H0: 70.0 + Om0: 0.3 +tables: + galaxy: + z, sm: !skypy.galaxies.schechter_smf + redshift: $z_range + m_star: $m_star + phi_star: $phi_star + alpha: $alpha_val + m_min: $m_min + m_max: $m_max + sky_area: $sky_area diff --git a/skypy/halos/sham.py b/skypy/halos/sham.py new file mode 100644 index 000000000..b3af5a0dd --- /dev/null +++ b/skypy/halos/sham.py @@ -0,0 +1,1107 @@ +'''Subhalo Abundance Matching module + +This module generates halos and galaxies using SkyPy from user defined +parameters, decides which halos will be assigned which galaxies and then +outputs the matched arrays + +Models +====== +.. autosummary:: + :nosignatures: + :toctree: ../api/ + + quenching_funct + find_min + run_file + gen_sub_cat + assignment + scatter_proxy + sham_plots + run_sham +''' + +# Imports +import numpy as np +from skypy.pipeline import Pipeline +from skypy.halos import mass # Vale and Ostriker 2004 num of subhalos +from skypy.galaxies import schechter_smf +from time import time +from scipy.special import erf # Error function for quenching +from scipy.integrate import trapezoid as trap # Integration of galaxies +from astropy import units as u + +try: + import colossus # noqa F401 +except ImportError: + HAS_COLOSSUS = False +else: + HAS_COLOSSUS = True + +__all__ = [ + 'quenching_funct', + 'find_min', + 'run_file', + 'gen_sub_cat', + 'assignment', + 'scatter_proxy', + 'sham_plots', + 'run_sham', + ] + +# General functions + + +# Quenching function +def quenching_funct(mass, M_mu, sigma, baseline=0): + r'''Quenching function applied to halos + This function computes the fraction of quenched halos (halos assigned + quenched galaxies) as a function of mass using the parameterisation found + in [1]_ and returns a truth list. + + Parameters + ----------- + mass : (nm,) array_like + Array for the halo mass, in units of solar mass. + M_mu : float + The mean of the error function, or the mass at which half of the halos + are quenched, in units of solar mass. + sigma : float + The standard deviation of the error function, controls the width of + the transition from blue to red galaxies + baseline : float + The initial value of the error function, or the fraction quenched at + their lowest mass. For central halos this should be zero but for + subhalos this will be non zero. + + Returns + -------- + truth_array: (nm,) array_like + Array which sets which halos are quenched + + Examples + --------- + >>> import numpy as np + >>> from skypy.halos import sham + + This example will generate a central halo catalogue from a yaml file and + then pass it through the quenching function to find which halos are + quenched, following the mass distribution of the quenching function: + + >>> #Parameters + >>> M_mu, sigma, base = 10**(12), 0.4, 0.3 + >>> h_mass = sham.run_file('halo.yaml', table1 = 'halo', info1 = 'mass') + >>> h_quench = sham.quenching_funct(h_mass, M_mu, sigma, base) + + + References + ---------- + .. [1] Peng Y.-j., et al., 2010, Astrophysical Journal, 721, 193 + ''' + mass = np.atleast_1d(mass) + + # Errors + if M_mu <= 0: + raise Exception('M_mu cannot be negative or 0 solar masses') + if sigma <= 0: + raise Exception('sigma cannot be negative or 0') + if baseline < 0: + baseline = -baseline + raise Warning('baseline cannot be negative, set_baseline = -baseline') + if baseline > 1: + raise Exception('baseline cannot be more than 1') + if len(mass) == 1: + raise TypeError('Mass must be array-like') + + # Only works with an increasing array due to normalisation? + qu_id = np.arange(0, len(mass)) # Random order elements (ie halo original order) + + # Elements of sorted halos + stack1 = np.stack((mass, qu_id), axis=1) + order1 = stack1[np.argsort(stack1[:, 0])] # Order on the halos + mass = order1[:, 0] + order_qu_id = order1[:, 1] + + calc = np.log10(mass/M_mu)/sigma + old_prob = 0.5 * (1.0 + erf(calc/np.sqrt(2))) # Base error funct, 0 -> 1 + add_val = baseline/(1-baseline) + prob_sub = old_prob + add_val # Add value to plot to get baseline + prob_array = prob_sub/prob_sub[-1] # Normalise, base -> 1 + rand_array = np.random.uniform(size=len(mass)) + qu_arr = rand_array < prob_array + + # Reorder truth to match original order + stack2 = np.stack((order_qu_id, qu_arr), axis=1) + order2 = stack2[np.argsort(stack2[:, 0])] + return order2[:, 1] + + +# Find minimum galaxy mass for given halo population through the integral +def find_min(m_star, phi_star, alpha, cosmology, z_range, skyarea, max_mass, no_halo, + print_out=False, run_anyway=False): + r'''Function to find a minimum mass for a galaxy type for a given number + of halos + + This function computes a look up table of the integral of the galaxy + Schechter function and interpolates on the number of halos to find the + minimum mass. + + Parameters + ----------- + m_star : float + Exponential tail off of Schechter function, units of solar mass + phi_star : float + Normalisation for the Schechter function, units of Mpc^{-3} + alpha : float + Power law of Schechter function + cosmology : Cosmology + Astropy cosmology object for calculating comoving volume + z_range : (2, ), array_like + Minimum and maximum redshift of galaxies + skyarea : float + Sky area galaxies are 'observed' from, units of deg^2 + max_mass : float + Maximum mass to integrate Schechter function + no_halo : float + Number of halos to match to given galaxy function + run_anyway (debug): Boolean + True/False statement to determine if function runs outside of integral + look up table + + Returns + -------- + min_mass_log: float + Log10 of the minimum mass returned from interpolation + + Examples + --------- + >>> from astropy.cosmology import FlatLambdaCDM + >>> from skypy.halos import sham + + This example uses the Schechter parameters for red centrals from [1]_ and + a number of halos these might be paired with to generate a minimum mass: + + >>> #Define input variables + >>> m_star, phi_star, alpha = 10**(10.75), 10**(-2.37), -0.18 + >>> cosmology = FlatLambdaCDM(H0=70, Om0=0.3, name='FlatLambdaCDM') + >>> z_range = [0.01, 0.1] + >>> skyarea = 600 + >>> max_mass = 10**(14) + >>> no_halo = 10124 + >>> #Run function + >>> min_mass = sham.find_min(m_star, phi_star, alpha, cosmology, + ... z_range, skyarea, max_mass, no_halo) + + + References + ---------- + .. [1] Weigel A. K., Schawinski K., Bruderer C., 2016, Monthly Notices of the + Royal Astronomical Society, 459, 2150 + ''' + z_range = np.atleast_1d(z_range) + + # Errors + if z_range.shape != (2, ): + raise Exception('The wrong number of redshifts were given') + if z_range[0] < 0 or z_range[1] < 0: + raise Exception('Redshift cannot be negative') + if z_range[1] <= z_range[0]: + raise Exception('The second redshift should be more than the first') + if alpha >= 0: + alpha = -1*alpha + raise Warning('Schechter function defined so alpha < 0, set_alpha = -alpha') + if m_star <= 0 or phi_star <= 0: + raise Exception('M* and phi* must be positive and non-zero numbers') + if skyarea <= 0: + raise Exception('The skyarea must be a positive non-zero number') + if max_mass <= 10**6: + raise Exception('The maximum mass must be more than 10^6 solar masses') + + skyarea = skyarea*u.deg**2 + z = np.linspace(z_range[0], z_range[1], 1000) + dVdz = (cosmology.differential_comoving_volume(z)*skyarea).to_value('Mpc3') + + # Other + mass_mins = 10**np.arange(6, 10, 0.1) # Minimums to integrate, linear in logspace + phi_m = phi_star/m_star + # Want all integral bins to have same resolution in log space, should give accurate integral + res = 10**(-2.3) + + # Integrate + int_vals = [] + for ii in mass_mins: + mass_bin = 10**np.arange(np.log10(ii), np.log10(max_mass), res) # Masses to integrate over + m_m = mass_bin/m_star + dndmdV = phi_m*np.e**(-m_m)*(m_m)**alpha # Mass function + dndV = trap(dndmdV, mass_bin) + dn = dndV*dVdz + int_vals.append(trap(dn, z)) # Integral per mass bin + + int_vals = np.array(int_vals) + # Integral must be increasing so flip arrays + min_mass = np.interp(no_halo, np.flip(int_vals), np.flip(mass_mins)) + + # Region within Poisson error + poisson_min = min(int_vals) - np.sqrt(min(int_vals)) + poisson_max = max(int_vals) + np.sqrt(max(int_vals)) + + if (no_halo < poisson_min or no_halo > poisson_max) and not run_anyway: + # Interpolation falls outside of created range + raise Exception('Number of halos not within reachable bounds') + + elif (no_halo < poisson_min or no_halo > poisson_max) and run_anyway: + if print_out: + print('Outside of interpolation, but running anyway') + return 10**(7) + + return min_mass + + +# Generate catalogues +def run_file(file_name, table1, info1, info2=None): + r'''Function that runs a yaml file to generate a catalogue + + This function runs a yaml file using the SkyPy pipeline and + produces a catalogue from the file. + + Parameters + ----------- + file_name : str + String of file name to be run + table1 : str + String of table to access from file + info1 : str + String of variable to access mass, this must be provided + info2 : str + String of variable to access redshift + + Returns + -------- + catalogue: (nm, ) array_like + List of masses + redshifts: (nm, ) array_like + List of redshifts generated for each mass, if requested + + Examples + --------- + >>> from skypy.halos import sham + + This example uses a yaml file for a galaxy type and returns a + galaxy catalogue: + + >>> galaxy_cat = sham.run_file('red_central.yaml', 'galaxy', 'sm') + + This example uses a yaml file for halos and returns a + halo catalogue and their generated redshifts: + + >>> halo_cat, h_z = sham.run_file('halo.yaml', 'halo', 'mass', 'z') + + References + ---------- + .. + ''' + import os + # Errors + if type(file_name) is not str: + raise Exception('File name must be a string') + if not os.path.exists(file_name): + raise Exception('File does not exist') + + # Use pipeline + pipe = Pipeline.read(file_name) + pipe.execute() + + # Get information + info = pipe[table1] + cat = info[info1] + if info2 is not None: + z = info[info2] + return cat, z + + return cat + + +# Generate subhalo catalogues +def gen_sub_cat(parent_halo, z_halo, sub_alpha, sub_beta, sub_gamma, sub_x): + r'''Function that generates a catalogue of subhalos from a population + of parent halos + + This function uses the conditional mass function and occupation number + from [1]_ to generate a catalogue of subhalos for each parent halo + in the halo catalogue + + Parameters + ----------- + parent_halo : (nm, ), array_like + Parent halo catalogue used to generate subhalos, units + of solar mass + z_halo : (nm, ), array_like + Redshifts of parent halos + sub_alpha : float + Low mass power law slope of conditional mass function + sub_beta : float + Exponential cutoff of subhalo masses and a fraction of the + parent halo + sub_gamma : float + Present day mass fraction of parent halo in sum of generated + subhalo mass + sub_x : float + Factor that allows for stripping of mass for current subhalos, + x = 1 is the current day stripped mass, x > 1 for unstripped + subhalo mass + + Returns + -------- + ID_halo: (nm, ), array_like + ID values for parent halos + sub_masses: (ns, ), array_like + Subhalo masses for parent halos + ID_sub: (ns, ), array_like + ID values for subhalos to assign them to their parent halos + z_sub: (ns, ), array_like + Redshifts of subhalos (same as their parent halo) + + Examples + --------- + >>> from skypy.halos import sham + + This example generates a halo catalogue and then creates a + subhalo catalogue: + + >>> halo_cat, h_z = sham.run_file('halo.yaml', 'halo', 'mass', 'z') + >>> #Variables + >>> a, b, c, x = 1.91, 0.39, 0.1, 3 + >>> ID_halo, sub_masses, ID_sub, z_sub = sham.gen_sub_cat(halo_cat, h_z, + ... a, b, c, x) + + + References + ---------- + .. [1] Vale A., Ostriker J. P., 2004, Monthly Notices of the Royal + Astronomical Society, 353, 189 + ''' + parent_halo = np.atleast_1d(parent_halo) + z_halo = np.atleast_1d(z_halo) + + # Errors + if len(parent_halo) != len(z_halo): + raise Exception('Catalogue of halos and redshifts must be the same length') + if len(np.where(parent_halo <= 0)[0]) > 0: + raise Exception('Masses in catalogue should be positive and non-zero') + if len(np.where(z_halo < 0)[0]) > 0: + raise Exception('Redshifts in catalogue should be positive') + if sub_alpha < 0: + sub_alpha = -1*sub_alpha + raise Warning('Subhalo mass function defined alpha > 0, set_alpha = -alpha') + if sub_alpha >= 2: + raise Exception('Subhalo alpha must be less than 2') + if sub_x < 1: + raise Exception('Subhalo x cannot be less than 1') + if sub_beta <= 0 or sub_beta > 1: + raise Exception('Subhalo beta must be between 0 and 1') + if sub_gamma < 0 or sub_gamma > 1: + raise Exception('Subhalo gamma must be between 0 and 1') + + if parent_halo.size == 1: # ie only one halo is provided + m_min = 10**(10) + else: + m_min = min(parent_halo) # Min mass of subhalo to generate (RESOLUTION OF PARENT HALOS) + ID_halo = -1*np.arange(1, len(parent_halo)+1, dtype=int) # Halo IDs + + # Get list of halos that will have subhalos + halo_to_sub = parent_halo[np.where(parent_halo >= 10**(10))] # TODO change to a larger number? + ID_to_sub = -1*ID_halo[np.where(parent_halo >= 10**(10))] + z_to_sub = z_halo[np.where(parent_halo >= 10**(10))] + + # Get subhalos + no_sub = mass.number_subhalos(halo_to_sub, sub_alpha, sub_beta, + sub_gamma, sub_x, m_min, noise=True) + sub_masses = mass.subhalo_mass_sampler(halo_to_sub, no_sub, sub_alpha, sub_beta, sub_x, m_min) + + # Assign subhalos to parents by ID value + ID_sub = [] + z_sub = [] + ID_count = 0 + for ii in no_sub: + ID_set = ID_to_sub[ID_count] + ID_sub.extend(ID_set*np.ones(ii, dtype=int)) # Positive ID for subhalo + z_sub.extend(z_to_sub[ID_count]*np.ones(ii)) # Same redshift as the parent + ID_count += 1 + + # Delete any generated subhalos smaller than the resolution + del_sub = np.where(sub_masses < m_min)[0] + sub_masses = np.delete(sub_masses, del_sub) + ID_sub = np.delete(ID_sub, del_sub) + z_sub = np.delete(z_sub, del_sub) + + return ID_halo, sub_masses, ID_sub, z_sub + + +# Assignment +def assignment(hs_order, rc_order, rs_order, bc_order, bs_order, qu_order, id_order, z_order, + scatter=False): + r'''Function that assigns galaxies to halos based on type and the + quenching function + + This function runs through the halo catalogue, assesses whether it has a + quenched galaxy and if it is a central or satellite, and assigns the next + galaxy in the appropriate list. This therefore creates an ordered array + of halos and galaxies which are assigned to each other. Unassigned halos + and galaxies are deleted or ignored. + + Parameters + ----------- + hs_order : (nm, ) array_like + Ordered (most massive first) array of central and subhalo masses, + in units of solar mass + rc_order : (ng1, ) array_like + Ordered (most massive first) array of red central galaxy masses, + in units of solar mass + rs_order : (ng2, ) array_like + Ordered (most massive first) array of red satellite galaxy masses, + in units of solar mass + bc_order : (ng3, ) array_like + Ordered (most massive first) array of blue central galaxy masses, + in units of solar mass + bs_order : (ng4, ) array_like + Ordered (most massive first) array of blue satellite galaxy masses, + in units of solar mass + qu_order : (nm, ) array_like + Truth array of which halos should be assigned quenched (red) galaxies + id_order : (nm, ) array_like + Array of IDs to determine if the halo is central or satellite + z_order : (nm, ) array_like + Redshifts of generated halos + scatter : Boolean + True/false if the galaxies lists have undergone scatter + + Returns + -------- + hs_fin: (nh, ) array_like + List of assigned halo masses, in units of solar masses + gal_fin: (nh, ) array_like + List of assigned galaxies, in units of solar mass + id_fin: (nh, ) array_like + List of ID values for halos + z_fin: (nh, ) array_like + List of redshifts for halos + gal_type_fin: (nh, ) array_like + List of assigned galaxy types: 'red central', 'red satellite', + 'blue central' and 'blue satellite' + + Examples + --------- + >>> from skypy.halos import sham + >>> from skypy.galaxies import schechter_smf + >>> from astropy.cosmology import FlatLambdaCDM + + This example generates the required catalogues, assigns which halos will + be given a quenched galaxy and then assigns them. Assume the parameters + are the same as previous examples or come from the same work. + + >>> #Generate the catalogues + >>> halo_cat, h_z = sham.run_file('halo.yaml', 'halo', 'mass', 'z') + >>> a, b, c, x = 1.91, 0.39, 0.1, 3 + >>> ID_halo, sub_masses, ID_sub, z_sub = sham.gen_sub_cat(halo_cat, h_z, + ... a, b, c, x) + >>> #Sky parameters + >>> cosmology = FlatLambdaCDM(H0=70, Om0=0.3, name='FlatLambdaCDM') + >>> z_range = [min(h_z), max(h_z)] + >>> skyarea = 100 + >>> #Galaxy parameters + >>> m_star1, phi_star1, alpha1 = 10**(10.75), 10**(-2.37), -0.18 + >>> m_star2, phi_star2, alpha2 = 10**(10.72), 10**(-2.66), -0.71 + >>> m_star3, phi_star3, alpha3 = 10**(10.59), 10**(-2.52), -1.15 + >>> m_star4, phi_star4, alpha4 = 10**(10.59), 10**(-3.09), -1.31 + >>> min_mass = 10**(7) + >>> max_mass = 10**(14) + >>> #Generate the galaxies + >>> rc_cat = schechter_smf(z_range, m_star1, phi_star1, alpha1, min_mass, + max_mass, skyarea*u.deg**2, cosmology)[1] + >>> rs_cat = schechter_smf(z_range, m_star2, phi_star2, alpha2, min_mass, + max_mass, skyarea*u.deg**2, cosmology)[1] + >>> bc_cat = schechter_smf(z_range, m_star3, phi_star3, alpha3, min_mass, + max_mass, skyarea*u.deg**2, cosmology)[1] + >>> bs_cat = schechter_smf(z_range, m_star4, phi_star4, alpha4, min_mass, + max_mass, skyarea*u.deg**2, cosmology)[1] + >>> #Quench the halos + >>> m_mu, sigma, base = 10**(12), 0.4, 0.4 + >>> h_quench = sham.quenching_funct(halo_cat, m_mu, sigma) + >>> s_quench = sham.quenching_funct(sub_masses, m_mu, sigma, base) + >>> #Order the arrays + >>> halo_subhalo = np.concatenate((halo_cat, sub_masses), axis=0) + >>> ID_list = np.concatenate((ID_halo, ID_sub), axis=0) + >>> z_list = np.concatenate((h_z, z_sub), axis=0) + >>> q_list = np.concatenate((h_quench, s_quench), axis=0) + >>> stack = np.stack((halo_subhalo, ID_list, z_list, q_list), axis=1) + >>> order1 = stack[np.argsort(stack[:,0])] + >>> hs_order = np.flip(order1[:,0]) + >>> id_hs = np.flip(order1[:,1]) + >>> z_hs = np.flip(order1[:,2]) + >>> qu_hs = np.flip(order1[:,3]) + >>> rc_order = np.flip(np.sort(rc_cat)) #Galaxies + >>> rs_order = np.flip(np.sort(rs_cat)) + >>> bc_order = np.flip(np.sort(bc_cat)) + >>> bs_order = np.flip(np.sort(bs_cat)) + >>> #Call function + >>> hs, gal, ID, z, gal_type = sham.assignment(hs_order, rc_order, + ... rs_order, bc_order, + ... bs_order, qu_hs, + ... id_hs, z_hs) + + References + ---------- + .. + ''' + hs_order = np.atleast_1d(hs_order) + rc_order = np.atleast_1d(rc_order) + rs_order = np.atleast_1d(rs_order) + bc_order = np.atleast_1d(bc_order) + bs_order = np.atleast_1d(bs_order) + qu_order = np.atleast_1d(qu_order) + id_order = np.atleast_1d(id_order) + z_order = np.atleast_1d(z_order) + + # Order check + hs_order_check = np.diff(hs_order) + rc_order_check = np.diff(rc_order) + rs_order_check = np.diff(rs_order) + bc_order_check = np.diff(bc_order) + bs_order_check = np.diff(bs_order) + + # Shape check + hs_shape = hs_order.shape + qu_shape = qu_order.shape + id_shape = id_order.shape + z_shape = z_order.shape + + # Errors + if ((hs_order <= 0)).any(): + raise Exception('Halo masses must be positive and non-zero') + if ((rc_order <= 0)).any() or ((rs_order <= 0)).any(): + raise Exception('Galaxy masses must be positive and non-zero') + if ((bc_order <= 0)).any() or ((bs_order <= 0)).any(): + raise Exception('Galaxy masses must be positive and non-zero') + if ((hs_order_check > 0)).all(): + hs_order = np.flip(hs_order) + raise Warning('Halo masses were in the wrong order and have been correct') + elif ((hs_order_check > 0)).any(): + raise Exception('Halos are not in a sorted order') + if hs_shape != qu_shape or qu_shape != id_shape or id_shape != z_shape: + raise Exception('All arrays pertaining to halos must be the same shape') + + # Only check ordering if not scattered + if not scatter: + if ((rc_order_check > 0)).all(): + rc_order = np.flip(rc_order) + raise Warning('Red central galaxies were in wrong order now correct') + elif ((rc_order_check > 0)).any(): + raise Exception('Red central galaxies are not in a sorted order') + if ((rs_order_check > 0)).all(): + rs_order = np.flip(rs_order) + raise Warning('Red satellite galaxies were in wrong order now correct') + elif ((rs_order_check > 0)).any(): + raise Exception('Red satellite galaxies are not in a sorted order') + if ((bc_order_check > 0)).all(): + bc_order = np.flip(bc_order) + raise Warning('Blue central galaxies were in wrong order and now correct') + elif ((bc_order_check > 0)).any(): + raise Exception('Blue central galaxies are not in a sorted order') + if ((bs_order_check > 0)).all(): + bs_order = np.flip(bs_order) + raise Warning('Blue satellite galaxies were in wrong order and now correct') + elif ((bs_order_check > 0)).any(): + raise Exception('Blue satellite galaxies are not in a sorted order') + + # Assign galaxies to halos + del_ele = [] # Elements to delete later + gal_assigned = [] # Assigned galaxies + gal_type_A = [] # Population assigned + + # Total number of galaxies + gal_num = len(rc_order) + len(rs_order) + len(bc_order) + len(bs_order) + + rc_counter = 0 # Counters for each population + rs_counter = 0 + bc_counter = 0 + bs_counter = 0 + + for ii in range(len(hs_order)): + qu = qu_order[ii] + total_counter = rc_counter + rs_counter + bc_counter + bs_counter + ID_A = id_order[ii] + + if total_counter == gal_num: # All galaxies assigned + del_array = np.arange(ii, len(hs_order)) + del_ele.extend(del_array) + break + + if qu == 1: # Halo is quenched + if ID_A < 0 and rc_counter != len(rc_order): # Halo assigned mass quenched + gal_assigned.append(rc_order[rc_counter]) + gal_type_A.append('red central') + rc_counter += 1 + + elif ID_A > 0 and rs_counter != len(rs_order): # Subhalo assigned environment quenched + gal_assigned.append(rs_order[rs_counter]) + gal_type_A.append('red satellite') + rs_counter += 1 + + else: # No red to assign + del_ele.append(ii) # Delete unassigned halos + + else: # Halo not quenched + if ID_A < 0 and bc_counter != len(bc_order): # Halo assigned blue central + gal_assigned.append(bc_order[bc_counter]) + gal_type_A.append('blue central') + bc_counter += 1 + + elif ID_A > 0 and bs_counter != len(bs_order): # Subhalo assigned blue satellite + gal_assigned.append(bs_order[bs_counter]) + gal_type_A.append('blue satellite') + bs_counter += 1 + + else: # No blue to assign + del_ele.append(ii) + + # Delete and array final lists + hs_fin = np.delete(hs_order, del_ele) + id_fin = np.delete(id_order, del_ele) + z_fin = np.delete(z_order, del_ele) + gal_fin = np.array(gal_assigned) + gal_type_fin = np.array(gal_type_A) + + return hs_fin, gal_fin, id_fin, z_fin, gal_type_fin + + +def scatter_proxy(gal): + r'''Function to add scatter to galaxies using a proxy mass + + This function takes the generated galaxy masses from the catalogues + and generates a gaussian proxy mass for each, then reorders based on + this proxy mass. This adds scatter to the galaxies. + + Parameters + ----------- + gal : (nm, ) array_like + List of galaxy masses to be scattered, in units of solar mass + + Returns + -------- + new_gal : (nm, ) array_like + Reordered galaxy masses according to generated proxy mass, highest + to lowest mass, in units of solar mass + + Examples + --------- + >>> from skypy.halos import sham + + This example generates a galaxy catalogue and then reorders it + + >>> gal_cat = sham.run_file('galaxies.yaml', galaxy, mass) + >>> new_gal = sham.scatter_proxy(gal_cat) + + References + ---------- + .. + ''' + # Errors + gal = np.atleast_1d(gal) + if ((gal <= 0)).any(): + raise Exception('Galaxy masses must be positive and non-zero') + + # Draw one sample for each galaxy + # Use galaxy mass as mean and sigma of one order less + proxy_gal = np.random.normal(loc=gal, scale=gal/10) + + # Order by proxy mass + proxy_stack = np.stack((proxy_gal, gal), axis=1) + gal_stack = proxy_stack[np.argsort(proxy_stack[:, 0])] + return np.flip(gal_stack[:, 1]) # Return new galaxies, highest mass first + + +# Separate populations +def sham_plots(hs_fin, gal_fin, gal_type_fin, print_out=False): + r'''Function that pulls assigned galaxies together into groups for + plotting + + This function takes in the final assigned halos and galaxies and splits + them into each type of galaxy for plotting and comparing to data. It also + computes the separate central and satellite populations with no colour + split. + + Parameters + ----------- + hs_fin : (nm, ) array_like + List of assigned halos and subhalos, in units of solar mass. The + subhalos are assumed to be their current stripped mass + gal_fin : (nm, ) array_like + List of assigned galaxies, in units of solar mass + gal_type_fin : (nm, ) array_like + List of assigned galaxies types with the tags + 'red central', 'red satellite', 'blue central', 'blue satellite' + print_out : Boolean + True/false whether to output print statements for progress and + timings + + Returns + -------- + sham_rc: (nh1, 2) array_like + Stacked array of halos ([:,0]) and galaxies ([:,1]) containing + red centrals + sham_rs: (nh2, 2) array_like + Stacked array of halos ([:,0]) and galaxies ([:,1]) containing + red satellites + sham_bc: (nh3, 2) array_like + Stacked array of halos ([:,0]) and galaxies ([:,1]) containing + blue centrals + sham_bs: (nh4, 2) array_like + Stacked array of halos ([:,0]) and galaxies ([:,1]) containing + blue satellites + sham_cen: (nh5, 2) array_like + Stacked array of halos ([:,0]) and galaxies ([:,1]) containing + centrals + sham_sub: (nh6, 2) array_like + Stacked array of halos ([:,0]) and galaxies ([:,1]) containing + satellites + + Examples + --------- + >>> from skypy.halos import sham + + This example uses assigned halos and galaxies and then finds the + separated SHAM arrays. + + >>> hs, gal, id, z, gal_type = sham.assignment(hs_order, rc_order, + ... rs_order, bc_order, + ... bs_order, qu_order, + ... id_order, z_order) + >>> rc, rs, bc, bs, cen, sub = sham.sham_plots(hs, gal, gal_type) + + References + ---------- + .. + ''' + # Errors + if hs_fin.shape != gal_fin.shape or gal_fin.shape != gal_type_fin.shape: + raise Exception('All arrays must be the same shape') + if ((hs_fin <= 0)).any(): + raise Exception('Halo masses must be positive and non-zero') + if ((gal_fin <= 0)).any(): + raise Exception('Galaxy masses must be positive and non-zero') + + rc_dots = np.where(gal_type_fin == 'red central')[0] + rs_dots = np.where(gal_type_fin == 'red satellite')[0] + bc_dots = np.where(gal_type_fin == 'blue central')[0] + bs_dots = np.where(gal_type_fin == 'blue satellite')[0] + + # SHAMs by galaxy population + sham_rc = np.stack((hs_fin[rc_dots], gal_fin[rc_dots]), axis=1) + sham_rs = np.stack((hs_fin[rs_dots], gal_fin[rs_dots]), axis=1) + sham_bc = np.stack((hs_fin[bc_dots], gal_fin[bc_dots]), axis=1) + sham_bs = np.stack((hs_fin[bs_dots], gal_fin[bs_dots]), axis=1) + + # Combine both sets of centrals + halo_all = np.concatenate((sham_rc[:, 0], sham_bc[:, 0]), axis=0) # All halos and galaxies + gals_all = np.concatenate((sham_rc[:, 1], sham_bc[:, 1]), axis=0) + + sham_concat = np.stack((halo_all, gals_all), axis=1) + sham_order = sham_concat[np.argsort(sham_concat[:, 0])] # Order by halo mass + sham_cen = np.flip(sham_order, 0) + + # Combine both sets of satellites + halo_sub = np.concatenate((sham_rs[:, 0], sham_bs[:, 0]), axis=0) # All halos and galaxies + gals_sub = np.concatenate((sham_rs[:, 1], sham_bs[:, 1]), axis=0) + + sham_concats = np.stack((halo_sub, gals_sub), axis=1) + sham_orders = sham_concats[np.argsort(sham_concats[:, 0])] # Order by halo mass + sham_sub = np.flip(sham_orders, 0) + + if print_out: + print('Created SHAMs') + + return sham_rc, sham_rs, sham_bc, sham_bs, sham_cen, sham_sub + + +# Called SHAM function +def run_sham(h_file, gal_param, cosmology, z_range, skyarea, qu_param, + sub_param=[1.91, 0.39, 0.1, 3], gal_max_h=10**(14), gal_max_s=10**(13), + print_out=False, run_anyway=False, scatter_prox=False): + r'''Function that takes all inputs for the halos and galaxies and runs + a SHAM over them + + This function runs all the functions required to make a complete SHAM. + It generates the catalogues from the input values, quenches the halos + and assigns the galaxies. Then it outputs the sorted SHAM arrays by + galaxy type. + + Parameters + ----------- + h_file : str + String for YAML file of halos to generate catalogue + gal_param : (4, 3) array_like + Parameters (M* (solar mass), phi* (Mpc^{-3}), alpha) for Schechter + functions of four galaxy types (red centrals, red satellites, + blue centrals, blue satellites) + cosmology : Cosmology + Astropy cosmology object, must specify H0 (km/s/Mpc), Om0 and astropy + name. Should be the same as specified in the halo YAML file + z_range : (2, ) array_like + Minimum and maximum redshifts to generate for galaxies. Should be the + same as specified in the halo YAML file + skyarea : float + Sky area galaxies are 'observed' from, units of deg^2. Should be the + same as specified in the halo YAML file + qu_param : (3, ) array_like + The mean, standard deviation and baseline of the error function + to quench the halos (see 'quenching_funct' for more detail) + sub_param : (4, ) array_like + alpha, beta, gamma and x parameters for subhalo generation. Defaults + from [1]_ (see 'gen_sub_cat' for more detail) + gal_max_h : float + Maximum mass central galaxy to generate, units of solar mass + gal_max_s : float + Maximum mass satellite galaxy to generate, units of solar mass + print_out : Boolean + True/false whether to output print statements for progress and timings + run_anyway : Boolean + True/false whether to run the SHAM code even if the galaxy abundance + does not match the number of halos + scatter_prox : Boolean + True/false whether to scatter the galaxies by a proxy mass when the + catalogues are generated + + Returns + -------- + sham_dict: (nm, 5) dictionary + A dictionary containing the necessary outputs stored as numpy arrays. + - 'Halo mass' parent and subhalo masses, in units of solar mass, + the subhalos are their present stripped mass + - 'Galaxy mass' assigned galaxy masses, in units of solar mass + - 'Galaxy type' type of galaxy assigned, written as 'red central', + 'red satellite', 'blue central' and 'blue satellite' + - 'ID value' ID number of halos and subhalos, negative for halos, + positive for subhalos + - 'Redshift' Redshift of halos/galaxies + + Examples + --------- + >>> from skypy.halos import sham + >>> from astropy.cosmology import FlatLambdaCDM + >>> import numpy as np + + This example shows the format for each of the inputs and receives the + outputs. Subhalo parameters are from [1]_, galaxy parameters are from [2]_ + and quenching parameters are approximately from [3]_ + + >>> #Parameters + >>> h_file = 'halo.yaml' + >>> rc_p = [10**(10.75), 10**(-2.37), -0.18] + >>> rs_p = [10**(10.72), 10**(-2.66), -0.71] + >>> bc_p = [10**(10.59), 10**(-2.52), -1.15] + >>> bs_p = [10**(10.59), 10**(-3.09), -1.31] + >>> gal_param = [rc_p, rs_p, bc_p, bs_p] + >>> cosmology = FlatLambdaCDM(H0=70, Om0=0.3, name='FlatLambdaCDM') + >>> z_range = np.array([0.01, 0.1]) + >>> skyarea = 600. + >>> qu = np.array([10**(11.9), 0.4, 0.5]) + >>> #Run function + >>> sham_dict = run_sham(h_file, gal_param, cosmology, z_range, skyarea, + ... qu, sub_param = [1.91, 0.39, 0.1, 3], + ... gal_max_h = 10**(14), gal_max_s = 10**(13), + ... print_out=False, run_anyway=True, + scatter_prox=True) + + References + ---------- + .. [1] Vale A., Ostriker J. P., 2004, Monthly Notices of the Royal + Astronomical Society, 353, 189 + .. [2] Weigel A. K., Schawinski K., Bruderer C., 2016, Monthly Notices of the + Royal Astronomical Society, 459, 2150 + .. [3] Peng Y.-j., et al., 2010, Astrophysical Journal, 721, 193 + ''' + sham_st = time() + + gal_param = np.atleast_1d(gal_param) + qu_param = np.atleast_1d(qu_param) + + # Check that all inputs are of the correct type and size + if type(h_file) is not str: + raise Exception('Halo YAML file must be provided as a string') + if gal_param.shape != (4, 3): + if gal_param.shape[0] != 4: + raise Exception('The wrong number of galaxies are in galaxy parameters') + elif gal_param.shape[1] != 3: + raise Exception('The wrong number of galaxy parameters have been provided') + else: + raise Exception('Supplied galaxy parameters are not the correct shape') + if qu_param.shape != (3,): + raise Exception('Provided incorrect number of quenching parameters') + if z_range.shape != (2, ): + raise Exception('The wrong number of redshifts were given') + if z_range[0] < 0 or z_range[1] < 0: + raise Exception('Redshift cannot be negative') + if z_range[1] <= z_range[0]: + raise Exception('The second redshift should be more than the first') + + # Check galaxy parameters are correct sign + if np.any(np.array([gal_param[0][0], gal_param[1][0], gal_param[2][0], gal_param[3][0]]) <= 0): + raise Warning('M* values must be positive and non-zero') + if np.any(np.array([gal_param[0][1], gal_param[1][1], gal_param[2][1], gal_param[3][1]]) <= 0): + raise Exception('phi* values must be positive and non-zero') + if np.any(np.array([gal_param[0][2], gal_param[1][2], gal_param[2][2], gal_param[3][2]]) > 0): + raise Warning('Galaxy Schechter function alphas must be < 0') + if skyarea <= 0: + raise Exception('The skyarea must be a positive non-zero number') + if cosmology.name is None: + raise Exception('Cosmology object must have an astropy cosmology name') + + # Generate parent halos from YAML file + # TODO remove hard coded table/variable names + h_st = time() + parent_halo, z_halo = run_file(h_file, 'halo', 'mass', 'z') + if print_out: + print('Halo catalogues generated in', round((time() - h_st), 2), 's') + + # Generate subhalos and IDs + sub_alpha = sub_param[0] # Vale and Ostriker 2004 mostly + sub_beta = sub_param[1] + sub_gamma = sub_param[2] + sub_x = sub_param[3] + + sub_tim = time() + ID_halo, subhalo_m, ID_sub, z_sub = gen_sub_cat(parent_halo, z_halo, sub_alpha, + sub_beta, sub_gamma, sub_x) + + if print_out: + print('Generated subhalos and IDs in', round((time() - sub_tim), 2), 's') + print('') + + # Quench halos and subhalos + M_mu = qu_param[0] # Parameters + sigma = qu_param[1] + baseline = qu_param[2] + + h_quench = quenching_funct(parent_halo, M_mu, sigma, 0) # Quenched halos + sub_quench = quenching_funct(subhalo_m, M_mu, sigma, baseline) # Quenched subhalos + + # Galaxy Schechter function parameters + rc_param = gal_param[0] + rs_param = gal_param[1] + bc_param = gal_param[2] + bs_param = gal_param[3] + + # Number of halos for each population + rc_halo = len(np.where(h_quench == 1)[0]) + rs_halo = len(np.where(sub_quench == 1)[0]) + bc_halo = len(np.where(h_quench == 0)[0]) + bs_halo = len(np.where(sub_quench == 0)[0]) + + # Find galaxy mass range (m_star, phi_star, alpha, tag) + # TODO Add a way to import a look up table + range1 = time() + rc_min = find_min(rc_param[0], rc_param[1], rc_param[2], cosmology, z_range, skyarea, + gal_max_h, rc_halo, print_out, run_anyway) + if print_out: + print('Red central log(minimum mass)', np.log10(rc_min), 'in', + round((time() - range1), 4), 's') + + range2 = time() + rs_min = find_min(rs_param[0], rs_param[1], rs_param[2], cosmology, z_range, skyarea, + gal_max_s, rs_halo, print_out, run_anyway) + if print_out: + print('Red satellite log(minimum mass)', np.log10(rs_min), 'in', + round((time() - range2), 4), 's') + + range3 = time() + bc_min = find_min(bc_param[0], bc_param[1], bc_param[2], cosmology, z_range, skyarea, + gal_max_h, bc_halo, print_out, run_anyway) + if print_out: + print('Blue central log(minimum mass)', np.log10(bc_min), 'in', + round((time() - range3), 4), 's') + + range4 = time() + bs_min = find_min(bs_param[0], bs_param[1], bs_param[2], cosmology, z_range, skyarea, + gal_max_s, bs_halo, print_out, run_anyway) + if print_out: + print('Blue satellite log(minimum mass)', np.log10(bs_min), 'in', + round((time() - range4), 4), 's') + print('') + + # Get a catalogue for each population + cat_time = time() + + redshift = np.linspace(z_range[0], z_range[1], 100) # Redshift volume + + rc_cat = schechter_smf(redshift, rc_param[0], rc_param[1], rc_param[2], rc_min, + gal_max_h, skyarea*u.deg**2, cosmology)[1] + rs_cat = schechter_smf(redshift, rs_param[0], rs_param[1], rs_param[2], rs_min, + gal_max_s, skyarea*u.deg**2, cosmology)[1] + bc_cat = schechter_smf(redshift, bc_param[0], bc_param[1], bc_param[2], bc_min, + gal_max_h, skyarea*u.deg**2, cosmology)[1] + bs_cat = schechter_smf(redshift, bs_param[0], bs_param[1], bs_param[2], bs_min, + gal_max_s, skyarea*u.deg**2, cosmology)[1] + + if scatter_prox: + if print_out: + print('Using scattering') + print('') + rc_order = scatter_proxy(rc_cat) + rs_order = scatter_proxy(rs_cat) + bc_order = scatter_proxy(bc_cat) + bs_order = scatter_proxy(bs_cat) + + if print_out: + print('Galaxy catalogues generated in', round((time() - cat_time), 2), 's') + + # Order and process DM and galaxies + # Concatenate halos and subhalos + halo_subhalo = np.concatenate((parent_halo, subhalo_m), axis=0) + ID_list = np.concatenate((ID_halo, ID_sub), axis=0) + z_list = np.concatenate((z_halo, z_sub), axis=0) + q_list = np.concatenate((h_quench, sub_quench), axis=0) + + # Sort lists by halo mass + stack = np.stack((halo_subhalo, ID_list, z_list, q_list), axis=1) + order1 = stack[np.argsort(stack[:, 0])] + hs_order = np.flip(order1[:, 0]) + id_hs = np.flip(order1[:, 1]) + z_hs = np.flip(order1[:, 2]) + qu_hs = np.flip(order1[:, 3]) + + # List galaxies if not already sorted by scattering + if not scatter_prox: + if print_out: + print('No scatter applied') + print('') + rc_order = np.flip(np.sort(rc_cat)) + rs_order = np.flip(np.sort(rs_cat)) + bc_order = np.flip(np.sort(bc_cat)) + bs_order = np.flip(np.sort(bs_cat)) + + # Assignment of galaxies + assign_st = time() + hs_fin, gal_fin, id_fin, z_fin, gal_type_fin = assignment(hs_order, rc_order, rs_order, + bc_order, bs_order, qu_hs, id_hs, + z_hs, scatter=scatter_prox) + + if print_out: + print('Galaxies assigned in', round((time() - assign_st), 2), 's') + print('') + + # Strip subhalos + sub_loc = np.where(id_fin > 0) + hs_fin[sub_loc] = hs_fin[sub_loc]/sub_x + + # Create output dictionary + sham_dict = { + 'Halo mass': hs_fin, + 'Galaxy mass': gal_fin, + 'Galaxy type': gal_type_fin, + 'ID value': id_fin, + 'Redshift': z_fin + } + + if print_out: + print('SHAM run in', round(((time() - sham_st)/60), 2), 'min') + + return sham_dict diff --git a/skypy/halos/test_gal.yaml b/skypy/halos/test_gal.yaml new file mode 100644 index 000000000..e1ce329d5 --- /dev/null +++ b/skypy/halos/test_gal.yaml @@ -0,0 +1,20 @@ +m_star: !numpy.power [10, 10.58] +phi_star: !numpy.power [10, -2.77] +alpha_val: -0.33 +m_min: !numpy.power [10, 7.5] +m_max: !numpy.power [10, 14.0] +sky_area: 600.0 deg2 +z_range: !numpy.linspace [0.01, 0.1, 100] +cosmology: !astropy.cosmology.FlatLambdaCDM + H0: 70.0 + Om0: 0.3 +tables: + galaxy: + z, sm: !skypy.galaxies.schechter_smf + redshift: $z_range + m_star: $m_star + phi_star: $phi_star + alpha: $alpha_val + m_min: $m_min + m_max: $m_max + sky_area: $sky_area \ No newline at end of file diff --git a/skypy/halos/tests/test_gal.yaml b/skypy/halos/tests/test_gal.yaml new file mode 100644 index 000000000..8d422900c --- /dev/null +++ b/skypy/halos/tests/test_gal.yaml @@ -0,0 +1,20 @@ +m_star: !numpy.power [10, 10.75] +phi_star: !numpy.power [10, -2.37] +alpha_val: -0.18 +m_min: !numpy.power [10, 7.0] +m_max: !numpy.power [10, 14.0] +sky_area: 600.0 deg2 +z_range: !numpy.linspace [0.01, 0.1, 100] +cosmology: !astropy.cosmology.FlatLambdaCDM + H0: 70.0 + Om0: 0.3 +tables: + galaxy: + z, sm: !skypy.galaxies.schechter_smf + redshift: $z_range + m_star: $m_star + phi_star: $phi_star + alpha: $alpha_val + m_min: $m_min + m_max: $m_max + sky_area: $sky_area diff --git a/skypy/halos/tests/test_halo.yaml b/skypy/halos/tests/test_halo.yaml new file mode 100644 index 000000000..4288ad416 --- /dev/null +++ b/skypy/halos/tests/test_halo.yaml @@ -0,0 +1,26 @@ +parameters: + model: 'sheth99' + mdef: 'fof' + m_min: 1.e+9 + m_max: 1.e+15 + sky_area: 600. deg2 # Same as galaxies + sigma_8: 0.8 + ns: 0.96 +cosmology: !astropy.cosmology.FlatLambdaCDM + H0: 70 + Om0: 0.3 + name: 'FlatLambdaCDM' # Requires user defined name + Ob0: 0.045 + Tcmb0: 2.7 +z_range: !numpy.linspace [0.01, 0.1, 100] # Same as galaxies +tables: + halo: + z, mass: !skypy.halos.colossus_mf + redshift: $z_range + model: $model + mdef: $mdef + m_min: $m_min + m_max: $m_max + sky_area: $sky_area + sigma8: $sigma_8 + ns: $ns \ No newline at end of file diff --git a/skypy/halos/tests/test_sham.py b/skypy/halos/tests/test_sham.py new file mode 100644 index 000000000..cbe7a0eb7 --- /dev/null +++ b/skypy/halos/tests/test_sham.py @@ -0,0 +1,1043 @@ +# Imports +import numpy as np +import pytest +from pytest import approx +from skypy.halos._colossus import HAS_COLOSSUS + + +@pytest.mark.skipif(not HAS_COLOSSUS, reason='test requires colossus') +@pytest.mark.flaky +def test_quenching_funct(): + # Run a test on the quenching function to check it follows the expected function values + + # Create halo catalogue + from astropy.cosmology import WMAP9 # Cannot be FlatLambdaCDM + from skypy.halos.mass import colossus_mass_sampler + from skypy.halos.sham import quenching_funct + m_min, m_max, size = 1e+10, 1e+16, 1000000 + parent_halo = colossus_mass_sampler(redshift=0.1, model='sheth99', + mdef='fof', m_min=m_min, m_max=m_max, + cosmology=WMAP9, sigma8=0.8, ns=1., + size=size, resolution=1000) + + # Parameters + m_mu, sigma, base = 10**(12.1), 0.4, 0.2 + + h_quench = quenching_funct(parent_halo, m_mu, sigma, base) # Quenched halos + + assert ((h_quench == 0) | (h_quench == 1)).all() # Check result is binary + assert h_quench.shape == parent_halo.shape # Check the shapes are the same + + # Make histogram of results + qu_stack = np.stack((parent_halo, h_quench), axis=1) + order = qu_stack[np.argsort(qu_stack[:, 0])] + ha_order = order[:, 0] + qu_order = order[:, 1] + + mass_bins = np.geomspace(min(parent_halo), max(parent_halo), 50) + mass_mid = (mass_bins[:-1] + mass_bins[1:])/2 + + fract = [] + for ii in range(0, len(mass_bins) - 1): + l, r = mass_bins[ii], mass_bins[ii+1] + # Halos in bin + h_bin = np.where(ha_order[np.where(ha_order < r)] >= l)[0] + q_bin = qu_order[h_bin] + + if len(q_bin) != 0: + fract.append(len(np.where(q_bin == 1)[0])/len(q_bin)) + else: + fract.append(0) # Some bins near end won't have any data in + + fract = np.array(fract) + + # Original PDF + from scipy.special import erf + mass = np.geomspace(min(parent_halo), max(parent_halo), 100) + calc = np.log10(mass/m_mu)/sigma + old_prob = 0.5 * (1.0 + erf(calc/np.sqrt(2))) # Base error funct, 0 -> 1 + add_val = base/(1-base) + prob_sub = old_prob + add_val # Add value to plot to get baseline + prob_array = prob_sub/prob_sub[-1] # Normalise, base -> 1 + + # Interp to get same masses + prob = np.interp(mass_mid, mass, prob_array) + + # Compare sampled to expected distribution + dist = (prob - fract)**2 + dist = np.delete(dist, np.where(dist > 10**(-1))[0]) # Delete values where there was no data + assert sum(dist)/len(dist) < 1e-4 # Is average difference small enough + + # Check the first value is around baseline + assert fract[0] == approx(base, rel=10**(-1)) + + # Check the mean is correct + assert np.interp(m_mu, mass_mid, fract) == approx((1 + base)/2, rel=10**(-1)) + + # Check the end is correct + fract_one = fract[np.where(fract > 0)] # Remove any zero values ie where there is no data + assert fract_one[-1] == approx(1, rel=10**(-1)) + + # Check the errors trigger + with pytest.raises(Exception) as excinfo: + m_mu, sigma, base = -10**(12), 0.4, 0.2 + quenching_funct(parent_halo, m_mu, sigma, base) + assert str(excinfo.value) == 'M_mu cannot be negative or 0 solar masses' + with pytest.raises(Exception) as excinfo: + m_mu, sigma, base = 0, 0.4, 0.2 + quenching_funct(parent_halo, m_mu, sigma, base) + assert str(excinfo.value) == 'M_mu cannot be negative or 0 solar masses' + with pytest.raises(Exception) as excinfo: + m_mu, sigma, base = 10**(12), -0.3, 0.2 + quenching_funct(parent_halo, m_mu, sigma, base) + assert str(excinfo.value) == 'sigma cannot be negative or 0' + with pytest.raises(Exception) as excinfo: + m_mu, sigma, base = 10**(12), 0, 0.2 + quenching_funct(parent_halo, m_mu, sigma, base) + assert str(excinfo.value) == 'sigma cannot be negative or 0' + with pytest.raises(Exception) as excinfo: + m_mu, sigma, base = 10**(12), 0.3, -0.7 + quenching_funct(parent_halo, m_mu, sigma, base) + assert str(excinfo.value) == 'baseline cannot be negative, set_baseline = -baseline' + with pytest.raises(Exception) as excinfo: + m_mu, sigma, base = 10**(12), 0.3, 1.3 + quenching_funct(parent_halo, m_mu, sigma, base) + assert str(excinfo.value) == 'baseline cannot be more than 1' + with pytest.raises(Exception) as excinfo: + m_mu, sigma, base = 10**(12), 0.3, 0.3 + quenching_funct(10**(11), m_mu, sigma, base) + assert str(excinfo.value) == 'Mass must be array-like' + + +@pytest.mark.flaky +def test_find_min(): + # Run a test on the galaxy's find_min function to check it gives expected values + + from scipy.integrate import trapezoid as trap + from astropy import units as u + from astropy.cosmology import FlatLambdaCDM + from skypy.halos.sham import find_min + + # Parameters + m_star, phi_star, alpha = 10**(10.58), 10**(-2.77), -0.33 + cosmology = FlatLambdaCDM(H0=70, Om0=0.3, name='FlatLambdaCDM') + z_range = [0.01, 0.1] + skyarea = 600 + max_mass = 10**(14) + no_halo_y, no_halo_min, no_halo_max = 10000, 600, 50000 + + # Run function when interpolation should work + min_mass = find_min(m_star, phi_star, alpha, cosmology, z_range, skyarea, max_mass, no_halo_y) + + # Check type and expected number + assert type(min_mass) is np.float64 + assert np.log10(min_mass) == approx(7.9976, rel=10**(-4)) + + # Check that the expected integral gives the number of halos + skyarea = skyarea*u.deg**2 + z = np.linspace(z_range[0], z_range[1], 1000) + dVdz = (cosmology.differential_comoving_volume(z)*skyarea).to_value('Mpc3') + + mass_range = np.geomspace(min_mass, max_mass, 1000) + phi_m = phi_star/m_star + m_m = mass_range/m_star + + dndmdV = phi_m*np.e**(-m_m)*(m_m)**alpha # Mass function + dndV = trap(dndmdV, mass_range) + dn = dndV*dVdz + + assert trap(dn, z) == approx(no_halo_y, rel=10**(-1)) + + # Run function with a non-optimal number of halos, check for exceptions + skyarea = 600 + with pytest.raises(Exception) as excinfo: + find_min(m_star, phi_star, alpha, cosmology, z_range, skyarea, max_mass, no_halo_min) + assert str(excinfo.value) == 'Number of halos not within reachable bounds' + + with pytest.raises(Exception) as excinfo: + find_min(m_star, phi_star, alpha, cosmology, z_range, skyarea, max_mass, no_halo_max) + assert str(excinfo.value) == 'Number of halos not within reachable bounds' + + # Check function gives expected output using 'run_anyway' + test1 = find_min(m_star, phi_star, alpha, cosmology, z_range, skyarea, + max_mass, no_halo_min, run_anyway=True) + test2 = find_min(m_star, phi_star, alpha, cosmology, z_range, skyarea, + max_mass, no_halo_max, run_anyway=True) + assert test1 == approx(10**(7)) + assert test2 == approx(10**(7)) + + # Check errors trigger + with pytest.raises(Exception) as excinfo: + m_star, phi_star, alpha = 10**(10.58), 10**(-2.77), -0.33 + cosmology = FlatLambdaCDM(H0=70, Om0=0.3, name='FlatLambdaCDM') + z_range = [0.01, 0.1, 10] + skyarea = 600 + max_mass = 10**(14) + no_halo_y, no_halo_min, no_halo_max = 10000, 600, 50000 + find_min(m_star, phi_star, alpha, cosmology, z_range, skyarea, max_mass, no_halo_y) + assert str(excinfo.value) == 'The wrong number of redshifts were given' + + with pytest.raises(Exception) as excinfo: + m_star, phi_star, alpha = 10**(10.58), 10**(-2.77), -0.33 + cosmology = FlatLambdaCDM(H0=70, Om0=0.3, name='FlatLambdaCDM') + z_range = [-0.01, 0.1] + skyarea = 600 + max_mass = 10**(14) + no_halo_y, no_halo_min, no_halo_max = 10000, 600, 50000 + find_min(m_star, phi_star, alpha, cosmology, z_range, skyarea, max_mass, no_halo_y) + assert str(excinfo.value) == 'Redshift cannot be negative' + + with pytest.raises(Exception) as excinfo: + m_star, phi_star, alpha = 10**(10.58), 10**(-2.77), -0.33 + cosmology = FlatLambdaCDM(H0=70, Om0=0.3, name='FlatLambdaCDM') + z_range = [0.01, -0.1] + skyarea = 600 + max_mass = 10**(14) + no_halo_y, no_halo_min, no_halo_max = 10000, 600, 50000 + find_min(m_star, phi_star, alpha, cosmology, z_range, skyarea, max_mass, no_halo_y) + assert str(excinfo.value) == 'Redshift cannot be negative' + + with pytest.raises(Exception) as excinfo: + m_star, phi_star, alpha = 10**(10.58), 10**(-2.77), -0.33 + cosmology = FlatLambdaCDM(H0=70, Om0=0.3, name='FlatLambdaCDM') + z_range = [-0.01, -0.1] + skyarea = 600 + max_mass = 10**(14) + no_halo_y, no_halo_min, no_halo_max = 10000, 600, 50000 + find_min(m_star, phi_star, alpha, cosmology, z_range, skyarea, max_mass, no_halo_y) + assert str(excinfo.value) == 'Redshift cannot be negative' + + with pytest.raises(Exception) as excinfo: + m_star, phi_star, alpha = 10**(10.58), 10**(-2.77), -0.33 + cosmology = FlatLambdaCDM(H0=70, Om0=0.3, name='FlatLambdaCDM') + z_range = [0.1, 0.01] + skyarea = 600 + max_mass = 10**(14) + no_halo_y, no_halo_min, no_halo_max = 10000, 600, 50000 + find_min(m_star, phi_star, alpha, cosmology, z_range, skyarea, max_mass, no_halo_y) + assert str(excinfo.value) == 'The second redshift should be more than the first' + + with pytest.raises(Exception) as excinfo: + m_star, phi_star, alpha = 10**(10.58), 10**(-2.77), 0.33 + cosmology = FlatLambdaCDM(H0=70, Om0=0.3, name='FlatLambdaCDM') + z_range = [0.01, 0.1] + skyarea = 600 + max_mass = 10**(14) + no_halo_y, no_halo_min, no_halo_max = 10000, 600, 50000 + find_min(m_star, phi_star, alpha, cosmology, z_range, skyarea, max_mass, no_halo_y) + assert str(excinfo.value) == 'Schechter function defined so alpha < 0, set_alpha = -alpha' + + with pytest.raises(Exception) as excinfo: + m_star, phi_star, alpha = -10**(10.58), 10**(-2.77), -0.33 + cosmology = FlatLambdaCDM(H0=70, Om0=0.3, name='FlatLambdaCDM') + z_range = [0.01, 0.1] + skyarea = 600 + max_mass = 10**(14) + no_halo_y, no_halo_min, no_halo_max = 10000, 600, 50000 + find_min(m_star, phi_star, alpha, cosmology, z_range, skyarea, max_mass, no_halo_y) + assert str(excinfo.value) == 'M* and phi* must be positive and non-zero numbers' + + with pytest.raises(Exception) as excinfo: + m_star, phi_star, alpha = -10**(10.58), 0, -0.33 + cosmology = FlatLambdaCDM(H0=70, Om0=0.3, name='FlatLambdaCDM') + z_range = [0.01, 0.1] + skyarea = 600 + max_mass = 10**(14) + no_halo_y, no_halo_min, no_halo_max = 10000, 600, 50000 + find_min(m_star, phi_star, alpha, cosmology, z_range, skyarea, max_mass, no_halo_y) + assert str(excinfo.value) == 'M* and phi* must be positive and non-zero numbers' + + with pytest.raises(Exception) as excinfo: + m_star, phi_star, alpha = 10**(10.58), 10**(-2.77), -0.33 + cosmology = FlatLambdaCDM(H0=70, Om0=0.3, name='FlatLambdaCDM') + z_range = [0.01, 0.1] + skyarea = -60 + max_mass = 10**(14) + no_halo_y, no_halo_min, no_halo_max = 10000, 600, 50000 + find_min(m_star, phi_star, alpha, cosmology, z_range, skyarea, max_mass, no_halo_y) + assert str(excinfo.value) == 'The skyarea must be a positive non-zero number' + + with pytest.raises(Exception) as excinfo: + m_star, phi_star, alpha = 10**(10.58), 10**(-2.77), -0.33 + cosmology = FlatLambdaCDM(H0=70, Om0=0.3, name='FlatLambdaCDM') + z_range = [0.1, 0.01] + skyarea = 60 + max_mass = 10**(4) + no_halo_y, no_halo_min, no_halo_max = 10000, 600, 50000 + find_min(m_star, phi_star, alpha, cosmology, z_range, skyarea, max_mass, no_halo_y) + assert str(excinfo.value) == 'The maximum mass must be more than 10^6 solar masses' + + +@pytest.mark.flaky +def test_run_file(): + # Run a test on run_file to check it correctly produces a catalogue from a yaml file + + import astropy + from scipy.integrate import trapezoid as trap + from astropy import units as u + from astropy.cosmology import FlatLambdaCDM + from skypy.halos.sham import run_file + + # Test catalogue is the expected length + m_star, phi_star, alpha = 10**(10.58), 10**(-2.77), -0.33 + cosmology = FlatLambdaCDM(H0=70, Om0=0.3, name='FlatLambdaCDM') + z_range = [0.01, 0.1] + skyarea = 600 + min_mass = 10**(7.5) + max_mass = 10**(14) + + file_test = 'test_gal.yaml' + test_cat, test_z = run_file(file_test, 'galaxy', 'sm', 'z') + + skyarea = skyarea*u.deg**2 + z = np.linspace(z_range[0], z_range[1], 1000) + dVdz = (cosmology.differential_comoving_volume(z)*skyarea).to_value('Mpc3') + + mass_range = np.geomspace(min_mass, max_mass, 1000) + phi_m = phi_star/m_star + m_m = mass_range/m_star + + dndmdV = phi_m*np.e**(-m_m)*(m_m)**alpha # Mass function + dndV = trap(dndmdV, mass_range) + dn = dndV*dVdz + + assert type(test_cat) is astropy.table.column.Column + assert type(test_z) is astropy.table.column.Column + assert test_cat.ndim == 1 + assert test_z.ndim == 1 + + # Length same within Poisson sampling error + assert trap(dn, z) == approx(len(test_cat), rel=np.sqrt(trap(dn, z))) + assert len(test_cat) == len(test_z) + + # Check errors trigger + with pytest.raises(Exception) as excinfo: + run_file(82, 'galaxy', 'sm', 'z') + assert str(excinfo.value) == 'File name must be a string' + + with pytest.raises(Exception) as excinfo: + run_file('GNDN.yaml', 'galaxy', 'sm', 'z') + assert str(excinfo.value) == 'File does not exist' + + +@pytest.mark.skipif(not HAS_COLOSSUS, reason='test requires colossus') +@pytest.mark.flaky +def test_gen_sub_cat(): + # Run a test on generating the subhalos to check it outputs correctly + + from astropy.cosmology import WMAP9 # Cannot be FlatLambdaCDM + from skypy.halos.mass import colossus_mass_sampler + from skypy.halos.sham import gen_sub_cat + + # Parameters + alpha, beta, gamma, x = 1.91, 0.39, 0.1, 3 + + # Catalogue + m_min, m_max, size = 1e+10, 1e+16, 1000 + parent_halo = colossus_mass_sampler(redshift=0.1, model='sheth99', + mdef='fof', m_min=m_min, m_max=m_max, + cosmology=WMAP9, sigma8=0.8, ns=1., + size=size, resolution=1000) + z_halo = np.linspace(0.01, 0.1, len(parent_halo)) + + # Function + ID_halo, sub_masses, ID_sub, z_sub = gen_sub_cat(parent_halo, z_halo, alpha, beta, gamma, x) + + # Check array sets + assert np.all(np.array([ID_halo.ndim, sub_masses.ndim, ID_sub.ndim, z_sub.ndim]) == 1) + assert len(ID_halo) == len(parent_halo) + assert len(sub_masses) == len(ID_sub) + assert len(sub_masses) == len(z_sub) + + # Check the data inside is correct + assert min(sub_masses) >= min(parent_halo) + assert max(sub_masses) <= max(parent_halo)/2 + + assert np.all(np.isin(z_sub, z_halo)) + assert ((ID_halo < 0)).all() + assert ((ID_sub > 0)).all() + + i_type = [] + id_list = np.append(ID_halo, ID_sub) + for ii in id_list: + if type(ii) is np.int64: + i_type.append(True) + else: + i_type.append(False) + assert np.all(i_type) # Check ids are all integers + + # Check a single value works + # Parameters + alpha, beta, gamma, x = 1.91, 0.39, 0.1, 3 + + # Catalogue + parent_halo = 1e+16 + z_halo = 0.01 + + # Function + ID_halo, sub_masses, ID_sub, z_sub = gen_sub_cat(parent_halo, z_halo, alpha, beta, gamma, x) + + assert sub_masses.size # Check array is not empty, i.e. it worked + + # Check errors trigger + with pytest.raises(Exception) as excinfo: + # Parameters + alpha, beta, gamma, x = 1.91, 0.39, 0.1, 3 + + # Catalogue + m_min, m_max, size = 1e+10, 1e+16, 1000 + parent_halo = colossus_mass_sampler(redshift=0.1, model='sheth99', + mdef='fof', m_min=m_min, m_max=m_max, + cosmology=WMAP9, sigma8=0.8, ns=1., + size=size, resolution=1000) + z_halo = np.linspace(0.01, 0.1, 500) + gen_sub_cat(parent_halo, z_halo, alpha, beta, gamma, x) + assert str(excinfo.value) == 'Catalogue of halos and redshifts must be the same length' + + with pytest.raises(Exception) as excinfo: + # Parameters + alpha, beta, gamma, x = 1.91, 0.39, 0.1, 3 + + # Catalogue + parent_halo = [-10, 1, 0] + z_halo = np.linspace(0.01, 0.1, 3) + gen_sub_cat(parent_halo, z_halo, alpha, beta, gamma, x) + assert str(excinfo.value) == 'Masses in catalogue should be positive and non-zero' + + with pytest.raises(Exception) as excinfo: + # Parameters + alpha, beta, gamma, x = 1.91, 0.39, 0.1, 3 + + # Catalogue + parent_halo = 10**(np.array([10, 12, 11])) + z_halo = [-0.05, 0.1, 0.01] + gen_sub_cat(parent_halo, z_halo, alpha, beta, gamma, x) + assert str(excinfo.value) == 'Redshifts in catalogue should be positive' + + with pytest.raises(Exception) as excinfo: + # Parameters + alpha, beta, gamma, x = -1.91, 0.39, 0.1, 3 + + # Catalogue + parent_halo = 10**(np.array([10, 12, 11])) + z_halo = [0.05, 0.1, 0.01] + gen_sub_cat(parent_halo, z_halo, alpha, beta, gamma, x) + assert str(excinfo.value) == 'Subhalo mass function defined alpha > 0, set_alpha = -alpha' + + with pytest.raises(Exception) as excinfo: + # Parameters + alpha, beta, gamma, x = 2, 0.39, 0.1, 3 + + # Catalogue + parent_halo = 10**(np.array([10, 12, 11])) + z_halo = [0.05, 0.1, 0.01] + gen_sub_cat(parent_halo, z_halo, alpha, beta, gamma, x) + assert str(excinfo.value) == 'Subhalo alpha must be less than 2' + + with pytest.raises(Exception) as excinfo: + # Parameters + alpha, beta, gamma, x = 1.91, 0.39, 0.1, 0.5 + + # Catalogue + parent_halo = 10**(np.array([10, 12, 11])) + z_halo = [0.05, 0.1, 0.01] + gen_sub_cat(parent_halo, z_halo, alpha, beta, gamma, x) + assert str(excinfo.value) == 'Subhalo x cannot be less than 1' + + with pytest.raises(Exception) as excinfo: + # Parameters + alpha, beta, gamma, x = 1.91, -0.39, 0.1, 3 + + # Catalogue + parent_halo = 10**(np.array([10, 12, 11])) + z_halo = [0.05, 0.1, 0.01] + gen_sub_cat(parent_halo, z_halo, alpha, beta, gamma, x) + assert str(excinfo.value) == 'Subhalo beta must be between 0 and 1' + + with pytest.raises(Exception) as excinfo: + # Parameters + alpha, beta, gamma, x = 1.91, 1.39, -0.1, 3 + + # Catalogue + parent_halo = 10**(np.array([10, 12, 11])) + z_halo = [0.05, 0.1, 0.01] + gen_sub_cat(parent_halo, z_halo, alpha, beta, gamma, x) + assert str(excinfo.value) == 'Subhalo beta must be between 0 and 1' + + with pytest.raises(Exception) as excinfo: + # Parameters + alpha, beta, gamma, x = 1.91, 0.39, -0.1, 3 + + # Catalogue + parent_halo = 10**(np.array([10, 12, 11])) + z_halo = [0.05, 0.1, 0.01] + gen_sub_cat(parent_halo, z_halo, alpha, beta, gamma, x) + assert str(excinfo.value) == 'Subhalo gamma must be between 0 and 1' + + with pytest.raises(Exception) as excinfo: + # Parameters + alpha, beta, gamma, x = 1.91, 1.39, 1.1, 3 + + # Catalogue + parent_halo = 10**(np.array([10, 12, 11])) + z_halo = [0.05, 0.1, 0.01] + gen_sub_cat(parent_halo, z_halo, alpha, beta, gamma, x) + assert str(excinfo.value) == 'Subhalo gamma must be between 0 and 1' + + +@pytest.mark.skipif(not HAS_COLOSSUS, reason='test requires colossus') +@pytest.mark.flaky +def test_assignment(): + # Run a test on the assignment function to check catalogues are matched correctly + + # Generate the catalogues + from astropy.cosmology import WMAP9 + from astropy.cosmology import FlatLambdaCDM + from astropy import units as u + from skypy.halos.mass import colossus_mass_sampler + from skypy.halos.sham import assignment + from skypy.halos.sham import gen_sub_cat + from skypy.halos.sham import quenching_funct + from skypy.galaxies._schechter import schechter_smf + + # Catalogues + cosmo = FlatLambdaCDM(H0=70, Om0=0.3, name='FlatLambdaCDM') + m_min, m_max, size = 1e+10, 1e+16, 1000 + halo_cat = colossus_mass_sampler(redshift=0.1, model='sheth99', + mdef='fof', m_min=m_min, m_max=m_max, + cosmology=WMAP9, sigma8=0.8, ns=1., + size=size, resolution=1000) + + h_z = np.linspace(0.01, 0.1, len(halo_cat)) + z_range = [min(h_z), max(h_z)] + skyarea = 600*u.deg**2 + + alpha, beta, gamma, x = 1.91, 0.39, 0.1, 3 + ID_halo, sub_masses, ID_sub, z_sub = gen_sub_cat(halo_cat, h_z, alpha, beta, gamma, x) + + m_star1, phi_star1, alpha1 = 10**(10.58), 10**(-2.77), -0.33 + m_star2, phi_star2, alpha2 = 10**(10.64), 10**(-4.24), -1.54 + m_star3, phi_star3, alpha3 = 10**(10.65), 10**(-2.98), -1.48 + m_star4, phi_star4, alpha4 = 10**(10.55), 10**(-3.96), -1.53 + + m_min1 = 10**(6.7) + m_min2 = 10**(6.6) + m_min3 = 10**(7.0) + m_min4 = 10**(7.0) + + m_max1 = 10**(14) + m_max2 = 10**(13) + + rc_cat = schechter_smf(z_range, m_star1, phi_star1, alpha1, m_min1, m_max1, + skyarea, cosmo, noise=False)[1] + rs_cat = schechter_smf(z_range, m_star2, phi_star2, alpha2, m_min2, m_max2, + skyarea, cosmo, noise=False)[1] + bc_cat = schechter_smf(z_range, m_star3, phi_star3, alpha3, m_min3, m_max1, + skyarea, cosmo, noise=False)[1] + bs_cat = schechter_smf(z_range, m_star4, phi_star4, alpha4, m_min4, m_max2, + skyarea, cosmo, noise=False)[1] + + # Quench the halos + h_quench = quenching_funct(halo_cat, 10**(12), 0.4) + s_quench = quenching_funct(sub_masses, 10**(12), 0.4, 0.4) + + # Order the arrays + halo_subhalo = np.concatenate((halo_cat, sub_masses), axis=0) # Halos + ID_list = np.concatenate((ID_halo, ID_sub), axis=0) # IDs + z_list = np.concatenate((h_z, z_sub), axis=0) # Redshifts + q_list = np.concatenate((h_quench, s_quench), axis=0) # Quenching + + stack = np.stack((halo_subhalo, ID_list, z_list, q_list), axis=1) + order1 = stack[np.argsort(stack[:, 0])] + hs_order = np.flip(order1[:, 0]) + id_hs = np.flip(order1[:, 1]) + z_hs = np.flip(order1[:, 2]) + qu_hs = np.flip(order1[:, 3]) + rc_order = np.flip(np.sort(rc_cat)) # Galaxies + rs_order = np.flip(np.sort(rs_cat)) + bc_order = np.flip(np.sort(bc_cat)) + bs_order = np.flip(np.sort(bs_cat)) + + # Function + hs_fin, gal_fin, id_fin, z_fin, gal_type = assignment(hs_order, rc_order, rs_order, + bc_order, bs_order, qu_hs, id_hs, z_hs) + + # Check arrays shape/type + assert hs_fin.shape == gal_fin.shape == id_fin.shape == z_fin.shape == gal_type.shape + assert hs_fin.ndim == 1 + + # Check arrays contain expected values + assert ((gal_type == 'red central') | (gal_type == 'red satellite') | + (gal_type == 'blue central') | (gal_type == 'blue satellite')).all() + assert len(np.where(gal_type == 'red central')[0]) > 1 # Check some of each were assigned + assert len(np.where(gal_type == 'red satellite')[0]) > 1 + assert len(np.where(gal_type == 'blue central')[0]) > 1 + assert len(np.where(gal_type == 'blue satellite')[0]) > 1 + + assert min(z_fin) >= min(z_range) # z within range + assert max(z_fin) <= max(z_range) + + i_type = [] + for ii in id_fin: + if abs(ii - int(ii)) > 0: + i_type.append(False) + else: + i_type.append(True) + + assert np.all(i_type) # Check ids are all integers + assert len(np.where(id_fin > 0)[0]) > 1 # Check there are some subhalos + assert len(np.where(id_fin < 0)[0]) > 1 # Check there are some parents + + assert max(gal_fin) <= 10**(15) # Check galaxies in correct range + assert min(gal_fin) >= 10**(4) + assert max(hs_fin) <= 10**(17) # Check halos in correct range + assert min(hs_fin) >= 10**(7) + + # Check galaxies are assigned to correct halos + check_parent_t = gal_type[np.where(id_fin < 0)] # Finished arrays + check_parent_m = gal_fin[np.where(id_fin < 0)] + check_subhal_t = gal_type[np.where(id_fin > 0)] + check_subhal_m = gal_fin[np.where(id_fin > 0)] + check_blu_t = gal_type[np.where(qu_hs == 0)] + check_blu_m = gal_fin[np.where(qu_hs == 0)] + check_red_t = gal_type[np.where(qu_hs == 1)] + check_red_m = gal_fin[np.where(qu_hs == 1)] + + # Check correct tags are applied + assert ((check_parent_t == 'red central') | (check_parent_t == 'blue central')).all() + assert ((check_subhal_t == 'red satellite') | (check_subhal_t == 'blue satellite')).all() + assert ((check_blu_t == 'blue central') | (check_blu_t == 'blue satellite')).all() + assert ((check_red_t == 'red central') | (check_red_t == 'red satellite')).all() + + cen_append = np.append(rc_order, bc_order) # Catalogue + sub_append = np.append(rs_order, bs_order) + blu_append = np.append(bc_order, bs_order) + red_append = np.append(rc_order, rs_order) + + assert np.all(np.isin(check_parent_m, cen_append)) # Check masses are from correct lists + assert np.all(np.isin(check_subhal_m, sub_append)) + assert np.all(np.isin(check_blu_m, blu_append)) + assert np.all(np.isin(check_red_m, red_append)) + + # Check errors trigger + with pytest.raises(Exception) as excinfo: + hs_test = -1*hs_order + hs_fin, gal_fin, id_fin, z_fin, gal_type = assignment(hs_test, rc_order, rs_order, + bc_order, bs_order, qu_hs, id_hs, + z_hs) + assert str(excinfo.value) == 'Halo masses must be positive and non-zero' + + with pytest.raises(Exception) as excinfo: + rs_test = -1*rs_order + hs_fin, gal_fin, id_fin, z_fin, gal_type = assignment(hs_order, rc_order, rs_test, + bc_order, bs_order, qu_hs, id_hs, + z_hs) + assert str(excinfo.value) == 'Galaxy masses must be positive and non-zero' + + with pytest.raises(Exception) as excinfo: + hs_test = np.flip(hs_order) + hs_fin, gal_fin, id_fin, z_fin, gal_type = assignment(hs_test, rc_order, rs_order, + bc_order, bs_order, qu_hs, id_hs, + z_hs) + assert str(excinfo.value) == 'Halo masses were in the wrong order and have been corrected' + + with pytest.raises(Exception) as excinfo: + hs_test = np.random.randint(min(hs_order), max(hs_order), len(hs_order)) + hs_fin, gal_fin, id_fin, z_fin, gal_type = assignment(hs_test, rc_order, rs_order, + bc_order, bs_order, qu_hs, id_hs, + z_hs) + assert str(excinfo.value) == 'Halos are not in a sorted order' + + with pytest.raises(Exception) as excinfo: + rc_test = np.flip(rc_order) + hs_fin, gal_fin, id_fin, z_fin, gal_type = assignment(hs_order, rc_test, rs_order, + bc_order, bs_order, qu_hs, id_hs, + z_hs) + assert str(excinfo.value) == 'Red central galaxies were in wrong order now correct' + + with pytest.raises(Exception) as excinfo: + rc_test = np.random.randint(min(rc_order), max(rc_order), len(rc_order)) + hs_fin, gal_fin, id_fin, z_fin, gal_type = assignment(hs_order, rc_test, rs_order, + bc_order, bs_order, qu_hs, id_hs, + z_hs) + assert str(excinfo.value) == 'Red central galaxies are not in a sorted order' + + with pytest.raises(Exception) as excinfo: + rs_test = np.flip(rs_order) + hs_fin, gal_fin, id_fin, z_fin, gal_type = assignment(hs_order, rc_order, rs_test, + bc_order, bs_order, qu_hs, id_hs, + z_hs) + assert str(excinfo.value) == 'Red satellite galaxies were in wrong order now correct' + + with pytest.raises(Exception) as excinfo: + rs_test = np.random.randint(min(rs_order), max(rs_order), len(rs_order)) + hs_fin, gal_fin, id_fin, z_fin, gal_type = assignment(hs_order, rc_order, rs_test, + bc_order, bs_order, qu_hs, id_hs, + z_hs) + assert str(excinfo.value) == 'Red satellite galaxies are not in a sorted order' + + with pytest.raises(Exception) as excinfo: + bc_test = np.flip(bc_order) + hs_fin, gal_fin, id_fin, z_fin, gal_type = assignment(hs_order, rc_order, rs_order, + bc_test, bs_order, qu_hs, id_hs, + z_hs) + assert str(excinfo.value) == 'Blue central galaxies were in wrong order and now correct' + + with pytest.raises(Exception) as excinfo: + bc_test = np.random.randint(min(bc_order), max(bc_order), len(bc_order)) + hs_fin, gal_fin, id_fin, z_fin, gal_type = assignment(hs_order, rc_order, rs_order, + bc_test, bs_order, qu_hs, id_hs, + z_hs) + assert str(excinfo.value) == 'Blue central galaxies are not in a sorted order' + + with pytest.raises(Exception) as excinfo: + bs_test = np.flip(bs_order) + hs_fin, gal_fin, id_fin, z_fin, gal_type = assignment(hs_order, rc_order, rs_order, + bc_order, bs_test, qu_hs, id_hs, + z_hs) + assert str(excinfo.value) == 'Blue satellite galaxies were in wrong order and now correct' + + with pytest.raises(Exception) as excinfo: + bs_test = np.random.randint(min(bs_order), max(bs_order), len(bs_order)) + hs_fin, gal_fin, id_fin, z_fin, gal_type = assignment(hs_order, rc_order, rs_order, + bc_order, bs_test, qu_hs, id_hs, + z_hs) + assert str(excinfo.value) == 'Blue satellite galaxies are not in a sorted order' + + with pytest.raises(Exception) as excinfo: + hs_test = np.geomspace(max(hs_order), min(hs_order), 10) + hs_fin, gal_fin, id_fin, z_fin, gal_type = assignment(hs_test, rc_order, rs_order, + bc_order, bs_order, qu_hs, id_hs, + z_hs) + assert str(excinfo.value) == 'All arrays pertaining to halos must be the same shape' + + with pytest.raises(Exception) as excinfo: + id_test = np.geomspace(min(id_hs), max(id_hs), 10) + hs_fin, gal_fin, id_fin, z_fin, gal_type = assignment(hs_test, rc_order, rs_order, + bc_order, bs_order, qu_hs, id_test, + z_hs) + assert str(excinfo.value) == 'All arrays pertaining to halos must be the same shape' + + with pytest.raises(Exception) as excinfo: + z_test = np.geomspace(min(z_hs), max(z_hs), 10) + hs_fin, gal_fin, id_fin, z_fin, gal_type = assignment(hs_order, rc_order, rs_order, + bc_order, bs_order, qu_hs, id_hs, + z_test) + assert str(excinfo.value) == 'All arrays pertaining to halos must be the same shape' + + +@pytest.mark.skipif(not HAS_COLOSSUS, reason='test requires colossus') +@pytest.mark.flaky +def test_scatter_proxy(): + # Run a test on the scattering function to check outputs are correct + + from skypy.halos.sham import run_file + from skypy.halos.sham import scatter_proxy + + # Generate the galaxy catalogue and scatter + gal_cat = run_file('test_gal.yaml', 'galaxy', 'sm') + new_gal = scatter_proxy(gal_cat) + + # Check new galaxies are not in direct mass order + assert np.any(np.diff(new_gal) > 0) + + # Check input galaxies are just rearranged + assert np.all(np.isin(gal_cat, new_gal)) + + # Check that the general trend is reducing mass + # By difference in mass means + sets = 500 + ii = 0 + mean = [] + while ii < len(new_gal): + gal_set = new_gal[ii:ii+sets] + mean.append(np.mean(gal_set)) + ii += sets + + assert np.all(np.diff(mean) < 0) + + # By start and end values + assert new_gal[0] > new_gal[-1] + + with pytest.raises(Exception) as excinfo: + gal_cat = np.array([-1, 2, 3, 0, -10, -20]) + scatter_proxy(gal_cat) + assert str(excinfo.value) == 'Galaxy masses must be positive and non-zero' + + +@pytest.mark.skipif(not HAS_COLOSSUS, reason='test requires colossus') +@pytest.mark.flaky +def test_sham_plots(): + # Run a test on the plotting function to check outputs are correctly labelled + + # Generate the catalogues + from astropy.cosmology import WMAP9 + from astropy.cosmology import FlatLambdaCDM + from astropy import units as u + from skypy.halos.mass import colossus_mass_sampler + from skypy.halos.sham import assignment + from skypy.halos.sham import gen_sub_cat + from skypy.halos.sham import quenching_funct + from skypy.halos.sham import sham_plots + from skypy.galaxies._schechter import schechter_smf + + # Catalogues + cosmo = FlatLambdaCDM(H0=70, Om0=0.3, name='FlatLambdaCDM') + m_min, m_max, size = 1e+10, 1e+16, 1000 + halo_cat = colossus_mass_sampler(redshift=0.1, model='sheth99', + mdef='fof', m_min=m_min, m_max=m_max, + cosmology=WMAP9, sigma8=0.8, ns=1., + size=size, resolution=1000) + h_z = np.linspace(0.01, 0.1, len(halo_cat)) + z_range = [min(h_z), max(h_z)] + skyarea = 600*u.deg**2 + + alpha, beta, gamma, x = 1.91, 0.39, 0.1, 3 + ID_halo, sub_masses, ID_sub, z_sub = gen_sub_cat(halo_cat, h_z, alpha, beta, gamma, x) + + m_star1, phi_star1, alpha1 = 10**(10.58), 10**(-2.77), -0.33 + m_star2, phi_star2, alpha2 = 10**(10.64), 10**(-4.24), -1.54 + m_star3, phi_star3, alpha3 = 10**(10.65), 10**(-2.98), -1.48 + m_star4, phi_star4, alpha4 = 10**(10.55), 10**(-3.96), -1.53 + + m_min1 = 10**(6.7) + m_min2 = 10**(6.6) + m_min3 = 10**(7.0) + m_min4 = 10**(7.0) + + m_max1 = 10**(14) + m_max2 = 10**(13) + + rc_cat = schechter_smf(z_range, m_star1, phi_star1, alpha1, m_min1, m_max1, + skyarea, cosmo, noise=False)[1] + rs_cat = schechter_smf(z_range, m_star2, phi_star2, alpha2, m_min2, m_max2, + skyarea, cosmo, noise=False)[1] + bc_cat = schechter_smf(z_range, m_star3, phi_star3, alpha3, m_min3, m_max1, + skyarea, cosmo, noise=False)[1] + bs_cat = schechter_smf(z_range, m_star4, phi_star4, alpha4, m_min4, m_max2, + skyarea, cosmo, noise=False)[1] + + # Quench the halos + h_quench = quenching_funct(halo_cat, 10**(12), 0.4) + s_quench = quenching_funct(sub_masses, 10**(12), 0.4, 0.4) + + # Order the arrays + halo_subhalo = np.concatenate((halo_cat, sub_masses), axis=0) # Halos + ID_list = np.concatenate((ID_halo, ID_sub), axis=0) # IDs + z_list = np.concatenate((h_z, z_sub), axis=0) # Redshifts + q_list = np.concatenate((h_quench, s_quench), axis=0) # Quenching + + stack = np.stack((halo_subhalo, ID_list, z_list, q_list), axis=1) + order1 = stack[np.argsort(stack[:, 0])] + hs_order = np.flip(order1[:, 0]) + id_hs = np.flip(order1[:, 1]) + z_hs = np.flip(order1[:, 2]) + qu_hs = np.flip(order1[:, 3]) + rc_order = np.flip(np.sort(rc_cat)) # Galaxies + rs_order = np.flip(np.sort(rs_cat)) + bc_order = np.flip(np.sort(bc_cat)) + bs_order = np.flip(np.sort(bs_cat)) + + # Assign the galaxies + hs_fin, gal_fin, id_fin, z_fin, gal_type_fin = assignment(hs_order, rc_order, rs_order, + bc_order, bs_order, qu_hs, id_hs, + z_hs) + + # Function + sham_rc, sham_rs, sham_bc, sham_bs, sham_cen, sham_sub = sham_plots(hs_fin, gal_fin, + gal_type_fin) + + # Arrays correct type and size + assert type(sham_rc) is type(sham_rs) is type(sham_bc) is type(sham_bs) is np.ndarray + assert sham_rc.ndim == sham_rs.ndim == sham_bc.ndim == sham_bs.ndim == 2 + assert sham_cen.ndim == sham_sub.ndim == 2 + assert sham_rc.size # Check lists are not empty + assert sham_rs.size + assert sham_bc.size + assert sham_bs.size + assert sham_cen.size + assert sham_sub.size + + assert z_fin.size + + # Check number of centrals and satellites correct + no_cen = len(np.where(id_fin < 0)[0]) + no_sub = len(np.where(id_fin > 0)[0]) + assert len(sham_cen[:, 0]) == no_cen + assert len(sham_sub[:, 0]) == no_sub + + # Check values correct + assert np.all(np.isin(sham_rc[:, 0], hs_fin)) + assert np.all(np.isin(sham_rs[:, 0], hs_fin)) + assert np.all(np.isin(sham_bc[:, 0], hs_fin)) + assert np.all(np.isin(sham_bs[:, 0], hs_fin)) + + assert np.all(np.isin(sham_rc[:, 1], gal_fin)) + assert np.all(np.isin(sham_rs[:, 1], gal_fin)) + assert np.all(np.isin(sham_bc[:, 1], gal_fin)) + assert np.all(np.isin(sham_bs[:, 1], gal_fin)) + + assert np.all(np.isin(sham_cen[:, 0], hs_fin)) + assert np.all(np.isin(sham_sub[:, 0], hs_fin)) + assert np.all(np.isin(sham_cen[:, 1], gal_fin)) + assert np.all(np.isin(sham_sub[:, 1], gal_fin)) + + # Check errors trigger + with pytest.raises(Exception) as excinfo: + hs_test = np.linspace(min(hs_order), max(hs_order), 10) + sham_plots(hs_test, gal_fin, gal_type_fin) + assert str(excinfo.value) == 'All arrays must be the same shape' + + with pytest.raises(Exception) as excinfo: + gal_type_test = np.random.randint(1, 5, 10) + sham_plots(hs_fin, gal_fin, gal_type_test) + assert str(excinfo.value) == 'All arrays must be the same shape' + + with pytest.raises(Exception) as excinfo: + hs_test = -1*(hs_fin) + sham_plots(hs_test, gal_fin, gal_type_fin) + assert str(excinfo.value) == 'Halo masses must be positive and non-zero' + + with pytest.raises(Exception) as excinfo: + gal_test = -1*(gal_fin) + sham_plots(hs_fin, gal_test, gal_type_fin) + assert str(excinfo.value) == 'Galaxy masses must be positive and non-zero' + + +@pytest.mark.skipif(not HAS_COLOSSUS, reason='test requires colossus') +@pytest.mark.flaky +def test_run_sham(): + # Run a test on the SHAM function to check outputs are in correct format + + from astropy.cosmology import FlatLambdaCDM + from skypy.halos.sham import run_sham + + # Parameters + h_file = 'test_halo.yaml' + gal_param = np.array([[10**(10.58), 10**(-2.77), -0.33], [10**(10.64), 10**(-4.24), -1.54], + [10**(10.65), 10**(-2.98), -1.48], [10**(10.55), 10**(-3.96), -1.53]]) + cosmology = FlatLambdaCDM(H0=70, Om0=0.3, name='FlatLambdaCDM') + z_range = [0.01, 0.1] + skyarea = 60 + qu = np.array([10**(12.16), 0.45, 0.5]) + + # Function + sham = run_sham(h_file, gal_param, cosmology, z_range, skyarea, qu, + sub_param=[1.91, 0.39, 0.1, 3], gal_max_h=10**(14), gal_max_s=10**(13), + print_out=False, run_anyway=True) + + # Check arrays + assert type(sham['Halo mass']) is np.ndarray + assert type(sham['Galaxy mass']) is np.ndarray + assert type(sham['Galaxy type']) is np.ndarray + assert type(sham['ID value']) is np.ndarray + assert type(sham['Redshift']) is np.ndarray + + assert len(sham['Halo mass']) == len(sham['Galaxy mass']) + assert len(sham['Galaxy mass']) == len(sham['Galaxy type']) + assert len(sham['Galaxy type']) == len(sham['ID value']) + assert len(sham['ID value']) == len(sham['Redshift']) + + # Check errors trigger + with pytest.raises(Exception) as excinfo: + h_f = 3 + run_sham(h_f, gal_param, cosmology, z_range, skyarea, qu, + sub_param=[1.91, 0.39, 0.1, 3], gal_max_h=10**(14), gal_max_s=10**(13), + print_out=False, run_anyway=True) + assert str(excinfo.value) == 'Halo YAML file must be provided as a string' + + with pytest.raises(Exception) as excinfo: + gal_p = [[1, 2, 3], [4, 5, 6]] + run_sham(h_file, gal_p, cosmology, z_range, skyarea, qu, + sub_param=[1.91, 0.39, 0.1, 3], gal_max_h=10**(14), gal_max_s=10**(13), + print_out=False, run_anyway=True) + assert str(excinfo.value) == 'The wrong number of galaxies are in galaxy parameters' + + with pytest.raises(Exception) as excinfo: + gal_p = [[1, 2], [4, 5], [7, 8], [3, 4]] + run_sham(h_file, gal_p, cosmology, z_range, skyarea, qu, + sub_param=[1.91, 0.39, 0.1, 3], gal_max_h=10**(14), gal_max_s=10**(13), + print_out=False, run_anyway=True) + assert str(excinfo.value) == 'The wrong number of galaxy parameters have been provided' + + with pytest.raises(Exception) as excinfo: + gal_p = [[1, 2], [4, 5], [7, 8]] + run_sham(h_file, gal_p, cosmology, z_range, skyarea, qu, + sub_param=[1.91, 0.39, 0.1, 3], gal_max_h=10**(14), gal_max_s=10**(13), + print_out=False, run_anyway=True) + assert str(excinfo.value) == 'Supplied galaxy parameters are not the correct shape' + + with pytest.raises(Exception) as excinfo: + qu_t = [1] + run_sham(h_file, gal_param, cosmology, z_range, skyarea, qu_t, + sub_param=[1.91, 0.39, 0.1, 3], gal_max_h=10**(14), gal_max_s=10**(13), + print_out=False, run_anyway=True) + assert str(excinfo.value) == 'Provided incorrect number of quenching parameters' + + with pytest.raises(Exception) as excinfo: + z_t = [1] + run_sham(h_file, gal_param, cosmology, z_t, skyarea, qu, + sub_param=[1.91, 0.39, 0.1, 3], gal_max_h=10**(14), gal_max_s=10**(13), + print_out=False, run_anyway=True) + assert str(excinfo.value) == 'The wrong number of redshifts were given' + + with pytest.raises(Exception) as excinfo: + z_t = [-0.01, 0.1] + run_sham(h_file, gal_param, cosmology, z_t, skyarea, qu, + sub_param=[1.91, 0.39, 0.1, 3], gal_max_h=10**(14), gal_max_s=10**(13), + print_out=False, run_anyway=True) + assert str(excinfo.value) == 'Redshift cannot be negative' + + with pytest.raises(Exception) as excinfo: + z_t = [0.1, 0.01] + run_sham(h_file, gal_param, cosmology, z_t, skyarea, qu, + sub_param=[1.91, 0.39, 0.1, 3], gal_max_h=10**(14), gal_max_s=10**(13), + print_out=False, run_anyway=True) + assert str(excinfo.value) == 'The second redshift should be more than the first' + + with pytest.raises(Exception) as excinfo: + rc_p = [-10**(10.58), 10**(-2.77), -0.33] + rs_p = [10**(10.64), 10**(-4.24), -1.54] + bc_p = [0, 10**(-2.98), -1.48] + bs_p = [10**(10.55), 10**(-3.96), -1.53] + gal_t = [rc_p, rs_p, bc_p, bs_p] + run_sham(h_file, gal_t, cosmology, z_range, skyarea, qu, + sub_param=[1.91, 0.39, 0.1, 3], gal_max_h=10**(14), gal_max_s=10**(13), + print_out=False, run_anyway=True) + assert str(excinfo.value) == 'M* values must be positive and non-zero' + + with pytest.raises(Exception) as excinfo: + rc_p = [10**(10.58), -10**(-2.77), -0.33] + rs_p = [10**(10.64), 10**(-4.24), -1.54] + bc_p = [10**(10.65), 0, -1.48] + bs_p = [10**(10.55), 10**(-3.96), -1.53] + gal_t = [rc_p, rs_p, bc_p, bs_p] + run_sham(h_file, gal_t, cosmology, z_range, skyarea, qu, + sub_param=[1.91, 0.39, 0.1, 3], gal_max_h=10**(14), gal_max_s=10**(13), + print_out=False, run_anyway=True) + assert str(excinfo.value) == 'phi* values must be positive and non-zero' + + with pytest.raises(Exception) as excinfo: + rc_p = [10**(10.58), 10**(-2.77), 0.33] + rs_p = [10**(10.64), 10**(-4.24), -1.54] + bc_p = [10**(10.65), 10**(-2.98), -1.48] + bs_p = [10**(10.55), 10**(-3.96), -1.53] + gal_t = [rc_p, rs_p, bc_p, bs_p] + run_sham(h_file, gal_t, cosmology, z_range, skyarea, qu, + sub_param=[1.91, 0.39, 0.1, 3], gal_max_h=10**(14), gal_max_s=10**(13), + print_out=False, run_anyway=True) + assert str(excinfo.value) == 'Galaxy Schechter function alphas must be < 0' + + with pytest.raises(Exception) as excinfo: + sky_t = -10 + run_sham(h_file, gal_param, cosmology, z_range, sky_t, qu, + sub_param=[1.91, 0.39, 0.1, 3], gal_max_h=10**(14), gal_max_s=10**(13), + print_out=False, run_anyway=True) + assert str(excinfo.value) == 'The skyarea must be a positive non-zero number' + + with pytest.raises(Exception) as excinfo: + cosmo_t = FlatLambdaCDM(H0=70, Om0=0.3) + run_sham(h_file, gal_param, cosmo_t, z_range, skyarea, qu, + sub_param=[1.91, 0.39, 0.1, 3], gal_max_h=10**(14), gal_max_s=10**(13), + print_out=False, run_anyway=True) + assert str(excinfo.value) == 'Cosmology object must have an astropy cosmology name' diff --git a/skypy/power_spectrum/tests/test_growth.py b/skypy/power_spectrum/tests/test_growth.py index efc65fc5f..53a4c6e59 100644 --- a/skypy/power_spectrum/tests/test_growth.py +++ b/skypy/power_spectrum/tests/test_growth.py @@ -57,15 +57,15 @@ def test_growth(): Dzprime = growth_function_derivative(redshift, cosmology_flat) # Test growth factor - assert redshift.shape == fz.shape,\ + assert redshift.shape == fz.shape, \ "Length of redshift array and growth rate array do not match" - assert isclose(fz[0], 1.0),\ + assert isclose(fz[0], 1.0), \ "Growth factor at redshift 0 is not close to 1.0" # Test growth function - assert redshift.shape == Dz.shape,\ + assert redshift.shape == Dz.shape, \ "Length of redshift array and growth function array do not match" - assert allclose(Dz, 1. / (1. + redshift)),\ + assert allclose(Dz, 1. / (1. + redshift)), \ "Growth function is not close to the scale factor" # make sure that growth_function with scalar returns scalar @@ -74,9 +74,9 @@ def test_growth(): assert Dz2 == Dz[2], 'growth function with scalar produced inconsistent result' # Test growth function derivative - assert redshift.shape == Dzprime.shape,\ + assert redshift.shape == Dzprime.shape, \ "Length of redshift array and growth function array do not match" - assert isclose(Dzprime[0], -1.0),\ + assert isclose(Dzprime[0], -1.0), \ "Derivative of growth function at redshift 0 is not close to -1.0" # Test against precomputed values using Planck15 cosmology