text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Nov  8 11:20:37 2021

@author: lillycummins

Motive: count number of genes per COG cat from eggnog annotation files and show
#COG cat % as heatmap for each ST (optimised for visualisation in spyder v5.0.5)
"""

import os
import csv
import argparse
from collections import defaultdict
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib.patches import Patch

parser = argparse.ArgumentParser(description="Create functional distribution heatmap")
parser.add_argument("--path", "-p", help="Path to directory containing eggnog annotations", type=str, required=True) #path to directory of eggnog annotation tsv files
parser.add_argument("--output", "-o", help="Specify output file name", type=str, required=True) #path to output matrix
args = parser.parse_args()

possiblecategories=['A','B','C','D','E','F','G','H','I','J','K','L','M','N','O',
                    'P','Q','S','T','U','V','W','Y','Z']
df=pd.DataFrame()

for filename in os.listdir(args.path):
    COG_categories=list()
    ST=filename.split('_')[0]
    with open(args.path+filename) as csvfile:
        reader = csv.reader(csvfile, delimiter='\t')
        next(reader)
        for column in reader: #column 6 for narrowest or 9 for best COG category
            if '#' in column[0]:
                continue
            else:
                category = column[9]
                COG_categories.extend([category])

###COUNT INSTANCES OF EACH COG CATEGORY
    d=defaultdict(int)
    for x in possiblecategories:
        d[x]=0

    for x in COG_categories:
        if x=='-':
            d['?']+=1
        for n in list(x):
            if n in possiblecategories: #check character is category and not a comma or space
                d[n]+=1

    ditems=sorted(d.items(), key=lambda x:x[0]) #order categories alphabetically

    x, function = list(zip(*ditems)) #x is d.values, function is d.keys

    percent=np.empty(len(ditems))
    for i in range(0,len(ditems)):
        percent[i]=function[i]/len(COG_categories)*100
    percents=percent.tolist()
    row_df = pd.DataFrame([percents], index = [ST], columns=x)
    df=pd.concat([df,row_df])

##REMOVE 'EMPTY' COG CATEGORIES
def find_empty_cats(data_set):
    redundant_cols = list()
    for i in range(0,data_set.shape[1]):
        column = data_set.iloc[:,i] #examine every column
        if all(column.iloc[x] == 0 for x in range(0,data_set.shape[0])) == True: #find empty cat
            redundant_cols.append(column.name)
        else:
            continue
    return redundant_cols

df=df.drop(columns = find_empty_cats(df))

df=df.drop(columns = ['S','?']) #remove groups of unknown function

##PHYLOGROUP ASSIGNMENT
phylogroup = {'ST10':'A','ST11':'E','ST117':'F','ST12':'B2','ST127':'B2','ST131':'B2',
              'ST14':'B2','ST141':'B2','ST144':'B2','ST167':'A','ST17':'B1','ST21':'B1',
              'ST28':'B2','ST3':'B1','ST372':'B2','ST38':'D','ST410':'A','ST648':'F',
              'ST69':'D','ST73':'B2','ST95':'B2'}

phylo_ser=pd.Series(data=phylogroup,index=['ST10','ST11','ST117','ST12','ST127','ST131',
                                           'ST14','ST141','ST144','ST167','ST17','ST21',
                                           'ST28','ST3','ST372','ST38','ST410','ST648',
                                           'ST69','ST73','ST95'])

cols=['mediumpurple','salmon','pink','orange','lightgreen','c']
lut = dict(zip(phylo_ser.unique(), cols))
row_colors = phylo_ser.map(lut)

##PLOT CLUSTERMAP
plt.figure()
sns.set(font_scale = 1.2)
cmap = sns.color_palette("Blues", as_cmap=True)
cm = sns.clustermap(df, xticklabels=1, tree_kws=dict(linewidths=1.5), cmap=cmap,
                    row_colors=row_colors, cbar_pos=(.03, .6, .03, .3))
ax=cm.ax_heatmap
ax.set_xticklabels(ax.get_xticklabels(), rotation=0, horizontalalignment='right')
ax.set(ylabel=None)
handles = [Patch(facecolor=lut[name]) for name in sorted(lut)]
plt.legend(handles, sorted(lut), title='Phylogroup', prop={'size': 11},
            bbox_to_anchor=(1, 1), bbox_transform=plt.gcf().transFigure, loc = 'upper right')

plt.show()

df.to_csv(args.output+'.csv')

max_stars_repo_name: lillycummins/PaperScripts

max_issues_repo_name: lillycummins/PaperScripts

max_forks_repo_name: lillycummins/PaperScripts


text:
import torch

import numpy as np
import random

import model.meta.learner as Learner
import model.meta.modelfactory as mf
import sys

if not sys.warnoptions:
    import warnings
    warnings.simplefilter("once")

"""
Multi task
    big batch size, set increment 100 so that it is treated as 1 task with all classes in the dataset
    inference time for acc eval, use offsets
