amplitude_model.py 3.87 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
import numpy as np

from ..domains import PowerSpace, UnstructuredDomain
from ..field import Field
from ..multi import MultiField
from ..sugar import makeOp, sqrt


def _ceps_kernel(dof_space, k, a, k0):
    return a**2/(1+(k/(k0*dof_space.bindistances[0]))**2)**2


def make_amplitude_model(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
                         keys=['tau', 'phi']):
    '''
    Method for construction of amplitude model
    Computes a smooth power spectrum.
    Output lives in PowerSpace.

    Parameters
    ----------

    Npixdof : #pix in dof_space

    ceps_a, ceps_k0 : Smoothnessparameters in ceps_kernel
                        eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2
                        a = ceps_a,  k0 = ceps_k0

    sm, sv : slope_mean = expected exponent of powerlaw (e.g. -4),
                slope_variance (default=1)

    im, iv : y-intercept_mean, y-intercept_variance  of power_slope
    '''
    from ..operators import (ExpTransform, QHTOperator, SlopeOperator,
                             SymmetrizingOperator)
36 37
    from ..models import Variable, Constant, PointwiseExponential

38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
    h_space = s_space.get_default_codomain()
    p_space = PowerSpace(h_space)
    exp_transform = ExpTransform(p_space, Npixdof)
    logk_space = exp_transform.domain[0]
    dof_space = logk_space.get_default_codomain()
    param_space = UnstructuredDomain(2)
    qht = QHTOperator(dof_space, logk_space)
    sym = SymmetrizingOperator(logk_space)

    phi_mean = np.array([sm, im])
    phi_sig = np.array([sv, iv])

    slope = SlopeOperator(param_space, logk_space, phi_sig)
    norm_phi_mean = Field(param_space, val=phi_mean/phi_sig)

    fields = {keys[0]: Field.from_random('normal', dof_space),
              keys[1]: Field.from_random('normal', param_space)}

    position = MultiField(fields)

    dof_space = position[keys[0]].domain[0]
    kern = lambda k: _ceps_kernel(dof_space, k, ceps_a, ceps_k)
    cepstrum = create_cepstrum_amplitude_field(dof_space, kern)

    ceps = makeOp(sqrt(cepstrum))
    smooth_op = sym * qht * ceps
    smooth_spec = smooth_op(Variable(position)[keys[0]])
65
    phi = Variable(position)[keys[1]] + Constant(position, norm_phi_mean)
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123
    linear_spec = slope(phi)
    loglog_spec = smooth_spec + linear_spec
    xlog_ampl = PointwiseExponential(0.5*loglog_spec)

    internals = {'loglog_spec': loglog_spec,
                 'qht': qht,
                 'ceps': ceps,
                 'norm_phi_mean': norm_phi_mean}
    return exp_transform(xlog_ampl), internals


def create_cepstrum_amplitude_field(domain, cepstrum):
    """Creates a ...
    Writes the sum of all modes into the zero-mode.

    Parameters
    ----------
    domain: ???
        ???
    cepstrum: Callable
        ???
    """

    dim = len(domain.shape)
    dist = domain.bindistances
    shape = domain.shape

    # Prepare q_array
    q_array = np.zeros((dim,) + shape)
    if dim == 1:
        ks = domain.get_k_length_array().val
        q_array = np.array([ks])
    else:
        for i in range(dim):
            ks = np.minimum(shape[i] - np.arange(shape[i]) +
                            1, np.arange(shape[i])) * dist[i]
            fst_dims = (1,) * i
            lst_dims = (1,) * (dim - i - 1)
            q_array[i] += ks.reshape(fst_dims + (shape[i],) + lst_dims)

    # Fill cepstrum field (all non-zero modes)
    no_zero_modes = (slice(1, None),) * dim
    ks = q_array[(slice(None),) + no_zero_modes]
    cepstrum_field = np.zeros(shape)
    cepstrum_field[no_zero_modes] = cepstrum(ks)

    # Fill cepstrum field (zero-mode subspaces)
    for i in range(dim):
        # Prepare indices
        fst_dims = (slice(None),) * i
        lst_dims = (slice(None),) * (dim - i - 1)
        sl = fst_dims + (slice(1, None),) + lst_dims
        sl2 = fst_dims + (0,) + lst_dims

        # Do summation
        cepstrum_field[sl2] = np.sum(cepstrum_field[sl], axis=i)

    return Field(domain, val=cepstrum_field)