# Path: ugali/utils/healpix.py
# def read_map(filename, nest=False, hdu=None, h=False, verbose=True):
#     """Read a healpix map from a fits file.  Partial-sky files,
#     if properly identified, are expanded to full size and filled with UNSEEN.
#     Uses fitsio to mirror much (but not all) of the functionality of healpy.read_map
#     
#     Parameters:
#     -----------
#     filename : str 
#       the fits file name
#     nest : bool, optional
#       If True return the map in NEST ordering, otherwise in RING ordering;
#       use fits keyword ORDERING to decide whether conversion is needed or not
#       If None, no conversion is performed.
#     hdu : int, optional
#       the header number to look at (start at 0)
#     h : bool, optional
#       If True, return also the header. Default: False.
#     verbose : bool, optional
#       If True, print a number of diagnostic messages
#     
#     Returns
#     -------
#     m [, header] : array, optionally with header appended
#       The map read from the file, and the header if *h* is True.
#     """
#     import fitsio
#     data,hdr = fitsio.read(filename,header=True,ext=hdu)
# 
#     nside = int(hdr.get('NSIDE'))
#     if verbose: print('NSIDE = {0:d}'.format(nside))
# 
#     if not hp.isnsideok(nside):
#         raise ValueError('Wrong nside parameter.')
#     sz=hp.nside2npix(nside)
# 
#     ordering = hdr.get('ORDERING','UNDEF').strip()
#     if verbose: print('ORDERING = {0:s} in fits file'.format(ordering))
# 
#     schm = hdr.get('INDXSCHM', 'UNDEF').strip()
#     if verbose: print('INDXSCHM = {0:s}'.format(schm))
#     if schm == 'EXPLICIT':
#         partial = True
#     elif schm == 'IMPLICIT':
#         partial = False
# 
#     # monkey patch on a field method
#     fields = data.dtype.names
# 
#     # Could be done more efficiently (but complicated) by reordering first
#     if hdr['INDXSCHM'] == 'EXPLICIT':
#         m = hp.UNSEEN*np.ones(sz,dtype=data[fields[1]].dtype)
#         m[data[fields[0]]] = data[fields[1]]
#     else:
#         m = data[fields[0]].ravel()
# 
#     if (not hp.isnpixok(m.size) or (sz>0 and sz != m.size)) and verbose:
#         print('nside={0:d}, sz={1:d}, m.size={2:d}'.format(nside,sz,m.size))
#         raise ValueError('Wrong nside parameter.')
#     if nest is not None:
#         if nest and ordering.startswith('RING'):
#             idx = hp.nest2ring(nside,np.arange(m.size,dtype=np.int32))
#             m = m[idx]
#             if verbose: print('Ordering converted to NEST')
#         elif (not nest) and ordering.startswith('NESTED'):
#             idx = hp.ring2nest(nside,np.arange(m.size,dtype=np.int32))
#             m = m[idx]
#             if verbose: print('Ordering converted to RING')
# 
#     if h: return m, hdr
#     else: return m


import os
import copy
import collections
import yaml
import numpy as np
import scipy.interpolate
import healpy as hp
import fitsio
import astropy.io.fits as pyfits
import ugali.utils.projector
import ugali.utils.healpix
import ugali.analysis.source
import ugali.analysis.kernel
import ugali.analysis.imf
import ugali.analysis.results
import ugali.simulation.population
        import pylab
    import argparse
from collections import OrderedDict as odict
from ugali.isochrone import factory as isochrone_factory
from ugali.utils.healpix import read_map
    assert mc_source_id_start >= 1, "Starting mc_source_id must be >= 1" 
    assert n % n_chunk == 0, "Total number of satellites must be divisible by the chunk size"
    nside_pix = 256 # NSIDE = 128 -> 27.5 arcmin, NSIDE = 256 -> 13.7 arcmin 
    
    if not os.path.exists(tag): os.makedirs(tag)

    if isinstance(config,str): config = yaml.load(open(config))
    assert config['survey'] in ['des', 'ps1', 'lsst']

    infile_ebv = config['ebv']
    infile_fracdet = config['fracdet']
    infile_maglim_g = config['maglim_g']
    infile_maglim_r = config['maglim_r']
    infile_density = config['stellar_density']

    range_distance      = config.get('range_distance',[5., 500.])
    range_stellar_mass  = config.get('range_stellar_mass',[1.e1, 1.e6])
    range_r_physical    = config.get('range_r_physical',[1.e-3, 2.0])
    range_ellipticity   = config.get('range_ellipticity',[0.1, 0.8])
    range_position_angle= config.get('range_position_angle',[0.0, 180.0])
    choice_age          = config.get('choice_age',[10., 12.0, 13.5])
    choice_metal        = config.get('choice_metal',[0.00010,0.00020])
    dwarf_file          = config.get('known_dwarfs',None)

    m_density = np.load(infile_density)
    nside_density = hp.npix2nside(len(m_density))
    m_fracdet = read_map(infile_fracdet, nest=False) #.astype(np.float16)
    nside_fracdet = hp.npix2nside(len(m_fracdet))

    m_maglim_g = read_map(infile_maglim_g, nest=False) #.astype(np.float16)

