How to use PyGLM - 10 common examples

def plot_connectivity(dataset, run, algs,):
    # Load the data and results
    train, test, true_model = load_data(dataset)
    res_dir = os.path.join("results", dataset, "run%03d" % run)
    results = load_results(dataset, run, algs)

    # samples = results["gibbs"][0]

    ###########################################################
    # Get the average connectivity
    ###########################################################
    if "bfgs" in algs:
        W_mean = results["bfgs"].W.sum(2)
    elif "gibbs" in algs:
        W_samples = [smpl.weight_model.W_effective
                         for smpl in results["gibbs"][0]]

        offset = len(W_samples) // 2
        W_samples = np.array(W_samples[offset:])

def plot_impulse_responses(dataset, run):
        # Load the data and results
    from pyglm.models import BernoulliEigenmodelPopulation
    model = BernoulliEigenmodelPopulation(
        N=27, dt=1.0, dt_max=10.0, B=5)
    res_dir = os.path.join("results", dataset, "run%03d" % run)
    results = load_results(dataset, run, ["gibbs"])
    samples = results["gibbs"][0]

    ###########################################################
    # Get the average connectivity
    ###########################################################
    W_samples = [smpl.weight_model.W_effective for smpl in samples]
    offset = len(W_samples) // 2
    W_samples = np.array(W_samples[offset:])
    W_lim = np.amax(abs(W_samples))
    imp_lim = W_lim * np.amax(abs(model.basis.basis))

    syn_sort = np.argsort(W_samples.mean(0).sum(2).ravel())
    i_synapses = zip(*np.unravel_index(syn_sort[:3], (27,27)))
    e_synapses = zip(*np.unravel_index(syn_sort[-3:], (27,27)))
    synapses = i_synapses + e_synapses

###########################################################
    #   Network hyperparameters
    ###########################################################
    # c = np.arange(C).repeat((N // C))                       # Neuron to cluster assignments
    # p = 0.5 * np.ones((C,C))                                      # Probability of connection for each pair of clusters
    # # p = 0.9 * np.eye(C) + 0.05 * (1-np.eye(C))              # Probability of connection for each pair of clusters
    # mu = np.zeros((C,C,B))                                  # Mean weight for each pair of clusters
    # Sigma = np.tile( 3**2 * np.eye(B)[None,None,:,:], (C,C,1,1))    # Covariance of weight for each pair of clusters
    # network_hypers = {'C': C, 'c': c, 'p': p, 'mu': mu, 'Sigma': Sigma}

    ###########################################################
    # Create the model with these parameters
    ###########################################################
    network_hypers = {"Sigma_0": 10*np.eye(B)}
    true_model = Population(N=N, dt=dt, dt_max=dt_max, B=B,
                            bias_hypers=bias_hypers,
                            network_hypers=network_hypers)

    ###########################################################
    # Sample from the true model
    ###########################################################
    S, Psi_hat = true_model.generate(T=T, keep=True, return_Psi=True)
    print "Number of generated spikes:\t", S.sum(0)

    ###########################################################
    #  Plot the network, the spike train and mean rate
    ###########################################################
    plt.figure()
    plt.imshow(true_model.weight_model.W_effective.sum(2), vmin=-1.0, vmax=1.0, interpolation="none", cmap="RdGy")

    plt.figure()

def run_synth_test():
    """ Run a test with synthetic data and MAP inference with cross validation
    """
    options, popn, data, popn_true, x_true = initialize_test_harness()
    
    # Get the list of models for cross validation
    base_model = make_model(options.model, N=data['N'], dt=0.001)
    models = get_xv_models(base_model)

    # TODO Segment data into training and cross validation sets
    train_frac = 0.75
    T_split = data['T'] * train_frac
    train_data = segment_data(data, (0,T_split))
    xv_data = segment_data(data, (T_split,data['T']))

    # Preprocess the data sequences
    train_data = popn.preprocess_data(train_data)
    xv_data = popn.preprocess_data(xv_data)

    # Sample random initial state
    x0 = popn.sample()

    # Track the best model and parameters

