Python autograd.numpy.cumsum() Examples

def simple_admixture_3pop(x=None):
    if x is None:
        x = np.random.normal(size=7)
    t = np.cumsum(np.exp(x[:5]))
    p = 1.0 / (1.0 + np.exp(x[5:]))

    model = momi.DemographicModel(1., .25)
    model.add_leaf("b")
    model.add_leaf("a")
    model.add_leaf("c")
    model.move_lineages("a", "c", t[1], p=1.-p[1])
    model.move_lineages("a", "d", t[0], p=1.-p[0])
    model.move_lineages("c", "d", t[2])
    model.move_lineages("d", "b", t[3])
    model.move_lineages("a", "b", t[4])
    return model 

def simple_admixture_demo(x=np.random.normal(size=7)):
    t = np.cumsum(np.exp(x[:5]))
    p = 1.0 / (1.0 + np.exp(x[5:]))

    ret = momi.DemographicModel(1., .25)
    ret.add_leaf("b")
    ret.add_leaf("a")
    ret.move_lineages("a", 2, t[1], p=1.-p[1])
    ret.move_lineages("a", 3, t[0], p=1.-p[0])
    ret.move_lineages(2, 3, t[2])
    ret.move_lineages(3, "b", t[3])
    ret.move_lineages("a", "b", t[4])
    return ret 

def resampling(w, rs):
    """
    Stratified resampling with "nograd_primitive" to ensure autograd 
    takes no derivatives through it.
    """
    N = w.shape[0]
    bins = np.cumsum(w)
    ind = np.arange(N)
    u = (ind  + rs.rand(N))/N
    
    return np.digitize(u, bins) 

def get_downsampled_mcl(self, mcl_fractions):
        # Returns the mean camber line in downsampled form

        mcl = self.mcl_coordinates
        # Find distances along mcl, assuming linear interpolation
        mcl_distances_between_points = np.sqrt(
            np.power(mcl[:-1, 0] - mcl[1:, 0], 2) +
            np.power(mcl[:-1, 1] - mcl[1:, 1], 2)
        )
        mcl_distances_cumulative = np.hstack((0, np.cumsum(mcl_distances_between_points)))
        mcl_distances_cumulative_normalized = mcl_distances_cumulative / mcl_distances_cumulative[-1]

        mcl_downsampled_x = np.interp(
            x=mcl_fractions,
            xp=mcl_distances_cumulative_normalized,
            fp=mcl[:, 0]
        )
        mcl_downsampled_y = np.interp(
            x=mcl_fractions,
            xp=mcl_distances_cumulative_normalized,
            fp=mcl[:, 1]
        )

        mcl_downsampled = np.column_stack((mcl_downsampled_x, mcl_downsampled_y))

        return mcl_downsampled 

def repanel_current_airfoil(self, n_points_per_side=100):
        # Returns a repaneled version of the airfoil with cosine-spaced coordinates on the upper and lower surfaces.
        # Inputs:
        #   # n_points_per_side is the number of points PER SIDE (upper and lower) of the airfoil. 100 is a good number.
        # Notes: The number of points defining the final airfoil will be n_points_per_side*2-1,
        # since one point (the leading edge point) is shared by both the upper and lower surfaces.

        upper_original_coors = self.upper_coordinates()  # Note: includes leading edge point, be careful about duplicates
        lower_original_coors = self.lower_coordinates()  # Note: includes leading edge point, be careful about duplicates

        # Find distances between coordinates, assuming linear interpolation
        upper_distances_between_points = np.sqrt(
            np.power(upper_original_coors[:-1, 0] - upper_original_coors[1:, 0], 2) +
            np.power(upper_original_coors[:-1, 1] - upper_original_coors[1:, 1], 2)
        )
        lower_distances_between_points = np.sqrt(
            np.power(lower_original_coors[:-1, 0] - lower_original_coors[1:, 0], 2) +
            np.power(lower_original_coors[:-1, 1] - lower_original_coors[1:, 1], 2)
        )
        upper_distances_from_TE = np.hstack((0, np.cumsum(upper_distances_between_points)))
        lower_distances_from_LE = np.hstack((0, np.cumsum(lower_distances_between_points)))
        upper_distances_from_TE_normalized = upper_distances_from_TE / upper_distances_from_TE[-1]
        lower_distances_from_LE_normalized = lower_distances_from_LE / lower_distances_from_LE[-1]

        # Generate a cosine-spaced list of points from 0 to 1
        s = cosspace(n_points=n_points_per_side)

        x_upper_func = sp_interp.PchipInterpolator(x=upper_distances_from_TE_normalized, y=upper_original_coors[:, 0])
        y_upper_func = sp_interp.PchipInterpolator(x=upper_distances_from_TE_normalized, y=upper_original_coors[:, 1])
        x_lower_func = sp_interp.PchipInterpolator(x=lower_distances_from_LE_normalized, y=lower_original_coors[:, 0])
        y_lower_func = sp_interp.PchipInterpolator(x=lower_distances_from_LE_normalized, y=lower_original_coors[:, 1])

        x_coors = np.hstack((x_upper_func(s), x_lower_func(s)[1:]))
        y_coors = np.hstack((y_upper_func(s), y_lower_func(s)[1:]))

        coordinates = np.column_stack((x_coors, y_coors))
        self.coordinates = coordinates 

