Python RUtil类代码示例

def get_r_tikz_corr_plot(nsels, time_stats):
    """
    @param time_stats: a list of stats for each time point
    @return: tikz code corresponding to an R plot
    """
    out = StringIO()
    time_stats_trans = zip(*time_stats)
    y_low = -1
    y_high = 1
    ylim = RUtil.mk_call_str("c", y_low, y_high)
    print >> out, RUtil.mk_call_str(
        "plot",
        "my.table$t",
        "my.table$corr.mi.diag.approx",
        type='"n"',
        ylim=ylim,
        xlab='"time"',
        ylab='"correlation"',
        main='"correlation with mutual information"',
    )
    colors = ("red", "orange", "green", "blue", "black")
    plot_indices = (7, 8, 9, 10, 11)
    for c, plot_index in zip(colors, plot_indices):
        header = g_time_stats_headers[plot_index]
        print >> out, RUtil.mk_call_str("lines", "my.table$t", "my.table$%s" % header, col='"%s"' % c)
    return out.getvalue()

def get_r_tikz_mi_plot_script(nsels, time_stats):
    """
    At each time point plot mutual information for all matrices.
    @param time_stats: a list of stats for each time point
    @return: tikz code corresponding to an R plot
    """
    out = StringIO()
    time_stats_trans = zip(*time_stats)
    mi_mut = time_stats_trans[1]
    mi_min_sels = time_stats_trans[6]
    mi_max_sels = time_stats_trans[2]
    y_low = min(mi_min_sels + mi_mut)
    y_high = max(mi_max_sels + mi_mut)
    ylim = RUtil.mk_call_str("c", y_low, y_high)
    print >> out, RUtil.mk_call_str(
        "plot",
        "my.table$t",
        "my.table$mut",
        type='"n"',
        ylim=ylim,
        xlab='"time"',
        ylab='"MI"',
        main='"MI for mut process and %d mut.sel processes"' % nsels,
    )
    colors = ("red", "blue", "green", "black", "green", "blue")
    plot_indices = (1, 2, 3, 4, 5, 6)
    for c, plot_index in zip(colors, plot_indices):
        header = g_time_stats_headers[plot_index]
        print >> out, RUtil.mk_call_str("lines", "my.table$t", "my.table$%s" % header, col='"%s"' % c)
    return out.getvalue()

def get_r_tikz_info_plot(nsels, time_stats):
    """
    @param time_stats: a list of stats for each time point
    @return: tikz code corresponding to an R plot
    """
    out = StringIO()
    time_stats_trans = zip(*time_stats)
    y_low = 0
    y_high = math.log(2)
    ylim = RUtil.mk_call_str("c", y_low, y_high)
    print >> out, RUtil.mk_call_str(
        "plot",
        "my.table$t",
        "my.table$info.mi.diag.approx",
        type='"n"',
        ylim=ylim,
        xlab='"time"',
        ylab='"info"',
        main='"informativeness with respect to MI"',
    )
    colors = ("red", "orange", "green", "blue", "black")
    plot_indices = (17, 18, 19, 20, 21)
    for c, plot_index in zip(colors, plot_indices):
        header = g_time_stats_headers[plot_index]
        print >> out, RUtil.mk_call_str("lines", "my.table$t", "my.table$%s" % header, col='"%s"' % c)
    return out.getvalue()

def get_latex_documentbody(fs):
    """
    This is obsolete because I am now using pure R output.
    The latex documentbody should have a bunch of tikz pieces in it.
    Each tikz piece should have been generated from R.
    """
    Q_mut, Q_sels = get_qmut_qsels(fs)
    # compute the statistics
    ER_ratios, NSR_ratios, ER_NSR_ratios  = get_statistic_ratios(Q_mut, Q_sels)
    M = zip(*(ER_ratios, NSR_ratios, ER_NSR_ratios))
    column_headers = ('ER.ratio', 'NSR.ratio', 'ER.times.NSR.ratio')
    table_string = RUtil.get_table_string(M, column_headers)
    nsels = len(Q_sels)
    # define the R scripts
    scripts = []
    for name in column_headers:
        scripts.append(get_r_tikz_script(nsels, name))
    # get the tikz codes from R, for each histogram
    retcode, r_out, r_err, tikz_code_list = RUtil.run_plotter_multiple_scripts(
            table_string, scripts, 'tikz',
            width=3, height=2)
    if retcode:
        raise RUtil.RError(r_err)
    #
    # show some timings
    print 'R did not fail, but here is its stderr:'
    print r_err
    #
    # write the latex code
    out = StringIO()
    #print >> out, '\\pagestyle{empty}'
    for tikz_code in tikz_code_list:
        print >> out, tikz_code
    # return the latex code, consisting mainly of a bunch of tikz plots
    return out.getvalue()

