text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Oct 21 11:38:47 2021

@author: mariaolaru
SVM implementation
"""

import pandas as pd
import numpy as np
import glob
import os
from matplotlib import pyplot as plt
pd.options.display.float_format = '{:.2f}'.format
#pd.reset_option('display.float_format')

def import_pkg_files(pkg_dir):
    paths = glob.glob(pkg_dir + 'scores*')
    gp = os.path.abspath(os.path.join(paths[0], "../.."))
    subj_id = os.path.basename(gp)

    md_fp = pkg_dir + subj_id + '_pkg_data_combined.csv'
    md = pd.DataFrame([])

    if (os.path.exists(md_fp)):
        md = pd.read_csv(md_fp, header=0)
    else:
        for i in range(0, len(paths)):
            p = paths[i]
            df = pd.read_csv(p)
            df['timestamp'] = pd.to_datetime(df['Date_Time']).astype(int) / 10**9
            md = pd.concat([md, df])

        md = md.sort_values('timestamp')
        md = md.reset_index(drop=True)
        md.to_csv(md_fp, index = False, header=True)
        md['timestamp'] = md['timestamp'].astype(int)
        md['timestamp'] = md['timestamp']*1000
    return md

def normalize_data(x, min_x, max_x):
    x_norm = (x-min_x)/(max_x-min_x)
    return x_norm

def remove_outliers(col):
    x = col.copy()
    y = x.median() + 5*x.std()
    indx = x[x > y].index
    x[indx] = np.nan
    return x

def preproc_pkg(pkg_dir):
    pkg_df = import_pkg_files(pkg_dir)
    pkg_df['timestamp'] = pkg_df['timestamp'].astype(int)*1000
    pkg_df['timestamp'] = pkg_df['timestamp'] +  60000 #add 1min to end time
    pkg_df['timestamp'] = pkg_df['timestamp'] + 25200000 #add 7hr to end time (b/c Datetime in PT, not GMT)

    #Include data where watch is on wrist, from even indices
    pkg_df = pkg_df[(pkg_df['Off_Wrist'] == 0)]

    start_time = pkg_df['timestamp'].head(1).values[0]
    stop_time = pkg_df['timestamp'].tail(1).values[0]

    return [pkg_df, start_time, stop_time]

def preproc_phs(phs_fp, start_time, stop_time):
    phs_df = pd.read_csv(phs_fp)
    phs_df = phs_df.rename({'timestamp_end': 'timestamp'}, axis = 1)
    phs_df['timestamp'] = phs_df['timestamp'].astype(int)
    #subset data
    phs_df = phs_df[(phs_df['timestamp'] > start_time) & (phs_df['timestamp'] < stop_time)]

    #create wide phs df
    phs_dfw = phs_df.pivot_table(index = ['timestamp', 'contacts'], values = ['max_amp'], columns = ['band'])
    phs_dfw.columns = [f'{x}_{y}' for x,y in phs_dfw.columns]
    phs_dfw = phs_dfw.reset_index(level=['timestamp', 'contacts'])
    phs_dfw = phs_dfw.sort_values('timestamp')

    phs_dfwb = phs_dfw[phs_dfw['contacts'] == '+3-1']
    phs_dfwb = pd.concat([phs_dfwb['timestamp'], phs_dfwb['max_amp_beta']], axis = 1)

    phs_dfwg = phs_dfw[phs_dfw['contacts'] == '+11-9']
    phs_dfwg = pd.concat([phs_dfwg['timestamp'], phs_dfwg['max_amp_gamma']], axis = 1)

    #Merge phs tables
    phs_merged = pd.merge(phs_dfwb, phs_dfwg, how = 'inner', on = 'timestamp')
    return phs_merged

def preproc_psd(fp_psd, start_time, stop_time):
    psd_df = pd.read_csv(fp_psd)
    psd_df = psd_df.rename({'timestamp_end': 'timestamp'}, axis = 1)
#    psd_df['timestamp'] = psd_df['timestamp'].astype(int)*1000
    #subset data
    psd_df = psd_df[(psd_df['timestamp'] > start_time) & (psd_df['timestamp'] < stop_time)]

    #create wide phs df
    psd_dfw = psd_df.pivot_table(index = ['timestamp'], values = ['spectra'], columns = ['f_0', 'contacts'])
    psd_dfw = psd_dfw['spectra']

    return psd_dfw

def merge_pkg_df(df_pkg, df):
    #Both dfs must have 'timestamp' column to merge on
    df_merged = pd.merge(df_pkg, df, how = 'inner', on = 'timestamp')
    df_merged = df_merged.sort_values('timestamp')
    df_merged = df_merged.reset_index(drop=True)
    return df_merged

def process_pkg(df_merged):
    df_merged['BK_rev'] = df_merged['BK']*-1
    df_merged['BK_rev'] = df_merged['BK_rev'] + (df_merged['BK_rev'].min()*-1)

    #Remove outliers from each max_amp column separately
    df_merged['DK'] = remove_outliers(df_merged['DK'])
    df_merged['BK_rev'] = remove_outliers(df_merged['BK_rev'])

    #Normalize data from 0-1
    df_merged['BK_norm'] = normalize_data(df_merged['BK_rev'], 0, np.nanmax(df_merged['BK_rev']))
    df_merged['DK_norm'] = normalize_data(df_merged['DK'], 0, np.nanmax(df_merged['DK']))
    return df_merged

def find_nan_chunks(col, num_nans):
    #find nan values
    x = col.isnull().astype(int).groupby(col.notnull().astype(int).cumsum()).cumsum()
    #find indices of nan vals >= num_nan
    indx_remove = x[x >= num_nans].index
    return indx_remove

def process_psd(df_merged):
    #df_merged['max_amp_beta_norm'] = normalize_data(df_merged['max_amp_beta'], np.nanmin(df_merged['max_amp_beta']), np.nanmax(df_merged['max_amp_beta']))
    #df_merged['max_amp_gamma_norm'] = normalize_data(df_merged['max_amp_gamma'], np.nanmin(df_merged['max_amp_gamma']), np.nanmax(df_merged['max_amp_gamma']))
    #df_merged['max_amp_diff_norm'] = df_merged['max_amp_beta_norm'] - df_merged['max_amp_gamma_norm']
    #df_merged['max_amp_diff_norm'] = normalize_data(df_merged['max_amp_diff'], 0, np.nanmax(df_merged['max_amp_diff']))
    cols = np.array(df_merged.columns)
    indices = [i for i, s in enumerate(list(cols)) if '+' in s]
    for i in indices:
        col_name = cols[i]
        df_merged[col_name] = normalize_data(df_merged[col_name], np.nanmin(df_merged[col_name]), np.nanmax(df_merged[col_name]))
    return df_merged

def process_dfs(df_pkg, df_lfp):
    df_merged = merge_pkg_df(df_pkg, df_lfp)

    df_merged = process_pkg(df_merged)
    df_merged = process_psd(df_merged)
    #Remove data with large nan blocks
    #indx_nan_chunks = find_nan_chunks(df_merged['+10-8'], 3)
    #indx_nan_chunks = find_nan_chunks(df_merged.max_amp_gamma, 3)
    #df_merged_trim = df_merged.drop(indx_nan_chunks, axis = 0)

    #indx_nan_chunks_2 = find_nan_chunks(df_merged_trim.DK, 6)
    #df_merged_trim = df_merged_trim.drop(indx_nan_chunks_2, axis = 0)
    #Fill remaining nan values
    #df_mergedf = df_merged_trim.ffill().reset_index(drop=True)
    #df_mergedfb = df_mergedf.bfill().reset_index(drop=True)

    return df_merged

def plot_pkg_sync(df_merged, freq_band):
    #Plotting
    #plt.plot(np.arange(1, len(df_mergedfb)+1, 1), df_mergedfb['BK_norm'], alpha = 0.7, label = 'PKG-BK', markersize = 1, color = 'indianred')

    #plt.plot(np.arange(1, len(df_mergedfb)+1, 1), df_mergedfb['max_amp_gamma_norm'], alpha = 0.7, label = 'RCS-gamma', markersize = 1, color = 'darkorange')
    #plt.plot(np.arange(1, len(df_mergedfb)+1, 1), df_mergedfb['max_amp_beta_norm'], alpha = 0.7, label = 'RCS-beta', markersize = 1, color = 'mediumpurple')
    title = ("freq_band: " + str(freq_band) + "Hz")
    plt.title(title)
    plt.rcParams["figure.figsize"] = (30,3)
    #plt.plot(np.arange(1, len(df_merged)+1, 1), df_merged[(contacts[0]+'_norm')], alpha = 0.7, label = contacts[0], markersize = 1, color = 'orchid')
    plt.plot(np.arange(1, len(df_merged)+1, 1), df_merged[(contacts[1] + '_norm')], alpha = 0.7, label = contacts[1], markersize = 1, color = 'mediumpurple')
    #plt.plot(np.arange(1, len(df_merged)+1, 1), df_merged[(contacts[2] + '_norm')], alpha = 0.7, label = contacts[2], markersize = 1, color = 'darkkhaki')
    plt.plot(np.arange(1, len(df_merged)+1, 1), df_merged[(contacts[3]+'_norm')], alpha = 0.7, label = contacts[3], markersize = 1, color = 'darkorange')

    plt.plot(np.arange(1, len(df_merged)+1, 1), df_merged['DK_norm'], alpha = 0.7, label = 'PKG-DK', markersize = 1, color = 'steelblue')

    #plt.plot(df_mergedfb['timestamp'], df_mergedfb['max_amp_diff_norm'], alpha = 0.7, markersize = 1, label = 'RCS (beta-gamma)')
    plt.legend(ncol = 5, loc = 'upper right')
    plt.ylabel('scores (normalized)')
    plt.xlabel('time (samples)')

    out_dir = '/Users/mariaolaru/Documents/temp/RCS02/RCS02L_pkg_rcs/plots'
    plt.savefig(out_dir + '/' + 'psd_' + 'freq' + str(freq_band) + '.pdf')
    plt.close()

def plot_corrs(df_corr, PKG_measure):
    df_corr_s = df_corr.loc[PKG_measure, :]
    plt.rcParams["figure.figsize"] = (5,5)
    plt.plot(df_corr_s['freq_band'], df_corr_s['+2-0'], label = '+2-0', color = 'orchid')
    plt.plot(df_corr_s['freq_band'], df_corr_s['+3-1'], label = '+3-1', color = 'mediumpurple')
    plt.plot(df_corr_s['freq_band'], df_corr_s['+10-8'], label = '+10-8', color = 'darkkhaki')
    plt.plot(df_corr_s['freq_band'], df_corr_s['+11-9'], label = '+11-9', color = 'darkorange')
    plt.title(PKG_measure)
    plt.legend(ncol = 2, loc = 'upper right')
    plt.ylim(-0.3, 0.3)
    plt.ylabel('Pearson coef (r)')
    plt.xlabel('Frequency (Hz)')

    out_dir = '/Users/mariaolaru/Documents/temp/RCS02/RCS02L_pkg_rcs/plots'
    plt.savefig(out_dir + '/' + 'psd_corr_' + 'PKG_' + PKG_measure + '.pdf')
    plt.close()

pkg_dir = '/Users/mariaolaru/Documents/temp/RCS02/RCS02L/RCS02L_ge_onmed/RCS02R_pkg_data/'
[df_pkg, start_time, stop_time] = preproc_pkg(pkg_dir)

#fp_phs = '/Users/mariaolaru/Documents/temp/RCS02/RCS02L/RCS02L_pre-stim/RCS02L_pre-stim_phs.csv'
#df_phs = preproc_phs(fp_phs, start_time, stop_time)

fp_psd = '/Users/mariaolaru/Documents/temp/RCS02/RCS02L/RCS02L_pre-stim/RCS02L_pre-stim_psd_partial_2m_wholemin.csv'
df_psd = preproc_psd(fp_psd, start_time, stop_time)
#df_psd.columns = [''.join(str(col)) for col in df_psd.columns]
#df_psd = df_psd.reset_index()
#df_merged = process_dfs(df_pkg, df_psd)
contacts = ['+2-0', '+3-1', '+10-8', '+11-9']

#Processing
freq_bands = np.array(df_psd.columns.levels[0])
df_corr = pd.DataFrame([])
for i in range(len(freq_bands)):
    freq_band = freq_bands[i]
    print("freq band: " + str(freq_band) + "hz")

    df_merged = process_dfs(df_pkg, df_psd[freq_band])

    df_merged['timestamp_dt'] = pd.to_datetime(df_merged['Date_Time'])
    df_merged['timestamp_dt_h'] = [i.hour for i in df_merged['timestamp_dt']]
    df_merged['asleep'] = (df_merged['timestamp_dt_h'] >= 20) | (df_merged['timestamp_dt_h'] <= 8)
    df_merged['asleep'] = df_merged['asleep'].astype(int)

    df_merged = df_merged[df_merged['asleep'] == 0]

    plot_pkg_sync(df_merged, freq_band)

    #get corr coefs
    df_vars = df_merged.loc[:, ['BK_norm', 'DK_norm', (contacts[3]), (contacts[2]), (contacts[1]), (contacts[0])]]

    #df_vars = df_mergedfb.loc[:, ['BK_norm', 'DK_norm', 'max_amp_beta', 'max_amp_gamma']]
    df_corr_ind = df_vars.corr()

    out_fp = '/Users/mariaolaru/Documents/temp/RCS02/RCS02L_pkg_rcs' + '/RCS02_corrs_' + str(freq_band) + 'hz' + '.csv'

    df_corr_ind = df_corr_ind.round(3)
    df_corr_ind['freq_band'] = freq_band
    df_corr = pd.concat([df_corr, df_corr_ind])
df_corr.to_csv(out_fp)
print(df_corr)
print('/n')

#get correlation plots
plot_corrs(df_corr, 'DK_norm')
plot_corrs(df_corr, 'BK_norm')

"""
#Run an SVM
from sklearn.utils import Bunch
from sklearn.model_selection import train_test_split
from sklearn import svm
from sklearn import metrics

