"""
   bpz: Bayesian Photo-Z estimation
   Reference: Benitez 2000, ApJ, 536, p.571
   Usage:
   python bpz.py catalog.cat
   Needs a catalog.columns file which describes the contents of catalog.cat
"""
from __future__ import print_function
from __future__ import division

from builtins import str
from builtins import map
from builtins import input
from builtins import range
from past.utils import old_div
from useful import *
rolex = watch()
rolex.set()

#from Numeric import *
from numpy import *
from bpz_tools import *
from string import *
import os, glob, sys
import time
import pickle
import shelve

from coetools import pause, params_cl

class Printer():
    """Print things to stdout on one line dynamically"""
    def __init__(self, data):
        sys.stdout.write("\r\x1b[K" + data.__str__())
        sys.stdout.flush()

def seglist(vals, mask=None):
    """Split vals into lists based on mask > 0"""
    if mask == None:
        mask = greater(vals, 0)
    lists = []
    i = 0
    lastgood = False
    list1 = []
    for i in range(len(vals)):
        if mask[i] == False:
            if lastgood:
                lists.append(list1)
                list1 = []
            lastgood = False
        if mask[i]:
            list1.append(vals[i])
            lastgood = True

    if lastgood:
        lists.append(list1)
    return lists

# Initialization and definitions#

#Current directory
homedir = os.getcwd()

#Parameter definition
pars = params()

pars.d = {
    'SPECTRA': 'CWWSB4.list',  # template list
    #'PRIOR':   'hdfn_SB',      # prior name
    'PRIOR': 'hdfn_gen',  # prior name
    'NTYPES':
    None,  # Number of Elliptical, Spiral, and Starburst/Irregular templates  Default: 1,2,n-3
    'DZ': 0.01,  # redshift resolution
    'ZMIN': 0.01,  # minimum redshift
    'ZMAX': 10.,  # maximum redshift
    'MAG': 'yes',  # Data in magnitudes?
    'MIN_MAGERR': 0.001,  # minimum magnitude uncertainty --DC
    'ODDS': 0.95,  # Odds threshold: affects confidence limits definition
    'INTERP':
    0,  # Number of interpolated templates between each of the original ones
    'EXCLUDE': 'none',  # Filters to be excluded from the estimation
    'NEW_AB': 'no',  # If yes, generate new AB files even if they already exist
    'CHECK':
    'yes',  # Perform some checks, compare observed colors with templates, etc.
    'VERBOSE': 'yes',  # Print estimated redshifts to the standard output
    'PROBS':
    'no',  # Save all the galaxy probability distributions (it will create a very large file)
    'PROBS2':
    'no',  # Save all the galaxy probability distributions P(z,t) (but not priors) -- Compact
    'PROBS_LITE': 'yes',  # Save only the final probability distribution
    'GET_Z': 'yes',  # Actually obtain photo-z
    'ONLY_TYPE': 'no',  # Use spectroscopic redshifts instead of photo-z
    'MADAU': 'yes',  #Apply Madau correction to spectra
    'Z_THR': 0,  #Integrate probability for z>z_thr
    'COLOR': 'no',  #Use colors instead of fluxes
    'PLOTS': 'no',  #Don't produce plots
    'INTERACTIVE': 'yes',  #Don't query the user
    'PHOTO_ERRORS':
    'no',  #Define the confidence interval using only the photometric errors
    'MIN_RMS':
    0.05,  #"Intrinsic"  photo-z rms in dz /(1+z) (Change to 0.05 for templates from Benitez et al. 2004
    'N_PEAKS': 1,
    'MERGE_PEAKS': 'no',
    'CONVOLVE_P': 'yes',
    'P_MIN': 1e-2,
    'SED_DIR': sed_dir,
    'AB_DIR': ab_dir,
    'FILTER_DIR': fil_dir,
    'DELTA_M_0': 0.,
    'ZP_OFFSETS': 0.,
    'ZC': None,
    'FC': None,
    "ADD_SPEC_PROB": None,
    "ADD_CONTINUOUS_PROB": None,
    "NMAX": None  # Useful for testing
}

if pars.d['PLOTS'] == 'no': plots = 0

if plots:
    # If pylab installed show plots
    plots = 'pylab'
    try:
        import matplotlib
        matplotlib.use('TkAgg')
        from pylab import *
        # from coeplot2a import *
        plot([1])
        title('KILL THIS WINDOW!')
        show()
        ioff()
    except:
        try:
            from biggles import *
            plots = 'biggles'
        except:
            plots = 0

#Define the default values of the parameters
pars.d['INPUT'] = sys.argv[1]  # catalog with the photometry
obs_file = pars.d['INPUT']
root = os.path.splitext(pars.d['INPUT'])[0]
pars.d[
    'COLUMNS'] = root + '.columns'  # column information for the input catalog
pars.d['OUTPUT'] = root + '.bpz'  # output

nargs = len(sys.argv)

ipar = 2

if nargs > 2:  #Check for parameter file and update parameters
    if sys.argv[2] == '-P':
        pars.fromfile(sys.argv[3])
        ipar = 4
# Update the parameters using command line additions
#pars.fromcommandline(sys.argv[ipar:])
#for key in pars.d:
#    print key, pars.d[key]
#pause()
pars.d.update(
    params_cl())  # allows for flag only (no value after), e.g., -CHECK


def updateblank(var, ext):
    global pars
    if pars.d[var] in [None, 'yes']:
        pars.d[var] = root + '.' + ext


updateblank('CHECK', 'flux_comparison')
updateblank('PROBS_LITE', 'probs')
updateblank('PROBS', 'full_probs')
updateblank('PROBS2', 'chisq')

#if pars.d['CHECK'] in [None, 'yes']:
#    pars.d['CHECK'] = root+'.flux_comparison'

#This allows to change the auxiliary directories used by BPZ
if pars.d['SED_DIR'] != sed_dir:
    print("Changing sed_dir to ", pars.d['SED_DIR'])
    sed_dir = pars.d['SED_DIR']
    if sed_dir[-1] != '/': sed_dir += '/'
if pars.d['AB_DIR'] != ab_dir:
    print("Changing ab_dir to ", pars.d['AB_DIR'])
    ab_dir = pars.d['AB_DIR']
    if ab_dir[-1] != '/': ab_dir += '/'
if pars.d['FILTER_DIR'] != fil_dir:
    print("Changing fil_dir to ", pars.d['FILTER_DIR'])
    fil_dir = pars.d['FILTER_DIR']
    if fil_dir[-1] != '/': fil_dir += '/'

#Better safe than sorry
if pars.d['OUTPUT'] == obs_file or pars.d['PROBS'] == obs_file or pars.d[
        'PROBS2'] == obs_file or pars.d['PROBS_LITE'] == obs_file:
    print("This would delete the input file!")
    sys.exit()
if pars.d['OUTPUT'] == pars.d['COLUMNS'] or pars.d['PROBS_LITE'] == pars.d[
        'COLUMNS'] or pars.d['PROBS'] == pars.d['COLUMNS']:
    print("This would delete the .columns file!")
    sys.exit()

#Assign the intrinsin rms
if pars.d['SPECTRA'] == 'CWWSB.list':
    print('Setting the intrinsic rms to 0.067(1+z)')
    pars.d['MIN_RMS'] = 0.067

pars.d['MIN_RMS'] = float(pars.d['MIN_RMS'])
pars.d['MIN_MAGERR'] = float(pars.d['MIN_MAGERR'])
if pars.d['INTERACTIVE'] == 'no': interactive = 0
else: interactive = 1
if pars.d['VERBOSE'] == 'yes':
    print("Current parameters")
    view_keys(pars.d)