def get_response_content(fs):
    M, R = get_input_matrices(fs)
    # create the R table string and scripts
    headers = [
            't',
            'mi.true.mut',
            'mi.true.mutsel',
            'mi.analog.mut',
            'mi.analog.mutsel']
    npoints = 100
    t_low = 0.0
    t_high = 5.0
    t_incr = (t_high - t_low) / (npoints - 1)
    t_values = [t_low + t_incr*i for i in range(npoints)]
    # get the data for the R table
    arr = []
    for t in t_values:
        mi_mut = ctmcmi.get_mutual_information(M, t)
        mi_mutsel = ctmcmi.get_mutual_information(R, t)
        mi_analog_mut = ctmcmi.get_ll_ratio_wrong(M, t)
        mi_analog_mutsel = ctmcmi.get_ll_ratio_wrong(R, t)
        row = [t, mi_mut, mi_mutsel, mi_analog_mut, mi_analog_mutsel]
        arr.append(row)
    # get the R table
    table_string = RUtil.get_table_string(arr, headers)
    # get the R script
    script = get_ggplot()
    # create the R plot image
    device_name = Form.g_imageformat_to_r_function[fs.imageformat]
    retcode, r_out, r_err, image_data = RUtil.run_plotter(
            table_string, script, device_name)
    if retcode:
        raise RUtil.RError(r_err)
    return image_data

def get_response_content(fs):
    # precompute some transition matrices
    P_drift_selection = pgmsinglesite.create_drift_selection_transition_matrix(
            fs.npop, fs.selection_ratio)
    MatrixUtil.assert_transition_matrix(P_drift_selection)
    P_mutation = pgmsinglesite.create_mutation_transition_matrix(
            fs.npop, fs.mutation_ab, fs.mutation_ba)
    MatrixUtil.assert_transition_matrix(P_mutation)
    # define the R table headers
    headers = ['generation', 'number.of.mutants']
    # compute the path samples
    P = np.dot(P_drift_selection, P_mutation)
    mypath = PathSampler.sample_endpoint_conditioned_path(
            fs.nmutants_initial, fs.nmutants_final, fs.ngenerations, P)
    arr = [[i, nmutants] for i, nmutants in enumerate(mypath)]
    # create the R table string and scripts
    # get the R table
    table_string = RUtil.get_table_string(arr, headers)
    # get the R script
    script = get_ggplot()
    # create the R plot image
    device_name = Form.g_imageformat_to_r_function[fs.imageformat]
    retcode, r_out, r_err, image_data = RUtil.run_plotter(
            table_string, script, device_name)
    if retcode:
        raise RUtil.RError(r_err)
    return image_data

def get_r_tikz_prop_plot(nsels, time_stats):
    """
    @param time_stats: a list of stats for each time point
    @return: tikz code corresponding to an R plot
    """
    out = StringIO()
    time_stats_trans = zip(*time_stats)
    y_low = 0
    y_high = 1
    ylim = RUtil.mk_call_str("c", y_low, y_high)
    print >> out, RUtil.mk_call_str(
        "plot",
        "my.table$t",
        "my.table$prop.mi.diag.approx",
        type='"n"',
        ylim=ylim,
        xlab='"time"',
        ylab='"proportion"',
        main='"proportion of same sign difference as MI"',
    )
    colors = ("red", "orange", "green", "blue", "black")
    plot_indices = (12, 13, 14, 15, 16)
    for c, plot_index in zip(colors, plot_indices):
        header = g_time_stats_headers[plot_index]
        print >> out, RUtil.mk_call_str("lines", "my.table$t", "my.table$%s" % header, col='"%s"' % c)
    return out.getvalue()

def get_response_content(fs):
    f_info = ctmcmi.get_mutual_info_known_distn
    # define the R table headers
    headers = ['log.probability.ratio', 'mutual.information']
    # make the array
    arr = []
    for x in np.linspace(fs.x_min, fs.x_max, 101):
        row = [x]
        proc = evozoo.AlternatingHypercube_d_1(3)
        X = np.array([x])
        distn = proc.get_distn(X)
        Q = proc.get_rate_matrix(X)
        info = f_info(Q, distn, fs.t)
        row.append(info)
        arr.append(row)
    # create the R table string and scripts
    # get the R table
    table_string = RUtil.get_table_string(arr, headers)
    # get the R script
    script = get_ggplot()
    # create the R plot image
    device_name = Form.g_imageformat_to_r_function[fs.imageformat]
    retcode, r_out, r_err, image_data = RUtil.run_plotter(
            table_string, script, device_name)
    if retcode:
        raise RUtil.RError(r_err)
    return image_data

