Entry 7 (First Place)

../../_images/entry7.png

Authors

  • Edward Taylor

Caption

Demonstrating the quality of our fits to the joint colour–mass distributions.— The histograms in each panel show the incompleteness-corrected distribution of colours of z < 0.12 galaxies, computed in bins of log M∗ centred on log M∗ = 8.8, 9.0, ..., 11.4. The left panel is for the effective (g-i) colour; the right panel is for the intrinsic, dust-corrected, stellar colour, (g-i)*. The errorbars show the statistical uncertainties on each these distributions, derived by bootstrap resampling. The smooth curves show the results of our modelling: the blue and red curves show the fit distributions for the B- and R-populations; the black curves shows the net B+R distributions. Underlaid beneath all this, the grey points show the data themselves; the size of each point reflects the relative 1/Vmax scaling used to account for selection effects. Note that we have not binned the data in the course of the fits; the binning in this Figure is for illustrative purposes only. It is clear that the fit model provides an excellent description of the observed data.

Description

These two panels are the showpiece of a new study that i have spent the last 3 years (!!) working on, looking into the idea of ‘bimodality’ among galaxies. I have just (finally!) submitted this work to MNRAS; the paper is actually the first in a series of 5 in which i am developing a phenomenological description of the main stages of galaxy evolution.

The observational result of ‘bimodality’ has played a crucial role in the development of our theoretical understanding of the process of galaxy formation and evolution. Let me explain how, as briefly as i can.

It seems that there are two and only two kinds of galaxies: red ones, and blue ones. (This is, of course, an over simplification – in fact, it is exactly this kind of oversimplification that i am trying to bust! But it is fair to say that this is how a lot of people in the field talk about things, even if they don’t really believe it.) In simple terms, a blue colour implies that a galaxy has formed at least some stars in the past few hundred million years. A red colour can mean one of two things: it can mean either a very old stellar population, with no star formation over the past billion years or so, or it can mean that there is a lot of dust in the galaxy. (Dust makes things red: think about a sunset on a smoggy evening.) So - modulo dust - the difference between red and blue is the difference between ‘star-forming’ and ‘dead’

The weird thing is that the most massive galaxies tend to be red and dead. This is weird because gravity is what drives star formation. Naively, you might think that more mass means more gravity, and so more star formation. But that’s not what happens: apparently, there is something that kills or ‘quenches’ the process of star formation in massive galaxies. Nobody knows what this might be, even though this quenching problem is one of the most intensively studied areas of extragalactic astronomy over the last ten or twelve years.

Models of galaxy evolution have to include some ad hoc prescription for when star formation is quenched, but pretty much every model has its own prescription – no one really understands what the real physical mechanism is. As ridiculous as it sounds, the best way to constrain these models is simply to measure the relative and absolute numbers of ‘red’ and ‘blue’ galaxies.

Even more ridiculously, even after ten or twelve years, no one has really measured this well. Let me explain what i mean by that.

What people have done in the part (with only one notable exception, in 2004) is simply to draw a line that divides ‘red’ galaxies from ‘blue’ ones. The problem with this is that if you draw a different line, you get a different answer. And the decision about precisely where to draw the line is always arguable, hence arbitrary. This is even worse than it sounds: different choices for the line give you qualitatively different measurements.

So what I have done is to develop a mixture modelling approach. Instead of drawing a line to separate ‘red’ and ‘blue’ galaxies, I model the data as a mixture of two distinct but overlapping populations. This lets me figure out the number of galaxies in the ‘red’ and ‘blue’ populations without ever specifying whether any particular galaxy belongs to which population. This lets me get around the problem of having to draw a line – in fact, it gives a way to derive an objective, phenomenological classification scheme that is based on the data.

These plots show just how good a description of the data the mixture modelling provides.

