text:
# -*- coding: utf-8 -*-
"""
This file analyze the random samples from the enzyme-constrained models
created by Eduard Kerkhoven by using the Raven Toolbox.

This file is more a collection of snippets used to create figures for the manuscript than anything else. It is poorly commented / documented.

Author: Snorre Sulheim
Created: 12.03.2019
Modified: 30.09.2019
Email: snorre.sulheim@sintef.no
"""
import pandas as pd
from pathlib import Path
import cobra
from matplotlib import pyplot as plt
import matplotlib
import numpy as np
import seaborn as sns

pd.set_option("display.max_rows", 150)
pd.set_option("display.max_columns", 101)

# MAT
matplotlib.rcParams.update({'font.size': 10, 'legend.loc':'upper right'})

# Definition of different groups of pathways used om plotting
SELECTED_PATHWAYS_M145 = ["Glycolysis/Gluconeogenesis", "Pyruvate metabolism", "Oxidative phosphorylation", "Glycine, serine and threonine metabolism",
                          "Citric Acid Cycle", "Alanine, aspartate and glutamate metabolism", "Pentose phosphate pathway",
                          "Fatty acid biosynthesis", "Valine, leucine and isoleucine degradation", "Nucleotide biosynthesis",
                          "Undecylprodigiosin Biosynthesis", "Calcium-Dependent Antibiotics Biosynthesis", "Actinorhodin Biosynthesis", "Coelimycin biosynthesis"]#,
                          # "Fatty acid biosynthesis", "Glycerophospholipid metabolism"]

SELECTED_PATHWAYS_M1152 = ["Glycolysis/Gluconeogenesis", "Pyruvate metabolism", "Oxidative phosphorylation", "Glycine, serine and threonine metabolism",
                          "Citric Acid Cycle", "Alanine, aspartate and glutamate metabolism", "Pentose phosphate pathway", "Glyoxylate and dicarboxylate metabolism",
                          "Fatty acid biosynthesis", "Valine, leucine and isoleucine degradation", "Nucleotide biosynthesis"]

NUCLEOTIDE_METABOLISM = ["Purine metabolism", "Pyrimidine metabolism"]

CENTRAL_CARBON =  ["Glycolysis/Gluconeogenesis",
                  "Citric Acid Cycle",
                  "Pentose phosphate pathway",
                  "Fructose and mannose metabolism",
                  "Alternate Carbon metabolism",
                  "Pyruvate metabolism",
                  "Pentose and glucuronate interconversions",
                  "Glycogen metabolism",
                  "Peptidoglycan biosynthesis",
                  "Murein Recycling",
                  "Murein Biosynthesis",
                  "Glyoxylate and dicarboxylate metabolism",
                  "Starch and sucrose metabolism",
                  "Butanoate metabolism",
                  "Propanoate metabolism",
                  "Carbon fixation pathways in prokaryotes",
                  "Amino sugar and nucleotide sugar metabolism",
                  "C5-Branched dibasic acid metabolism",
                  "Inositol phosphate metabolism",
                ]

AMINO_ACIDS = ["Alanine, aspartate and glutamate metabolism",
               "Glycine, serine and threonine metabolism",
               "Valine, leucine and isoleucine degradation",
               "Phenylalanine, tyrosine and tryptophan biosynthesis",
               "Valine, leucine and isoleucine biosynthesis",
               "Arginine and proline metabolism",
               "Cysteine and methionine metabolism",
               "Lysine degradation",
               "Arginine biosynthesis",
               "Cysteine metabolism",
               "Lysine biosynthesis",
               "Histidine metabolism",
               "Phenylalanine metabolism",
               "Methionine metabolism",
               "Thiamine metabolism",
               "Tyrosine metabolism",
               "Tryptophan metabolism",
               "D-Alanine metabolism",
               "Ascorbate and aldarate metabolism",
               "D-Glutamine and D-glutamate metabolism",
               "Nitrogen metabolism"]

LIPIDS =  ["Fatty acid biosynthesis",
          "Glycerolipid metabolism",
          "Fatty acid degradation",
          "Glycerophospholipid metabolism"]

VITAMINS_COFACTORS =  ["Vitamin B6 metabolism",
                      "Pantothenate and CoA biosynthesis",
                      "beta-Alanine metabolism",
                      "Nicotinate and nicotinamide metabolism",
                      "Folate biosynthesis",
                      "Porphyrin and chlorophyll metabolism",
                      "Riboflavin metabolism",
                      "Folate metabolism",
                      "One carbon pool by folate",
                      "Lipoic acid metabolism",
                      "Biotin metabolism"]

TOXIC_COMPOUND_DEGRADATION = ["Atrazine degradation",
                              "Benzoate degradation",
                              "Aminobenzoate degradation",
                              ]

ALL_BGC =  ["Terpenoid backbone biosynthesis",
           "Desferrioxamine biosynthesis",
           "Germicidin Biosynthesis",
           "Sesquiterpenoid and triterpenoid biosynthesis",
           "Ubiquinone and other terpenoid-quinone biosynthesis",
           "Albaflavenol biosynthesis",
           "Geosmin synthase",
           "Biflaviolin synthase",
           "2-methylisoborneol synthase"]

REMOVED_BGC = ["Undecylprodigiosin Biosynthesis",
               "Actinorhodin Biosynthesis",
               "Calcium-Dependent Antibiotics Biosynthesis",
               "Coelimycin biosynthesis"]

OXIDATIVE_STRESS = ["Mycothiol metabolism",
                    "Catalase",
                    "Superoxide dismutase"]

IGNORE = ["D-Arginine and D-ornithine metabolism",
          "Toluene degradation",
          "Polyketide sugar unit biosynthesis",
          "Galactose metabolism",
          "Biosynthesis of type II polyketide backbone",
          "Glutathione metabolism"
          ]

OTHER = ["Nucleotide biosynthesis",
         "Methane metabolism",
         "Sulfur metabolism",
         "Oxidative phosphorylation",
         "NADH repair",
         "tRNA Charging",
         ]

def pathway_analysis_plot(model_fn, random_samples_fn, selected_pwys, row_order = None, labels = None, mask_rows = None, sep = "\t", key = "pathway", absolute_values = True, folder = None):
    """
    Plot random sampling heatmap for both M145 and M1152 in idividual plots
    """
    if not folder:
      folder = "."

    model = cobra.io.read_sbml_model(model_fn)
    df = get_random_samples(random_samples_fn, model, key = key, sep = sep, column_key = "MEAN", absolute_values = absolute_values)
    sub_df = df.groupby(key).sum()

    # Remove all zero rows
    sub_df = sub_df.loc[(sub_df != 0).any(axis = 1), :]

    # Standardize
    standardized_df = sub_df.subtract(sub_df.mean(axis = 1), axis = 0).divide(sub_df.std(axis = 1), axis = 0)

    #M145
    M145_columns = [x for x in standardized_df.columns if "M145" in x]
    M145_sub_df = standardized_df.loc[selected_pwys, M145_columns]
    M145_sub_df.index.name = None
    # M145_sub_df.columns = [x.replace("_", " - ")+"h" for x in M145_sub_df.columns.values]
    M145_sub_df.columns = [x.split("_")[-1] for x in M145_sub_df.columns.values]

    # M1152
    M1152_columns = [x for x in standardized_df.columns if "M1152" in x]
    M1152_sub_df = standardized_df.loc[selected_pwys, M1152_columns]
    M1152_sub_df.index.name = None
    # M1152_sub_df.columns = [x.replace("_", " - ")+"h" for x in M1152_sub_df.columns.values]
    M1152_sub_df.columns = [x.split("_")[-1] for x in M1152_sub_df.columns.values]

    if mask_rows:
        M1152_sub_df.iloc[mask_rows, :] = None

    if absolute_values:
        abs_string = "abs"
    else:
        abs_string = "real"

    # Plot Settings
    # fig, [ax_M145, ax_M1152] = plt.subplots(1, 2, sharey = True)
    figsize = (14, 14)
    cmap = sns.cm.vlag
    cmap.set_bad("gray", alpha = 0)

    g_M145 = sns.clustermap(M145_sub_df, cmap = cmap, col_cluster = False, vmin = -1.3, vmax = 1.3, figsize = figsize)
    print("Row order:", g_M145.dendrogram_row.reordered_ind)
    # plt.setp(g_M145.ax_heatmap.xaxis.get_majorticklabels(), rotation=90)
    plt.subplots_adjust(left = 0.03, right = 0.5)

    path_M145 = folder + "/all_pathways_uptake_M145_{0}.svg".format(abs_string)
    plt.savefig(path_M145)
    # plt.savefig("C:/Users/snorres/OneDrive - SINTEF/SINTEF projects/INBioPharm/scoGEM/random sampling/Eduard random sampling/M145/all_pathways.png")
    # plt.show()
    plt.close()

    if row_order:
        M1152_sub_df = M1152_sub_df.iloc[g_M145.dendrogram_row.reordered_ind, :]
    # Set hatches
    g_M1152 = sns.clustermap(M1152_sub_df, cmap = cmap, col_cluster = False, row_cluster = False, vmin = -1.3, vmax = 1.3, figsize = figsize)
    g_M1152.ax_heatmap.patch.set(hatch='//', edgecolor='black')
    # plt.setp(g_M1152.ax_heatmap.xaxis.get_majorticklabels(), rotation=45)
    plt.subplots_adjust(left = 0.03, right = 0.5)
    # plt.savefig("C:/Users/snorres/OneDrive - SINTEF/SINTEF projects/INBioPharm/scoGEM/random sampling/Eduard random sampling/M1152/all_pathways.png")
    path_M1152 = folder + "/all_pathways_uptake_M1152_{0}.eps".format(abs_string)
    plt.savefig(path_M1152)
    # plt.show()
    plt.close()