#Normalized DK scores > 0.2 == Dyskinesia
df_merged['DK_class'] = 0
indxs = np.where(df_merged['DK_norm'] > 0.1)[0]
df_merged.loc[indxs, 'DK_class'] = 1

cols = np.array(df_merged.columns)
feature_col_indxs = [i for i, s in enumerate(list(cols)) if '+' in s]
df_data = df_merged.iloc[:, feature_col_indxs]

data = df_data.to_numpy()
feature_names = df_data.columns
target = df_merged['DK_class'].to_numpy()
target_names = np.array(['present', 'absent'])

dyskinesia = Bunch(data = data, feature_names = feature_names, target = target, target_names = target_names)

accuracy = np.array([])
precision = np.array([])
recall = np.array([])
for i in range(100):
    print(str(i))
    # Split dataset into training set and test set
    X_train, X_test, y_train, y_test = train_test_split(dyskinesia.data, dyskinesia.target, test_size=0.2)

    #Create classifier
    clf = svm.SVC(kernel='linear') # Linear Kernel

    #Train the model using the training sets
    clf.fit(X_train, y_train)

    #Predict the response for test dataset
    y_pred = clf.predict(X_test)

    accuracy = np.append(accuracy, metrics.accuracy_score(y_test, y_pred))
    precision = np.append(precision, metrics.precision_score(y_test, y_pred))
    recall = np.append(recall, metrics.recall_score(y_test, y_pred))

print("Accuracy:", str(accuracy.mean()))
print("Precision:", str(precision.mean()))
print("Recall:", str(recall.mean()))

plt.scatter(np.arange(0,len(y_pred)), y_pred, alpha = 0.7, label = 'predict')
plt.scatter(np.arange(0,len(y_pred)), y_test, alpha = 0.7, label = 'test')
plt.legend()

#Smooth data
df_mergedfb['DK_smoo'] = df_mergedfb['DK_norm'].rolling(window = 2, win_type = 'gaussian', center = True).mean(std=0.5)
df_mergedfb['max_amp_gamma_smoo'] = df_mergedfb['max_amp_gamma_norm'].rolling(window = 8, win_type = 'gaussian', center = True).mean(std=0.5)