pars.d['N_PEAKS'] = int(pars.d['N_PEAKS'])
if pars.d["ADD_SPEC_PROB"] != None:
    specprob = 1
    specfile = pars.d["ADD_SPEC_PROB"]
    spec = get_2Darray(specfile)
    ns = spec.shape[1]
    if old_div(ns, 2) != (old_div(ns, 2.)):
        print("Number of columns in SPEC_PROB is odd")
        sys.exit()
    z_spec = spec[:, :old_div(ns, 2)]
    p_spec = spec[:, old_div(ns, 2):]
    # Write output file header
    header = "#ID "
    header += ns / 2 * " z_spec%i"
    header += ns / 2 * " p_spec%i"
    header += "\n"
    header = header % tuple(list(range(old_div(ns, 2))) + list(range(old_div(
        ns, 2))))
    specout = open(specfile.split()[0] + ".p_spec", "w")
    specout.write(header)
else:
    specprob = 0
pars.d['DELTA_M_0'] = float(pars.d['DELTA_M_0'])

#Some misc. initialization info useful for the .columns file
#nofilters=['M_0','OTHER','ID','Z_S','X','Y']
nofilters = ['M_0', 'OTHER', 'ID', 'Z_S']

#Numerical codes for nondetection, etc. in the photometric catalog
unobs = -99.  #Objects not observed
undet = 99.  #Objects not detected

#Define the z-grid
zmin = float(pars.d['ZMIN'])
zmax = float(pars.d['ZMAX'])
if zmin > zmax: raise 'zmin < zmax !'
dz = float(pars.d['DZ'])

linear = 1
if linear:
    z = arange(zmin, zmax + dz, dz)
else:
    if zmax != 0.:
        zi = zmin
        z = []
        while zi <= zmax:
            z.append(zi)
            zi = zi + dz * (1. + zi)
        z = array(z)
    else:
        z = array([0.])

#Now check the contents of the FILTERS,SED and A diBrectories

#Get the filters in stock
filters_db = []
filters_db = glob.glob(fil_dir + '*.res')
for i in range(len(filters_db)):
    filters_db[i] = os.path.basename(filters_db[i])
    filters_db[i] = filters_db[i][:-4]

    #Get the SEDs in stock
sed_db = []
sed_db = glob.glob(sed_dir + '*.sed')
for i in range(len(sed_db)):
    sed_db[i] = os.path.basename(sed_db[i])
    sed_db[i] = sed_db[i][:-4]

#Get the ABflux files in stock
ab_db = []
ab_db = glob.glob(ab_dir + '*.AB')
for i in range(len(ab_db)):
    ab_db[i] = os.path.basename(ab_db[i])
    ab_db[i] = ab_db[i][:-3]

#Get a list with the filter names and check whether they are in stock
col_file = pars.d['COLUMNS']
filters = get_str(col_file, 0)

for cosa in nofilters:
    if filters.count(cosa): filters.remove(cosa)

if pars.d['EXCLUDE'] != 'none':
    if type(pars.d['EXCLUDE']) == type(' '):
        pars.d['EXCLUDE'] = [pars.d['EXCLUDE']]
    for cosa in pars.d['EXCLUDE']:
        if filters.count(cosa): filters.remove(cosa)

for filter in filters:
    if filter[-4:] == '.res': filter = filter[:-4]
    if filter not in filters_db:
        print('filter ', filter, 'not in database at', fil_dir, ':')
        if ask('Print filters in database?'):
            for line in filters_db:
                print(line)
        sys.exit()

#Get a list with the spectrum names and check whether they're in stock
#Look for the list in the home directory first,
#if it's not there, look in the SED directory
spectra_file = os.path.join(homedir, pars.d['SPECTRA'])
if not os.path.exists(spectra_file):
    spectra_file = os.path.join(sed_dir, pars.d['SPECTRA'])

spectra = get_str(spectra_file, 0)
for i in range(len(spectra)):
    if spectra[i][-4:] == '.sed': spectra[i] = spectra[i][:-4]

nf = len(filters)
nt = len(spectra)
nz = len(z)

#Get the model fluxes
f_mod = zeros((nz, nt, nf)) * 0.
abfiles = []

for it in range(nt):
    for jf in range(nf):
        if filters[jf][-4:] == '.res': filtro = filters[jf][:-4]
        else: filtro = filters[jf]
        #model = join([spectra[it], filtro, 'AB'], '.')
        model = '.'.join([spectra[it], filtro, 'AB'])
        model_path = os.path.join(ab_dir, model)
        abfiles.append(model)
        #Generate new ABflux files if not present
        # or if new_ab flag on
        if pars.d['NEW_AB'] == 'yes' or model[:-3] not in ab_db:
            if spectra[it] not in sed_db:
                print('SED ', spectra[it], 'not in database at', sed_dir)
                #		for line in sed_db:
                #                    print line
                sys.exit()
        #print spectra[it],filters[jf]
            print('     Generating ', model, '....')
            ABflux(spectra[it], filtro, madau=pars.d['MADAU'])
            #z_ab=arange(0.,zmax_ab,dz_ab) #zmax_ab and dz_ab are def. in bpz_tools
            # abflux=f_z_sed(spectra[it],filters[jf], z_ab,units='nu',madau=pars.d['MADAU'])
            # abflux=clip(abflux,0.,1e400)
            # buffer=join(['#',spectra[it],filters[jf], 'AB','\n'])
            #for i in range(len(z_ab)):
            #	 buffer=buffer+join([`z_ab[i]`,`abflux[i]`,'\n'])
            #open(model_path,'w').write(buffer)
            #zo=z_ab
            #f_mod_0=abflux
            #else:
            #Read the data

        zo, f_mod_0 = get_data(model_path, (0, 1))
        #Rebin the data to the required redshift resolution
        f_mod[:, it, jf] = match_resol(zo, f_mod_0, z)
        #if sometrue(less(f_mod[:,it,jf],0.)):
        if less(f_mod[:, it, jf], 0.).any():
            print('Warning: some values of the model AB fluxes are <0')
            print('due to the interpolation ')
            print('Clipping them to f>=0 values')
            #To avoid rounding errors in the calculation of the likelihood
            f_mod[:, it, jf] = clip(f_mod[:, it, jf], 0., 1e300)

        #We forbid f_mod to take values in the (0,1e-100) interval
        #f_mod[:,it,jf]=where(less(f_mod[:,it,jf],1e-100)*greater(f_mod[:,it,jf],0.),0.,f_mod[:,it,jf])

        #Here goes the interpolacion between the colors
ninterp = int(pars.d['INTERP'])

ntypes = pars.d['NTYPES']
if ntypes == None:
    nt0 = nt
else:
    nt0 = list(ntypes)
    for i, nt1 in enumerate(nt0):
        print(i, nt1)
        nt0[i] = int(nt1)
    if (len(nt0) != 3) or (sum(nt0) != nt):
        print()
        print('%d ellipticals + %d spirals + %d ellipticals' % tuple(nt0))
        print('does not add up to %d templates' % nt)
        print('USAGE: -NTYPES nell,nsp,nsb')
        print('nell = # of elliptical templates')
        print('nsp  = # of spiral templates')
        print('nsb  = # of starburst templates')
        print(
            'These must add up to the number of templates in the SPECTRA list')
        print('Quitting BPZ.')
        sys.exit()