def pathway_analysis(model_fn, random_samples_fn, selected_pwys, strain = "M145", row_order = None, labels = None,
                     mask_rows = None, store = True, key = "pathway", sep = ",", absolute_values = True,
                     show_cbar = False, mask_cells = None):
    """
    Plot pathway heatmap for either M145, M1152 or both in one and same plot
    """
    model = cobra.io.read_sbml_model(model_fn)
    df = get_random_samples(random_samples_fn, model, key = key, sep = sep, column_key = "MEAN", absolute_values = absolute_values)
    sub_df = df.groupby(key).sum()

    # Remove all zero rows
    sub_df = sub_df.loc[(sub_df.abs() > 1e-8).any(axis = 1), :]
    print(sub_df.loc[(sub_df.abs() < 1e-8).any(axis = 1), :].abs().max(axis = 1).sort_values(ascending = False))

    # Standardize
    standardized_df = sub_df.subtract(sub_df.mean(axis = 1), axis = 0).divide(sub_df.std(axis = 1), axis = 0)

    # print(sub_df.mean(axis = 1).sort_values(ascending = False))
    std_df = sub_df.iloc[:, 1:].std(axis = 1).sort_values(ascending = False)
    # rstd_df = sub_df.std(axis = 1).divide(sub_df.mean(axis = 1), axis = 0).sort_values(ascending = False)
    print(std_df)
    # print(rstd_df)

    if store:
        store_df = sub_df.copy()
        store_df["std"] = sub_df.std(axis = 1)
        norm_type = random_samples_fn.rsplit("_")[-1].split(".")[0]
        csv_fn = "C:/Users/snorres/OneDrive - SINTEF/SINTEF projects/INBioPharm/scoGEM/random sampling/{0}/sub_df_{1}.svg".format(strain, norm_type)
        store_df.to_csv(csv_fn)

    strain_columns = [x for x in standardized_df.columns if strain in x]
    strain_df = standardized_df.loc[:, strain_columns]

    # Select top X
    # selected_pwys = list(std_df.index)[:70]
    if not selected_pwys:
        selected_pwys = strain_df.index
    selected_df = strain_df.loc[selected_pwys, :]
    selected_df.index.names = ["Pathways"]

    if labels:
        print(labels)
        selected_df = selected_df.loc[:, labels]

    # if isinstance(mask_cells, np.ndarray):
    #     print(selected_df.shape, mask_cells.shape)
    #     selected_df[mask_cells] = None

    selected_df.columns = ["-".join(x.split("_")[1:]) for x in selected_df.columns.values]

    # Settings
    figsize = (10, len(selected_pwys)/3)
    cmap = sns.cm.vlag
    cmap.set_bad("gray", alpha = 0)
    print(selected_df)
    if row_order:
        selected_df = selected_df.iloc[row_order, :]
        g = sns.clustermap(selected_df, cmap = cmap, col_cluster = False, row_cluster = False, figsize = figsize, vmin = -1.2, vmax = 1.2, yticklabels = True)
    else:
        g = sns.clustermap(selected_df, cmap = cmap, col_cluster = False, figsize = figsize, vmin = -1.2, vmax = 1.2, yticklabels = True, mask=mask_cells)
        print("Row order:", g.dendrogram_row.reordered_ind)

    # Remove labels
    ax = g.ax_heatmap
    ax.set_xlabel("")
    ax.set_ylabel("")

    # Show colorbar
    g.cax.set_visible(show_cbar)

    # Set hatches
    g.ax_heatmap.patch.set(hatch='//', edgecolor='black')
    # .xticks(rotation=45)
    #plt.setp(g.ax_heatmap.xaxis.get_majorticklabels(), rotation=45, ha = "center")
    plt.subplots_adjust(left = 0.03, right = 0.60, bottom = 0.27, top = 0.97)
    plt.show()

def get_random_samples(fn, model, key = "subsystem", sep = "\t", column_key = "MEAN", skipcolumns = 0, absolute_values = True):
    """
    Reads the random samples and performs some data preprocessing.
    """
    df = pd.read_csv(fn, index_col = 0, header = 0, sep = sep)
    df_data = df.iloc[:, skipcolumns:]

    # A few reactions are changed sign of to make the total flux even out with a correlated "looped" reaction
    # PGM and PGK is just defined in the opposite direction in the model

    df_data.loc[["SUCOAS", "PGM", "PGK"], :] = df_data.loc[["SUCOAS", "PGM", "PGK"], :].abs()
    df_data.loc["HACD1", :] *=-1

    if absolute_values:
        df_data = df_data.abs()
    if column_key:
        columns = [x for x in df.columns.values if column_key in x]
        df_data = df_data.loc[:, columns]

    subsystem_df = add_subsystem_to_df(df_data, model, key)
    return subsystem_df

def get_random_samples_raw(fn, sep = "\t", column_key = "MEAN", skipcolumns = 0):
    """
    Reads the random sampling data without the preprocessing
    """
    df = pd.read_csv(fn, index_col = 0, header = 0, sep = sep)
    df_data = df.iloc[:, skipcolumns:]

    if column_key:
        columns = [x for x in df.columns.values if column_key in x]
        df_data = df_data.loc[:, columns]
    return df_data

def get_random_samples_and_split(fn, model, key = "subsystem"):
    df = pd.read_csv(fn, index_col = 0, header = 0, sep = "\t")
    df_data = df.iloc[:, 2:].abs()
    subsystem_df = add_subsystem_to_df(df_data, model, key)

    M1152_columns = [x for x in df.columns.values if "MEAN_M1152" in x] + [key]
    M1152_df = subsystem_df[M1152_columns]
    M145_columns = [x for x in df.columns.values if "MEAN_M145" in x] + [key]
    M145_df = subsystem_df[M145_columns]
    return M145_df, M1152_df

def add_subsystem_to_df(df, model, key = "subsystem"):
    """
    df is a dataframe where the indexes are reaction ids and the columns are different timepoints / strains
    the values are the different fluxes
    """
    df[key] = None
    keep_columns = list(df.columns.values)
    add_rows = []
    df["Keep"] = True
    for r_id, row in df.iterrows():
        try:
            r = model.reactions.get_by_id(r_id)
        except KeyError:
            df.loc[r_id, "Keep"] = False
            # print(r_id)
            continue
        try:
            subsystem = r.annotation[key]
        except KeyError:
            df.loc[r_id, "Keep"] = False
            print("Missing subsystem: {0}".format(r_id))
            continue

        if isinstance(subsystem, str) and len(subsystem):
            if "Valine, leucine and isoleucine" in subsystem:

                # Currently there are Valine, leucine and isoleucine metabolism, biosynthesis and degradation,
                # but it makes sense to combine these
                # NOTE FOR REVIEW: Undid this merging
                # subsystem = "Valine, leucine and isoleucine metabolism"
                pass
            # elif subsystem in OTHER_AMINO_ACIDS:
            #     subsystem = "Metabolism of other amino acids"
            elif subsystem == "FERI metabolism":
                # FERI metabolism is only one reaction, FNOR, which recycles NADP with ferredoxin used by AKGDH2 reaction
                subsystem = "Citric Acid Cycle"
            elif subsystem in NUCLEOTIDE_METABOLISM:
                subsystem = "Nucleotide biosynthesis"
            elif r_id == "CDAS11":
                subsystem = "Fatty acid degradation"
            elif r.id == "MALTHIK":
                subsystem = "Glyoxylate and dicarboxylate metabolism"
            elif r.id == "HEX1":
                subsystem = "Glycolysis/Gluconeogenesis"
            elif r.id == "HACD1":
                subsystem = "Butanoate metabolism"

            df.loc[r_id, key] = subsystem
            # print(r_id, subsystem)
        else:
            print("Multiple subsystem annotations for: ", r.id)
            df.loc[r_id, "Keep"] = False

    df = df.loc[df["Keep"], :]
    df_to_keep = df.loc[:, keep_columns]
    # df_to_keep.loc[:, keep_columns[:-1]] = df_to_keep.loc[:, keep_columns[:-1]].abs()/df_to_keep.loc[:, keep_columns[:-1]].abs().sum()
    return df_to_keep

def make_contribution_subsystems_plots(model_fn, random_samples_fn, sep = "\t", key = "pathway", absolute_values = False):
    """
    Used for data QC. Plots histogram to see which reactions that have the largest contribution to each pathway
    """
    model = cobra.io.read_sbml_model(model_fn)
    df = get_random_samples(random_samples_fn, model, key = key, sep = sep, column_key = "MEAN", absolute_values = absolute_values)

    if absolute_values:
        abs_string = "abs"
    else:
        abs_string = "real"

    for pathway, pwy_df in df.groupby(key):
        fig, ax = plt.subplots(1)
        pwy_df.drop(key, axis = 1, inplace = True)
        # pwy_df = pwy_df.abs()
        pwy_df.plot(kind = "bar", logy = absolute_values, title = pathway, figsize = (20, 14), ax = ax).legend(loc = 1)
        ax.set_ylim(1.1*pwy_df.min().min(), 1.1*pwy_df.max().max())
        pathway = pathway.replace("/", "-")
        fn = "C:/Users/snorres/OneDrive - SINTEF/SINTEF projects/INBioPharm/scoGEM/random sampling/Random sampling july/pathway contribution plot/{0}_{1}.svg".format(pathway, abs_string)
        fig.savefig(fn)
        plt.close()

def analyze_accoa_consumption(model_fn, random_samples_fn, sep = "\t", subplot = True):
    metabolite_id = "accoa_c"
    model = cobra.io.read_sbml_model(model_fn)
    df = get_random_samples_raw(random_samples_fn, sep = sep, column_key = "MEAN")

    metabolite = model.metabolites.get_by_id(metabolite_id)

    reaction_ids = []
    reaction_multiplier = []
    for r in metabolite.reactions:
        reaction_ids.append(r.id)
        reaction_multiplier.append(r.get_coefficient(metabolite_id))

    df_selected = (df.loc[reaction_ids, :].T*reaction_multiplier).T
    strain_annotation = [x.split("_")[1] for x in df_selected.columns.values]
    timepoint_annotation = [1,2,3,4,5,6,7,8,9,1,2,3,4,5,6,7,8]
    print(len(strain_annotation), len(timepoint_annotation))

    # Sum ACTCOALIG anf PTAr
    print(df_selected)

    df_all_0 = (df_selected.abs() < 1e-4).all(axis = 1)
    df_selected = df_selected.loc[~df_all_0, :].T

    # ACTOALIG and PTAr is creating a loop
    # df_selected["ACTCOALIG - PTAr"] = df_selected["ACTCOALIG"] + df_selected["PTAr"]
    # df_selected.drop(["PTAr", "ACTCOALIG"], axis = 1, inplace = True)

    df_selected["Strain"] = strain_annotation
    df_selected["Time point"] = timepoint_annotation
    df_long_form = pd.melt(df_selected, id_vars = ["Strain", "Time point"], value_vars = df_selected.columns.values[:-2], value_name = "Flux")

    # Production
    df_production = df_long_form.loc[df_long_form["Flux"] > 1e-8, :]
    df_consumption = df_long_form.loc[df_long_form["Flux"] < -1e-8, :]
    df_consumption["Flux"] *= -1

    if subplot:
        fig, [ax1, ax2] = plt.subplots(2)
    else:
        fig1, ax1 = plt.subplots(1)
        fig2, ax2 = plt.subplots(1)

    sns.lineplot(data = df_production, x = "Time point", y = "Flux", hue = "rxns", style = "Strain", ax = ax1)
    ax1.set_title("Acetyl-CoA production")
    ax1.set_yscale("log")

    sns.lineplot(data = df_consumption, x = "Time point", y = "Flux", hue = "rxns", style = "Strain", ax = ax2)
    ax2.set_title("Acetyl-CoA consumption")
    ax2.set_yscale("log")
    plt.show()