plt.plot(np.arange(1, len(df_mergedfb)+1, 1), df_mergedfb['max_amp_gamma_norm'], alpha = 0.7, label = 'RCS-gamma', markersize = 1, color = 'darkorange')
plt.plot(np.arange(1, len(df_mergedfb)+1, 1), df_mergedfb['max_amp_gamma_smoo'], alpha = 0.7, label = 'smoo', markersize = 1, color = 'red')
plt.legend(ncol = 4, loc = 'upper right')
plt.ylabel('scores (normalized)')
plt.xlabel('time (samples)')
"""


text:
import os
import sys
import random
import math
import re
import time
import numpy as np
import cv2

ROOT_DIR = os.path.abspath("../../")

sys.path.append(ROOT_DIR)  # To find local version of the library
from mrcnn.config import Config
from mrcnn import utils
import mrcnn.model as modellib
from mrcnn import visualize
from mrcnn.model import log

# Directory to save logs and trained model
MODEL_DIR = os.path.join(ROOT_DIR, "logs")

class YCBVDataset(utils.Dataset):
    """ Loads YCB-Video dataset for Mask RCNN training
    """
    def __init__(self, data_dir, split=None, name='ycb'):
        super(YCBVDataset, self).__init__()
        self.data_dir = data_dir
        self.split = split
        self.name = name
        self.mask_files = self.set_all_mask_files()

    def set_all_mask_files(self):
        mask_files = {}
        masks_dir = os.path.join(self.data_dir, self.split + '/mask')
        for file in os.listdir(masks_dir):
            mask_file = masks_dir + '/' + file
            IMID, class_id = file.split('_')
            class_id, _ = class_id.split('.')
            mask_files.setdefault(IMID, {}).update({class_id: mask_file})
        return mask_files

    def load_ycbv(self):
        class_names = self.get_class_names()
        for class_id in range(1, len(class_names)):
            self.add_class('ycb', class_id, class_names[class_id])
        assert self.split in ['train', 'val', 'test']
        images = {}
        images_dir = os.path.join(self.data_dir, self.split + '/rgb')
        for filename in os.listdir(images_dir):
            IMID, _ = filename.split('.')
            images[IMID] = images_dir + '/' + filename
        images = dict(sorted(images.items(), key=lambda x: x[1]))
        for image_id in list(images.keys()):
            self.add_image('ycb', image_id, path=images[image_id])
        self.image_info = sorted(self.image_info, key=lambda k: int(k['id']))

    def get_class_names(self):
        class_names = ['background', '002_master_chef_can', '003_cracker_box',
                       '004_sugar_box', '005_tomato_soup_can', '006_mustard_bottle',
                       '007_tuna_fish_can', '008_pudding_box', '009_gelatin_box',
                       '010_potted_meat_can', '011_banana', '019_pitcher_base',
                       '021_bleach_cleanser', '024_bowl',
                       '025_mug', '035_power_drill',
                       '036_wood_block', '037_scissors',
                       '040_large_marker', '051_large_clamp',
                       '052_extra_large_clamp', '061_foam_brick']
        return class_names

    def image_reference(self, image_id):
        info = self.image_info[image_id]
        if info['source'] == 'ycb':
            return info['ycb']
        else:
            super(self.__class__).image_reference(self, image_id)

    def load_mask(self, image_id):
        image_masks = self.mask_files[str(image_id)]
        instance_masks = []
        class_ids = []
        for class_id in list(image_masks.keys()):
            bgr_mask = cv2.imread(image_masks[class_id])
            binary_mask = np.clip(bgr_mask, 0, 1)
            instance_masks.append(binary_mask[:, :, 0])
            class_ids.append(class_id)
        class_ids = np.array(class_ids, dtype=np.int32)
        mask = np.stack(instance_masks, axis=2).astype(np.bool)
        return mask, class_ids

    def load_rgb_mask(self, image_id):
        image_masks = self.mask_files[str(image_id)]
        instance_masks = []
        class_ids = []
        for class_id in list(image_masks.keys()):
            bgr_mask = cv2.imread(image_masks[class_id])
            rgb_mask = cv2.cvtColor(bgr_mask, cv2.COLOR_BGR2RGB)
            instance_masks.append(rgb_mask)
            class_ids.append(class_id)
        class_ids = np.array(class_ids, dtype=np.int32)
        mask = np.stack(instance_masks, axis=2)
        mask = np.reshape(mask, (320, 320, 3*len(class_ids)))
        mask = mask / 255.0
        return mask, class_ids

ext: py

max_stars_repo_name: SushmaDG/MaskRCNN

max_issues_repo_name: SushmaDG/MaskRCNN

max_forks_repo_name: SushmaDG/MaskRCNN

reason: import numpy


text:
import argparse
import tensorflow as tf
import numpy as np

def main(params):
    cell_numbers = params.cells_number
    fashion_mnist = tf.keras.datasets.fashion_mnist
    (train_images, train_labels), (test_images, test_labels) = fashion_mnist.load_data()
    train_images = train_images / 255.0
    test_images = test_images / 255.0
    model = tf.keras.Sequential([
    tf.keras.layers.Flatten(input_shape=(28, 28)),
    tf.keras.layers.Dense(128, activation='relu'),
    tf.keras.layers.Dense(10)])

    model.compile(optimizer='adam',
             loss=tf.keras.losses.SparseCategoricalCrossentropy(from_logits=True),
             metrics=['accuracy'])
    model.fit(train_images, train_labels, epochs=10)
    test_acc = model.evaluate(test_images,  test_labels, verbose=2)
    print('accuracy:', test_acc)
    model.save('my_model')

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description='katib experiment')
    parser.add_argument('--cells_number', type=int, default=128)
    params = parser.parse_args()
    main(params)

max_stars_repo_name: datastacktv/kubeflow-introduction


text:
import sys
import numpy as np
import os
import time
import math
from PIL import Image
import cv2
from datetime import datetime
from pynq import Xlnk
from pynq import Overlay
import pynq
import struct
from multiprocessing import Process, Pipe, Queue, Event, Manager
from IoU import Average_IoU

anchor = [1.4940052559648322, 2.3598481287086823, 4.0113013115312155, 5.760873975661669]
bbox_m = [52., 48., 28., 30., 124., 47., 52., 23., 23., 125.]
qm = 131072.0
w = 40
h = 20
IMG_DIR = '../sample1000/'
BATCH_SIZE = 4

def sigmoid(x):
    return 1/(1+np.exp(-x))

def get_image_names():
    names_temp = [f for f in os.listdir(IMG_DIR) if f.endswith('.jpg')]
    names_temp.sort(key= lambda x:int(x[:-4]))
    return names_temp

def get_image_batch():
    image_list = get_image_names()
    batches = list()
    for i in range(0, len(image_list), BATCH_SIZE):
        batches.append(image_list[i:i+BATCH_SIZE])
    return batches

IMAGE_NAMES  = get_image_batch()
IAMGE_NAMES_F = np.array(IMAGE_NAMES).reshape(-1).tolist()
IMAGE_NAMES_LEN = len(IMAGE_NAMES)
IMAGE_BUFF  = []

def stitch(batch):
    image   = np.zeros((4,160,320,4),np.uint8)
    image[0] = np.array(Image.open(IMG_DIR+batch[0]).convert('RGBA').resize((320, 160)))
    image[1] = np.array(Image.open(IMG_DIR+batch[1]).convert('RGBA').resize((320, 160)))
    image[2] = np.array(Image.open(IMG_DIR+batch[2]).convert('RGBA').resize((320, 160)))
    image[3] = np.array(Image.open(IMG_DIR+batch[3]).convert('RGBA').resize((320, 160)))
    return image

def compute_bounding_box(bbox, bbox_origin, batch, result):
    for b in range(4):
        if(bbox_origin[b,4]>0):
            xs = bbox_origin[b][0]*bbox_m[5]/qm
            ys = bbox_origin[b][1]*bbox_m[6]/qm
            ws = bbox_origin[b][2]*bbox_m[7]/qm
            hs = bbox_origin[b][3]*bbox_m[8]/qm
            ws_inb = math.exp(ws)*anchor[2]
            hs_inb = math.exp(hs)*anchor[3]
        else:
            xs = bbox_origin[b][0]*bbox_m[0]/qm
            ys = bbox_origin[b][1]*bbox_m[1]/qm
            ws = bbox_origin[b][2]*bbox_m[2]/qm
            hs = bbox_origin[b][3]*bbox_m[3]/qm
            ws_inb = np.exp(ws)*anchor[0]
            hs_inb = np.exp(hs)*anchor[1]
        xs_inb = sigmoid(xs) + bbox_origin[b][5]
        ys_inb = sigmoid(ys) + bbox_origin[b][6]
        bcx = xs_inb/w
        bcy = ys_inb/h
        bw = ws_inb/w
        bh = hs_inb/h
        bbox[b][0] = bcx - bw/2.0
        bbox[b][1] = bcy - bh/2.0
        bbox[b][2] = bcx + bw/2.0
        bbox[b][3] = bcy + bh/2.0

        x1 = int(round(bbox[b][0] * 640))
        y1 = int(round(bbox[b][1] * 360))
        x2 = int(round(bbox[b][2] * 640))
        y2 = int(round(bbox[b][3] * 360))
        x1 = np.clip(x1,1,640)
        y1 = np.clip(y1,1,360)
        x2 = np.clip(x2,1,640)
        y2 = np.clip(y2,1,360)

        result.write(batch[b].split('.')[0].zfill(3) + '.jpg' +' '+str([x1, x2, y1, y2])+'\n')
        print(batch[b], [x1, x2, y1, y2])

################################## Init FPGA ##################################
xlnk = Xlnk()
xlnk.xlnk_reset()

img   = xlnk.cma_array(shape=[4,160,320,4], dtype=np.uint8)
fm= xlnk.cma_array(shape=(628115*32), dtype=np.uint8)
weight = xlnk.cma_array(shape=(220672),  dtype=np.int16)
biasm  = xlnk.cma_array(shape=(432*16),  dtype=np.int16)
print("Allocating memory done")

parameter = np.fromfile("./weight/SkyNet.bin", dtype=np.int16)
np.copyto(weight, parameter[0:220672])
np.copyto(biasm[0:428*16], parameter[220672:])
print("Parameters loading done")

overlay = Overlay("SkyNet.bit")
print("Bitstream loaded")

SkyNet = overlay.SkyNet
SkyNet.write(0x10, img.physical_address)
SkyNet.write(0x1c, fm.physical_address)
SkyNet.write(0x28, weight.physical_address)
SkyNet.write(0x34, biasm.physical_address)

rails = pynq.get_rails()
recorder = pynq.DataRecorder(rails["power1"].power)

################################## Main process ##################################
bbox_origin = np.empty(64, dtype=np.int16)
bbox  = np.zeros((4,4),dtype=np.float32)
result= open('predict.txt','w+')
batch_buff  = None

image_buff  = np.zeros((4,160,320,4),np.uint8)

print("\n**** Start to detect")
start = time.time()
with recorder.record(0.05):
    for idx, batch in enumerate(IMAGE_NAMES):

        if (idx == 0):
            np.copyto(img, stitch(IMAGE_NAMES[idx]))
        else:
            np.copyto(img, image_buff)

        acc_start = time.time()

        SkyNet.write(0x00, 1)
        if (0 <= idx < IMAGE_NAMES_LEN - 1):
            # pass
            image_buff = stitch(IMAGE_NAMES[idx+1])
        if (0 < idx <= IMAGE_NAMES_LEN - 1):
            compute_bounding_box(bbox ,bbox_origin.reshape(4,16), batch_buff, result)

        isready = SkyNet.read(0x00)
        while( isready == 1 ):
            isready = SkyNet.read(0x00)

        np.copyto(bbox_origin, biasm[428*16:])
        batch_buff = batch

        if (idx == IMAGE_NAMES_LEN - 1):
            compute_bounding_box(bbox ,bbox_origin.reshape(4,16), batch, result)

        acc_end = time.time()
        total_time = acc_end - acc_start

result.close()

end = time.time()
total_time = end - start
print("**** Detection finished\n")
energy = recorder.frame["power1_power"].mean() * total_time
print(str(energy)+'J')

IoU = Average_IoU(IMG_DIR+'ground_truth.txt', 'predict.txt')
print('Total time: ' + str(total_time) + ' s')


text:
""" Code for the MAML algorithm and network definitions. """
from __future__ import print_function
import numpy as np
import sys
import tensorflow as tf
try:
    import special_grads
except KeyError as e:
    print('WARN: Cannot define MaxPoolGrad, likely already defined for this version of tensorflow: %s' % e,
          file=sys.stderr)

from tensorflow.python.platform import flags
from utils import mse, xent, conv_block, normalize

FLAGS = flags.FLAGS

class R2D2:
    def __init__(self, dim_input=1, dim_output=1, test_num_updates=5):
        #tf.reset_default_graph()

        """ must call construct_model() after initializing MAML! """
        self.dim_input = dim_input
        self.dim_output = dim_output
        self.update_lr = FLAGS.update_lr
        self.meta_lr = tf.placeholder_with_default(FLAGS.meta_lr, ())
        self.classification = False
        self.test_num_updates = test_num_updates
        if FLAGS.datasource == 'sinusoid':
            self.dim_hidden = [40, 40]
            self.loss_func = mse
            self.forward = self.forward_fc
            self.construct_weights = self.construct_fc_weights
        elif FLAGS.datasource == 'omniglot' or FLAGS.datasource == 'miniimagenet' or FLAGS.datasource == 'cifarfs':
            self.loss_func = xent
            self.classification = True

            if FLAGS.conv:
                self.dim_hidden = FLAGS.num_filters
                self.forward = self.forward_conv
                self.construct_weights = self.construct_conv_weights
            else:
                self.dim_hidden = [256, 128, 64, 64]
                self.forward=self.forward_fc
                self.construct_weights = self.construct_fc_weights

            # Determine amount of channels to use
            if FLAGS.datasource == 'miniimagenet' or FLAGS.datasource == 'cifarfs':
                self.channels = 3
            else:
                self.channels = 1

            # Compute image width (=height)
            self.img_size = int(np.sqrt(self.dim_input/self.channels)) # dim input is length of totally flattened image
        else:
            raise ValueError('Unrecognized data source.')

    def construct_model(self, input_tensors=None, prefix='metatrain_'):
        # This function constructs the model, and defines the ops. The ops are not called yet! That happens in session.run(...)

        # a: training data for inner gradient, b: test data for meta gradient
        if input_tensors is None:
            self.inputa = tf.placeholder(tf.float32)
            self.inputb = tf.placeholder(tf.float32)
            self.labela = tf.placeholder(tf.float32)
            self.labelb = tf.placeholder(tf.float32)
        else: # Directly couple input tensors from tf queue to object variables
            self.inputa = input_tensors['inputa']
            self.inputb = input_tensors['inputb']
            self.labela = input_tensors['labela']
            self.labelb = input_tensors['labelb']

        with tf.variable_scope('model', reuse=None) as training_scope:
            if 'weights' in dir(self):
                # weights were already initialized during some training, reuse those
                training_scope.reuse_variables()
                weights = self.weights
            else:
                # Define the weights
                # this is done when construct_model is called
                self.weights = weights = self.construct_weights()

            # outputbs[i] and lossesb[i] is the output and loss after i+1 gradient updates
            lossesa, outputas, lossesb, outputbs, labelas, labelbs = [], [], [], [], [], []
            accuraciesa, accuraciesb = [], []
            num_updates = max(self.test_num_updates, FLAGS.num_updates)
            outputbs = [[]]*num_updates
            lossesb = [[]]*num_updates
            accuraciesb = [[]]*num_updates

            def task_metalearn(inp, reuse=True):
                """ Perform gradient descent for one task in the meta-batch. """
                inputa, inputb, labela, labelb = inp
                task_outputbs, task_lossesb = [], []

                if self.classification:
                    task_accuraciesb = []

                task_outputa = self.forward(inputa, weights, reuse=reuse)  # only reuse on the first iter
                task_lossa = self.loss_func(task_outputa, labela)

                ## Pass through CNN
                x = self.forward_conv_CNN(inputa, weights, reuse=True)

                ## Linear Regression with Woodbury Identity
                # using training set (a) to determine new weights for linear regressor
                xT = tf.transpose(x)
                xxT = tf.matmul(x,xT)

                # Calculate new LINEAR REGRESSION weights on train set, using the Woodbury identity
                fast_weights = dict(zip(weights.keys(), [weights[key] for key in weights.keys()]))
                fast_weights['stop_w5'] = tf.stop_gradient(tf.matmul(tf.matmul(xT,tf.linalg.inv(xxT + weights['lr_lambda'] * tf.eye(tf.shape(xxT)[0],tf.shape(xxT)[1]))),labela))

                # MAML line 8: calculate output/loss on test set (b), internally does LR conversion with scale alpha and bias beta
                output = self.forward(inputb, fast_weights, reuse=True)
                task_outputbs.append(output)
                task_lossesb.append(self.loss_func(output, labelb))

                # NO (further) INNER STEPS REQUIRED FOR LINEAR REGRESSION
                """
                for j in range(num_updates - 1):
                    loss = self.loss_func(self.forward(inputa, fast_weights, reuse=True), labela)

                    # MAML line 5: evaluate grads on train set (a)
                    grads = tf.gradients(loss, list(fast_weights.values()))
                    if FLAGS.stop_grad:
                        grads = [tf.stop_gradient(grad) for grad in grads]
                    gradients = dict(zip(fast_weights.keys(), grads))

                    # MAML line 6: compute updates (adapted parameters)
                    fast_weights = dict(zip(fast_weights.keys(), [fast_weights[key] - self.update_lr*gradients[key] for key in fast_weights.keys()]))

                    # MAML line 8: calculate output/loss on test set (b)
                    output = self.forward(inputb, fast_weights, reuse=True)
                    task_outputbs.append(output)
                    task_lossesb.append(self.loss_func(output, labelb))
                """
                task_output = [task_outputa, task_outputbs, task_lossa, task_lossesb]

                if self.classification:
                    task_accuracya = tf.contrib.metrics.accuracy(tf.argmax(tf.nn.softmax(task_outputa), 1), tf.argmax(labela, 1))
                    for j in range(num_updates):
                        task_accuraciesb.append(tf.contrib.metrics.accuracy(tf.argmax(tf.nn.softmax(task_outputbs[j]), 1), tf.argmax(labelb, 1)))
                    task_output.extend([task_accuracya, task_accuraciesb])

                return task_output

            if FLAGS.norm is not 'None': # to initialize batch norm variables
                # to initialize the batch norm vars, might want to combine this, and not run idx 0 twice.
                unused = task_metalearn((self.inputa[0], self.inputb[0], self.labela[0], self.labelb[0]), False)

            out_dtype = [tf.float32, [tf.float32]*num_updates, tf.float32, [tf.float32]*num_updates]
            if self.classification: # accuracies are also stored
                out_dtype.extend([tf.float32, [tf.float32]*num_updates])
            # THE REAL LEARNING CONSTRUCTION OCCURS HERE
            # IMPORTANT: executes in parallel for ALL TASKS in batch I guess? The inputs are formatted in a special way to contain multiple tasks?
            result = tf.map_fn(task_metalearn, elems=(self.inputa, self.inputb, self.labela, self.labelb), dtype=out_dtype, parallel_iterations=FLAGS.meta_batch_size)
            if self.classification:
                outputas, outputbs, lossesa, lossesb, accuraciesa, accuraciesb = result
            else:
                outputas, outputbs, lossesa, lossesb  = result

        ## Performance & Optimization
        if 'train' in prefix:
            self.total_loss1 = total_loss1 = tf.reduce_sum(lossesa) / tf.to_float(FLAGS.meta_batch_size)
            self.total_losses2 = total_losses2 = [tf.reduce_sum(lossesb[j]) / tf.to_float(FLAGS.meta_batch_size) for j in range(num_updates)]
            # after the map_fn
            self.outputas, self.outputbs = outputas, outputbs
            #self.labelas, self.labelbs = outputas, outputbs
            if self.classification:
                self.total_accuracy1 = total_accuracy1 = tf.reduce_sum(accuraciesa) / tf.to_float(FLAGS.meta_batch_size)
                self.total_accuracies2 = total_accuracies2 = [tf.reduce_sum(accuraciesb[j]) / tf.to_float(FLAGS.meta_batch_size) for j in range(num_updates)]
            self.pretrain_op = tf.train.AdamOptimizer(self.meta_lr).minimize(total_loss1)

            if FLAGS.metatrain_iterations > 0: # FLAGS.metatrain_iterations = how many times to execute
                # This is the meta optimizer
                optimizer = tf.train.AdamOptimizer(self.meta_lr)

                # Compute gradients after num_updates
                self.gvs = gvs = optimizer.compute_gradients(self.total_losses2[FLAGS.num_updates-1])

                #grads = tf.gradients(loss, list(fast_weights.values()))
                #grads = [tf.stop_gradient(grad) for grad in gvs]

                # Gradients are clipped by [-10,10] to avoid explosion?
                if FLAGS.datasource == 'miniimagenet' or FLAGS.datasource == 'cifarfs':
                    gvs = [(tf.clip_by_value(grad, -10, 10), var) for grad, var in gvs if grad is not None]

                # update parameters
                self.metatrain_op = optimizer.apply_gradients(gvs)
        else:
            self.metaval_total_loss1 = total_loss1 = tf.reduce_sum(lossesa) / tf.to_float(FLAGS.meta_batch_size)
            self.metaval_total_losses2 = total_losses2 = [tf.reduce_sum(lossesb[j]) / tf.to_float(FLAGS.meta_batch_size) for j in range(num_updates)]
            if self.classification:
                self.metaval_total_accuracy1 = total_accuracy1 = tf.reduce_sum(accuraciesa) / tf.to_float(FLAGS.meta_batch_size)
                self.metaval_total_accuracies2 = total_accuracies2 =[tf.reduce_sum(accuraciesb[j]) / tf.to_float(FLAGS.meta_batch_size) for j in range(num_updates)]

                # For diagnostic purposes
                self.test_accuraciesa = accuraciesa
                self.test_accuraciesb = [accuraciesb[j] for j in range(num_updates)]

                self.test_outputas = outputas
                self.test_outputbs = outputbs

                self.labelas = self.labela
                self.labelbs = self.labelb

        ## Summaries
        tf.summary.scalar(prefix+'Pre-update loss', total_loss1)
        if self.classification:
            tf.summary.scalar(prefix+'Pre-update accuracy', total_accuracy1)

        for j in range(num_updates):
            tf.summary.scalar(prefix+'Post-update loss, step ' + str(j+1), total_losses2[j])
            if self.classification:
                tf.summary.scalar(prefix+'Post-update accuracy, step ' + str(j+1), total_accuracies2[j])

    ### Network construction functions
    ## CNN
    # initialize and return weights for CNN
    def construct_conv_weights(self):
        weights = {}

        dtype = tf.float32
        conv_initializer =  tf.contrib.layers.xavier_initializer_conv2d(dtype=dtype)
        fc_initializer =  tf.contrib.layers.xavier_initializer(dtype=dtype)
        k = 3

        # CNN weights
        weights['conv1'] = tf.get_variable('conv1', [k, k, self.channels, self.dim_hidden], initializer=conv_initializer, dtype=dtype)
        weights['b1'] = tf.Variable(tf.zeros([self.dim_hidden]))
        weights['conv2'] = tf.get_variable('conv2', [k, k, self.dim_hidden, self.dim_hidden], initializer=conv_initializer, dtype=dtype)
        weights['b2'] = tf.Variable(tf.zeros([self.dim_hidden]))
        weights['conv3'] = tf.get_variable('conv3', [k, k, self.dim_hidden, self.dim_hidden], initializer=conv_initializer, dtype=dtype)
        weights['b3'] = tf.Variable(tf.zeros([self.dim_hidden]))
        weights['conv4'] = tf.get_variable('conv4', [k, k, self.dim_hidden, self.dim_hidden], initializer=conv_initializer, dtype=dtype)
        weights['b4'] = tf.Variable(tf.zeros([self.dim_hidden]))

        # RR weights
        # assumes max pooling, flat_dim is concatenated flattened output of layer 3 and 4
        flat_dim = 640
        if FLAGS.datasource == 'miniimagenet': # 84x84 * (1/2 + 1/2/2/2)
            flat_dim = 4000
        else:# cifarfs 32x32 * (1/2 + 1/2/2/2) = 640
            flat_dim = 640

        weights['stop_w5'] = tf.get_variable('stop_w5', [flat_dim, self.dim_output], initializer=fc_initializer)

        # hyper parameters of base learner, to be learnt in outer loop together with CNN parameters
        #weights['lr_lambda'] = tf.get_variable('lr_lambda', initializer=tf.constant(1., dtype=dtype), dtype=dtype)
        #weights['lr_alpha'] = tf.get_variable('lr_alpha',initializer=tf.constant(1., dtype=dtype), dtype=dtype)
        #weights['lr_beta'] = tf.get_variable('lr_beta', initializer=tf.constant(1., dtype=dtype), dtype=dtype)

        weights['lr_lambda'] = tf.Variable(tf.zeros(1, dtype = dtype))
        weights['lr_alpha'] = tf.Variable(tf.zeros(1, dtype = dtype))
        weights['lr_beta'] = tf.Variable(tf.zeros(1, dtype = dtype))

        #weights['b5'] = tf.Variable(tf.zeros([self.dim_output]), name='b5')

        return weights

    # return output of input image, with weights given as argument!
    # This is only to be used in the meta-learning step, during base training the direct solution for LR is used!
    def forward_conv(self, inp, weights, reuse=False, scope=''):
        out = self.forward_conv_CNN(inp, weights, reuse=reuse, scope=scope)
        out = self.forward_conv_lr(out, weights, reuse=reuse, scope=scope)

        return out

    def forward_conv_CNN(self, inp, weights, reuse=False, scope=''):
        # reuse is for the normalization parameters.
        channels = self.channels
        inp = tf.reshape(inp, [-1, self.img_size, self.img_size, channels])

        hidden1 = conv_block(inp, weights['conv1'], weights['b1'], reuse, scope+'0')
        hidden2 = conv_block(hidden1, weights['conv2'], weights['b2'], reuse, scope+'1')
        hidden3 = conv_block(hidden2, weights['conv3'], weights['b3'], reuse, scope+'2')
        hidden4 = conv_block(hidden3, weights['conv4'], weights['b4'], reuse, scope+'3')

        # Flattening of blocks 3 and 4
        hidden3 = tf.reshape(hidden3, [-1, np.prod([int(dim) for dim in hidden3.get_shape()[1:]])])
        hidden4 = tf.reshape(hidden4, [-1, np.prod([int(dim) for dim in hidden4.get_shape()[1:]])])

        # Concatenate
        flatconcat34 = tf.concat([hidden3, hidden4], axis=1) # keep batched (axis 0), concatenate columns (axis 1)

        return flatconcat34

    def forward_conv_lr(self, inp, weights, reuse=False, scope=''):
        W = tf.stop_gradient(weights['stop_w5']) # stop backpropagation to CNN weights through RR calculation
        return tf.multiply(weights['lr_alpha'],tf.matmul(inp, W)) + tf.multiply(weights['lr_beta'],tf.ones(shape=[inp.get_shape()[0], self.dim_output], dtype=tf.float32))


text:
#!/usr/bin/python
# -*- coding: utf-8 -*-
##
# rb.py: Models for accelerated randomized benchmarking.
##
# © 2017, Chris Ferrie (csferrie@gmail.com) and
#         Christopher Granade (cgranade@cgranade.com).
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
#     1. Redistributions of source code must retain the above copyright
#        notice, this list of conditions and the following disclaimer.
#
#     2. Redistributions in binary form must reproduce the above copyright
#        notice, this list of conditions and the following disclaimer in the
#        documentation and/or other materials provided with the distribution.
#
#     3. Neither the name of the copyright holder nor the names of its
#        contributors may be used to endorse or promote products derived from
#        this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
##

## FEATURES ###################################################################

from __future__ import absolute_import
from __future__ import division

## ALL ########################################################################

__all__ = [
    'RandomizedBenchmarkingModel'
]

## IMPORTS ####################################################################

from itertools import starmap

import numpy as np
from qinfer._due import due, Doi
from qinfer.abstract_model import FiniteOutcomeModel, DifferentiableModel

from operator import mul

## FUNCTIONS ##################################################################

def p(F, d=2):
    """
    Given the fidelity of a gate in :math:`d` dimensions, returns the
    depolarizating probability of the twirled channel.

    :param float F: Fidelity of a gate.
    :param int d: Dimensionality of the Hilbert space on which the gate acts.
    """
    return (d * F - 1) / (d - 1)

def F(p, d=2):
    """
    Given the depolarizating probabilty of a twirled channel in :math:`d`
    dimensions, returns the fidelity of the original gate.

    :param float p: Depolarizing parameter for the twirled channel.
    :param int d: Dimensionality of the Hilbert space on which the gate acts.
    """
    return 1 - (1 - p) * (d - 1) / d

## CLASSES ####################################################################

class RandomizedBenchmarkingModel(FiniteOutcomeModel, DifferentiableModel):
    r"""
    Implements the randomized benchmarking or interleaved randomized
    benchmarking protocol, such that the depolarizing strength :math:`p`
    of the twirled channel is a parameter to be estimated, given a sequnce
    length :math:`m` as an experimental control. In addition, the zeroth-order
    "fitting"-parameters :math:`A` and :math:`B` are represented as model
    parameters to be estimated.

    :param bool interleaved: If `True`, the model implements the interleaved
        protocol, with :math:`\tilde{p}` being the depolarizing parameter for
        the interleaved gate and with :math:`p_{\text{ref}}` being the reference
        parameter.

    :modelparam p: Fidelity of the twirled error channel :math:`\Lambda`, represented as
        a decay rate :math:`p = (d F - 1) / (d - 1)`, where :math:`F`
        is the fidelity and :math:`d` is the dimension of the Hilbert space.
    :modelparam A: Scale of the randomized benchmarking decay, defined as
        :math:`\Tr[Q \Lambda(\rho - \ident / d)]`, where :math:`Q` is the final
        measurement, and where :math:`\ident` is the initial preparation.
    :modelparam B: Offset of the randomized benchmarking decay, defined as
        :math:`\Tr[Q \Lambda(\ident / d)]`.

    :expparam int m: Length of the randomized benchmarking sequence
        that was measured.
    """
    # TODO: add citations to the above docstring.

    @due.dcite(
        Doi("10.1088/1367-2630/17/1/013042"),
        description="Accelerated randomized benchmarking",
        tags=["implementation"]
    )
    def __init__(self, interleaved=False, order=0):
        self._il = interleaved
        if order != 0:
            raise NotImplementedError(
                "Only zeroth-order is currently implemented."
            )
        super(RandomizedBenchmarkingModel, self).__init__()

    @property
    def n_modelparams(self):
        return 3 + (1 if self._il else 0)

    @property
    def modelparam_names(self):
        return (
            # We want to know \tilde{p} := p_C / p, and so we make it
            # a model parameter directly. This means that later, we'll
            # need to extract p_C = p \tilde{p}.
            [r'\tilde{p}', 'p', 'A', 'B']
            if self._il else
            ['p', 'A', 'B']
        )

    @property
    def is_n_outcomes_constant(self):
        return True
    @property
    def expparams_dtype(self):
        return [('m', 'uint')] + (
            [('reference', bool)] if self._il else []
        )

    def n_outcomes(self, expparams):
        return 2

    def are_models_valid(self, modelparams):
        if self._il:
            p_C, p, A, B = modelparams.T
            return np.all([
                0 <= p,
                p <= 1,
                0 <= p_C,
                p_C <= 1,
                0 <= A,
                A <= 1,
                0 <= B,
                B <= 1,
                A + B <= 1,
                A * p + B <= 1,
                A * p_C + B <= 1
            ], axis=0)
        else:
            p, A, B = modelparams.T
            return np.all([
                0 <= p,
                p <= 1,
                0 <= A,
                A <= 1,
                0 <= B,
                B <= 1,
                A + B <= 1,
                A * p + B <= 1
            ], axis=0)

    def likelihood(self, outcomes, modelparams, expparams):
        super(RandomizedBenchmarkingModel, self).likelihood(outcomes, modelparams, expparams)

        if self._il:
            p_tilde, p, A, B = modelparams.T[:, :, np.newaxis]

            p_C = p_tilde * p

            p = np.where(expparams['reference'][np.newaxis, :], p, p_C)
        else:
            p, A, B = modelparams.T[:, :, np.newaxis]

        m = expparams['m'][np.newaxis, :]

        pr0 = np.zeros((modelparams.shape[0], expparams.shape[0]))
        pr0[:, :] = 1 - (A * (p ** m) + B)

        return FiniteOutcomeModel.pr0_to_likelihood_array(outcomes, pr0)

    def score(self, outcomes, modelparams, expparams, return_L=False):

        na = np.newaxis
        n_m = modelparams.shape[0]
        n_e = expparams.shape[0]
        n_o = outcomes.shape[0]
        n_p = self.n_modelparams

        m = expparams['m'].reshape((1, 1, 1, n_e))

        L = self.likelihood(outcomes, modelparams, expparams)[na, ...]
        outcomes = outcomes.reshape((1, n_o, 1, 1))

        if not self._il:

            p, A, B = modelparams.T[:, :, np.newaxis]
            p = p.reshape((1, 1, n_m, 1))
            A = A.reshape((1, 1, n_m, 1))
            B = B.reshape((1, 1, n_m, 1))

            q = (-1)**(1-outcomes) * np.concatenate(np.broadcast_arrays(
                A * m * (p ** (m-1)), p**m, np.ones_like(p),
            ), axis=0) / L

        else:

            p_tilde, p_ref, A, B = modelparams.T[:, :, np.newaxis]
            p_C = p_tilde * p_ref

            mode = expparams['reference'][np.newaxis, :]

            p = np.where(mode, p_ref, p_C)

            p = p.reshape((1, 1, n_m, n_e))
            A = A.reshape((1, 1, n_m, 1))
            B = B.reshape((1, 1, n_m, 1))

            q = (-1)**(1-outcomes) * np.concatenate(np.broadcast_arrays(
                np.where(mode, 0, A * m * (p_tilde ** (m - 1)) * (p_ref ** m)),
                np.where(mode,
                    A * m * (p_ref ** (m - 1)),
                    A * m * (p_ref ** (m - 1)) * (p_tilde ** m)
                ),
                p**m, np.ones_like(p)
            ), axis=0) / L

        if return_L:
            # Need to strip off the extra axis we added for broadcasting to q.
            return q, L[0, ...]
        else:
            return q

max_stars_repo_name: mikedeltalima/python-qinfer

max_issues_repo_name: mikedeltalima/python-qinfer

max_forks_repo_name: mikedeltalima/python-qinfer


text:
# -*- coding: utf-8 -*-
"""Kelvin-Helmholtz test problem
Initial conditions taken from
"A Well Posed Kelvin-Helmholtz Instability Test Problem", C.P. McNally, et. al
See https://arxiv.org/pdf/1111.1764.pdf or
https://iopscience.iop.org/article/10.1088/0067-0049/201/2/18/meta"""
try:
    import matplotlib.pyplot as plt