"""

class Net(torch.nn.Module):
    def __init__(self,
                 n_inputs,
                 n_outputs,
                 n_tasks,
                 args):
        super(Net, self).__init__()
        self.args = args
        self.nt = n_tasks

        self.n_feat = n_outputs
        self.n_classes = n_outputs

        arch = args.arch
        nl, nh = args.n_layers, args.n_hiddens
        config = mf.ModelFactory.get_model(model_type = arch, sizes = [n_inputs] + [nh] * nl + [n_outputs],
                                                dataset = args.dataset, args=args)
        self.net = Learner.Learner(config, args)

        # setup optimizer
        self.opt = torch.optim.SGD(self.parameters(), lr=args.lr)

        # setup losses
        self.loss = torch.nn.CrossEntropyLoss()

        self.gpu = args.cuda
        self.nc_per_task = int(n_outputs / n_tasks)
        self.n_outputs = n_outputs

    def compute_offsets(self, task):
        offset1 = task * self.nc_per_task
        offset2 = (task + 1) * self.nc_per_task
        return int(offset1), int(offset2)

    def take_multitask_loss(self, bt, logits, y):
        loss = 0.0
        for i, ti in enumerate(bt):
            offset1, offset2 = self.compute_offsets(ti)
            loss += self.loss(logits[i, offset1:offset2].unsqueeze(0), y[i].unsqueeze(0)-offset1)
        return loss/len(bt)

    def forward(self, x, t):

        output = self.net.forward(x)

        # make sure we predict classes within the current task
        if torch.unique(t).shape[0] == 1:
            offset1, offset2 = self.compute_offsets(t[0].item())
            if offset1 > 0:
                output[:, :offset1].data.fill_(-10e10)
            if offset2 < self.n_outputs:
                output[:, int(offset2):self.n_outputs].data.fill_(-10e10)
        else:
            for i in range(len(t)):
                offset1, offset2 = self.compute_offsets(t[i])
                if offset1 > 0:
                    output[i, :offset1].data.fill_(-10e10)
                if offset2 < self.n_outputs:
                    output[i, int(offset2):self.n_outputs].data.fill_(-10e10)
        return output

    def observe(self, x, y, t):
        self.net.train()

        self.net.zero_grad()
        logits = self.net.forward(x)
        loss = self.take_multitask_loss(t, logits, y)
        loss.backward()
        self.opt.step()

        return loss.item()

ext: py

max_stars_repo_name: Zotkin/La-MAML

max_issues_repo_name: Zotkin/La-MAML

max_forks_repo_name: Zotkin/La-MAML

reason: import numpy


text:
import numpy as np
from unittest import mock
import joblib

EXPECTED_RESPONSE = np.array([3, 3])

def test_model(*args, **kwargs):
    data = ["text 1", "text 2"]

    m = joblib.load("model.joblib")
    result = m.predict(data)
    assert all(result == EXPECTED_RESPONSE)

max_stars_repo_name: karthiklinux/crescendo-jenkins-x-seldon

max_issues_repo_name: karthiklinux/crescendo-jenkins-x-seldon

max_forks_repo_name: karthiklinux/crescendo-jenkins-x-seldon


text:
import unittest
import filterdesigner.IIRDesign as IIRDesign
import scipy.signal as signal
import numpy as np

class TestEllip(unittest.TestCase):

    def setUp(self):
        self.n = 3
        self.Rp = 1
        self.Rs = 120
        self.Wp1 = 0.5
        self.Wp2 = [0.25, 0.75]
        self.Wps = np.pi/2

    def test_ellip_1(self):
        # Test case for lowpass filter with default
        IIR = IIRDesign.ellip(self.n, self.Rp, self.Rs, self.Wp1)
        iir = signal.ellip(self.n, self.Rp, self.Rs, self.Wp1, btype='lowpass', fs=2)
        self.assertTrue((IIR[0] == iir[0]).all() and (IIR[1] == iir[1]).all())

    def test_ellip_2(self):
        # Test case for lowpass filter without default
        IIR = IIRDesign.ellip(self.n, self.Rp, self.Rs, self.Wp1, ftype='low')
        iir = signal.ellip(self.n, self.Rp, self.Rs, self.Wp1, btype='lowpass', fs=2)
        self.assertTrue((IIR[0] == iir[0]).all() and (IIR[1] == iir[1]).all())

    def test_ellip_3(self):
        # Test case for highpass filter
        IIR = IIRDesign.ellip(self.n, self.Rp, self.Rs, self.Wp1, ftype='high')
        iir = signal.ellip(self.n, self.Rp, self.Rs, self.Wp1, btype='highpass', fs=2)
        self.assertTrue((IIR[0] == iir[0]).all() and (IIR[1] == iir[1]).all())

    def test_ellip_4(self):
        # Test case for bandpass filter with default
        IIR = IIRDesign.ellip(self.n, self.Rp, self.Rs, self.Wp2)
        iir = signal.ellip(self.n, self.Rp, self.Rs, self.Wp2, btype='bandpass', fs=2)
        self.assertTrue((IIR[0] == iir[0]).all() and (IIR[1] == iir[1]).all())

    def test_ellip_5(self):
        # Test case for bandpass filter without default
        IIR = IIRDesign.ellip(self.n, self.Rp, self.Rs, self.Wp2, ftype='bandpass')
        iir = signal.ellip(self.n, self.Rp, self.Rs, self.Wp2, btype='bandpass', fs=2)
        self.assertTrue((IIR[0] == iir[0]).all() and (IIR[1] == iir[1]).all())

    def test_ellip_6(self):
        # Test case for bandstop filter
        IIR = IIRDesign.ellip(self.n, self.Rp, self.Rs, self.Wp2, ftype='stop')
        iir = signal.ellip(self.n, self.Rp, self.Rs, self.Wp2, btype='bandstop', fs=2)
        self.assertTrue((IIR[0] == iir[0]).all() and (IIR[1] == iir[1]).all())

    def test_ellip_7(self):
        # Test case for analog filter
        IIR = IIRDesign.ellip(self.n, self.Rp, self.Rs, self.Wps, zs='s')
        iir = signal.ellip(self.n, self.Rp, self.Rs, self.Wps, btype='lowpass', analog=True)
        self.assertTrue((IIR[0] == iir[0]).all() and (IIR[1] == iir[1]).all())

    def test_ellip_8(self):
        # test case for Exception 1
        with self.assertRaises(ValueError):
            IIRDesign.ellip(self.n, self.Rp, self.Rs, self.Wps, zs='x')

    def test_ellip_9(self):
        # test case for Exception 2
        with self.assertRaises(ValueError):
            IIRDesign.ellip(2.5, self.Rp, self.Rs, self.Wp1)

    def test_ellip_10(self):
        # test case for Exception 3
        with self.assertRaises(ValueError):
            IIRDesign.ellip(self.n, self.Rp, self.Rs, self.Wp1, ftype='x')

    def test_ellip_11(self):
        # test case for Exception 4
        with self.assertRaises(ValueError):
            IIRDesign.ellip(self.n, self.Rp, self.Rs, 1.4, zs='z')

    def test_ellip_12(self):
        # test case for Exception 5
        with self.assertRaises(ValueError):
            IIRDesign.ellip(self.n, self.Rp, self.Rs, [0.5, 1.5], zs='z')

    def test_ellip_13(self):
        # test case for Exception 6
        with self.assertRaises(ValueError):
            IIRDesign.ellip(self.n, self.Rp, self.Rs, self.Wp2, ftype='low')

    def test_ellip_14(self):
        # test case for Exception 7
        with self.assertRaises(ValueError):
            IIRDesign.ellip(self.n, self.Rp, self.Rs, self.Wp2, ftype='high')

    def test_ellip_15(self):
        # test case for Exception 8
        with self.assertRaises(ValueError):
            IIRDesign.ellip(self.n, self.Rp, self.Rs, self.Wp1, ftype='stop')

    def test_ellip_16(self):
        # test case for Exception 9
        with self.assertRaises(ValueError):
            IIRDesign.ellip(self.n, self.Rp, self.Rs, self.Wp1, ftype='bandpass')

max_stars_repo_name: Yuki-F-HCU/filterdesigner

max_issues_repo_name: Yuki-F-HCU/filterdesigner

max_forks_repo_name: Yuki-F-HCU/filterdesigner

reason: import numpy,import scipy


text:
import tensorflow as tf
from tensorflow.python.framework import ops
from tensorflow.python.framework import dtypes
import os
import numpy as np
import scipy
from scipy import misc

def get_filename_list(path, config):
    fd = open(path)
    image_filenames = []
    label_filenames = []
    for i in fd:
        i = i.strip().split(" ")
        image_filenames.append(i[0])
        label_filenames.append(i[1])

    image_filenames = [config["IMG_PREFIX"] + name for name in image_filenames]
    label_filenames = [config["LABEL_PREFIX"] + name for name in label_filenames]
    return image_filenames, label_filenames

def dataset_reader(filename_queue, config):  # prev name: CamVid_reader

    image_filename = filename_queue[0]  # tensor of type string
    label_filename = filename_queue[1]  # tensor of type string

    # get png encoded image
    imageValue = tf.read_file(image_filename)
    labelValue = tf.read_file(label_filename)

    # decodes a png image into a uint8 or uint16 tensor
    # returns a tensor of type dtype with shape [height, width, depth]
    image_bytes = tf.image.decode_png(imageValue)
    label_bytes = tf.image.decode_png(labelValue)  # Labels are png, not jpeg

    image = tf.reshape(image_bytes, (config["INPUT_HEIGHT"], config["INPUT_WIDTH"], config["INPUT_CHANNELS"]))
    label = tf.reshape(label_bytes, (config["INPUT_HEIGHT"], config["INPUT_WIDTH"], 1))

    return image, label

def dataset_inputs(image_filenames, label_filenames, batch_size, config):
    images = ops.convert_to_tensor(image_filenames, dtype=dtypes.string)
    labels = ops.convert_to_tensor(label_filenames, dtype=dtypes.string)

    filename_queue = tf.train.slice_input_producer([images, labels], shuffle=True)

    image, label = dataset_reader(filename_queue, config)
    reshaped_image = tf.cast(image, tf.float32)
    min_queue_examples = 300
    print('Filling queue with %d input images before starting to train. '
          'This may take some time.' % min_queue_examples)

    # Generate a batch of images and labels by building up a queue of examples.
    return _generate_image_and_label_batch(reshaped_image, label,
                                           min_queue_examples, batch_size,
                                           shuffle=True)

def _generate_image_and_label_batch(image, label, min_queue_examples,
                                    batch_size, shuffle):
    """Construct a queued batch of images and labels.
    Args:
        image: 3-D Tensor of [height, width, 3] of type.float32.
        label: 3-D Tensor of [height, width, 1] type.int32
        min_queue_examples: int32, minimum number of samples to retain
        in the queue that provides of batches of examples.
        batch_size: Number of images per batch.
        shuffle: boolean indicating whether to use a shuffling queue.
    Returns:
        images: Images. 4D tensor of [batch_size, height, width, 3] size.
        labels: Labels. 3D tensor of [batch_size, height, width ,1] size.
    """
    # Create a queue that shuffles the examples, and then
    # read 'batch_size' images + labels from the example queue.

    # TODO: test if setting threads to higher number!
    num_preprocess_threads = 1
    if shuffle:
        images, label_batch = tf.train.shuffle_batch(
            [image, label],
            batch_size=batch_size,
            num_threads=num_preprocess_threads,
            capacity=min_queue_examples + 3 * batch_size,
            min_after_dequeue=min_queue_examples)
    else:
        images, label_batch = tf.train.batch(
            [image, label],
            batch_size=batch_size,
            num_threads=num_preprocess_threads,
            capacity=min_queue_examples + 3 * batch_size)

    # Display the training images in the visualizer.
    tf.summary.image('training_images', images)
    print('generating image and label batch:')
    return images, label_batch

def get_all_test_data(im_list, la_list):
    images = []
    labels = []
    index = 0
    for im_filename, la_filename in zip(im_list, la_list):
        im = scipy.misc.imread(im_filename)
        la = scipy.misc.imread(la_filename)
        images.append(im)
        labels.append(la)
        index = index + 1

    print('%d images are loaded' % index)
    return images, labels

reason: import numpy,import scipy,from scipy


text:
import numpy as np
import gizeh
from ..types import Matrix4, Sequence3
from .. import v4hb
from . import Segment

def make_gizeh(seg: Segment, matrix: Matrix4, fill:Sequence3=(1, 0, 0), num_points:int=16) -> gizeh.Element:
    """Make gizeh element.

    The transformation from coordinate system of segment to 2D scene is given by matrix.
    """
    assert all(surface is not None for surface in seg.surfaces), 'Segment must be bounded by surfaces'
    assert np.prod(seg.beam.mode.shape) == 1, 'Only one PBG allowed (for now)'
    length = np.asscalar(seg.length)  # Extra dimensions confuse things.
    origin = seg.beam.mode.line_base.origin.reshape((4,))
    vector = seg.beam.mode.line_base.vector.reshape((4,))
    rhat = v4hb.cross(matrix[2, :], vector)

    def calc_point(d, sign):
        """Calculate beam edge point.

        Args:
            d (scalar): Distance along beam.
            sign (scalar): Direction along rhat.

        Returns:
            ...x4 array: Beam edge point(s).
        """
        r = seg.beam.calc_edge_distance(d, sign * rhat)
        point = origin + d * vector + sign * rhat * r
        return point

    dss = []
    for sign in (-1, 1):
        dss.append([])
        for (surface, d_guess) in zip(seg.surfaces, (0, length)):
            d_root = surface.intersect_global_curve(lambda d: calc_point(d, sign).reshape((4,)), d_guess)
            dss[-1].append(d_root)

    dm = np.linspace(dss[0][0], dss[0][1], num_points)[:, None]
    pointsm = calc_point(dm, -1).reshape((-1, 4))
    dp = np.linspace(dss[1][1], dss[1][0], num_points)[:, None]
    pointsp = calc_point(dp, 1).reshape((-1, 4))
    to_section = np.linalg.inv(matrix)[:, :2]
    points_section = np.concatenate((pointsm, pointsp)).dot(to_section)
    return gizeh.polyline(points_section, close_path=True, fill=fill)

max_stars_repo_name: draustin/otk

max_issues_repo_name: uamhforever/otk

max_forks_repo_name: uamhforever/otk


text:
import argparse
import os

from utils.utils import directory_for_dataset, directory_for_model

def parse_args():
    """
    Imitation Learning - Extension to an arbitrary number of agents

    usage: task1_extensions.py [--args]

    :return args:
                  --help                Show this help message and exit
                  --gui                 Run simulation using the gui (default: False)
                  --n-simulations N     Number of runs for each simulation (default: 1000)
                  --task TASK
                                        Choose the task to perform in the current execution
                                        between task1 and task2 (default: task1)
                  --myt-quantity MYT_QUANTITY
                                        Number of thymios for the simulation (default: variable)
                  --avg-gap AVG_GAP
                                        Average gap distance between thymios (default: variable)
                  --generate-dataset    Generate the dataset containing the simulations
                                        (default: False)
                  --generate-split      Generate the indices for the split of the dataset
                                        (default: False)
                  --plots-dataset       Generate the plots of regarding the dataset
                                        (default: False)
                  --check-dataset       Generate the plots that check the dataset conformity
                                        (default: False)
                  --compare-all         Generate plots that compare all the experiments in
                                        terms of distance from goal (default: False)
                  --generate-animations
                                        Generate animations that compare the controllers
                                        (default: False)
                  --controller CONTROLLER
                                        Choose the controller for the current execution.
                                        Usually between all, learned, manual and omniscient
                                        (default: all)
                  --dataset-folder DATASET_FOLDER
                                        Name of the directory containing the datasets
                                        (default: datasets)
                  --dataset DATASET     Choose the datasets to use in the current execution
                                        (default: all)
                  --models-folder MODELS_FOLDER
                                        Name of the directory containing the models
                                        (default: models)
                  --model-type MODEL_TYPE
                                        Name of the sub-directory containing the models
                                        (default: distributed)
                  --model MODEL         Name of the model (default: net1)
                  --train-net           Train the model (default: False)
                  --save-net            Save the model in onnx format (default: False)
                  --net-input SENSING
                                        Choose the input of the net between prox_values,
                                        prox_comm or all_sensors (default: all_sensors)
                  --plots-net           Generate the plots of regarding the model
                                        (default: False)

    """
    parser = argparse.ArgumentParser(description='Imitation Learning - Extension to an arbitrary number of agents')

    parser.add_argument('--gui', action="store_true",
                        help='Run simulation using the gui (default: False)')
    parser.add_argument('--n-simulations', type=int, default=1000, metavar='N',
                        help='Number of runs for each simulation (default: 1000)')
    parser.add_argument('--task', default='task1', choices=['task1', 'task2'],
                        help='Choose the task to perform in the current execution between task1 and task2 (default: task1)')

    parser.add_argument('--myt-quantity', default='variable', type=str,
                        help='Number of thymios for the simulation (default: variable)')
    parser.add_argument('--avg-gap', default='variable', type=str,
                        help='Average gap distance between thymios (default: variable)')

    parser.add_argument('--generate-dataset', action="store_true",
                        help='Generate the dataset containing the simulations (default: False)')
    parser.add_argument('--generate-split', action="store_true",
                        help='Generate the indices for the split of the dataset (default: False)')

    parser.add_argument('--plots-dataset', action="store_true",
                        help='Generate the plots of regarding the dataset (default: False)')
    parser.add_argument('--check-dataset', action="store_true",
                        help='Generate the plots that check the dataset conformity (default: False)')
    parser.add_argument('--compare-all', action="store_true",
                        help='Generate plots that compare all the experiments in terms of distance from goal ('
                             'default: False)')
    parser.add_argument('--generate-animations', action="store_true",
                        help='Generate animations that compare the controllers (default: False)')

    parser.add_argument('--controller', default='all', type=str,
                        help='Choose the controller for the current execution. Usually between all, learned, '
                             'manual and omniscient (default: all)')

    parser.add_argument('--dataset-folder', default='datasets', type=str,
                        help='Name of the directory containing the datasets (default: datasets)')
    parser.add_argument('--dataset', default='all', type=str,
                        help='Choose the datasets to use in the current execution (default: all)')

    parser.add_argument('--models-folder', default='models', type=str,
                        help='Name of the directory containing the models (default: models)')
    parser.add_argument('--model-type', default='distributed', type=str,
                        help='Name of the sub-directory containing the models (default: distributed)')
    parser.add_argument('--model', default='net1', type=str, help='Name of the model (default: net1)')

    parser.add_argument('--train-net', action='store_true', help='Train the model (default: False)')
    parser.add_argument('--save-net', action='store_true', help='Save the model in onnx format (default: False)')

    parser.add_argument('--net-input', default='all_sensors', choices=['prox_values', 'prox_comm', 'all_sensors'],
                        help='Choose the input of the net between prox_values and prox_comm_events (default: '
                             'all_sensors)')
    parser.add_argument('--plots-net', action="store_true",
                        help='Generate the plots of regarding the model (default: False)')

    args = parser.parse_args()

    return args

if __name__ == '__main__':

    args = parse_args()

    d = os.path.join(args.dataset_folder, args.task, args.net_input, args.dataset)

    if args.controller == 'all':
        controllers = ['omniscient', 'learned_distributed', 'learned_communication', 'manual']
    else:
        controllers = [args.controller]

    runs_dir_omniscient = os.path.join(d, 'omniscient')
    runs_dir_manual = os.path.join(d, 'manual')
    runs_dir_learned = os.path.join(d, 'learned')
    runs_dir_learned_dist = os.path.join(d, 'learned_distributed')
    runs_dir_learned_comm = os.path.join(d, 'learned_communication')

    for c in controllers:
        run_dir, run_img_dir, run_video_dir = directory_for_dataset(d, c)
        model_dir, model_img_dir, model_video_dir, metrics_path = directory_for_model(args)

        if args.model_type == 'communication':
            communication = True
        elif args.model_type == 'distributed':
            communication = False
        else:
            raise ValueError('Invalid value for model_type')

        if args.generate_dataset:
            from generate_simulation_data import GenerateSimulationData as sim

            print('Generating %s simulations for %s %s controller…' % (args.n_simulations, d, c))

            if c == 'learned':
                sim.generate_simulation(run_dir=run_dir, n_simulations=args.n_simulations, controller=c,
                                        args=args, model_dir=model_dir, model=args.model, communication=communication)
            else:
                sim.generate_simulation(run_dir=run_dir, n_simulations=args.n_simulations, controller=c,
                                        args=args, communication=communication)

        if args.plots_dataset:
            from utils.my_plots import visualise_simulation_over_time_all_sensors, plot_compared_distance_compressed, \
                                       thymio_quantity_distribution, visualise_position_over_time, \
                                       visualise_control_over_time, plot_predicted_message

            print('Generating plots for %s %s controller…' % (d, c))
            runs = [0, 1, 2, 3, 4, 5, 251, 634]
            # for i in runs:
            #     visualise_simulation_over_time_all_sensors(run_dir, run_img_dir, i,
            #                                                'Thymio positions over time - Simulation run %d - %s %s'
            #                                                % (i, args.net_input, c))
            #
            # # FIXME substitute controllers with datasets
            # datasets = ['omniscient']
            # plot_compared_distance_compressed([run_dir], run_img_dir, controllers,
            #                                   'Robot distances from goal - %s' % (args.net_input),
            #                                   'distances-from-goal-compressed')
            #
            # plot_compared_distance_compressed([run_dir], run_img_dir, controllers,
            #                                   'Robot distances from goal - %s' % (args.net_input),
            #                                   'distances-from-goal-absolute-compressed', absolute=False)
            if args.myt_quantity != 'variable':
                visualise_position_over_time(run_dir, run_img_dir, 'position-overtime-%s' % c)
            # visualise_control_over_time(run_dir, run_img_dir, 'control-overtime-%s' % c)

            thymio_quantity_distribution(run_dir, run_img_dir, 'Thymio quantity distribution - %s' % (args.net_input),
                                         'thymio-quantity-distribution')

            if communication:
                from utils.my_plots import visualise_communication_simulation, visualise_communication_vs_control, \
                                           visualise_communication_vs_distance
                for i in range(5):
                    visualise_communication_simulation(run_dir, run_img_dir, i,
                                                       'Simulation run %d - %s %s - communication' %
                                                       (i, args.net_input, c))

                visualise_communication_vs_control(run_dir, run_img_dir,
                                                   'Communication vs Control - %s %s' % (args.net_input, c))

                visualise_communication_vs_distance(run_dir, run_img_dir,
                                                    'Communication vs Distance from goal - %s %s' % (args.net_input, c))
                    plot_predicted_message(run_dir, run_img_dir, i)

        if args.check_dataset:
            from generate_simulation_data import GenerateSimulationData as sim
            print('\nChecking conformity of %s %s dataset…' % (d, c))
            sim.check_dataset_conformity(run_dir, run_img_dir, 'Dataset - %s %s' % (args.net_input, c), c,
                                         net_input=args.net_input, communication=communication)

        if args.train_net or args.plots_net:
            from utils.utils import prepare_dataset

            indices = prepare_dataset(run_dir, args.generate_split, args.n_simulations)
            file_losses = os.path.join(model_dir, 'losses.npy')

            if args.train_net:
                from network_training_distribute import network_train
                network_train(indices, file_losses, runs_dir_omniscient, model_dir, args.model, communication,
                              net_input=args.net_input, save_net=args.save_net)

            if args.plots_net:
                from network_evaluation import network_evaluation

                network_evaluation(indices, file_losses, runs_dir_omniscient, model_dir, args.model, model_img_dir,
                                   'omniscient', 'manual', communication, net_input=args.net_input, task=args.task)

    if args.compare_all:
        from utils.my_plots import plot_compared_distance_compressed, test_controller_given_init_positions

        print('\nGenerating comparison plots among all datasets of type %s…' % (args.net_input))
        runs_img_dir = os.path.join(d, 'images')

        dataset_folders_comm = [runs_dir_omniscient, runs_dir_manual, runs_dir_learned_dist, runs_dir_learned_comm]
        datasets_comm = ['omniscient', 'manual', 'distributed', 'communication']

        dataset_folders_dist = [runs_dir_omniscient, runs_dir_manual, runs_dir_learned_dist]
        datasets_dist = ['omniscient', 'manual', 'distributed']

        plot_compared_distance_compressed(dataset_folders_dist, runs_img_dir, datasets_dist,
                                         'Robot distances from goal', 'distances-from-goal-compressed-distributed')

        plot_compared_distance_compressed(dataset_folders_comm, runs_img_dir, datasets_comm,
                                         'Robot distances from goal', 'distances-from-goal-compressed-communication')

        # Evaluate the learned controllers by passing a specific initial position configuration and compare them with
        # the omniscient and the manual controllers
        test_controller_given_init_positions(runs_img_dir, args.model, args.net_input)

    if args.generate_animations:
        from utils.my_plots import animate_simulation, plot_simulations, visualise_position_over_time
        from utils.utils import generate_fake_simulations, generate_init_positions, check_dir
        import numpy as np
        animations_dir = os.path.join(d, 'animations')
        check_dir(animations_dir)

        plots_dir = os.path.join(d, 'plots')
        check_dir(plots_dir)

        # Create a simulation for each of the controller using the same initial position
        if args.myt_quantity == 'variable':
            myt_quantities = np.arange(5, 11)
        else:
            myt_quantities = [int(args.myt_quantity)]

        for N in myt_quantities:
            dir = os.path.join(plots_dir, 'N%d' % N)
            check_dir(dir)
            _ = generate_fake_simulations(dir, args.model, N, simulations=100)
            out_dirs = [os.path.join(dir, 'omniscient'), os.path.join(dir, 'manual'), os.path.join(dir, 'distributed'), os.path.join(dir, 'communication')]

            for c in out_dirs:
                controller = os.path.basename(os.path.normpath(c))
                img_dir = os.path.join(c, 'images')
                check_dir(img_dir)
                visualise_position_over_time(c, img_dir, 'position-overtime-%s' % controller)

            # FIXME
            # dir = os.path.join(animations_dir, 'N%d' % N)
            # check_dir(dir)
            # out_dirs = generate_fake_simulations(dir, args.model, N, simulations=1)
            #
            # for c in out_dirs:
            #     controller = os.path.basename(os.path.normpath(c))
            #     img_dir = os.path.join(c, 'images')
            #     check_dir(img_dir)
            #
            #     animate_simulation(out_dirs, N)
            #     plot_simulations(out_dirs, N)


text:
# Copyright (c) Microsoft Corporation. All rights reserved.
# Licensed under the MIT license.

import time
import argparse
from os.path import join as pjoin

import numpy as np

import textworld
from textworld import g_rng

from textworld.generator.game import GameOptions

def generate_never_ending_game(args):
    g_rng.set_seed(args.seed)

    msg = "--max-steps {} --nb-objects {} --nb-rooms {} --quest-length {} --quest-breadth {} --seed {}"
    print(msg.format(args.max_steps, args.nb_objects, args.nb_rooms, args.quest_length, args.quest_breadth, g_rng.seed))
    print("Generating game...")

    options = GameOptions()
    options.seeds = g_rng.seed
    options.nb_rooms = args.nb_rooms
    options.nb_objects = args.nb_objects
    options.quest_length = args.quest_length
    options.quest_breadth = args.quest_breadth

    game = textworld.generator.make_game(options)
    if args.no_quest:
        game.quests = []

    game_name = "neverending"
    path = pjoin(args.output, game_name + ".ulx")
    options = textworld.GameOptions()
    options.path = path
    options.force_recompile = True
    game_file = textworld.generator.compile_game(game, options)
    return game_file

def benchmark(game_file, args):
    infos = textworld.EnvInfos()
    if args.activate_state_tracking or args.mode == "random-cmd":
        infos.admissible_commands = True

    if args.compute_intermediate_reward:
        infos.intermediate_reward = True

    env = textworld.start(game_file, infos)
    print("Using {}".format(env))

    if args.mode == "random":
        agent = textworld.agents.NaiveAgent()
    elif args.mode == "random-cmd":
        agent = textworld.agents.RandomCommandAgent(seed=args.agent_seed)
    elif args.mode == "walkthrough":
        agent = textworld.agents.WalkthroughAgent()

    agent.reset(env)
    game_state = env.reset()

    if args.verbose:
        env.render()

    reward = 0
    done = False
    nb_resets = 1
    start_time = time.time()
    for _ in range(args.max_steps):
        command = agent.act(game_state, reward, done)
        game_state, reward, done = env.step(command)

        if done:
            env.reset()
            nb_resets += 1
            done = False

        if args.verbose:
            env.render()

    duration = time.time() - start_time
    speed = args.max_steps / duration
    msg = "Done {:,} steps in {:.2f} secs ({:,.1f} steps/sec) with {} resets."
    print(msg.format(args.max_steps, duration, speed, nb_resets))
    return speed

def parse_args():
    parser = argparse.ArgumentParser()
    parser.add_argument("--game",
                        help="Path to an existing game.")
    parser.add_argument("--nb-rooms", type=int, default=20,
                        help="Nb. of rooms in the world. Default: %(default)s")
    parser.add_argument("--nb-objects", type=int, default=50,
                        help="Nb. of objects in the world. Default: %(default)s")
    parser.add_argument("--quest-length", type=int, default=5,
                        help="Minimum nb. of actions the quest requires to be completed. Default: %(default)s")
    parser.add_argument("--quest-breadth", type=int, default=1,
                        help="Control how non-linear a quest can be. Default: %(default)s")
    parser.add_argument("--max-steps", type=int, default=1000,
                        help="Stop the game after that many steps. Default: %(default)s")
    parser.add_argument("--output", default="./gen_games/",
                        help="Output folder to save generated game files.")
    parser.add_argument("--mode", default="random-cmd", choices=["random", "random-cmd", "walkthrough"])
    parser.add_argument("--batch-size", type=int, default=1)
    parser.add_argument("--no-quest", action="store_true")
    parser.add_argument("--compute_intermediate_reward", action="store_true")
    parser.add_argument("--activate_state_tracking", action="store_true")
    parser.add_argument("--seed", type=int)
    parser.add_argument("--agent-seed", type=int, default=2018)
    parser.add_argument("-v", "--verbose", action="store_true")
    return parser.parse_args()

if __name__ == "__main__":
    args = parse_args()
    game_file = args.game
    if args.game is None:
        game_file = generate_never_ending_game(args)

    speeds = []
    for _ in range(3):
        speed = benchmark(game_file, args)
        speeds.append(speed)
        args.agent_seed = args.agent_seed + 1

    print("-----\nAverage: {:,.1f} steps/sec".format(np.mean(speeds)))


text:
import matplotlib
# Force matplotlib to not use any Xwindows backend.
matplotlib.use('Agg')
import numpy as np, pandas as pd, matplotlib.pyplot as plt; plt.rc('text', usetex=True)
import seaborn as sns; sns.set(style="white", context="talk")
from  matplotlib import colors, ticker, cm
from astroML.utils import convert_2D_cov
from astroML.plotting.tools import draw_ellipse
from astroML.plotting import setup_text_plots
setup_text_plots(fontsize=8, usetex=False)
from scipy import stats
from sklearn.cluster import KMeans
import scipy.cluster.hierarchy as sch
from sklearn.preprocessing import scale
import rpy2.robjects as robjects
from rpy2.robjects.packages import importr
circularlib=importr('circular')
circular=robjects.r('circular')
corcircular=robjects.r('cor.circular')
import pickle
from LIMBR import imputation

def plot_cormat(df,fname,ptitle):
    fig, ax = plt.subplots(nrows=1, ncols=1)
    cors = np.corrcoef(df.values.T)
    z_min, z_max = -np.abs(cors).max(), np.abs(cors).max()
    pc = ax.pcolor(cors, cmap='seismic', vmin=z_min, vmax=z_max)
    for i in range(0,cors.shape[1],3):
        ax.axvline(x=i, linestyle='-', color='black', linewidth=0.7,    zorder=3)
    for i in range(0,cors.shape[0],3):
        ax.axhline(y=i, linestyle='-', color='black', linewidth=0.7,    zorder=3)
    tpts = [i*2 for i in range(1,(int(cors.shape[1]/3)+1))]
    ax.set_xticks([i+1.5 for i in range(0,cors.shape[1],3)])
    ax.set_xticklabels([(i-12) % 22 for i in tpts], fontsize=18)
    ax.set_yticks([i+1.5 for i in range(0,cors.shape[1],3)])
    ax.set_yticklabels([(i-12) % 22 for i in tpts], fontsize=18)
    ax.set_title(ptitle, fontsize=16)
    fig.subplots_adjust(right=0.85)
    cbar_ax = fig.add_axes([0.87, 0.15, 0.02, 0.7])
    cb = fig.colorbar(pc, cax=cbar_ax)
    cb.set_label('Correlation', fontsize=12)
    plt.savefig(fname+'.pdf')
    plt.close()

def plot_cormat_prep(df,fname,ptitle):
    fig, ax = plt.subplots(nrows=1, ncols=1)
    cors = np.corrcoef(df.values.T)
    z_min, z_max = -np.abs(cors).max(), np.abs(cors).max()
    pc = ax.pcolor(cors, cmap='seismic', vmin=z_min, vmax=z_max)
    for i in prep_breaks:
        ax.axvline(x=i, linestyle='-', color='black', linewidth=0.7,    zorder=3)
    for i in prep_breaks:
        ax.axhline(y=i, linestyle='-', color='black', linewidth=0.7,    zorder=3)
    ax.set_xticks([prep_breaks[i] + prep_diffs[i]/2 for i in range(len(prep_breaks))])
    ax.set_xticklabels([i+1 for i in range(16)], fontsize=18)
    ax.set_yticks([prep_breaks[i] + prep_diffs[i]/2 for i in range(len(prep_breaks))])
    ax.set_yticklabels([i+1 for i in range(16)], fontsize=18)
    ax.set_title(ptitle, fontsize=16)
    fig.subplots_adjust(right=0.85)
    cbar_ax = fig.add_axes([0.87, 0.15, 0.02, 0.7])
    cb = fig.colorbar(pc, cax=cbar_ax)
    cb.set_label('Correlation', fontsize=18)
    plt.savefig(fname+'.pdf')
    plt.close()

def scale_df(df):
    sc = df.apply(lambda x: x-np.mean(x),axis=1)
    sc = sc.apply(lambda x: x/np.std(x),axis=1)
    return sc

def reorder_heatmap(arr):
    Y = sch.linkage(arr, method='complete', metric='euclidean')
    Z = sch.dendrogram(Y, orientation='left',no_plot=True)
    idx = Z['leaves']
    return arr[idx,:]

def plot_heatmap(df,fname,ptitle):
    z_min, z_max = -np.abs(df).max(), np.abs(df).max()
    random_state = 42
    y_pred = KMeans(n_clusters=2, random_state=random_state).fit_predict(df)
    fig, ax = plt.subplots(nrows=1, ncols=1)
    ordered = reorder_heatmap(df)
    pc = ax.pcolor(ordered, cmap='seismic', vmin=z_min*.8, vmax=z_max*.8)
    for i in range(0,ordered.shape[1],3):
        ax.axvline(x=i, linestyle='-', color='black', linewidth=0.4,    zorder=3)
    ax.set_xticks([i+1.5 for i in range(0,ordered.shape[1],3)])
    ax.set_xticklabels([i*2 for i in range(1, (int(ordered.shape[1]/3)+1))], fontsize=18)
    ax.set_yticks([])
    ax.set_title(ptitle, fontsize=16)
    fig.subplots_adjust(right=0.85)
    cbar_ax = fig.add_axes([0.87, 0.15, 0.02, 0.7])
    cb = fig.colorbar(pc, cax=cbar_ax)
    cbar_ax.text(2.5,.49,r'$\sigma$', fontsize=12)
    plt.savefig(fname+'.pdf')
    plt.close()

def autocorr(l,shift):
    return np.dot(l, np.roll(l, shift)) / np.dot(l, l)

def autocorr(x,shift):
    result = np.correlate(x, x, mode='full')
    return result[result.size//2:][shift]

def get_tpoints(df):
    tpoints = [i.replace('CT','') for i in df.columns.values]
    tpoints = [int(i.split('.')[0]) if '.' in i else int(i.split('_')[0]) for i in tpoints]
    return np.asarray(tpoints)

def circ_cor(df):
    per = 6
    cors = []
    tpoints = get_tpoints(df)
    for row in df.values:
        ave = []
        #might eventually need to account for case where all replicates of a timepoint are missing (in this case the experiment is probably irreparably broken anyway though)
        for k in set(tpoints):
            ave.append((np.mean([row[i] for i, j in enumerate(tpoints) if j == k])))
        cors.append(autocorr(ave,(per*2))-autocorr(ave,per))
    return np.asarray(cors)

def get_fstats(ts,cs,ctypes):
    results = {}
    results2 = {}
    for class_val in ctypes:
        groups = []
        for i in list(set(cs[class_val].values)):
            groups.append(list(cs[cs[class_val]==i]['Sample ID'].values))
        fortest = {}
        for j in np.sort(list(set(ts['Trend'].values))):
            templist = []
            for i in range(len(groups)):
                templist.append(list(ts[ts['variable'].isin(groups[i]) & (ts['Trend'] == j)]['value'].values))
            fortest[j] = templist
        tempdict = {}
        tempicc = {}
        for k, v in fortest.items():
            tempdict[k] = stats.f_oneway(*v)[1]
            tempicc[k] = icc(v)
        results[class_val] = tempdict
        results2[class_val] = tempicc
    return pd.DataFrame.from_dict(results), pd.DataFrame.from_dict(results2)

def gen_biasmap(df1,df2,fname,ptitle):
    df1.sort_index(ascending=False,inplace=True)
    df2.sort_index(ascending=False,inplace=True)
    fig, ax = plt.subplots(nrows=1, ncols=1)
    pc = ax.pcolor(df1.values, cmap='viridis')
    ax.set_xticks([i+.5 for i in range(len(df1.columns.values))])
    ax.set_xticklabels(df1.columns.values, fontsize=18)
    ax.set_yticks([i+.5 for i in range(len(df1))])
    ylabs = [i for i in range(1,len(df1)+1)]
    ylabs.reverse()
    ax.set_yticklabels(ylabs, fontsize=18)
    ax.xaxis.tick_top()
    ax.set_ylabel('Bias Trend')
    ax.set_title(ptitle, fontsize=16, y=1.08)
    for y in range(df2.values.shape[0]):
        for x in range(df2.values.shape[1]):
            if df2.values[y,x] < .001/df2.values.shape[1]:
                plt.text(x + 0.5, y + 0.5, '*', horizontalalignment='center',verticalalignment='center', fontsize=16)
    fig.subplots_adjust(right=0.85)
    cbar_ax = fig.add_axes([0.87, 0.15, 0.02, 0.7])
    cb = fig.colorbar(pc, cax=cbar_ax)
    cb.set_label('ICC', fontsize=12)
    plt.savefig(fname+'.pdf')
    plt.close()

def icc(arr):
    m = np.mean([np.mean(i) for i in arr])
    alphas = [np.mean(i) for i in arr]-m
    es = [list(arr[i]-m-alphas[i]) for i in range(len(arr))]
    es = [item for sublist in es for item in sublist]
    return (np.std(alphas)**2)/(np.std(alphas)**2+np.std(es)**2)

def scale_peps(df):
    sc = df.apply(lambda x: x-np.mean(x),axis=0)
    sc = sc.apply(lambda x: x/np.std(x),axis=0)
    return sc

def plot_cormat_pep(arr,fname,ptitle):
    fig, ax = plt.subplots(nrows=1, ncols=1)
    #z_min, z_max = -np.abs(arr).max(), np.abs(arr).max()
    z_min, z_max = arr.min(), arr.max()
    pc = ax.pcolor(arr, cmap='seismic', vmin=z_min, vmax=z_max)
    ax.set_xticks([i+.5 for i in range(0,arr.shape[1],1)])
    ax.set_yticks([i+.5 for i in range(0,arr.shape[1],1)])
    ax.set_xticklabels(range(1, arr.shape[1]+1), fontsize=18)
    ax.set_yticklabels(range(1, arr.shape[1]+1), fontsize=18)
    ax.set_title(ptitle, fontsize=16)
    fig.subplots_adjust(right=0.85)
    cbar_ax = fig.add_axes([0.87, 0.15, 0.02, 0.7])
    cb = fig.colorbar(pc, cax=cbar_ax)
    cb.set_label('Correlation', fontsize=12)
    plt.savefig(fname+'.pdf')
    plt.close()

def phase_hist(phases,ptitle):
    N = 11
    phases = [i%22 for i in phases]
    radii = np.histogram(phases,N)[0]
    bottom = 0
    theta = np.linspace(0.0, 2 * np.pi, N, endpoint=False)
    width = (2*np.pi) / N
    ax = plt.subplot(111, polar=True)
    ax.set_theta_offset(0.5*np.pi)
    ax.set_theta_direction(-1)
    ax.set_rlabel_position(180)
    ax.set_xticks(np.pi/180. * np.linspace(0,  360, N, endpoint=False))
    ax.set_xticklabels([str(i)+' hrs' for i in range(0, 2*N+1, 2)], fontsize=18)
    bars = ax.bar(theta, radii, width=width, bottom=bottom)
    # Use custom colors and opacity
    for r, bar in zip(radii, bars):
        bar.set_facecolor('black')
        bar.set_alpha(0.8)
    plt.title(ptitle, y=1.08, fontsize=16)
    ax.text(0.85, 0.05, r'$\rho_{cc}$ = %f'"\n"r'p = %e'%(float(circorr[0].r_repr()), float(circorr[2].r_repr())), transform=ax.transAxes, fontsize=8, verticalalignment='top', bbox=dict(facecolor='white'))
    plt.tight_layout()
    plt.savefig('output/figs/WT_CSP_phasehist.pdf')
    plt.close()

wt_p = pd.read_csv('output/actual/wt_lowess_normed__jtkout_GammaP.txt',sep='\t',index_col=0)
csp_p = pd.read_csv('output/actual/csp_lowess_normed__jtkout_GammaP.txt',sep='\t',index_col=0)
rna_p = pd.read_csv('output/actual/rna_lowess_normed__jtkout_GammaP.txt',sep='\t',index_col=0)
wt_old_p = pd.read_csv('output/actual/wt_old_normed__jtkout_GammaP.txt',sep='\t',index_col=0)
csp_old_p = pd.read_csv('output/actual/csp_old_normed__jtkout_GammaP.txt',sep='\t',index_col=0)
rna_old_p = pd.read_csv('output/actual/rna_old_normed__jtkout_GammaP.txt',sep='\t',index_col=0)

wt = pd.read_csv('output/actual/wt_lowess_normed.txt',sep='\t',index_col=0)
csp = pd.read_csv('output/actual/csp_lowess_normed.txt',sep='\t',index_col=0)
rna = pd.read_csv('output/actual/rna_lowess_normed.txt',sep='\t',index_col=0)
wt_old = pd.read_csv('output/actual/wt_old_normed.txt',sep='\t',index_col=0)
csp_old = pd.read_csv('output/actual/csp_old_normed.txt',sep='\t',index_col=0)
rna_old = pd.read_csv('output/actual/rna_old_normed.txt',sep='\t',index_col=0)

#Figure 1A
wt_class = pd.read_csv('data/wt_classes.csv')
csp_class = pd.read_csv('data/csp_classes.csv')

plot_cormat(scale_df(wt_old),'output/figs/initial_cormat','Correlation Matrix of Traditionally Normalized WT Samples Ordered by Timepoint')
new_cols = [i.replace('_1','') for i in wt_class[~wt_class['Sample ID'].isin(['CT00_3', 'CT00_2', 'CT00_1'])]['Sample ID'].values]
new_cols = [i.replace('_2','.1') for i in new_cols]
new_cols = [i.replace('_3','.2') for i in new_cols]
#new_cols = [i for i in new_cols if i not in ['CT00.2', 'CT00.1', 'CT00']]
prep = wt_class[~wt_class['Sample ID'].isin(['CT00_3', 'CT00_2', 'CT00_1'])]['Prep Set'].values
prep_breaks = []
curr_class = prep[0]
for i in range(len(prep)):
    if prep[i] != curr_class:
        prep_breaks.append(i)
    curr_class = prep[i]

prep_diffs = [4] + [prep_breaks[i+1]-prep_breaks[i] for i in range(len(prep_breaks)-1)] +[2]
prep_breaks = [0] + prep_breaks

plot_cormat_prep(scale_df(wt_old[new_cols]),'output/figs/initial_cormat_prep','Correlation Matrix of Traditionally Normalized WT Samples Ordered by Prep Set')

#Figure 1B
rawdata = pd.read_csv('data/Formatted_input.txt',sep='\t')
dist = rawdata.isnull().sum(axis=1).values.tolist()
density = stats.gaussian_kde(dist)
xs = np.linspace(0,max(dist),200)
density.covariance_factor = lambda : .1
density._compute_covariance()
values, base = np.histogram(dist, bins=100)
cumulative = np.cumsum(values)

fig, ax1 = plt.subplots()
ax1.plot(xs,density(xs),'b')
ax1.set_xlabel('Number of Missing Values')
# Make the y-axis label, ticks and tick labels match the line color.
ax1.set_ylabel('Density of Peptides', color='b')
ax1.tick_params('y', colors='b')

ax2 = ax1.twinx()
ax2.plot(base[:-1], cumulative, c='r')
ax2.set_ylabel('Cumulative Number of Peptides \n at or Below Missingness Level', color='r')
ax2.tick_params('y', colors='r')

fig.tight_layout()
plt.savefig('output/figs/missingness.pdf')
plt.close()

#Figure 2A
simres = pd.read_csv('output/simdata/simdata_mb.csv',index_col=0)
plt.scatter(simres['Base_FPR'].values,simres['Base_TPR'].values,marker='x', label='Baseline classification without noise',color='g')
plt.scatter(simres['Circ_FPR'].values,simres['Circ_TPR'].values,marker='x', label='Circadian denoised classification',color='b')
plt.scatter(simres['Noise_FPR'].values,simres['Noise_TPR'].values,marker='x', label='Classification with noise',color='r')
plt.axvline(x=0.05,color='black',ls='dashed')
plt.legend(loc="lower right")
plt.xlabel('1 - Specificity')
plt.ylabel('Sensitivity')
plt.xlim(0,1)
plt.ylim(0,1)
plt.savefig('output/figs/sim_ROC_mb.pdf')
plt.close()

#Figure 2B
simres0 = pd.read_csv('output/simdata/simdata_mb.csv')
melted_simres = pd.melt(simres0, id_vars=['Iteration'], value_vars=[r'Base_auc', r'Noise_auc',r'Circ_auc'])
ax = sns.boxplot(x='variable',y='value',data=melted_simres)
ax = sns.swarmplot(x='variable',y='value',data=melted_simres)
plt.savefig('output/figs/sim_ROC_auc_mb.pdf')
plt.close()

#Figure 2C
key = pd.read_csv('output/simdata/mb_simulated_data_key_1.txt',sep='\t',index_col=0)
simcirc_base = pd.read_csv('output/simdata/mb_simulated_data_baseline_1.txt',sep='\t',index_col=0)[key['circ']==1]
simcirc_noise = pd.read_csv('output/simdata/mb_simulated_data_with_noise_1.txt',sep='\t',index_col=0)[key['circ']==1]
simcirc_denoise = pd.read_csv('output/simdata/mb_denoised_circ_lowess_1.txt',sep='\t',index_col=0)[key['circ']==1]

plot_heatmap(simcirc_base.values,r'output/figs/simulated_circ_genes_mb','Expression as a Function of Timepoint in Simulated Circadian Genes')
plot_heatmap(simcirc_noise.values,r'output/figs/simulated_noise_circ_genes_mb','Expression as a Function of Timepoint in Simulated Circadian Genes with Noise')
plot_heatmap(simcirc_denoise.values,r'output/figs/simulated_denoise_circ_genes_mb','Expression as a Function of Timepoint in Denoised Simulated Circadian Genes')

#Figure S1
acorr = key.copy()
acorr['cors'] = circ_cor(pd.read_csv('output/simdata/mb_simulated_data_baseline_1.txt',sep='\t',index_col=0))
acorr['Circadian'] = acorr['circ'].apply(lambda x: 'Circ' if x == 1 else 'Non-Circ')
acorr['included'] = [(i<(np.percentile(acorr['cors'].values,25))) for i in acorr['cors'].values]
plt.axhline(y=np.percentile(acorr['cors'].values,25),color='black',ls='dashed')
sns.swarmplot(x='circ', y='cors', hue='included', data = acorr)
plt.ylabel(r'$\Delta$ Autocorrelation')
plt.savefig('output/figs/acorr_val_mb.pdf')
plt.close()

#Fig 3 B and C
def get_point(fname_jtk, fname_key):
    jtk = pd.read_csv(fname_jtk,sep='\t', index_col=0)
    key = pd.read_csv(fname_key,sep='\t', index_col=0)
    merged = pd.merge(jtk, key,left_index=True,right_index=True)
    return ((1-len(merged[(merged['GammaP']>.05)&(merged['circ']==0)])/len(merged[merged['circ']==0]))),(len(merged[(merged['GammaP']<.05)&(merged['circ']==1)])/len(merged[merged['circ']==1]))

def get_auc(fname_jtk, fname_key):
    jtk = pd.read_csv(fname_jtk,sep='\t', index_col=0)
    key = pd.read_csv(fname_key,sep='\t', index_col=0)
    merged = pd.merge(jtk, key,left_index=True,right_index=True)
    scoresn = np.asarray([1-i for i in merged['GammaP'].values.astype(float)])
    yn = merged['circ'].values.astype(int)
    fprn, tprn, thresholdsn = metrics.roc_curve(yn, scoresn)
    return metrics.auc(fprn, tprn)

baseline = []
noise = []
eigenms = []
limbr = []

for i in range(1,20):
    baseline.append(get_point('output/simdata/mb_simulated_data_baseline_'+str(i)+'__jtkout_GammaP.txt','output/simdata/mb_simulated_data_key_'+str(i)+'.txt'))
    noise.append(get_point('output/simdata/simulated_data_with_noise_'+str(i)+'__jtkout_GammaP.txt','output/simdata/mb_simulated_data_key_'+str(i)+'.txt'))
    eigenms.append(get_point('output/simdata/mb_simdata_eigenMS_'+str(i)+'__jtkout_GammaP.txt','output/simdata/mb_simulated_data_key_'+str(i)+'.txt'))
    limbr.append(get_point('output/simdata/mb_denoised_circ_lowess_'+str(i)+'__jtkout_GammaP.txt','output/simdata/mb_simulated_data_key_'+str(i)+'.txt'))

plt.scatter(*zip(*baseline),color='b',marker='x', label='Baseline ROC curve')

plt.scatter(*zip(*noise),color='r',marker='x', label='Noise ROC curve')

plt.scatter(*zip(*eigenms),color='y',marker='x', label='EigenMS ROC curve')

plt.scatter(*zip(*limbr),color='g',marker='x', label='LIMBR ROC curve')

plt.axvline(x=0.05,color='black',ls='dashed')
plt.legend(loc="lower right")
plt.axis([0, 1, 0, 1])
plt.xlabel('1 - Specificity')
plt.ylabel('Sensitivity')
plt.savefig('output/figs/sim_ROC_mb.pdf')
plt.close()

baseline_auc = []
noise_auc = []
eigenms_auc = []
limbr_auc = []

for i in range(1,20):
    baseline_auc.append(get_auc('output/simdata/mb_simulated_data_baseline_'+str(i)+'__jtkout_GammaP.txt','output/simdata/mb_simulated_data_key_'+str(i)+'.txt'))
    noise_auc.append(get_auc('output/simdata/simulated_data_with_noise_'+str(i)+'__jtkout_GammaP.txt','output/simdata/mb_simulated_data_key_'+str(i)+'.txt'))
    eigenms_auc.append(get_auc('output/simdata/mb_simdata_eigenMS_'+str(i)+'__jtkout_GammaP.txt','output/simdata/mb_simulated_data_key_'+str(i)+'.txt'))
    limbr_auc.append(get_auc('output/simdata/mb_denoised_circ_lowess_'+str(i)+'__jtkout_GammaP.txt','output/simdata/mb_simulated_data_key_'+str(i)+'.txt'))

auc = pd.DataFrame({'baseline' : baseline_auc, 'noise' : noise_auc, 'eigenms' : eigenms_auc, 'LIMBR' : limbr_auc})

auc = pd.melt(auc)

auc.columns = ['Processing','Area Under the Curve']

sns.boxplot(x="Processing", y="Area Under the Curve", data=auc, order=['baseline','noise','eigenms','LIMBR'], palette=['b','r','y','g'])

plt.savefig('output/figs/sim_AUC_mb.pdf')
plt.close()

#Figure 4A

wt_trends = pd.read_csv('output/actual/wt_lowess_normed_trends.txt',sep='\t')
wt_trends = wt_trends.drop('Unnamed: 0',axis=1)
wt_trends['Trend'] = range(1,len(wt_trends)+1)
wt_trends = pd.melt(wt_trends, id_vars=['Trend'], value_vars=wt_trends.columns.values.tolist()[:-1])

wt_bias_f, wt_bias_icc = get_fstats(wt_trends,wt_class,['TMT Set', 'Prep Set'])

gen_biasmap(wt_bias_icc,wt_bias_f,'output/figs/wt_bcor',r'Correlation of Bias Trends with Groupings in Experimental Procedure (WT)')

#Figure 4B
csp_trends = pd.read_csv('output/actual/csp_lowess_normed_trends.txt',sep='\t')
csp_trends = csp_trends.drop('Unnamed: 0',axis=1)
csp_trends['Trend'] = range(1,len(csp_trends)+1)
csp_trends = pd.melt(csp_trends, id_vars=['Trend'], value_vars=csp_trends.columns.values.tolist()[:-1])

csp_bias_f, csp_bias_icc = get_fstats(csp_trends,csp_class,['TMT Set', 'Prep Set'])

gen_biasmap(csp_bias_icc,csp_bias_f,'output/figs/csp_bcor',r'Correlation of Bias Trends with Groupings in Experimental Procedure ($\Delta$CSP-1)')

#Figure 4C
wt_pepts =  pd.read_csv('output/actual/wt_lowess_normed_pep_bias.txt',sep='\t')
csp_pepts =  pd.read_csv('output/actual/csp_lowess_normed_pep_bias.txt',sep='\t')
merge_pepts = pd.merge(wt_pepts.drop('Protein',axis=1), csp_pepts.drop('Protein',axis=1), on='Peptide').set_index('Peptide')
merge_pepts = scale(merge_pepts)

pep_corr = np.corrcoef(merge_pepts.T)[(len(wt_pepts.columns)-2):((len(wt_pepts.columns)-2)+5),0:5]

plot_cormat_pep(np.abs(pep_corr),'output/figs/peptide_btrend_cormat','Correlation Matrix of Peptide Level Bias Trends')

#Figure 5A
plot_cormat(scale_df(wt),'output/figs/final_cormat','Correlation Matrix of LIMBR Normalized WT Samples Ordered by Timepoint')

prot = pd.merge(wt_p[['GammaP','Phase']], csp_p[['GammaP','Phase']], right_index=True, left_index=True,suffixes=('_wt','_csp'))
prot['zscore_wt'] = stats.norm.ppf(1-prot['GammaP_wt'])
prot['zscore_csp'] = stats.norm.ppf(1-prot['GammaP_csp'])
prot['delta_z'] = prot['zscore_wt'] - prot['zscore_csp']
prot['delta_p'] = prot['Phase_csp'] - prot['Phase_wt']

plot_cormat(scale_df(wt[wt.index.isin(prot[prot['GammaP_wt']<.05].index.tolist())]),'output/figs/final_cormat_circ','Correlation Matrix of LIMBR Normalized WT Circadian Proteins Ordered by Timepoint')

#Figure 5B

defcirc = prot[(prot['GammaP_wt']<.05)&(prot['GammaP_csp']<.05)]

cpw = circular(robjects.FloatVector([j*2*np.pi/22 for j in defcirc['Phase_wt'].values]),units="radians", zero=0, rotation='clock')
cpc = circular(robjects.FloatVector([j*2*np.pi/22 for j in defcirc['Phase_csp'].values]),units="radians", zero=0, rotation='clock')

circorr = corcircular(cpw,cpc,test='TRUE')

phase_hist(defcirc['delta_p'].values,r'Histogram of Phase Lags Between'"\n"r'Expression in WT and $\Delta$CSP-1 Genotypes'"\n"r'For Proteins Circadian in Both Genotypes')

#Figure 1d
to_impute = imputation.imputable("data/Formatted_input.txt",.3)
to_impute.deduplicate()
miss = []
blocks = []
for i in range(1,20):
    miss.append(len((to_impute.data[to_impute.data.columns[0:(i*10)]]).dropna()))
    blocks.append(i)

def fit_exp_linear(t, y, C=0):
    y = y - C
    y = np.log(y)
    K, A_log = np.polyfit(t, y, 1)
    A = np.exp(A_log)
    return A, K

def model_func(t, tau):
    return np.multiply((np.max(miss)-np.min(miss)), np.exp(-t / tau)) + np.min(miss)

def fit_exp_nonlinear(x, y):
    opt_parms, parm_cov = sp.optimize.curve_fit(model_func, x, y, maxfev=10000)
    tau = opt_parms
    return tau

tau = fit_exp_nonlinear(blocks, miss)

fig, ax1 = plt.subplots()
ax1.plot(blocks, [model_func(i, tau) for i in blocks], color='r',ls='--',lw=0.5)[0]
ax1.plot(blocks, miss, color='b', marker='.', lw=0, ms=5)[0]
ax1.set_ylabel('Number of peptides without any missing values', color='k')
ax1.set_xlabel('Number of MS runs', color='k')
ax1.set_xticks(range(1,20))
ax1.set_ylim(bottom=0)
plt.savefig('output/figs/missingness_vs_experiments.pdf')
plt.close()

#Figure S1
noisea = pd.read_csv("grid/melted_aucs.txt", sep='\t').sort_values(by=['Noise'])

fig, ax1 = plt.subplots()
ax1.plot(noisea[noisea['Method'] == 'LIMBR']['Noise'].values, noisea[noisea['Method'] == 'LIMBR']['AUC'].values, color='b', marker='.', lw=0.5, ms=2)[0]
ax1.plot(noisea[noisea['Method'] == 'baseline']['Noise'].values, noisea[noisea['Method'] == 'baseline']['AUC'].values, color='g', marker='.', lw=0.5, ms=2)[0]
ax1.plot(noisea[noisea['Method'] == 'eigenMS']['Noise'].values, noisea[noisea['Method'] == 'eigenMS']['AUC'].values, color='y', marker='.', lw=0.5, ms=2)[0]
ax1.plot(noisea[noisea['Method'] == 'traditional']['Noise'].values, noisea[noisea['Method'] == 'traditional']['AUC'].values, color='r', marker='.', lw=0.5, ms=2)[0]
ax1.set_ylabel('Area Under the Receiver Operating Curve', color='k')
ax1.set_xlabel('Relative Noise Level', color='k')
ax1.set_xticks([0.1,0.5,1,2])
ax1.set_ylim(bottom=0)
plt.savefig('output/figs/AUC_vs_noise.pdf')
plt.close()

#Figure S2
grid = pd.read_csv("grid/grid.txt", sep='\t')

fig, ax1 = plt.subplots()
ax1.plot(grid[grid['NN'] == 5]['Missing'].values, grid[grid['NN'] == 5]['AUC'].values, color='b', marker='.', lw=0.5, ms=2)[0]
ax1.set_ylabel('Area Under the Receiver Operating Curve', color='k')
ax1.set_xlabel('Relative Noise Level', color='k')
ax1.set_xticks([0.2, 0.3, 0.4, 0.5,0.6,0.7])
ax1.set_ylim(bottom=0,top=1)
plt.savefig('output/figs/AUC_vs_NN.pdf')
plt.close()

#Figure ???
#ax = plt.subplot(111)
#hb = ax.hexbin(prot['zscore_wt'].values, prot['zscore_csp'].values, gridsize=48, extent=(-2,4,-2,4))
#ax.set_xlabel(r"WT Z Score", fontsize=12)
#ax.set_ylabel(r'$\Delta$CSP-1 Z Score', fontsize=12)
#cb = fig.colorbar(hb, ax=ax, ticks=np.linspace(0,20,11))
#cb.set_label('count', fontsize=12)
#ax.set_title('Distribution of Z Scores from eJTK of Proteomics Data by Genotype', fontsize=16)
#ax.text(0.015, 0.985, "Pearsons r = %f, p = %e"%(r, p_value), transform=ax.transAxes, fontsize=8, verticalalignment='top', bbox=dict(facecolor='white'))
#plt.savefig('prot_hex.pdf')
#plt.close()

#Figure ???
#sns.jointplot('zscore_wt','zscore_csp',data=prot,kind='kde',joint_kws=dict(bw=.1),marginal_kws=dict(bw=.1))

#Figure ???

max_stars_repo_name: aleccrowell/LIMBR

max_issues_repo_name: aleccrowell/LIMBR

max_forks_repo_name: aleccrowell/LIMBR


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Copyright (c) 2020 by Inria
Authoried by Xiaoyu BIE (xiaoyu.bie@inrai.fr)
License agreement in LICENSE.txt
"""