def run_synth_test():
    """ Run a test with synthetic data and MAP inference with cross validation
    """
    options, popn, data, popn_true, x_true = initialize_test_harness()
    
    # Get the list of models for cross validation
    base_model = make_model(options.model, N=data['N'], dt=0.001)
    models = get_xv_models(base_model)

    # TODO Segment data into training and cross validation sets
    train_frac = 0.75
    T_split = data['T'] * train_frac
    train_data = segment_data(data, (0,T_split))
    xv_data = segment_data(data, (T_split,data['T']))

    # Preprocess the data sequences
    train_data = popn.preprocess_data(train_data)
    xv_data = popn.preprocess_data(xv_data)

    # Sample random initial state
    x0 = popn.sample()

    # Track the best model and parameters
    best_ind = -1
    best_xv_ll = -np.Inf
    best_x = x0
    best_model = None

    # Fit each model using the optimum of the previous models
    train_lls = np.zeros(len(models))

def run_synth_test():
    """ Run a test with synthetic data and MAP inference with cross validation
    """
    options, popn, data, popn_true, x_true = initialize_test_harness()
    
    # Get the list of models for cross validation
    base_model = make_model(options.model, N=data['N'], dt=0.001)
    models = get_xv_models(base_model)

    # TODO Segment data into training and cross validation sets
    train_frac = 0.75
    T_split = data['T'] * train_frac
    train_data = segment_data(data, (0,T_split))
    xv_data = segment_data(data, (T_split,data['T']))

    # Preprocess the data sequences
    train_data = popn.preprocess_data(train_data)
    xv_data = popn.preprocess_data(xv_data)

    # Sample random initial state
    x0 = popn.sample()

    # Track the best model and parameters
    best_ind = -1
    best_xv_ll = -np.Inf
    best_x = x0
    best_model = None

    # Fit each model using the optimum of the previous models

def run_synth_test():
    """ Run a test with synthetic data and MCMC inference
    """
    options, popn, data, popn_true, x_true = initialize_test_harness()

    # Sample random initial state
    x0 = popn.sample()
    ll0 = popn.compute_log_p(x0)
    print "LL0: %f" % ll0

    # Perform inference
    x_inf = coord_descent(popn, x0=x0, maxiter=1)
    ll_inf = popn.compute_log_p(x_inf)
    print "LL_inf: %f" % ll_inf

    # Save results
    results_file = os.path.join(options.resultsDir, 'results.pkl')
    print "Saving results to %s" % results_file
    with open(results_file, 'w') as f:
        cPickle.dump(x_inf, f, protocol=-1)

    # Plot results
    plot_results(popn, x_inf, popn_true, x_true, resdir=options.resultsDir)

def resample_x(z):
    # Resample x from its (scaled) Bernoulli lkhd
    p = logistic(z / T)
    x_smpl = np.random.rand() < p
    return x_smpl

def update(model):
    model.resample_model()
    z_inf = model.states_list[0].stateseq
    C_inf = model.C
    psi_inf = z_inf.dot(C_inf.T)
    p_inf = logistic(psi_inf)
    return model.log_likelihood(), p_inf

B      = true_model.B
    dt     = true_model.dt
    dt_max = true_model.dt_max

    ###########################################################
    # Create and fit a standard model for initialization
    ###########################################################
    with gzip.open(init_path, 'r') as f:
        init_model = cPickle.load(f)

    ###########################################################
    # Create a test spike-and-slab model
    ###########################################################

    # Copy the network hypers.
    test_model = NegativeBinomialEigenmodelPopulation(N=N, dt=dt, dt_max=dt_max, B=B,
                            basis_hypers=true_model.basis_hypers,
                            observation_hypers=true_model.observation_hypers,
                            activation_hypers=true_model.activation_hypers,
                            weight_hypers=true_model.weight_hypers,
                            bias_hypers=true_model.bias_hypers,
                            network_hypers=true_model.network_hypers)
    test_model.add_data(S)

    # Initialize with the standard model
    test_model.initialize_with_standard_model(init_model)

    # Convolve the test data for fast heldout likelihood calculations
    F_test = test_model.basis.convolve_with_basis(S_test)

    ###########################################################
    # Fit the test model with Gibbs sampling