def test_cumsum_axis0():
    def fun(x): return np.cumsum(x, axis=0)
    mat = npr.randn(10, 11)
    check_grads(fun)(mat) 

def test_cumsum_axis1():
    def fun(x): return np.cumsum(x, axis=1)
    mat = npr.randn(10, 11)
    check_grads(fun)(mat) 

def test_cumsum_1d():
    def fun(x): return np.cumsum(x)
    mat = npr.randn(10)
    check_grads(fun)(mat) 

def test_cumsum_no_axis():
    def fun(x): return np.cumsum(x)
    mat = npr.randn(10, 11)
    check_grads(fun)(mat) 

def _split(x, shapes):
    x_split = np.split(x, np.cumsum([np.prod(shape) for shape in shapes[:-1]]))
    optim_vars = [var.reshape(*shape) for (var, shape) in zip(x_split, shapes)]
    return optim_vars 

def simple_five_pop_demo(x=np.random.normal(size=30)):
    assert len(x) == 30
    # make all params positive
    x = np.exp(x)

    # # allow negative growth rates
    # for i in range(15,20):
    #     x[i] = np.log(x[i])
    # # make times increasing
    # for i in range(1,15):
    #     x[i] = x[i] + x[i-1]

    t = np.cumsum(x[:15])
    # allow negative growth rates
    g = np.log(x[15:20])

    model = momi.DemographicModel(1.0, .25)
    for pop in range(1, 6):
        model.add_leaf(pop)
    model.set_size(5, t[0], g=g[0])
    model.set_size(4, t[1], g=g[1])
    model.set_size(3, t[2], g=g[2])
    model.set_size(2, t[3], g=g[3])
    model.set_size(1, t[4], g=g[4])
    model.move_lineages(5, 4, t=t[5], N=x[20])
    model.set_size(3, t=t[6], N=x[21])
    model.set_size(2, t=t[7], N=x[22])
    model.set_size(1, t[8], N=x[23])
    model.move_lineages(4, 3, t[9], N=x[24])
    model.set_size(2, t[10], N=x[25])
    model.set_size(1, t[11], N=x[26])
    model.move_lineages(3, 2, t[12], N=x[27])
    model.set_size(1, t[13], N=x[28])
    model.move_lineages(2, 1, t[14], N=x[29])
    return model
    ## number of edges is 2n-1
    #events_list = [('-eg', t[0], 5, g[0]),
    #               ('-eg', t[1], 4, g[1]),
    #               ('-eg', t[2], 3, g[2]),
    #               ('-eg', t[3], 2, g[3]),
    #               ('-eg', t[4], 1, g[4]),
    #               ('-ej', t[5], 5, 4), ('-en', t[5], 4, x[20]),
    #               ('-en', t[6], 3, x[21]),
    #               ('-en', t[7], 2, x[22]),
    #               ('-en', t[8], 1, x[23]),
    #               ('-ej', t[9], 4, 3), ('-en', t[9], 3, x[24]),
    #               ('-en', t[10], 2, x[25]),
    #               ('-en', t[11], 1, x[26]),
    #               ('-ej', t[12], 3, 2), ('-en', t[12], 2, x[27]),
    #               ('-en', t[13], 1, x[28]),
    #               ('-ej', t[14], 2, 1), ('-en', t[14], 1, x[29])]
    #demo = make_demo_hist(events_list, sampled_pops=list(
    #    range(1, len(n_lins) + 1)), sampled_n=n_lins)
    ##demo = make_demography(events_list, sampled_pops = list(range(1,len(n_lins)+1)), sampled_n = n_lins)
    ##demo = demo.rescaled(1e4)
    #return demo 

def sim_q(prop_params, model_params, y, smc_obj, rs, verbose=False):
    """
    Simulates a single sample from the VSMC approximation.

    Requires an SMC object with 2 member functions:
    -- sim_prop(t, x_{t-1}, y, prop_params, model_params, rs)
    -- log_weights(t, x_t, x_{t-1}, y, prop_params, model_params)
    """
    # Extract constants
    T = y.shape[0]
    Dx = smc_obj.Dx
    N = smc_obj.N

    # Initialize SMC
    X = np.zeros((N,T,Dx))
    logW = np.zeros(N)
    W = np.zeros((N,T))
    ESS = np.zeros(T)

    for t in range(T):
        # Resampling
        if t > 0:
            ancestors = resampling(W[:,t-1], rs)
            X[:,:t,:] = X[ancestors,:t,:]

        # Propagation
        X[:,t,:] = smc_obj.sim_prop(t, X[:,t-1,:], y, prop_params, model_params, rs)

        # Weighting
        logW = smc_obj.log_weights(t, X[:,t,:], X[:,t-1,:], y, prop_params, model_params)
        max_logW = np.max(logW)
        W[:,t] = np.exp(logW-max_logW)
        W[:,t] /= np.sum(W[:,t])
        ESS[t] = 1./np.sum(W[:,t]**2)

    # Sample from the empirical approximation
    bins = np.cumsum(W[:,-1])
    u = rs.rand()
    B = np.digitize(u,bins)

    if verbose:
        print('Mean ESS', np.mean(ESS)/N)
        print('Min ESS', np.min(ESS))
        
    return X[B,:,:] 

