text:
import numpy as np
import pandas as pd
import scipy.sparse as sparse
import random
import copy
from multiprocessing import Pool
from functools import partial
from sklearn.metrics.pairwise import cosine_similarity

# метрика на 1 карту
def fuzzy_crd(rec_ids, fact_ids, item_factors, plu_category_dict, weight_by_price=True, recs_prices_list=None):
    if (len(rec_ids) != 0) & (len(fact_ids) != 0):
        similarity_matrix = get_cosine_similarity(fact_ids, rec_ids, item_factors)
        heaviside_matrix = get_heaviside_matrix(fact_ids, rec_ids, plu_category_dict)
        fzy_col = np.amax(similarity_matrix * heaviside_matrix.T, axis=1)

        if weight_by_price:
            fzy_col_w = [f*price for f, price in zip(fzy_col, recs_prices_list)]
            fzy = sum(fzy_col_w) / sum(recs_prices_list)
        else:
            fzy = sum(fzy_col) / len(fzy_col)
        return fzy
    else:
        return np.nan

def precision_crd(rec_ids, fact_ids, weight_by_price=True, recs_prices_list=None):
    if (len(rec_ids) != 0) & (len(fact_ids) != 0):
        if weight_by_price:
            guessed = [rec in fact_ids for rec in rec_ids]
            guessed_w = [rec*price for rec, price in zip(guessed, recs_prices_list)]
            precision_score = sum(guessed_w) / sum(recs_prices_list)
        else:
            guessed = [rec in fact_ids for rec in rec_ids]
            precision_score = sum(guessed) / len(fact_ids)
        return precision_score
    else:
        return np.nan

def get_heaviside_matrix(fact_ids, rec_ids, plu_category_dict):
    heaviside_matrix = []
    for fct in fact_ids:
        dim = []
        for rcm in rec_ids:
            dim.append(heaviside(plu_category_dict[fct], plu_category_dict[rcm]))
        heaviside_matrix.append(dim)
    return np.array(heaviside_matrix)

def heaviside(cat1, cat2):
    if cat1 == cat2:
        return 1
    else:
        return 0

def get_cosine_similarity(id1, id2, item_factors):
    factors_1_full = []
    factors_2_full = []
    for id_ in id1:
        factors_1 = item_factors.loc[id_].tolist()
        factors_1_full.append(factors_1)
    for id_ in id2:
        factors_2 = item_factors.loc[id_].tolist()
        factors_2_full.append(factors_2)
    return cosine_similarity(factors_2_full, factors_1_full)

def get_recommendations(model, user_items_sparse_matrix, index_crd_dict, filter_already_liked, num_recommendations):
    recs = model.recommend_all(user_items=sparse.csr_matrix(user_items_sparse_matrix), show_progress=False,
                               filter_already_liked_items=filter_already_liked, N=num_recommendations)
    recs = pd.DataFrame(recs).reset_index().rename(columns={'index': 'crd_no'})
    recs['crd_no'] = recs['crd_no'].map(lambda x: index_crd_dict[x])
    recs['recs'] = recs.apply(lambda row_df: [int(row_df[x]) for x in range(num_recommendations)], axis=1)
    recs_dict = recs[['crd_no', 'recs']].set_index('crd_no').to_dict()['recs']
    return recs_dict

# метрики на все карты
def fuzzy(user_item_weighted_matrix, model, facts_dict, plu_category_dict, method='mean',
          num_recommendations=5, num_threads=16,
          filter_already_liked=False,
          weight_by_price=True, plu_price=None,
          score_baseline=None):
    """
    Calculates fuzzy precision (cosine_similarity*presence_in_same_category/num_predictions)
    Args:
        user_item_weighted_matrix: pd.core.frame.DataFrame user_item matrix
        model: implicit.als.AlternatingLeastSquares fitted instance or similar
        facts_dict: dictionary of fact sells: {'crd_no': '[plu_1, plu_2, ..., plu_n]'}
        plu_category_dict: dictionary {'plu': 'категория_lvl_2'}
        method: 'mean' or 'max' closeness of recommended item embeddings to true (bought) item embeddings
        num_recommendations: num items to recommend while checking fuzzy
        filter_already_liked: boolean for filtering already presented items
        weight_by_price: boolean. If you need to weight prices.
        plu_price: dict. Dict format {'plu':'price'}
        score_baseline:
        num_threads: Num threads to parallel
    Returns:
        Fuzzy score
    """
    index_crd_dict = user_item_weighted_matrix.reset_index()['index'].to_dict()
    item_factors = pd.DataFrame(model.item_factors, index=user_item_weighted_matrix.columns.tolist())
    user_items_sparse_matrix = sparse.csr_matrix(user_item_weighted_matrix)
    recs_dict = get_recommendations(model, user_items_sparse_matrix, index_crd_dict,
                                    filter_already_liked, num_recommendations)
    id_plu_dict = user_item_weighted_matrix.T.reset_index()['index'].to_dict()
    p = Pool(num_threads)
    facts_split = np.array_split(np.array(list(facts_dict.keys())), num_threads)
    fz_list_of_lists = p.map(partial(_fuzzy_iter,
                                     recs_dict=recs_dict, plu_category_dict=plu_category_dict, facts_dict=facts_dict,
                                     item_factors=item_factors, id_plu_dict=id_plu_dict,
                                     weight_by_price=weight_by_price, plu_price=plu_price,
                                     score_baseline=score_baseline), facts_split)
    if method == 'mean':
        return np.mean([item for sublist in fz_list_of_lists for item in sublist])
    elif method == 'max':
        return np.max([item for sublist in fz_list_of_lists for item in sublist])
    else:
        raise NotImplementedError('Method: %s Not implemented' % method)

def _fuzzy_iter(facts_subset, facts_dict, recs_dict, plu_category_dict,
                item_factors, weight_by_price, id_plu_dict,
                plu_price, score_baseline):
    fz_list = []
    recs_prices_list = []
    for crd in facts_subset:
        recs_list = recs_dict[crd]
        # бейзлайн возвращает сразу plu_id. Поэтому для бейзлайна эта функция не применяется
        if not score_baseline:
            recs_list = [id_plu_dict[x] for x in recs_list]
        if (len(recs_list) > 0) & weight_by_price:
            recs_prices_list = [plu_price[str(rec)] for rec in recs_list]
        facts_list = facts_dict[crd]
        fz = fuzzy_crd(recs_list, facts_list, item_factors,
                       plu_category_dict, weight_by_price, recs_prices_list)
        fz_list.append(fz)
    return fz_list

def precision(user_item_weighted_matrix, model, facts_dict,
              num_recommendations=5,
              filter_already_liked=False, weight_by_price=True, plu_price=None,
              score_baseline=None):
    """
    Calculates fuzzy precision (cosine_similarity*presence_in_same_category/num_predictions)
    Args:
        user_item_weighted_matrix: pd.core.frame.DataFrame user_item matrix
        model: implicit.als.AlternatingLeastSquares fitted instance or similar
        facts_dict: dictionary of fact sells: {'crd_no': '[plu_1, plu_2, ..., plu_n]'}
        num_recommendations: num items to recommend while checking fuzzy
        filter_already_liked: boolean for filtering already presented items
        weight_by_price: boolean. If you need to weight prices.
        plu_price: dict. Dict format {'plu':'price'}
        score_baseline:
    Returns:
        Fuzzy score
    """
    index_crd_dict = user_item_weighted_matrix.reset_index()['index'].to_dict()
    user_items_sparse_matrix = sparse.csr_matrix(user_item_weighted_matrix)
    recs_dict = get_recommendations(model, user_items_sparse_matrix, index_crd_dict,
                                    filter_already_liked, num_recommendations)
    id_plu_dict = user_item_weighted_matrix.T.reset_index()['index'].to_dict()

    recs_prices_list = []
    prec_list = []
    for crd in facts_dict.keys():
        recs_list = recs_dict[crd]
        # бейзлайн возвращает сразу plu_id. Поэтому для бейзлайна эта функция не применяется
        if not score_baseline:
            recs_list = [id_plu_dict[x] for x in recs_list]
        if (len(recs_list) > 0) & weight_by_price:
            recs_prices_list = [plu_price[str(rec)] for rec in recs_list]
        facts_list = facts_dict[crd]
        # Удаляем из словаря фактов для того, чтобы не попадал туда заново
        precision_score = precision_crd(recs_list, facts_list, weight_by_price, recs_prices_list)
        prec_list.append(precision_score)

    return np.mean(prec_list)

max_stars_repo_name: geangohn/RecSys

max_issues_repo_name: geangohn/RecSys

max_forks_repo_name: geangohn/RecSys

reason: import numpy,import scipy


text:
import numpy as np
import argparse
import glob
import os
from functools import partial
import vispy
import scipy.misc as misc
from tqdm import tqdm
import yaml
import json
import time
import sys
import shutil
from mesh import write_ply, read_ply, output_3d_photo
from utils import get_MiDaS_samples, read_MiDaS_depth
import torch
import cv2
from skimage.transform import resize
import imageio
import copy
from networks import Inpaint_Color_Net, Inpaint_Depth_Net, Inpaint_Edge_Net
from MiDaS.run import run_depth
from MiDaS.monodepth_net import MonoDepthNet
import MiDaS.MiDaS_utils as MiDaS_utils
from bilateral_filtering import sparse_bilateral_filtering
import boto3

parser = argparse.ArgumentParser()
parser.add_argument('--config', type=str, default='argument.yml',help='Configure of post processing')
args = parser.parse_args()
config = yaml.load(open(args.config, 'r'))
if config['offscreen_rendering'] is True:
    vispy.app.use_app('osmesa')
os.makedirs(config['mesh_folder'], exist_ok=True)
os.makedirs(config['video_folder'], exist_ok=True)
os.makedirs(config['depth_folder'], exist_ok=True)
if isinstance(config["gpu_ids"], int) and (config["gpu_ids"] >= 0):
    device = config["gpu_ids"]
else:
    device = "cpu"
print(f"running on device {device}")

sqs = boto3.client("sqs")
s3 = boto3.resource("s3")
s3_uploader = boto3.client("s3")
url = os.environ["AMENA_INPUT_QUEUE_URL"]
input_bucket_name = os.environ["AMENA_INPUT_BUCKET_NAME"]
output_bucket_name = os.environ["AMENA_OUTPUT_BUCKET_NAME"]
s3_input_image_path = "{account_id}/{input_id}/input.jpg"
s3_output_image_path = "{account_id}/{input_id}/output.mp4"
s3_output_error_path = "{account_id}/{input_id}/error.json"
local_input_image_path = config['src_folder'] + "/{account_id}_{input_id}.jpg"
dequeue_wait_time = 20
input_bucket = s3.Bucket(input_bucket_name)