def analyze_malcoa_consumption(model_fn, random_samples_fn, sep = "\t", subplot = True):
    """
    Make two plots: One of the reactions producing malonyl-CoA and on for the reactions consuming Malonyl-CoA
    """
    metabolite_id = "malcoa_c"
    model = cobra.io.read_sbml_model(model_fn)
    df = get_random_samples_raw(random_samples_fn, sep = sep, column_key = "MEAN")

    metabolite = model.metabolites.get_by_id(metabolite_id)

    reaction_ids = []
    reaction_multiplier = []
    for r in metabolite.reactions:
        reaction_ids.append(r.id)
        reaction_multiplier.append(r.get_coefficient(metabolite_id))

    df_selected = (df.loc[reaction_ids, :].T*reaction_multiplier).T
    strain_annotation = [x.split("_")[1] for x in df_selected.columns.values]
    timepoint_annotation = [1,2,3,4,5,6,7,8,9,1,2,3,4,5,6,7,8]
    print(len(strain_annotation), len(timepoint_annotation))

    df_all_0 = (df_selected.abs() < 1e-10).all(axis = 1)
    df_selected = df_selected.loc[~df_all_0, :].T

    # Sum RED, ACT
    red_columns = [x for x in df_selected.columns if x[:4] == "REDS"]

    df_selected["Red"] = df_selected.loc[:, red_columns].sum(axis = 1)
    df_selected.drop(red_columns, axis = 1, inplace = True)

    cpk_columns = [x for x in df_selected.columns if x[:4] == "CPKS"]
    df_selected["Cpk"] = df_selected.loc[:, cpk_columns].sum(axis = 1)

    df_selected.drop(cpk_columns, axis = 1, inplace = True)

    cda_columns = [x for x in df_selected.columns if x[:4] == "CDAS"]
    df_selected["CDA"] = df_selected.loc[:, cda_columns].sum(axis = 1)
    df_selected.drop(cda_columns, axis = 1, inplace = True)

    act_columns = [x for x in df_selected.columns if x[:4] == "ACTS"]
    df_selected["Act"] = df_selected.loc[:, act_columns].sum(axis = 1)
    df_selected.drop(act_columns, axis = 1, inplace = True)

    # SUM ACCOAC and MCOATA
    # df_selected["ACCOAC + MCOATA + ACCOAC_1"] = df_selected["ACCOAC"] + df_selected["MCOATA"] + df_selected["ACCOAC"]
    # df_selected.drop(["ACCOAC", "MCOATA"], axis = 1, inplace = True)

    df_selected["Strain"] = strain_annotation
    df_selected["Time point"] = timepoint_annotation
    print(df_selected)
    df_selected = df_selected.reindex(["ACCOAC", "ACCOAC_1", "Act", "MCOATA", "Cpk", "Red", "CDA", "THYDNAPS", "Strain", "Time point"], axis = 1)
    df_long_form = pd.melt(df_selected, id_vars = ["Strain", "Time point"], value_vars = df_selected.columns.values[:-2], value_name = "Normalized flux")

    # Production
    df_production = df_long_form.loc[df_long_form["Normalized flux"] > 1e-8, :]
    df_consumption = df_long_form.loc[df_long_form["Normalized flux"] < -1e-8, :]
    df_consumption["Normalized flux"] *= -1
    print(df_consumption)

    if subplot:
        fig, [ax1, ax2] = plt.subplots(1, 2, sharey = True)
    else:
        fig1, ax1 = plt.subplots(1)
        fig2, ax2 = plt.subplots(1)

    sns.lineplot(data = df_production, x = "Time point", y = "Normalized flux", hue = "rxns", style = "Strain", ax = ax1, lw = 2)
    ax1.set_title("Malonyl-CoA production")
    ax1.set_yscale("log")
    ax1.legend(loc='lower left')#, bbox_to_anchor=(1.01, 0.5), ncol=1)
    # ax1.spines['right'].set_visible(False)
    # ax1.spines['top'].set_visible(False)

    sns.lineplot(data = df_consumption, x = "Time point", y = "Normalized flux", hue = "rxns", style = "Strain", ax = ax2, lw = 2)
    ax2.set_title("Malonyl-CoA consumption")
    ax2.set_yscale("log")
    ax2.legend(loc='lower left')
    # ax2.spines['right'].set_visible(False)
    # ax2.spines['top'].set_visible(False)
    sns.despine()

    plt.show()

    print(df_consumption[df_consumption["rxns"]!="MCOATA"].groupby(["Strain", "Time point"]).sum())
    print(df_consumption[df_consumption["rxns"]=="MCOATA"].groupby(["Strain", "Time point"]).sum())

def analyze_met_consumption(model_fn, random_samples_fn, metabolite_id = "nadh_c", metabolite_name = "NADH", sep = "\t", lim = 1e-3):
    """
    Make two plots: One of the reactions producing the chosen metabolite and one for the reactions consuming the metabolite
    """
    model = cobra.io.read_sbml_model(model_fn)
    df = get_random_samples_raw(random_samples_fn, sep = sep, column_key = "MEAN")

    metabolite = model.metabolites.get_by_id(metabolite_id)

    reaction_ids = []
    reaction_multiplier = []
    for r in metabolite.reactions:
        reaction_ids.append(r.id)
        reaction_multiplier.append(r.get_coefficient(metabolite_id))

    df_selected = (df.loc[reaction_ids, :].T*reaction_multiplier).T
    strain_annotation = [x.split("_")[1] for x in df_selected.columns.values]
    timepoint_annotation = [1,2,3,4,5,6,7,8,9,1,2,3,4,5,6,7,8]
    print(strain_annotation, timepoint_annotation)

    df_all_0 = (df_selected.abs() < lim).all(axis = 1)
    df_selected = df_selected.loc[~df_all_0, :].T

    df_selected["Strain"] = strain_annotation
    df_selected["Time point"] = timepoint_annotation
    df_long_form = pd.melt(df_selected, id_vars = ["Strain", "Time point"], value_vars = df_selected.columns.values[:-2], value_name = "Flux")

    print(df_long_form)

    # Production
    df_production = df_long_form[df_long_form["Flux"] > 1e-8]
    df_consumption = df_long_form[df_long_form["Flux"] < -1e-8]
    df_consumption["Flux"] *= -1

    print(df_production.groupby(["Strain", "Time point"]).sum())

    fig, [ax1, ax2] = plt.subplots(2)

    sns.lineplot(data = df_production, x = "Time point", y = "Flux", hue = "rxns", style = "Strain", ax = ax1, style_order = ["M145", "M1152"])
    ax1.set_title("{0} production".format(metabolite_name))
    ax1.set_yscale("log")

    sns.lineplot(data = df_consumption, x = "Time point", y = "Flux", hue = "rxns", style = "Strain", ax = ax2, style_order = ["M145", "M1152"])
    ax2.set_title("{0} consumption".format(metabolite_name))
    ax2.set_yscale("log")
    plt.show()

def plot_key_metabolic_reactions(random_samples_fn, title = "", sep = ","):
    """
    Plot the mean predicted flux through the following reactions selected by Tjasa Kumelj
    - PFK
    - G6PDH2r
    - TAGS140
    - PDH
    - MCOATA
    - ACTS1
    - ICDHyr
    - AKGDH2
    - CYO2a
    - CYO2b
    - NADH17b
    - GLUN
    - GLUDxi
    """
    key_reactions = ["TAGS140","MCOATA","PFK", "PDH","ACTS1","AKGDH2","ICDHyr","CYO2a","CYO2b","G6PDH2r","NADH17b","GLUN","GLUDxi"]
    df = get_random_samples_raw(random_samples_fn, sep = sep, column_key = "MEAN")
    df_selected = df.loc[key_reactions, :]
    df_selected.columns = ["-".join(x.split("_")[1:]) for x in df_selected.columns.values]
    # print(df_selected)
    # Standardize
    print(df_selected.subtract(df_selected.mean(axis = 1), 0))
    df_centered = df_selected.subtract(df_selected.mean(axis = 1), 0)
    df_z_score = df_centered.div(df_centered.std(axis = 1), 0)
    # strain_annotation = [x.split("_")[1] for x in df_selected.columns.values]
    # timepoint_annotation = [1,2,3,4,5,6,7,8,9,1,2,3,4,5,6,7,8]
    cmap = sns.cm.vlag
    # g = sns.clustermap(df_selected, row_cluster = False, col_cluster = False, z_score = 0, cmap = cmap, vmin = -1.5, vmax = 1.5)
    sns.heatmap(df_z_score, cmap = cmap, vmin = -1.5, vmax = 1.5, cbar_kws={'label': 'Flux Z-score'})
    plt.title(title)
    plt.subplots_adjust(bottom = 0.2)
    save_key = random_samples_fn.split("/")[-1].split(".")[0]
    plt.savefig("C:/Users/snorres/OneDrive - SINTEF/SINTEF projects/INBioPharm/scoGEM/random sampling/Random sampling july/key_metabolic_reactions_{0}.svg".format(save_key))
    # plt.show()
    plt.close()

def plot_selected_reactions(random_samples_fn, reactions, title = "", sep = ",", strain = "M145"):
    """
    Plot the flux for a list of chosen reactions
    """
    df = get_random_samples_raw(random_samples_fn, sep = sep, column_key = "MEAN")
    df_selected = df.loc[reactions, :]
    df_selected.columns = ["-".join(x.split("_")[1:]) for x in df_selected.columns.values]
    if strain:
        selected_columns = [x for x in df_selected.columns if strain in x]
        df_selected = df_selected.loc[:, selected_columns]
        df_selected.columns = [x.split("-")[1] for x in df_selected.columns]

    print(df_selected)

    ax = df_selected.T.plot(alpha = 0.9, logy = False)
    ax.set_ylabel("Normalized mean flux")
    ax.set_xlabel("Hours")
    plt.show()

def print_genes_in_subsystem(model_fn, pathway):
    model = cobra.io.read_sbml_model(model_fn)

    pathway_genes = []
    for r in model.reactions:
        try:
            p_r = r.annotation["pathway"]
        except:
            continue
        if p_r == pathway:
            pathway_genes += [g.id for g in r.genes]

    print("Genes related to {0};".format(pathway))
    print(list(set(pathway_genes)))