What’s more, in the ‘intrinsic’ plot, I have actually performed a correction for dust, so that I really am distinguishing galaxies based on their stars, not just their apparent colours. And even better than that, although this point isn’t visually emphasised in the plots I’ve submitted, i get the same answer in both cases. In other words, if you do the mixture modelling well, you get the right answer, no matter whether you do the dust correction or not – I have definitely gotten around the dust problem. And that’s pretty neat.

The only other thing that I’ll say is that I had to learn a lot of statistics in order to do all this. The fits shown in each panel have 35-40 parameters (having done the full and proper Bayesian model selection to make sure that these really are the best and simplest descriptions of the data). This is a level of statistical rigor that is very, very rare in the field. But it is, in my opinion the way forward in the new age of ‘big data’ astronomy.

That was a lot of background, but i hope it wasn’t too painful to read!

Instructions

There are two json archives that contain the data for each of the two panels. There is also the script bimix.py that should use these data to make the two plots. You can either (i hope!) just python bimix.py at the command line, or run bimix.makePlots() in interactive mode.

Source

#! /usr/bin/env python

import numpy as np
from matplotlib import pyplot as p

import json

import os, sys, random, time

bluecolor = ( 0, .185, 1, 1 )
redcolor =  ( 1, .185, 0, 1 )

masslimit = 8.7

def plotFitDistros( parVec, argList, figure=1, figsize=(4,9), figstring=None,
                    nbins=64, scaleby=0.765, massbin=0.2, massbins=None, 
                    showdata=True, showbad=True, scalepoints=True, 
                    applyerrors=True,
                    nowhite=False, bootstraps=200, printthings=False,
                    xlim=[8.675, 11.525], ylim=(-0.25, 1.45), 
                    xlabel='stellar mass, log $M_*$', 
                    ylabel='intrinsic stellar colour, $(g-i)_*$     ',
                    bestfit=False, logyblue=False, 
                    redcolor=redcolor, bluecolor=bluecolor,
                    ):

    # compute the best fit (or whatever) model from parVec
    ( xobs, yobs, xerr, yerr, xycovar, weights, redlike, bluelike, badlike,
      bluem, pblue, redm, pred, bluem1, pblue1, redm1, pred1, 
      bluem2, pblue2, redm2, pred2, bluey, bluedy, redy, reddy, result ) \
      = bivarModel( parVec, *argList, returnFullModel=True, plotprob=0 )
    goodlike = redlike + bluelike

    if yerr.shape != yobs.shape : applyerrors = False        

    # open the figure
    p.figure( figure, figsize=figsize ) ; p.clf()
    p.subplots_adjust( left=0.22, right=0.95, top=0.995, bottom=0.065 )
    
    # define point sizes, colors, etc.
    if scalepoints  :
        psize = 2. + 1./weights * 2.
    else :
        psize = 8
    pcolor= p.cm.gray_r( 0.2 )
    palpha= 1.

    bsize = 68. * np.clip( badlike / ( badlike + goodlike ), 0., 1. )**2.
 
    # plot the datapoints
    if showdata :
        if showbad :
            p.scatter( yobs, xobs, bsize, 'k', marker='x', zorder=-10 )
        p.scatter( yobs, xobs, psize, pcolor, alpha=palpha, 
                   edgecolors='none', zorder=-9 )

#    p.xlim( ylim ) ; p.draw()
    
    # define the mass bins in which histograms will be plotted
    if massbins == None :
        massbins = np.arange( masslimit, 11.5+massbin, massbin )[ ::-1 ]
#    print massbins

    # setup for showing the fit curves
    ytp = np.linspace( ylim[0], ylim[1], nbins*4 )

    densenorm = np.sum( 1 / weights )

    # now, for each bin ...
    for mi, m0 in enumerate( massbins ):
      # ... select the galaxies that fall in this bin
      inbin = np.where( np.abs( xobs - m0 - massbin/2. ) < massbin/2. ) 

      # histograms for the observed data
      hist, cedges = np.histogram(
          yobs[inbin], weights=1./weights[inbin], range=ylim, bins=nbins )
      pbins = ( cedges[:-1] + cedges[1:] ) / 2.
      cbinwidth = cedges[1]-cedges[0]
      hist = hist.astype( 'f' ) / densenorm / cbinwidth * scaleby
      # will be plotted later

      # bootstrapping to find statistical errors, if necessary
      histerr = np.zeros( (2, pbins.shape[0]) )
      if bootstraps <= 0 : 
          hihist = hist
          lohist = hist
      else :
          bootvalues = np.zeros( (pbins.shape[0], bootstraps) )
          binpop = len( inbin[0] )