except ImportError:
    pass

from configparser import ConfigParser
import numpy as np
import sys
import os

sys.path.append(os.path.abspath("../../.."))
from pycato import *

# Make the empty grid
domain = make_uniform_grid(
    n_cells=(200, 200), xrange=(0.0, 1.0), yrange=(0.0, 1.0), input_file="input.ini",
)

# Set the initial conditions
domain["p"] = domain["p"] * 2.5
rho = domain["rho"].m
u = domain["u"].m
v = domain["v"].m
rho_1 = 1.0
rho_2 = 2.0
L = 0.025
U_1 = 0.5
U_2 = -0.5
U_m = (U_1 - U_2) / 2.0
rho_m = (rho_1 - rho_2) / 2.0

for i in range(domain["xc"].shape[0]):
    for j in range(domain["xc"].shape[1]):
        x = domain["xc"][i, j].m
        y = domain["yc"][i, j].m

        # Perturbed velocity
        v[i, j] = 0.01 * np.sin(4.0 * np.pi * x)

        if 0.25 > y >= 0.0:
            rho[i, j] = rho_1 - rho_m * np.exp((y - 0.25) / L)
            u[i, j] = U_1 - U_m * np.exp((y - 0.25) / L)
        elif 0.5 > y >= 0.25:
            rho[i, j] = rho_2 + rho_m * np.exp((-y + 0.25) / L)
            u[i, j] = U_2 + U_m * np.exp((-y + 0.25) / L)
        elif 0.75 > y >= 0.5:
            rho[i, j] = rho_2 + rho_m * np.exp(-(0.75 - y) / L)
            u[i, j] = U_2 + U_m * np.exp(-(0.75 - y) / L)
        elif 1.0 > y >= 0.75:
            rho[i, j] = rho_1 - rho_m * np.exp(-(y - 0.75) / L)
            u[i, j] = U_1 - U_m * np.exp(-(y - 0.75) / L)

