Resampling mass histograms

gistfile1.py

'''
Created on 19/05/2014

@author: andre
'''

import numpy as np 
import matplotlib.pyplot as plt

###################################################################################################

def bin_edges(bin_center):
    '''
    Computes the bin edges as the bissection of the bins,
    expanding the borders accordingly.

        * Units: :math:`[\log age / yr]`
        * Shape: ``(len(bin_center) + 1)``
    '''
    bin_edges = np.empty(len(bin_center) + 1)
    bin_edges[1:-1] = (bin_center[:-1] + bin_center[1:]) / 2.0
    bin_edges[0] = bin_center[0] - (bin_center[1] - bin_center[0]) / 2.0
    bin_edges[-1] = bin_center[-1] + (bin_center[-1] - bin_center[-2]) / 2.0
    return bin_edges

###################################################################################################

def hist_resample(bins_o, bins_r, v):
    n_v_r = len(bins_r) - 1
    v_r = np.zeros(n_v_r, dtype=v.dtype)
    i_0 = 0
    j = 0
    lbo = bins_o[i_0]
    rbo = bins_o[i_0+1]
    lbr = bins_r[j]
    rbr = bins_r[j+1]
    
    if lbo < lbr:
        # Skip leading original bins that do not overlap the resampled bins.
        while lbr > rbo:
            i_0 += 1
            lbo = bins_o[i_0]
            rbo = bins_o[i_0+1]
            dbo = rbo - lbo
        frac = (rbo - lbr) / dbo
        v_r[j] += frac * v[i_0]
        i_0 += 1
        last_edge = rbo
    else:
        # Skip leading resampled bins that do not overlap the original bins.
        while rbr < lbo:
            j += 1
            lbr = bins_r[j]
            rbr = bins_r[j+1]
        last_edge = rbr
        
    for i in xrange(i_0, len(v)):
        rbo = bins_o[i+1]
        if rbo < rbr:
            # This original bin is entirely contained in the resampled bin.
            v_r[j] += v[i]
            last_edge = rbo
            continue
        lbo = bins_o[i]
        dbo = rbo - lbo
        while rbr < rbo:
            # Compute the slices of the original bin going to the resampled bins.
            frac = (rbr - last_edge) / dbo
            v_r[j] += frac * v[i]
            last_edge = rbr
            j += 1
            if j >= n_v_r:
                break
            rbr = bins_r[j+1]
        if j >= n_v_r:
            break
        # Compute the last slice of the original bin.
        frac = (rbo - last_edge) / dbo
        v_r[j] += frac * v[i]
        last_edge = rbo
    return v_r

###################################################################################################

# Initial bins spaced logarithmically.
logtc = np.arange(7, 11, 0.5)
logtc_bins = bin_edges(logtc)
tc_bins = 10**logtc_bins
dtc = tc_bins[1:] - tc_bins[:-1]

# Linear resampled bins.
dt = 0.5e9
tl = np.arange(tc_bins.min(), tc_bins.max()+dt, dt)
tl_bins = bin_edges(tl)
if tl_bins[0] < 0.0:
    tl_bins[0] = 0.0
dtl = tl_bins[1:] - tl_bins[:-1]

# Mass is proportional to bin width in linear time (equal to a constant SFR)
M = 10 * dtc

# This function treats M as an histogram.
M_r = hist_resample(tc_bins, tl_bins, M)

# Try not to create/destroy mass.
print 'M = %f, M_r = %f, M/M_r = %f' % (M.sum(), M_r.sum(), M.sum() / M_r.sum())

# plot the rates (M / dt), should be the same in both cases.
plt.ioff()
plt.bar(tc_bins[:-1], M/dtc, dtc, color='r')
plt.bar(tl_bins[:-1], M_r/dtl, dtl, color='b')
plt.show()