#          print binpop
          for bi in xrange( bootstraps ):
              if bi == 0 :
                  choices = inbin
              else :
                  scores = ( np.random.rand( binpop ) * binpop ).astype( 'i' )
                  choices = inbin[0][ scores.ravel() ]
                  
              boothist, cedges = np.histogram(
                  yobs[choices], weights=1./weights[choices], 
                  range=ylim, bins=nbins )
              bootvalues[ :, bi ] = boothist

              nz = np.where( boothist != 0 )

          bootvalues = bootvalues.astype( 'f' ) / densenorm / cbinwidth * scaleby
          counter = np.linspace( 0., 1., bootstraps )
          percentiles = [ 0.158655, 0.841345 ]
          for pi in range( pbins.shape[0] ):
              sorted = np.argsort( bootvalues[ pi, : ] )
              histerr[ :, pi ] = np.interp( percentiles, counter, 
                                           bootvalues[ pi, sorted ] )
          hihist = histerr[ 1, : ]
          lohist = histerr[ 0, : ]

      if applyerrors :
          # figure out how to deal with the observational errors 
          redweight = ( redlike / yerr*2 / weights )[ inbin ]
          reddymean = np.sum( redweight * yerr[inbin] ) \
              / np.sum( redweight )

          blueweight = ( redlike / yerr*2 / weights )[ inbin ]
          bluedymean = np.sum( redweight * yerr[inbin] ) \
              / np.sum( redweight )

      # figure out how to show the fits themselves
      sumover = np.where( np.abs( redm - m0 - massbin/2. ) < massbin/2. )[0]
      rxtp = np.zeros( ytp.shape )
      bxtp = np.zeros( ytp.shape )
      for mj in sumover:
          rednorm= pred[ mj ]
          redcen = redy[ mj ]
          redvar = reddy[mj ]**2 
          if applyerrors : redvar += reddymean**2.

          rxtp += rednorm / np.sqrt( 2. * np.pi * redvar ) \
              * np.exp( -0.5*(ytp - redcen)**2 / redvar ) 

          bluenorm= pblue[ mj ]
          bluecen = bluey[ mj ]
          bluevar = bluedy[mj ]**2.
          if applyerrors : bluevar += bluedymean**2.
              
          bxtp +=bluenorm / np.sqrt( 2. * np.pi *bluevar ) \
                  * np.exp( -0.5*(ytp - bluecen)**2 /bluevar ) 

#          print redm[ mj ]
#          print 'red  :', np.sum(rxtp)/np.sum(rxtp+bxtp), redcen,np.sqrt( redvar)
#          print 'blue :', np.sum(bxtp)/np.sum(rxtp+bxtp),bluecen,np.sqrt(bluevar)

      rxtp *= massbin * scaleby / sumover.shape[0]
      bxtp *= massbin * scaleby / sumover.shape[0]