def plot_all_reactions_for_metabolite(model_fn, metabolite_id, random_samples_fn, sep = "\t"):
    model = cobra.io.read_sbml_model(model_fn)
    df = get_random_samples_raw(random_samples_fn, sep = sep, column_key = "MEAN")

    metabolite = model.metabolites.get_by_id(metabolite_id)
    reactions = [r.id for r in metabolite.reactions]

    coeff = [r.metabolites[metabolite] for r in metabolite.reactions]
    df_scaled = df.loc[reactions, :].multiply(coeff, axis = 0)
    # print(df_scaled)

    std_df = df_scaled.std(axis = 1).sort_values(ascending = False)
    print(std_df)
    # print(df_scaled[reactions, :].head())
    df_selected = df_scaled.loc[std_df > 1e-3, :]

    # Row colors
    lut = dict(zip([1,-1], "gr"))
    row_sign = [np.sign(r.metabolites[metabolite]) for r in metabolite.reactions if r.id in df_selected.index]
    row_colors = [lut[x] for x in row_sign]

    g = sns.clustermap(df_selected, z_score = None, cmap = sns.cm.vlag, col_cluster = False, vmin = -0.1, vmax = 0.1)#row_colors = row_colors)
    plt.show()

if __name__ == '__main__':
    co2_normalized_random_samples_fn = "C:/Users/snorres/Google Drive/scoGEM community model/Supporting information/Model/randomsampling_july/ec-RandSampComb_proteomics_CO2norm.tsv"
    gluglc_normalized_random_samples_fn = "C:/Users/snorres/Google Drive/scoGEM community model/Supporting information/Model/randomsampling_july/ec-RandSampComb_proteomics_GlcGluNorm.tsv"
    growth_normalized_random_samples_fn = "C:/Users/snorres/Google Drive/scoGEM community model/Supporting information/Model/randomsampling_july/ec-RandSampComb_proteomics_growthNorm.tsv"
    all_flux_normalized_random_samples_fn = "C:/Users/snorres/Google Drive/scoGEM community model/Supporting information/Model/randomsampling_july/ec-RandSampComb_proteomics_AllFluxNorm.tsv"

    model_fn = "../../ModelFiles/xml/Sco-GEM.xml"
    if 0:
        labels = ["M145_21", "M145_29", "M145_33", "M145_37", "M145_41", "M145_45", "M145_49", "M145_53", "M145_57"]  #
        subsystem_analysis(model_fn, co2_normalized_random_samples_fn, strain = "M145", labels = labels)
    if 0:
        labels = ["M1152_33", "M1152_41", "M1152_45", "M1152_49", "M1152_53", "M1152_57", "M1152_61", "M1152_65"] #
        # row_order = [2, 8, 10, 6, 13, 4, 14, 1, 9, 5, 7, 0, 12, 3, 11]
        # row_order = [6, 13, 2, 10, 0, 5, 7, 9, 14, 1, 4, 8, 11, 3, 12]
        row_order = None
        subsystem_analysis(model_fn, co2_normalized_random_samples_fn, strain = "M1152", labels = labels, row_order = row_order)

    if 0:
        labels = ["M145_21", "M145_29", "M145_33", "M145_37", "M145_41", "M145_45", "M145_49", "M145_53", "M145_57"]  #,
        subsystem_analysis(model_fn, uptake_normalized_random_samples_fn, strain = "M145", labels = labels)
    if 0:
        labels = ["M1152_33", "M1152_41", "M1152_45", "M1152_49", "M1152_53", "M1152_57", "M1152_61", "M1152_65"] #
        # row_order = [3, 0, 14, 4, 1, 9, 5, 7, 8, 12, 6, 13, 2, 10, 11] # all x33
        row_order = [3, 12, 8, 0, 14, 4, 1, 7, 5, 9, 6, 13, 2, 10, 11]
        subsystem_analysis(model_fn, uptake_normalized_random_samples_fn, strain = "M1152", labels = labels, row_order = row_order)

    if 0:
        matplotlib.rcParams.update({'font.size': 6, 'legend.loc':'upper right'})
        row_order = [4, 58, 19, 40, 26, 66, 23, 24, 69, 72, 46, 79, 18, 54, 14, 27, 7, 5, 67, 53, 29, 21, 74, 51, 48, 8, 52, 42, 49, 63, 70, 59, 81, 76, 44, 47, 61, 22, 30, 35, 56, 64, 80, 32, 45, 38, 65, 62, 43, 13, 28, 9, 71, 73, 31, 57, 36, 37, 15, 1, 78, 0, 12, 34, 17, 77, 20, 6, 11, 2, 55, 16, 25, 60, 75, 3, 10, 41, 50, 39, 33, 68]
        # subplots adjust left = 0.03, right = 0.35, bottom = 0.27, top = 0.97
        pathway_analysis(model_fn, co2_normalized_random_samples_fn, SELECTED_PATHWAYS_M145, strain = "MEAN", sep = "\t", store = False, absolute_values = False)
        #pathway_analysis(model_fn, gluglc_normalized_random_samples_fn, SELECTED_PATHWAYS_M145, strain = "MEAN", sep = "\t", store = False, absolute_values = False, row_order = row_order)

        # row_order = [7, 6, 2, 9, 12, 3, 10, 13, 11, 0, 1, 5, 4, 8]
        # mask_rows = [10,11,12,13]
        # pathway_analysis(model_fn, co2_normalized_random_samples_fn, SELECTED_PATHWAYS_M145, strain = "M1152", row_order = row_order, mask_rows = mask_rows)

    if 1:

        matplotlib.rcParams.update({'font.size': 16, 'legend.loc':'upper right'})
        row_order = [12, 13, 10, 11, 3, 2, 8, 6, 7, 9, 5, 4, 0, 1]
        mask_rows = [10,11,12,13]
        pathway_analysis_plot(model_fn, co2_normalized_random_samples_fn, SELECTED_PATHWAYS_M145, row_order = True, mask_rows = mask_rows, absolute_values = False)
        # pathway_analysis_plot(model_fn, uptake_normalized_random_samples_fn, SELECTED_PATHWAYS_M145, row_order = row_order, mask_rows = mask_rows)
    if 0:
        analyze_accoa_consumption(model_fn, co2_normalized_random_samples_fn, subplot = False)

    if 0:
        analyze_malcoa_consumption(model_fn, co2_normalized_random_samples_fn, subplot = True)

    if 0:
        analyze_pyruvate_consumption(model_fn, co2_normalized_random_samples_fn)
    if 0:
        analyze_glutamate_consumption(model_fn, co2_normalized_random_samples_fn)
    if 0:
        analyze_met_consumption(model_fn, co2_normalized_random_samples_fn, "akg_c", "Glutamate", lim = 1e-2)
    if 0:
        make_contribution_subsystems_plots(model_fn, co2_normalized_random_samples_fn, absolute_values = False)
    if 0:
        plot_key_metabolic_reactions(co2_normalized_random_samples_fn, "CO2-normalized ", "\t")
        plot_key_metabolic_reactions(gluglc_normalized_random_samples_fn, "Glucose and glutamate carbon uptake normalized","\t")
        plot_key_metabolic_reactions(growth_normalized_random_samples_fn, "Growth rate normalized", "\t")
        plot_key_metabolic_reactions(all_flux_normalized_random_samples_fn, "Normalized by sum of fluxes", "\t")

    if 0:
        # plot_selected_reactions(co2_normalized_random_samples_fn, ["CS", "ICDHyr"], "Aceyl-CoA consumption into TCA cycle", "\t")
        # plot_selected_reactions(co2_normalized_random_samples_fn, ["FBA", "PFK", "ENO", "G6PDH2r"], "Aceyl-CoA consumption into TCA cycle", "\t", strain = "M145")
        plot_selected_reactions(co2_normalized_random_samples_fn, ["ILETA", "LEUTA", "VALTA", "ILEDHr", "VALDHr", "LLEUDr"], "Branched-chain amino acids", "\t", strain = "M145")
    if 0:
        # Plot Supplementary figure 3 as several subpanels
        matplotlib.rcParams.update({'font.size': 10, 'legend.loc':'upper right'})

        # Glycolysis
        #pathway_analysis(model_fn, co2_normalized_random_samples_fn, CENTRAL_CARBON, strain = "MEAN", sep = "\t", store = False, absolute_values = False)

        # Amino acids
        # pathway_analysis(model_fn, co2_normalized_random_samples_fn, AMINO_ACIDS, strain = "MEAN", sep = "\t", store = False, absolute_values = False)

        # for key in [OXIDATIVE_STRESS]:
        #     # Lipids
        #     pathway_analysis(model_fn, co2_normalized_random_samples_fn, key, strain = "MEAN", sep = "\t", store = False, absolute_values = False, show_cbar = False)

        # BGC
        mask_cells = np.zeros((len(REMOVED_BGC+ALL_BGC), 17), dtype = bool)
        mask_cells[0:4, 9:] = 1
        # mask_cells = mask_cells.T
        pathway_analysis(model_fn, co2_normalized_random_samples_fn, REMOVED_BGC + ALL_BGC, strain = "MEAN", sep = "\t", store = False, absolute_values = False, mask_cells = mask_cells)
    if 0:
        # print_genes_in_subsystem(model_fn, "Alanine, aspartate and glutamate metabolism")
        print_genes_in_subsystem(model_fn, "D-Glutamine and D-glutamate metabolism")
    if 0:
        plot_all_reactions_for_metabolite(model_fn, "glu__L_c", co2_normalized_random_samples_fn)

max_stars_repo_name: edkerk/scoGEM

max_issues_repo_name: edkerk/scoGEM

max_forks_repo_name: edkerk/scoGEM


text:
import tensorflow as tf
import numpy as np
from pyrr import Quaternion
import os

from tensorflow_graphics.geometry.transformation import quaternion

from models.dense_net import relu, conv_costum_init_layer
from models.diff_renderer import MyDiffRendererStages
from models.dense_net import UDenseNet

FLAGS = tf.app.flags.FLAGS

crop_dimensions = {
    "x_crop_dim": np.array([-76, 76]),
    "y_crop_dim": np.array([-76, 76])
}

image_params = {
    "img_size": [640, 480],
}

dense_net_params = {
    "grow_rate": 16,
}

model_params = {
    "grid_size": 38,
}

model_names = ['ape','cam','cat','duck','glue','iron','phone',
               'benchvise','can','driller','eggbox','holepuncher','lamp']

scaling_2_mm = 1000.0
far_near_cut = [10000.0, 100.0]

R_init = np.array([
    [1.0, 0.0, 0.0],
    [0.0, 1.0, 0.0],
    [0.0, 0.0, 1.0],
])

