Commit f7c80363 authored by Médéric Boquien's avatar Médéric Boquien

Merge branch 'release/v2018.0'

parents 0769fddf f7623fbf
pcigale authors list pcigale authors list
==================== ====================
This document lists alphabetically the various authors who wrote the pcigale This document lists alphabetically the core team who wrote the pcigale code
code with their current email address and affiliation. with their current email address and affiliation.
* Médéric Boquien <mederic.boquien@uantof.cl>, * Médéric Boquien <mederic.boquien@uantof.cl>,
Universidad de Antofagasta, Chile Universidad de Antofagasta, Chile
...@@ -10,4 +10,8 @@ code with their current email address and affiliation. ...@@ -10,4 +10,8 @@ code with their current email address and affiliation.
Laboratoire d'Astrophysique de Marseille, France Laboratoire d'Astrophysique de Marseille, France
* Laure Ciesla <ciesla@lam.fr>, * Laure Ciesla <ciesla@lam.fr>,
Laboratoire d'Astrophysique de Marseille, France Laboratoire d'Astrophysique de Marseille, France
* David Corre <david.corre@lam.fr>
Laboratoire d'Astrophysique de Marseille, France
* Yannick Roehlly <yannick@iaora.eu> * Yannick Roehlly <yannick@iaora.eu>
* Héctor Salas Olave <hector.salas.o@gmail.com>
Universidad de Antofagasta, Chile
# Change Log # Change Log
## 2018.0 (2018-11-06)
### Added
- It is now possible to optionally indicate the distance in Mpc in the input file. If present it will be used in lieu of the distance computed from the redshift. This is especially useful in the nearby universe where the redshift is a very poor indicator of the actual distance. (Médéric Boquien)
- It is now possible to fit any physical property indicated by the code (e.g. equivalent width, dust luminosity, etc.). For this the physical property needs to be given in the input file and the properties to be fitted must be given in the properties filed in pcigale.ini. (Héctor Salas & Médéric Boquien)
- It is now possible to fit emission lines. For this the line has to be indicated in the same way as any other band both in the input flux file (in units of W/m²) and in the list of bands in `pcigale.ini`. Lines are prefixed with `line.` followed by the name of the line, for instance `line.H-alpha` for Hɑ. The following lines are supported at the moment: `Ly-alpha`, `CII-133.5`, `SiIV-139.7`, `CIV-154.9`, `HeII-164.0`, `OIII-166.5`, `CIII-190.9`, `CII-232.6`, `MgII-279.8`, `OII-372.7`, `H-10`, `H-9`, `NeIII-386.9` `HeI-388.9`, `H-epsilon`, `SII-407.0`, `H-delta`, `H-gamma`, `H-beta`, `OIII-495.9`, `OIII-500.7`, `OI-630.0`, `NII-654.8`, `H-alpha`, `NII-658.4`, `SII-671.6`, `SII-673.1`. (Médéric Boquien)
- When emission lines are not corrected for absorption lines (e.g., in the case of very low resolution spectroscopy) the previous method, which computes the theoretical line fluxes, is not optimal. Rather we offer the possibility to measure the fluxes through special filters that are used like regular filters. The idea is to define filters with a positive part on the line, a negative part on the continuum, and a zero-valued integral. In such case the integration over the spectrum directly gives the flux of the integral. So this works at all redshifts, the filter is automatically redshifted at runtime. (Médéric Boquien & David Corre)
- Two new dust attenuation modules have been added: `dustatt\_modified\_CF00` and `dustatt\_modified\_starburst`. The former implements a modified 2-component Charlot & Fall (2000) model whereas the latter implements a modified starburst law with the continuum attenuated with a Calzetti (2000) curve and the lines extincted with a Milky Way or a Magellanic Cloud law. The previous models `dustatt\_powerlaw`, `dustatt\_2powerlaws`, and `dustatt\_calzleit` are still available but are deprecated. (Médéric Boquien & David Corre)
- In addition to the physical properties, the fluxes are now also estimated through a Bayesian analysis. (Médéric Boquien)
- The module `sfhdelayedbq` has been added. It implements a delayed SFH with a burst/quench. It is fully described in Ciesla et al. (2017).
### Changed
- The `sfhdelayed` module has been extended to optionally include an exponential burst to model the latest episode of star formation. (Médéric Boquien & Barbara Lo Faro)
- On Linux platforms the method to start the parallel processes has been changed from "spawn" to "fork". This allows for a much faster startup. On other platforms is remains unchanged as Windows does not support "fork" and MacOS is bugged when using "fork", resulting in a free of cigale. (Médéric Boquien)
- The list of modules has been made more explicit in the `pcigale.ini` file. (Médéric Boquien)
### Fixed
- The histogram bin width was not computed optimally when some models were invalid. (David Corre & Médéric Boquien)
- Missing import in the `m2005` module. (Médéric Boquien, reported by Dominika Wylezalek)
- The plot of the PDF could not be generated for physical properties estimated in log (Médéric Boquien)
- We do not attempt anymore to estimate the physical properties of galaxies with insanely large χ² that lead to an underflow in the computation of the likelihood. (Médéric Boquien)
- The best fit is now plotted at the exact redshift rather than at the rounded redshift. (Médéric Boquien)
- It is now possible to plot the best fit obtained in redshifting mode. (Médéric Boquien)
### Optimised
- The estimation of the physical properties is made a bit faster when all the models are valid. (Médéric Boquien)
- The access to the module cache has been made faster and the model cache has been made much simpler, avoiding plenty of complex computations. This results in a speedup of at least ~6% in the computation of the models. The speedup can be higher when using few photometric bands. At the same time it considerably reduces the number of page faults seen in some rare circumstances. (Médéric Boquien)
- The models counter was a bottleneck when using many cores as updating it could stall other parallel processes. Now the internal counter is updated much less frequently. The speedup goes from between negligible (few cores) up to a factor of a few (many cores). The downside is the the updates on the screen may be a bit irregular. (Médéric Boquien)
- It turns out that elevating an array to some power is an especially slow operation in python. The `dustatt_calzleit` module has been optimised leading to a massive speed improvement. This speedup is especially large for models that do not include dust emission. (Médéric Boquien)
- Making copies of partially computed SED when storing them to the cache can be slow. Now we avoid making copies of the redshifted SED. The speedup should be especially noticeable when computing a set of models with numerous redshifts. (Médéric Boquien)
## 0.12.1 (2018-02-27) ## 0.12.1 (2018-02-27)
### Fixed ### Fixed
- The best fit could not be computed in photo-z mode because the redshift was negative. (Médéric Boquien) - The best fit could not be computed in photo-z mode because the redshift was negative. (Médéric Boquien)
......
...@@ -168,7 +168,8 @@ def build_filters(base): ...@@ -168,7 +168,8 @@ def build_filters(base):
# We normalise the filter and compute the pivot wavelength. If the # We normalise the filter and compute the pivot wavelength. If the
# filter is a pseudo-filter used to compute line fluxes, it should not # filter is a pseudo-filter used to compute line fluxes, it should not
# be normalised. # be normalised.
if not filter_name.startswith('PSEUDO'): if not (filter_name.startswith('PSEUDO') or
filter_name.startswith('linefilter')):
new_filter.normalise() new_filter.normalise()
else: else:
new_filter.pivot_wavelength = np.mean( new_filter.pivot_wavelength = np.mean(
...@@ -696,7 +697,10 @@ def build_nebular(base): ...@@ -696,7 +697,10 @@ def build_nebular(base):
print("Importing {}...".format(nebular_dir + 'lines.dat')) print("Importing {}...".format(nebular_dir + 'lines.dat'))
lines = np.genfromtxt(nebular_dir + 'lines.dat') lines = np.genfromtxt(nebular_dir + 'lines.dat')
wave_lines = np.genfromtxt(nebular_dir + 'line_wavelengths.dat') tmp = Table.read(nebular_dir + 'line_wavelengths.dat', format='ascii')
wave_lines = tmp['col1'].data
name_lines = tmp['col2'].data
print("Importing {}...".format(nebular_dir + 'continuum.dat')) print("Importing {}...".format(nebular_dir + 'continuum.dat'))
cont = np.genfromtxt(nebular_dir + 'continuum.dat') cont = np.genfromtxt(nebular_dir + 'continuum.dat')
...@@ -727,8 +731,8 @@ def build_nebular(base): ...@@ -727,8 +731,8 @@ def build_nebular(base):
spectra = lines[idx::6, :] spectra = lines[idx::6, :]
for logU, spectrum in zip(np.around(np.arange(-4., -.9, .1), 1), for logU, spectrum in zip(np.around(np.arange(-4., -.9, .1), 1),
spectra.T): spectra.T):
models_lines.append(NebularLines(metallicity, logU, wave_lines, models_lines.append(NebularLines(metallicity, logU, name_lines,
spectrum)) wave_lines, spectrum))
# Import continuum # Import continuum
for idx, metallicity in enumerate(metallicities): for idx, metallicity in enumerate(metallicities):
......
# linefilter.H-alpha-100
# energy
# H-alpha pseudo filter, R=100. The flux will not be corrected for any absorption line. Make sure the continuum part is not affected by emission lines. Feel free to adapt the filter if the resolution is not adequate.
6497.370000 -1.0
6530.185000 -1.0
6530.185001 1.0
6595.814999 1.0
6595.815000 -1.0
6628.630000 -1.0
# linefilter.H-alpha-50
# energy
# H-alpha pseudo filter, R=50. The flux will not be corrected for any absorption line. Make sure the continuum part is not affected by emission lines. Feel free to adapt the filter if the resolution is not adequate.
6431.740000 -1.0
6497.370000 -1.0
6497.370001 1.0
6628.629999 1.0
6628.630000 -1.0
6694.260000 -1.0
# linefilter.H-alpha-75
# energy
# H-alpha pseudo filter, R=75. The flux will not be corrected for any absorption line. Make sure the continuum part is not affected by emission lines. Feel free to adapt the filter if the resolution is not adequate.
6475.493333 -1.0
6519.246667 -1.0
6519.246668 1.0
6606.753332 1.0
6606.753333 -1.0
6650.506667 -1.0
# linefilter.H-beta-100
# energy
# H-beta pseudo filter, R=100. The flux will not be corrected for any absorption line. Make sure the continuum part is not affected by emission lines. Feel free to adapt the filter if the resolution is not adequate.
4812.390000 -1.0
4836.695000 -1.0
4836.695001 1.0
4885.304999 1.0
4885.305000 -1.0
4909.610000 -1.0
# linefilter.H-beta-125
# energy
# H-beta pseudo filter, R=125. The flux will not be corrected for any absorption line. Make sure the continuum part is not affected by emission lines. Feel free to adapt the filter if the resolution is not adequate.
4822.112000 -1.0
4841.556000 -1.0
4841.556001 1.0
4880.443999 1.0
4880.444000 -1.0
4899.888000 -1.0
# linefilter.H-beta-75
# energy
# H-beta pseudo filter, R=75. The flux will not be corrected for any absorption line. Make sure the continuum part is not affected by emission lines. Feel free to adapt the filter if the resolution is not adequate.
4796.186667 -1.0
4828.593333 -1.0
4828.593334 1.0
4893.406666 1.0
4893.406667 -1.0
4925.813333 -1.0
# linefilter.OIII-50
# energy
# OIII pseudo filter, R=50. The flux will not be corrected for any absorption line. Make sure the continuum part is not affected by emission lines. Feel free to adapt the filter if the resolution is not adequate.
4883.340000 -1.0
4933.170000 -1.0
4933.170001 1.0
5032.829999 1.0
5032.830000 -1.0
5082.660000 -1.0
# linefilter.OIII-75
# energy
# OIII pseudo filter, R=75. The flux will not be corrected for any absorption line. Make sure the continuum part is not affected by emission lines. Feel free to adapt the filter if the resolution is not adequate.
4916.560000 -1.0
4949.780000 -1.0
4949.780001 1.0
5016.219999 1.0
5016.220000 -1.0
5049.440000 -1.0
303.8 # He 2 303.8A 303.8,HeII-30.38
584.3 # He 1 584.3A 584.3,HeI-58.43
625.6 # He 1 625.6A 625.6,HeI-62.56
972.5 # H 1 972.5A Ly gamma (1-4) 972.5,Ly-gamma
977.0 # C 3 977.0A 977.0,CIII-97.70
1026 # H 1 1026A Ly beta (1-3) 1026,Ly-beta
1216 # H 1 1216A Ly alpha (1-2) 1216,Ly-alpha
1335 # C 2 1335A C II 1335,CII-133.5
1397 # TOTL 1397A Si IV 1394, 1397, 1398 1397,SiIV-139.7
1486 # TOTL 1486A N IV] 1486,NIV-148.6
1549 # TOTL 1549A C IV 1548, 1551 1549,CIV-154.9
1640 # He 2 1640A 1640,HeII-164.0
1665 # TOTL 1665A O III] 1661, 1666 1665,OIII-166.5
1750 # TOTL 1750A N III] 1750,NIII-175.0
1888 # TOTL 1888A Si III] 1883, 1892 1888,SiIII-188.8
1909 # TOTL 1909A C III] 1909,CIII-190.9
2326 # TOTL 2326A C II] 2326,CII-232.6
2400 # Fe 2 2400A 2400,FeII-240.0
2471 # O II 2471A 2471,OII-247.1
2798 # TOTL 2798A Mg II 2796, 2803 2798,MgII-279.8
2829 # He 1 2829A 2829,HeI-282.9
2836 # Fe 4 2836A 2836,FeIV-283.6
2945 # He 1 2945A 2945,HeI-294.5
3096 # Fe 4 3096A 3096,FeIV-309.6
3188 # He 1 3188A 3188,HeI-318.8
3203 # He 2 3203A 3203,HeII-320.3
3671 # Ca B 3671A H 24 (Case B) 3671,H-24
3674 # Ca B 3674A H 23 (Case B) 3674,H-23
3676 # Ca B 3676A H 22 (Case B) 3676,H-22
3679 # Ca B 3679A H 21 (Case B) 3679,H-21
3683 # Ca B 3683A H 20 (Case B) 3683,H-20
3687 # Ca B 3687A H 19 (Case B) 3687,H-19
3692 # Ca B 3692A H 18 (Case B) 3692,H-18
3697 # Ca B 3697A H 17 (Case B) 3697,H-17
3704 # Ca B 3704A H 16 (Case B) 3704,H-16
3712 # Ca B 3712A H 15 (Case B) 3712,H-15
3722 # Ca B 3722A H 14 (Case B) 3722,H-14
3727 # TOTL 3727A [OII] 3726, 3729 3727,OII-372.7
3734 # Ca B 3734A H 13 (Case B) 3734,H-13
3750 # Ca B 3750A H 12 (Case B) 3750,H-12
3771 # Ca B 3771A H 11 (Case B) 3771,H-11
3798 # H 1 3798A H 10 3798,H-10
3820 # He 1 3820A 3820,HeI-382.0
3835 # H 1 3835A H 9 3835,H-9
3869 # Ne 3 3869A 3869,NeIII-386.9
3889 # He 1 3889A 3889,HeI-388.9
3889 # H 1 3889A H 8 3889,H-8
3965 # He 1 3965A 3965,HeI-396.5
3968 # Ne 3 3968A 3968,NeIII-396.8
3970 # H 1 3970A H epsilon (2-7) 3970,H-epsilon
4026 # He 1 4026A 4026,HeI-402.6
4070 # S II 4070A 4070,SII-407.0
4074 # S 2 4074A 4074,SII-407.4
4102 # H 1 4102A H delta (2-6) 4102,H-delta
4300 # Fe 2 4300A 4300,FeII-430.0
4340 # H 1 4340A H gamma (2-5) 4340,H-gamma
4363 # TOTL 4363A 4363,OIII-436.3
4471 # He 1 4471A 4471,HeI-447.1
4659 # Fe 3 4659A 4659,FeIII-465.9
4686 # He 2 4686A 4686,HeII-468.6
4861 # H 1 4861A H beta (2-4) 4861,H-beta
4922 # He 1 4922A 4922,HeI-492.2
4959 # O 3 4959A 4959,OIII-495.9
5007 # O 3 5007A 5007,OIII-500.7
5016 # He 1 5016A 5016,HeI-501.6
5876 # He 1 5876A 5876,HeI-587.6
6300 # O 1 6300A 6300,OI-630.0
6312 # S 3 6312A 6312,SIII-631.2
6548 # N 2 6548A 6548,NII-654.8
6563 # H 1 6563A H alpha (2-3) 6563,H-alpha
6584 # N 2 6584A 6584,NII-658.4
6678 # He 1 6678A 6678,HeI-667.8
6716 # S II 6716A 6716,SII-671.6
6720 # S 2 6720A 6720,SII-672.0
6731 # S II 6731A 6731,SII-673.1
7065 # He 1 7065A 7065,HeI-706.5
7135 # Ar 3 7135A 7135,ArIII-713.5
7325 # TOTL 7325A [OII] 7323, 7332 7325,OII-732.5
7751 # Ar 3 7751A 7751,ArIII-775.1
8334 # Ca B 8334A Pa 24 (Case B) 8334,Pa-24
8346 # Ca B 8346A Pa 23 (Case B) 8346,Pa-23
8359 # Ca B 8359A Pa 22 (Case B) 8359,Pa-22
8374 # Ca B 8374A Pa 21 (Case B) 8374,Pa-21
8392 # Ca B 8392A Pa 20 (Case B) 8392,Pa-20
8413 # Ca B 8413A Pa 19 (Case B) 8413,Pa-19
8438 # Ca B 8438A Pa 18 (Case B) 8438,Pa-18
8467 # Ca B 8467A Pa 17 (Case B) 8467,Pa-17
8502 # Ca B 8502A Pa 16 (Case B) 8502,Pa-16
8545 # Ca B 8545A Pa 15 (Case B) 8545,Pa-15
8598 # Ca B 8598A Pa 14 (Case B) 8598,Pa-14
8665 # Ca B 8665A Pa 13 (Case B) 8665,Pa-13
8750 # Ca B 8750A Pa 12 8750,Pa-12
8863 # Ca B 8863A Pa 11 8863,Pa-11
9015 # H 1 9015A Pa 10 9015,Pa-10
9069 # S 3 9069A 9069,SIII-906.9
9229 # H 1 9229A Pa 9 9229,Pa9
9532 # S 3 9532A 9532,SIII-953.2
9546 # H 1 9546A Pa epsilon (3-8) 9546,Pa-epsilon
10050 # H 1 1.005m Pa delta (3-7) 10050,Pa-delta
10830 # TOTL 1.083m He I 10830,HeI-1.083
10940 # H 1 1.094m Pa gamma (3-6) 10940,Pa-gamma
12820 # H 1 1.282m Pa beta (3-5) 12820,Pa-beta
18750 # H 1 1.875m Pa alpha (3-4) 18750,Pa-alpha
21660 # H 1 2.166m Br gamma (4-7) 21660,Br-gamma
26250 # H 1 2.625m Br beta (4-6) 26250,Br-beta
40510 # H 1 4.051m Br alpha (4-5) 40510,Br-alpha
69800 # Ar 2 6.980m 69800,ArII-6.980
74580 # H 1 7.458m Pf alpha (5-6) 74580,Pf-alpha
90000 # Ar 3 9.000m 90000,ArIII-9.000
105100 # S 4 10.51m 105100,SIV-10.51
128100 # Ne 2 12.81m 128100,NeII-12.81
155500 # Ne 3 15.55m 155500,NeIII-15.55
186700 # S 3 18.67m 186700,SIII-18.67
334700 # S 3 33.47m 334700,SIII-33.47
348100 # Si 2 34.81m 348100,SiII-34.81
360100 # Ne 3 36.01m 360100,NeIII-36.01
518000 # O 3 51.80m 518000,OIII-51.80
572100 # N 3 57.21m 572100,NIII-57.21
631700 # O 1 63.17m 631700,OI-63.17
883300 # O 3 88.33m 883300,OIII-88.33
1217000 # N 2 121.7m 1217000,NII-121.7
1455000 # O 1 145.5m 1455000,OI-145.5
1576000 # C 2 157.6m 1576000,CII-157.6
2054000 # N 2 205.4m 2054000,NII-205.4
...@@ -64,8 +64,12 @@ def main(): ...@@ -64,8 +64,12 @@ def main():
# We set the sub processes start method to spawn because it solves # We set the sub processes start method to spawn because it solves
# deadlocks when a library cannot handle being used on two sides of a # deadlocks when a library cannot handle being used on two sides of a
# forked process. This happens on modern Macs with the Accelerate library # forked process. This happens on modern Macs with the Accelerate library
# for instance. # for instance. On Linux we should be pretty safe with a fork, which allows
mp.set_start_method('spawn') # to start processes much more rapidly.
if sys.platform.startswith('linux'):
mp.set_start_method('fork')
else:
mp.set_start_method('spawn')
parser = argparse.ArgumentParser() parser = argparse.ArgumentParser()
......
...@@ -27,11 +27,11 @@ reduced χ²) is given for each observation. ...@@ -27,11 +27,11 @@ reduced χ²) is given for each observation.
from collections import OrderedDict from collections import OrderedDict
import multiprocessing as mp import multiprocessing as mp
import time
import numpy as np import numpy as np
from .. import AnalysisModule from .. import AnalysisModule
from ..utils import Counter
from .workers import sed as worker_sed from .workers import sed as worker_sed
from .workers import init_sed as init_worker_sed from .workers import init_sed as init_worker_sed
from .workers import init_analysis as init_worker_analysis from .workers import init_analysis as init_worker_analysis
...@@ -95,28 +95,31 @@ class PdfAnalysis(AnalysisModule): ...@@ -95,28 +95,31 @@ class PdfAnalysis(AnalysisModule):
)) ))
]) ])
def _compute_models(self, conf, obs, params, iblock): def _compute_models(self, conf, obs, params, iblock):
models = ModelsManager(conf, obs, params, iblock) models = ModelsManager(conf, obs, params, iblock)
counter = Counter(len(params.blocks[iblock]), 50, 250)
initargs = (models, counter)
initargs = (models, time.time(), mp.Value('i', 0))
self._parallel_job(worker_sed, params.blocks[iblock], initargs, self._parallel_job(worker_sed, params.blocks[iblock], initargs,
init_worker_sed, conf['cores']) init_worker_sed, conf['cores'])
# Print the final value as it may not otherwise be printed
if counter.global_counter.value % 250 != 0:
counter.pprint(len(params.blocks[iblock]))
return models return models
def _compute_bayes(self, conf, obs, models): def _compute_bayes(self, conf, obs, models):
results = ResultsManager(models) results = ResultsManager(models)
initargs = (models, results, time.time(), mp.Value('i', 0)) initargs = (models, results, Counter(len(obs)))
self._parallel_job(worker_analysis, obs, initargs, self._parallel_job(worker_analysis, obs, initargs,
init_worker_analysis, conf['cores']) init_worker_analysis, conf['cores'])
return results return results
def _compute_best(self, conf, obs, params, results): def _compute_best(self, conf, obs, params, results):
initargs = (conf, params, obs, results, time.time(), initargs = (conf, params, obs, results, Counter(len(obs)))
mp.Value('i', 0))
self._parallel_job(worker_bestfit, obs, initargs, self._parallel_job(worker_bestfit, obs, initargs,
init_worker_bestfit, conf['cores']) init_worker_bestfit, conf['cores'])
...@@ -181,16 +184,16 @@ class PdfAnalysis(AnalysisModule): ...@@ -181,16 +184,16 @@ class PdfAnalysis(AnalysisModule):
# Rename the output directory if it exists # Rename the output directory if it exists
self.prepare_dirs() self.prepare_dirs()
# Store the grid of parameters in a manager to facilitate the
# computation of the models
params = ParametersManager(conf)
# Store the observations in a manager which sanitises the data, checks # Store the observations in a manager which sanitises the data, checks
# all the required fluxes are present, adding errors if needed, # all the required fluxes are present, adding errors if needed,
# discarding invalid fluxes, etc. # discarding invalid fluxes, etc.
obs = ObservationsManager(conf) obs = ObservationsManager(conf, params)
obs.save('observations') obs.save('observations')
# Store the grid of parameters in a manager to facilitate the
# computation of the models
params = ParametersManager(conf)
results = self._compute(conf, obs, params) results = self._compute(conf, obs, params)
results.best.analyse_chi2() results.best.analyse_chi2()
...@@ -201,7 +204,7 @@ class PdfAnalysis(AnalysisModule): ...@@ -201,7 +204,7 @@ class PdfAnalysis(AnalysisModule):
print("\nAnalysing the mock observations...") print("\nAnalysing the mock observations...")
# For the mock analysis we do not save the ancillary files. # For the mock analysis we do not save the ancillary files.
for k in ['best_sed', 'chi2', 'pdf']: for k in ['best_sed', 'chi2']:
conf['analysis_params']["save_{}".format(k)] = False conf['analysis_params']["save_{}".format(k)] = False
# We replace the observations with a mock catalogue.. # We replace the observations with a mock catalogue..
......
...@@ -17,9 +17,9 @@ The data file is used only to get the list of fluxes to be computed. ...@@ -17,9 +17,9 @@ The data file is used only to get the list of fluxes to be computed.
from collections import OrderedDict from collections import OrderedDict
import multiprocessing as mp import multiprocessing as mp
import time
from .. import AnalysisModule from .. import AnalysisModule
from ..utils import Counter
from .workers import init_fluxes as init_worker_fluxes from .workers import init_fluxes as init_worker_fluxes
from .workers import fluxes as worker_fluxes from .workers import fluxes as worker_fluxes
from ...managers.models import ModelsManager from ...managers.models import ModelsManager
...@@ -75,11 +75,16 @@ class SaveFluxes(AnalysisModule): ...@@ -75,11 +75,16 @@ class SaveFluxes(AnalysisModule):
nblocks)) nblocks))
models = ModelsManager(conf, obs, params, iblock) models = ModelsManager(conf, obs, params, iblock)
counter = Counter(len(params.blocks[iblock]), 50, 250)
initargs = (models, time.time(), mp.Value('i', 0)) initargs = (models, counter)
self._parallel_job(worker_fluxes, params.blocks[iblock], initargs, self._parallel_job(worker_fluxes, params.blocks[iblock], initargs,
init_worker_fluxes, conf['cores']) init_worker_fluxes, conf['cores'])
# Print the final value as it may not otherwise be printed
if counter.global_counter.value % 250 != 0:
counter.pprint(len(params.blocks[iblock]))
print("Saving the models ....") print("Saving the models ....")
models.save('models-block-{}'.format(iblock)) models.save('models-block-{}'.format(iblock))
......
...@@ -5,15 +5,13 @@ ...@@ -5,15 +5,13 @@
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt # Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Author: Yannick Roehlly & Médéric Boquien # Author: Yannick Roehlly & Médéric Boquien
import time
import numpy as np import numpy as np
from ...warehouse import SedWarehouse from ...warehouse import SedWarehouse
from ..utils import nothread from ..utils import nothread
def init_fluxes(models, t0, ncomputed): def init_fluxes(models, counter):
"""Initializer of the pool of processes. It is mostly used to convert """Initializer of the pool of processes. It is mostly used to convert
RawArrays into numpy arrays. The latter are defined as global variables to RawArrays into numpy arrays. The latter are defined as global variables to
be accessible from the workers. be accessible from the workers.
...@@ -22,29 +20,22 @@ def init_fluxes(models, t0, ncomputed): ...@@ -22,29 +20,22 @@ def init_fluxes(models, t0, ncomputed):
---------- ----------
models: ModelsManagers models: ModelsManagers
Manages the storage of the computed models (fluxes and properties). Manages the storage of the computed models (fluxes and properties).
t0: float counter: Counter class object
Time of the beginning of the computation. Counter for the number of models computed
ncomputed: Value
Number of computed models. Shared among workers.
""" """
global gbl_previous_idx, gbl_warehouse, gbl_models, gbl_obs global gbl_warehouse, gbl_models, gbl_obs, gbl_save, gbl_counter
global gbl_properties, gbl_save, gbl_t0, gbl_ncomputed
# Limit the number of threads to 1 if we use MKL in order to limit the # Limit the number of threads to 1 if we use MKL in order to limit the
# oversubscription of the CPU/RAM. # oversubscription of the CPU/RAM.
nothread() nothread()
gbl_previous_idx = -1
gbl_warehouse = SedWarehouse() gbl_warehouse = SedWarehouse()
gbl_models = models gbl_models = models
gbl_obs = models.obs gbl_obs = models.obs
gbl_properties = [prop[:-4] if prop.endswith('_log') else prop for prop in
models.propertiesnames]
gbl_save = models.conf['analysis_params']['save_sed'] gbl_save = models.conf['analysis_params']['save_sed']
gbl_t0 = t0 gbl_counter = counter
gbl_ncomputed = ncomputed
def fluxes(idx, midx): def fluxes(idx, midx):
...@@ -57,33 +48,25 @@ def fluxes(idx, midx): ...@@ -57,33 +48,25 @@ def fluxes(idx, midx):
Index of the model within the current block of models. Index of the model within the current block of models.
""" """
global gbl_previous_idx
if gbl_previous_idx > -1:
gbl_warehouse.partial_clear_cache(
gbl_models.params.index_module_changed(gbl_previous_idx, midx))
gbl_previous_idx = midx
sed = gbl_warehouse.get_sed(gbl_models.params.modules, sed = gbl_warehouse.get_sed(gbl_models.params.modules,
gbl_models.params.from_index(midx)) gbl_models.params.from_index(midx))
if 'sfh.age' in sed.info and sed.info['sfh.age'] > sed.info['universe.age']: if 'sfh.age' in sed.info and sed.info['sfh.age'] > sed.info['universe.age']:
gbl_models.fluxes[:, idx] = np.full(len(gbl_obs.bands), np.nan) for band in gbl_models.flux:
gbl_models.properties[:, idx] = np.full(len(gbl_properties), np.nan) gbl_models.flux[band][idx] = np.nan
for prop in gbl_models.extprop:
gbl_models.extprop[prop][idx] = np.nan
for prop in gbl_models.intprop:
gbl_models.intprop[prop][idx] = np.nan
else: else:
gbl_models.fluxes[:, idx] = [sed.compute_fnu(filter_) for band in gbl_models.flux:
for filter_ in gbl_obs.bands] gbl_models.flux[band][idx] = sed.compute_fnu(band)
gbl_models.properties[:, idx] = [sed.info[name] for prop in gbl_models.extprop:
for name in gbl_properties]