#      print 'red  :', np.sum(rxtp)/np.sum(rxtp+bxtp), redcen,np.sqrt( redvar)
#      print 'blue :', np.sum(bxtp)/np.sum(rxtp+bxtp),bluecen,np.sqrt(bluevar)

      if not nowhite :
          # fill in white under curves for this bin
          p.fill_between( ytp, m0+0*ytp, m0+rxtp+bxtp, color='w', zorder=mi-0.8 )
          p.fill_between( pbins, m0+0*hihist, m0+hihist, color='w', zorder=mi-0.8 )

          # replot data points as necessary
          limit1 = np.interp( yobs, ytp, m0+rxtp+bxtp )
          limit2 = np.interp( yobs, pbins, m0+hihist )
          limit = np.where( limit1 > limit2, limit1, limit2 )
          replot = np.where( ( xobs > m0 ) & ( xobs < limit ) )
          if len( replot[0] ) and showdata :
              p.scatter( yobs[replot], xobs[replot], psize[replot], pcolor, 
                         alpha=palpha, edgecolors='none', zorder=mi-0.7 )

      # fill in under the fit curves
      for ( xtp, ctp ) in zip( ( rxtp, bxtp ),  ( redcolor, bluecolor ) ):
          sig = np.where( xtp > 1e-4 )
          p.fill_between( ytp[ sig ], m0 + 0.*xtp[ sig ],
                          m0+xtp[ sig ], color=ctp, alpha=0.2, zorder=mi-0.6  )
                  
      # plot the statistical uncertainties from bootstrapping
      if bootstraps >= 0 :
          for pi in np.where( hist > 0 )[0] :
              p.plot( pbins[pi] * np.ones(2), [ m0+lohist[pi], m0+hihist[pi] ], 
                      'w-', lw=2.4, zorder=mi-0.5 )
              p.plot( pbins[pi] * np.ones(2), [ m0+lohist[pi], m0+hihist[pi] ], 
                      'k-', lw=0.8, zorder=mi )

      # plot the histograms for the data
      p.plot( pbins, m0+hist, 'w-', lw=4, drawstyle='steps-mid', zorder=mi-0.5 )
      p.plot( pbins, m0+hist, 'k-',lw=1.65,drawstyle='steps-mid',zorder=mi )

      # actually show the fits
      for ( xtp, ctp, lw ) in zip( ( rxtp + bxtp, bxtp, rxtp ), 
                                   ( 'k', bluecolor, redcolor ), 
                                   ( 1.5, 1, 1 ) ) :
          sig = np.where( xtp > 1e-4 )
          p.plot( ytp[ sig ], m0+xtp[ sig ], '-', color='w', lw=lw+1, zorder=mi )
          p.plot( ytp[ sig ], m0+xtp[ sig ], '-', color=ctp, lw=lw, zorder=mi+.1 )
                  
      # draw out the bin limits
      p.plot( ylim, np.ones(2) * m0, '-', 
              color=p.cm.gray( 0.7 ), lw=lw, zorder=mi+.3 )

    # limits, labels, and draw the fucker
    p.yticks( np.arange( 8.7, 11.6, 0.2 ) ) ; p.ylim( xlim ) 

    
    ytv = np.arange( -0.4, 2.0, 0.1 )
    ytn = ytv.tolist()
    for i in np.arange( 1, ytv.shape[0], 1 ):
        if ( (ytv[i] % 0.2 ) > .001 ) and ( (ytv[i] % 0.2 ) < .199 ) :
            ytn[ i ] = ' '
        else :
            ytn[ i ] = '%.1f' % ytv[i]
    p.xticks( ytv, ytn ) ; p.xlim( ylim )
        
    p.ylabel( xlabel, fontsize='xx-large' )
    p.xlabel( ylabel, fontsize='xx-large' )

    if figstring != None :
        figstr = '%s - %s' % ( figstring, time.ctime()[ 4:16 ] )
        p.figtext( 0.99, 0.995, figstr, fontsize='small', color='gray', 
                   va='top', ha='right', rotation='vertical' )

    p.draw()


# _____________________________________________________________________________
#
# THE BASIC DEFINITION OF THE MODEL
# _____________________________________________________________________________