import os
import sys
import argparse
from tqdm import tqdm
import torch
import numpy as np
import librosa
import soundfile as sf
import matplotlib.pyplot as plt
from dvae.learning_algo import LearningAlgorithm
from dvae.learning_algo_ss import LearningAlgorithm_ss
from dvae.utils.eval_metric import compute_median, EvalMetrics

torch.manual_seed(0)
np.random.seed(0)

class Options:
    def __init__(self):
        self.parser = argparse.ArgumentParser()
        self.opt = None

    def _initial(self):
        # Basic config file
        self.parser.add_argument('--ss', action='store_true', help='schedule sampling')
        self.parser.add_argument('--cfg', type=str, default=None, help='config path')
        self.parser.add_argument('--saved_dict', type=str, default=None, help='trained model dict')
        # Dataset
        self.parser.add_argument('--test_dir', type=str, default='./data/clean_speech/wsj0_si_et_05', help='test dataset')
        # Restuls directory
        self.parser.add_argument('--ret_dir', type=str, default='./data/tmp', help='tmp dir for audio reconstruction')
    def get_params(self):
        self._initial()
        self.opt = self.parser.parse_args()
        params = vars(self.opt)
        return params

params = Options().get_params()

if params['ss']:
    learning_algo = LearningAlgorithm_ss(params=params)