write_initial_hdf5(filename="initial_conditions", initial_condition_dict=domain)

# Plot the results
try:
    fig, (ax1, ax2, ax3) = plt.subplots(figsize=(18, 8), nrows=1, ncols=3)

    dens = ax1.pcolormesh(
        domain["x"],
        domain["y"],
        domain["rho"],
        edgecolor="k",
        lw=0.001,
        antialiased=True,
    )
    fig.colorbar(dens, ax=ax1, label="Density")

    uax = ax2.pcolormesh(
        domain["x"], domain["y"], domain["u"], edgecolor="k", lw=0.001, antialiased=True
    )
    fig.colorbar(uax, ax=ax2, label="X Velocity")

    vax = ax3.pcolormesh(
        domain["x"], domain["y"], domain["v"], edgecolor="k", lw=0.001, antialiased=True
    )
    fig.colorbar(vax, ax=ax3, label="Y Velocity")

    for ax in [ax1, ax2, ax3]:
        ax.set_xlabel("X")
        ax.set_ylabel("Y")
        ax.axis("equal")

    plt.tight_layout()
    plt.show()
except Exception:
    pass

max_issues_repo_name: smillerc/cato

max_forks_repo_name: smillerc/cato


text:
#  -*- coding: utf-8 -*-
# *****************************************************************************
# NICOS, the Networked Instrument Control System of the MLZ
# Copyright (c) 2009-2022 by the NICOS contributors (see AUTHORS)
#
# This program is free software; you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free Software
# Foundation; either version 2 of the License, or (at your option) any later
# version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# this program; if not, write to the Free Software Foundation, Inc.,
# 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
#
# Module authors:
#   Michele Brambilla <michele.brambilla@psi.ch>
#
# *****************************************************************************