def bivarModel( parVec, parList, 
                xobs, xunc, yobs, yunc, xycovar, weights,
                bluePrior, redPrior, badPrior, isinsample, 
                blueyfunc, bluedyfunc, redyfunc, reddyfunc, 
                mbinsize=0.08, 
                logyblue=False, redNormalMF=False, usePriors=False, 
                plotprob=1e-3 + 0.0000, printprob=.02 + 1e-3, 
                returnFullModel=False, computelogp=True, showdata=True ) :

        plotscore = random.random()
        plotme = plotscore < plotprob
        printscore = random.random()
        printme = printscore < printprob

        if plotme : printme = True
        if printme :
            startTime = time.time()

        par = unpackPars( parVec, parList )

        mmod = np.arange( masslimit, 11.75, mbinsize ) + mbinsize/2.

        # enforce priors as needed
        if ( (par['redfrac'] < 0) 
             | (par['redfrac'] < 0) | (par['redfrac'] > 1 )
             | (par['red2frac'] < 0) | (par['red2frac'] > 1 )
             | (par['blue2frac'] < 0) | (par['blue2frac'] > 1 )
             | (par['badrate'] < 0) | (par['badrate'] > 1)
             | (par['yerrscale'] < 0) 
             | (par['xscatterbad'] < 0) | (par['yscatterbad'] < 0) 
             | (par['blueq0']-par['blueq2'] < 9.)
             | (par['redq0']-par['redq2'] < 9.)
             | (par['blueq0']+par['blueq2'] > 11.5 )
             | (par['redq0']+par['redq2'] > 11.5 )
#             | (par['bluet0']-par['bluet2'] < 9.)
#             | (par['redt0']-par['redt2'] < 9.)
#             | (par['bluet0']+par['bluet2'] > 11.5 )
#             | (par['redt0']+par['redt2'] > 11.5 )
             | (par['blueq2'] < 0) | (par['redq2'] < 0) 
             | (par['blueq2'] > 3) | (par['redq2'] > 3) 
#             | (par['bluet2'] < 0) | (par['redt2'] < 0) 
#             | (par['bluet2'] > 3) | (par['redt2'] > 3) 
             | (abs( par['bluelogm0'] - 10 ) > 2 )
             | (abs( par['bluelogm02'] - 10 ) > 2 )
             | (abs( par['redlogm0'] - 10 ) > 2 )
             | (abs( par['redlogm02'] - 10 ) > 2 )
             ) and not returnFullModel :
#            print '_',
            return -np.inf

        # compute blue mass function at discrete mmods
        bluem1, pblue1 \
            = SchechterFunction( mmod, par['bluelogm0'], par['bluealpha'] )
        if par['blue2frac'] > 0 :
            pblue1 *= ( 1. - par['blue2frac'] ) * ( 1. - par[ 'redfrac' ] )
            bluem2, pblue2 \
                = SchechterFunction( mmod,par['bluelogm02'],par['bluealpha2'] )
            pblue2 *= par['blue2frac' ] * ( 1. - par[ 'redfrac' ] )
            bluem = ( pblue1*bluem1 + pblue2*bluem2 ) / ( pblue1 + pblue2 )
            pblue = pblue1 + pblue2
        else :
            bluem2, pblue2 = 0, 0
            pblue1 *= ( 1. - par[ 'redfrac' ] )
            bluem, pblue = bluem1, pblue1

        # compute red mass function at discrete mmods
        redm1, pred1 \
            = SchechterFunction( mmod, par['redlogm0'], par['redalpha'] )
        if par['red2frac'] > 0 :
            pred1 *= ( 1. - par['red2frac'] ) * par[ 'redfrac' ]
            redm2, pred2 \
                = SchechterFunction( mmod, par['redlogm02'], par['redalpha2'] )
            pred2 *= par['red2frac' ] * par[ 'redfrac' ]
            redm = ( pred1*redm1 + pred2*redm2 ) / ( pred1 + pred2 )
            pred = pred1 + pred2
        else :
            redm2, pred2 = 0, 0.
            pred1 *= par[ 'redfrac' ]
            redm, pred = redm1, pred1

#        print np.sum( pred1 )

        if redNormalMF :
            redm1 = mmod
            pred1 = par['redfrac']/np.sqrt( 2. * np.pi * par['redalpha']**2. ) \
                * np.exp( -0.5 * ( redm1 - par['redlogm0'] )**2 
                          / par['redalpha']**2. )
            redm2 = mmod
            pred2 = 0 * pred1
            pred = pred1