def dequeueInputImage():
    received_queues = []

    while True:
        received_queues = sqs.receive_message(
            QueueUrl=url,
            MaxNumberOfMessages=1,
            VisibilityTimeout=5,
            WaitTimeSeconds=dequeue_wait_time
        )
        if "Messages" in received_queues:
            break

    received_queue = received_queues["Messages"][0]

    queue_id = received_queue["MessageId"]
    receipt_handle = received_queue["ReceiptHandle"]
    queue_body = json.loads(received_queue["Body"])
    queue_type = queue_body["type"]
    account_id = queue_body["account_id"]

    s3_object_key = s3_input_image_path.replace("{account_id}", account_id).replace("{input_id}", queue_id)
    s3_object = input_bucket.Object(s3_object_key)
    img_byte = s3_object.get().get('Body').read()
    img = cv2.imdecode(np.asarray(bytearray(img_byte)), cv2.IMREAD_COLOR)

    return {
        "img": img,
        "account_id": account_id,
        "input_id": queue_id,
        "s3_object": s3_object,
        "sqs_receipt_handle": receipt_handle
    }

while True:
    # Clear work folders.
    if os.path.exists(config['src_folder']):
        shutil.rmtree(config['src_folder'])
    if os.path.exists(config['video_folder']):
        shutil.rmtree(config['video_folder'])
    if os.path.exists(config['depth_folder']):
        shutil.rmtree(config['depth_folder'])
    if os.path.exists(config['mesh_folder']):
        shutil.rmtree(config['mesh_folder'])

    os.mkdir(config['src_folder'])
    os.mkdir(config['video_folder'])
    os.mkdir(config['depth_folder'])
    os.mkdir(config['mesh_folder'])

    # Dequeue request from sqs
    # if failed, continue next
    try:
        dequeued = dequeueInputImage()
    except Exception as e:
        continue

    sqs.change_message_visibility(
        QueueUrl=url,
        ReceiptHandle=dequeued["sqs_receipt_handle"],
        VisibilityTimeout=900
    )

    src_path = local_input_image_path.replace("{account_id}", dequeued["account_id"]).replace("{input_id}", dequeued["input_id"])
    cv2.imwrite(src_path, dequeued["img"])
    print(f"Dequeued: {src_path}")

    s3_output_key = s3_output_image_path.replace("{account_id}", dequeued["account_id"]).replace("{input_id}", dequeued["input_id"])
    s3_error_key = s3_output_error_path.replace("{account_id}", dequeued["account_id"]).replace("{input_id}", dequeued["input_id"])

    try:
        # Get dequeued local input images.
        sample_list = get_MiDaS_samples(config['src_folder'], config['depth_folder'], config, config['specific'])
        normal_canvas, all_canvas = None, None

        depth = None
        sample = sample_list[0]
        print("Current Source ==> ", sample['src_pair_name'])
        mesh_fi = os.path.join(config['mesh_folder'], sample['src_pair_name'] +'.ply')
        image = imageio.imread(sample['ref_img_fi'])

        print(f"Running depth extraction at {time.time()}")
        if config['require_midas'] is True:
            run_depth([sample['ref_img_fi']], config['src_folder'], config['depth_folder'],
                    config['MiDaS_model_ckpt'], MonoDepthNet, MiDaS_utils, target_w=640)
        if 'npy' in config['depth_format']:
            config['output_h'], config['output_w'] = np.load(sample['depth_fi']).shape[:2]
        else:
            config['output_h'], config['output_w'] = imageio.imread(sample['depth_fi']).shape[:2]
        frac = config['longer_side_len'] / max(config['output_h'], config['output_w'])
        config['output_h'], config['output_w'] = int(config['output_h'] * frac), int(config['output_w'] * frac)
        config['original_h'], config['original_w'] = config['output_h'], config['output_w']
        if image.ndim == 2:
            image = image[..., None].repeat(3, -1)
        if np.sum(np.abs(image[..., 0] - image[..., 1])) == 0 and np.sum(np.abs(image[..., 1] - image[..., 2])) == 0:
            config['gray_image'] = True
        else:
            config['gray_image'] = False
        image = cv2.resize(image, (config['output_w'], config['output_h']), interpolation=cv2.INTER_AREA)
        depth = read_MiDaS_depth(sample['depth_fi'], 3.0, config['output_h'], config['output_w'])
        mean_loc_depth = depth[depth.shape[0]//2, depth.shape[1]//2]
        if not(config['load_ply'] is True and os.path.exists(mesh_fi)):
            vis_photos, vis_depths = sparse_bilateral_filtering(depth.copy(), image.copy(), config, num_iter=config['sparse_iter'], spdb=False)
            depth = vis_depths[-1]
            model = None
            torch.cuda.empty_cache()
            print("Start Running 3D_Photo ...")
            print(f"Loading edge model at {time.time()}")
            depth_edge_model = Inpaint_Edge_Net(init_weights=True)
            depth_edge_weight = torch.load(config['depth_edge_model_ckpt'],
                                        map_location=torch.device(device))
            depth_edge_model.load_state_dict(depth_edge_weight)
            depth_edge_model = depth_edge_model.to(device)
            depth_edge_model.eval()

            print(f"Loading depth model at {time.time()}")
            depth_feat_model = Inpaint_Depth_Net()
            depth_feat_weight = torch.load(config['depth_feat_model_ckpt'],
                                        map_location=torch.device(device))
            depth_feat_model.load_state_dict(depth_feat_weight, strict=True)
            depth_feat_model = depth_feat_model.to(device)
            depth_feat_model.eval()
            depth_feat_model = depth_feat_model.to(device)
            print(f"Loading rgb model at {time.time()}")
            rgb_model = Inpaint_Color_Net()
            rgb_feat_weight = torch.load(config['rgb_feat_model_ckpt'],
                                        map_location=torch.device(device))
            rgb_model.load_state_dict(rgb_feat_weight)
            rgb_model.eval()
            rgb_model = rgb_model.to(device)
            graph = None

            print(f"Writing depth ply (and basically doing everything) at {time.time()}")
            rt_info = write_ply(image,
                                depth,
                                sample['int_mtx'],
                                mesh_fi,
                                config,
                                rgb_model,
                                depth_edge_model,
                                depth_edge_model,
                                depth_feat_model)

            if rt_info is False:
                continue
            rgb_model = None
            color_feat_model = None
            depth_edge_model = None
            depth_feat_model = None
            torch.cuda.empty_cache()
        if config['save_ply'] is True or config['load_ply'] is True:
            verts, colors, faces, Height, Width, hFov, vFov = read_ply(mesh_fi)
        else:
            verts, colors, faces, Height, Width, hFov, vFov = rt_info

        print(f"Making video at {time.time()}")
        videos_poses, video_basename = copy.deepcopy(sample['tgts_poses']), sample['tgt_name']
        top = (config.get('original_h') // 2 - sample['int_mtx'][1, 2] * config['output_h'])
        left = (config.get('original_w') // 2 - sample['int_mtx'][0, 2] * config['output_w'])
        down, right = top + config['output_h'], left + config['output_w']
        border = [int(xx) for xx in [top, down, left, right]]
        normal_canvas, all_canvas, video_path = output_3d_photo(verts.copy(), colors.copy(), faces.copy(), copy.deepcopy(Height), copy.deepcopy(Width), copy.deepcopy(hFov), copy.deepcopy(vFov),
                            copy.deepcopy(sample['tgt_pose']), sample['video_postfix'], copy.deepcopy(sample['ref_pose']), copy.deepcopy(config['video_folder']),
                            image.copy(), copy.deepcopy(sample['int_mtx']), config, image,
                            videos_poses, video_basename, config.get('original_h'), config.get('original_w'), border=border, depth=depth, normal_canvas=normal_canvas, all_canvas=all_canvas,
                            mean_loc_depth=mean_loc_depth)

        # Put artifacts to s3
        s3_uploader.upload_file(video_path, output_bucket_name, s3_output_key)
        print(f"Put artifacts: {s3_output_key}")

    except Exception as e:
        print("Error: " + repr(e))

        error = { "message": "Unknown error." }

        s3_uploader.put_object(
            Bucket=output_bucket_name,
            Key=s3_error_key,
            Body=json.dumps(error)
        )

        print(f"Put error log: {s3_error_key}")

    finally:
        # delete input source
        sqs.delete_message(
            QueueUrl=url,
            ReceiptHandle=dequeued["sqs_receipt_handle"]
        )

        dequeued["s3_object"].delete()

ext: py

max_stars_repo_name: uehr/amena

max_issues_repo_name: uehr/amena

max_forks_repo_name: uehr/amena


text:
'''
Function to generate capital

By Dr. Raymond Hoogendoorn
Copyright 2020
'''

from numpy.random import choice
from random import randint

def generateCapital(prob_distribution, age):
    elements = range(0,12)
    min_range = [-50000, -5000, 0, 1000, 5000, 10000, 20000, 50000, 100000, 200000, 500000, 1000000]
    max_range = [-5000, 0, 1000, 5000, 10000, 20000, 50000, 100000, 200000, 500000, 1000000, 5000000]
    if age >= 18:
        capitalclass = choice(elements, p = prob_distribution)
        min_capital = min_range[capitalclass]
        max_capital = max_range[capitalclass]
        capital = randint(min_capital, max_capital)
        return capital
    else:
        return 0

max_stars_repo_name: virtualsociety/simulation-py

max_issues_repo_name: virtualsociety/simulation-py

max_forks_repo_name: virtualsociety/simulation-py

reason: from numpy


text:
import os
import chainer
from chainer.training import extension
from chainer import Variable, cuda
import chainer.functions as F
import numpy as np
import os
import cv2
from .save_images import save_images_grid

def gan_sampling(gen, eval_folder, gpu, rows=6, cols=6, latent_len=128):
    @chainer.training.make_extension()
    def samples_generation(trainer):
        if not os.path.exists(eval_folder):
            os.makedirs(eval_folder)
        z = np.random.normal(size=(rows*cols, latent_len)).astype("f")
        if gpu>=0:
            z = cuda.to_gpu(z)
        z = Variable(z, volatile=True)
        imgs = gen(z, test=True)
        save_images_grid(imgs, path=eval_folder+"/iter_"+str(trainer.updater.iteration)+".jpg",
            grid_w=rows, grid_h=cols)

    return samples_generation

def gan_sampling_tags(gen, eval_folder, gpu, rows=6, cols=6, latent_len=128, attr_len=38, threshold=0.25):
    @chainer.training.make_extension()
    def get_fake_tag():
        prob2 = np.random.rand(attr_len)
        tags = np.zeros((attr_len)).astype("f")
        tags[:] = -1.0
        tags[np.argmax(prob2[0:13])]=1.0
        tags[27 + np.argmax(prob2[27:])] = 1.0
        prob2[prob2<threshold] = -1.0
        prob2[prob2>=threshold] = 1.0
        for i in range(13, 27):
            tags[i] = prob2[i]
        return tags

    def get_fake_tag_batch():
        xp = gen.xp
        batch = rows*cols
        tags = xp.zeros((batch, attr_len)).astype("f")
        for i in range(batch):
            tags[i] = xp.asarray(get_fake_tag())
        return tags

    def samples_generation(trainer):
        if not os.path.exists(eval_folder):
            os.makedirs(eval_folder)
        z = np.random.normal(size=(rows*cols, latent_len)).astype("f")
        if gpu>=0:
            z = cuda.to_gpu(z)
        tags =get_fake_tag_batch()
        z = Variable(z, volatile=True)
        tags = Variable(tags, volatile=True)
        imgs = gen(F.concat([z,tags]), test=True)
        save_images_grid(imgs, path=eval_folder+"/iter_"+str(trainer.updater.iteration)+".jpg",
            grid_w=rows, grid_h=cols)

    return samples_generation

def ae_reconstruction(enc, dec, eval_folder, gpu, data_iter, batch_size=32, img_chan=3, img_size=64):
    @chainer.training.make_extension()
    def sample_reconstruction(trainer):
        xp = enc.xp
        batch = data_iter.next()
        d_real = xp.zeros((batch_size, img_chan, img_size, img_size)).astype("f")
        for i in range(batch_size):
            d_real[i, :] = xp.asarray(batch[i])
        x = Variable(d_real, volatile=True)
        imgs = dec(enc(x, test=True), test=True)
        save_images_grid(imgs, path=eval_folder+"/iter_"+str(trainer.updater.iteration)+".rec.jpg",
            grid_w=batch_size//8, grid_h=8)
        save_images_grid(d_real, path=eval_folder+"/iter_"+str(trainer.updater.iteration)+".real.jpg",
            grid_w=batch_size//8, grid_h=8)

    return sample_reconstruction

def analogy(gen, output, samples=12, latent_len=128, points=10):
    xp = gen.xp
    z0 = xp.random.normal(size=(samples, latent_len)).astype("f")
    z1 = xp.random.normal(size=(samples, latent_len)).astype("f")
    values = np.linspace(0.0, 1.0, num=points)
    results = []
    for i in range(points):
        z = (values[i]*z0 + (1.0-values[i])*z1).astype("f")
        z = Variable(z, volatile=True)
        imgs = gen(z, test=True)
        results.append(imgs.data)
    results = xp.concatenate(results, axis=0)
    save_images_grid(results, path=output, grid_w=points, grid_h=samples, transposed=True)

max_stars_repo_name: Aixile/chainer-gan-experiments

max_issues_repo_name: Aixile/chainer-gan-experiments

max_forks_repo_name: Aixile/chainer-gan-experiments

reason: import numpy


text:
import numpy as np
import matplotlib.pyplot as plt

test = np.load("/Users/robin/PycharmProjects/DeepCCA/examples/new_features.npy", allow_pickle=True)

test_data, test_pred, test_label = test[2]

plt.plot(test_data[:,0], test_pred[:,0], 'ro')
plt.show()

max_stars_repo_name: robinthibaut/deep-cca

max_issues_repo_name: robinthibaut/deep-cca

max_forks_repo_name: robinthibaut/deep-cca


text:
import pandas as pd
import numpy as np
import time
import sklearn
import pickle

from collections import namedtuple
from sklearn.linear_model import LinearRegression, Ridge
from sklearn import ensemble
from sklearn.neural_network import MLPRegressor
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.pipeline import make_pipeline

import xgboost as xgb

# This is for querying the estimated area of a given spatial node

class ReportStats(namedtuple('ReportStats', ['mse', 'mae', 'r2', 'abs_error', 'rel_error', \
                                             'failed_node_names', 'failed_node_gold', 'failed_node_pred'
                                             ])):
    pass

class IRTrainingSession(namedtuple('IRTrainingSession', ['node_type', 'input_feature_list'])):
    pass

class ModelMetaData(namedtuple('ModelMetaData', ['best_model_name', 'best_mae', 'model_pickle', 'prediction_task'])):
    pass

XilinxFPGAResourceTargets = [
    'LUTs', 'FFs', 'SRLs', 'RAMB32', 'RAMB18', 'DSPs'
]

type = "par"
input_file = "regression_full." + type + ".csv"
df = pd.read_csv(input_file, header=0)
trained_model_dir = "./trained_models/"
df.sort_values(by='nodetype', axis=0, inplace=True)
names = df['nodetype'].unique().tolist()
subframes = [df.loc[df.nodetype == name].dropna(axis=1, how='all') for name in names]

ensemble_params = {
    "n_estimators": 500,
    "max_depth": 4,
    "min_samples_split": 2,
    "learning_rate": 0.01,
    "loss": "ls"
}
model_zoo = [
    #     (
    #         'PolynomialFeatures(1)+xgb.XGBRegressor',
    #         make_pipeline(PolynomialFeatures(1), xgb.XGBRegressor(objective="reg:linear", random_state=42))
    #     ),
    #     (
    #         'PolynomialFeatures(2)+xgb.XGBRegressor',
    #         make_pipeline(PolynomialFeatures(2), xgb.XGBRegressor(objective="reg:linear", random_state=42))
    #     ),
    #     (
    #         'LinearRegression',
    #         make_pipeline(LinearRegression())
    #     ),

    #     (
    #         'PolynomialFeatures(2)+Ridge',
    #         make_pipeline(PolynomialFeatures(2), Ridge())
    #     ),
    #     (
    #         'StandardScaler+MLPRegressor',
    #         make_pipeline(StandardScaler(), MLPRegressor(hidden_layer_sizes=(128,), max_iter=10000))
    #     ),
    #     (
    #         'PolynomialFeatures(2)+MLPRegressor',
    #         make_pipeline(PolynomialFeatures(2), MLPRegressor(hidden_layer_sizes=(256,), max_iter=10000))
    #     ),
    #     (
    #         'MLPRegressor',
    #         make_pipeline(MLPRegressor(hidden_layer_sizes=(128,), max_iter=10000))
    #     ),
    #     (
    #         'MLPRegressor(100, 100)',
    #         make_pipeline(MLPRegressor(hidden_layer_sizes=(100, 100), max_iter=100000))
    #     ),
    (
        'RandomForestRegressor',
        make_pipeline(ensemble.RandomForestRegressor())
    ),
    (
        'ensemble_GradientBoostingRegressor',
        make_pipeline(ensemble.GradientBoostingRegressor(**ensemble_params))
    ),
    (
        'XGBRegressor',
        make_pipeline(xgb.XGBRegressor(objective="reg:linear", random_state=42))
    ),
    (
        'PolynomialFeatures-2-RandomForestRegressor',
        make_pipeline(PolynomialFeatures(2), ensemble.RandomForestRegressor())
    ),
    (
        'PolynomialFeatures-2-ensemble_GradientBoostingRegressor',
        make_pipeline(PolynomialFeatures(2), ensemble.GradientBoostingRegressor(**ensemble_params))
    ),

    (
        'PolynomialFeatures-2-XGBRegressor',
        make_pipeline(PolynomialFeatures(2), xgb.XGBRegressor(objective="reg:linear", random_state=42))
    )
]

# clean up the dataframe where LUTs = 0
def predict(model, x):
    start = time.time()
    result = model.predict(x).astype(int).clip(min=0)
    end = time.time()
    # print("prediction takes " + str(end - start))
    return result

def get_shuffled_data(data_frame, X_column_names, Y_column_names, split_ratio=0.80):
    X = data_frame.loc[:, X_column_names]
    X = X.fillna(0)
    Y = data_frame.loc[:, Y_column_names]
    n_samples, n_features = X.shape
    _, n_targets = Y.shape
    n_split_index = int(split_ratio * n_samples)

    XX = np.array(X).reshape(n_samples, n_features)
    Y = np.array(Y).reshape(n_samples, n_targets)
    XY = np.concatenate([X, Y], axis=1)
    # clean up invalid entries
    #     invalid_rows = np.where(Y == 0)
    #     XY = np.delete(XY, invalid_rows, axis=0)

    np.random.shuffle(XY)

    train_fullname = XY[:n_split_index, 0]
    test_fullname = XY[n_split_index:, 0]
    train_set = XY[:n_split_index, 1:]
    train_set_X = train_set[:, :-1]
    train_set_Y = train_set[:, -1]
    test_set = XY[n_split_index:, 1:]
    test_set_X = test_set[:, :-1]
    test_set_Y = test_set[:, -1]

    return (train_fullname, test_fullname, train_set_X, train_set_Y, test_set_X, test_set_Y)

def get_metrics(gold, pred, names, fail_thresh=0.2):
    abs_error = np.abs(gold - pred).flatten()
    rel_error = abs_error  # if (gold == 0) else abs_error / gold; cannot think of a good way to do this...
    failed_idx = np.where(rel_error > fail_thresh)
    failed_node_names = names[failed_idx]
    failed_node_gold = gold[failed_idx]
    failed_node_pred = pred[failed_idx]

    return ReportStats(
        mse=mean_squared_error(gold, pred),
        mae=mean_absolute_error(gold, pred),
        r2=r2_score(gold, pred),
        abs_error=np.amax(abs_error),
        rel_error=np.amax(rel_error),
        failed_node_names=failed_node_names,
        failed_node_gold=failed_node_gold,
        failed_node_pred=failed_node_pred
    )

def eval_model(model_template, x_train, y_train, x_test, y_test, train_names, test_names):
    model = sklearn.base.clone(model_template)
    model.fit(x_train, y_train)
    y_test_pred = predict(model, x_test)
    y_train_pred = predict(model, x_train)

    train_stats = get_metrics(y_train, y_train_pred, train_names)
    test_stats = get_metrics(y_test, y_test_pred, test_names)
    return (train_stats, test_stats, model)

def test_and_store_pickled_model(df, sess, predict_type, model_name, pickled_model):
    test_model = pickle.loads(pickled_model)
    file_name = sess.node_type + "_" + predict_type + ".pkl"
    pickle.dump(test_model, open(file_name, "wb"))
    ff = df.loc[df.nodetype == sess.node_type]

    Y_column_names = [predict_type]
    X_column_names = sess.input_feature_list

    _, _, \
    x_test_0, y_test_0, \
    x_test_1, y_test_1 = \
        get_shuffled_data(
            ff, X_column_names, Y_column_names, split_ratio=0.80
        )

    # This is a rough check to estimate how well the model works on the full data set...
    # Cannot be taken as a valid statistical measure...
    # Would like to see what's the largest deviation...
    x_test = np.concatenate([x_test_0, x_test_1], axis=0)
    y_test = np.concatenate([y_test_0, y_test_1], axis=0)

    y_pred = test_model.predict(x_test)

    print("Testing the pickled model on " + predict_type + ". mae = ", mean_absolute_error(y_pred, y_test))

def retrain_models_on_session(df, sess):
    """
    Retrain a session specification with all the models in the model zoo.
    Report the one that achieves the lowest MAE.
    """
    ff = df.loc[df.nodetype == sess.node_type]

    def get_single_sess_regression(predict_type):
        Y_column_names = [predict_type]
        X_column_names = sess.input_feature_list
        train_full_name, test_full_name, \
        x_train, y_train, \
        x_test, y_test = \
            get_shuffled_data(
                ff, X_column_names, Y_column_names, split_ratio=0.80
            )
        print("Training Set contains ", x_train.shape[0], " points")
        print("Test Set contains ", x_test.shape[0], " points")

        stats = [
            (
                name,
                eval_model(
                    f, x_train, y_train, \
                    x_test, y_test, \
                    train_full_name, test_full_name
                )
            ) for (name, f) in model_zoo
        ]
        model_names = np.array([ts[0] for ts in stats])
        maes = np.array([ts[1][1].mae for ts in stats])
        model_copies = np.array([ts[1][2] for ts in stats])

        best_mae = np.amin(maes)
        best_mae_idx = np.argmin(maes)
        best_model_name = model_names[best_mae_idx]
        best_model_pickle = pickle.dumps(model_copies[best_mae_idx])

        selected_model = ModelMetaData(best_model_name, best_mae, best_model_pickle, predict_type)

        return selected_model

    best_models = [get_single_sess_regression(predict_type) for predict_type in XilinxFPGAResourceTargets]
    return best_models

FIFOSess = IRTrainingSession(
    "FIFONew",
    [
        "fullname",
        "mem_B0",  # not related?
        "mem_N0",
        "mem_a0",
        "mem_bitwidth",
        "mem_dim0",
        "mem_hist0muxwidth",
        "mem_hist0rlanes",
        "mem_hist0wlanes",
        "mem_nbufs",  # not related?
        "mem_p0"
    ]
)
# [1,1,1,1,1,8,1,1,1,1,1,1,1,1,1,32,2048,1,1,1,1,1,8,8,8,2,1,1,1,1,1,8,1,1,1,1]
SRAMSess = IRTrainingSession(
    "SRAMNew",
    [
        "fullname",
        "mem_B0",
        "mem_B1",
        "mem_B2",
        "mem_B3",
        "mem_B4",
        "mem_N0",
        "mem_N1",
        "mem_N2",
        "mem_N3",
        "mem_N4",
        "mem_a0",
        "mem_a1",
        "mem_a2",
        "mem_a3",
        "mem_a4",
        "mem_bitwidth",
        "mem_dim0",
        "mem_dim1",
        "mem_dim2",
        "mem_dim3",
        "mem_dim4",
        "mem_hist0muxwidth",
        "mem_hist0rlanes",
        "mem_hist0wlanes",
        "mem_hist1muxwidth",
        "mem_hist1rlanes",
        "mem_hist1wlanes",
        "mem_hist2muxwidth",
        "mem_hist2rlanes",
        "mem_hist2wlanes",
        "mem_nbufs",
        "mem_p0",
        "mem_p1",
        "mem_p2",
        "mem_p3",
        "mem_p4"
    ]
)

CounterChainSess = IRTrainingSession(
    "CounterChainNew",
    [
        "fullname",
        "cchain_ctr0par",
        "cchain_ctr0start",
        "cchain_ctr0step",
        "cchain_ctr0stop",
        "cchain_ctr1par",
        "cchain_ctr1start",
        "cchain_ctr1step",
        "cchain_ctr1stop",
        "cchain_ctr2par",
        "cchain_ctr2start",
        "cchain_ctr2step"
    ]
)

FixAddSess = IRTrainingSession(
    "FixAdd",
    [
        "fullname",
        "op_consta",
        "op_constb",
        "tp_dec",
        "tp_frac",
        "tp_sgn"
    ]
)

FixDivSess = IRTrainingSession(
    "FixDiv",
    [
        "fullname",
        "op_consta",
        "op_constb",
        "tp_dec",
        "tp_frac",
        "tp_sgn"
    ]
)

FixFMASess = IRTrainingSession(
    "FixFMA",
    [
        "fullname",
        "tp_dec",
        "tp_frac",
        "tp_sign"
    ]
)

FixMulSess = IRTrainingSession(
    "FixMul",
    [
        "op_consta",
        "op_constb",
        "tp_dec",
        "tp_frac",
        "tp_sign"
    ]
)

FixSubSess = IRTrainingSession(
    "FixSub",
    [
        "op_consta",
        "op_constb",
        "tp_dec",
        "tp_frac",
        "tp_sign"
    ]
)

FixToFixSess = IRTrainingSession(
    "FixToFix",
    [
        "tp_dec",
        "tp_frac",
        "tp_sgn",
        "tp_srcdec",
        "tp_srcfrac",
        "tp_srcsgn"
    ]
)

LIFONewSess = IRTrainingSession(
    "LIFONew",
    [
        "mem_B0",
        "mem_N0",
        "mem_a0",
        "mem_bitwidth",
        "mem_dim0",
        "mem_hist0muxwidth",
        "mem_hist0rlanes",
        "mem_hist0wlanes",
        "mem_nbufs",
        "mem_p0"

    ]
)

LineBufferNewSess = IRTrainingSession(
    "LineBufferNew",
    [
        "mem_B0",
        "mem_B1",
        "mem_N0",
        "mem_N1",
        "mem_a0",
        "mem_a1",
        "mem_bitwidth",
        "mem_dim0",
        "mem_dim1",
        "mem_hist0muxwidth",
        "mem_hist0rlanes",
        "mem_hist0wlanes",
        "mem_hist1muxwidth",
        "mem_hist1rlanes",
        "mem_hist1wlanes",
        "mem_nbufs",
        "mem_p0",
        "mem_p1"
    ]
)

ParallelPipeSess = IRTrainingSession(
    "ParallelPipe",
    [
        "ctrl_children",
        "ctrl_ii",
        "ctrl_lat",
        "ctrl_level",
        "ctrl_style"
    ]
)

RegFileNewSess = IRTrainingSession(
    "RegFileNew",
    [
        "mem_B0",
        "mem_B1",
        "mem_N0",
        "mem_N1",
        "mem_a0",
        "mem_a1",
        "mem_bitwidth",
        "mem_dim0",
        "mem_dim1",
        "mem_hist0muxwidth",
        "mem_hist0rlanes",
        "mem_hist0wlanes",
        "mem_hist1muxwidth",
        "mem_hist1rlanes",
        "mem_hist1wlanes",
        "mem_nbufs",
        "mem_p0",
        "mem_p1"
    ]
)

RegNewSess = IRTrainingSession(
    "RegNew",
    [
        "mem_bitwidth",
        "mem_nbufs"
    ]
)

StateMachineSess = IRTrainingSession(
    "StateMachine",
    [
        "ctrl_children",
        "ctrl_ii",
        "ctrl_lat",
        "ctrl_level",
        "ctrl_style"
    ]
)

UnbMulSess = IRTrainingSession(
    "UnbMul",
    [
        "op_consta",
        "tp_dec",
        "tp_frac",
        "tp_sgn"
    ]
)

UnbSatMulSess = IRTrainingSession(
    "UnbSatMul",
    [
        "tp_dec",
        "tp_frac",
        "tp_sgn"
    ]
)

UnitPipeSess = IRTrainingSession(
    "UnitPipe",
    [
        "ctrl_children",
        "ctrl_ii",
        "ctrl_lat",
        "ctrl_level",
        "ctrl_style"
    ]
)

UnrolledForeachSess = IRTrainingSession(
    "UnrolledForeach",
    [
        "ctrl_children",
        "ctrl_ii",
        "ctrl_lat",
        "ctrl_level",
        "ctrl_style"

    ]
)

UnrolledReduceSess = IRTrainingSession(
    "UnrolledReduce",
    [
        "ctrl_children",
        "ctrl_ii",
        "ctrl_lat",
        "ctrl_level",
        "ctrl_style"

    ]
)

sessions = [
    FIFOSess,
    SRAMSess,
    CounterChainSess,
    FixAddSess,
    FixDivSess,
    FixFMASess,
    FixMulSess,
    FixSubSess,
    FixToFixSess,
    LIFONewSess,
    LineBufferNewSess,
    ParallelPipeSess,
    RegFileNewSess,
    RegNewSess,
    StateMachineSess,
    UnbMulSess,
    UnbSatMulSess,
    UnitPipeSess,
    UnrolledForeachSess,
    UnrolledReduceSess
]

with open("test_mae.csv", "w+") as f:
    f.write(" ,")
    for rsc in XilinxFPGAResourceTargets:
        f.write(rsc + ", ")
    f.write("\n")

    for sess in sessions:
        f.write(sess.node_type + ", ")
        print(sess.node_type)
        model_metas = retrain_models_on_session(df, sess)
        for mm in model_metas:
            f.write(str(mm.best_mae) + ", ")
            # print(mm.best_mae, mm.best_model_name, mm.prediction_task)
        f.write("\n")

        for mm in model_metas:
            test_and_store_pickled_model(df, sess, mm.prediction_task, mm.best_model_name, mm.model_pickle)

max_stars_repo_name: Hecmay/spatial

max_issues_repo_name: Hecmay/spatial

max_forks_repo_name: Hecmay/spatial


text:
from mayavi import mlab
import h5py
import numpy as np
import vigra
import time

import trainingcore
import visCell
import hdf5io

from traits.api import *
from traitsui.api import *

from mayavi.tools.mlab_scene_model import MlabSceneModel
from mayavi.core.ui.mayavi_scene import MayaviScene
from tvtk.pyface.scene_editor import SceneEditor

class msgdialog(HasTraits):
    message = Label("")

    def __init__(self, message):
        self.message = message

    view = View( Item('message',show_label=False,style='readonly'),
                 buttons = ['OK'],
                 title='Message')

def showmessage(message):
    msg = msgdialog(message)
    msg.configure_traits()

class traininggui(HasTraits):

    """ Group: File selection """
    fileIndex = Int(0)
    btnLoadFile = Button()
    btnLoadNextFile = Button()
    btnLoadPreviousFile = Button()
    btnNextHere = Button()
    btnPrevHere = Button()

    """ Info """
    numFiles = Int(0)
    numObjects = Int(0)
    numLabels = Int(0)
    numClasses = Int(0)
    classCount = Array(dtype=np.int32)

    """ Group: Jump navigation """
    btnJump = Button()
    btnJumpRandom = Button()
    index = Int(1)

    """ Group: Forward navigation """
    btnShowNext = Button()
    btnShowPrevious = Button()
    btnShowNextLabel = Button()
    btnShowPreviousLabel = Button()
    nextLabel = Int(0)
    btnShowNextUncertain = Button()

    """ Group: Labeling"""
    label = Int(-1)
    btnZero = Button()
    btnOne = Button()
    btnLblAs = Button()
    labelAs = Int(1)

    """ Group: Settings """
    borderSize = Int(30)
    datasetTraining = Str('training')
    datasetPrediction = Str('prediction')
    useRandomForest = Bool(False)
    randomForestTrees = Int(255)

    """ Group: Train&Predict """
    btnTrainIntern = Button()
    extLabelSource = Str('training')
    btnTrainExtern = Button()

    predDestination = Str('prediction')
    predType = Enum('Class', 'Class Probabilities')
    btnPredict = Button()

    btnSaveRF = Button()
    btnLoadRF = Button()
    RFFile = Str('RandomForest.h5')

    """ Scene views """
    scene3d = Instance(MlabSceneModel, ())
    scene2d = Instance(MlabSceneModel, ())

    def __init__(self, files):
        self.files = files
        self.fileIndex = 0
        self.t = trainingcore.trainingcore(self.files[self.fileIndex], scene3d=self.scene3d, scene2d=self.scene2d)
        self.t.jump_to_cell(self.index)
        self.prepareRF()
        self.update_info()
        self.seriousRF = None

    def _btnLoadFile_fired(self):
        if self.fileIndex >= 0 and self.fileIndex < len(self.files):
            self.t = trainingcore.trainingcore(self.files[self.fileIndex], self.index)
            self.t.jump_to_cell(self.index)
            self.update_info()
            self.t.borderSize = self.borderSize
        else:
            showmessage("Error: Cannot view this time step - index out of bounds.")

    def _btnLoadNextFile_fired(self):
        if self.fileIndex+1 < len(self.files):
            self.fileIndex = self.fileIndex +1
            self.t = trainingcore.trainingcore(self.files[self.fileIndex], self.index)
            self.t.jump_to_cell(self.index)
            self.update_info()
            self.t.borderSize = self.borderSize
        else:
            showmessage("Error: Cannot load next time step - index out of bounds.")

    def _btnLoadPreviousFile_fired(self):
        if self.fileIndex > 0:
            self.fileIndex = self.fileIndex -1
            self.t = trainingcore.trainingcore(self.files[self.fileIndex], self.index)
            self.t.jump_to_cell(self.index)
            self.update_info()
            self.t.borderSize = self.borderSize
        else:
            showmessage("Error: Cannot load previous time step - index out of bounds.")

    def _btnNextHere_fired(self):
        if self.fileIndex+1 < len(self.files):
            currentBBox = self.t.currentBBox
            glanceFile = self.files[self.fileIndex+1]
            t = trainingcore.trainingcore(glanceFile)
            t.borderSize = self.borderSize
            t.show_volume(currentBBox)
            t.currentBBox = currentBBox
        else:
            showmessage("Error: Cannot view next time step - already reached last file.")

    def _btnPrevHere_fired(self):
        if self.fileIndex-1 >= 0:
            currentBBox = self.t.currentBBox
            glanceFile = self.files[self.fileIndex-1]
            t = trainingcore.trainingcore(glanceFile)
            t.borderSize = self.borderSize
            t.show_volume(currentBBox)
            t.currentBBox = currentBBox
        else:
            showmessage("Error: Cannot view previous time step - already reached first file.")

    def _btnShowNext_fired(self):
        if self.useRandomForest:
            self.t.RF = self.RF
        else:
            self.t.RF = None
        self.index = self.t.show_next_cell()
        self.label = self.find_label(self.index)

    def _btnShowPrevious_fired(self):
        if self.useRandomForest:
            self.t.RF = self.RF
        else:
            self.t.RF = None
        self.index = self.t.show_previous_cell()
        self.label = self.find_label(self.index)

    def _btnShowNextLabel_fired(self):
        if self.useRandomForest:
            self.t.RF = self.RF
        else:
            self.t.RF = None
        self.index = self.t.show_next_label_cell(self.datasetPrediction,self.nextLabel)
        self.label = self.find_label(self.index)

    def _btnShowPreviousLabel_fired(self):
        if self.useRandomForest:
            self.t.RF = self.RF
        else:
            self.t.RF = None
        self.index = self.t.show_previous_label_cell(self.datasetPrediction,self.nextLabel)
        self.label = self.find_label(self.index)

    def _btnShowNextUncertain_fired(self):
        if self.useRandomForest and self.RF != None:
            self.t.RF = self.RF
            self.index = self.t.show_next_uncertain_cell()
            self.label = self.find_label(self.index)
        else:
            self.t.RF = None
            showmessage("Error: Cannot find cells with high uncertainty. Need a pretrained Random Forest.")

    def _btnOne_fired(self):
        with h5py.File(self.files[self.fileIndex], 'a') as f:
            hdf5io.set_one_label(f,self.datasetTraining,self.index,1)
            newf = hdf5io.load_features_one_object(f,self.index)

        self.labels = np.append(self.labels,[[1]],axis=0).astype(np.int32)

        self.features = np.append(self.features, newf, axis=0).astype(np.float32)

        if self.useRandomForest and self.RF != None:
            self.RF.learnRF(self.features, self.labels.astype(np.uint32))
            self.t.RF = self.RF
        else:
            self.t.RF = None

        self.index = self.t.show_next_cell()
        self.label = self.find_label(self.index)
        self.update_info()

    def _btnZero_fired(self):
        with h5py.File(self.files[self.fileIndex], 'a') as f:
            hdf5io.set_one_label(f,self.datasetTraining,self.index,0)
            newf = hdf5io.load_features_one_object(f,self.index)

        self.labels = np.append(self.labels,[[0]],axis=0).astype(np.int32)

        self.features = np.append(self.features, newf, axis=0).astype(np.float32)

        if self.useRandomForest and self.RF != None:
            self.RF.learnRF(self.features, self.labels.astype(np.uint32))
            self.t.RF = self.RF
        else:
            self.t.RF = None

        self.index = self.t.show_next_cell()
        self.label = self.find_label(self.index)
        self.update_info()

    def _btnLblAs_fired(self):
        with h5py.File(self.files[self.fileIndex], 'a') as f:
            hdf5io.set_one_label(f,self.datasetTraining,self.index,self.labelAs)

        if self.labelAs >= 0:
            self.labels = np.append(self.labels,[[self.labelAs]],axis=0).astype(np.int32)
            with h5py.File(self.files[self.fileIndex], 'r') as f:
                newf = hdf5io.load_features_one_object(f,self.index)
            self.features = np.append(self.features, newf, axis=0).astype(np.float32)

        if self.useRandomForest and self.RF != None:
            self.RF.learnRF(self.features, self.labels.astype(np.uint32))
            self.t.RF = self.RF
        else:
            self.t.RF = None

        self.index = self.t.show_next_cell()
        self.label = self.find_label(self.index)
        self.update_info()

    def _btnJump_fired(self):
        if self.useRandomForest:
            self.t.RF = self.RF
        else:
            self.t.RF = None
        self.index = self.t.jump_to_cell(self.index)
        self.label = self.find_label(self.index)

    def _btnJumpRandom_fired(self):
        if self.useRandomForest:
            self.t.RF = self.RF
        else:
            self.t.RF = None
        self.index = self.t.jump_to_cell(np.random.randint(1,self.t.maxLabel+1))
        self.label = self.find_label(self.index)

    def _btnTrainIntern_fired(self):
        """
        """
        if np.unique(self.labels).shape[0] >= 2:
            self.seriousRF = vigra.learning.RandomForest(treeCount=self.randomForestTrees)
            self.seriousRF.learnRF(self.features, self.labels.astype(np.uint32))
            showmessage("Training successful. Used %i trees and %i training samples of %i classes."%(self.randomForestTrees,self.labels.shape[0],np.unique(self.labels).shape[0]))
        else:
            showmessage("Error: Not enough labels provided. Provide labels of at least two classes!")

    def _btnTrainExtern_fired(self):
        """
        """
        with h5py.File(self.files[self.fileIndex], 'r') as f:
            l = hdf5io.get_labels(f,self.extLabelSource).astype(np.int32)
        indices = l >= 0
        l = l[indices]

        if np.unique(l).shape[0] >= 2:
            with h5py.File(self.files[self.fileIndex], 'r') as fl:
                f = hdf5io.load_object_features(fl)[indices,:]
            self.seriousRF = vigra.learning.RandomForest(treeCount=self.randomForestTrees)
            self.seriousRF.learnRF(f, l.astype(np.uint32))
            showmessage("Training successful. Used %i trees and %i training samples of %i classes."%(self.randomForestTrees,l.shape[0],np.unique(l).shape[0]))
        else:
            showmessage("Error: Not enough labels found in dataset. Provide labels of at least two classes!")

    def _btnSaveRF_fired(self):
        """
        """
        if self.seriousRF != None and self.RFFile != "":
            self.seriousRF.writeHDF5(str(self.RFFile))
            showmessage('Successfully wrote Random Forest to file %s'%self.RFFile)
        else:
            showmessage('Error: train a Random Forest first and provide a filename.')

    def _btnLoadRF_fired(self):
        """
        """
        if self.RFFile != "":
            self.seriousRF = vigra.learning.RandomForest(str(self.RFFile))
            showmessage('Successfully loaded Random Forest from file %s'%self.RFFile)
        else:
            showmessage('Error: provide a filename.')

    def _btnPredict_fired(self):
        """
        """
        if self.seriousRF != None:
            with h5py.File(self.files[self.fileIndex], 'a') as fl:
                f = hdf5io.load_object_features(fl)
            if self.predType == 'Class':
                l = self.seriousRF.predictLabels(f)
                with h5py.File(self.files[self.fileIndex], 'a') as fl:
                    hdf5io.set_labels(fl,self.predDestination,l)
                showmessage('Successfully predicted labels of %i objects.'%f.shape[0])
            if self.predType == 'Class Probabilities':
                l = self.seriousRF.predictProbabilities(f)
                with h5py.File(self.files[self.fileIndex], 'a') as fl:
                    hdf5io.set_probabilities(fl,self.predDestination,l)
                showmessage('Successfully predicted probabilities of %i objects.'%f.shape[0])
        else:
            showmessage('Error: train the Random Forest first.')

    def find_label(self, index):
        """
        Returns the predicted label. If not available, it returns the
        trained label. If no label is found, -1 is returned.
        """
        with h5py.File(self.files[self.fileIndex], 'r') as f:
            l = hdf5io.get_labels(f,self.datasetPrediction)[index-1,0]
            if l == -1:
                l = hdf5io.get_labels(f,self.datasetTraining)[index-1,0]
        return l

    def update_info(self):
        """
        Update the values in the info field.
        """
        uniquelabels = np.unique(self.labels)
        count = np.zeros((uniquelabels.shape[0],2),dtype=np.int32)
        count[:,0] = uniquelabels[:]
        for i in range(uniquelabels.shape[0]):
            count[i,1] = np.sum(self.labels == uniquelabels[i])

        self.numFiles = len(self.files)
        self.numObjects = self.t.maxLabel
        self.numLabels = self.labels.shape[0]
        self.numClasses = len(uniquelabels)
        self.classCount = count

    def _borderSize_changed(self):
        """
        Apply a changed border size immediately
        """
        self.t.borderSize = self.borderSize

    def prepareRF(self):
        # prepare the random forest based on the features in h5file
        self.RF = vigra.learning.RandomForest(treeCount=80)
        with h5py.File(self.files[self.fileIndex], 'r') as f:
            nfeatures = hdf5io.load_features_one_object(f, 1).size
        self.features = np.zeros([0,nfeatures], dtype=np.float32)
        self.labels = np.zeros([0,1], dtype=np.int32)

    view = View(VSplit(
                    HSplit(
                        Item(name='scene3d',
                            editor=SceneEditor(scene_class=MayaviScene),
                            show_label=False,
                            height = 500
                        ),
                        Item(name='scene2d',
                            editor=SceneEditor(scene_class=MayaviScene),
                            show_label=False,
                            height = 500
                        ),
                    ),
                    HGroup(
                        Tabbed(
                            Group(
                                Item('fileIndex', label='Timestep' ),
                                '_',
                                HGroup(
                                    Item('btnLoadPreviousFile', label='<<<', show_label=False, tooltip="Load the previous time step." ),
                                    Item('btnLoadFile', label='Load This File', show_label=False, tooltip="Load/switch back to the selected time step." ),
                                    Item('btnLoadNextFile', label='>>>', show_label=False, tooltip="Loat the next time step." ),

                                ),
                                HGroup(
                                    Item('btnPrevHere', label='Glance previous', show_label=False, tooltip="Show the same position in the previous time step." ),
                                    Item('btnNextHere', label='Glance next', show_label=False, tooltip="Show the same position in the next time step." ),

                                ),
                                '_',
                                Label('- Info -'),
                                Item('numFiles', label='Number of Files', style='readonly'),
                                Item('numObjects', label='Number of Objects', style='readonly'),
                                Item('numLabels', label='Number of Labels', style='readonly'),
                                Item('numClasses', label='Number of Classes', style='readonly'),
                                Item('classCount', label='Labels per Class', style='readonly'),

                                label = 'File'
                            ),
                            Group(
                                Item('borderSize', label='Neighborhood Size', tooltip="Size of the neighborhood displayed around the object." ),
                                Item('datasetTraining', label='Write Labels' ),
                                Item('datasetPrediction', label='Read Labels' ),
                                Item('useRandomForest', label='Interactive Labeling', tooltip="Enable interactive labeling using a Random Forest classifier. This requires that labels for at least two different classes are provided."  ),
                                Item('randomForestTrees', label='RF Trees', tooltip="Number of trees used for the serious Random Forests." ),
                                label = "Settings"
                            ),
                        ),
                        Tabbed(
                            HGroup(
                                Group(
                                    HGroup(
                                        Item('btnShowPrevious', label='<<<', show_label=False, tooltip="Decrease object index by one."),
                                        Item('index', label='Selected index', show_label=False, tooltip="Object index."),
                                        Item('btnShowNext', label='>>>', show_label=False, tooltip="Increase object index by one."),
                                    ),
                                    HGroup(
                                        Item('btnJump', label='Jump to index', show_label=False, tooltip="Jump to object index given in the field above." ),
                                        Item('btnJumpRandom', label='Random index', show_label=False, tooltip="Jump to some ranodm object index." ),
                                    ),
                                    '_',
                                    Label('Find object with label:'),
                                    HGroup(
                                        Item('btnShowPreviousLabel',label='<<<',show_label=False, tooltip="Search for label in backward direction."),
                                        Item('nextLabel', label='Next label',show_label=False, tooltip="A label name to search for, e.g. '0'. Searches in the dataset specified in 'Read Labels'."),
                                        Item('btnShowNextLabel',label='>>>',show_label=False, tooltip="Search for label in forward direction."),
                                    ),
                                    '_',
                                    Label('Suggest:'),
                                    Item('btnShowNextUncertain',label='Suggest object',show_label=False, tooltip="Let the Random Forest decide which object to present next."),
                                ),
                                '_',
                                Group(
                                    HGroup(
                                        Label('Current label: '),
                                        Item('label', style='readonly', show_label=False, tooltip="Label assigned to selected object."),
                                    ),
                                    HGroup(
                                        Item('btnZero', label='Label as 0', show_label=False, tooltip="Assign label 0 to object." ),
                                        Item('btnOne', label='Label as 1', show_label=False, tooltip="Assign label 1 to object." ),
                                    ),
                                    HGroup(
                                        Item('btnLblAs', label='Label as ...', show_label=False, tooltip = "Assign a custom label to the object." ),
                                        Item('labelAs',label='Label',show_label=False, tooltip="Provide a custom label >= 0."),
                                    ),
                                ),
                                label = "Controls"
                            ),
                            Group(
                                Item('btnTrainIntern',label='Train with recorded labels',show_label=False, tooltip="Use the labels recorded in this session.", width=150),
                                HGroup(
                                    Item('btnTrainExtern',label='Train with label dataset',show_label=False, tooltip="Use labels loaded from the dataset.", width=150),
                                    Item('extLabelSource', label='Label dataset', tooltip = "Dataset with labels. Stored in '/objects/meta'." ),

                                ),
                                '_',
                                HGroup(
                                    Item('predDestination', label='Write prediction to', tooltip = "Write the prediction in this dataset in '/objects/meta'." ),
                                    Item('predType', label='Prediction type', tooltip = "'Class' predicts only the class for each object, 'Probabilities' predicts the probabilities for each class." ),
                                ),
                                Item('btnPredict',label='Predict',show_label=False, tooltip="Calculate predictions for the opened file."),
                                '_',
                                Item('RFFile', label='Random Forest filename', tooltip = "HDF5 filename of Random Forest file. Random Forest must be stored in '/'." ),
                                HGroup(
                                    Item('btnSaveRF', label='Save Random Forest to file',show_label=False, tooltip = "Write Random Forest to HDF5 file." ),
                                    Item('btnLoadRF', label='Load Random Forest from file',show_label=False, tooltip = "Load Random Forest from HDF5 file." ),
                                ),
                                label = "Train&Predict"
                            ),
                        ),
                     ),
                ),
                title = 'Training GUI',
                resizable=True,
                height=0.9,
                width=0.9,
                )

max_stars_repo_name: m-novikov/hytra

max_issues_repo_name: m-novikov/hytra

max_forks_repo_name: m-novikov/hytra


text:
import itertools

import numpy as np

from sklearn.linear_model import SGDClassifier, SGDRanking
from sklearn import metrics
from minirank.compat import RankSVM as MinirankSVM
from scipy import stats

def transform_pairwise(X, y):
    """Transforms data into pairs with balanced labels for ranking
    Transforms a n-class ranking problem into a two-class classification
    problem. Subclasses implementing particular strategies for choosing
    pairs should override this method.
    In this method, all pairs are choosen, except for those that have the
    same target value. The output is an array of balanced classes, i.e.
    there are the same number of -1 as +1
    Parameters
    ----------
    X : array, shape (n_samples, n_features)
        The data
    y : array, shape (n_samples,) or (n_samples, 2)
        Target labels. If it's a 2D array, the second column represents
        the grouping of samples, i.e., samples with different groups will
        not be considered.
    Returns
    -------
    X_trans : array, shape (k, n_feaures)
        Data as pairs
    y_trans : array, shape (k,)
        Output class labels, where classes have values {-1, +1}
    """
    X_new = []
    y_new = []
    y = np.asarray(y)
    if y.ndim == 1:
        y = np.c_[y, np.ones(y.shape[0])]
    comb = itertools.combinations(range(X.shape[0]), 2)
    for k, (i, j) in enumerate(comb):
        if y[i, 0] == y[j, 0] or y[i, 1] != y[j, 1]:
            # skip if same target or different group
            continue
        X_new.append(X[i] - X[j])
        y_new.append(np.sign(y[i, 0] - y[j, 0]))
        # output balanced classes
        if y_new[-1] != (-1) ** k:
            y_new[-1] = - y_new[-1]
            X_new[-1] = - X_new[-1]
    return np.asarray(X_new), np.asarray(y_new).ravel()

class RankSVM(SGDClassifier):
    """Performs pairwise ranking with an underlying SGDClassifer model
    Input should be a n-class ranking problem, this object will convert it
    into a two-class classification problem, a setting known as
    `pairwise ranking`.
    Authors: Fabian Pedregosa <fabian@fseoane.net>
             Alexandre Gramfort <alexandre.gramfort@inria.fr>
    https://gist.github.com/2071994
    """

    def fit(self, X, y):
        """
        Fit a pairwise ranking model.
        Parameters
        ----------
        X : array, shape (n_samples, n_features)
        y : array, shape (n_samples,) or (n_samples, 2)
        Returns
        -------
        self
        """
        X_trans, y_trans = transform_pairwise(X, y)
        super(RankSVM, self).fit(X_trans, y_trans)
        return self

    def predict(self, X):
        pred = super(RankSVM, self).predict(X)
        # preds are mapped to {-1,1}
        # FIXME only works in this example!!!
        pred[pred == -1] = 0
        return pred

    def score(self, X, y):
        """
        Because we transformed into a pairwise problem, chance level is at 0.5
        """
        X_trans, y_trans = transform_pairwise(X, y)
        return np.mean(super(RankSVM, self).predict(X_trans) == y_trans)

def rank(clf,X):
    if clf.coef_.shape[0] == 1:
        coef = clf.coef_[0]
    else:
        coef = clf.coef_
    order = np.argsort(np.dot(X,coef))
    order_inv = np.zeros_like(order)
    order_inv[order] = np.arange(len(order))
    return order_inv

def kendalltau(clf,X,y):
    if clf.coef_.shape[0] == 1:
        coef = clf.coef_[0]
    else:
        coef = clf.coef_
    tau, _ = stats.kendalltau(np.dot(X, coef), y)
    return np.abs(tau)

if __name__=="__main__":
    rs = np.random.RandomState(0)
    n_samples_1 = 10000
    n_samples_2 = 100
    X = np.r_[1.5 * rs.randn(n_samples_1, 2),
              0.5 * rs.randn(n_samples_2, 2) + [2, 2]]
    y = np.array([0] * (n_samples_1) + [1] * (n_samples_2))
    idx = np.arange(y.shape[0])
    rs.shuffle(idx)
    X = X[idx]
    y = y[idx]
    mean = X.mean(axis=0)
    std = X.std(axis=0)
    X = (X - mean) / std

    for clf, name in ((SGDClassifier(n_iter=100, alpha=0.01), "plain sgd"),
                      (SGDClassifier(n_iter=100, alpha=0.01,
                                     class_weight={1: 10}),"weighted sgd"),
                      (SGDRanking(n_iter=1000, alpha=0.01,
                                     loss='roc_pairwise_ranking'), "pairwise sgd"),
                      (RankSVM(n_iter=100, alpha=0.01, loss='hinge'), 'RankSVM'),
                      ):
        clf.fit(X, y)
        print clf
        pred = clf.predict(X)

        print "ACC: %.4f" % metrics.zero_one_score(y, pred)
        print "AUC: %.4f" % metrics.auc_score(y, pred)
        print "CONFUSION MATRIX: "
        print metrics.confusion_matrix(y, pred)
        print "Kendall Tau: %.4f" % kendalltau(clf,X,y)
        print 80*'='

    clf = MinirankSVM(max_iter=100, alpha=0.01).fit(X,y)
    print clf
    scores = np.dot(X,clf.coef_)
    pred = (scores > 0).astype(np.int)
    print "ACC: %.4f" % metrics.zero_one_score(y, pred)
    print "AUC: %.4f" % metrics.auc_score(y, pred)
    print "CONFUSION MATRIX: "
    print metrics.confusion_matrix(y, pred)
    print "Kendall Tau: %.4f" % kendalltau(clf,X,y)
    print 80*'='

reason: import numpy,from scipy


text:
import matplotlib
matplotlib.use('Agg')
from matplotlib import rcParams
import matplotlib.pyplot as plt
from config import PLOT_PATH, TAB_PATH, DATABASE_PATH, normal_dict, dark_dict
from config import order_caller, order_model, order_strategy, order_mutation, order_filter
from utils import dfs_to_sheet, noparam_wrm, get_pal, nanlog10
import seaborn as sns
import pandas as pd
import numpy as np
import os
import sys

patient = 'NA12877'

if len(sys.argv) > 1:
    patient = sys.argv[1]

# setting plot parameters
snsdeep = sns.color_palette('deep', 10)
snsdark = sns.color_palette('dark', 10)
snsset2 = sns.color_palette('Set2', 10)

sns.set_theme(
        style='ticks',
        context='talk',
        palette='deep')
params = {
    # 'font.family':'serif',
    # 'font.serif':'Times New Roman',
    'font.family': 'sans-serif',
    'font.sans-serif': 'Arial',
    'font.style':'normal',
    'font.weight':'normal', #or 'blod'
    'pdf.fonttype':42
    # 'font.size':'medium',#or large,small
}
rcParams.update(params)

# for leave-one-patient-out setting
alias_dict = {**order_caller,**order_model,**order_strategy,**order_mutation,**order_filter}

## metr table
metr_df = pd.read_pickle(os.path.join(DATABASE_PATH, 'rawfeat_filter_metr_df_additional.pkl'))
select_row = metr_df['patient']==patient
metr_df = metr_df.loc[select_row]
metr_df.insert(metr_df.columns.to_list().index('TNR')+1, 'TPR', 1 - metr_df['FNR'])
metr_df.insert(metr_df.columns.to_list().index('OFO')+1, 'log OFO', [nanlog10(x) for x in metr_df['OFO'].to_numpy()])
metr_df[['Caller','Mutation','Filter']] = metr_df[['software','mutation','filter']].applymap(lambda x: alias_dict[x])
metr_df_pp = metr_df.groupby(by=['depth','Caller','Mutation','Filter','patient']).mean()
metr_df_mn = metr_df.groupby(by=['depth','Caller','Mutation','Filter']).mean()
metr_df_sd = metr_df.groupby(by=['depth','Caller','Mutation','Filter']).std()
dfs_to_sheet([metr_df, metr_df_pp, metr_df_mn, metr_df_sd], ['raw', 'per patient', 'mean', 'std'], filename='rawfeat_filter_{}_mutation.xlsx'.format(patient), root_path=TAB_PATH)

## roc plot
roc_df = pd.read_pickle(os.path.join(DATABASE_PATH, 'rawfeat_filter_roc_df_additional.pkl'))
select_row = roc_df['patient']==patient
roc_df = roc_df.loc[select_row]
roc_df[['Caller','Mutation','Filter']] = roc_df[['software','mutation','filter']].applymap(lambda x: alias_dict[x])

for depth in ['50X']:
    g = sns.relplot(data=metr_df.loc[metr_df['depth']==depth], row='Mutation', col='Caller', kind='scatter', height=5, aspect=1,
        x='FPR', y='TPR', style='Filter', markers=['x']*len(set(metr_df['Filter'])),
        hue='Filter', hue_order=['FNVC', 'VEF', 'Frequency', 'Garfield', 'VQSR', 'Hard Filter'],
        palette=[snsdeep[i] for i in [0,9,1,3,4,8]])
    g.set(xlim=(-0.1,1.1), ylim=(-0.1,1.1))
    plt.savefig(os.path.join(PLOT_PATH, 'rawfeat_filter_{}_mutation_roc_scatter_{}.pdf'.format(patient, depth)))
    # plt.show()

    g = sns.relplot(data=roc_df.loc[roc_df['depth']==depth], row='Mutation', col='Caller', kind='line', height=5, aspect=1,
        x='fpr', y='tpr',
        hue='Filter', hue_order=['FNVC', 'VEF', 'Frequency', 'Garfield', 'VQSR', 'Hard Filter'],
        palette=[snsdeep[i] for i in [0,9,1,3,4,8]])
    g.set(xlim=(-0.1,1.1), ylim=(-0.1,1.1))
    plt.savefig(os.path.join(PLOT_PATH, 'rawfeat_filter_{}_mutation_roc_line_{}.pdf'.format(patient, depth)))
    # plt.show()

## bar plot
for depth in ['50X']:
    for y_metrname in ['AUC', 'G1-score', 'MCC', 'log OFO']:
        sns.catplot(
            data=metr_df.loc[metr_df['depth']==depth], x='Caller', y=y_metrname,
            hue='Filter', row='Mutation',
            order=['GATK', 'Varscan2', 'Mutect2', 'DeepVariant'],
            hue_order=['FNVC', 'VEF', 'Frequency', 'Garfield', 'VQSR', 'Hard Filter'],
            palette=[snsdeep[i] for i in [0,1,3,4,8,9]],
            kind='bar', aspect=1.6, height=4, errwidth=1, capsize=0.05)
        plt.savefig(os.path.join(PLOT_PATH, 'rawfeat_filter_{}_mutation_bar_{}-{}.pdf'.format(patient, depth, y_metrname.replace(' ', '.'))))
        # plt.show()

## pvalue
test_tab_df_list = []
test_tab_name_list = []
for depth in ['50X']:
    for y_metrname in ['AUC', 'G1-score', 'MCC', 'log OFO']:
        tab_list = []
        for mutation in set(metr_df['Mutation']):
            for caller in set(metr_df['Caller']):
                for filter in set(metr_df['Filter']).difference(['FNVC']):
                    x = metr_df.loc[
                        (metr_df['depth']==depth)&
                        (metr_df['Mutation']==mutation)&
                        (metr_df['Caller']==caller)&
                        (metr_df['Filter']=='FNVC')][y_metrname].values
                    y = metr_df.loc[
                        (metr_df['depth']==depth)&
                        (metr_df['Mutation']==mutation)&
                        (metr_df['Caller']==caller)&
                        (metr_df['Filter']==filter)][y_metrname].values
                    print(x.shape, y.shape)
                    if x.shape[0] == y.shape[0]:
                        tab_list.append([mutation, caller, filter, *noparam_wrm(x,y)])
        test_tab_df_list.append(pd.DataFrame(np.array(tab_list), columns=['Mutation','Caller', 'Filter', 'FNVC Ave','FNVC SD','Filter Ave','Filter SD','statistic','pvalue']))
        test_tab_name_list.append('{}-{}'.format(depth,y_metrname))
dfs_to_sheet(test_tab_df_list, test_tab_name_list, 'rawfeat_filter_{}_mutation_pvalue.xlsx'.format(patient), TAB_PATH)

max_stars_repo_path: Code for Reproducing Results/Plot_figures/11.rawfeat_aditional_mutation.py

max_stars_repo_name: yyren/FVC

max_issues_repo_name: yyren/FVC

max_forks_repo_name: yyren/FVC


text:
#!/usr/bin/env python
from __future__ import print_function
from optparse import OptionParser
import copy, os, pdb, subprocess, sys

import h5py
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
import seaborn as sns

import basset_sat
import dna_io
from seq_logo import seq_logo
import bvcf

################################################################################
# basset_sat_vcf.py
#
# Perform an in silico saturated mutagenesis of the regions surrounding a list
# of SNPs given in VCF format using the given model.
################################################################################

################################################################################
# main
############################s####################################################
def main():
    usage = 'usage: %prog [options] <model_file> <vcf_file>'
    parser = OptionParser(usage)
    parser.add_option('-d', dest='model_hdf5_file', default=None, help='Pre-computed model output as HDF5 [Default: %default]')
    parser.add_option('-f', dest='genome_fasta', default='%s/data/genomes/hg19.fa'%os.environ['BASSETDIR'], help='Genome FASTA from which sequences will be drawn [Default: %default]')
    parser.add_option('--f1', dest='genome1_fasta', default=None, help='Genome FASTA which which major allele sequences will be drawn')
    parser.add_option('--f2', dest='genome2_fasta', default=None, help='Genome FASTA which which minor allele sequences will be drawn')
    parser.add_option('-g', dest='gain_height', default=False, action='store_true', help='Nucleotide heights determined by the max of loss and gain [Default: %default]')
    parser.add_option('-l', dest='seq_len', type='int', default=600, help='Sequence length provided to the model [Default: %default]')
    parser.add_option('-m', dest='min_limit', default=0.1, type='float', help='Minimum heatmap limit [Default: %default]')
    parser.add_option('-n', dest='center_nt', default=200, type='int', help='Nt around the SNP to mutate and plot in the heatmap [Default: %default]')
    parser.add_option('-o', dest='out_dir', default='heat', help='Output directory [Default: %default]')
    parser.add_option('-t', dest='targets', default='-1', help='Comma-separated list of target indexes to plot (or -1 for all) [Default: %default]')
    (options,args) = parser.parse_args()

    if len(args) != 2:
        parser.error('Must provide Basset model file and input SNPs in VCF format')
    else:
        model_file = args[0]
        vcf_file = args[1]

    if not os.path.isdir(options.out_dir):
        os.mkdir(options.out_dir)

    #################################################################
    # prep SNP sequences
    #################################################################
    # load SNPs
    snps = bvcf.vcf_snps(vcf_file, pos2=(options.genome2_fasta is not None))

    # get one hot coded input sequences
    if not options.genome1_fasta or not options.genome2_fasta:
        seqs_1hot, seq_headers, snps, seqs = bvcf.snps_seq1(snps, options.seq_len, options.genome_fasta, return_seqs=True)
    else:
        seqs_1hot, seq_headers, snps, seqs = bvcf.snps2_seq1(snps, options.seq_len, options.genome1_fasta, options.genome2_fasta, return_seqs=True)

    # reshape sequences for torch
    seqs_1hot = seqs_1hot.reshape((seqs_1hot.shape[0],4,1,seqs_1hot.shape[1]//4))

    # write to HDF5
    model_input_hdf5 = '%s/model_in.h5'%options.out_dir
    h5f = h5py.File(model_input_hdf5, 'w')
    h5f.create_dataset('test_in', data=seqs_1hot)
    h5f.close()

    #################################################################
    # Torch predict modifications
    #################################################################
    if options.model_hdf5_file is None:
        options.model_hdf5_file = '%s/model_out.h5' % options.out_dir
        torch_cmd = '%s/src/basset_sat_predict.lua -center_nt %d %s %s %s' % (os.environ['BASSETDIR'], options.center_nt, model_file, model_input_hdf5, options.model_hdf5_file)
        subprocess.call(torch_cmd, shell=True)

    #################################################################
    # load modification predictions
    #################################################################
    hdf5_in = h5py.File(options.model_hdf5_file, 'r')
    seq_mod_preds = np.array(hdf5_in['seq_mod_preds'])
    hdf5_in.close()

    # trim seqs to match seq_mod_preds length
    delta_start = 0
    delta_len = seq_mod_preds.shape[2]
    if delta_len < options.seq_len:
        delta_start = (options.seq_len - delta_len)//2
        for si in range(len(seqs)):
            seqs[si] = seqs[si][delta_start:delta_start+delta_len]

    # decide which cells to plot
    if options.targets == '-1':
        plot_cells = xrange(seq_mod_preds.shape[3])
    else:
        plot_cells = [int(ci) for ci in options.targets.split(',')]

    #################################################################
    # plot
    #################################################################
    table_out = open('%s/table.txt' % options.out_dir, 'w')

    rdbu = sns.color_palette("RdBu_r", 10)

    nts = 'ACGT'
    for si in range(seq_mod_preds.shape[0]):
        header = seq_headers[si]
        header_fs = fs_clean(header)
        seq = seqs[si]

        # plot some descriptive heatmaps for each individual cell type
        for ci in plot_cells:
            seq_mod_preds_cell = seq_mod_preds[si,:,:,ci]
            real_pred_cell = basset_sat.get_real_pred(seq_mod_preds_cell, seq)

            # compute matrices
            norm_matrix = seq_mod_preds_cell - real_pred_cell
            min_scores = seq_mod_preds_cell.min(axis=0)
            max_scores = seq_mod_preds_cell.max(axis=0)
            minmax_matrix = np.vstack([min_scores - real_pred_cell, max_scores - real_pred_cell])

            # prepare figure
            sns.set(style='white', font_scale=0.5)
            sns.axes_style({'axes.linewidth':1})
            spp = basset_sat.subplot_params(seq_mod_preds_cell.shape[1])
            fig = plt.figure(figsize=(20,3))
            ax_logo = plt.subplot2grid((3,spp['heat_cols']), (0,spp['logo_start']), colspan=(spp['logo_end']-spp['logo_start']))
            ax_sad = plt.subplot2grid((3,spp['heat_cols']), (1,spp['sad_start']), colspan=(spp['sad_end']-spp['sad_start']))
            ax_heat = plt.subplot2grid((3,spp['heat_cols']), (2,0), colspan=spp['heat_cols'])

            # print a WebLogo of the sequence
            vlim = max(options.min_limit, abs(minmax_matrix).max())
            if options.gain_height:
                seq_heights = 0.25 + 1.75/vlim*(abs(minmax_matrix).max(axis=0))
            else:
                seq_heights = 0.25 + 1.75/vlim*(-minmax_matrix[0])
            logo_eps = '%s/%s_c%d_seq.eps' % (options.out_dir, header_fs, ci)
            seq_logo(seq, seq_heights, logo_eps, color_mode='meme')

            # add to figure
            logo_png = '%s.png' % logo_eps[:-4]
            subprocess.call('convert -density 300 %s %s' % (logo_eps, logo_png), shell=True)
            logo = Image.open(logo_png)
            ax_logo.imshow(logo)
            ax_logo.set_axis_off()

            # plot loss and gain SAD scores
            ax_sad.plot(-minmax_matrix[0], c=rdbu[0], label='loss', linewidth=1)
            ax_sad.plot(minmax_matrix[1], c=rdbu[-1], label='gain', linewidth=1)
            ax_sad.set_xlim(0,minmax_matrix.shape[1])
            ax_sad.legend()
            # ax_sad.grid(True, linestyle=':')
            for axis in ['top','bottom','left','right']:
                ax_sad.spines[axis].set_linewidth(0.5)

            # plot real-normalized scores
            vlim = max(options.min_limit, abs(norm_matrix).max())
            sns.heatmap(norm_matrix, linewidths=0, cmap='RdBu_r', vmin=-vlim, vmax=vlim, xticklabels=False, ax=ax_heat)
            ax_heat.yaxis.set_ticklabels('TGCA', rotation='horizontal') # , size=10)

            # save final figure
            plt.tight_layout()
            plt.savefig('%s/%s_c%d_heat.pdf' % (options.out_dir, header_fs, ci), dpi=300)
            plt.close()

        #################################################################
        # print table of nt variability for each cell
        #################################################################
        for ci in range(seq_mod_preds.shape[3]):
            seq_mod_preds_cell = seq_mod_preds[si,:,:,ci]
            real_pred_cell = basset_sat.get_real_pred(seq_mod_preds_cell, seq)

            min_scores = seq_mod_preds_cell.min(axis=0)
            max_scores = seq_mod_preds_cell.max(axis=0)

            loss_matrix = real_pred_cell - seq_mod_preds_cell.min(axis=0)
            gain_matrix = seq_mod_preds_cell.max(axis=0) - real_pred_cell

            for pos in range(seq_mod_preds_cell.shape[1]):
                cols = [header, delta_start+pos, ci, loss_matrix[pos], gain_matrix[pos]]
                print('\t'.join([str(c) for c in cols]), file=table_out)

    table_out.close()

def fs_clean(header):
    ''' Clean up the headers to valid filenames. '''
    header = header.replace(':','_')
    header = header.replace('>','_')
    return header

################################################################################
# __main__
################################################################################
if __name__ == '__main__':
    main()
    #pdb.runcall(main)


text:
import numpy as module_np

def function_tsig(
    list_indicator
):

    '''
    This function accepts a list of
    indicator values and converts it
    to a trade signal series as a list.
    '''

    # Transform the entries of the
    # indicator series to create the
    # trade signal series 'list_sign'
    list_sign = [

        # Pass through 'sign' function
        module_np.sign(

            # Pass through 'maximum' function
            module_np.maximum(
                float_value,
                0
            )
        )

        # Iterate for all values in 'list_indicator'
        for float_value
        in list_indicator

    ]

    # Output 'list_sign'
    return list_sign

def function_action(
    list_sign
):

    '''
    This function accepts a list
    representing a trade signal series
    and transforms it to an action
    series.
    '''

    # Shift 'list_sign' by one step
    # and fill the first value with zero
    list_shift = [0] + list_sign[:-1]

    # Subtract entries 'list_shift'
    # from 'list_sign' element-wise
    list_output = [

        list_sign[int_i] -
        list_shift[int_i]

        for int_i
        in range(
            0,
            len(list_sign)
        )
    ]

    # Return final action series
    # as a list
    return list_output

def function_backtester(
    list_prices,
    list_indicator,
    float_ratio
):

    '''
    Backtester for either a long or short position.
    Long positions are tested with positive 'float_ratio'.
    Short positions are tested with negative 'float_ratio'.
    '''

    # Derive the trade signal series
    # from the indicator series in the form
    # of 'list_inidcator'
    list_tsig = function_tsig(
        list_indicator = list_indicator
    )

    # Derive the action series from
    # the previously derived trade
    # signal series
    list_action = function_action(
        list_sign = list_tsig
    )

    # Create a list for storing the
    # alpha series of the trading
    # strategy
    list_alpha = []

    # Set the length of the list of
    # prices
    int_length = len(list_prices)

    # Carry out the algorithm to
    # generate the alpha series
    for int_i in range(0, int_length):

        if list_action[int_i] == 1:

            list_alpha.append(1)
            float_inprice = list_prices[int_i]

        if list_action[int_i] == 0:

            if list_tsig[int_i] == 1:

                list_alpha.append(
                    (float_inprice + (list_prices[int_i] - float_inprice) * float_ratio) /
                    (float_inprice + (list_prices[int_i - 1] - float_inprice) * float_ratio)
                )

            if list_tsig[int_i] == 0:

                list_alpha.append(1)

        if list_action[int_i] == -1:

            list_alpha.append(
                (float_inprice + (list_prices[int_i] - float_inprice) * float_ratio) /
                (float_inprice + (list_prices[int_i - 1] - float_inprice) * float_ratio)
            )

    # Return the final alpha series
    return list_alpha

def function_heikin_ashi_gen(
    list_values,
    float_initial,
    float_strength
):

    list_output = [
        float_initial
    ]

    int_length = len(
        list_values
    )

    for int_i in range(0, int_length - 1):

        float_new_ha = list_output[int_i] * float_strength + list_values[int_i + 1] - list_values[int_i]

        list_output.append(
            float_new_ha
        )

    return list_output