def generate_smooth_gp_re_a(out_fname='data.csv', country_variation=True):
    """ Generate random data based on a nested gaussian process random
    effects model with age, with covariates that vary smoothly over
    time (where unexplained variation in time does not interact with
    unexplained variation in age)

    This function generates data for all countries in all regions, and
    all age groups based on the model::

        Y_r,c,t = beta * X_r,c,t + f_r(t) + g_r(a) + f_c(t)

        beta = [30., -.5, .1, .1, -.1, 0., 0., 0., 0., 0.]
        f_r ~ GP(0, C(3.))
        g_r ~ GP(0, C(2.))
        f_c ~ GP(0, C(1.)) or 0 depending on country_variation flag
        C(amp) = Matern(amp, scale=20., diff_degree=2)

        X_r,c,t[0] = 1
        X_r,c,t[1] = t - 1990.
        X_r,c,t[k] ~ GP(t; 0, C(1)) for k >= 2
    """
    c4 = countries_by_region()

    data = col_names()

    beta = [30., -.5, .1, .1, -.1, 0., 0., 0., 0., 0.]
    C0 = gp.matern.euclidean(time_range, time_range, amp=1., scale=25., diff_degree=2)
    C1 = gp.matern.euclidean(age_range, age_range, amp=1., scale=25., diff_degree=2)
    C2 = gp.matern.euclidean(time_range, time_range, amp=.1, scale=25., diff_degree=2)
    C3 = gp.matern.euclidean(time_range, time_range, amp=1., scale=25., diff_degree=2)

    g = mc.rmv_normal_cov(pl.zeros_like(age_range), C1)
    for r in c4:
        f_r = mc.rmv_normal_cov(pl.zeros_like(time_range), C0)
        g_r = mc.rmv_normal_cov(g, C1)
        for c in c4[r]:
            f_c = mc.rmv_normal_cov(pl.zeros_like(time_range), C2)

            x_gp = {}
            for k in range(2,10):
                x_gp[k] = mc.rmv_normal_cov(pl.zeros_like(time_range), C3)

            for j, t in enumerate(time_range):
                for i, a in enumerate(age_range):
                    x = [1] + [j] + [x_gp[k][j] for k in range(2,10)]
                    y = float(pl.dot(beta, x)) + f_r[j] + g_r[i]
                    if country_variation:
                        y += f_c[j]
                    se = 0.
                    data.append([r, c, t, a, y, se] + list(x))
    write(data, out_fname)

def show_grey_channels(I):
    K = average(I, axis=2)
    for i in range(3):
        J = zeros_like(I)
        J[:, :, i] = K
        figure(i+10)
        imshow(J)

def plot_count(fname, dpi=70):
    # Load data
    date, libxc, c, code, test, doc = np.loadtxt(fname, unpack=True)
    zero = pl.zeros_like(date)

    fig = pl.figure(1, figsize=(10, 5), dpi=dpi)
    ax = fig.add_subplot(111)
    polygon(date, c + libxc + code + test, c + libxc + code + test + doc, facecolor="m", label="Documentation")
    polygon(date, c + libxc + code, c + libxc + code + test, facecolor="y", label="Tests")
    polygon(date, c + libxc, c + libxc + code, facecolor="g", label="Python-code")
    polygon(date, c, c + libxc, facecolor="c", label="LibXC-code")
    polygon(date, zero, c, facecolor="r", label="C-code")
    polygon(date, zero, zero, facecolor="b", label="Fortran-code")

    months = pl.MonthLocator()
    months4 = pl.MonthLocator(interval=4)
    month_year_fmt = pl.DateFormatter("%b '%y")

    ax.xaxis.set_major_locator(months4)
    ax.xaxis.set_minor_locator(months)
    ax.xaxis.set_major_formatter(month_year_fmt)
    labels = ax.get_xticklabels()
    pl.setp(labels, rotation=30)
    pl.axis("tight")
    pl.legend(loc="upper left")
    pl.title("Number of lines")
    pl.savefig(fname.split(".")[0] + ".png", dpi=dpi)

 def plot(self, plot_range):
     u = pb.linspace(-plot_range, plot_range, 1000)
     z = pb.zeros_like(u)
     for i, j in enumerate(u):
         z[i] = self(j)
     pb.plot(u, z)
     pb.show()

