def run_comparison(data_path, test_path, output_path, T_train=None, seed=None):
    if seed is None:
        seed = np.random.randint(2**32)

    print "Setting seed to ", seed
    np.random.seed(seed)

    assert os.path.exists(os.path.dirname(output_path)), "Output directory does not exist!"

    if data_path.endswith(".gz"):
        with gzip.open(data_path, 'r') as f:
            S, true_model = cPickle.load(f)
    else:
        with open(data_path, 'r') as f:
            S, true_model = cPickle.load(f)

    # If T_train is given, only use a fraction of the dataset
    if T_train is not None:
        S = S[:T_train,:]

    if test_path.endswith(".gz"):
        with gzip.open(test_path, 'r') as f:
            S_test, test_model = cPickle.load(f)
    else:
        with open(test_path, 'r') as f:
            S_test, test_model = cPickle.load(f)

    K      = true_model.K
    C      = true_model.C
    B      = true_model.B
    dt     = true_model.dt
    dt_max = true_model.dt_max

    # Fix the sparsity for the spike and slab model
    p      = 0.4 * np.eye(C) + 0.01 * (1-np.eye(C))

    use_parse_results = False
    if use_parse_results and  os.path.exists(output_path + ".parsed_results.pkl"):
        with open(output_path + ".parsed_results.pkl") as f:
            auc_rocs, auc_prcs, plls, timestamps = cPickle.load(f)
            timestamps['svi'] = np.array(timestamps['svi'])
        import pdb; pdb.set_trace()

    else:
        # Compute the cross correlation to estimate the connectivity
        W_xcorr = infer_net_from_xcorr(S, dtmax=true_model.dt_max // true_model.dt)

        # Fit a standard Hawkes model on subset of data with BFGS
        bfgs_model, bfgs_time = fit_standard_hawkes_model_bfgs(S, K, B, dt, dt_max,
                                                               output_path=output_path)

        # Fit a standard Hawkes model with SGD
        # standard_models, timestamps = fit_standard_hawkes_model_sgd(S, K, B, dt, dt_max,
        #                                                         init_model=init_model)
        #
        # # Save the models
        # with open(output_path + ".sgd.pkl", 'w') as f:
        #     print "Saving SGD results to ", (output_path + ".sgd.pkl")
        #     cPickle.dump((standard_models, timestamps), f, protocol=-1)

        # Fit a network Hawkes model with Gibbs
        gibbs_samples = gibbs_timestamps = None
        # gibbs_samples, gibbs_timestamps = fit_network_hawkes_gibbs(S, K, C, B, dt, dt_max,
        #                                                      output_path=output_path,
        #                                                      standard_model=bfgs_model)

        # Fit a spike and slab network Hawkes model with Gibbs
        gibbs_ss_samples, gibbs_ss_timestamps = fit_network_hawkes_gibbs_ss(S, K, C, B, dt, dt_max,
                                                                 output_path=output_path,
                                                                 p=p,
                                                                 standard_model=bfgs_model)

        # Fit a network Hawkes model with Batch VB
        vb_models, vb_timestamps = fit_network_hawkes_vb(S, K, C, B, dt, dt_max,
                                                      output_path=output_path,
                                                      standard_model=bfgs_model)

        # Fit a network Hawkes model with SVI
        # svi_models = svi_timestamps = None
        svi_models, svi_timestamps = fit_network_hawkes_svi(S, K, C, B, dt, dt_max,
                                                        output_path,
                                                        standard_model=bfgs_model)

        # Combine timestamps into a dict
        timestamps = {}
        timestamps['bfgs'] = bfgs_time
        timestamps['gibbs'] = gibbs_timestamps
        timestamps['gibbs_ss'] = gibbs_ss_timestamps
        timestamps['svi'] = svi_timestamps
        timestamps['vb'] = vb_timestamps

        auc_rocs = compute_auc(true_model,
                           W_xcorr=W_xcorr,
                           bfgs_model=bfgs_model,
                           gibbs_samples=gibbs_samples,
                           gibbs_ss_samples=gibbs_ss_samples,
                           svi_models=svi_models,
                           vb_models=vb_models)
        print "AUC-ROC"
        pprint.pprint(auc_rocs)

        # Compute area under precisino recall curve of inferred network
        auc_prcs = compute_auc_prc(true_model,
                                   W_xcorr=W_xcorr,
                                   bfgs_model=bfgs_model,
                                   gibbs_samples=gibbs_samples,
                                   gibbs_ss_samples=gibbs_ss_samples,
                                   svi_models=svi_models,
                                   vb_models=vb_models)
        print "AUC-PRC"
        pprint.pprint(auc_prcs)


        plls = compute_predictive_ll(S_test, S,
                                     true_model=true_model,
                                     bfgs_model=bfgs_model,
                                     gibbs_samples=gibbs_samples,
                                     gibbs_ss_samples=gibbs_ss_samples,
                                     svi_models=svi_models,
                                     vb_models=vb_models)

        with open(output_path + ".parsed_results.pkl", 'w') as f:
            print "Saving parsed results to ", output_path + ".parsed_results.pkl"
            cPickle.dump((auc_rocs, auc_prcs, plls, timestamps), f, protocol=-1)

    plot_pred_ll_vs_time(plls, timestamps, Z=float(S.size), T_train=T_train)
Run the comparison on the given data file
dataset/ETHPy150Open slinderman/pyhawkes/experiments/synthetic_comparison.py/run_comparison