def predict_the_pose(pos_init, quat_init, xy_pred, z_pred, quat_pred, img_scaling, cam_mat):

    img_center = tf.matmul(tf.expand_dims(pos_init, axis=1), cam_mat, transpose_b=True)
    img_center = tf.reduce_sum(img_center, axis=1)
    img_center = tf.div(img_center, tf.expand_dims(img_center[:, -1], axis=1))

    xy_image = img_center[:, :2] + xy_pred / tf.expand_dims(img_scaling, axis=-1)

    z_model = pos_init[:, 2]

    x_model_cam_frame = (xy_image[:, 0] - cam_mat[:, 0, 2]) / cam_mat[:, 0, 0]
    x_model_cam_frame = tf.expand_dims(tf.multiply(x_model_cam_frame, z_model), axis=-1)
    y_model_cam_frame = (xy_image[:, 1] - cam_mat[:, 1, 2]) / cam_mat[:, 1, 1]
    y_model_cam_frame = tf.expand_dims(tf.multiply(y_model_cam_frame, z_model), axis=-1)

    z_model = tf.expand_dims(z_model, axis=-1)
    pose_model_lateral = tf.concat((x_model_cam_frame, y_model_cam_frame, z_model), axis=-1)

    z_out = z_pred / 3000.0
    pos_out = tf.multiply(tf.exp(z_out), pose_model_lateral)

    quat = quat_pred / 10.0
    quat_len = tf.expand_dims(tf.norm(quat, axis=-1), axis=-1)
    quat_out = tf.div(quat, quat_len)
    quat_out = quaternion.multiply(quat_out, quat_init)

    return pos_out, quat_out

class ModelUDenseNetAttentionStages(tf.keras.layers.Layer):

    def __init__(self,  mode, weight_decay, dense_net_params, size_params, mesh_folder):
        super(ModelUDenseNetAttentionStages, self).__init__()

        if mode == tf.estimator.ModeKeys.TRAIN:
            training = True
        else:
            training = False
        self.training = training

        self.stages = []
        for it_stage in range(4):
            stage = ModelUDenseNetAttentionLayer(weight_decay, dense_net_params, size_params, mesh_folder, self.training)
            self.stages.append(stage)

    def call(self, input):

        img_input = input[0]
        quat_pos_ind_cam_mat = input[1]

        pose_prediction = []
        quat_prediction = []

        quat_init, pos_init, obj_ind, cam_mat = tf.split(quat_pos_ind_cam_mat, (4, 3, 1, 9), axis=-1)
        cam_mat_3x3 = tf.reshape(cam_mat, [-1, 3, 3])
        quat_pos_ind_cam_mat_pos_init = tf.concat((quat_init, pos_init, obj_ind, cam_mat, pos_init, quat_init), axis=-1)

        for it_stage in self.stages:

            quat, pos, obj_ind, cam_mat, pos_init, quat_init = tf.split(quat_pos_ind_cam_mat_pos_init,
                                                                        (4, 3, 1, 9, 3, 4), axis=-1)

            xy_out, z_out, quat_out, im_center, image_scaling = it_stage([img_input, quat_pos_ind_cam_mat_pos_init])
            pos_out, quat_out = predict_the_pose(pos, quat, xy_out, z_out, quat_out, image_scaling, cam_mat_3x3)

            pose_prediction.append(pos_out)
            quat_prediction.append(quat_out)

            quat_pos_ind_cam_mat_pos_init = tf.concat((quat_out, pos_out, obj_ind, cam_mat, pos_init, quat_init), axis=-1)

        return pose_prediction, quat_prediction

class ModelUDenseNetAttentionLayer(tf.keras.layers.Layer):

    def __init__(self, weight_decay, dense_net_params, size_params, mesh_folder, training):
        super(ModelUDenseNetAttentionLayer, self).__init__()

        self.batch_size = size_params["batch_size"]
        self.grid_size = size_params["grid_size"]
        weight_decay = (weight_decay, 0)
        dense_net_params["weight_decay"] = weight_decay

        self.training = training

        models_path = []
        for it_model in model_names:
            models_path.append(os.path.join(mesh_folder, it_model + '_color_4000.ply'))

        self.diff_rend = MyDiffRendererStages(far_near_cut, models_path, R_init, scaling_2_mm, use_texture=True)
        self.concat_join_in = tf.keras.layers.Concatenate(axis=-1)

        self.u_dense_net = UDenseNet('dense_net', dense_net_params, training)
        self.batch_norm_layer = tf.keras.layers.BatchNormalization(name='dense_net_out_bn')

        # xy stream
        self.xy_conv_layer = conv_costum_init_layer(256, 1, 'quat_pos_conv_layer', weight_decay,
                                                    1.0 / self.u_dense_net.features_count)
        self.xy_att = tf.keras.layers.Dense(256, name="quat_pos_att")
        self.xy_conv_out = tf.keras.layers.Dense(2, name="quat_pos_out")
        self.xy_out_batch_norm = tf.keras.layers.BatchNormalization(name='quat_pose_out_bn')

        # quat stream
        self.quat_conv_layer = conv_costum_init_layer(256, 1, 'quat_conv_layer', weight_decay,
                                                      1.0 / self.u_dense_net.features_count)
        self.quat_att = tf.keras.layers.Dense(256, name="quat_att")
        self.quat_out = tf.keras.layers.Dense(4, name="quat_out")
        self.quat_out_batch_norm = tf.keras.layers.BatchNormalization(name='quat_out_bn')

        # z stream
        self.z_conv_layer = conv_costum_init_layer(256, 1, 'z_conv_layer', weight_decay,
                                                   1.0 / self.u_dense_net.features_count)
        self.z_out_att = tf.keras.layers.Dense(256, name="z_out_att")
        self.z_out = tf.keras.layers.Dense(1, name="z_out")
        self.z_out_batch_norm = tf.keras.layers.BatchNormalization(name='z_out_bn')

        self.attention_out_batch_norm_layer = tf.keras.layers.BatchNormalization(name='attention_out_bn')

        # att xy
        self.attention_conv_layer_xy = conv_costum_init_layer(256, 1, 'attention_conv_layer', weight_decay,
                                                              1.0 / self.u_dense_net.features_count)
        self.attention_out_xy = conv_costum_init_layer(1, 1, 'attention_out', weight_decay,
                                                       1.0 / 256)
        self.attention_out_xy_batch_norm = tf.keras.layers.BatchNormalization(name='attention_out_bn')

        # att z
        self.attention_conv_layer_z = conv_costum_init_layer(256, 1, 'attention_conv_layer_z', weight_decay,
                                                             1.0 / self.u_dense_net.features_count)
        self.attention_out_z = conv_costum_init_layer(1, 1, 'attention_out_z', weight_decay,
                                                      1.0 / 256)
        self.attention_out_z_batch_norm = tf.keras.layers.BatchNormalization(name='attention_z_out_bn')

        # att quat
        self.attention_conv_layer_q = conv_costum_init_layer(256, 1, 'attention_conv_layer_quat', weight_decay,
                                                             1.0 / self.u_dense_net.features_count)
        self.attention_out_q = conv_costum_init_layer(1, 1, 'attention_out_quat', weight_decay,
                                                      1.0 / 256)
        self.attention_out_q_batch_norm = tf.keras.layers.BatchNormalization(name='attention_q_out_bn')

    def call(self, input):
        img_input = input[0]
        quat_pos_input = input[1]

        # diff renderer
        diff_renderer_out, cube_vertices_clip, im_center, image_scaling = self.diff_rend(quat_pos_input)
        full_image = tf.zeros((self.batch_size, 2, 2))
        bb_xy_bounds = tf.zeros((self.batch_size, 2, 2))
        diff_renderer_out_zero_mean = tf.subtract(diff_renderer_out, 0.5)

        # dense net 38x38 grid predictions
        joint_feature = self.concat_join_in([img_input, diff_renderer_out_zero_mean])
        ae_out, bottleneck_out, decoder_2_out, decoder_1_out = self.u_dense_net(joint_feature)

        x = self.batch_norm_layer(ae_out, training=self.training)
        x = relu(x)

        # attention out
        attention_bn = self.attention_out_batch_norm_layer(ae_out, training=self.training)
        attention = relu(attention_bn)
        attention = self.attention_conv_layer_xy(attention)
        attention = self.attention_out_xy_batch_norm(attention, training=self.training)
        attention = relu(attention)
        attention_xy = self.attention_out_xy(attention)

        attention = relu(attention_bn)
        attention = self.attention_conv_layer_z(attention)
        attention = self.attention_out_z_batch_norm(attention, training=self.training)
        attention = relu(attention)
        attention_z = self.attention_out_z(attention)

        attention = relu(attention_bn)
        attention = self.attention_conv_layer_q(attention)
        attention = self.attention_out_q_batch_norm(attention, training=self.training)
        attention = relu(attention)
        attention_quat = self.attention_out_q(attention)

        # z out
        attention_map_z = tf.reduce_sum(attention_z, axis=-1)
        attention_map_z = tf.reshape(attention_map_z, shape=[self.batch_size, self.grid_size * self.grid_size])
        attention_map_z = tf.nn.softmax(attention_map_z, axis=-1)
        attention_map_z = tf.reshape(attention_map_z, shape=[self.batch_size, self.grid_size, self.grid_size])
        attention_map_z = tf.expand_dims(attention_map_z, axis=-1)

        z_out = self.z_conv_layer(x)
        z_out = self.z_out_batch_norm(z_out, training=self.training)
        z_out = relu(z_out)
        z_out_weighted = tf.multiply(z_out, attention_map_z)
        z_out_sum = tf.reduce_sum(z_out_weighted, axis=[1, 2])
        z_out = self.z_out_att(z_out_sum)
        z_out = self.z_out(z_out)

        # xy out
        attention_map_xy = tf.reduce_sum(attention_xy, axis=-1)
        attention_map_xy = tf.reshape(attention_map_xy, shape=[self.batch_size, self.grid_size * self.grid_size])
        attention_map_xy = tf.nn.softmax(attention_map_xy, axis=-1)
        attention_map_xy = tf.reshape(attention_map_xy, shape=[self.batch_size, self.grid_size, self.grid_size])
        attention_map_xy = tf.expand_dims(attention_map_xy, axis=-1)

        xy_out = self.xy_conv_layer(x)
        xy_out = self.xy_out_batch_norm(xy_out, training=self.training)
        xy_out = relu(xy_out)
        xy_out_weighted = tf.multiply(xy_out, attention_map_xy)
        xy_out_sum = tf.reduce_sum(xy_out_weighted, axis=[1, 2])
        xy_out = self.xy_att(xy_out_sum)
        xy_out = self.xy_conv_out(xy_out)

        #  quat out
        attention_map_quat = tf.reduce_sum(attention_quat, axis=-1)
        attention_map_quat = tf.reshape(attention_map_quat, shape=[self.batch_size, self.grid_size * self.grid_size])
        attention_map_quat = tf.nn.softmax(attention_map_quat, axis=-1)
        attention_map_quat = tf.reshape(attention_map_quat, shape=[self.batch_size, self.grid_size, self.grid_size])
        attention_map_quat = tf.expand_dims(attention_map_quat, axis=-1)

        quat_out = self.quat_conv_layer(x)
        quat_out = self.quat_out_batch_norm(quat_out, training=self.training)
        quat_out = relu(quat_out)
        qaut_out_weighted = tf.multiply(quat_out, attention_map_quat)
        quat_out_sum = tf.reduce_sum(qaut_out_weighted, axis=[1, 2])
        quat_out = self.quat_att(quat_out_sum)
        quat_out = self.quat_out(quat_out)

        return xy_out, z_out, quat_out, im_center, image_scaling,