def bad_model(X):
    """ Results in a matrix with shape matching X, but all rows sum to 1"""
    N, T, J = X.shape
    Y = pl.zeros_like(X)
    for t in range(T):
        Y[:,t,:] = X[:,t,:] / pl.outer(pl.array(X[:,t,:]).sum(axis=1), pl.ones(J))
    return Y.view(pl.recarray) 

def dS_dX(x0, PR, h_mag = .0005):
    """
    calculates the Jacobian of the SLIP at the given point x0,
    with PR beeing the parameters for that step
    coordinates under consideration are:
        y
        vx
        vz
    only for a single step!
    """
    df = []
    for dim in range(len(x0)):
        delta = zeros_like(x0)
        delta[dim] = 1.            
        h = h_mag * delta      
        # in positive direction           
        resRp = sl.SLIP_step3D(x0 + h, PR)
        SRp = array([resRp['y'][-1], resRp['vx'][-1], resRp['vz'][-1]])
        #fhp = array(SR2 - x0)
        # in negative direction
        resRn = sl.SLIP_step3D(x0 - h, PR)
        SRn = array([resRn['y'][-1], resRn['vx'][-1], resRn['vz'][-1]])
        #fhn = array(SR2 - x0)
        # derivative: difference quotient
        df.append( (SRp - SRn)/(2.*h_mag) )
    
    return vstack(df).T

def dist2(x):
  R, GSIGMAS = pylab.meshgrid(r[r<fit_rcutoff], gsigmas)
  g = pylab.zeros_like(GSIGMAS)
  g = evalg(x, GSIGMAS, R)
  gfit = pylab.reshape(g, len(eta)*len(r[r<fit_rcutoff]))
  return gfit - pylab.reshape([g[r<fit_rcutoff] for g in ghs],
                              len(gsigmas)*len(r[r<fit_rcutoff]))

def dist2(x):
  R, ETA = pylab.meshgrid(r[r<fit_rcutoff], eta)
  g = pylab.zeros_like(ETA)
  g = evalg(x, ETA, R)
  gfit = pylab.reshape(g, len(eta)*len(r[r<fit_rcutoff]))
  return gfit - pylab.reshape([g[r<fit_rcutoff] for g in ghs],
                              len(eta)*len(r[r<fit_rcutoff]))

def plot_count(fname, dpi=70):
    # Load data
    date, libxc, c, code, test, doc = np.loadtxt(fname, unpack=True)
    zero = pl.zeros_like(date)

    fig = pl.figure(1, figsize=(10, 5), dpi=dpi)
    ax = fig.add_subplot(111)
    polygon(date, c + code + test, c + code + test + doc,
             facecolor='m', label='Documentation')
    polygon(date, c + code, c + code + test,
             facecolor='y', label='Tests')
    polygon(date, c, c + code,
            facecolor='g', label='Python-code')
    polygon(date, zero, c,
            facecolor='r', label='C-code')
    polygon(date, zero, zero,
            facecolor='b', label='Fortran-code')

    months = pl.MonthLocator()
    months4 = pl.MonthLocator(interval=4)
    month_year_fmt = pl.DateFormatter("%b '%y")

    ax.xaxis.set_major_locator(months4)
    ax.xaxis.set_minor_locator(months)
    ax.xaxis.set_major_formatter(month_year_fmt)
    labels = ax.get_xticklabels()
    pl.setp(labels, rotation=30)
    pl.axis('tight')
    pl.legend(loc='upper left')
    pl.title('Number of lines')
    pl.savefig(fname.split('.')[0] + '.png', dpi=dpi)