import kafka
import numpy

from nicos.core import Attach, status
from nicos.core.constants import LIVE
from nicos.devices.generic.detector import ImageChannelMixin, PassiveChannel

from nicos_ess.devices.kafka.consumer import KafkaSubscriber
from nicos_sinq.devices.epics.area_detector import ADKafkaPlugin

try:
    from nicos_ess.devices.fbschemas.hs00 import Array, ArrayUInt, \
        EventHistogram

except ImportError:
    EventHistogram = None
    Array = None
    ArrayUInt = None

class ADKafkaImageDetector(KafkaSubscriber, ImageChannelMixin, PassiveChannel):
    """
    Class for reading images from the stream associated with the ADPluginKafka
    """
    attached_devices = {
        'kafka_plugin': Attach('NICOS device for the ADPluginKafka',
                               ADKafkaPlugin, optional=False)
    }

    def doPreinit(self, mode):
        broker = getattr(self._attached_kafka_plugin, 'broker')
        if not broker:
            raise Exception('Can\'t find broker address in ADPluginKafka')
        self._consumer = kafka.KafkaConsumer(
            bootstrap_servers=[broker],
            auto_offset_reset='latest'  # start at latest offset
        )

        # Settings for thread to fetch new message
        self._stoprequest = True
        self._updater_thread = None
        self._lastmessage = None

    def doInit(self, mode):
        topic = getattr(self._attached_kafka_plugin, 'topic')
        if not topic:
            raise Exception('Can\'t find topic in ADPluginKafka')
        self._consumer.subscribe([topic])

    def doReadArray(self, quality=LIVE):
        deserializer = HistogramFlatbuffersDeserializer()
        try:
            data, _, _ = deserializer.decode(self._lastmessage[1])
        except Exception as e:
            self.log.error(e)
            return []
        self.readresult = []
        return data.tolist()

    def new_messages_callback(self, messages):
        if not messages:
            return
        self._lastmessage = sorted(messages, key=lambda m: m[0])[~0]

    def doStatus(self, maxage=0):
        st = self._attached_kafka_plugin.doStatus()
        if st[0] != status.OK:
            return st
        if not self._consumer:
            return status.ERROR, 'Broker failure'
        if not self._consumer.subscription():
            return status.WARN, 'No topic subscribed'
        return status.OK, ''

    def valueInfo(self):
        return ()

class HistogramFlatbuffersDeserializer:
    """
    Decode the histogram using the flatbuffers schema hs00
    """
    file_identifier = "hs00"

    @staticmethod
    def _read_as_numpy(histogram, dtype=ArrayUInt.ArrayUInt):
        result = dtype()
        result.Init(histogram.Data().Bytes, histogram.Data().Pos)
        return result.ValueAsNumpy()

    @staticmethod
    def _read_current_shape(histogram):
        result = ArrayUInt.ArrayUInt()
        result.Init(histogram.CurrentShape().Bytes,
                    histogram.CurrentShape().Pos)
        return result.ValueAsNumpy()

    def decode(self, buff):
        histogram = EventHistogram.EventHistogram().GetRootAsEventHistogram(
            buff, 0)
        data_type = histogram.DataType()
        value = []
        if data_type == Array.Array.ArrayUInt:
            value = self._read_as_numpy(histogram=histogram)
        if data_type == Array.Array.ArrayULong:
            raise NotImplementedError('This data type is not yet implemented.')
            # value = self._read_as_numpy(histogram=histogram,
            # dtype=ArrayULong.ArrayULong)
        if data_type == Array.Array.ArrayFloat:
            raise NotImplementedError('This data type is not yet implemented.')
            # value = self._read_as_numpy(histogram=histogram,
            # dtype=ArrayFloat.ArrayFloat)
        if data_type == Array.Array.ArrayDouble:
            raise NotImplementedError('This data type is not yet implemented.')
            # value = self._read_as_numpy(histogram=histogram,
            # dtype=ArrayDouble.ArrayDouble)

        value = numpy.reshape(value, histogram.CurrentShapeAsNumpy())
        return value, histogram.Timestamp(), histogram.Source().decode()


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*

#/// DEPENDENCIES
import discord                    #python3.7 -m pip install -U discord.py
import logging
import math
from numexpr import evaluate as evl
from util import embedify
from sympy.simplify.simplify import nsimplify as nsimp
from discord.ext import commands
from discord.ext.commands import Bot, MissingPermissions, has_permissions
from chk.enbl import enbl

##///---------------------///##
##///    BOT  COMMANDS    ///##
##///---------------------///##

def sqrt(n):
    return n ** 0.5

@commands.command(
    aliases = [],
    help = 'math',
    brief = 'Solves the quadratic formula',
    usage = ';]quad {A} {B} {C}',
    description = '''\
A, B, C [NUMBER] - Corrosponds to "Ax^2 + Bx + C"
'''
)
@commands.check(enbl)
async def quad(ctx, A: float, B: float, C: float):
    root = str(nsimp(f'sqrt({B**2}-4*{A}*{C})'))
    imaginary, K, D = False, 1, 1
    if '*I' in root:
        root, imaginary = root.replace('*I', ''), True
    if '*' not in root and 'sqrt(' not in root:
        K, D = float(eval(root)), 1
    elif '*' not in root and 'sqrt(' in root:
        K, D = 1, float(eval(root))
    else:
        K, D = int(root.split('*')[0]), float(eval(root))
    if imaginary:
        K = K*1j
    sol1 = evl(f"({-B}+{K}*sqrt({D}))/(2*{A})")
    sol2 = evl(f"({-B}-{K}*sqrt({D}))/(2*{A})")
    await ctx.send(embed=embedify.embedify(desc=f'''```md
#]QUADRATICS``````
{-B}+-{K}√{D}
-----------------
{2*A}``````diff
+] [{-B/(2*A)} + {K/(2*A)}√{D}] {sol1}
-] [{-B/(2*A)} - {K/(2*A)}√{D}] {sol2}```'''))

##///---------------------///##
##///     OTHER STUFF     ///##
##///---------------------///##
def setup(bot):
    print('+COM')
    bot.add_command(quad)
    print('GOOD')

def teardown(bot):
    print('-COM')
    bot.remove_command('quad')
    print('GOOD')

max_stars_repo_name: VoxelPrismatic/prizai

max_issues_repo_name: VoxelPrismatic/prizai