def loss_stages(labels, net_out, keypoints_3D_model):

    loss_list = []
    loss_sum = tf.constant(0.0)

    for it_pred_pose, it_pred_quat in zip(net_out['pos'], net_out['quat']) :

        loss_it = 0.5 * stage_loss_keypoints(it_pred_pose, it_pred_quat, labels['pos'], labels['quat'], keypoints_3D_model)
        loss_it = tf.reduce_mean(loss_it)
        loss_list.append(loss_it)
        loss_sum += loss_it

    loss = loss_sum/len(loss_list)

    return loss, loss_list

def stage_loss_keypoints(pos_out, quat_out, pos_gt, quat_gt, keypoints_3D_model):

    keypoints_3D = quaternion.rotate(keypoints_3D_model, tf.expand_dims(quat_gt, axis=-2))
    keypoints_3D_predicted = quaternion.rotate(keypoints_3D_model, tf.expand_dims(quat_out, axis=1))

    keypoints_3D = keypoints_3D + tf.expand_dims(pos_gt, axis=1)
    keypoints_3D_predicted = keypoints_3D_predicted + tf.expand_dims(pos_out, axis=1)

    keypoint_distance = keypoints_3D_predicted - keypoints_3D
    keypoint_distance = tf.norm(keypoint_distance, axis=-1)
    loss_keypoint = tf.reduce_mean(0.5 * tf.square(keypoint_distance), axis=-1)

    return loss_keypoint

def model_fn(features, labels, mode, hparams):

    weight_decay = hparams.weight_decay
    batch_size = tf.shape(features["pos_init"])[0]

    # init pos and quaternion
    pos_init = features["pos_init"]
    quat_init = features["quat_init"]
    obj_ind = features["object_ind"]
    keypoints_3D_model = tf.convert_to_tensor(hparams.keypoints_list)

    keypoints_3D_model_batch = tf.tile(tf.expand_dims(keypoints_3D_model, axis=0), [batch_size, 1, 1, 1])
    keypoints_3D_model = tf.gather_nd(keypoints_3D_model_batch, obj_ind, batch_dims=1)
    keypoints_3D_model = tf.cast(keypoints_3D_model, dtype=tf.float32)

    grid_size = model_params["grid_size"]
    cam_mat = features["cam_mat"]

    quat_pos_indx_cam_mat_input = tf.concat([quat_init, pos_init, tf.cast(obj_ind, dtype=tf.float32),
                                             tf.reshape(cam_mat, [-1, 3*3])], axis=-1)

    size_params = {
        "batch_size": batch_size,
        "grid_size": grid_size,
    }

    model = ModelUDenseNetAttentionStages(mode, weight_decay, dense_net_params, size_params, hparams.mesh_folder)

    img_input = features['img']
    img_input_scaled = tf.subtract(tf.div(tf.cast(img_input, dtype=tf.float32), 256.0), 0.5) # this is ok, cropping does not scale values

    ind = tf.cast(features["object_ind"], dtype=tf.float32)
    ind = tf.expand_dims(ind, axis=1)
    ind = tf.expand_dims(ind, axis=1)
    x_img_dim = crop_dimensions["x_crop_dim"][1] - crop_dimensions["x_crop_dim"][0]
    y_img_dim = crop_dimensions["y_crop_dim"][1] - crop_dimensions["y_crop_dim"][0]
    cls_img_in = tf.ones((batch_size, y_img_dim, x_img_dim, 1)) * ind / 13.0 - 0.5

    output_pos, output_quat = model([img_input_scaled, quat_pos_indx_cam_mat_input, cls_img_in, cam_mat])

    if mode == tf.estimator.ModeKeys.PREDICT:
        labels = {
            "pos": features["pos"],
            "quat": features["quat"],
        }

    net_out = {
        "pos": output_pos,
        "quat": output_quat,
    }

    loss, loss_list = loss_stages(labels, net_out, keypoints_3D_model)

    # summaries
    tf.summary.scalar('loss_stage_1', loss_list[0])
    tf.summary.scalar('loss_stage_2', loss_list[1])
    tf.summary.scalar('loss_stage_3', loss_list[2])
    tf.summary.scalar('loss_stage_4', loss_list[3])

    return {
        "loss": loss,
        "predictions": {
            "pos_out_0": output_pos[0],
            "pos_out_1": output_pos[1],
            "pos_out_2": output_pos[2],
            "pos_out_3": output_pos[3],
            "quat_out_0": output_quat[0],
            "quat_out_1": output_quat[1],
            "quat_out_2": output_quat[2],
            "quat_out_3": output_quat[3],
            "input": features["img"],
            "pos_init": features["pos_init"],
            "quat_init": features["quat_init"],
            "pos_gt": features["pos"],
            "quat_gt": features["quat"],
            "cam_mat": cam_mat,
        },
        "eval_metric_ops": {
        },
        "model_keras": model,
    }

def estimator_model_fn(features, labels, mode, params):
    ret = model_fn(features, labels, mode, params)
    train_op = None
    training_hooks = []

    return tf.estimator.EstimatorSpec(
        mode=mode,
        predictions=ret["predictions"],
        loss=ret["loss"],
        train_op=train_op,
        eval_metric_ops=ret["eval_metric_ops"],
        training_hooks=training_hooks)

def predict_eval_and_save_external(steps, hparams, sync_replicas, input_fn, external_options=None):

    object_to_test = external_options["obj_name"]
    out_folder = external_options["out_folder"]
    ckpt_path = external_options["ckpt_path"]
    dataset_folder = external_options["dataset_folder"]
    mesh_folder = external_options["mesh_folder"]

    run_config = tf.estimator.RunConfig(
        keep_checkpoint_every_n_hours=0.5,
        save_checkpoints_secs=180,
        save_summary_steps=50)

    estimator = tf.estimator.Estimator(
        model_fn=estimator_model_fn,
        params=hparams, config=run_config)

    print("loading model: ", ckpt_path)

    input_fn = input_fn(dataset=object_to_test, folder_tfrecord=dataset_folder, mesh_folder=mesh_folder, batch_size=1)
    predictions = estimator.predict(input_fn=input_fn, checkpoint_path=ckpt_path)

    if not os.path.exists(out_folder):
        os.mkdir(out_folder)
    out_folder_object = os.path.join(out_folder, object_to_test)
    if not os.path.exists(out_folder_object):
        os.mkdir(out_folder_object)

    it_img = 0

    for item in predictions:

        pos_gt = item["pos_gt"]
        quat_gt = item["quat_gt"]

        pos_init = item["pos_init"]
        quat_init = item["quat_init"]

        pred_0 = item["pos_out_0"]
        pred_1 = item["pos_out_1"]
        pred_2 = item["pos_out_2"]
        pred_3 = item["pos_out_3"]

        quat_0 = item["quat_out_0"]
        quat_1 = item["quat_out_1"]
        quat_2 = item["quat_out_2"]
        quat_3 = item["quat_out_3"]

        pose_pred_it = np.zeros((3, 4))
        pose_pred_it[:, 3] = pred_3
        quat_pred = Quaternion(quat_3)
        R_pred = np.array(quat_pred.matrix33)
        pose_pred_it[:3, :3] = R_pred

        pose_pred_it_0 = np.zeros((3, 4))
        quat_pred = Quaternion(quat_0)
        R_pred = np.array(quat_pred.matrix33)
        pose_pred_it_0[:3, :3] = R_pred
        pose_pred_it_0[:, 3] = pred_0

        pose_pred_it_1 = np.zeros((3, 4))
        quat_pred = Quaternion(quat_1)
        R_pred = np.array(quat_pred.matrix33)
        pose_pred_it_1[:3, :3] = R_pred
        pose_pred_it_1[:, 3] = pred_1

        pose_pred_it_2 = np.zeros((3, 4))
        quat_pred = Quaternion(quat_2)
        R_pred = np.array(quat_pred.matrix33)
        pose_pred_it_2[:3, :3] = R_pred
        pose_pred_it_2[:, 3] = pred_2

        pose_gt_it = np.zeros((3, 4))
        pose_gt_it[:, 3] = pos_gt
        quat_gt = Quaternion(quat_gt)
        R_gt = np.array(quat_gt.matrix33)
        pose_gt_it[:3, :3] = R_gt

        if it_img % 100 == 0:
            print(it_img)

        object_to_test = "%06d" % it_img
        out_file = os.path.join(out_folder_object, object_to_test + "_predict.npy")
        out_dict = {
            'pose_pred': pose_pred_it,
            'pose_gt': pose_gt_it,
            'pose_pred_0': pose_pred_it_0,
            'pose_pred_1': pose_pred_it_1,
            'pose_pred_2': pose_pred_it_2,
        }
        np.save(out_file, out_dict)

        it_img += 1

ext: py

reason: import numpy


text:
"""
checkBonds.py

Brief:  Check that all D atoms are bonded to one
other D atom in d144.

Author:  Ian S. Dunn, Lawrence Livermore National Laboratory
"""

# Import modules.
import sys
import numpy as np

# Load .xyz movie file.
with open(sys.argv[1], 'r') as f:
    lines = f.readlines()

ang2bohr=1.8897269
bohr2ang=1./ang2bohr

box = np.array([42.4813, 21.2406, 21.2406])  * bohr2ang

current_timestep = []
num_timesteps = 0

for i in range(int(sys.argv[2]), len(lines)):

    if 'H' in lines[i]:
        coords = lines[i].rstrip('\n').split()[1:4]
        coords = np.array([float(coord) for coord in coords])
        current_timestep.append(coords)

    elif ' ' not in lines[i] and lines[i] != '\n':
        if len(current_timestep) == 0:
            continue

        current_timestep = np.array(current_timestep)

        # Check for points with no neighbors.
        neighborless = 0

        for k in range(len(current_timestep)):

            smallest_distance = 100000000

            for j in range(len(current_timestep)):

                if j == k:
                    continue

                distance = current_timestep[k] - current_timestep[j]

                for l in range(3):
                    distance[l] = min(abs(distance[l]), abs(abs(distance[l]) - box[l]))

                distance = np.linalg.norm(distance)

                if distance < smallest_distance:
                    smallest_distance = distance

            if smallest_distance > 4.0:
                neighborless += 1
                print k

        print "Time step "+str(np.ceil(i / 218.0))+": "+str(neighborless)
        num_timesteps += 1
        current_timestep = []


text:
"""
Display two vectors layers ontop of a 4-D image layer. One of the vectors
layers is 3D and "sliced" with a different set of vectors appearing on
different 3D slices. Another is 2D and "broadcast" with the same vectors
apprearing on each slice.
"""

import numpy as np
from skimage import data
import napari