else:
    learning_algo = LearningAlgorithm(params=params)
learning_algo.build_model()
dvae = learning_algo.model
dvae.load_state_dict(torch.load(params['saved_dict'], map_location='cpu'))
eval_metrics = EvalMetrics(metric='all')
dvae.eval()
cfg = learning_algo.cfg
print('Total params: %.2fM' % (sum(p.numel() for p in dvae.parameters()) / 1000000.0))

list_rmse = []
list_sisdr = []
list_pesq = []
list_estoi = []

file_list = librosa.util.find_files(test_dir, ext='wav')

if not os.path.isdir(params['ret_dir']):
    os.makedirs(params['ret_dir'])

for audio_file in tqdm(file_list):

    root, file = os.path.split(audio_file)
    filename, _ = os.path.splitext(file)
    recon_audio = os.path.join(params['ret_dir'], 'recon_{}.wav'.format(filename))
    orig_audio = os.path.join(params['ret_dir'], 'orig_{}.wav'.format(filename))

    # STFT parameters
    wlen_sec = cfg.getfloat('STFT', 'wlen_sec')
    hop_percent = cfg.getfloat('STFT', 'hop_percent')
    fs = cfg.getint('STFT', 'fs')
    zp_percent = cfg.getint('STFT', 'zp_percent')
    wlen = wlen_sec * fs
    wlen = np.int(np.power(2, np.ceil(np.log2(wlen)))) # pwoer of 2
    hop = np.int(hop_percent * wlen)
    nfft = wlen + zp_percent * wlen
    win = np.sin(np.arange(0.5, wlen+0.5) / wlen * np.pi)
    trim = cfg.getboolean('STFT', 'trim')

    x, fs_x = sf.read(audio_file)

    if trim:
        x, _ = librosa.effects.trim(x, top_db=30)

    #####################
    # Scaling on waveform
    #####################
    # scale = 1
    scale = np.max(np.abs(x)) # normalized by Max(|x|)
    x = x / scale

    # STFT
    X = librosa.stft(x, n_fft=nfft, hop_length=hop, win_length=wlen, window=win)

    # Prepare data input
    data_orig = np.abs(X) ** 2 # (x_dim, seq_len)
    data_orig = torch.from_numpy(data_orig.astype(np.float32)).to(dvae.device)
    data_orig = data_orig.permute(1,0).unsqueeze(1) #  (x_dim, seq_len) => (seq_len, 1, x_dim)

    # Reconstruction
    with torch.no_grad():
        data_recon = torch.exp(dvae(data_orig))

    data_recon = data_recon.to('cpu').detach().squeeze().permute(1,0).numpy()

    # Re-synthesis
    X_recon = np.sqrt(data_recon) * np.exp(1j * np.angle(X))
    x_recon = librosa.istft(X_recon, hop_length=hop, win_length=wlen, window=win)

    # Wrtie audio file
    scale_norm = 1 / (np.maximum(np.max(np.abs(x_recon)), np.max(np.abs(x)))) * 0.9
    sf.write(recon_audio, scale_norm*x_recon, fs_x)
    sf.write(orig_audio, scale_norm*x, fs_x)

    rmse, sisdr, pesq, estoi = eval_metrics.eval(audio_est=recon_audio, audio_ref=orig_audio)

    # print('File: {}, rmse: {:.4f}, pesq: {:.4f}, estoi: {:.4f}'.format(filename, rmse, pesq, estoi))

    list_rmse.append(rmse)
    list_sisdr.append(sisdr)
    list_pesq.append(pesq)
    list_estoi.append(estoi)