max_forks_repo_name: VoxelPrismatic/prizai

reason: from sympy,from numexpr


text:
# Import system services
import os
import sys

# The main package behind the graph data structures.
import networkx as nx
import random

class forest:
    """
    This class is created by an Election to track delegations that have happened. It stores delegations in a forest
    structure.
    """

    # The constructors for this class
    def __init__(self, n):
        # create empty graph with n nodes, add edges between them as we add delegations
        self.G = nx.DiGraph()
        self.G.add_nodes_from(list(range(n)))
        self.guru_tree = None   # arbitrary constructed tree in which all nodes are gurus

    def add_delegation(self, i, j):
        """
        Add a delegation from voter i to voter j.
        #TODO: This method could check for cycles and not add an edge that is not allowed
        :param i:
        :param j:
        :return: True or False depending on whether the delegation was made successfully
        """
        if self.get_guru_of_voter(j) == i:
            # i cannot delegate to someone who has i as a guru, would cause a cycle
            return False

        # remove any other delegations i is making
        i_edges = list(self.G.out_edges(i))
        self.G.remove_edges_from(i_edges)

        if i == j:
            # voter is becoming a guru, stop after removing all existing delegations
            return True

        # add new edge from i to j
        self.G.add_edge(i, j)

        # check if cycle exists
        try:
            cycle = nx.find_cycle(self.G, source=i)
        except nx.NetworkXNoCycle:
            cycle = False

        if cycle:
            # remove newly added edge
            # add back old delegations
            self.G.remove_edge(i, j)
            self.G.add_edges_from(i_edges)
            return False

        return True

    def make_random_delegation(self, delegator, voter_network, allow_guru_changes=False):
        """
        Have the given node make a random delegation. Based on the parameter, this may be a delegation to itself which
        causes the node to become (or remain) a guru. No guarantee that the result_accuracy will be different than the
        current delegation.
        :param delegator: index of the voter changing their delegation
        :param voter_network: The network the voters exist upon
        :param allow_guru_changes:
        :return:
        """
        neighbors = list(voter_network.neighbors(delegator))
        if allow_guru_changes:
            neighbors.append(delegator)
        # if allow_guru_changes and random.random() < 0.3:
        #     neighbors = [delegator]

        new_delegate = random.choice(neighbors)
        return self.add_delegation(delegator, new_delegate), new_delegate

    def get_gurus(self):
        """
        Return a list of the node indices that are gurus (have no outgoing edges). The result_accuracy has no ordering guarantee.
        :return:
        """
        sink_nodes = [node for node in self.G.nodes() if self.G.out_degree(node) == 0]
        return sink_nodes

    def get_delegators(self):
        """
        Return a list of the node indices that delegate (have an outgoing edge). The result_accuracy has no ordering guarantee.
        :return:
        """
        delegating_nodes = [node for node in self.G.nodes() if self.G.out_degree(node) == 1]
        return delegating_nodes

    def get_guru_of_voter(self, voter):
        """
        Return the index of the voter that is this voter's guru. Possibly the voter themself.
        :param voter:
        :return:
        """
        sinks = self.get_gurus()

        for sink in sinks:
            paths = []
            for vote_path in nx.all_simple_paths(self.G, voter, sink):
                paths.append(vote_path)
            if len(paths) > 0:
                break

        return paths

    def get_guru_weights(self):
        """
        Find the number of voters delegating to each guru, including itself.
        Return a dictionary mapping each guru index to their weight.
        # TODO: There may be a much faster way to do this
        :return:
        """
        # weights = []
        weights = dict()

        U = self.G.to_undirected()
        for node in self.get_gurus():
            w = nx.node_connected_component(U, node)
            weights[node] = len(w)
            # weights.append(len(w))

        return weights

    def get_guru_accuracies(self, all_accuracies):
        """
        Return a dictionary mapping the index of each guru to their accuracy.
        :param all_accuracies: List of accuracy of each voter
        :return:
        """
        return {i: all_accuracies[i] for i in self.get_gurus()}

    def get_delegation_list(self):
        """
        Return a list of length n where element i is the voter that voter i delegates to, or i if i is a guru.
        :return:
        """
        # TODO: Consider making this a parameter that is updated whenever a delegation changes
        delegations = []
        for i in range(len(self.G.nodes)):
            if self.G.out_degree(i) > 0:
                delegations.append(list(self.G.successors(i))[0])
            else:
                delegations.append(i)
        return delegations

    def create_guru_tree(self):
        """
        Create a binary tree containing all the current gurus as nodes.
        Pay no attention to the particular location of each nodes.
        :return:
        """

        gurus = self.get_gurus()
        self.guru_tree = nx.generators.classic.full_rary_tree(r=2, n=len(gurus), create_using=nx.DiGraph)

        # relabel nodes so indices match the guru indices
        mapping = {i: gurus[i] for i in range(len(gurus))}
        nx.relabel_nodes(self.guru_tree, mapping, copy=False)

    def accuracy_of_forest(self, all_accuracies):
        """
        Given accuracy of all voters in the entire network, determine the accuracy of system with the current
        delegations.
        :param all_accuracies:
        :return:
        """
        self.create_guru_tree()

        root = None
        for node in self.guru_tree:
            if len(list(self.guru_tree.predecessors(node))) == 0:
                root = node
                break
        guru_accuracies = self.get_guru_accuracies(all_accuracies)

        guru_weights = self.get_guru_weights()

        accs = self.subtree_accuracy(left_accs=None,
                                    right_accs=None,
                                    curr_weight=guru_weights[root],
                                    curr_idx=root,
                                    guru_weights=guru_weights,
                                    guru_accuracies=guru_accuracies)

        majority_probability = sum(accs[len(accs) - len(accs) // 2:])
        return majority_probability

    def subtree_accuracy(self, left_accs, right_accs, curr_weight, curr_idx, guru_weights, guru_accuracies):
        """
        Calculate the chance that this guru, along with the gurus underneath it on the guru tree, will result_accuracy in exactly
        i correct votes.
        In left_accs, right_accs, and the returned list, index i is the probability that the guru creating that list
        will result_accuracy in exactly i correct votes. Sum for i > n/2 to get the probability of a majority.
        :param left_accs: list of accuracies of left child; [0] if no left child; None if not yet set
        :param right_accs: list of accuracies of right child or empty list if no right child
        :param curr_weight: weight of the guru currently adding to this
        :param curr_idx: index/label of the guru currently calculating new_accuracies
        :param guru_weights: dictionary mapping each guru index to it's weight
        :param guru_accuracies: dictionary mapping each guru index to it's accuracy
        :return:
        """

        # base cases
        curr_successors = list(self.guru_tree.successors(curr_idx))
        if left_accs is None and len(curr_successors) > 0:
            # check if left subtree exists
                # lowest-valued (or only) child/successor of the current node
            ls = min(curr_successors)
            left_accs = self.subtree_accuracy(left_accs=None,
                                  right_accs=None,
                                  curr_weight=guru_weights[ls],
                                  curr_idx=ls,
                                  guru_weights=guru_weights,
                                  guru_accuracies=guru_accuracies)
        elif left_accs is None:
            left_accs = [1]     # base case: if no children, they have probability 1 of exactly 0 correct votes
        if right_accs is None and len(curr_successors) == 2:
            # check if right subtree exists
                # highest-valued child/successor of the current node, only if two children exist
            rs = max(curr_successors)
            right_accs = self.subtree_accuracy(left_accs=None,
                                  right_accs=None,
                                  curr_weight=guru_weights[rs],
                                  curr_idx=rs,
                                  guru_weights=guru_weights,
                                  guru_accuracies=guru_accuracies)
        elif right_accs is None:
            right_accs = [1]

        max_child_weight = max(len(left_accs)-1, 0) + max(len(right_accs)-1, 0)
        child_accuracies = [0 for _ in range(max_child_weight+curr_weight+1)]   # weight of both subtrees and this plus 1 for 0 weight

        # merge the accuracy of both child subtrees
        # TODO: Doesn't account for one of the subtrees being empty (can we make it [0]?)
        for l_idx, l_acc in enumerate(left_accs):
            for r_idx, r_acc in enumerate(right_accs):

                child_accuracies[l_idx + r_idx] += l_acc * r_acc

        # add this guru's accuracy to that of the combined children's accuracies
        curr_acc = guru_accuracies[curr_idx]
        new_accuracies = child_accuracies.copy()
        for w in range(0, len(child_accuracies)):

            if w >= curr_weight:
                # chance of being incorrect and children getting w correct + chance of children getting w-curr_weight correct
                new_accuracies[w] = child_accuracies[w]*(1-curr_acc) + child_accuracies[w-curr_weight]*curr_acc
            else:
                # weights below curr_weight can only be reached if this guru is incorrect
                new_accuracies[w] = child_accuracies[w] * (1 - curr_acc)

        return new_accuracies

def delegation_forest_from_delegations(delegations):
    """
    Construct a delegation forest from the given list of delegations
    :param delegations: List where value at index i indicates where voter i delegates
    :return:
    """
    f = forest(n=len(delegations))
    for delegator, recipient in enumerate(delegations):
        f.add_delegation(delegator, recipient)

    return f

reason: import networkx


text:
import os.path
from amuse.test.amusetest import get_path_to_results
try:
    from matplotlib import pyplot
    HAS_MATPLOTLIB = True
    from amuse.plot import plot, semilogy, xlabel, ylabel, loglog
except ImportError:
    HAS_MATPLOTLIB = False

from amuse.units import units
from amuse.units import constants
from amuse.units.generic_unit_converter import ConvertBetweenGenericAndSiUnits
from amuse.support.exceptions import AmuseException
from amuse.community.mesa.interface import MESA
from amuse.community.gadget2.interface import Gadget2
from amuse.community.fi.interface import Fi
from amuse.ext.star_to_sph import convert_stellar_model_to_SPH

import numpy

from amuse.datamodel import Particles
from amuse.datamodel import ParticlesSuperset
from amuse.datamodel import Grid
def head_on_stellar_merger(
        masses = [0.3, 3.0] | units.MSun,
        star_age = 310.0 | units.Myr,
        initial_separation = 4.0 | units.RSun,
        angle = numpy.pi / 3,
        initial_speed = 3000.0 | units.km / units.s,
        initial_speed_perpendicular = 30.0 | units.km / units.s,
        number_of_sph_particles = 50000,
        t_end = 1.0e4 | units.s,
        sph_code = Fi,
    ):
    """
    masses: Mass of the two stars
    star_age: Initial age of the stars
    number_of_sph_particles: Total number of particles of both stars, divided according to their masses
    t_end: (Physical, not computational) duration of the hydrodynamics simulation
    sph_code: Code to use for the hydrodynamics simulation
    """

    # Convert some of the input parameters to string, for use in output file names:
    n_string = "n" + ("%1.0e"%(number_of_sph_particles)).replace("+0","").replace("+","")
    t_end_string = "t" + ("%1.0e"%(t_end.value_in(units.s))).replace("+0","").replace("+","")

    base_output_file_name = os.path.join(get_path_to_results(), "stellar_merger_"+n_string+"_"+t_end_string)

    stars =  Particles(2)
    stars.mass = masses
    try:
        stellar_evolution = MESA()
        stellar_evolution.initialize_code()
    except:
        print "MESA was not built. Returning."
        return
    stellar_evolution.commit_parameters()
    stellar_evolution.particles.add_particles(stars)
    stellar_evolution.commit_particles()
    print "Evolving stars with MESA..."
    stellar_evolution.evolve_model(star_age)

    number_of_sph_particles_1 = int(round(number_of_sph_particles *
        (stellar_evolution.particles[0].mass / stellar_evolution.particles.mass.sum())))
    number_of_sph_particles_2 = number_of_sph_particles - number_of_sph_particles_1
    print "Creating initial conditions from a MESA stellar evolution model:"
    print stellar_evolution.particles[0].mass, "star consisting of", number_of_sph_particles_1, "particles."
    sph_particles_1 = convert_stellar_model_to_SPH(
        stellar_evolution.particles[0],
        number_of_sph_particles_1,
        seed=12345
    ).gas_particles
    print stellar_evolution.particles[1].mass, "star consisting of", number_of_sph_particles_2, "particles."
    sph_particles_2 = convert_stellar_model_to_SPH(
        stellar_evolution.particles[1],
        number_of_sph_particles_2
    ).gas_particles

    initial_separation += stellar_evolution.particles.radius.sum()
    sph_particles_2.x  += numpy.cos(angle) * initial_separation
    sph_particles_2.y  += numpy.sin(angle) * initial_separation
    sph_particles_1.vx += numpy.cos(angle) * initial_speed - numpy.sin(angle) * initial_speed_perpendicular
    sph_particles_1.vy += numpy.cos(angle) * initial_speed_perpendicular + numpy.sin(angle) * initial_speed
    view = [-0.5, 0.5, -0.5, 0.5] * (initial_separation + stellar_evolution.particles.radius.sum())
    stellar_evolution.stop()

    all_sph_particles = ParticlesSuperset([sph_particles_1, sph_particles_2])
    all_sph_particles.move_to_center()

    unit_converter = ConvertBetweenGenericAndSiUnits(1.0 | units.RSun, constants.G, t_end)
    hydro_legacy_code = sph_code(unit_converter)
    n_steps = 100
    hydro_legacy_code.parameters.n_smooth = 96
    try:
        hydro_legacy_code.parameters.timestep = t_end / n_steps
    except Exception as exc:
        if not "parameter is read-only" in str(exc): raise
    hydro_legacy_code.gas_particles.add_particles(all_sph_particles)

    print "Evolving to t =", t_end, " (using", sph_code.__name__, "SPH code)."
    for time, i_step in [(i*t_end/n_steps, i) for i in range(1, n_steps+1)]:
        hydro_legacy_code.evolve_model(time)
        if not i_step % 4:
            hydro_plot(
                view,
                hydro_legacy_code,
                (300, 300),
                base_output_file_name + "_hydro_image{0:=03}.png".format(i_step)
            )
    hydro_legacy_code.stop()
    print "All done!\n"

def hydro_plot(view, hydro_code, image_size, figname):
    """
    view: the (physical) region to plot [xmin, xmax, ymin, ymax]
    hydro_code: hydrodynamics code in which the gas to be plotted is defined
    image_size: size of the output image in pixels (x, y)
    """
    shape = (image_size[0], image_size[1], 1)
    size = image_size[0] * image_size[1]
    axis_lengths = [0.0, 0.0, 0.0] | units.m
    axis_lengths[0] = view[1] - view[0]
    axis_lengths[1] = view[3] - view[2]
    grid = Grid.create(shape, axis_lengths)
    grid.x += view[0]
    grid.y += view[2]
    speed = grid.z.reshape(size) * (0 | 1/units.s)
    rho, rhovx, rhovy, rhovz, rhoe = hydro_code.get_hydro_state_at_point(grid.x.reshape(size),
        grid.y.reshape(size), grid.z.reshape(size), speed, speed, speed)

    min_v =  800.0 | units.km / units.s
    max_v = 3000.0 | units.km / units.s
    min_rho = 3.0e-4 | units.g / units.cm**3
    max_rho = 0.1 | units.g / units.cm**3
    min_E = 1.0e11 | units.J / units.kg
    max_E = 1.0e13 | units.J / units.kg

    v_sqr = (rhovx**2 + rhovy**2 + rhovz**2) / rho**2
    E = rhoe / rho
    log_v = numpy.log((v_sqr / min_v**2)) / numpy.log((max_v**2 / min_v**2))
    log_rho = numpy.log((rho / min_rho)) / numpy.log((max_rho / min_rho))
    log_E = numpy.log((E / min_E)) / numpy.log((max_E / min_E))

    red   = numpy.minimum(numpy.ones_like(rho.number), numpy.maximum(numpy.zeros_like(rho.number), log_rho)).reshape(shape)
    green = numpy.minimum(numpy.ones_like(rho.number), numpy.maximum(numpy.zeros_like(rho.number), log_v)).reshape(shape)
    blue  = numpy.minimum(numpy.ones_like(rho.number), numpy.maximum(numpy.zeros_like(rho.number), log_E)).reshape(shape)
    alpha = numpy.minimum(numpy.ones_like(log_v), numpy.maximum(numpy.zeros_like(log_v),
        numpy.log((rho / (10*min_rho))))).reshape(shape)

    rgba = numpy.concatenate((red, green, blue, alpha), axis = 2)

    pyplot.figure(figsize = (image_size[0]/100.0, image_size[1]/100.0), dpi = 100)
    im = pyplot.figimage(rgba, origin='lower')

    pyplot.savefig(figname, transparent=True, dpi = 100)
    print "\nHydroplot was saved to: ", figname
    pyplot.close()

if __name__  == "__main__":
    print "Running the simulation that formed the basis of the christmas card of Leiden Observatory of 2010."
    print
    print "Details:"
    print "The ornaments are the result of a smoothed particle simulation " \
        "with 50000 equal mass particles of a 310 Myr star of 0.3 solar mass, " \
        "which is ejected from a distance of 4 solar radii (left) " \
        "with a velocity of 3000 km/s into a 3.0 solar mass star at an " \
        "age of 310 Myr. The calculation was performed using the AMUSE " \
        "(amusecode.org) software environment in which the stars were evolved " \
        "using MESA to an age of 310 Myr before the encounter was performed " \
        "using Fi. Each ornament, generated using pyplot, is a snapshot from the " \
        "simulation, from top left to bottom right. The peak is created from a " \
        "blend of all snapshots. The colors of the ornaments are, red: log " \
        "of the density, green: log of the speed and for blue we used " \
        "the log of the specific internal energy."
    print
    if HAS_MATPLOTLIB:
        head_on_stellar_merger()
    else:
        print "matplotlib is not installed. Install it in the site-packages folder of your Python installation. Returning."

max_stars_repo_name: egpbos/amuse

max_issues_repo_name: egpbos/amuse

max_forks_repo_name: egpbos/amuse


text:
import cv2
import sys
import os
import tensorflow.contrib.tensorrt as trt
import tensorflow as tf
import numpy as np
import time
from tf_trt_models.detection import download_detection_model, build_detection_graph

IMAGE_PATH = './COCO_val2014_000000581929.jpg'

def open_pb(PB_PATH):
    detection_graph = tf.Graph()
    with detection_graph.as_default():
        od_graph_def = tf.GraphDef()
        with tf.gfile.GFile(PB_PATH, 'rb') as fid:
            serialized_graph = fid.read()
            od_graph_def.ParseFromString(serialized_graph)
            tf.import_graph_def(od_graph_def, name='')
    return detection_graph

trt_graph = open_pb('./data/faster_rcnn_resnet101_coco_trt_16.pb')

tf_config = tf.ConfigProto()
tf_config.gpu_options.allow_growth = True

with trt_graph.as_default():
    tf_sess = tf.Session(config=tf_config)

    #tf.import_graph_def(trt_graph, name='')

    tf_input = tf.get_default_graph().get_tensor_by_name('image_tensor:0')
    tf_scores = tf_sess.graph.get_tensor_by_name('detection_scores:0')
    tf_boxes = tf_sess.graph.get_tensor_by_name('detection_boxes:0')
    tf_classes = tf_sess.graph.get_tensor_by_name('detection_classes:0')
    tf_num_detections = tf_sess.graph.get_tensor_by_name('num_detections:0')

    image_bgr = cv2.imread(IMAGE_PATH)
    image_origin = np.copy(image_bgr)
    image = image_bgr[..., ::-1]

    #image_resized = cv2.resize(image, (300, 300))

    scores, boxes, classes, num_detections = tf_sess.run([tf_scores, tf_boxes, tf_classes, tf_num_detections], feed_dict={
        tf_input: image[None, ...]
    })

    boxes = boxes[0] # index by 0 to remove batch dimension
    scores = scores[0]
    classes = classes[0]
    num_detections = int(num_detections[0])

    rows = image_bgr.shape[0]
    cols = image_bgr.shape[1]
    # plot boxes exceeding score threshold
    for i in range(num_detections):
        classId = classes[i]
        score = scores[i]
        bbox = [float(v) for v in boxes[i]]
        if score > 0.1:
            x = bbox[1] * cols
            y = bbox[0] * rows
            right = bbox[3] * cols
            bottom = bbox[2] * rows
            cv2.rectangle(image_bgr, (int(x), int(y)), (int(right), int(bottom)), (100, 255, 10), 5)

    cv2.imshow("origin", image_origin)
    cv2.imshow("result", image_bgr)
    cv2.waitKey()