if ninterp:
    nti = nt + (nt - 1) * ninterp
    buffer = zeros((nz, nti, nf)) * 1.
    tipos = arange(0., float(nti), float(ninterp) + 1.)
    xtipos = arange(float(nti))
    for iz in arange(nz):
        for jf in range(nf):
            buffer[iz, :, jf] = match_resol(tipos, f_mod[iz, :, jf], xtipos)
    nt = nti
    f_mod = buffer

#for j in range(nf):
#    plot=FramedPlot()
#    for i in range(nt): plot.add(Curve(z,log(f_mod[:,i,j]+1e-40)))
#    plot.show()
#    ask('More?')

#Load all the parameters in the columns file to a dictionary
col_pars = params()
col_pars.fromfile(col_file)

# Read which filters are in which columns
flux_cols = []
eflux_cols = []
cals = []
zp_errors = []
zp_offsets = []
for filter in filters:
    datos = col_pars.d[filter]
    flux_cols.append(int(datos[0]) - 1)
    eflux_cols.append(int(datos[1]) - 1)
    cals.append(datos[2])
    zp_errors.append(datos[3])
    zp_offsets.append(datos[4])
zp_offsets = array(list(map(float, zp_offsets)))
if pars.d['ZP_OFFSETS']:
    zp_offsets += array(list(map(float, pars.d['ZP_OFFSETS'])))

flux_cols = tuple(flux_cols)
eflux_cols = tuple(eflux_cols)

#READ the flux and errors from obs_file
f_obs = get_2Darray(obs_file, flux_cols)
ef_obs = get_2Darray(obs_file, eflux_cols)

#Convert them to arbitrary fluxes if they are in magnitudes
if pars.d['MAG'] == 'yes':
    seen = greater(f_obs, 0.) * less(f_obs, undet)
    no_seen = equal(f_obs, undet)
    no_observed = equal(f_obs, unobs)
    todo = seen + no_seen + no_observed
    #The minimum photometric error is 0.01
    #ef_obs=ef_obs+seen*equal(ef_obs,0.)*0.001
    ef_obs = where(
        greater_equal(ef_obs, 0.), clip(ef_obs, pars.d['MIN_MAGERR'], 1e10),
        ef_obs)
    if add.reduce(add.reduce(todo)) != todo.shape[0] * todo.shape[1]:
        print('Objects with unexpected magnitudes!')
        print("""Allowed values for magnitudes are
	0<m<""" + repr(undet) + " m=" + repr(undet) + "(non detection), m=" + repr(
            unobs) + "(not observed)")
        for i in range(len(todo)):
            if not alltrue(todo[i, :]):
                print(i + 1, f_obs[i, :], ef_obs[i, :])
        sys.exit()

#Detected objects
    try:
        f_obs = where(seen, 10.**(-.4 * f_obs), f_obs)
    except OverflowError:
        print(
            'Some of the input magnitudes have values which are >700 or <-700')
        print('Purge the input photometric catalog')
        print('Minimum value', min(f_obs))
        print('Maximum value', max(f_obs))
        print('Indexes for minimum values', argmin(f_obs, 0.))
        print('Indexes for maximum values', argmax(f_obs, 0.))
        print('Bye.')
        sys.exit()

    try:
        ef_obs = where(seen, (10.**(.4 * ef_obs) - 1.) * f_obs, ef_obs)
    except OverflowError:
        print(
            'Some of the input magnitude errors have values which are >700 or <-700')
        print('Purge the input photometric catalog')
        print('Minimum value', min(ef_obs))
        print('Maximum value', max(ef_obs))
        print('Indexes for minimum values', argmin(ef_obs, 0.))
        print('Indexes for maximum values', argmax(ef_obs, 0.))
        print('Bye.')
        sys.exit()

    #print 'ef', ef_obs[0,:nf]
    #print 'f',  f_obs[1,:nf]
    #print 'ef', ef_obs[1,:nf]

    #Looked at, but not detected objects (mag=99.)
    #We take the flux equal to zero, and the error in the flux equal to the 1-sigma detection error.
    #If m=99, the corresponding error magnitude column in supposed to be dm=m_1sigma, to avoid errors
    #with the sign we take the absolute value of dm
    f_obs = where(no_seen, 0., f_obs)
    ef_obs = where(no_seen, 10.**(-.4 * abs(ef_obs)), ef_obs)

    #Objects not looked at (mag=-99.)
    f_obs = where(no_observed, 0., f_obs)
    ef_obs = where(no_observed, 0., ef_obs)

#Flux codes:
# If f>0 and ef>0 : normal objects
# If f==0 and ef>0 :object not detected
# If f==0 and ef==0: object not observed
#Everything else will crash the program

#Check that the observed error fluxes are reasonable
#if sometrue(less(ef_obs,0.)): raise 'Negative input flux errors'
if less(ef_obs, 0.).any():
    raise ValueError('Negative input flux errors')

f_obs = where(less(f_obs, 0.), 0., f_obs)  #Put non-detections to 0
ef_obs = where(
    less(f_obs, 0.), maximum(1e-100, f_obs + ef_obs),
    ef_obs)  # Error equivalent to 1 sigma upper limit

#if sometrue(less(f_obs,0.)) : raise 'Negative input fluxes'
seen = greater(f_obs, 0.) * greater(ef_obs, 0.)
no_seen = equal(f_obs, 0.) * greater(ef_obs, 0.)
no_observed = equal(f_obs, 0.) * equal(ef_obs, 0.)

todo = seen + no_seen + no_observed
if add.reduce(add.reduce(todo)) != todo.shape[0] * todo.shape[1]:
    print('Objects with unexpected fluxes/errors')

#Convert (internally) objects with zero flux and zero error(non observed)
#to objects with almost infinite (~1e108) error and still zero flux
#This will yield reasonable likelihoods (flat ones) for these objects
ef_obs = where(no_observed, 1e108, ef_obs)

#Include the zero point errors
zp_errors = array(list(map(float, zp_errors)))
zp_frac = e_mag2frac(zp_errors)
#zp_frac=10.**(.4*zp_errors)-1.
ef_obs = where(seen, sqrt(ef_obs * ef_obs + (zp_frac * f_obs)**2), ef_obs)
ef_obs = where(no_seen,
               sqrt(ef_obs * ef_obs + (zp_frac * (old_div(ef_obs, 2.)))**2),
               ef_obs)

#Add the zero-points offset
#The offsets are defined as m_new-m_old
zp_offsets = array(list(map(float, zp_offsets)))
zp_offsets = where(not_equal(zp_offsets, 0.), 10.**(-.4 * zp_offsets), 1.)
f_obs = f_obs * zp_offsets
ef_obs = ef_obs * zp_offsets

#Convert fluxes to AB if needed
for i in range(f_obs.shape[1]):
    if cals[i] == 'Vega':
        const = mag2flux(VegatoAB(0., filters[i]))
        f_obs[:, i] = f_obs[:, i] * const
        ef_obs[:, i] = ef_obs[:, i] * const
    elif cals[i] == 'AB':
        continue
    else:
        print('AB or Vega?. Check ' + col_file + ' file')
        sys.exit()

        #Get m_0 (if present)
if 'M_0' in col_pars.d:
    m_0_col = int(col_pars.d['M_0']) - 1
    m_0 = get_data(obs_file, m_0_col)
    m_0 += pars.d['DELTA_M_0']

#Get the objects ID (as a string)
if 'ID' in col_pars.d:
    #    print col_pars.d['ID']
    id_col = int(col_pars.d['ID']) - 1
    id = get_str(obs_file, id_col)