def get_repaneled_airfoil(self, n_points_per_side=100):
        # Returns a repaneled version of the airfoil with cosine-spaced coordinates on the upper and lower surfaces.
        # Inputs:
        #   # n_points_per_side is the number of points PER SIDE (upper and lower) of the airfoil. 100 is a good number.
        # Notes: The number of points defining the final airfoil will be n_points_per_side*2-1,
        # since one point (the leading edge point) is shared by both the upper and lower surfaces.

        upper_original_coors = self.upper_coordinates()  # Note: includes leading edge point, be careful about duplicates
        lower_original_coors = self.lower_coordinates()  # Note: includes leading edge point, be careful about duplicates

        # Find distances between coordinates, assuming linear interpolation
        upper_distances_between_points = np.sqrt(
            np.power(upper_original_coors[:-1, 0] - upper_original_coors[1:, 0], 2) +
            np.power(upper_original_coors[:-1, 1] - upper_original_coors[1:, 1], 2)
        )
        lower_distances_between_points = np.sqrt(
            np.power(lower_original_coors[:-1, 0] - lower_original_coors[1:, 0], 2) +
            np.power(lower_original_coors[:-1, 1] - lower_original_coors[1:, 1], 2)
        )
        upper_distances_from_TE = np.hstack((0, np.cumsum(upper_distances_between_points)))
        lower_distances_from_LE = np.hstack((0, np.cumsum(lower_distances_between_points)))
        upper_distances_from_TE_normalized = upper_distances_from_TE / upper_distances_from_TE[-1]
        lower_distances_from_LE_normalized = lower_distances_from_LE / lower_distances_from_LE[-1]

        # Generate a cosine-spaced list of points from 0 to 1
        s = cosspace(n_points=n_points_per_side)

        x_upper_func = sp_interp.PchipInterpolator(x=upper_distances_from_TE_normalized, y=upper_original_coors[:, 0])
        y_upper_func = sp_interp.PchipInterpolator(x=upper_distances_from_TE_normalized, y=upper_original_coors[:, 1])
        x_lower_func = sp_interp.PchipInterpolator(x=lower_distances_from_LE_normalized, y=lower_original_coors[:, 0])
        y_lower_func = sp_interp.PchipInterpolator(x=lower_distances_from_LE_normalized, y=lower_original_coors[:, 1])

        x_coors = np.hstack((x_upper_func(s), x_lower_func(s)[1:]))
        y_coors = np.hstack((y_upper_func(s), y_lower_func(s)[1:]))

        coordinates = np.column_stack((x_coors, y_coors))

        # Make a new airfoil with the coordinates
        name = self.name + ", repaneled to " + str(n_points_per_side) + " pts"
        new_airfoil = Airfoil(name=name, coordinates=coordinates, repanel=False)

        return new_airfoil 

def sample(self, n_samples=2000, observed_states=None, random_state=None):
        """Generate random samples from the self.

        Parameters
        ----------
        n : int
            Number of samples to generate.

        observed_states : array
            If provided, states are not sampled.

        random_state: RandomState or an int seed
            A random number generator instance. If None is given, the
            object's random_state is used

        Returns
        -------
        samples : array_like, length (``n_samples``, ``n_features``)
                  List of samples

        states : array_like, shape (``n_samples``)
                 List of hidden states (accounting for tied states by giving
                 them the same index)
        """

        if random_state is None:
            random_state = self.random_state
        random_state = check_random_state(random_state)

        samples = np.zeros((n_samples, self.n_features))
        states = np.zeros(n_samples)

        if observed_states is None:
            startprob_pdf = np.exp(np.copy(self._log_startprob))
            startdist = stats.rv_discrete(name='custm',
                                      values=(np.arange(startprob_pdf.shape[0]),
                                                        startprob_pdf),
                                      seed=random_state)
            states[0] = startdist.rvs(size=1)[0]

            transmat_pdf = np.exp(np.copy(self._log_transmat))
            transmat_cdf = np.cumsum(transmat_pdf, 1)

            nrand = random_state.rand(n_samples)
            for idx in range(1,n_samples):
                newstate = (transmat_cdf[int(states[idx-1])] > nrand[idx-1]).argmax()
                states[idx] = newstate

        else:
            states = observed_states

        mu = np.copy(self._mu_)
        precision = np.copy(self._precision_)
        for idx in range(n_samples):
            mean_ = mu[states[idx]]

            covar_ = np.linalg.inv(precision[states[idx]])
            samples[idx] = multivariate_normal.rvs(mean=mean_, cov=covar_,
                                                   random_state=random_state)
        states = self._process_sequence(states)
        return samples, states