with napari.gui_qt():
    blobs = np.stack(
        [
            data.binary_blobs(
                length=128, blob_size_fraction=0.05, n_dim=3, volume_fraction=f
            )
            for f in np.linspace(0.05, 0.5, 10)
        ],
        axis=0,
    )
    viewer = napari.view(blobs.astype(float))

    # sample vector coord-like data
    n = 200
    pos = np.zeros((n, 2, 2), dtype=np.float32)
    phi_space = np.linspace(0, 4 * np.pi, n)
    radius_space = np.linspace(0, 20, n)

    # assign x-y position
    pos[:, 0, 0] = radius_space * np.cos(phi_space) + 64
    pos[:, 0, 1] = radius_space * np.sin(phi_space) + 64

    # assign x-y projection
    pos[:, 1, 0] = 2 * radius_space * np.cos(phi_space)
    pos[:, 1, 1] = 2 * radius_space * np.sin(phi_space)

    planes = np.round(np.linspace(0, 128, n)).astype(int)
    planes = np.concatenate(
        (planes.reshape((n, 1, 1)), np.zeros((n, 1, 1))), axis=1
    )
    vectors = np.concatenate((planes, pos), axis=2)

    # add the sliced vectors
    layer = viewer.add_vectors(
        vectors, edge_width=0.4, name='sliced vectors', edge_color='blue'
    )

    viewer.dims.ndisplay = 3

max_stars_repo_name: arokem/napari

max_issues_repo_name: arokem/napari

max_forks_repo_name: arokem/napari


text:
import sys

from time import sleep

import numpy as np

from scipy import interpolate

from beamngpy import BeamNGpy, Scenario, Road, Vehicle, setup_logging
import json

SIZE = 1024

def main():
    pristine_filename = 'C:/Users/merie/Documents/BeamNG.research/vehicles/pickup/pristine.save.json'
    filename = 'C:/Users/merie/Documents/BeamNG.research/vehicles/pickup/poppedtires.save.json'
    with open(filename, 'r') as f:
        file_json = f.read().replace('\n', '')
        j = json.loads(file_json)
        print(j)
        print()
        print(j.keys())
        for k in j.keys():
            try:
                print("j[{}].keys() = {}".format(k, j[k].keys()))
                for k2 in j[k].keys():
                    #if 'tire' in k2:
                    print("\tj[{}][{}] = {}".format(k, k2, j[k][k2]))
                        #print(j[k][k2])
                print()
            except:
                if 'luaState' in k:
                    print("j[{}] =".format(k))
                    parts = j[k].split(',')
                    for p in parts:
                        print("\t{}".format(p))
                    print()
                else:
                    print("j[{}] = {}\n".format(k, j[k]))

if __name__ == '__main__':
    main()

max_stars_repo_name: MissMeriel/BeamNGpy

max_issues_repo_name: MissMeriel/BeamNGpy

max_forks_repo_name: MissMeriel/BeamNGpy

reason: import numpy,from scipy


text:
# -*- coding: utf-8 -*-
"""Vehicle_Navigation-Shortest_Path.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/103unhaid5NTquk2vVaEhzRrYFdc7e2Y5
"""

!pip install osmnx==0.16.2

!pip install contextily

# Commented out IPython magic to ensure Python compatibility.
#Import necessary libraries.
import numpy as np
import math
import folium
import contextily as ctx
import shapely
import osmnx as ox
import geopandas as gpd
from pyproj import Transformer
from shapely.geometry import  Point, Polygon, LineString
from matplotlib.ticker import FormatStrFormatter
from matplotlib.patches import Circle
from osgeo import ogr
import matplotlib.pyplot as plt
# %matplotlib inline

def pointreproject(point_x,point_y,inp_CRS,out_CRS):

  """Converts from (latitude, longitude) to native map projection (x,y) coordinates and vice versa.

    Args:
        point_x: input CRS coordinate (x or latitude).
        point_y: input CRS coordinate (y or longitude).
        inp_CRS: reference system of input coordinates: 'epsg:4326','epsg:3857'...
        out_CRS: reference system of output coordinates: 'epsg:4326','epsg:3857'...

    Returns:
        A tuple (latitude, longitude) or (x,y) of the reprojected coordinates of the input point in the output CRS.

  """
  transformer = Transformer.from_crs(inp_CRS, out_CRS)
  out_x,out_y = transformer.transform(point_x, point_y)
  transformer = []

  return out_x,out_y

def rgraphfrompoints(starting_point,ending_point):

  """Creates a graph, with edges (roads) and nodes (road intersections), based on the two given points.
     The two given points must be in the WGS84 coordinate reference system and in the form of (latitude, longitude).

    Args:
        starting_point: the route's starting point. -> A tuple of (latitude, longitude).
        ending_point: the route's ending point. -> A tuple of (latitude, longitude).

    Returns:
        A reprojected (from EPSG:4326 to EPSG:3857) MultiDiGraph from the osmnx package.

  """
  #Input points geographical coordinates.
  lat1 = starting_point[0]
  lat2 = ending_point[0]
  long1 = starting_point[1]
  long2 = ending_point[1]

  #Find the function's argument values.
  north = lat2
  south = lat1
  if  lat1 > lat2:
    north = lat1
    south = lat2

  east = long2
  west = long1
  if  long1 > long2:
    east = long1
    west = long2

  #Construct the road graph.
  roadgraph = ox.graph_from_bbox(north+0.01, south-0.01, east+0.01, west-0.01, network_type = 'drive', simplify=False )

  #Reproject the road graph.
  projected_roadgraph = ox.projection.project_graph(roadgraph, to_crs='epsg:3857')

  return projected_roadgraph

def nearest_edge_node_id(point,roadgraph):

  """Finds the nearest network's node id for the given point, which lies in the nearest network edge.
     Input Point -> Nearest Edge Point -> Nearest Network Node (Route Starting/Ending Point).
     The given point's coordinates must be in the EPSG:4326.

    Args:
        point: The point:(latitude, longitude) that needs to be transfered to the nearest network node..
        roadgraph: A reprojected (from EPSG:4326 to EPSG:3857) MultiDiGraph from the osmnx package.

    Returns:
        The nearest network's node id for the given point, which is in the nearest network edge.

  """
  #Reproject point coordinates from EPSG:4326 to the EPSG:3857.
  point_x, point_y = pointreproject(point[0], point[1],'epsg:4326','epsg:3857')

  #Find the nearest road.
  _, _, _, nearest_road = ox.distance.get_nearest_edge(roadgraph, (point_y,point_x), return_geom=True)

  #Find all point coordinates that are constructing the specific road.
  nearest_road_coords = np.array(nearest_road.xy)

  #Calculate all distances between the input point and all the points constructing the nearest road line.
  e_distances = []
  for i in range(0,nearest_road_coords[0].size):
    e_dist = ox.distance.euclidean_dist_vec(point_y,point_x,nearest_road_coords[1][i],nearest_road_coords[0][i])
    e_distances.append(e_dist)
  mindist = min(e_distances)
  for i in range(0,len(e_distances)):
    if e_distances[i]==mindist: pos=i

  #Nearest edge node coordinates.
  nearest_edge_node_x = nearest_road_coords[0][pos]
  nearest_edge_node_y = nearest_road_coords[1][pos]

  #Find the nearest network's node id, which has the minimum distance from the road point that has the minimum distance from the input point.
  node_id= ox.distance.get_nearest_node(roadgraph, (nearest_edge_node_y,nearest_edge_node_x),  method='euclidean')

  return node_id

def pointreproject2edge(point,roadgraph):

  """
    Finds the nearest network's edge for the given point, and reprojects it perpendicularly to its nearest edge.
    If the orthogonal projection of the point lies outside the nearest edge's nodes, then the point isn't reprojected
    but it is transferred to the nearest network node instead.

     The given point's coordinates must be in the EPSG:4326, and following the form of (latitude, longitude).

    Args:
        point: The point:(latitude, longitude) that needs to be projected to the nearest network edge.
        roadgraph: A MultiDiGraph (EPSG:3857) from the osmnx package.

    Returns:
        A tuple (latitude, longitude) of the reprojected coordinates of the input point.

  """

  #Reproject point coordinates from EPSG:4326 to the EPSG:3857.
  point_x, point_y = pointreproject(point[0], point[1],'epsg:4326','epsg:3857')

  #Find the nearest road.
  _, _, _, nearest_road, dist__pp = ox.distance.get_nearest_edge(roadgraph, (point_y,point_x), return_geom=True, return_dist=True)

  #Extract the nearest road coordinates.
  nearest_road_coords = np.array(nearest_road.xy)

  #Calculate the euclidean distance between the input point and the nearest edge's first node.
  dist__pv = ox.distance.euclidean_dist_vec(point_y,point_x,nearest_road_coords[1][0],nearest_road_coords[0][0])

  #Calculate the hypotenuse.
  dist = np.sqrt(dist__pv**2-dist__pp**2)

  #Create the corresponding circle.
  circle = Point(nearest_road_coords[0][0], nearest_road_coords[1][0]).buffer(dist).boundary

  #Find the circle's intersection point with the given edge.
  intersection_points = circle.intersection(nearest_road)

  lies_in = False

  #If the orthogonal projection of the point lies inside the nearest edge's nodes.
  if type(intersection_points) == shapely.geometry.point.Point:
    lies_in = True
    intersection_point_x = intersection_points.x
    intersection_point_y = intersection_points.y

  #If the orthogonal projection of the point lies outside the nearest edge's nodes.
  else:
    gdf = ox.utils_graph.graph_to_gdfs(roadgraph)
    nodes_gdf = gdf[0].reset_index(drop=True)
    nearest_node_id = nearest_edge_node_id(point,roadgraph)
    intersection_point_x = nodes_gdf.loc[(nodes_gdf['osmid']==nearest_node_id)].x.values[0]
    intersection_point_y = nodes_gdf.loc[(nodes_gdf['osmid']==nearest_node_id)].y.values[0]

  #Find the nearest node id from the intersection point.
  node_id = ox.distance.get_nearest_node(roadgraph, (intersection_point_y,intersection_point_x),  method='euclidean')

  #Reproject the calculated grid point from epsg:3857 to epsg:4326.
  point_lat, point_long = pointreproject(intersection_point_x, intersection_point_y,'epsg:3857','epsg:4326')

  return point_lat, point_long, lies_in

