Source code for pyraeus.perfphotoz

## Copyright (C) 2013 APC CNRS Universite Paris Diderot 
## <cecile.cavet@apc.univ-paris7.fr> 
## <lejeune@apc.univ-paris7.fr>  
## <cecile.portello-roucelle@apc.univ-paris7.fr>
## 
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, see http://www.gnu.org/licenses/gpl.html

import numpy as np
#from pylab import *

#### Warning: formula to compute outliers are changing from author to author. 
####

[docs]def error(z_cat,z_reco):
    """ Comput the redshift error
      
    1. Formula from Hildebrandt et al., AA, 523, A31 (2010):
        delta_z  = z_cat - z_reco
    
    2. Formula from Ilbert et al., AA, 457, 841 (2006):
        delta_z = delta_z / (1. + z_cat)

    **Examples**

    >>> import numpy as np
    >>> print error(np.array([1.0,1.4]),np.array([1.1,1.5]))
    (array([-0.1, -0.1]), array([-0.05      , -0.04166667]))
    """
    
    delta_z      = z_cat - z_reco
    delta_z_norm = delta_z / (1. + z_cat)
    
    return delta_z, delta_z_norm 
    
[docs]def outliers(delta_z,rate_outliers):
    """ Comput the rate of outliers in %
    
    outliers = abs(delta_z) > 0.10 - 0.15
    eta = outliers * 100 / nb_outliers 
    """
    
    delta_z_abs  = np.abs(delta_z)  
    outliers     = len(delta_z_abs[delta_z_abs>rate_outliers])
    eta          = 100. * float(outliers) / float(len(delta_z_abs)) 
    
    return outliers, eta
        
[docs]def reject_outliers(z_cat,delta_z,rate_outliers):
    """ Reject the outlier objects from the data
    
    See the function outliers(arg1,arg2) for formula.
    """

    delta_z_abs = np.abs(delta_z) 
    z_cat_wo    = z_cat[delta_z_abs<=rate_outliers]
    delta_z_wo  = delta_z[delta_z_abs<=rate_outliers]

    return z_cat_wo, delta_z_wo
    
[docs]def bias(delta_z_wo):
    """ Comput the redshift bias on the data without outliers
    
    e_z = abs(median(delta_z_wo))
    """
    
    e_z  = np.abs(np.median(delta_z_wo)) 
    
    return e_z
    
[docs]def acuracy(delta_z_wo):
    """ Comput the redshift accuracy on the data without outliers (?)
    
    1. Root-mean-square (RMS) scatter:
    sigma_rms  = standard_deviation(delta_z_wo)
    
    2. Normalized median absolute deviation (NMAD):
    sigma_nmad = 1.48 * median(abs(delta_z_norm))
    """
    
    sigma_rms  = delta_z_wo.std()

    sigma_nmad = 1.48 * np.median(np.abs(delta_z_wo)) 
    
    return sigma_rms, sigma_nmad