def get_response_content(fs):
    # validate and store user input
    if fs.x_max <= fs.x_min:
        raise ValueError('check the min and max logs')
    f_info = divtime.get_fisher_info_known_distn_fast
    # define the R table headers
    headers = ['log.probability.ratio', 'fisher.information']
    # make the array
    arr = []
    for x in np.linspace(fs.x_min, fs.x_max, 101):
        row = [x]
        proc = evozoo.DistinguishedCornerPairHypercube_d_1(3)
        X = np.array([x])
        distn = proc.get_distn(X)
        Q = proc.get_rate_matrix(X)
        info = f_info(Q, distn, fs.t)
        row.append(info)
        arr.append(row)
    # create the R table string and scripts
    # get the R table
    table_string = RUtil.get_table_string(arr, headers)
    # get the R script
    script = get_ggplot()
    # create the R plot image
    device_name = Form.g_imageformat_to_r_function[fs.imageformat]
    retcode, r_out, r_err, image_data = RUtil.run_plotter(
            table_string, script, device_name)
    if retcode:
        raise RUtil.RError(r_err)
    return image_data

def get_response_content(fs):
    # precompute some transition matrices
    P_drift_selection = pgmsinglesite.create_drift_selection_transition_matrix(
            fs.npop, fs.selection_ratio)
    MatrixUtil.assert_transition_matrix(P_drift_selection)
    P_mutation = pgmsinglesite.create_mutation_transition_matrix(
            fs.npop, fs.mutation_ab, fs.mutation_ba)
    MatrixUtil.assert_transition_matrix(P_mutation)
    # define the R table headers
    headers = [
            'generation',
            'number.of.mutants',
            'probability',
            'log.prob',
            ]
    # compute the transition matrix
    P = np.dot(P_drift_selection, P_mutation)
    # Compute the endpoint conditional probabilities for various states
    # along the unobserved path.
    nstates = fs.npop + 1
    M = np.zeros((nstates, fs.ngenerations))
    M[fs.nmutants_initial, 0] = 1.0
    M[fs.nmutants_final, fs.ngenerations-1] = 1.0
    for i in range(fs.ngenerations-2):
        A_exponent = i + 1
        B_exponent = fs.ngenerations - 1 - A_exponent
        A = np.linalg.matrix_power(P, A_exponent)
        B = np.linalg.matrix_power(P, B_exponent)
        weights = np.zeros(nstates)
        for k in range(nstates):
            weights[k] = A[fs.nmutants_initial, k] * B[k, fs.nmutants_final]
        weights /= np.sum(weights)
        for k, p in enumerate(weights):
            M[k, i+1] = p
    arr = []
    for g in range(fs.ngenerations):
        for k in range(nstates):
            p = M[k, g]
            if p:
                logp = math.log(p)
            else:
                logp = float('-inf')
            row = [g, k, p, logp]
            arr.append(row)
    # create the R table string and scripts
    # get the R table
    table_string = RUtil.get_table_string(arr, headers)
    # get the R script
    script = get_ggplot()
    # create the R plot image
    device_name = Form.g_imageformat_to_r_function[fs.imageformat]
    retcode, r_out, r_err, image_data = RUtil.run_plotter(
            table_string, script, device_name)
    if retcode:
        raise RUtil.RError(r_err)
    return image_data

def get_response_content(fs):
    f_info = divtime.get_fisher_info_known_distn_fast
    requested_triples = []
    for triple in g_process_triples:
        name, desc, zoo_obj = triple
        if getattr(fs, name):
            requested_triples.append(triple)
    if not requested_triples:
        raise ValueError('nothing to plot')
    # define the R table headers
    r_names = [a.replace('_', '.') for a, b, c in requested_triples]
    headers = ['t'] + r_names
    # Spend a lot of time doing the optimizations
    # to construct the points for the R table.
    arr = []
    for t in cbreaker.throttled(
            progrid.gen_binary(fs.start_time, fs.stop_time),
            nseconds=5, ncount=200):
        row = [t]
        for python_name, desc, zoo_class in requested_triples:
            zoo_obj = zoo_class(fs.d)
            df = zoo_obj.get_df()
            opt_dep = OptDep(zoo_obj, t, f_info)
            if df:
                X0 = np.random.randn(df)
                xopt = scipy.optimize.fmin(
                        opt_dep, X0, maxiter=10000, maxfun=10000)
                # I would like to use scipy.optimize.minimize
                # except that this requires a newer version of
                # scipy than is packaged for ubuntu right now.
                # fmin_bfgs seems to have problems sometimes
                # either hanging or maxiter=10K is too big.
                """
                xopt = scipy.optimize.fmin_bfgs(opt_dep, X0,
                        gtol=1e-8, maxiter=10000)
                """
            else:
                xopt = np.array([])
            info_value = -opt_dep(xopt)
            row.append(info_value)
        arr.append(row)
    arr.sort()
    npoints = len(arr)
    # create the R table string and scripts
    # get the R table
    table_string = RUtil.get_table_string(arr, headers)
    # get the R script
    script = get_ggplot()
    # create the R plot image
    device_name = Form.g_imageformat_to_r_function[fs.imageformat]
    retcode, r_out, r_err, image_data = RUtil.run_plotter(
            table_string, script, device_name)
    if retcode:
        raise RUtil.RError(r_err)
    return image_data