np_rmse = np.array(list_rmse)
np_sisdr = np.array(list_sisdr)
np_pesq = np.array(list_pesq)
np_estoi = np.array(list_estoi)

print('Re-synthesis finished')
print('RMSE: {:.4f}'.format(np.mean(np_rmse)))
print('SI-SDR: {:.4f}'.format(np.mean(np_sisdr)))
print('PESQ: {:.4f}'.format(np.mean(np_pesq)))
print('ESTOI: {:.4f}'.format(np.mean(np_estoi)))

rmse_median, rmse_ci = compute_median(np_rmse)
sisdr_median, sisdr_ci = compute_median(np_sisdr)
pesq_median, pesq_ci = compute_median(np_pesq)
estoi_median, estoi_ci = compute_median(np_estoi)

print("Median evaluation")
print('median rmse score: {:.4f} +/- {:.4f}'.format(rmse_median, rmse_ci))
print('median sisdr score: {:.4f} +/- {:.4f}'.format(sisdr_median, sisdr_ci))
print('median pesq score: {:.4f} +/- {:.4f}'.format(pesq_median, pesq_ci))
print('median estoi score: {:.4f} +/- {:.4f}'.format(estoi_median, estoi_ci))

max_stars_repo_name: hcost/DVAE

max_issues_repo_name: hcost/DVAE