else:
    id = list(map(str, list(range(1, len(f_obs[:, 0]) + 1))))

#Get spectroscopic redshifts (if present)
if 'Z_S' in col_pars.d:
    z_s_col = int(col_pars.d['Z_S']) - 1
    z_s = get_data(obs_file, z_s_col)

#Get the X,Y coordinates
if 'X' in col_pars.d:
    datos = col_pars.d['X']
    if len(datos) == 1:  # OTHERWISE IT'S A FILTER!
        x_col = int(col_pars.d['X']) - 1
        x = get_data(obs_file, x_col)
if 'Y' in col_pars.d:
    datos = col_pars.d['Y']
    if len(datos) == 1:  # OTHERWISE IT'S A FILTER!
        y_col = int(datos) - 1
        y = get_data(obs_file, y_col)

#If 'check' on, initialize some variables
check = pars.d['CHECK']

# This generates a file with m,z,T and observed/expected colors
#if check=='yes': pars.d['FLUX_COMPARISON']=root+'.flux_comparison'

checkSED = check != 'no'

ng = f_obs.shape[0]
if checkSED:
    # PHOTOMETRIC CALIBRATION CHECK
    #r=zeros((ng,nf),float)+1.
    #dm=zeros((ng,nf),float)+1.
    #w=r*0.
    # Defaults: r=1, dm=1, w=0
    frat = ones((ng, nf), float)
    dmag = ones((ng, nf), float)
    fw = zeros((ng, nf), float)

#Visualize the colors of the galaxies and the templates

#When there are spectroscopic redshifts available
if interactive and 'Z_S' in col_pars.d and plots and checkSED and ask(
        'Plot colors vs spectroscopic redshifts?'):
    color_m = zeros((nz, nt, nf - 1)) * 1.
    if plots == 'pylab':
        figure(1)
    nrows = 2
    ncols = old_div((nf - 1), nrows)
    if (nf - 1) % nrows: ncols += 1
    for i in range(nf - 1):
        ##plot=FramedPlot()
        # Check for overflows
        fmu = f_obs[:, i + 1]
        fml = f_obs[:, i]
        good = greater(fml, 1e-100) * greater(fmu, 1e-100)
        zz, fmu, fml = multicompress(good, (z_s, fmu, fml))
        colour = old_div(fmu, fml)
        colour = clip(colour, 1e-5, 1e5)
        colour = 2.5 * log10(colour)
        if plots == 'pylab':
            subplot(nrows, ncols, i + 1)
            plot(zz, colour, "bo")
        elif plots == 'biggles':
            d = Points(zz, colour, color='blue')
            plot.add(d)
        for it in range(nt):
            #Prevent overflows
            fmu = f_mod[:, it, i + 1]
            fml = f_mod[:, it, i]
            good = greater(fml, 1e-100)
            zz, fmu, fml = multicompress(good, (z, fmu, fml))
            colour = old_div(fmu, fml)
            colour = clip(colour, 1e-5, 1e5)
            colour = 2.5 * log10(colour)
            if plots == 'pylab':
                plot(zz, colour, "r")
            elif plots == 'biggles':
                d = Curve(zz, colour, color='red')
                plot.add(d)
        if plots == 'pylab':
            xlabel(r'$z$')
            ylabel('%s - %s' % (filters[i], filters[i + 1]))
        elif plots == 'biggles':
            plot.xlabel = r'$z$'
            plot.ylabel = '%s - %s' % (filters[i], filters[i + 1])
            plot.save_as_eps('%s-%s.eps' % (filters[i], filters[i + 1]))
            plot.show()
    if plots == 'pylab':
        show()
        inp = eval(input('Hit Enter to continue.'))

#Get other information which will go in the output file (as strings)
if 'OTHER' in col_pars.d:
    if col_pars.d['OTHER'] != 'all':
        other_cols = col_pars.d['OTHER']
        if type(other_cols) == type((2, )):
            other_cols = tuple(map(int, other_cols))
        else:
            other_cols = (int(other_cols), )
        other_cols = [x - 1 for x in other_cols]
        n_other = len(other_cols)
    else:
        n_other = get_2Darray(obs_file, cols='all', nrows=1).shape[1]
        other_cols = list(range(n_other))

    others = get_str(obs_file, other_cols)

    if len(other_cols) > 1:
        other = []
        for j in range(len(others[0])):
            lista = []
            for i in range(len(others)):
                lista.append(others[i][j])
            other.append(join(lista))
    else:
        other = others

if pars.d['GET_Z'] == 'no': get_z = 0
else: get_z = 1

#Prepare the output file
out_name = pars.d['OUTPUT']
if get_z:
    if os.path.exists(out_name):
        os.system('cp %s %s.bak' % (out_name, out_name))
        print("File %s exists. Copying it to %s.bak" % (out_name, out_name))
    output = open(out_name, 'w')

if pars.d['PROBS_LITE'] == 'no': save_probs = 0
else: save_probs = 1

if pars.d['PROBS'] == 'no': save_full_probs = 0
else: save_full_probs = 1

if pars.d['PROBS2'] == 'no': save_probs2 = 0
else: save_probs2 = 1

#Include some header information

#   File name and the date...
time_stamp = time.ctime(time.time())
if get_z: output.write('## File ' + out_name + '  ' + time_stamp + '\n')

#and also the parameters used to run bpz...
if get_z: output.write("""##
##Parameters used to run BPZ:
##
""")
claves = list(pars.d.keys())
claves.sort()
for key in claves:
    if type(pars.d[key]) == type((1, )):
        cosa = join(list(pars.d[key]), ',')
    else:
        cosa = str(pars.d[key])
    if get_z: output.write('##' + key.upper() + '=' + cosa + '\n')

if save_full_probs:
    #Shelve some info on the run
    full_probs = shelve.open(pars.d['PROBS'])
    full_probs['TIME'] = time_stamp
    full_probs['PARS'] = pars.d

if save_probs:
    probs = open(pars.d['PROBS_LITE'], 'w')
    probs.write('# ID  p_bayes(z)  where z=arange(%.4f,%.4f,%.4f) \n' %
                (zmin, zmax + dz, dz))

if save_probs2:
    probs2 = open(pars.d['PROBS2'], 'w')
    probs2.write(
        '# id t  z1    P(z1) P(z1+dz) P(z1+2*dz) ...  where dz = %.4f\n' % dz)
    #probs2.write('# ID\n')
    #probs2.write('# t  z1  P(z1)  P(z1+dz)  P(z1+2*dz) ...  where dz = %.4f\n' % dz)

    #Use a empirical prior?
tipo_prior = pars.d['PRIOR']
useprior = 0
if 'M_0' in col_pars.d:
    has_mags = 1
else:
    has_mags = 0
if has_mags and tipo_prior != 'none' and tipo_prior != 'flat':
    useprior = 1

#Add cluster 'spikes' to the prior?
cluster_prior = 0.
if pars.d['ZC']:
    cluster_prior = 1
    if type(pars.d['ZC']) == type(""): zc = array([float(pars.d['ZC'])])
    else: zc = array(list(map(float, pars.d['ZC'])))
    if type(pars.d['FC']) == type(""): fc = array([float(pars.d['FC'])])
    else: fc = array(list(map(float, pars.d['FC'])))

    fcc = add.reduce(fc)
    if fcc > 1.:
        print(ftc)
        raise 'Too many galaxies in clusters!'
    pi_c = zeros((nz, nt)) * 1.
    #Go over the different cluster spikes
    for i in range(len(zc)):
        #We define the cluster within dz=0.01 limits
        cluster_range = less_equal(abs(z - zc[i]), .01) * 1.
        #Clip values to avoid overflow
        exponente = clip(-(z - zc[i])**2 / 2. / (0.00333)**2, -700., 0.)
        #Outside the cluster range g is 0
        g = exp(exponente) * cluster_range
        norm = add.reduce(g)
        pi_c[:, 0] = pi_c[:, 0] + g / norm * fc[i]

    #Go over the different types
    print('We only apply the cluster prior to the early type galaxies')
    for i in range(1, 3 + 2 * ninterp):
        pi_c[:, i] = pi_c[:, i] + pi_c[:, 0]