def edge_props2(sheet, cut_node=None, trail=None, dead_end=None):
    counts_file = '%s/reformated_counts%s.csv' % (DATASETS_DIR, sheet)
    df = pd.read_csv(counts_file, names = ['source', 'dest', 'time'], skipinitialspace=True)
    df['time'] = pd.to_datetime(df['time'])
    df.sort_values(by='time', inplace=True)
    
    times = list(df['time'])
    deltas = []
    starttime = times[0]
    for time in times:
        deltas.append((time - starttime) / pylab.timedelta64(1, 's'))

    sources = list(df['source'])
    dests = list(df['dest'])
    counts = {}
    
    delta_edges = defaultdict(list)
    for i in xrange(len(deltas)):
        delta = deltas[i]
        source = sources[i]
        dest = dests[i]
        #edge = (source, dest)
        edge = tuple(sorted((source, dest)))
        if cut_node == None or cut_node in edge:
            delta_edges[delta].append(edge)
            counts[edge] = []

    for delta in sorted(delta_edges.keys()):
        for edge in counts:
            count = 0
            if len(counts[edge]) > 0:
                count = counts[edge][- 1]
            if edge in delta_edges[delta]:
                count += 1
            counts[edge].append(count)

    step_times = pylab.arange(1, len(delta_edges.keys()) + 1, dtype=pylab.float64)
    norms = pylab.zeros_like(step_times)
    for edge in counts:
        counts[edge] /= step_times
        norms += counts[edge]
 
    pylab.figure()
    for edge in counts:
        counts[edge] /= norms
        if (trail == None and dead_end == None) or ((edge == trail) or (edge == dead_end)):
            label = edge
            if trail != None and edge == trail:
                label = 'trail'
            elif dead_end != None and edge == dead_end:
                label = 'dead end'
            pylab.plot(sorted(delta_edges.keys()), counts[edge], label=label)

    pylab.legend()
    pylab.xlabel('time (seconds)')
    pylab.ylabel('proportion of choices on edge')
    pylab.savefig('cut_edge_props%s.pdf' % sheet, format='pdf')
    pylab.close()