#        print np.sum( pred1 )

        bluey = GeneralLines( 
            bluem, par['bluep0' ], par['bluep1' ],
            par['blueq0' ], par['blueq1' ], par['blueq2' ], par['blueq3' ] )


        if bluedyfunc != '' and bluedyfunc != None :
            bluedy= bluedyfunc( ( par['bluet3'], par['bluet2'], par['bluet1'],
                                 par['bluet0'], par['blues1'], par['blues0'],
                                 ), bluem - 10.3 )
        else :
            bluedy= GeneralLines(
                bluem, par['blues0' ], par['blues1' ], par['bluet0' ], 
                par['bluet1' ], par['bluet2' ], par['bluet3' ] )

        redy  = GeneralLines(
            redm , par[ 'redp0' ], par[ 'redp1' ],
            par[ 'redq0' ], par[ 'redq1' ], par[ 'redq2' ], par['redq3' ] ) 

        if reddyfunc != '' and reddyfunc != None :
            reddy= reddyfunc( ( par['redt3'], par['redt2'], par['redt1'],
                                par['redt0'], par['reds1'], par['reds0'],
                                ), redm - 10.3 )
        else :
            reddy = GeneralLines( 
                redm , par[ 'reds0' ], par[ 'reds1' ], par[ 'redt0' ], 
                par[ 'redt1' ], par[ 'redt2' ], par['redt3' ] ) 

#        bluedy = bluey - bluedy
#        reddy = reddy - redy

#        if np.any( bluedy < 0 ) or np.any( reddy < 0 ):
#            return -np.inf

        xvar = xunc**2.

        # scale y uncertainties as required
        yerr = np.sqrt( ( par['yerrscale']*yunc )**2.
                        + ( par['yerradjust'] )**2 )

        rhoxy = xycovar * xunc * yerr

        if returnFullModel and not computelogp :
            return ( xobs, yobs, xunc, yerr, xycovar, weights,
                     0., 0., 0., 
                     bluem, pblue, redm, pred, 
                     bluem1, pblue1, redm1, pred1, 
                     bluem2, pblue2, redm2, pred2, 
                     bluey, bluedy, redy, reddy, 
                     0. )

        if ( ( reddy.min() < 0 ) or ( bluedy.min() < 0 ) ) and not returnFullModel:
#            print '|',
            return -np.inf

        redlike, bluelike, badlike, like = 0., 0., 0., 0.
        for i, (mr, mb, ry, by, rs, bs, pr, pb) in enumerate( 
            zip( redm, bluem, redy, bluey, reddy, bluedy, pred, pblue ) ):

            xoff = ( xobs - mr )
            yoff = ( yobs - ry )

            yvar = yerr**2 + rs**2

            red = pr * ( 1. - bluePrior ) * ( 1. - badPrior ) \
                * Gaussian( xoff, yoff, xvar, yvar, rhoxy )

            badred = pr * ( 1. - bluePrior ) * badPrior \
                * Gaussian( xoff, yoff, xvar + par['xscatterbad']**2, 
                                yvar + par['yscatterbad']**2, rhoxy )

            xoff = ( xobs - mb )
            if not logyblue :
                yoff = ( yobs - by )
                yvar = yerr**2 + bs**2
            else :
                yoff = ( np.log10( yobs ) - by )
                yvar = (yerr/yobs/np.log(10.))**2. + bs**2