#Output format
format = '%' + repr(maximum(5, len(id[0]))) + 's'  #ID format
format = format + pars.d[
    'N_PEAKS'] * ' %.3f %.3f  %.3f %.3f %.5f' + ' %.3f %.3f %10.3f'

#Add header with variable names to the output file
sxhdr = """##
##Column information
##
# 1 ID"""
k = 1

if pars.d['N_PEAKS'] > 1:
    for j in range(pars.d['N_PEAKS']):
        sxhdr += """
# %i Z_B_%i
# %i Z_B_MIN_%i
# %i Z_B_MAX_%i
# %i T_B_%i
# %i ODDS_%i""" % (k + 1, j + 1, k + 2, j + 1, k + 3, j + 1, k + 4, j + 1,
                   k + 5, j + 1)
        k += 5
else:
    sxhdr += """
# %i Z_B
# %i Z_B_MIN
# %i Z_B_MAX
# %i T_B
# %i ODDS""" % (k + 1, k + 2, k + 3, k + 4, k + 5)
    k += 5

sxhdr += """
# %i Z_ML
# %i T_ML
# %i CHI-SQUARED\n""" % (k + 1, k + 2, k + 3)

nh = k + 4
if 'Z_S' in col_pars.d:
    sxhdr = sxhdr + '# %i Z_S\n' % nh
    format = format + '  %.3f'
    nh += 1
if has_mags:
    format = format + '  %.3f'
    sxhdr = sxhdr + '# %i M_0\n' % nh
    nh += 1
if 'OTHER' in col_pars.d:
    sxhdr = sxhdr + '# %i OTHER\n' % nh
    format = format + ' %s'
    nh += n_other

#print sxhdr

if get_z: output.write(sxhdr + '##\n')

odds_i = float(pars.d['ODDS'])
oi = inv_gauss_int(odds_i)

print(odds_i, oi)

#Proceed to redshift estimation

if checkSED: buffer_flux_comparison = ""

if pars.d['CONVOLVE_P'] == 'yes':
    # Will Convolve with a dz=0.03 gaussian to make probabilities smoother
    # This is necessary; if not there are too many close peaks
    sigma_g = 0.03
    x = arange(-3. * sigma_g, 3. * sigma_g + old_div(dz, 10.),
               dz)  # made symmetric --DC
    gaus = exp(-(old_div(x, sigma_g))**2)