# Path: ugali/utils/healpix.py
# def read_map(filename, nest=False, hdu=None, h=False, verbose=True):
#     """Read a healpix map from a fits file.  Partial-sky files,
#     if properly identified, are expanded to full size and filled with UNSEEN.
#     Uses fitsio to mirror much (but not all) of the functionality of healpy.read_map
#     
#     Parameters:
#     -----------
#     filename : str 
#       the fits file name
#     nest : bool, optional
#       If True return the map in NEST ordering, otherwise in RING ordering;
#       use fits keyword ORDERING to decide whether conversion is needed or not
#       If None, no conversion is performed.
#     hdu : int, optional
#       the header number to look at (start at 0)
#     h : bool, optional
#       If True, return also the header. Default: False.
#     verbose : bool, optional
#       If True, print a number of diagnostic messages
#     
#     Returns
#     -------
#     m [, header] : array, optionally with header appended
#       The map read from the file, and the header if *h* is True.
#     """
#     import fitsio
#     data,hdr = fitsio.read(filename,header=True,ext=hdu)
# 
#     nside = int(hdr.get('NSIDE'))
#     if verbose: print('NSIDE = {0:d}'.format(nside))
# 
#     if not hp.isnsideok(nside):
#         raise ValueError('Wrong nside parameter.')
#     sz=hp.nside2npix(nside)
# 
#     ordering = hdr.get('ORDERING','UNDEF').strip()
#     if verbose: print('ORDERING = {0:s} in fits file'.format(ordering))
# 
#     schm = hdr.get('INDXSCHM', 'UNDEF').strip()
#     if verbose: print('INDXSCHM = {0:s}'.format(schm))
#     if schm == 'EXPLICIT':
#         partial = True
#     elif schm == 'IMPLICIT':
#         partial = False
# 
#     # monkey patch on a field method
#     fields = data.dtype.names
# 
#     # Could be done more efficiently (but complicated) by reordering first
#     if hdr['INDXSCHM'] == 'EXPLICIT':
#         m = hp.UNSEEN*np.ones(sz,dtype=data[fields[1]].dtype)
#         m[data[fields[0]]] = data[fields[1]]
#     else:
#         m = data[fields[0]].ravel()
# 
#     if (not hp.isnpixok(m.size) or (sz>0 and sz != m.size)) and verbose:
#         print('nside={0:d}, sz={1:d}, m.size={2:d}'.format(nside,sz,m.size))
#         raise ValueError('Wrong nside parameter.')
#     if nest is not None:
#         if nest and ordering.startswith('RING'):
#             idx = hp.nest2ring(nside,np.arange(m.size,dtype=np.int32))
#             m = m[idx]
#             if verbose: print('Ordering converted to NEST')
#         elif (not nest) and ordering.startswith('NESTED'):
#             idx = hp.ring2nest(nside,np.arange(m.size,dtype=np.int32))
#             m = m[idx]
#             if verbose: print('Ordering converted to RING')
# 
#     if h: return m, hdr
#     else: return m


# Path: ugali/scratch/simulation/simulate_population.py
import os
import copy
import collections
import yaml
import numpy as np
import scipy.interpolate
import healpy as hp
import fitsio
import astropy.io.fits as pyfits
import ugali.utils.projector
import ugali.utils.healpix
import ugali.analysis.source
import ugali.analysis.kernel
import ugali.analysis.imf
import ugali.analysis.results
import ugali.simulation.population
        import pylab
    import argparse
from collections import OrderedDict as odict
from ugali.isochrone import factory as isochrone_factory
from ugali.utils.healpix import read_map

    assert mc_source_id_start >= 1, "Starting mc_source_id must be >= 1" 
    assert n % n_chunk == 0, "Total number of satellites must be divisible by the chunk size"
    nside_pix = 256 # NSIDE = 128 -> 27.5 arcmin, NSIDE = 256 -> 13.7 arcmin 
    
    if not os.path.exists(tag): os.makedirs(tag)

    if isinstance(config,str): config = yaml.load(open(config))
    assert config['survey'] in ['des', 'ps1', 'lsst']

    infile_ebv = config['ebv']
    infile_fracdet = config['fracdet']
    infile_maglim_g = config['maglim_g']
    infile_maglim_r = config['maglim_r']
    infile_density = config['stellar_density']

    range_distance      = config.get('range_distance',[5., 500.])
    range_stellar_mass  = config.get('range_stellar_mass',[1.e1, 1.e6])
    range_r_physical    = config.get('range_r_physical',[1.e-3, 2.0])
    range_ellipticity   = config.get('range_ellipticity',[0.1, 0.8])
    range_position_angle= config.get('range_position_angle',[0.0, 180.0])
    choice_age          = config.get('choice_age',[10., 12.0, 13.5])
    choice_metal        = config.get('choice_metal',[0.00010,0.00020])
    dwarf_file          = config.get('known_dwarfs',None)

    m_density = np.load(infile_density)
    nside_density = hp.npix2nside(len(m_density))
    m_fracdet = read_map(infile_fracdet, nest=False) #.astype(np.float16)
    nside_fracdet = hp.npix2nside(len(m_fracdet))

    m_maglim_g = read_map(infile_maglim_g, nest=False) #.astype(np.float16)

    m_maglim_r = read_map(infile_maglim_r, nest=False) #.astype(np.float16)