def routeplot(starting_point,ending_point):

  """
     Finds and plots the sortest route -based on Dijkstra's pathfinding algorithm- between the
     two given points in the given roadgraph.The given points coordinates must be in the EPSG:4326
     and they must be following the form of (latitude, longitude).

    Args:
        starting_point: the route's starting point. -> A tuple of (latitude, longitude).
        ending_point: the route's ending point. -> A tuple of (latitude, longitude).

    Returns:
        A Folium map with the given starting and ending points and the shortest path between them.

  """

  #Calculate the roadgraph between the starting and ending points.
  roadgraph = rgraphfrompoints(starting_point,ending_point)

  #Find the routes starting end ending point which lies inside the given roadgraph.
  p_lat1,p_long1,lies_in1 = pointreproject2edge(starting_point,roadgraph)
  p_lat2,p_long2,lies_in2 = pointreproject2edge(ending_point,roadgraph)

  #Find the nearest edge node id.
  node1_id = nearest_edge_node_id(starting_point,roadgraph)
  node2_id = nearest_edge_node_id(ending_point,roadgraph)

  #Calculate the sortest path between the two given nodes.
  nodeids = ox.shortest_path(roadgraph, node1_id, node2_id)

  #Create the necessary geodataframes.
  gdf = ox.utils_graph.graph_to_gdfs(roadgraph)
  nodes_gdf = gdf[0].reset_index(drop=True)

  #Map's center point.
  center_lat = (starting_point[0]+ending_point[0])/2
  center_long = (starting_point[1]+ending_point[1])/2
  center = [center_lat,center_long]

  #Define the map.
  m = folium.Map(location=center, zoom_start=13, crs='EPSG3857')

  #Plot Original points.
  folium.Marker((starting_point[0], starting_point[1]),
                 popup="Original Starting Point", tooltip="Original Starting Point",
                 icon=folium.Icon(color='pink', icon="fas fa-car", prefix='fa')).add_to(m);

  folium.Marker((ending_point[0], ending_point[1]),
                 popup="Original Destination", tooltip="Original Destination",
                 icon=folium.Icon(color='pink', icon="fas fa-crosshairs", prefix='fa')).add_to(m);

  #Plot Starting and Destination points.
  folium.Marker((p_lat1,p_long1),
                popup="Starting Point", tooltip="Starting Point",
                icon=folium.Icon(color='darkblue', icon="fas fa-car", prefix='fa')).add_to(m);

  folium.Marker((p_lat2,p_long2),
                popup="Destination", tooltip="Destination",
                icon=folium.Icon(color='darkblue', icon="fas fa-crosshairs", prefix='fa')).add_to(m);

  #Create a list of tuples (locations) for each of route's node coordinates.
  locations = []
  locations.append((p_lat1,p_long1)) #Add the starting point to the route.
  for node_id in nodeids:
    p_lat = nodes_gdf.loc[(nodes_gdf['osmid']==node_id)].lat.values[0]
    p_long = nodes_gdf.loc[(nodes_gdf['osmid']==node_id)].lon.values[0]
    locations.append((p_lat,p_long))
  locations.append((p_lat2,p_long2)) #Add the ending point to the route.

  #Cut the final route segment if the starting point lies in it.
  if lies_in1:
    r1 = pointreproject(locations[1][0],locations[1][1],'epsg:4326','epsg:3857')
    r2 = pointreproject(locations[2][0],locations[2][1],'epsg:4326','epsg:3857')
    road1 = LineString([r1,r2])
    point1_x, point1_y = pointreproject(p_lat1,p_long1,'epsg:4326','epsg:3857')
    point1 = Point(point1_x, point1_y)
    if road1.distance(point1) < 1e-8 :del locations[1]

  #Cut the final route segment if the ending point lies in it.
  if lies_in2:
    r3 = pointreproject(locations[-3][0],locations[-3][1],'epsg:4326','epsg:3857')
    r4 = pointreproject(locations[-2][0],locations[-2][1],'epsg:4326','epsg:3857')
    road2 = LineString([r3,r4])
    point2_x, point2_y = pointreproject(p_lat2,p_long2,'epsg:4326','epsg:3857')
    point2 = Point(point2_x, point2_y)
    if road2.distance(point2) < 1e-8 :del locations[-2]

  #Create a folium polyline from the list of "locations".
  route = folium.PolyLine(locations=locations,tooltip='Route', smooth_factor=0.3, color="purple", weight=2.5, opacity=1)

  #Plot route.
  route.add_to(m);

  #Save map to html file.
  m.save('route.html') #saves the .html file in the local directory.

  #Plot map.
  return m

#Random points to WGS84 between latitude of 37.9000-38.0000 degrees and longitude of 23.7000-23.8000 degrees. (Athens, Attiki, Greece)
starting_point = (round(37.9+np.random.random()/10,4),round(23.7+np.random.random()/10,4))
ending_point = (round(37.9+np.random.random()/10,4),round(23.7+np.random.random()/10,4))

#Plot the route between starting and ending point.
routeplot(starting_point,ending_point)

#Print a .txt file with the code requirements.
!pip freeze > /content/drive/MyDrive/Projects/Python_Projects/Vehicle_Navigation-Sortest_Path/requirements.txt


text:
import numpy as np
import os
import pandas as pd
import pickle
from utils import *

class bern_emb_data():
    def __init__(self, cs, ns, n_minibatch, L):
        assert cs%2 == 0
        self.cs = cs
        self.ns = ns
        self.n_minibatch = n_minibatch
        self.L = L
        url = 'http://mattmahoney.net/dc/'
        filename = maybe_download(url, 'text8.zip', 31344016)
        words = read_data(filename)
        self.build_dataset(words)
        self.batch = self.batch_generator()
        self.N = len(self.data)

    def build_dataset(self, words):
        count = [['UNK', -1]]
        count.extend(collections.Counter(words).most_common(self.L - 1))
        dictionary = dict()
        for word, _ in count:
            dictionary[word] = len(dictionary)
        data = list()
        unk_count = 0
        for word in words:
            if word in dictionary:
              index = dictionary[word]
            else:
              index = 0
              unk_count += 1
            data.append(index)
        count[0][1] = unk_count
        reverse_dictionary = dict(zip(dictionary.values(), dictionary.keys()))
        self.data = np.array(data)
        self.count = count
        self.dictionary = dictionary
        self.labels = [reverse_dictionary[x] for x in range(len(reverse_dictionary))]
        unigram_dist = np.array([1.0*i for ii, i in count])
        unigram_dist = (unigram_dist/unigram_dist.sum())**(3.0/4)
        self.unigram = unigram_dist/unigram_dist.sum()
        with open('fits/vocab.tsv', 'w') as txt:
            for word in self.labels:
                txt.write(word+'\n')

    def batch_generator(self):
        batch_size = self.n_minibatch + self.cs
        data = self.data
        while True:
            if data.shape[0] < batch_size:
                data = np.hstack([data, self.data])
                if data.shape[0] < batch_size:
                    continue
            words = data[:batch_size]
            data = data[batch_size:]
            yield words

    def feed(self, placeholder):
        return {placeholder: self.batch.next()}


text:
import weldnumpy as wn
import numpy as np

def assert_correct(np_a, z):
    '''
    common part of the check.
    '''
    shape = []
    for s in z.shape:
        shape.append(s)
    shape = tuple(shape)

    np_a = np.reshape(np_a, shape)

    for i in range((z.shape[0])):
        for j in range(z.shape[1]):
            print('{}, {}'.format(i,j))
            assert np_a[i][j] == z[i][j]

def test_view():
    '''
    Adding the iter code to a view.
    In general, do we deal correctly with wn.array(view)??
    '''
    orig = np.random.rand(20,20)
    a = orig[3:15:1,4:20:2]

    # orig = np.random.rand(20,20,20)
    # a = orig[3:15:3,:,:]

    print(a.flags)
    assert not a.flags.contiguous
    a = wn.array(a)
    z = np.copy(a)

    shapes = a.weldobj.update(np.array(list(a.shape)))

    strides = []
    for s in a.strides:
        strides.append(s/z.itemsize)

    strides = a.weldobj.update(np.array(strides))

    end = 1
    for s in a.shape:
        end = end*s
    end = end

    iter_code = 'result(for(nditer({arr}, 0L, {end}L, 1L, {shapes}, {strides}), appender, \
    |b, i, e| merge(b,exp(e))))'.format(shapes=shapes, strides=strides, end=str(end), arr=a.name)

    a.weldobj.weld_code = iter_code
    z = np.exp(z)
    # convert the data represented by weldarray 'a', to a multi-dimensional numpy array of shape as z,
    # and then compare the values.
    np_a = a._eval()
    assert_correct(np_a, z)

def test_start():
    '''
    Has a different start from the base array.
    '''
    orig = np.random.rand(20,20)
    a = orig[0:20:1,0:20:1]

    start = (wn.addr(a) - wn.addr(orig)) / a.itemsize
    orig = wn.array(orig)
    z = np.copy(a)
    z = np.exp(z)

    shapes = orig.weldobj.update(np.array(list(a.shape)))

    strides = []
    for s in a.strides:
        strides.append(s/8)

    strides = orig.weldobj.update(np.array(strides))

    end = 1
    for s in a.shape:
        end = end*s
    end = end + start

    iter_code = 'result(for(nditer({arr}, {start}L, {end}L, 1L, {shapes}, {strides}), appender, \
    |b, i, e| merge(b,exp(e))))'.format(shapes=shapes, strides=strides, end=str(end),
            start=str(start), arr=orig.name)
    orig.weldobj.weld_code = iter_code
    np_a = orig._eval()

    assert_correct(np_a, z)

def test_zip():
    '''
    Has a different start from the base array.
    '''
    orig = np.random.rand(20,20)
    orig2 = np.random.rand(20,20)
    a = orig[5:20:1,3:20:2]
    b = orig2[5:20:1,3:20:2]
    start = (wn.addr(a) - wn.addr(orig)) / a.itemsize
    orig = wn.array(orig)

    # copying so we can test them later.
    z = np.copy(a)
    z2 = np.copy(b)
    # added orig2 to orig's weldobject.
    orig_2_name = orig.weldobj.update(orig2)

    shapes = orig.weldobj.update(np.array(list(a.shape)))

    strides = []
    for s in a.strides:
        strides.append(s/8)

    strides = orig.weldobj.update(np.array(strides))

    end = 1
    for s in a.shape:
        end = end*s
    end = end + start

    iter_code = 'result(for(zip(nditer({arr}, {start}l, {end}l, 1l, {shapes}, {strides}),  \
nditer({arr2}, {start}l, {end}l, 1l, {shapes}, {strides})), \
appender, |b, i, e| merge(b,e.$0+e.$1)))'.format(shapes=shapes, strides=strides, end=str(end),
start=str(start), arr=orig.name, arr2=orig_2_name)
    orig.weldobj.weld_code = iter_code
    # gives us a numpy array after evaluating the nditer code above.
    np_a = orig._eval()
    # update the copied array.
    z3 = z+z2;
    # test values are equal.
    assert_correct(np_a, z3)

# few different tests.
test_view()
test_start()
test_zip()

max_stars_repo_name: Bathtor/weld

max_issues_repo_name: Bathtor/weld

max_forks_repo_name: Bathtor/weld