if pars.d["NMAX"] != None: ng = int(pars.d["NMAX"])
for ig in range(ng):
    currentPercent = ig / ng * 100
    status = "{:.3f}% of {} completed.".format(currentPercent, ng)
    Printer(status)

    #Don't run BPZ on galaxies with have z_s > z_max
    #if col_pars.d.has_key('Z_S'):
    #    if z_s[ig]<9.9 and z_s[ig]>zmax : continue
    if not get_z: continue
    if pars.d['COLOR'] == 'yes':
        likelihood = p_c_z_t_color(f_obs[ig, :nf], ef_obs[ig, :nf],
                                   f_mod[:nz, :nt, :nf])
    else:
        likelihood = p_c_z_t(f_obs[ig, :nf], ef_obs[ig, :nf],
                             f_mod[:nz, :nt, :nf])

    if 0:
        print(f_obs[ig, :nf])
        print(ef_obs[ig, :nf])

    iz_ml = likelihood.i_z_ml
    t_ml = likelihood.i_t_ml
    red_chi2 = old_div(likelihood.min_chi2, float(nf - 1.))
    #p=likelihood.Bayes_likelihood
    #likelihood.various_plots()
    #print 'FULL BAYESAIN LIKELIHOOD'
    p = likelihood.likelihood
    if not ig:
        print('ML * prior -- NOT QUITE BAYESIAN')

    if pars.d[
            'ONLY_TYPE'] == 'yes':  #Use only the redshift information, no priors
        p_i = zeros((nz, nt)) * 1.
        j = searchsorted(z, z_s[ig])
        #print j,nt,z_s[ig]
        try:
            p_i[j, :] = old_div(1., float(nt))
        except IndexError:
            pass
    else:
        if useprior:
            if pars.d['PRIOR'] == 'lensing':
                p_i = prior(z, m_0[ig], tipo_prior, nt0, ninterp, x[ig], y[ig])
            else:
                p_i = prior(z, m_0[ig], tipo_prior, nt0, ninterp)
        else:
            p_i = old_div(ones((nz, nt), float), float(nz * nt))
        if cluster_prior: p_i = (1. - fcc) * p_i + pi_c

    if save_full_probs:
        full_probs[id[ig]] = [z, p_i[:nz, :nt], p[:nz, :nt], red_chi2]

    #Multiply the prior by the likelihood to find the final probability
    pb = p_i[:nz, :nt] * p[:nz, :nt]

    #plo=FramedPlot()
    #for i in range(p.shape[1]):
    #    plo.add(Curve(z,p_i[:nz,i]/sum(sum(p_i[:nz,:]))))
    #for i in range(p.shape[1]):
    #    plo.add(Curve(z,p[:nz,i]/sum(sum(p[:nz,:])),color='red'))
    #plo.add(Curve(z,pb[:nz,-1]/sum(pb[:nz,-1]),color='blue'))
    #plo.show()
    #ask('More?')

    #Convolve with a gaussian of width \sigma(1+z) to take into
    #accout the intrinsic scatter in the redshift estimation 0.06*(1+z)
    #(to be done)

    #Estimate the bayesian quantities
    p_bayes = add.reduce(pb[:nz, :nt], -1)
    #print p_bayes.shape
    #print argmax(p_bayes)
    #print p_bayes[300:310]

    #Convolve with a gaussian
    if pars.d['CONVOLVE_P'] == 'yes' and pars.d['ONLY_TYPE'] == 'no':
        #print 'GAUSS CONV'
        p_bayes = convolve(p_bayes, gaus, 1)
        #print 'gaus', gaus
        #print p_bayes.shape
        #print argmax(p_bayes)
        #print p_bayes[300:310]

        # Eliminate all low level features in the prob. distribution
    pmax = max(p_bayes)
    p_bayes = where(
        greater(p_bayes, pmax * float(pars.d['P_MIN'])), p_bayes, 0.)

    norm = add.reduce(p_bayes)
    p_bayes = old_div(p_bayes, norm)

    if specprob:
        p_spec[ig, :] = match_resol(z, p_bayes, z_spec[ig, :]) * p_spec[ig, :]
        norma = add.reduce(p_spec[ig, :])
        if norma == 0.: norma = 1.
        p_spec[ig, :] /= norma
        #vyjod=tuple([id[ig]]+list(z_spec[ig,:])+list(p_spec[ig,:])+[z_s[ig],
        #                int(float(other[ig]))])
        vyjod = tuple([id[ig]] + list(z_spec[ig, :]) + list(p_spec[ig, :]))
        formato = "%s " + 5 * " %.4f"
        formato += 5 * " %.3f"
        #formato+="  %4f %i"
        formato += "\n"
        print(formato % vyjod)
        specout.write(formato % vyjod)

    if pars.d['N_PEAKS'] > 1:
        # Identify  maxima and minima in the final probability
        g_max = less(p_bayes[2:], p_bayes[1:-1]) * less(p_bayes[:-2],
                                                        p_bayes[1:-1])
        g_min = greater(p_bayes[2:], p_bayes[1:-1]) * greater(p_bayes[:-2],
                                                              p_bayes[1:-1])

        g_min += equal(p_bayes[1:-1], 0.) * greater(p_bayes[2:], 0.)
        g_min += equal(p_bayes[1:-1], 0.) * greater(p_bayes[:-2], 0.)

        i_max = compress(g_max, arange(nz - 2)) + 1
        i_min = compress(g_min, arange(nz - 2)) + 1

        # Check that the first point and the last one are not minima or maxima,
        # if they are, add them to the index arrays

        if p_bayes[0] > p_bayes[1]:
            i_max = concatenate([[0], i_max])
            i_min = concatenate([[0], i_min])
        if p_bayes[-1] > p_bayes[-2]:
            i_max = concatenate([i_max, [nz - 1]])
            i_min = concatenate([i_min, [nz - 1]])
        if p_bayes[0] < p_bayes[1]:
            i_min = concatenate([[0], i_min])
        if p_bayes[-1] < p_bayes[-2]:
            i_min = concatenate([i_min, [nz - 1]])

        p_max = take(p_bayes, i_max)
        #p_min=take(p_bayes,i_min)
        p_tot = []
        z_peaks = []
        t_peaks = []
        # Sort them by probability values
        p_max, i_max = multisort(old_div(1., p_max), (p_max, i_max))
        # For each maximum, define the minima which sandwich it
        # Assign minima to each maximum
        jm = searchsorted(i_min, i_max)
        p_max = list(p_max)

        for i in range(len(i_max)):
            z_peaks.append([z[i_max[i]], z[i_min[jm[i] - 1]], z[i_min[jm[i]]]])
            t_peaks.append(argmax(pb[i_max[i], :nt]))
            p_tot.append(sum(p_bayes[i_min[jm[i] - 1]:i_min[jm[i]]]))
            # print z_peaks[-1][0],f_mod[i_max[i],t_peaks[-1]-1,:nf]

        if ninterp:
            t_peaks = list(old_div(array(t_peaks), (1. + ninterp)))

        if pars.d['MERGE_PEAKS'] == 'yes':
            # Merge peaks which are very close 0.03(1+z)
            merged = []
            for k in range(len(z_peaks)):
                for j in range(len(z_peaks)):
                    if j > k and k not in merged and j not in merged:
                        if abs(z_peaks[k][0] - z_peaks[j][0]) < 0.06 * (
                                1. + z_peaks[j][0]):
                            # Modify the element which receives the accretion
                            z_peaks[k][1] = minimum(z_peaks[k][1],
                                                    z_peaks[j][1])
                            z_peaks[k][2] = maximum(z_peaks[k][2],
                                                    z_peaks[j][2])
                            p_tot[k] += p_tot[j]
                            # Put the merged element in the list
                            merged.append(j)

                            #print merged
                            # Clean up
            copia = p_tot[:]
            for j in merged:
                p_tot.remove(copia[j])
            copia = z_peaks[:]
            for j in merged:
                z_peaks.remove(copia[j])
            copia = t_peaks[:]
            for j in merged:
                t_peaks.remove(copia[j])
            copia = p_max[:]
            for j in merged:
                p_max.remove(copia[j])

        if sum(array(p_tot)) != 1.:
            p_tot = old_div(array(p_tot), sum(array(p_tot)))

            # Define the peak
    iz_b = argmax(p_bayes)
    zb = z[iz_b]
    # OKAY, NOW THAT GAUSSIAN CONVOLUTION BUG IS FIXED
    # if pars.d['ONLY_TYPE']=='yes': zb=zb-dz/2. #This corrects a small bias
    # else: zb=zb-dz #This corrects another small bias --DC

    #Integrate within a ~ oi*sigma interval to estimate
    # the odds. (based on a sigma=pars.d['MIN_RMS']*(1+z))
    #Look for the number of sigma corresponding
    #to the odds_i confidence limit

    zo1 = zb - oi * pars.d['MIN_RMS'] * (1. + zb)
    zo2 = zb + oi * pars.d['MIN_RMS'] * (1. + zb)
    if pars.d['Z_THR'] > 0:
        zo1 = float(pars.d['Z_THR'])
        zo2 = float(pars.d['ZMAX'])
    o = odds(p_bayes[:nz], z, zo1, zo2)

    # Integrate within the same odds interval to find the type
    # izo1=maximum(0,searchsorted(z,zo1)-1)
    # izo2=minimum(nz,searchsorted(z,zo2))
    # t_b=argmax(add.reduce(p[izo1:izo2,:nt],0))

    it_b = argmax(pb[iz_b, :nt])
    t_b = it_b + 1

    if ninterp:
        tt_b = old_div(float(it_b), (1. + ninterp))
        tt_ml = old_div(float(t_ml), (1. + ninterp))
    else:
        tt_b = it_b
        tt_ml = t_ml

    if max(pb[iz_b, :]) < 1e-300:
        print('NO CLEAR BEST t_b; ALL PROBABILITIES ZERO')
        t_b = -1.
        tt_b = -1.

    #print it_b, t_b, tt_b, pb.shape

    if 0:
        print(f_mod[iz_b, it_b, :nf])

        print(min(ravel(p_i)), max(ravel(p_i)))
        print(min(ravel(p)), max(ravel(p)))
        print(p_i[iz_b, :])
        print(p[iz_b, :])
        print(p_i[iz_b, it_b])  # prior
        print(p[iz_b, it_b])  # chisq
        print(likelihood.likelihood[iz_b, it_b])
        print(likelihood.chi2[iz_b, it_b])
        print(likelihood.ftt[iz_b, it_b])
        print(likelihood.foo)

        print()
        print('t_b', t_b)
        print('iz_b', iz_b)
        print('nt', nt)
        print(max(ravel(pb)))
        impb = argmax(ravel(pb))
        impbz = old_div(impb, nt)
        impbt = impb % nt
        print(impb, impbz, impbt)
        print(ravel(pb)[impb])
        print(pb.shape, (nz, nt))
        print(pb[impbz, impbt])
        print(pb[iz_b, it_b])
        print('z, t', z[impbz], t_b)
        print(t_b)

    # Redshift confidence limits
    z1, z2 = interval(p_bayes[:nz], z, odds_i)
    if pars.d['PHOTO_ERRORS'] == 'no':
        zo1 = zb - oi * pars.d['MIN_RMS'] * (1. + zb)
        zo2 = zb + oi * pars.d['MIN_RMS'] * (1. + zb)
        if zo1 < z1: z1 = maximum(0., zo1)
        if zo2 > z2: z2 = zo2

    # Print output

    if pars.d['N_PEAKS'] == 1:
        salida = [id[ig], zb, z1, z2, tt_b + 1, o, z[iz_ml], tt_ml + 1,
                  red_chi2]
    else:
        salida = [id[ig]]
        for k in range(pars.d['N_PEAKS']):
            if k <= len(p_tot) - 1:
                salida = salida + list(z_peaks[k]) + [t_peaks[k] + 1, p_tot[k]]
            else:
                salida += [-1., -1., -1., -1., -1.]
        salida += [z[iz_ml], tt_ml + 1, red_chi2]

    if 'Z_S' in col_pars.d: salida.append(z_s[ig])
    if has_mags: salida.append(m_0[ig] - pars.d['DELTA_M_0'])
    if 'OTHER' in col_pars.d: salida.append(other[ig])

    if get_z: output.write(format % tuple(salida) + '\n')
    if pars.d['VERBOSE'] == 'yes': print(format % tuple(salida))

    #try:
    #    if sometrue(greater(z_peaks,7.5)):
    #        connect(z,p_bayes)
    #        ask('More?')
    #except:
    #    pass

    odd_check = odds_i

    if checkSED:
        ft = f_mod[iz_b, it_b, :]
        fo = f_obs[ig, :]
        efo = ef_obs[ig, :]
        dfosq = (old_div((ft - fo), efo))**2
        if 0:
            print(ft)
            print(fo)
            print(efo)
            print(dfosq)
            pause()
        factor = ft / efo / efo
        ftt = add.reduce(ft * factor)
        fot = add.reduce(fo * factor)
        am = old_div(fot, ftt)
        ft = ft * am
        if 0:
            print(factor)
            print(ftt)
            print(fot)
            print(am)
            print(ft)
            print()
            pause()

        flux_comparison = [id[ig], m_0[ig], z[iz_b], t_b, am] + list(
            concatenate([ft, fo, efo]))
        nfc = len(flux_comparison)

        format_fc = '%s  %.2f  %.2f   %i' + (nfc - 4) * '   %.3e' + '\n'
        buffer_flux_comparison = buffer_flux_comparison + format_fc % tuple(
            flux_comparison)
        if o >= odd_check:
            # PHOTOMETRIC CALIBRATION CHECK
            # Calculate flux ratios, but only for objects with ODDS >= odd_check
            #  (odd_check = 0.95 by default)
            # otherwise, leave weight w = 0 by default
            eps = 1e-10
            frat[ig, :] = divsafe(fo, ft, inf=eps, nan=eps)
            #fw[ig,:] = greater(fo, 0)
            fw[ig, :] = divsafe(fo, efo, inf=1e8, nan=0)
            fw[ig, :] = clip(fw[ig, :], 0, 100)
            #print fw[ig,:]
            #print

        if 0:
            bad = less_equal(ft, 0.)
            #Avoid overflow by setting r to 0.
            fo = where(bad, 0., fo)
            ft = where(bad, 1., ft)
            r[ig, :] = old_div(fo, ft)
            try:
                dm[ig, :] = -flux2mag(old_div(fo, ft))
            except:
                dm[ig, :] = -100
    # Clip ratio between 0.01 & 100
            r[ig, :] = where(greater(r[ig, :], 100.), 100., r[ig, :])
            r[ig, :] = where(less_equal(r[ig, :], 0.), 0.01, r[ig, :])
            #Weight by flux
            w[ig, :] = where(greater(fo, 0.), 1, 0.)
            #w[ig,:]=where(greater(fo,0.),fo,0.)
            #print fo
            #print r[ig,:]
            #print
            # This is no good becasue r is always > 0 (has been clipped that way)
            #w[ig,:]=where(greater(r[ig,:],0.),fo,0.)
            # The is bad because it would include non-detections:
            #w[ig,:]=where(greater(r[ig,:],0.),1.,0.)

    if save_probs:
        texto = '%s ' % str(id[ig])
        texto += len(p_bayes) * '%.3e ' + '\n'
        probs.write(texto % tuple(p_bayes))

    # pb[z,t] -> p_bayes[z]
    # 1. tb are summed over
    # 2. convolved with Gaussian if CONVOLVE_P
    # 3. Clipped above P_MIN * max(P), where P_MIN = 0.01 by default
    # 4. normalized such that sum(P(z)) = 1
    if save_probs2:  # P = exp(-chisq / 2)
        #probs2.write('%s\n' % id[ig])
        pmin = pmax * float(pars.d['P_MIN'])
        #pb = where(less(pb,pmin), 0, pb)
        chisq = -2 * log(pb)
        for itb in range(nt):
            chisqtb = chisq[:, itb]
            pqual = greater(pb[:, itb], pmin)
            chisqlists = seglist(chisqtb, pqual)
            if len(chisqlists) == 0:
                continue
            #print pb[:,itb]
            #print chisqlists
            zz = arange(zmin, zmax + dz, dz)
            zlists = seglist(zz, pqual)
            for i in range(len(zlists)):
                probs2.write('%s  %2d  %.3f  ' %
                             (id[ig], itb + 1, zlists[i][0]))
                fmt = len(chisqlists[i]) * '%4.2f ' + '\n'
                probs2.write(fmt % tuple(chisqlists[i]))
            #fmt = len(chisqtb) * '%4.2f '+'\n'
            #probs2.write('%d  ' % itb)
            #probs2.write(fmt % tuple(chisqtb))

            #if checkSED: open(pars.d['FLUX_COMPARISON'],'w').write(buffer_flux_comparison)