def hard_coded_analysis():
    branch_length = 5.0
    sequence_length = 1000
    nsequences = 1000
    estimate_triple_list = []
    column_headers = ('most.info', 'less.info', 'least.info')
    for i in range(nsequences):
        # sample sequence changes at three levels of informativeness
        sequence_changes = sample_sequence_changes(
                branch_length, sequence_length)
        # get a distance estimate for each level of informativeness
        estimate_triple = sample_distance(*sequence_changes)
        estimate_triple_list.append(estimate_triple)
    print RUtil.get_table_string(estimate_triple_list, column_headers)

def get_table_string_and_scripts(start_stop_pairs, nsamples):
    """
    Command-line only.
    """
    # build the array for the R table
    data_arr = []
    sequence_lengths = []
    midpoints = []
    for start_pos, stop_pos in start_stop_pairs:
        sequence_length = stop_pos - start_pos + 1
        means, variations, covs = get_value_lists(
                start_pos, stop_pos, nsamples)
        midpoint = (start_pos + stop_pos) / 2.0
        row = [sequence_length, midpoint]
        for values in means, variations, covs:
            corr_info = mcmc.Correlation()
            corr_info.analyze(values)
            hpd_low, hpd_high = mcmc.get_hpd_interval(0.95, values)
            row.extend([hpd_low, corr_info.mean, hpd_high])
        data_arr.append(row)
        sequence_lengths.append(sequence_length)
        midpoints.append(midpoint)
    # build the table string
    table_string = RUtil.get_table_string(data_arr, g_headers)
    # get the scripts
    scripts = get_ggplot2_scripts(nsamples, sequence_lengths, midpoints)
    # return the table string and scripts
    return table_string, scripts

def get_r_tikz_stub():
    user_script = RUtil.g_stub
    device_name = "tikz"
    retcode, r_out, r_err, tikz_code = RUtil.run_plotter_no_table(user_script, device_name)
    if retcode:
        raise RUtil.RError(r_err)
    return tikz_code

def get_response_content(fs):
    # define some fixed values
    N_diploid = 6
    N_hap = 2 * N_diploid
    plot_density = 8
    # define some mutation rates
    theta_values = [0.001, 0.01, 0.1, 1.0]
    # define some selection coefficients to plot
    Ns_low = 0.0
    Ns_high = 3.0
    Ns_values = np.linspace(Ns_low, Ns_high, 3 * plot_density + 1)
    # get the values for each h
    Nr_values = (0, 5)
    arr_0 = get_plot_array(N_diploid, Nr_values[0], theta_values, Ns_values)
    arr_1 = get_plot_array(N_diploid, Nr_values[1], theta_values, Ns_values)
    ylab = '"expected returns to AB"'
    # define x and y plot limits
    xlim = (Ns_low, Ns_high)
    ylim = (np.min((arr_0, arr_1)), np.max((arr_0, arr_1)))
    ylogstr = '""'
    # http://sphaerula.com/legacy/R/multiplePlotFigure.html
    out = StringIO()
    print >> out, mk_call_str("par", mfrow="c(1,2)", oma="c(0,0,2,0)")
    print >> out, get_plot("left", Nr_values[0], arr_0, theta_values, Ns_values, xlim, ylim, ylogstr, ylab)
    print >> out, get_plot("right", Nr_values[1], arr_1, theta_values, Ns_values, xlim, ylim, ylogstr, '""')
    print >> out, mk_call_str("title", '"expected number of returns to AB, 2N=%s"' % N_hap, outer="TRUE")
    script = out.getvalue().rstrip()
    # create the R plot image
    device_name = Form.g_imageformat_to_r_function[fs.imageformat]
    retcode, r_out, r_err, image_data = RUtil.run_plotter_no_table(script, device_name)
    if retcode:
        raise RUtil.RError(r_err)
    return image_data