def int_f(a, fs=1.):
    """
    A fourier-based integrator.

    ===========
    Parameters:
    ===========
    a : *array* (1D)
        The array which should be integrated
    fs : *float*
        sampling time of the data

    ========
    Returns:
    ========
    y : *array* (1D)
        The integrated array

    """

    if False:
    # version with "mirrored" code
        xp = hstack([a, a[::-1]])
        int_fluc = int_f0(xp, float(fs))[:len(a)]
        baseline = mean(a) * arange(len(a)) / float(fs)
        return int_fluc + baseline - int_fluc[0]
    
    # old version
    baseline = mean(a) * arange(len(a)) / float(fs)
    int_fluc = int_f0(a, float(fs))
    return int_fluc + baseline - int_fluc[0]

    # old code - remove eventually (comment on 02/2014)
    # periodify
    if False:
        baseline = linspace(a[0], a[-1], len(a))
        a0 = a - baseline
        m = a0[-1] - a0[-2]
        b2 = linspace(0, -.5 * m, len(a))
        baseline -= b2
        a0 += b2
        a2 = hstack([a0, -1. * a0[1:][::-1]]) # "smooth" periodic signal  

        dbase = baseline[1] - baseline[0]
        t_vec = arange(len(a)) / float(fs)
        baseint = baseline[0] * t_vec + .5 * dbase * t_vec ** 2
        
        # define frequencies
        T = len(a2) / float(fs)
        freqs = 1. / T * arange(len(a2))
        freqs[len(freqs) // 2 + 1 :] -= float(fs)

        spec = fft.fft(a2)
        spec_i = zeros_like(spec, dtype=complex)
        spec_i[1:] = spec[1:] / (2j * pi* freqs[1:])
        res_int = fft.ifft(spec_i).real[:len(a0)] + baseint
        return res_int - res_int[0]

def swap_subsample(I, k=1):
    for c, color in enumerate(colors):
        print "%s <-- %s" %(colors[c], colors[(c+k)%3])
    for i in range(3):
        J = zeros_like(I)
        for j in range(3):
            J[:, :, j] = I[:, :, (j+k)%3]
        J[:, :, i] = zoom(I[::4, ::4, (i+k)%3], 4)
        figure(i+10)
        title("%s channel subsampled" %colors[i])
        imshow(J)

def circle(center,radius,color):
	"""
		plot a circle with given center and radius

		Arguments
		----------
		center : matrix or ndarray
			it should be 2x1 ndarray or matrix
		radius:  float
			masses per mm
	"""
	u=pb.linspace(0,2*np.pi,200)
	x0=pb.zeros_like(u)
	y0=pb.zeros_like(u)
	center=pb.matrix(center)
	center.shape=2,1
	for i,j in enumerate(u):
		x0[i]=radius*pb.sin(j)+center[0,0]
		y0[i]=radius*pb.cos(j)+center[1,0]
	pb.plot(x0,y0,color)

def sub_mean(x, N):
    N = int(N)
    L = len(x)
    y = pl.zeros_like(x)
    ii = pl.arange(-N, N + 1)
    k = 1.0 / len(ii) # 1 / (2 * N + 1)
    for n in range(L):
        iii = pl.clip(ii + n, 0, L - 1)
        s = k * sum(x[iii])
        y[n] = x[n] - s
    print n, x[n], iii[0], iii[-1], s
    return y

def new_bad_model(F):
    """ Results in a matrix with shape matching X, but all rows sum to 1"""
    N, T, J = F.shape
    pi = pl.zeros_like(F)
    for t in range(T):
        u = F[:,t,:].var(axis=0)
        u /= pl.sqrt(pl.dot(u,u))
        F_t_par = pl.dot(pl.atleast_2d(pl.dot(F[:,t,:], u)).T, pl.atleast_2d(u))
        F_t_perp = F[:,t,:] - F_t_par
        for n in range(N):
            alpha = (1 - F_t_perp[n].sum()) / F_t_par[n].sum()
            pi[n,t,:] = F_t_perp[n,:] + alpha*F_t_par[n,:]
    return pi

def int_f0(x, fs=1.):
    """
    returns the 'basic' fourier integration of a signal
    
    """
    # define frequencies
    T = len(x) / float(fs)
    freqs = 1. / T * arange(len(x))
    freqs[len(freqs) // 2 + 1:] -= float(fs)

    spec = fft.fft(x)
    spec_i = zeros_like(spec, dtype=complex)
    # exclude frequency 0 - it cannot be integrated
    spec_i[1:] = spec[1:] / (2j * pi* freqs[1:])
    if mod(len(x), 2) == 0:
        spec_i[len(x) // 2] = 0.
    sig_d = fft.ifft(spec_i) 
    return sig_d.real

        def rates(S=data_sample,
                Xa=Xa, Xb=Xb,
                alpha=alpha, beta=beta, gamma=gamma,
                bounds_func=vars['bounds_func'],
                age_indices=ai,
                age_weights=aw):

            # calculate study-specific rate function
            shifts = pl.exp(pl.dot(Xa[S], alpha) + pl.dot(Xb[S], pl.atleast_1d(beta)))
            exp_gamma = pl.exp(gamma)
            mu = pl.zeros_like(shifts)
            for i,s in enumerate(S):
                mu[i] = pl.dot(age_weights[s], bounds_func(shifts[i] * exp_gamma[age_indices[s]], age_indices[s]))
                # TODO: evaluate speed increase and accuracy decrease of the following:
                #midpoint = age_indices[s][len(age_indices[s])/2]
                #mu[i] = bounds_func(shifts[i] * exp_gamma[midpoint], midpoint)
                # TODO: evaluate speed increase and accuracy decrease of the following: (to see speed increase, need to code this up using difference of running sums
                #mu[i] = pl.dot(pl.ones_like(age_weights[s]) / float(len(age_weights[s])),
                #               bounds_func(shifts[i] * exp_gamma[age_indices[s]], age_indices[s]))
            return mu

def connected_conponent(img):
    """
    input:  img cant be gray or color. 
            The shape of img is (height,width,..)
    output: list of components  and an image mask
    """
    img.resize((10,10))
    img = array(img)
    assert(img.ndim>=2)
    #mask = array([0]*img.size).reshape(img.shape)
    mask = zeros_like(img,dtype=uint64)  #[h,w]
    comp_list = []
    comp_index = 0   # the index of the first componet  is 1

    print mask.shape, mask.dtype
    print img.shape,img.dtype
    for y in range(img.shape[0]):
        for x in range(img.shape[1]):
            if mask[y,x]==0: # (y,x) will be a seed of component.
                comp_index+=1
                mask[y,x]=comp_index  
                #flood to ne
                flood(img,mask,(y,x))

    Itest = 5
    Vmax = 15

    #make grid
    plotMargin = 1.2
    sizex = Ca1.outer_radius*plotMargin
    sizez = Ca1.outer_thickness*plotMargin
    print sizex,sizez
    size = np.max([sizex,sizez])
    N=20.0    # number of grid points per axis
    i=j=pb.arange(N)
    n=float(N)
    x = ((i/(n-1))-.5)*size*2
    z = ((j/(n-1))-.5)*size*2
    X,Z = pb.meshgrid(x,z)
    Bx = pb.zeros_like(X)
    Bz = pb.zeros_like(X)
    Bnorm = pb.zeros_like(X)
    #get B field
    gaussPerTesla = 1e4    
    for ii in i:
        for jj in j:
            x = X[ii,jj]
            z = Z[ii,jj]
            if jj==0:
                print ii,jj,x,z
            B=(Ca1.Bvector([x,0,z])+Ca2.Bvector([x,0,z]))*gaussPerTesla
            Bnorm[ii,jj] = norm(B)
            Bx[ii,jj] = B[0]*np.log1p(np.fabs(1000*norm(B)))/norm(B)
            Bz[ii,jj] = B[2]*np.log1p(np.fabs(1000*norm(B)))/norm(B)
            

dm.description = 'With excess-mortality data removed'
dm.data = [d for d in dm.data if dm.relevant_to(d, 'all-cause_mortality', region, year, sex)] + \
    [d for d in dm.data if dm.relevant_to(d, 'prevalence_x_excess-mortality', region, year, sex)] + \
    [d for d in dm.data if d.get('country_iso3_code') == 'GBR']
dm.params['global_priors']['level_bounds']['excess_mortality'] = dict(lower=.2, upper=10.)
fit_posterior.fit_posterior(dm, region, sex, year, map_only=True, store_results=False)
models.append(dm)

dm = initialize_model()
dm.description = 'With excess-mortality data and priors removed'
dm.data = [d for d in dm.data if dm.relevant_to(d, 'all-cause_mortality', region, year, sex)] + \
    [d for d in dm.data if dm.relevant_to(d, 'prevalence_x_excess-mortality', region, year, sex)] + \
    [d for d in dm.data if d.get('country_iso3_code') == 'GBR']
dm.params['global_priors']['level_bounds']['excess_mortality'] = dict(lower=0., upper=10.)
dm.params['global_priors']['smoothness']['prevalence']['amount'] = 'No Prior'
dm.params['gamma_effect_excess-mortality'] = dict(mean=list(pl.zeros_like(dm.get_estimate_age_mesh())),
                                                  std=list(10.*pl.ones_like(dm.get_estimate_age_mesh())))
fit_posterior.fit_posterior(dm, region, sex, year, map_only=True, store_results=False)
models.append(dm)


dm = initialize_model()
dm.description = 'Without increasing prior on excess-mortality, but with "Very" smoothing'
dm.data = [d for d in dm.data if dm.relevant_to(d, 'all-cause_mortality', region, year, sex)] + \
    [d for d in dm.data if dm.relevant_to(d, 'prevalence_x_excess-mortality', region, year, sex)] + \
    [d for d in dm.data if d.get('country_iso3_code') == 'GBR'] + \
    [d for d in dm.data if dm.relevant_to(d, 'excess-mortality', 'all', 'all', 'all')]
dm.params['global_priors']['smoothness']['prevalence']['amount'] = 'Very'
fit_posterior.fit_posterior(dm, region, sex, year, map_only=True, store_results=False)
models.append(dm)