if checkSED: open(pars.d['CHECK'], 'w').write(buffer_flux_comparison)

if get_z: output.close()

#if checkSED and get_z:
if checkSED:
    #try:
    if 1:
        if interactive:
            print("")
            print("")
            print("PHOTOMETRIC CALIBRATION TESTS")
            # See PHOTOMETRIC CALIBRATION CHECK above
            #ratios=add.reduce(w*r,0)/add.reduce(w,0)
            #print "Average, weighted by flux ratios f_obs/f_model for objects with odds >= %g" % odd_check
            #print len(filters)*'  %s' % tuple(filters)
            #print  nf*' % 7.3f       ' % tuple(ratios)
            #print "Corresponding zero point shifts"
            #print  nf*' % 7.3f       ' % tuple(-flux2mag(ratios))
            #print

            fratavg = old_div(sum(fw * frat, axis=0), sum(fw, axis=0))
            dmavg = -flux2mag(fratavg)
            fnobj = sum(greater(fw, 0), axis=0)
            #print 'fratavg', fratavg
            #print 'dmavg', dmavg
            #print 'fnobj', fnobj
            #fnobj = sum(greater(w[:,i],0))
            print(
                "If the dmag are large, add them to the .columns file (zp_offset), then re-run BPZ.")
            print(
                "(For better results, first re-run with -ONLY_TYPE yes to fit SEDs to known spec-z.)")
            print()
            print('  fo/ft    dmag   nobj   filter')
            #print nf
            for i in range(nf):
                print('% 7.3f  % 7.3f %5d   %s'\
                    % (fratavg[i], dmavg[i], fnobj[i], filters[i]))
                #% (ratios[i], -flux2mag(ratios)[i], sum(greater(w[:,i],0)), filters[i])
                #print '  fo/ft    dmag    filter'
                #for i in range(nf):
                #    print '% 7.3f  % 7.3f   %s'  % (ratios[i], -flux2mag(ratios)[i], filters[i])
            print(
                "fo/ft = Average f_obs/f_model weighted by f_obs/ef_obs for objects with ODDS >= %g"
                % odd_check)
            print(
                "dmag = magnitude offset which should be applied (added) to the photometry (zp_offset)")
            print(
                "nobj = # of galaxies considered in that filter (detected and high ODDS >= %g)"
                % odd_check)
            # print r
            # print w
            #print
            #print "Number of galaxies considered (with ODDS >= %g):" % odd_check
            #print '  ', sum(greater(w,0)) / float(nf)
            #print '(Note a galaxy detected in only 5 / 6 filters counts as 5/6 = 0.833)'
            #print sum(greater(w,0))

            #This part is experimental and may not work in the general case
            #print "Median color offsets for objects with odds > "+`odd_check`+" (not weighted)"
            #print len(filters)*'  %s' % tuple(filters)
            #r=flux2mag(r)
            #print  nf*' %.3f       ' % tuple(-median(r))
            #print  nf*' %.3f       ' % tuple(median(dm))
            #rms=[]
            #efobs=[]

            #for j in range(nf):
            #    ee=where(greater(f_obs[:,j],0.),f_obs[:,j],2.)
            #    zz=e_frac2mag(ef_obs[:,j]/ee)
            #
            #    xer=arange(0.,1.,.02)
            #    hr=hist(abs(r[:,j]),xer)
            #    hee=hist(zz,xer)
            #    rms.append(std_log(compress(less_equal(r[:,j],1.),r[:,j])))
            #    zz=compress(less_equal(zz,1.),zz)
            #    efobs.append(sqrt(mean(zz*zz)))

            #print  nf*' %.3f       ' % tuple(rms)
            #print  nf*' %.3f       ' % tuple(efobs)
            #print  nf*' %.3f       ' % tuple(sqrt(abs(array(rms)**2-array(efobs)**2)))

            #except: pass

    if save_full_probs: full_probs.close()
    if save_probs: probs.close()
    if save_probs2: probs2.close()

