Python Population.eval_state方法代码示例

# 需要导入模块: from population import Population [as 别名]
# 或者: from population.Population import eval_state [as 别名]
def postprocess(popn, x_inf, popn_true, x_true, options):
    """ Compute an ROC curve from a set of inferred samples and the true state
    """
    true_state = popn_true.eval_state(x_true)

    # Make sure we have a list of x's
    if not isinstance(x_inf, list):
        x_inf = [x_inf]

    inf_state = [popn.eval_state(x) for x in x_inf]


    # Check if the inference model is a standard GLM or a network GLM
    if
    # Now compute the true and false positive rates for MAP
    (map_tpr, map_fpr) = compute_roc_from_std_glm(true_state, map_state)
    map_tprs.append(map_tpr)
    map_fprs.append(map_fpr)

    print "Loading mcmc estimate"
    x_mcmc = None
    with open(os.path.join(data_dir, d, 'mcmc', 'results.pkl'), 'r') as f:
        x_mcmc = cPickle.load(f)


    model_mcmc = make_model('sparse_weighted_model', N=data['N'])
    popn_mcmc = Population(model_mcmc)
    popn_mcmc.set_data(data)

    # Evaluate the state
    mcmc_state = []
    for x in x_mcmc:
        mcmc_state.append(popn_mcmc.eval_state(x))


    # Now compute the true and false positive rates for MCMC
    # For MCMC results, only consider the tail of the samples
    N_samples = len(mcmc_state)
    sample_frac = 0.2
    start_smpl = int(np.floor(N_samples - sample_frac*N_samples))
    (mcmc_tpr, mcmc_fpr) = compute_roc_from_sparse_glm_smpls(true_state, mcmc_state[start_smpl:])
    mcmc_tprs.append(mcmc_tpr)
    mcmc_fprs.append(mcmc_fpr)


    # Pickle the roc results
    with open(PKL_FNAME, 'w') as f:
        # TODO Dump the MCMC results too
        cPickle.dump({'map_tprs' : map_tprs,
                      'map_fprs' : map_fprs},
                     f,
                     protocol=-1)

    # Plot the actual ROC curve
    # Subsample to get about 10 errorbars
    subsample = N*N//10
    f = plt.figure()
    ax = f.add_subplot(111)
    plot_roc_curve(map_tprs, map_fprs, ax=ax, color='b', subsample=subsample)
    # plot_roc_curve(mcmc_tprs, mcmc_fprs, ax=ax, color='r', subsample=subsample)
    fname = os.path.join(PLOTDIR,'roc_N%dT%d.pdf' % (N,T))
    print "Saving ROC to %s" % fname
    f.savefig(fname)
    plt.close(f)