max_forks_repo_name: hcost/DVAE


text:
import torch as th
import numpy as np
from gluonnlp.torch.layers import SinusoidalPositionalEmbedding

def test_sinusoidal_pos_embed():
    embed1 = SinusoidalPositionalEmbedding(128, learnable=False)
    embed2 = SinusoidalPositionalEmbedding(128, learnable=True)
    assert len([(name, param) for name, param in embed1.named_parameters()
                if param.requires_grad]) == 0
    assert len([(name, param) for name, param in embed2.named_parameters()
                if param.requires_grad]) == 1
    inputs = th.randint(0, 128, (8, 4))
    np.testing.assert_allclose(embed1(inputs).detach().cpu().numpy(),
                               embed2(inputs).detach().cpu().numpy(), 1E-3, 1E-3)

max_stars_repo_name: leezu/gluon-nlp

max_issues_repo_name: leezu/gluon-nlp

max_forks_repo_name: leezu/gluon-nlp


text:
import numpy as np
import PIL.Image
from sklearn.cluster import KMeans
import colorsys

def make_image_array(img):

	# resize image to keep size relatively small and keep calculations running quickly
	basewidth = 100
	wpercent = (basewidth / float(img.size[0]))
	hsize = int((float(img.size[1]) * float(wpercent)))
	img = img.resize((basewidth, hsize), PIL.Image.ANTIALIAS)

	# get channels and turn colours into dataframe columns
	arr = np.array(img.convert('RGB')) # need .convert('RGB') as sometimes png may not be loaded as such

	r = arr[:,:,0].flatten()
	g = arr[:,:,1].flatten()
	b = arr[:,:,2].flatten()

	rgb = np.column_stack((r,g,b))

	return rgb

def prettify_colours(centroids):

	# for a more visually appealing palette, let's order the colours by hue
	colours = []
	for row in range(len(centroids)):

	    c = list(centroids[row,:]/255)

	    colours.append(list(colorsys.rgb_to_hsv(*c)))

	colours.sort()

	colours = [list(colorsys.hsv_to_rgb(*c)) for c in colours]

	return colours

def kmeans(img, n_colours=10):

	arr = make_image_array(img=img)

	# pick how many colours will populate our palette, then perform a KMeans, and finally get the centroids identifying the colours we will display
	n_clusters = n_colours

	kmeans = KMeans(n_clusters=n_clusters).fit(arr)

	centroids = kmeans.cluster_centers_

	colours = prettify_colours(centroids=centroids)

	return colours

max_stars_repo_name: pguarnieri/get-palette-app

max_issues_repo_name: pguarnieri/get-palette-app

max_forks_repo_name: pguarnieri/get-palette-app


text:
"""
Fibonacci Spreadsheet example for book.

Depends on having tkinter installed.
"""
try:
    import networkx as nx
except ImportError:
    import ch07.replacement as nx

from ch07.dependencies import tkinter_error

def fibonacci_example(spreadsheet):
    """Initialize Spreadsheet to small Fibonacci example for book."""
    entries = {
        'A1': 'N',
        'B1': 'Fibn',
        'C1': 'Sum',
        'A2': '0',
        'A3': '=(A2+1)',
        'A4': '=(A3+1)',
        'A5': '=(A4+1)',
#         'A6': '=(A5+1)',         ## Add these lines back in for full
#         'A7': '=(A6+1)',         ## example.
#         'A8': '=(A7+1)',
        'B2': '0',
        'B3': '1',
        'B4': '=(B2+B3)',
        'B5': '=(B3+B4)',
#         'B6': '=(B4+B5)',
#         'B7': '=(B5+B6)',
#         'B8': '=(B6+B7)',
        'C2': '=B2',
        'C3': '=(B3+C2)',
        'C4': '=(B4+C3)',
        'C5': '=(B5+C4)',
#         'C6': '=(B6+C5)',
#         'C7': '=(B7+C6)',
#         'C8': '=(B8+C7)'
    }
    for k in entries:
        spreadsheet.set(k,entries[k])

#######################################################################
if __name__ == '__main__':
    if tkinter_error:
        print('tkinter is not installed so unable to launch spreadsheet application')
    else:
        import tkinter
        from ch07.spreadsheet import Spreadsheet

        root = tkinter.Tk()
        ss = Spreadsheet(root, nx.DiGraph())
        fibonacci_example(ss)
        from ch07.digraph_search import topological_sort
        print(list(topological_sort(ss.digraph)))
        root.mainloop()

max_stars_repo_name: laszlokiraly/LearningAlgorithms

max_issues_repo_name: laszlokiraly/LearningAlgorithms