if plots and checkSED:
    zb, zm, zb1, zb2, o, tb = get_data(out_name, (1, 6, 2, 3, 5, 4))
    #Plot the comparison between z_spec and z_B

    if 'Z_S' in col_pars.d:
        if not interactive or ask('Compare z_B vs z_spec?'):
            good = less(z_s, 9.99)
            print(
                'Total initial number of objects with spectroscopic redshifts= ',
                sum(good))
            od_th = 0.
            if ask('Select for galaxy characteristics?\n'):
                od_th = eval(input('Odds threshold?\n'))
                good *= greater_equal(o, od_th)
                t_min = eval(input('Minimum spectral type\n'))
                t_max = eval(input('Maximum spectral type\n'))
                good *= less_equal(tb, t_max) * greater_equal(tb, t_min)
                if has_mags:
                    mg_min = eval(input('Bright magnitude limit?\n'))
                    mg_max = eval(input('Faint magnitude limit?\n'))
                    good = good * less_equal(m_0, mg_max) * greater_equal(
                        m_0, mg_min)

            zmo, zso, zbo, zb1o, zb2o, tb = multicompress(good, (zm, z_s, zb,
                                                                 zb1, zb2, tb))
            print('Number of objects with odds > %.2f= %i ' %
                  (od_th, len(zbo)))
            deltaz = old_div((zso - zbo), (1. + zso))
            sz = stat_robust(deltaz, 3., 3)
            sz.run()
            outliers = greater_equal(abs(deltaz), 3. * sz.rms)
            print('Number of outliers [dz >%.2f*(1+z)]=%i' %
                  (3. * sz.rms, add.reduce(outliers)))
            catastrophic = greater_equal(deltaz * (1. + zso), 1.)
            n_catast = sum(catastrophic)
            print('Number of catastrophic outliers [dz >1]=', n_catast)
            print('Delta z/(1+z) = %.4f +- %.4f' % (sz.median, sz.rms))
            if interactive and plots:
                if plots == 'pylab':
                    figure(2)
                    subplot(211)
                    plot(
                        arange(
                            min(zso), max(zso) + 0.01, 0.01), arange(
                                min(zso), max(zso) + 0.01, 0.01), "r")
                    errorbar(zso,
                             zbo, [abs(zbo - zb1o), abs(zb2o - zbo)],
                             fmt="bo")
                    xlabel(r'$z_{spec}$')
                    ylabel(r'$z_{bpz}$')
                    subplot(212)
                    plot(zso, zmo, "go", zso, zso, "r")
                    xlabel(r'$z_{spec}$')
                    ylabel(r'$z_{ML}$')
                    show()
                elif plots == 'biggles':
                    plot = FramedPlot()
                    if len(zso) > 2000: symbol = 'dot'
                    else: symbol = 'circle'
                    plot.add(Points(zso, zbo, symboltype=symbol, color='blue'))
                    plot.add(Curve(zso, zso, linewidth=2., color='red'))
                    plot.add(ErrorBarsY(zso, zb1o, zb2o))
                    plot.xlabel = r'$z_{spec}$'
                    plot.ylabel = r'$z_{bpz}$'
                    #	    plot.xrange=0.,1.5
                    #	    plot.yrange=0.,1.5
                    plot.show()
                    #
                    plot_ml = FramedPlot()
                    if len(zso) > 2000: symbol = 'dot'
                    else: symbol = 'circle'
                    plot_ml.add(Points(
                        zso, zmo, symboltype=symbol,
                        color='blue'))
                    plot_ml.add(Curve(zso, zso, linewidth=2., color='red'))
                    plot_ml.xlabel = r"$z_{spec}$"
                    plot_ml.ylabel = r"$z_{ML}$"
                    plot_ml.show()

    if interactive and plots and ask('Plot Bayesian photo-z histogram?'):
        if plots == 'biggles':
            dz = eval(input('Redshift interval?\n'))
            od_th = eval(input('Odds threshold?\n'))
            good = greater_equal(o, od_th)
            if has_mags:
                mg_min = eval(input('Bright magnitude limit?\n'))
                mg_max = eval(input('Faint magnitude limit?\n'))
                good = good * less_equal(m_0, mg_max) * greater_equal(m_0,
                                                                      mg_min)
            z = compress(good, zb)
            xz = arange(zmin, zmax, dz)
            hz = hist(z, xz)
            plot = FramedPlot()
            h = Histogram(hz, 0., dz, color='blue')
            plot.add(h)
            plot.xlabel = r'$z_{bpz}$'
            plot.ylabel = r'$N(z_{bpz})$'
            plot.show()
            if ask('Want to save plot as eps file?'):
                file = eval(input('File name?\n'))
                if file[-2:] != 'ps': file = file + '.eps'
                plot.save_as_eps(file)

    if interactive and plots and ask(
            'Compare colors with photometric redshifts?'):
        if plots == 'biggles':
            color_m = zeros((nz, nt, nf - 1)) * 1.
            for i in range(nf - 1):
                plot = FramedPlot()
                # Check for overflows
                fmu = f_obs[:, i + 1]
                fml = f_obs[:, i]
                good = greater(fml, 1e-100) * greater(fmu, 1e-100)
                zz, fmu, fml = multicompress(good, (zb, fmu, fml))
                colour = old_div(fmu, fml)
                colour = clip(colour, 1e-5, 1e5)
                colour = 2.5 * log10(colour)
                d = Points(zz, colour, color='blue')
                plot.add(d)
                for it in range(nt):
                    #Prevent overflows
                    fmu = f_mod[:, it, i + 1]
                    fml = f_mod[:, it, i]
                    good = greater(fml, 1e-100)
                    zz, fmu, fml = multicompress(good, (z, fmu, fml))
                    colour = old_div(fmu, fml)
                    colour = clip(colour, 1e-5, 1e5)
                    colour = 2.5 * log10(colour)
                    d = Curve(zz, colour, color='red')
                    plot.add(d)
                plot.xlabel = r'$z$'
                plot.ylabel = '%s - %s' % (filters[i], filters[i + 1])
                plot.save_as_eps('%s-%s.eps' % (filters[i], filters[i + 1]))
                plot.show()

rolex.check()