# 需要导入模块: from population import Population [as 别名]
# 或者: from population.Population import eval_state [as 别名]
def load_set_of_results(N, T, graph_model='er', sample_frac=0.1):
    data_dir = os.path.join('/group', 'hips', 'scott', 'pyglm', 'data', 'synth', graph_model, 'N%dT%d' % (N, T))

    # Evaluate the state for each of the parameter settings
    s_infs_mcmc = []
    s_infs_map = []
    s_trues = []

    # Enumerate the subdirectories containing the data
    subdirs = os.listdir(data_dir)
    subdirs = reduce(lambda sd, d: sd + [d] \
                                   if os.path.isdir(os.path.join(data_dir, d)) \
                                   else sd,
                     subdirs, [])

    # For each data subdirectory, load the true data, the MAP estimate, and the MCMC results
    print "WARNING: Make sure we sample all subdirs"
    # import pdb; pdb.set_trace()
    for d in subdirs:
        print "Loading data and results from %s" % d
        print "Loading true data"
        with open(os.path.join(data_dir, d, 'data.pkl'), 'r') as f:
            data = cPickle.load(f)

        print "Loading model"
        with open(os.path.join(data_dir, d, 'model.pkl'), 'r') as f:
            model_data = cPickle.load(f)
            #HACK
            if 'N_dims' not in model_data['network']['graph']:
                model_data['network']['graph']['N_dims'] = 1
            if 'location_prior' not in model_data['network']['graph']:
                model_data['network']['graph']['location_prior'] = \
                         {
                             'type' : 'gaussian',
                             'mu' : 0.0,
                             'sigma' : 1.0
                         }
            if 'L' in data['vars']['net']['graph']:
                data['vars']['net']['graph']['L'] = data['vars']['net']['graph']['L'].ravel()
        popn_data = Population(model_data)
        popn_data.set_data(data)
        s_trues.append(popn_data.eval_state(data['vars']))

        try:
            print "Loading map estimate"
            with open(os.path.join(data_dir, d, 'map', 'results.pkl'), 'r') as f:
                x_map = cPickle.load(f)

            model_map = make_model('standard_glm', N=data['N'])
            popn_map = Population(model_map)
            popn_map.set_data(data)
            print "Evaluating MAP state"
            s_infs_map.append(popn_map.eval_state(x_map))

        except Exception as e:
            print "ERROR: Failed to load MAP estimate"

        try:
            print "Loading mcmc estimate"
            with open(os.path.join(data_dir, d, 'mcmc', 'results.pkl'), 'r') as f:
                x_mcmc = cPickle.load(f)

            model_mcmc = make_model('sparse_weighted_model', N=data['N'])
            popn_mcmc = Population(model_mcmc)
            popn_mcmc.set_data(data)

            # Now compute the true and false positive rates for MCMC
            # For MCMC results, only consider the tail of the samples
            print "Evaluating MCMC states"
            N_samples = len(x_mcmc)
            start_smpl = int(np.floor(N_samples - sample_frac*N_samples))

            # Evaluate the state
            this_s_mcmc = []
            for i in range(start_smpl, N_samples):
                this_s_mcmc.append(popn_mcmc.eval_state(x_mcmc[i]))
            s_infs_mcmc.append(this_s_mcmc)
        except Exception as e:
            print "ERROR: Failed to load MCMC estimate"

    return s_trues, s_infs_map, s_infs_mcmc

# 需要导入模块: from population import Population [as 别名]
# 或者: from population.Population import eval_state [as 别名]
def run_gen_synth_data():
    """ Run a test with synthetic data and MCMC inference
    """
    options, args = parse_cmd_line_args()
    
    # Create the model
    model = make_model(options.model, N=options.N)
    # Set the sparsity level to minimize the risk of unstable networks
    stabilize_sparsity(model)

    print "Creating master population object"
    popn = Population(model)

    # Sample random parameters from the model
    x_true = popn.sample()

    # Check stability of matrix
    assert check_stability(model, x_true, options.N), "ERROR: Sampled network is unstable!"


    # Save the model so it can be loaded alongside the data
    fname_model = os.path.join(options.resultsDir, 'model.pkl')
    print "Saving data to %s" % fname_model
    with open(fname_model,'w') as f:
        cPickle.dump(model, f, protocol=-1)

    print "Generating synthetic data with %d neurons and %.2f seconds." % \
          (options.N, options.T_stop)

    # Set simulation parametrs
    dt = 0.001
    dt_stim = 0.1
    D_stim = 1
    stim = np.random.randn(options.T_stop/dt_stim, D_stim)

    data = gen_synth_data(options.N, options.T_stop, popn, x_true, dt, dt_stim, D_stim, stim)
    # Set the data so that the population state can be evaluated
    popn.set_data(data)
    
    # DEBUG Evaluate the firing rate and the simulated firing rate
    state = popn.eval_state(x_true)
    for n in np.arange(options.N):
        lam_true = state['glms'][n]['lam']
        lam_sim =  popn.glm.nlin_model.f_nlin(data['X'][:,n])
        assert np.allclose(lam_true, lam_sim)

    # Save the data for reuse
    #fname_mat = os.path.join(options.resultsDir, 'data.mat')
    #print "Saving data to %s" % fname_mat
    #scipy.io.savemat(fname_mat, data, oned_as='row')
        
    # Pickle the data so we can open it more easily
    fname_pkl = os.path.join(options.resultsDir, 'data.pkl')
    print "Saving data to %s" % fname_pkl
    with open(fname_pkl,'w') as f:
        cPickle.dump(data, f, protocol=-1)

    # Plot firing rates, stimulus responses, etc
    plot_results(popn, data['vars'],
                 resdir=options.resultsDir,
                 do_plot_stim_resp=False,
                 do_plot_imp_responses=False)