#                    yvar = bs**2

            blue = pb * ( 1. - redPrior ) * ( 1. - badPrior ) \
                * Gaussian( xoff, yoff, xvar, yvar, rhoxy )

            badblue = pb * ( 1. - redPrior ) * badPrior \
                * Gaussian( xoff, yoff, xvar + par['xscatterbad']**2, 
                            yvar + par['yscatterbad']**2, rhoxy )

            if logyblue :
                blue[ np.where( yobs <= 0 ) ] = 0.
                badblue[ np.where( yobs <= 0 ) ] = 0.

            if plotme or returnFullModel :
                redlike += ( 1. - par[ 'badrate' ] ) * red \
                    + ( par[ 'badrate' ] ) * badred
                bluelike += ( 1. - par[ 'badrate' ] ) * blue \
                    + ( par[ 'badrate' ] ) * badblue
                badlike += ( par[ 'badrate' ] ) * ( badred + badblue  )

            like += ( 1. - par['badrate'] ) * ( red + blue ) \
                + par['badrate'] * ( badred + badblue )

        use = np.where( isinsample == 1 )
        finalval = isinsample / weights * ( np.log( like ) - 1.837377 )
        result = np.sum( finalval[use] )
#        print len( use[0] ), result

        if plotme or printprob >= 1 :
            headstr = ''
            for param in parList :
                headstr += param[ -5: ]+' '
            print  headstr

        if plotme :
            plotModel( xobs, yobs, xvar, yvar, xycovar, weights,
                       redlike, bluelike, badlike,
                       bluem, pblue, redm, pred, 
                       bluem1, pblue1, redm1, pred1, 
                       bluem2, pblue2, redm2, pred2, 
                       bluey, bluedy, redy, reddy, 
                       result, logyblue=logyblue, showdata=showdata )

        if printme :
            printstr = ''
            for param in parList :
                printstr += ('%5.2f ' % par[ param ])
            print printstr, ('%.1f' % result),
            print ( ' (%.3fs)' ) % ( ( time.time() - startTime ) )
            sys.stdout.flush()

        if returnFullModel :
            return ( xobs, yobs, xunc, yerr, xycovar, weights,
                     redlike, bluelike, badlike,
                     bluem, pblue, redm, pred, 
                     bluem1, pblue1, redm1, pred1, 
                     bluem2, pblue2, redm2, pred2, 
                     bluey, bluedy, redy, reddy, 
                     result )

        if np.isfinite( result ) :
            return result
        else :
            return -np.inf

# _____________________________________________________________________________

def GeneralLines( x, intercept, theta, transition, delta, softening, 
                  dnorm=0., error=False ):

    line = np.tan( theta ) * ( x - transition ) + intercept

    deflect = np.tanh( ( x - transition ) / softening ) \
              * ( np.tan( delta ) * ( x - transition ) + dnorm )

    value = line + deflect

    if error :
        value = np.clip( value, 1e-8, value )

    return value

# _____________________________________________________________________________

def SchechterFunction( logm, logm0, alpha, 
                       cutoff=masslimit, justvalues=False, 
                       logmstep=0.005 ):

    logmmin, logmmax = masslimit-0.5, 12.8
    distlogm = np.arange( logmmax, logmmin, -logmstep )

    massfunc = ( 10**( distlogm - logm0 ) )**( alpha+1. ) \
        * np.exp( -10**( distlogm - logm0 ) )

    cumfunc = np.cumsum( massfunc ) * logmstep
    norm = np.interp( [-masslimit], -distlogm, cumfunc )[0]

    cumfunc /= norm
    massfunc /= norm

    if justvalues :
        return np.interp( -logm, -distlogm, massfunc )

    edges = np.zeros( logm.shape[0] + 1 )
    edges[ 0 ] = logm[0]
    edges[ 1:-1 ] = ( logm[ :-1 ] + logm[ 1: ] ) / 2.
    edges[ -1 ] = logm[-1]

    nivalues = - np.interp( -edges, -distlogm, -cumfunc )
    fvalues = - np.diff( nivalues ) / np.diff( edges )
    # now integrated number of galaxies per bin

    nmivalues = - np.interp( -edges, -distlogm, 
                                -np.cumsum( massfunc*10**distlogm ) * logmstep )
    newvalues = np.log10( np.diff( nmivalues ) / np.diff( nivalues ) )

#        nogood = np.where( fvalues <= 0 )
#        newvalues[ nogood ] = logm[ nogood ]

    return newvalues, fvalues