max_forks_repo_name: laszlokiraly/LearningAlgorithms

reason: import networkx


text:
import numpy as np
import pickle

"""
This class reads in the data necessary for the algorithm.

The theoretical frequency spectra need to be dumped to a
pickle file and must have the dimension (n_conformer,n_freq,2),
    where the first dimension is the number of conformer
    the second dimension is the number of normal modes, and
    the third dimension is the frequency (0) and dipol moment (1) OR ir intensity.
It must be specified, whether the IR intensity is saved or the dipole moment (default is IR intensity)

The energies of the conformer also must be dumped to a pickle file and must have the dimension (n_conformer).
The default energy assumed is Hartree (default in most QM programs), and is converted to kJ/mol.

"""
class Reader:
    def __init__(self):
        print("Reader init successfully")
    def read_pickle(self,file_name_E,file_name_f,threshold = 20):
        self.E = np.asarray(pickle.load(open(file_name_E,'rb')))
        self.freq = np.asarray(pickle.load(open(file_name_f,'rb')))
        ## SORT ENERGIES AND FREQUENCIES
        idx = np.argsort(self.E)
        self.E = self.E[idx]
        self.E -= np.min(self.E)
        self.freq = self.freq[idx]
        ## CONVERT HARTREE TO KJ/MOL
        self.E *= 2625.50
        #tmp = self.E <= threshold
        ##Only consider contributing conformers
        #self.E = self.E[0:1]
        #self.freq = self.freq[0:1]
    def get_freq(self):
        return self.freq
    def get_energy(self):
        return self.E
    def write(self,out,p,returnvalue,args):
        f_open = open(out+"results.txt","w+")
        f_open.write("PEARSON "+str(p)+"\n")
        f_open.write("ALIGNMENT_SCORE "+str(-returnvalue)+"\n")
        f_open.write("PRODUCT "+str(-returnvalue*p)+"\n")
        f_open.write("------------SETTINGS USED---------------\n")
        f_open.write(str(args))
        f_open.close()


text:
import numpy
import pandas as pd
from keras.models import Sequential
from keras.layers import Dense
from keras.wrappers.scikit_learn import KerasClassifier
from sklearn.model_selection import cross_val_score, cross_val_predict
from sklearn.preprocessing import LabelEncoder
from sklearn.model_selection import StratifiedKFold
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import confusion_matrix

from micro_influencer_utilities import *

# fix random seed for reproducibility
seed = 7
numpy.random.seed(seed)

#Reading the data set
pathToDataFolder = "Data"
parametersPath = "/03_users_parameters/"

if len(sys.argv)== 2:
	topic_selected = sys.argv[1]
	if not topic_selected.startswith('#'):
		topic_selected = "#"+topic_selected
else:
	topic_selected = topic_selection()
pathToTopic = pathToDataFolder+"/"+topic_selected
pathToUserParameters = pathToTopic + parametersPath

#filtering out not confoming data i.e max-value (oor) or no words in category (anomaly)
f = open(pathToUserParameters+"/table/"+topic_selected[1:]+".csv", "r")
f2 = open(pathToUserParameters+"/table/"+topic_selected[1:]+"2.csv", "w")
for line in f.readlines():
	if "max" in line:
		continue
	elements = line.split(",")
	if len(elements)<3040:
		continue
	else:
		f2.write(line)
f.close()
f2.close()

#read csv from the filtered and clean csv(just created), not the old one (dirt one)
dataset = pd.read_csv(pathToUserParameters+"/table/"+topic_selected[1:]+"2.csv")

# load dataset
# split into input (X) and output (Y) variables

#Pre-processing data set to data frame by choosing features (X) and output (y)
#print (dataset.head()) #debug
Y = dataset.y #micro influencer yes = 1, no = 0
X = dataset.drop(["user_screen_name", "Semb", "Srec", "Sint","y"],axis=1)# "num_of_words","distance", "y"], axis=1)

X = X.loc[:, ["big5_O", "big5_C", "big5_E", "big5_A", "big5_N","scoreselfdirection","scorestimulation", "scorehedonism",
"scoreachievement","scorepower" ,"scoresecurity","scoreconformity", "scoretradition", "scorebenevolence", "scoreuniversalism"]].astype(float)

# encode class values as integers
encoder = LabelEncoder()
encoder.fit(Y)
encoded_Y = encoder.transform(Y)

# baseline model
def create_baseline():
	# create model
	model = Sequential()
	model.add(Dense(15, input_dim=15, kernel_initializer='normal', activation='relu'))
	model.add(Dense(1, kernel_initializer='normal', activation='sigmoid'))
	# Compile model
	model.compile(loss='binary_crossentropy', optimizer='adam', metrics=['accuracy'])
	return model

# evaluate baseline model with standardized dataset
numpy.random.seed(seed)
estimators = []
estimators.append(('standardize', StandardScaler()))
estimators.append(('mlp', KerasClassifier(build_fn=create_baseline, epochs=100, batch_size=5, verbose=0)))
pipeline = Pipeline(estimators)
kfold = StratifiedKFold(n_splits=5, shuffle=True, random_state=seed)
results = cross_val_score(pipeline, X, encoded_Y, cv=kfold)
print("Standardized: %.2f%% (%.2f%%)" % (results.mean()*100, results.std()*100))

y_pred = cross_val_predict(pipeline, X, encoded_Y, cv=kfold)
y_pred = numpy.asarray(y_pred)
#print("y_pred", y_pred)
#print("y_true", encoded_Y)

tn, fp, fn, tp = confusion_matrix(encoded_Y, y_pred).ravel();
precision = tp/(tp+fp);
recall = tp/(tp+fn)
f1_score = 2*(precision*recall)/(precision+recall);
print("precision: ", precision)
print("recall: ", recall)
print("f1_score:", f1_score)
# import keras
# import keras_metrics as km
# import sklearn.metrics as sklm

# class Metrics(keras.callbacks.Callback):
#     def on_train_begin(self, logs={}):
#         self.confusion = []
#         self.precision = []
#         self.recall = []
#         self.f1s = []
#         self.kappa = []
#         self.auc = []

#     def on_epoch_end(self, epoch, logs={}):
#         score = numpy.asarray(self.model.predict(self.validation_data[0]))
#         predict = numpy.round(numpy.asarray(self.model.predict(self.validation_data[0])))
#         targ = self.validation_data[1]

#         self.auc.append(sklm.roc_auc_score(targ, score))
#         self.confusion.append(sklm.confusion_matrix(targ, predict))
#         self.precision.append(sklm.precision_score(targ, predict))
#         self.recall.append(sklm.recall_score(targ, predict))
#         self.f1s.append(sklm.f1_score(targ, predict))
#         self.kappa.append(sklm.cohen_kappa_score(targ, predict))

#         return

# #-----------------following doesn't work------------------
# def second_baseline():
# 	# create model
# 	model = Sequential()
# 	model.add(keras.layers.Dense(15, activation="sigmoid", input_dim=15))
# 	model.add(keras.layers.Dense(1, activation="softmax"))

# 	# Calculate precision for the label 1.
# 	precision = km.binary_precision(label=1)

# 	# Calculate recall for the label 1.
# 	recall = km.binary_recall(label=1)
# 	# Compile model
# 	model.compile(optimizer="sgd", loss="binary_crossentropy", metrics=[precision, recall])
# 	return model