# _____________________________________________________________________________

def Gaussian( x, y, xvar, yvar, rhoxy=0 ):
        # NB. rhoxy must be understood to be the off-diagonal element of the
        # error/covar matrix, NOT the Peason correlation coeff.

        norm = xvar * yvar - rhoxy**2

#        print xvar[ ::500 ]
#        print yvar[ ::500 ]
#        print rhoxy[ ::500 ]
#        print norm[ ::500 ]

        value = 1. / np.sqrt( norm ) \
            * np.exp( -0.5 * ( x**2 * yvar - 2. * x * y * rhoxy + y**2 * xvar )
                     / norm )

#        check = np.where( value <= 0 )[0] 
#        print len( check[0] ), x[check], y[check]

        return value


def unpackPars( parVec=[], parList=[], fittingfor='gminusi_stars' ):

# specify default parameters here
    parDict = { 'bluelogm0' : 10.5,
                'bluealpha' : -1.0,
                'redfrac'   : 0.2,
                'redlogm0'  : 10.9,
                'redalpha'  : 0.0,
                'blue2frac' : 0.00,
                'bluelogm02': 10.5,
                'bluealpha2': -1.2,
                'red2frac'  : 0.,
                'redlogm02' : 10.9,
                'redalpha2' : 0.0,
                'bluep0'    : 0.40,
                'bluep1'    : 0.00,
                'blueq0'    : 10.2,
                'blueq1'    : 0.00,
                'blueq2'    : 0.10,
                'blueq3'    : 0.00,
                'redp0'     : 1.00,
                'redp1'     : 0.00,
                'redq0'     : 10.2,
                'redq1'     : 0.00,
                'redq2'     : 0.10,
                'redq3'     : 0.00,
                'blues0'    : 0.15,
                'blues1'    : 0.00,
                'bluet0'    : 10.2,
                'bluet1'    : 0.00,
                'bluet2'    : 0.10,
                'bluet3'    : 0.00,
                'reds0'     : 0.15,
                'reds1'     : 0.00,
                'redt0'     : 10.2,
                'redt1'     : 0.00,
                'redt2'     : 0.10,
                'redt3'     : 0.00,
                'badrate'   : 0.00,
                'xscatterbad': 0.0,
                'yscatterbad': 0.5,
                'yerrscale' : 1.00,
                'yerradjust': 0.00 }

    if fittingfor == 'EWHalpha' :
        parDict[ 'bluep0' ] = 1.5
        parDict[ 'redp0' ] = -1.0
        parDict[ 'reds0' ] = 0.5
        parDict[ 'blues0' ] = 0.5
        parDict[ 'blueq0' ] = 9.5

    for parName, parValue in zip( parList, parVec ):
        parDict[ parName ] = parValue

    return parDict


def makePlots():

    savefile = json.loads( open( 'bimix_intrinsic.json', 'r' ).read() )
    parVec  = np.array( savefile[ u'parVec' ] )
    argList = savefile[ u'argList' ]
    for i in [ 1, 2, 3, 4, 5, 6, -7 ]:
        argList[ i ] = np.array( argList[ i ] )

    plotFitDistros( parVec, argList, ylim=( -0.25, 1.35 ), 
        ylabel='intrinsic stellar colour, $(g-i)_*$     ' )

    p.savefig( 'bimix_intrinsic.pdf' )


    savefile = json.loads( open( 'bimix_effective.json', 'r' ).read() )
    parVec  = np.array( savefile[ u'parVec' ] )
    argList = savefile[ u'argList' ]
    for i in [ 1, 2, 3, 4, 5, 6, -7 ]:
        argList[ i ] = np.array( argList[ i ] )

    plotFitDistros( parVec, argList, ylim=( -0.05, 1.55 ), 
        ylabel='effective stellar colour, $(g-i)_*$     '  )

    p.savefig( 'bimix_effective.pdf' )

if __name__ == '__main__' :

    makePlots()