# sec_estimators = []
# sec_estimators.append(('standardize', StandardScaler()))
# sec_estimators.append(('mlp', KerasClassifier(build_fn=second_baseline, epochs=100, batch_size=5, verbose=0)))
# sec_pipeline = Pipeline(sec_estimators)
# sec_kfold = StratifiedKFold(n_splits=10, shuffle=True, random_state=seed)
# sec_results = cross_val_score(sec_pipeline, X, encoded_Y, cv=sec_kfold)
# #print("Standardized: %.2f%% (%.2f%%)" % (results.mean()*100, results.std()*100))
# print("sec_results", sec_results)


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
@author: Frederik Kratzert. frederik.kratzert(at)boku.ac.at
"""

import os
import tensorflow as tf
import numpy as np

valid_modes = ['vgg', 'length', 'date', 'all']

class VGG16:
    """
    This class stores the model definition of the original VGG16 network,
    as well as the adapted versions for additional feature fusing into the
    second-to-last layer.
    """

    def __init__(self, mode = 'vgg', num_classes=7, keep_prob = 0.5):

        """
        Arguments:

          mode: one of ['vgg', 'length', 'date', 'all']. Creates the model
            graphs for the four different model versions.
          num_classes: number of classes. 7 for this experiment, 1000 for the
            original vgg16 model
          keep_prob: dropout probability. 0.5 is the value of the original model

        """

        if mode not in valid_modes:
            raise ValueError("Bad input argument. 'Mode' must be one of 'vgg', 'length', 'date' or 'all'.")

        self.mode = mode
        self.num_classes = num_classes
        self.keep_prob = keep_prob

    def build(self, images, length = None, weeks = None):
        """
        This function creates the model graph depending on the selected 'mode'.

        Inputs:
            - images: tf.placeholder of size (N, 224, 224, 3) with N number of
                      images per batch.
            - length: (optional) tf.placeholder of size (N,1)
            - weeks: (optional) tf.placeholder of size (N, 52)
        """

        # First convolutional block of VGG16
        conv1_1 = self.conv_layer(images, 3, 3, 64, 'conv1_1')
        conv1_2 = self.conv_layer(conv1_1, 3, 3, 64, 'conv1_2')
        pool1 = self.max_pool(conv1_2, 'pool1')

        # Second convolutional block of VGG16
        conv2_1 = self.conv_layer(pool1, 3, 3, 128, 'conv2_1')
        conv2_2 = self.conv_layer(conv2_1, 3, 3, 128, 'conv2_2')
        pool2 = self.max_pool(conv2_2, 'pool2')

        # Third convolutional block of VGG16
        conv3_1 = self.conv_layer(pool2, 3, 3, 256, 'conv3_1')
        conv3_2 = self.conv_layer(conv3_1, 3, 3, 256, 'conv3_2')
        conv3_3 = self.conv_layer(conv3_2, 3, 3, 256, 'conv3_3')
        pool3 = self.max_pool(conv3_3, 'pool3')

        # Fourth convolutional block of VGG16
        conv4_1 = self.conv_layer(pool3, 3, 3, 512, 'conv4_1')
        conv4_2 = self.conv_layer(conv4_1, 3, 3, 512, 'conv4_2')
        conv4_3 = self.conv_layer(conv4_2, 3, 3, 512, 'conv4_3')
        pool4 = self.max_pool(conv4_3, 'pool4')

        # Fifth convolutional block of VGG16
        conv5_1 = self.conv_layer(pool4, 3, 3, 512, 'conv5_1')
        conv5_2 = self.conv_layer(conv5_1, 3, 3, 512, 'conv5_2')
        conv5_3 = self.conv_layer(conv5_2, 3, 3, 512, 'conv5_3')
        pool5 = self.max_pool(conv5_3, 'pool5')

        # flatten pool5 into one long vector
        pool5_flat = tf.reshape(pool5, [-1, 512*7*7])

        # first fully connected layer of original VGG16 model + dropout
        fc6 = self.fc_layer(pool5_flat, 512*7*7, 4096, name = 'fc6')
        drpt6 = tf.nn.dropout(fc6, self.keep_prob)

        #apply batch norm
        bn6 = self.batch_norm(drpt6, 'bn6')

        # Depending on mode, select which top to but on the convolutional core
        if self.mode == 'vgg':

            self.fc8 = self.vgg_top(bn6)

        elif self.mode == 'length':

            self.fc8 = self.length_top(bn6, length)

        elif self.mode == 'date':

            self.fc8 = self.weeks_top(bn6, weeks)

        elif self.mode == 'all':

            self.fc8 = self.all_top(bn6, length, weeks)

    def conv_layer(self, x, filter_height, filter_width, num_filters, name):
        """
        Helper class to create a convolutional layer of arbitrary filter
        dimensions, input channels and number of filter
        """
        #number of channels of the input layer
        in_channels = int(x.shape[-1].value)

        with tf.variable_scope(name) as scope:

            #create variables for weights and biases
            weights = tf.get_variable('weights', shape = [filter_height, filter_width,
                                                          in_channels, num_filters])
            biases = tf.get_variable('biases', shape = [num_filters])

            #perform convolution, bias add and relu activation
            conv = tf.nn.conv2d(x, weights, strides = [1,1,1,1], padding = 'SAME')
            bias = tf.nn.bias_add(conv, biases)
            relu = tf.nn.relu(bias, name = scope.name)

            return relu

    def fc_layer(self, x, num_in, num_out, name, relu=True):
        """
        Helper class to create a fully connected layer of arbitrary number of
        input neurons and output neurons
        """
        with tf.variable_scope(name) as scope:

            #create variables for weights and biases
            weights = tf.get_variable('weights', shape=[num_in, num_out])
            biases = tf.get_variable('biases', shape=[num_out])

            #calculate activation
            act = tf.nn.xw_plus_b(x, weights, biases, name=scope.name)

            #apply non-linearity depending on set relu flag
            if relu == True:
                return tf.nn.relu(act)
            else:
                return act

    def max_pool(self, x, name):
        """
        Helper function for the max pooling layer
        """
        return tf.nn.max_pool(x, ksize = [1, 2, 2, 1], strides = [1, 2, 2, 1],
                            padding = 'SAME', name = name)

    def batch_norm(self, x, name, is_train = False):
        """
        Helper function for the batch_norm layer. Note that for the sake of
        testing the 'is_train' flag needs to be set to false. If you are
        interested in training a model with this function, you need to pass
        a 'is_train' = True during the trainings phase to correctly update the
        moving mean/std.
        """
        with tf.variable_scope(name) as scope:

            out = tf.contrib.layers.batch_norm(x, decay = 0.9, center = True,
                                            scale = True, updates_collections = None,
                                            is_training = is_train, fused = True,
                                            scope = scope.name)
            return out

    def vgg_top(self, bn6):
        """
        Helper function to create the upper part of the original VGG16 model
        starting from the fc7 layer.
        """
        #second fully connected layer (fc7) + batchnorm + dropout
        fc7 = self.fc_layer(bn6, 4096, 4096, name = 'fc7')
        bn7 = self.batch_norm(fc7, name = 'fc7')
        drpt7 = tf.nn.dropout(bn7, self.keep_prob)

        #output layer with self.num_class neurons and without non-linearity
        fc8 = self.fc_layer(drpt7, 4096, self.num_classes, name = 'fc8', relu = False)

        return fc8

    def length_top(self, bn6, length):
        """
        Helper function to create the upper part of the model, were the length
        feature is fused in as additional input.
        """
        #process length through seperate batch norm layer
        bn_length = self.batch_norm(length, 'input_length')

        #fuse image features and length feature together to one long vector
        fused_inputs = tf.concat([bn6, bn_length], axis = 1)

        #add fully connected layer (fc7) + batchnorm + dropout
        fc7 = self.fc_layer(fused_inputs, 4097, 4096, name = 'fc7')
        bn7 = self.batch_norm(fc7, name = 'fc7')
        drpt7 = tf.nn.dropout(bn7, self.keep_prob)

        #output layer with self.num_class neurons and without non-linearity
        fc8 = self.fc_layer(drpt7, 4096, self.num_classes, name = 'fc8', relu = False)

        return fc8

    def weeks_top(self, bn6, weeks):
        """
        Helper function to create the upper part of the model, were the one-hot
        week vector serves as additional feature.
        """

        #fuse week vector together with the bn6 output vector
        fused_inputs = tf.concat([bn6, weeks], axis = 1)

        #add fully connected layer (fc7) + batchnorm + dropout
        fc7 = self.fc_layer(fused_inputs, 4148, 4096, name = 'fc7')
        bn7 = self.batch_norm(fc7, name = 'fc7')
        drpt7 = tf.nn.dropout(bn7, self.keep_prob)

        #output layer with self.num_class neurons and without non-linearity
        fc8 = self.fc_layer(drpt7, 4096, self.num_classes, name = 'fc8', relu = False)

        return fc8

    def all_top(self, bn6, length, weeks):
        """
        Helper function to create the upper part of the model, were the one-hot
        week vector and the length feature are fused into the model.
        """
        #process length through seperate batch norm layer
        bn_length = self.batch_norm(length, 'input_length')

        #fuse image features and length feature together to one long vector
        fused_inputs = tf.concat([bn6, bn_length, weeks], axis = 1)

        #add fully connected layer (fc7) + batchnorm + dropout
        fc7 = self.fc_layer(fused_inputs, 4149, 4096, name = 'fc7')
        bn7 = self.batch_norm(fc7, name = 'fc7')
        drpt7 = tf.nn.dropout(bn7, self.keep_prob)

        #output layer with self.num_class neurons and without non-linearity
        fc8 = self.fc_layer(drpt7, 4096, self.num_classes, name = 'fc8', relu = False)

        return fc8


text:
# This script will generate the fingerprints for the mol objects
import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem import DataStructs

######################## generate_ftps() ##################################
#INPUT:
# an array or list of mol objects
#OUTPUT:
# a list of fingerprints (each fingerprint is an array)
# a list of corresponding SMILE strings (a list of char strings)
def generate_fpts(molList):
    pos_at=[1,6,7,8,9,11,15,16,17,19,35,53]
    pos_charges=[-1,0,1]
    # the array for the fingerprints
    fpt_arr=[]
    # The corresponding array for the smile strings
    smile_list=[]
    # mol_withH coresponds to moll, also corresponds to moli
    # Looks like you will need to iterate through molList witha for loop and call Chem.addHs(molList[i]), for instance
    for i in np.arange(len(molList)): # iterate through molList
        moli_withH=Chem.AddHs(molList[i])
        smile=Chem.MolToSmiles(moli_withH)
        # NOTE: some objects will be removed by this filter, so fpt_array and smile_list will be smaller than molList
        if moli_withH.GetNumAtoms()<=80 and all([x.GetAtomicNum() in pos_at for x in moli_withH.GetAtoms()]) and all([x.GetFormalCharge() in pos_charges for x in moli_withH.GetAtoms()]):
            # creating an empty numpy array
            arr = np.zeros((0,), dtype=np.float32)
            # Creating the fingerprints and also converting them to a numpy array, storing in arr5
            DataStructs.ConvertToNumpyArray(AllChem.GetMorganFingerprintAsBitVect(moli_withH,4, nBits=1024),arr)
            # adding the fingerprints (as np arrays) to list
            fpt_arr.append(arr)
            # creating a corresponding list of smile strings for the mols
            smile_list.append(smile)
    return (fpt_arr,smile_list)

# loading in the array of mol objects as generated by download_pubchem.py -> for Data Science Project
theMols=np.load("CS_numpy_objs/theMols.npy",allow_pickle=True)

# running the function on the pubChem data
fingerprints, smiles = generate_fpts(theMols)

################## saving the objects for later use ###############################
np.save("CS_numpy_objs/fingerprints.npy",fingerprints,allow_pickle=True)
np.save("CS_numpy_objs/smiles.npy",smiles,allow_pickle=True)


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Aug 13 16:19:11 2018

@author: madeleineernst
"""
import scipy
from scipy.cluster import hierarchy
import numpy as np
import pandas as pd
import obonet
import networkx
import os
import pkg_resources

# from https://stackoverflow.com/questions/28222179/save-dendrogram-to-newick-format
def getNewick(node, newick, parentdist, leaf_names):
    if node.is_leaf():
        return "%s:%.2f%s" % (leaf_names[node.id], parentdist - node.dist, newick)
    else:
        if len(newick) > 0:
            newick = "):%.2f%s" % (parentdist - node.dist, newick)
        else:
            newick = ");"
        newick = getNewick(node.get_left(), newick, node.dist, leaf_names)
        newick = getNewick(node.get_right(), ",%s" % (newick), node.dist, leaf_names)
        newick = "(%s" % (newick)
        return newick

def get_motiftrees(motifs, buckettable, method = "ward", metric = "euclidean", outputdir = "MotifTree/"):

    os.mkdir(outputdir)

    # select only features present in tree
    bt_sel = buckettable[buckettable['#OTU ID'].isin(list(set(buckettable['#OTU ID']) & set(motifs['scans'])))]
    bt_sel.to_csv(outputdir + 'Buckettable_Motifs.tsv', sep = '\t', index = False)

    motifs = motifs[pd.notnull(motifs['motif'])] # remove features, which do not contain any motif
    motifs = motifs.loc[:, (motifs != 0).any(axis=0)]
    motifs.index = motifs['scans']
    motifs = motifs.filter(like='motif_') # select all motif columns

    Z = scipy.cluster.hierarchy.linkage(motifs, method=method, metric=metric)
    leaf_names = motifs.index # remove white space from leaf labels
    tree = hierarchy.to_tree(Z,False)

    f = open(outputdir + 'Tree_Motifs.txt','w')
    f.write(getNewick(tree, "", tree.dist, leaf_names))
    f.close()

def get_classytrees(cf,buckettable, lev, method='average', metric='jaccard', outputdir = "ClassyTree/"):

    os.makedirs(outputdir, exist_ok=True)

    # remove inorganic compounds by default
    cf = cf[cf.CF_kingdom != 'Inorganic compounds']
    cf = cf.rename(columns = {'cluster index':'cluster.index'})

    graph_path = pkg_resources.resource_filename('MetaboDistTrees', 'data/ChemOnt_2_1.obo')
    graph = obonet.read_obo(graph_path)

    # Mapping from term ID to name
    id_to_name = {id_: data['name'] for id_, data in graph.nodes(data=True)}
    name_to_id = {data['name']: id_ for id_, data in graph.nodes(data=True) if 'name' in data}

    for i in range(len(lev)):

        if lev[i] == 'cluster.index':

            all_subclasses = list(set(cf['CF_Dparent']))
            all_subclasses = [x for x in all_subclasses  if str(x) != 'nan']
            all_subclasses = [x for x in all_subclasses  if str(x) != 'no matches']
            all_subclasses = [x for x in all_subclasses  if str(x) != ' ']
            all_subclasses = [x for x in all_subclasses  if str(x) != 'Non-metal phosphonates']

            a = []
            for j in all_subclasses:
                paths = networkx.all_simple_paths(
                        graph,source=name_to_id[j],
                        target=name_to_id['Organic compounds'])
                for path in paths:
                    a.append('#'.join(id_to_name[node] for node in path))

            notintree = list(set(cf['CF_Dparent']) - set(all_subclasses))
            cf = cf[~cf['CF_Dparent'].isin(notintree)]
            #b = list(set(cf[lev[i]]))
            b = list(cf['CF_Dparent'])
            for index,item in enumerate(b):
                b[index] = '#'.join([str(list(cf[lev[i]])[index]) ,a[all_subclasses.index(item)]])

            inames = []
            for index, item in enumerate(b):
                inames.append(b[index].split('#', 1)[0])
                b[index] = b[index].split('#', 1)[-1]

            df = pd.DataFrame(b)
            df = df.iloc[:,0].str.get_dummies(sep='#')
            df.index = inames

            df.to_csv(os.path.join(outputdir, 'Classlist_' + lev[i] + '.tsv'),sep='\t',index=True)

            # do hierarchical cluster analysis on classlist and convert to newick format
            Z = scipy.cluster.hierarchy.linkage(df, method=method, metric=metric)

            leaf_names = df.index # remove white space from leaf labels
            tree = hierarchy.to_tree(Z,False)

            f = open(os.path.join(outputdir, 'NewickTree_' + lev[i] + '.txt'),'w')
            f.write(getNewick(tree, "", tree.dist, leaf_names))
            f.close()

        else:
            all_subclasses = list(set(cf[lev[i]]))
            all_subclasses = [x for x in all_subclasses  if str(x) != 'nan']
            all_subclasses = [x for x in all_subclasses  if str(x) != 'no matches']
            all_subclasses = [x for x in all_subclasses  if str(x) != ' ']

            a = []
            for j in all_subclasses:
                paths = networkx.all_simple_paths(
                        graph,source=name_to_id[j],
                        target=name_to_id['Organic compounds'])
                for path in paths:
                    a.append('#'.join(id_to_name[node] for node in path))

            inames = []
            for index, item in enumerate(a):
                inames.append(a[index].split('#', 1)[0])
                a[index] = a[index].split('#', 1)[-1]

            df = pd.DataFrame(a)
            df = df.iloc[:,0].str.get_dummies(sep='#')
            df.index = [num.replace(' ', '.') for num in inames]
            df.index = [num.replace(',', '.') for num in df.index]
            df.index = [num.replace('(', '.') for num in df.index]
            df.index = [num.replace(')', '.') for num in df.index]
            df.index = [num.replace("'", '') for num in df.index]
            df.index = [num.replace('[', '.') for num in df.index]
            df.index = [num.replace(']', '.') for num in df.index]

            df.to_csv(os.path.join(outputdir, 'Classlist_' + lev[i] + '.tsv'),sep='\t',index=True)

            # do hierarchical cluster analysis on classlist and convert to newick format
            Z = scipy.cluster.hierarchy.linkage(df, method='average', metric='jaccard')

            leaf_names = df.index # remove white space from leaf labels
            tree = hierarchy.to_tree(Z,False)

            f = open(os.path.join(outputdir, 'NewickTree_' + lev[i] + '.txt'),'w')
            f.write(getNewick(tree, "", tree.dist, leaf_names))
            f.close()

    # select only features present in tree
    bt_sel = buckettable[buckettable['#OTU ID'].isin(list(set(buckettable['#OTU ID']) & set(cf['cluster.index'])))]
    bt_sel.to_csv(os.path.join(outputdir, 'Buckettable_ChemicalClasses.tsv'), sep = '\t', index = False)

reason: import numpy,import scipy,from scipy,import networkx