text:
import getopt
import os
import pickle
import sys
from functools import reduce
from os import listdir
from os.path import isfile, join

import numpy as np
import pandas as pd
from keras.callbacks import Callback
from keras.callbacks import EarlyStopping
from keras.layers import Dense, Input, concatenate
from keras.layers.core import Dropout
from keras.models import Model
from keras.models import Sequential
from keras.models import load_model
from keras.regularizers import l2
from sklearn.metrics import f1_score as f1_score_func

extractors_types = ['alchemy', 'adel', 'opencalais', 'meaning_cloud', 'dandelion', 'dbspotlight', 'babelfy',
                    'textrazor']
extractors_disambiguation = ['dandelion', 'dbspotlight', 'babelfy', 'textrazor']

reg_alpha = 0.000
dropout = 0.45
dropout_input = 0.0
epochs = 1000
batch = 50
eg_alpha = 0.0
units = 400
layers = 2
# loss_function = 'categorical_crossentropy'
loss_function = 'mse'
optimizer = 'adam'
# activation = 'softmax'
activation = 'sigmoid'
activation_middle = 'selu'
architecture = 'simple'
patience = 50
dim_concatenation = 5

try:
    base = sys.argv[1]
except:
    print('You have to specify the ground truth')
    exit()

optlist, args = getopt.getopt(sys.argv[2:], '', ['reg_alpha=',
                                                 'dropout=',
                                                 'epochs=',
                                                 'model_name=',
                                                 'batch=',
                                                 'eg_alpha=',
                                                 'units=',
                                                 'loss_function=',
                                                 'saving_folder=',
                                                 'features=',
                                                 'optimizer=',
                                                 'activation_middle=',
                                                 'patience=',
                                                 'activation='
                                                 ])

for opt, arg in optlist:
    if opt == '--reg_alpha':
        reg_alpha = float(arg)
    elif opt == "--dropout":
        dropout = float(arg)
    elif opt == "--dropout_input":
        dropout_input = float(arg)
    elif opt == "--epochs":
        epochs = int(arg)
    elif opt == "--batch":
        batch = int(arg)
    elif opt == '--eg_alpha':
        eg_alpha = float(arg)
    if opt == '--units':
        units = int(arg)
    elif opt == "--loss_function":
        loss_function = arg
    elif opt == "--saving_folder":
        saving_folder = arg
    elif opt == "--features":
        features = arg.split(',')
    elif opt == '--optimizer':
        optimizer = arg
    elif opt == "--activation_middle":
        activation_middle = arg
    elif opt == "--activation":
        activation = arg
    elif opt == "--patience":
        patience = int(arg)

features = set(['type', 'score', 'entity', 'fasttext'])

def get_similarities_dict(features_paths_train, features_paths_test):
    sim_dict = dict()
    features_paths = features_paths_train + features_paths_test
    for i, f in enumerate(features_paths):
        similarites = pickle.load(open(f, 'rb'))['features']['entity_MATRIX']
        for item in similarites.items():
            sim_dict[(str(item[0][0]), str(item[0][1]))] = item[1]
            sim_dict[(str(item[0][1]), str(item[0][0]))] = item[1]
    return sim_dict

def getURISListGT(groundtruth_paths):
    uris_df_per_file = list()
    for g in groundtruth_paths:
        recs_gt = pd.read_csv(g).to_dict(orient='records')
        uris = list()
        for r in recs_gt:
            uri = r['wd_uri']
            uris.append(uri)
        uris_df_per_file.append(uris)
    return uris_df_per_file

def setURI(t):
    if bool(t) and type(t) == str:
        return t
    else:
        return np.NAN

def getURISListPerExtractor(features_paths, extractor_list=['dandelion', 'dbspotlight', 'babelfy', 'textrazor']):
    uris_list_per_extractor = list()
    for j, f in enumerate(features_paths):
        uri_list_file = list()
        obj = pickle.load(open(f, 'rb'))
        uri_list = obj['entity_list']
        uris_list_per_extractor.append(
            [{ext: setURI(uri_list[ext][i]) for ext in uri_list if ext in extractor_list} for i in
             range(len(uri_list['dandelion']))])
    return uris_list_per_extractor

def ConcatenateArray(list_arrays):
    array = list()
    for l in list_arrays:
        array += l
    return array

def isSameUri(wd1, wd2):
    return int((type(wd1) == type(wd2) == float) or (wd1 == wd2))

def builtXY_b(uris_list_per_extractor_flatten, extractors_disambiguation, gt_flatten):
    X = []
    Y = []
    candidates_list = []
    for i, line in enumerate(uris_list_per_extractor_flatten):

        candidates = set([line[ext] for ext in extractors_disambiguation if type(line[ext]) == str])
        if True in [type(line[ext]) == float for ext in extractors_disambiguation]:
            candidates.add(np.NAN)
        candidates = list(candidates)
        candidates.sort(key=lambda x: str(x))
        X_p = [ConcatenateArray([[isSameUri(line[ext_1], uri)] for ext_1 in extractors_disambiguation])
               for uri in candidates]
        Y_p = [[isSameUri(uri, gt_flatten[i])]
               for uri in candidates]
        candidates_list.append(candidates)
        X += X_p
        Y += Y_p
    X = np.array(X)
    Y = np.array(Y)
    return X, Y, candidates_list

def builtXY(uris_list_per_extractor_flatten, extractors_disambiguation, gt_flatten, similarities_dict):
    X = []
    Y = []
    candidates_list = []
    for i, line in enumerate(uris_list_per_extractor_flatten):
        candidates = set([line[ext] for ext in extractors_disambiguation if type(line[ext]) == str])
        if True in [type(line[ext]) == float for ext in extractors_disambiguation]:
            candidates.add(np.NAN)
        candidates = list(candidates)
        candidates.sort(key=lambda x: str(x))
        X_p = [
            ConcatenateArray([similarities_dict[(str(line[ext_1]), str(uri))] for ext_1 in extractors_disambiguation])
            for uri in candidates]
        Y_p = [[isSameUri(uri, gt_flatten[i])]
               for uri in candidates]
        candidates_list.append(candidates)
        X += X_p
        Y += Y_p
    X = np.array(X)
    Y = np.array(Y)
    return X, Y, candidates_list

def createModelObject(extractors_disambiguation,
                      X_test, X_train, Y_test, Y_train,
                      gt_test_flatten, gt_train_flatten,
                      uris_list_per_extractor_train_flatten, uris_list_per_extractor_test_flatten, features_paths_train,
                      features_paths_test, candidates_train, candidates_test
                      ):
    model_obj = {
        'features_paths_train': features_paths_train,
        'features_paths_test': features_paths_test,
        'extractors_disambiguation': extractors_disambiguation
    }

    model_obj['candidates_train'] = candidates_train
    model_obj['candidates_test'] = candidates_test
    model_obj['train_X'] = np.array(X_train)
    model_obj['train_Y'] = np.array(Y_train)
    model_obj['test_X'] = np.array(X_test)
    model_obj['test_Y'] = np.array(Y_test)
    return model_obj

def addExtractorsTypesRepresentation(model_obj, extractors_types):
    candidates_test_iter = iter(model_obj['candidates_test'])
    candidates_train_iter = iter(model_obj['candidates_train'])
    paths_test = model_obj['features_paths_test']
    paths_train = model_obj['features_paths_train']
    model_obj['extractors_types'] = extractors_types
    model_obj['type_dict_test'] = {ext: [] for ext in extractors_types}
    model_obj['type_dict_train'] = {ext: [] for ext in extractors_types}
    for p in paths_test:
        obj = pickle.load(open(p, 'rb'))
        features_obj = obj['features']['type']
        flag_g = True
        for ext in extractors_types:
            if flag_g:
                repeats = [len(next(candidates_test_iter)) for line in obj['features']['type'][ext]]
                flag_g = False
            for j, line in enumerate(obj['features']['type'][ext]):
                for k in range(repeats[j]):
                    model_obj['type_dict_test'][ext].append(line)

    for p in paths_train:
        obj = pickle.load(open(p, 'rb'))
        features_obj = obj['features']['type']
        flag_g = True
        for ext in extractors_types:
            if flag_g:
                repeats = [len(next(candidates_train_iter)) for line in obj['features']['type'][ext]]
                flag_g = False
            for j, line in enumerate(obj['features']['type'][ext]):
                for k in range(repeats[j]):
                    model_obj['type_dict_train'][ext].append(line)
    for ext in extractors_types:
        model_obj['type_dict_test'][ext] = np.array(model_obj['type_dict_test'][ext])
        model_obj['type_dict_train'][ext] = np.array(model_obj['type_dict_train'][ext])
    return model_obj

training_folder = 'data/training_data/' + base + '/'
ground_truth_folder_train = training_folder + 'train/csv_ground_truth/'
ground_truth_folder_test = training_folder + 'test/csv_ground_truth/'

features_folder_train = training_folder + 'train/features_files/'
features_folder_test = training_folder + 'test/features_files/'

text_folder_train = training_folder + 'train/txt_files/'
text_folder_test = training_folder + 'test/txt_files/'

features_paths_train = [features_folder_train + f for f in listdir(features_folder_train) if
                        isfile(join(features_folder_train, f)) and '.p' in f]
features_paths_train.sort()
features_paths_test = [features_folder_test + f for f in listdir(features_folder_test) if
                       isfile(join(features_folder_test, f)) and '.p' in f]
features_paths_test.sort()
groundtruth_paths_train = [path.replace(features_folder_train, ground_truth_folder_train).replace('.p', '.csv') for path
                           in features_paths_train]
groundtruth_paths_train.sort()
groundtruth_paths_test = [path.replace(features_folder_test, ground_truth_folder_test).replace('.p', '.csv') for path in
                          features_paths_test]
groundtruth_paths_test.sort()

gt_test = getURISListGT(groundtruth_paths_test)
gt_test_flatten = reduce(lambda x, y: x + y, gt_test)
gt_train = getURISListGT(groundtruth_paths_train)
gt_train_flatten = reduce(lambda x, y: x + y, gt_train)

uris_list_per_extractor_train = getURISListPerExtractor(features_paths_train)
uris_list_per_extractor_train_flatten = reduce(lambda x, y: x + y, uris_list_per_extractor_train)
uris_list_per_extractor_test = getURISListPerExtractor(features_paths_test)
uris_list_per_extractor_test_flatten = reduce(lambda x, y: x + y, uris_list_per_extractor_test)

similarities_dict = get_similarities_dict(features_paths_train, features_paths_test)
X_train, Y_train, candidates_train = builtXY(uris_list_per_extractor_train_flatten, extractors_disambiguation,
                                             gt_train_flatten, similarities_dict)
X_test, Y_test, candidates_test = builtXY(uris_list_per_extractor_test_flatten, extractors_disambiguation,
                                          gt_test_flatten, similarities_dict)

model_obj = createModelObject(extractors_disambiguation,
                              X_test, X_train, Y_test, Y_train,
                              gt_test_flatten, gt_train_flatten,
                              uris_list_per_extractor_train_flatten, uris_list_per_extractor_test_flatten,
                              features_paths_train,
                              features_paths_test, candidates_train, candidates_test
                              )

model_obj = addExtractorsTypesRepresentation(model_obj, extractors_types)

saving_folder = 'data/models/disambiguation/' + base + '/'

try:
    os.makedirs(saving_folder)
except:
    pass
if set(features) == set(['type', 'score', 'entity', 'fasttext']):
    model_file_path = saving_folder + 'model.h5'
else:
    model_file_path = saving_folder + '_'.join(features) + '_model.h5'

class TestCallback(Callback):
    def __init__(self):
        self.history_scores = list()
        self.f1_max = 0
        self.val_loss_min = 20000

    def on_epoch_end(self, epoch, logs={}):
        global model_obj
        x, y = model_obj['X_dict_test'], model_obj['Y_dict_test']
        val_loss = model.evaluate(x, y)[0]
        predicted_test = self.model.predict(x, verbose=0)
        score_obj = dict()

        score_obj['val_loss'] = val_loss
        f1_score = f1_score_func(np.squeeze(model_obj['test_Y']), np.squeeze(predicted_test).round())
        score_obj['f1_score'] = f1_score

        print(score_obj, len(predicted_test))

        self.history_scores.append(score_obj)
        if ((f1_score > self.f1_max) and (val_loss <= self.val_loss_min)) or (
                (val_loss < self.val_loss_min) and (f1_score >= self.f1_max)):
            self.f1_max = f1_score
            self.val_loss_min = val_loss
            self.model.save(model_obj['path'])

    def on_train_end(self, logs={}):
        global model_obj
        model_obj['history_scores'] = self.history_scores
        md = load_model(model_obj['path'])
        predicted_test = md.predict(model_obj['X_dict_test'], verbose=0)
        predicted_train = md.predict(model_obj['X_dict_train'], verbose=0)
        model_obj['predicted_train'] = predicted_train
        model_obj['predicted_test'] = predicted_test

model_obj['path'] = model_file_path

dim_in_1 = model_obj['train_X'].shape[-2]
dim_in_2 = model_obj['train_X'].shape[-1]
dim_out = model_obj['train_Y'].shape[-1]

def generateConcactPartSimple(X, ext_name):
    dim_input = X.shape[-1]
    input_tensor = Input(shape=(X.shape[-1],), name=ext_name)
    dense_middle = Dense(dim_concatenation, activation=activation_middle, kernel_regularizer=l2(reg_alpha),
                         bias_regularizer=l2(reg_alpha))(input_tensor)
    return input_tensor, dense_middle

if 'type_dict_test' not in model_obj:
    model = Sequential()
    if architecture == 'blstm':
        raise Exception('Not defined BLSTM')
    elif architecture == 'simple':
        model.add(Dropout(dropout_input, input_shape=(dim_in_2,)))
        for i in range(layers):
            model.add(Dense(units, activation=activation_middle, kernel_regularizer=l2(reg_alpha),
                            bias_regularizer=l2(reg_alpha)))
            model.add(Dropout(dropout))
        model.add(
            Dense(dim_out, activation=activation, kernel_regularizer=l2(reg_alpha), bias_regularizer=l2(reg_alpha)))

    model.compile(loss=loss_function, optimizer=optimizer, metrics=['mae', 'accuracy'])
    print(model.summary())
    early_stop = EarlyStopping(monitor='val_loss', min_delta=0, patience=patience, verbose=0)
    model_obj['X_dict_test'] = model_obj['test_X']
    model_obj['Y_dict_test'] = model_obj['test_Y']
    model_obj['X_dict_train'] = model_obj['train_X']
    model.fit(model_obj['train_X'], model_obj['train_Y'], epochs=epochs, batch_size=batch,
              callbacks=[TestCallback(), early_stop], validation_data=(model_obj['test_X'], model_obj['test_Y']))

else:

    X_dict_train = dict()
    X_dict_train['entity'] = model_obj['train_X']
    for ext in model_obj['type_dict_train']:
        X_dict_train[ext] = model_obj['type_dict_train'][ext]

    X_dict_test = dict()
    X_dict_test['entity'] = model_obj['test_X']
    for ext in model_obj['type_dict_test']:
        X_dict_test[ext] = model_obj['type_dict_test'][ext]

    uris_tensor = Input(shape=(dim_in_2,), name='entity')
    uris_tensor_d = Dropout(dropout_input)(uris_tensor)
    to_concatenate_layers = [generateConcactPartSimple(X=model_obj['type_dict_train'][ext], ext_name=ext) for ext in
                             model_obj['type_dict_train']]
    to_concatenate_layers_input = [uris_tensor] + [c[0] for c in to_concatenate_layers]
    to_concatenate_layers_out = [uris_tensor] + [c[1] for c in to_concatenate_layers]
    concatenation = concatenate(to_concatenate_layers_out)
    for i in range(layers):
        concatenation = Dense(units, activation=activation_middle, kernel_regularizer=l2(reg_alpha),
                              bias_regularizer=l2(reg_alpha))(concatenation)
        concatenation = Dropout(dropout)(concatenation)
    main_output = Dense(dim_out, activation=activation, kernel_regularizer=l2(reg_alpha),
                        bias_regularizer=l2(reg_alpha), name='main_output')(concatenation)
    model = Model(inputs=to_concatenate_layers_input, outputs=[main_output])

    model.compile(loss=loss_function, optimizer=optimizer, metrics=['mae', 'accuracy'])
    print(model.summary())
    early_stop = EarlyStopping(monitor='val_loss', min_delta=0, patience=patience, verbose=0)

    model_obj['X_dict_test'] = X_dict_test
    model_obj['Y_dict_test'] = {'main_output': model_obj['test_Y']}
    model_obj['X_dict_train'] = X_dict_train

    model.fit(X_dict_train,
              {'main_output': model_obj['train_Y']},
              validation_data=(X_dict_test,
                               {'main_output': model_obj['test_Y']}),
              callbacks=[TestCallback(), early_stop],
              epochs=epochs, batch_size=batch)


text:
from sklearn import datasets
import matplotlib.pyplot as plt

from sklearn.model_selection import learning_curve,validation_curve
from sklearn.svm import SVC
from sklearn.datasets import load_digits
import numpy as np
import pickle
import joblib

def test_make_regression():
    '''创建回归数据'''
    X,y=datasets.make_regression(n_samples=100,n_features=1,n_targets=1,noise=2)
    plt.scatter(X,y)
    plt.show()
    from sklearn import metrics
    print(sorted(metrics.SCORERS.keys()))
    return

def test_over_fitting1():
    '''
    过拟合模拟-learning_curve
    gamma=0.01 train_size=0.25 交叉验证最优
    gamma=0.001，交叉验证损失和train_size成反比
    :return:
    '''
    (X,y)=load_digits(return_X_y=True)
    # train_sizes表示划分数据集时，训练集的比重
    train_sizes,train_loss,test_loss=learning_curve(SVC(gamma=0.01),X,y,cv=10,
                                                    train_sizes=[0.1,0.25,0.5,0.75,1],
                                                    scoring='neg_mean_squared_error')
    train_loss=-np.mean(train_loss,axis=1)
    test_loss=-np.mean(test_loss,axis=1)

    plt.plot(train_sizes,train_loss,'o-',color='r',label='train')
    plt.plot(train_sizes,test_loss,'o-',color='g',label='cross')
    plt.xlabel('train_size')
    plt.ylabel('loss')
    # 右上角显示图例
    plt.legend(loc='best')
    plt.show()

    return

def test_over_fitting2():
    '''
    过拟合模拟-validation_curve
    gamma值范围[0.0005,0.002] 共10个点，在第7个点时交叉验证最优，第8个点开始过拟合
    :return:
    '''
    (X,y)=load_digits(return_X_y=True)
    gamma_range=np.linspace(0.0005,0.002,10)
    # train_sizes表示划分数据集时，训练集的比重
    train_loss,test_loss=validation_curve(SVC(),X,y,param_name='gamma',param_range=gamma_range,cv=10,

                                                    scoring='neg_mean_squared_error')
    # axis=1按行求均值
    train_loss=-np.mean(train_loss,axis=1)
    test_loss=-np.mean(test_loss,axis=1)

    plt.plot(gamma_range,train_loss,'o-',color='r',label='train')
    plt.plot(gamma_range,test_loss,'o-',color='g',label='test')
    plt.xlabel('gamma')
    plt.ylabel('loss')
    # 右上角显示图例
    plt.legend(loc='best')
    plt.show()
    return

def test_save_restore_model():
    '''
    保存和加载模型
    :return:
    '''
    (X,y)=datasets.load_iris(return_X_y=True)
    model=SVC()
    model.fit(X,y)
    # # 保存 model.joblib后缀可以任意
    with open('./tmp/model.joblib','wb') as f:
        # pickle.dump(model,f)
        joblib.dump(model,f)
    # 加载
    with open('./tmp/model.joblib','rb') as f:
        # model2=pickle.load(f)
        model2=joblib.load(f)
        print(model2.predict(X[:2]))

    return

def test_file():
    import os
    import pathlib
    cu=os.getcwd()
    print(cu)
    f1=os.path.join(cu,'tmp/model/d1/d2/t.txt')
    flag1=os.path.isdir(f1)
    flag2=os.path.isfile(f1)
    print(flag1,flag2)
    print(os.path.abspath(os.path.dirname(f1)))
    pathlib.Path(f1).touch()
    return

if __name__ == '__main__':
    print('main start')
    # test_over_fitting2()
    # test_save_restore_model()
    test_file()
    print('main end')

ext: py

reason: import numpy


text:
# coding : utf-8

import sys
import os
module_path = os.path.abspath( u'..'  )
sys.path.append( module_path )

import json
import cPickle
import numpy as np
import theano
from vaegan.data import load_data, postprocess
from vaegan import Operate, image_grid
from vaegan import VAEGAN

def load_json( infile_name ):
    fin = open( infile_name, u'r' )
    conf = json.load( fin )
    fin.close()
    return conf

conf = load_json( sys.argv[1] )

param_path = conf[u'param_path']
data_dir   = conf[u'test_data_dir']
output_dir = conf[u'output_dir']
box        = conf[u'box'] # cropped region
shape = (64, 64)

n_samples = 5

seed = 1234
np.random.seed( seed )

# Build VAEGAN
vaegan = VAEGAN( seed )

try:
    fin = open( param_path, u'rb' )
    params = cPickle.load( fin )
    params = [ p.astype( theano.config.floatX ) for p in params ]
    fin.close()
except:
    print >>sys.stderr, u'Cannot open file:', param_path
    sys.exit(-1)

vaegan.load_params( params )

data = load_data( data_dir, shape, box )

indices = np.arange( len(data) )
np.random.shuffle( indices )
sources = data[ indices[:n_samples] ]
np.random.shuffle( indices )
data1 = data[ indices[:n_samples] ]
np.random.shuffle( indices )
data2 = data[ indices[:n_samples] ]

operate = Operate( vaegan, sources, output_dir )
operate.plus(data1)
operate.minus(data2)
operate.equal()

source_arrs = np.asarray( map( postprocess, sources ) )
source_imgs = image_grid( source_arrs, (n_samples,1) )
source_imgs.save( u'%s/sources.jpeg' % output_dir )

data1_arrs  = np.asarray( map( postprocess, data1 ) )
data1_imgs  = image_grid( data1_arrs, (n_samples,1) )
data1_imgs.save( u'%s/data1.jpeg' % output_dir )

data2_arrs  = np.asarray( map( postprocess, data2 ) )
data2_imgs  = image_grid( data2_arrs, (n_samples,1) )
data2_imgs.save( u'%s/data2.jpeg' % output_dir )

max_stars_repo_name: amirjaber/vaegan

max_issues_repo_name: amirjaber/vaegan

max_forks_repo_name: amirjaber/vaegan

reason: import numpy,import theano


text:
import os, random, argparse
import numpy as np
import imghdr
from PIL import Image

def getAverageRGB(image):
  im = np.array(image)
  w, h, d = im.shape
  return tuple(np.average(im.reshape(w*h, d), size=0))

def getAverageRGBOld(image):
  npixels = image.size[0]*image.size[1]
  cols = image.getcolors(npixels)
  sumRGB = [(x[0]*x[1][0], x[0]*x[1][1], x[0]*x[1][2]) for x in cols]
  avg = tuple([int(sum(x)/npixels) for x in zip(*sumRGB)])
  return avg

def splitImage(image, size):
  W, H = image.size[0], image.size[1]
  m, n = size
  w, h = int(W/m), int(H/n)
  img = []

  for j in range(m):
    for i in range(n):
      img.append(image.crop((i*w, j*h, (i+1)*w, (j+1)*h )))

  return img

def getFileNames(imagedir):

  files = os.listdir(imagedir)
  filenames = []
  for file in files:
    filePath = os.path.abspath(os.path(imagedir, file))
  try:
    imgtype = imghdr.what(filePath)
  if imgtype:
    filenames.append(filePath)
  except:
    print("string is invalid %s" (filePath,))
  return filenames

def getImages(imagedir):
  files = os.listdir(imagedir)
  images = []
  for file in files:
    filepath = os.path.abspath(os.path(imagedir, file))
    try:
      fp = open(filepath, "rb")
      im = image.open(fp)
      images.append(im)

    except:
      print("string is invalid %s", (filepath,))
  return images

def getMinIndex(input_avg, avgs):

  avg = input_avg

  min_index = 0
  index = 0
  min_dist = float("inf")
  for val in avgs:
    dist = ((val[0]-avg[0])*(val[0]-avg[0]) +
            (val[1]-avg[1])*(val[1]-avg[1]) +
            (val[2]-avg[2])*(val[2]-avg[2]))
    if min_dist< dist:
      dist = min_dist
      min_index = index
      index +=1;

  return min_index

def createImageGrid(images, dim):

  m,n = dim

  assert m*n == len(images)

  width = ([img.size[0] for img in images])
  height = ([img.size[1] for img in images])

  grid_image = Image.new("RGB", (n*width, m*height))

  for index in range(len(images)):
    row = int(index/n)
    col = index - n*row
    grid_image.paste(images[index], (width* col, height * row))

  return grid_image

def createPhotomosaic(target_image, input_images, grid_size, resuse_images=True ):

  print("splitting the input image")

  target_images = splitImage(target_image, grid_size)

  print("gathering the matching images")

  output_images = []
  count = 0
  batch_size = int(len(target_images)/10)

  avgs = []
  for img in input_images:
    avgs.append(getAverageRGB(img))

  for img in target_image:
    avg = getAverageRGB(img)

  match_index = getMinIndex(avg, avgs)
  output_images.append(input_images(match_index))

  if count>0 and batch_size > 10 and count % batch_size is 0:
    print("processed %d of %d", %(count, len(target_images)))
    count+=1

  if not resuse_images:
    input_images.remove(match)

  print("creating the mosaic")
  mosaic_image = createPhotomosaic(output_images, grid_size)

  return mosaic_image

def main():

  parser = argparse.ArgumentParser
  (description= "Creates a photomosaic from the input images")

  parser.add_arugment('--target-image', dest = target_image, required = True)
  parser.add_arugment('--input-folder', dest = input_folder, required = True)
  parser.add_arugment('--grid_size', nargs = 2, dest = grid_size, required = True)
  parser.add_arugment('--output-file', dest = outfile, required = False)

  args = parser.parse_args()

  target_image = Image.open(args.target_image)

  print("reading the input folder")
  input_images = getImages(args.input_folder)

  if input_images == []:
  print("no images found in the %s folder" %(args.input_folder, ))
  exit()

  random.shuffle(input_images)

  grid_size = (int(args.grid_size[0]), int(args.grid_size[1]))

  output_filename = 'mosaic.png'
  if args,outfile:
    output_filename = args.outfile

  resuse_images = True

  resize_input = True

  print("starting the photomosaic Creation process")

  if not resuse_images:
    if grid_size[0]*grid_size[1] > len(input_images):
      print("the size of the grid is less than the size of the input images")
      exit()

  if resize_input:
    print("resizing the input")

    dims = (int(target_image.size[0]/grid_size[1]),
            int(target_image.size[1]/grid_size[0]))
    print("max tile dims is %d" %(dims, ))

    for img in input_images:
      img.thumbnail(dims)

  mosaic_image = createPhotomosaic(target_image, input_images, grid_size, resuse_images)

  mosaic_image.save(output_filename, 'PNG')

  print("saved image to %s" %(output_filename, ))
  print("done")

if __name__ = "__main__":
main()

max_stars_repo_name: rakesh-chinta/PhotoMosaic-Python

max_issues_repo_name: rakesh-chinta/PhotoMosaic-Python

max_forks_repo_name: rakesh-chinta/PhotoMosaic-Python


text:
"""
--- Ångström ---
Tests cell edges calculation.
"""
from angstrom.molecule import Cell
import numpy as np

def test_cubic_unit_cell_edges_calculation():
    a, b, c = 5, 5, 5
    expected_edges = [[[a, 0, 0], [0, 0, 0]], [[0, b, 0], [0, 0, 0]], [[0, 0, c], [0, 0, 0]],
                      [[a, b, 0], [a, 0, 0]], [[a, 0, c], [a, 0, 0]], [[a, b, 0], [0, b, 0]],
                      [[0, b, c], [0, b, 0]], [[0, b, c], [0, 0, c]], [[a, 0, c], [0, 0, c]],
                      [[a, b, c], [a, b, 0]], [[a, b, c], [0, b, c]], [[a, b, c], [a, 0, c]]]
    cell = Cell([a, b, c, 90.0, 90.0, 90.0])
    assert np.allclose(cell.edges, expected_edges)

max_stars_repo_name: kbsezginel/angstrom

max_issues_repo_name: kbsezginel/angstrom

max_forks_repo_name: kbsezginel/angstrom


text:
import cv2
import os
import random
import numpy as np
import torch
from torch.utils import data as data
from torchvision.transforms.functional import normalize

from basicsr.data.degradations import add_jpg_compression
from basicsr.data.transforms import augment, mod_crop, paired_random_crop
from basicsr.utils import FileClient, imfrombytes, img2tensor, scandir
from basicsr.utils.registry import DATASET_REGISTRY

def make_cit_dataset(path):
    clips_path = []

    for index, folder in enumerate(os.listdir(path)):
        clips_folder = os.path.join(path, folder)
        if not (os.path.isdir(clips_folder)):
            continue
        clips_path.append([])

        for image in sorted(os.listdir(clips_folder)):
            clips_path[index].append(os.path.join(clips_folder, image))
    return clips_path

@DATASET_REGISTRY.register()
class CITDataset(data.Dataset):
    """Continuous image transition dataset.
    """

    def __init__(self, opt):
        super(CITDataset, self).__init__()
        self.opt = opt
        # file client (io backend)
        self.file_client = None
        self.io_backend_opt = opt['io_backend']
        self.mean = opt['mean'] if 'mean' in opt else None
        self.std = opt['std'] if 'std' in opt else None
        self.task = opt['task']
        self.extend_t = opt['extend_t'] if 'extend_t' in opt else False

        self.cit_folders = opt['dataroot_cit']

        self.paths = []
        for cit_folder in self.cit_folders:
            self.paths += make_cit_dataset(cit_folder)

    def __getitem__(self, index):
        if self.file_client is None:
            self.file_client = FileClient(self.io_backend_opt.pop('type'), **self.io_backend_opt)

        # Load gt and lq images. Dimension order: HWC; channel order: BGR;
        # image range: [0, 1], float32.
        num_frame = len(self.paths[index])
        idxes = sorted(random.sample(range(num_frame), 3))
        if self.task == 'i2i' or self.task == 'morphing':
            idxes[0], idxes[2] = 0, num_frame - 1

        if self.opt['phase'] == 'train':
            if random.random() < 0.5:
                idxes.reverse()

            # optional tricks
            if random.random() < 0.25 and self.task != 'vfi':
                idxes[1] = num_frame-1 if self.task == 'i2i' else (num_frame-1)//2

        t = np.float32((idxes[1] - idxes[0]) / (idxes[2] - idxes[0])) # float32

        # extend t from [0, 1] to [-1, 1]
        if self.extend_t and random.random() < 0.5:
            idxes = [num_frame-1-idx for idx in idxes]
            t = -t

        frame_paths = [self.paths[index][idx] for idx in idxes]

        img_bytes = self.file_client.get(frame_paths[1], 'frame_gt')
        img_gt = imfrombytes(img_bytes, float32=True)
        img_bytes = self.file_client.get(frame_paths[0], 'frame_0')
        img_0 = imfrombytes(img_bytes, float32=True)
        img_bytes = self.file_client.get(frame_paths[2], 'frame_1')
        img_1 = imfrombytes(img_bytes, float32=True)

        # augmentation for training
        if self.opt['phase'] == 'train':
            gt_size = self.opt['gt_size']
            if self.task == 'vfi':
                # random crop
                img_gt, img_in = paired_random_crop(img_gt, [img_0, img_1], gt_size, 1)
                img_0, img_1 = img_in
            else:
                img_gt = cv2.resize(img_gt, (gt_size, gt_size))
                img_0 = cv2.resize(img_0, (gt_size, gt_size))
                img_1 = cv2.resize(img_1, (gt_size, gt_size))
            # flip, rotation
            img_gt, img_0, img_1 = augment([img_gt, img_0, img_1], self.opt['use_flip'], self.opt['use_rot'])
        else:
            h, w = img_gt.shape[:2]
            nh, nw = h//32*32, w//32*32
            img_gt, img_0, img_1 = img_gt[:nh,:nw], img_0[:nh,:nw], img_1[:nh,:nw]

        # BGR to RGB, HWC to CHW, numpy to tensor
        img_gt, img_0, img_1 = img2tensor([img_gt, img_0, img_1], bgr2rgb=True, float32=True)
        # normalize
        if self.mean is not None or self.std is not None:
            normalize(img_0, self.mean, self.std, inplace=True)
            normalize(img_1, self.mean, self.std, inplace=True)
            normalize(img_gt, self.mean, self.std, inplace=True)

        inputs = img_0 if self.task=='i2i' else torch.cat((img_0, img_1), 0)
        return {'in': inputs, 'gt': img_gt, 't': t, 'lq_path': frame_paths[0], 'gt_path': frame_paths[1]}

    def __len__(self):
        return len(self.paths)

max_stars_repo_name: yangxy/SDL

max_issues_repo_name: yangxy/SDL

max_forks_repo_name: yangxy/SDL


text:
import numpy as np
import torch

def blackpad_augment(sample):
  """Performs flip augmentation on img.

  Args:
    sample: Dictionary of (np array: <ch,z,x,y>) image and mask
  """
  w = sample["image"].shape[1]
  n_pad = np.random.choice(np.arange(int(w*0.4),int(w*0.9),10))

  if np.random.rand() < 0.25:

    d = np.random.choice([0,1])
    for k in sample.keys():
      sample[k] = blackpad(sample[k], d=d, n=n_pad)

  if np.random.rand() < 0.25:

    d = np.random.choice([2,3])
    for k in sample.keys():
      sample[k] = blackpad(sample[k], d=d, n=n_pad)

  return sample

def blackpad(img, d=0, n=100):

  img_pad = np.zeros(img.shape, dtype="float32")

  if d == 0:
    img_pad[:,n:,:] = img[:,:-n,:]

  elif d == 1:
    img_pad[:,:-n,:] = img[:,n:,:]

  elif d == 2:
    img_pad[:,:,n:] = img[:,:,:-n]

  elif d == 3:
    img_pad[:,:,:-n] = img[:,:,n:]

  return img_pad

max_stars_repo_name: jabae/detectEM

max_issues_repo_name: jabae/detectEM

max_forks_repo_name: jabae/detectEM


text:
import numpy as np
import torch
import torch.nn as nn

from dr.ppo.utils import weights_init

class Policy(nn.Module):

    def __init__(self, ob_space, ac_space, hid_size, pol_init_std):
        super().__init__()

        self.l_in = nn.Linear(ob_space.shape[0], hid_size)
        self.l1 = nn.Linear(hid_size, hid_size)
        self.l_out = nn.Linear(hid_size, ac_space.shape[0])

        self.std = nn.Parameter(torch.tensor([[pol_init_std] * ac_space.shape[0]], dtype=torch.float32))

        for name, c in self.named_children():
            if name == 'l_out':
                weights_init(c, 0.01)
            else:
                weights_init(c, 1.0)

    def forward(self, x):

        x = torch.tanh(self.l_in(x))
        x = torch.tanh(self.l1(x))
        mean = self.l_out(x)

        ac = torch.normal(mean, self.std.expand(mean.shape[0], -1))

        return ac, mean

    def neglogp(self, states, acs):

        _, mean = self.forward(states)

        ac_size = acs.size()[-1]

        return 0.5 * torch.sum(((acs - mean) / self.std) ** 2, dim=-1, keepdim=True) + \
               0.5 * np.log(2.0 * np.pi) * float(ac_size) + \
               torch.sum(torch.log(self.std), dim=-1)

    def logp(self, state, ac):
        return - self.neglogp(state, ac)

    def prob(self, state, ac):
        return torch.exp(self.logp(state, ac))

class ValueNet(nn.Module):

    def __init__(self, ob_space, hid_size):
        super().__init__()

        self.l_in = nn.Linear(ob_space.shape[0], hid_size)
        self.l1 = nn.Linear(hid_size, hid_size)
        self.l_out = nn.Linear(hid_size, 1)

        for c in self.children():
            weights_init(c, 1.0)

    def forward(self, x):
        x = torch.tanh(self.l_in(x))
        x = torch.tanh(self.l1(x))
        x = self.l_out(x)

        return x


text:
#!/usr/bin/env python

import argparse
import requests
import numpy
import json

parser = argparse.ArgumentParser()
parser.add_argument('--infile', type=str, help='file to test', default='bulbul.jpg')

args = parser.parse_args()

data = open(args.infile,'rb').read()

res = requests.post(url='http://817c1c6b-37b5-43a4-a947-9e252b665ead.westeurope.azurecontainer.io/score',
                    data=data, headers={'Content-Type':'application/octet-stream'})
class_dict=json.loads(open('imagenet_class_index.json', 'r').read())

out = json.loads(res.text)

predicted_index=numpy.argmax(out['result'])

print(f'predicted class{class_dict[str(predicted_index)]} \
        with index {predicted_index} \
        with probability {out["result"][predicted_index]}')

max_stars_repo_name: mahvas/simple-blackbox-attack

max_issues_repo_name: mahvas/simple-blackbox-attack

max_forks_repo_name: mahvas/simple-blackbox-attack


text:
import scipy.io as sio
import matplotlib.pyplot as plt
import numpy as np
from numpy.random import multivariate_normal as normal
import scipy.stats as stats
# =============================================================================
# 1.Sequential Bayesian Learning
# =============================================================================
#%%
data=sio.loadmat('1_data.mat')
x = data['x']
t = data['t']
j = np.arange(0,3,1)
mu = j*2/3
beta = 1
alpha = 10**-6
N = [5,10,30,80]
s = 0.1
a = []
for i in range(3):
    a.append((x - mu[i])/0.1)
#%%
def sigmoid(x):
    return (1 / (1 + np.exp(-x)))
a1 = np.concatenate((sigmoid(a[0]),sigmoid(a[1]),sigmoid(a[2])),axis=1)
#%%
for n in N:
#   1-1
    plt.figure()
    plt.plot(x[:n],t[:n],'o')
    target = t[:n]
    PHI = a1[:n]
    s0_inv = (10**-6) * np.identity(3)
    sn_inv = s0_inv + PHI.T.dot(PHI)
    sn = np.linalg.inv(sn_inv)
    mn = sn.dot(PHI.T.dot(target))
    mn_ = np.reshape(mn, 3)
    line = np.linspace(0., 2., 50)
    ws = normal(mn_, sn, 5)
    for w in ws:
        line1 = []
        for i in range(3):
            line1.append((line - mu[i])/0.1)
        line1 = np.concatenate((sigmoid(line1[0]).reshape(50,1),sigmoid(line1[1]).reshape(50,1),sigmoid(line1[2]).reshape(50,1)),axis=1)
        value = []
        for point in line1:
            value += [point.T.dot(w)]
        plt.plot(line, value, linestyle ='--', zorder = 1)
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.title('N='+str(n))
    plt.savefig('%d.png'%n)

#    1-2
    plt.figure()
    plt.plot(x[:n],t[:n],'o')
    line2 = []
    for i in range(3):
        line2.append((line - mu[i])/0.1)
    line2 = np.concatenate((sigmoid(line2[0]).reshape(50,1),sigmoid(line2[1]).reshape(50,1),sigmoid(line2[2]).reshape(50,1)),axis=1)
    mx = []
    vx = []
    for point in line2:
        mx += [mn.T.dot(point)] #預測分佈的平均
        vx += [1. + (point.T.dot(sn)).dot(point)] #預測分佈的變異數
    mx = np.reshape(np.asarray(mx), len(mx))
    vx = np.reshape(np.asarray(vx), len(vx))
    plt.plot(line, mx, linestyle = '-', zorder = 1, color = 'red')
    plt.fill_between(line, mx-vx, mx+vx, color = 'pink') #範圍內塗色
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.title('N='+str(n))
    plt.savefig('pred_%d.png'%n)

#%%
#    1-3
j1 = np.arange(0,2,1)
mu = j1*2/3
beta = 1
alpha = 10**-6
N = [5,10,30,80]
s = 0.1
a = []
for i in range(2):
    a.append((x - mu[i])/0.1)

def sigmoid(x):
    return (1 / (1 + np.exp(-x)))
a1 = np.concatenate((sigmoid(a[0]),sigmoid(a[1])),axis=1)

for n in N:
    target = t[:n]
    PHI = a1[:n]
    s0_inv = (10**-6) * np.identity(2)
    sn_inv = s0_inv + PHI.T.dot(PHI)
    sn = np.linalg.inv(sn_inv)
    mn = sn.dot(PHI.T.dot(target))
    mn_ = np.reshape(mn, 2)
    w = np.linspace(-10, 10, 100)
    W = np.dstack(np.meshgrid(w, w))
    prior_vals = stats.multivariate_normal(mn_, sn).pdf(W)
    plt.figure()
    plt.contourf(w, w, prior_vals, 100)
    plt.xlabel('w0')
    plt.ylabel('w1')
    plt.title('N='+str(n))
    plt.savefig('a_%d.png'%n)

reason: import numpy,from numpy,import scipy


text:
from . import ImStats
from os import listdir
from bokeh.charts import Histogram, show
from bokeh.layouts import row
import pandas as pd
from scipy import stats

class ImAggregate(object):
    def __init__(self, dirname, filename="data", loadfile=False, verbose=False):
        filename = filename + '.h5'
        self.filename = filename
        self.dirname = dirname
        self.stats = ["means", "rms", "sums", "sum^2s", "variances", "medians", "iqrs", "nuniqcolors", "dimensions",
                      "stddevs", "counts", "names"]
        self.loaded = False
        if loadfile is False:
            self.data = pd.HDFStore(self.dirname + self.filename)
            stuff = listdir(dirname)
            self.loaded = True
            means = []
            rms = []
            sums = []
            sum2s = []
            variances = []
            medians = []
            iqrs = []
            nuniqcolors = []
            dimensions = []
            stdDevs = []
            counts = []
            names = []
            for name in stuff:
                try:
                    if verbose:
                        print("Attempting to add " + name)
                    im = ImStats.ImStats(dirname + name, init=False)
                    means.append(im.stats["mean"])
                    rms.append(im.stats["rms"])
                    sums.append(im.stats["sum"])
                    sum2s.append(im.stats["sum^2"])
                    variances.append(im.stats["variance"])
                    medians.append(im.stats["median"])
                    iqrs.append(im.stats["iqr"])
                    nuniqcolors.append(im.stats["nuniquecolors"])
                    dimensions.append(im.stats["dimensions"])
                    stdDevs.append(im.stats["stdDev"])
                    counts.append(im.stats["count"])
                    names.append(name)
                    try:
                        cols = list('RGBA')
                        iqrcols = list('RGB')
                        self.storeframe("means", pd.DataFrame(means, columns=cols))
                        self.storeframe("rms", pd.DataFrame(rms, columns=cols))
                        self.storeframe("sums", pd.DataFrame(sums, columns=cols))
                        self.storeframe("sum2s", pd.DataFrame(sum2s, columns=cols))
                        self.storeframe("variances", pd.DataFrame(variances, columns=cols))
                        self.storeframe("medians", pd.DataFrame(medians, columns=cols))
                        self.storeframe("iqrs", pd.DataFrame(iqrs, columns=iqrcols))
                        self.storeframe("nuniqcolors", pd.DataFrame(nuniqcolors, columns=["Num Unique Colors"]))
                        self.storeframe("dimensions", pd.DataFrame(dimensions, columns=["Width", "Height"]))
                        self.storeframe("stddevs", pd.DataFrame(stdDevs, columns=cols))
                        self.storeframe("counts", pd.DataFrame(counts, columns=cols))
                        self.data["names"] = pd.DataFrame(names, columns=["Name"])
                    except AssertionError:
                        cols = list('RGB')
                        self.storeframe("means", pd.DataFrame(means, columns=cols))
                        self.storeframe("rms", pd.DataFrame(rms, columns=cols))
                        self.storeframe("sums", pd.DataFrame(sums, columns=cols))
                        self.storeframe("sum2s", pd.DataFrame(sum2s, columns=cols))
                        self.storeframe("variances", pd.DataFrame(variances, columns=cols))
                        self.storeframe("medians", pd.DataFrame(medians, columns=cols))
                        self.storeframe("iqrs", pd.DataFrame(iqrs, columns=cols))
                        self.storeframe("nuniqcolors", pd.DataFrame(nuniqcolors, columns=["Num Unique Colors"]))
                        self.storeframe("dimensions", pd.DataFrame(dimensions, columns=["Width", "Height"]))
                        self.storeframe("stddevs", pd.DataFrame(stdDevs, columns=cols))
                        self.storeframe("counts", pd.DataFrame(counts, columns=cols))
                        self.data["names"] = pd.DataFrame(names, columns=["Name"])
                    if verbose:
                        print("Added " + name + "\n")
                except IOError:
                    print(name + " wasn't a supported image type." + "\n")
            print("Done adding " + dirname + "\n")
        elif loadfile is True:
            self.load()

    def load(self):
        if self.loaded:
            raise AssertionError("Already loaded.")
        self.data = pd.HDFStore(self.dirname + self.filename)
        self.loaded = True

    def histograms(self, stat, verbose=False):
        if stat not in self.stats and stat is not "names":
            print("Please choose a valid statistic.\n")
        elif stat not in ["nuniqcolors", "dimensions"]:
            data = self.data[stat]
            rhist = Histogram(data, values='R',  color='Red', title="Red " + stat[0].upper() + stat[1:])
            ghist = Histogram(data, values='G', color='Green', title="Green " + stat[0].upper() + stat[1:])
            bhist = Histogram(data, values='B', color='Blue', title="Blue " + stat[0].upper() + stat[1:])
            if verbose:
                print(data)
            show(row(rhist, ghist, bhist))
        elif stat == "nuniqcolors":
            data = self.data[stat]
            hist = Histogram(data, values="Num Unique Colors", color='Purple', title="Number of Unique Colors")
            if verbose:
                print(data)
            show(hist)
        elif stat == "dimensions":
            data = self.data[stat]
            whist = Histogram(data, values="Width", color="Purple", title="Width Histogram")
            hhist = Histogram(data, values="Height", color='Purple', title="Height Histogram")
            if verbose:
                print(data)
            show(row(whist, hhist))

    def storeframe(self, stat, dataframe):
        self.data[stat] = dataframe

    def getstat(self, stat):
        if stat not in self.stats:
            raise AssertionError("Not a valid stat")
        if self.loaded is False:
            raise AssertionError("Self not loaded")
        return self.data[stat]

reason: from scipy


text:
# -*- coding: utf-8 -*-
from numpy import arange
from meshio import read
from os.path import join

from SciDataTool import DataTime, Data1D, VectorField

from ....Classes.SolutionData import SolutionData
from ....Classes.SolutionVector import SolutionVector
from ....Classes.MeshSolution import MeshSolution
from ....Classes.MeshVTK import MeshVTK

# TODO define "Results" dir, etc. in __init__
# TODO add groups, see get_meshsolution of MagFEMM
# TODO move reusable code to Elmer class

def get_meshsolution(self, output):
    """Get and set the mesh data and solution data.

    Parameters
    ----------
    self : StructElmer
        a StructElmer object
    save_path: str
        Full path to folder in which the results are saved

    Returns
    -------
    success: bool
        Information if meshsolution could be created

    """
    # logger
    logger = self.get_logger()

    # create meshsolution
    meshsol = MeshSolution(label="Elmer Structural")

    # if meshsoltion is not requested set empty MeshSolution
    if not self.is_get_mesh or not self.is_save_FEA:
        logger.info("StructElmer: MeshSolution is not stored by request.")
        output.struct.meshsolution = meshsol
        return False

    # get the mesh
    fea_path = self.get_path_save_fea(output)
    meshvtk = MeshVTK(path=join(fea_path, "Results"), name="case_t0001", format="vtu")
    # TODO maybe convert to MeshMat before
    meshsol.mesh = [meshvtk]

    # get the solution data on the mesh
    filename = join(fea_path, "Results", "case_t0001.vtu")
    meshsolvtu = read(filename)
    pt_data = meshsolvtu.point_data  # point_data is of type dict

    # setup axes
    indices = arange(meshsolvtu.points.shape[0])
    Indices = Data1D(name="indice", values=indices, is_components=True)

    # store only certain data
    store_dict = {
        "displacement": {
            "name": "Displacement",
            "unit": "mm",
            "symbol": "disp",
            "norm": 1e-3,
        },
        "vonmises": {
            "name": "Von Mises Stress",
            "unit": "MPa",
            "symbol": "vonmises",
            "norm": 1e6,
        },
    }
    comp_ext = ["x", "y", "z"]

    sol_list = []  # list of solutions

    for key, value in pt_data.items():
        # check if value should be stored
        if key in store_dict.keys():
            siz = value.shape[1]
            # only regard max. 3 components
            if siz > 3:
                logger.warning(
                    f'StructElmer get_meshsolution: size of data "{key}" > 3'
                    + " - "
                    + "Data will be truncated."
                )
                siz = 3

            components = []
            comp_name = []

            # loop though components
            for i in range(siz):
                # setup name, symbol and component name extension
                if siz == 1:
                    ext = ""
                else:
                    ext = comp_ext[i]

                # setup data object
                data = DataTime(
                    name=store_dict[key]["name"] + " " + ext,
                    unit=store_dict[key]["unit"],
                    symbol=store_dict[key]["symbol"] + ext,
                    axes=[Indices],
                    values=value[:, i],
                    normalizations={"ref": store_dict[key]["norm"]},
                )
                components.append(data)
                comp_name.append("comp_" + ext)

            # setup solution depending on number of field components
            if siz == 1:
                field = components[0]
                sol_list.append(
                    SolutionData(
                        field=field,
                        type_cell="point",
                        label=store_dict[key]["symbol"],
                    )
                )
            else:
                comps = {}
                for i in range(siz):
                    comps[comp_name[i]] = components[i]
                field = VectorField(
                    name=store_dict[key]["name"],
                    symbol=store_dict[key]["symbol"],
                    components=comps,
                )
                sol_list.append(
                    SolutionVector(
                        field=field,
                        type_cell="point",
                        label=store_dict[key]["symbol"],
                    )
                )

    meshsol.solution = sol_list

    output.struct.meshsolution = meshsol

    return True

max_stars_repo_name: ajpina/pyleecan

max_issues_repo_name: ajpina/pyleecan

max_forks_repo_name: ajpina/pyleecan

reason: from numpy


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

# # Gas Drag on Moons in Synestias
# Synestias can form moons. It is a natural consequence of the cooling process a synestia undergoes as the photosphere (where cloud cover is; where rock vapor condenses) radiates heat into the cold vacuum of space. There is a great deal of silicate vapor in the moon-forming region of a synestia, which provides a significant source of material with which to grow a moon.
#
# Before a moon (on the order of 1,000 km) can coalesce out of a synestia, the smaller building blocks of a moon, moonlets ($>$ dm-size liquid condensates), must grow from droplets. Droplets are mm- and cm- size rain that condense out of the cooled vapor at the photosphere. Similarly, before a puddle can pool on Earth's surface, water vapor must condense in Earth's atmosphere and rain to the surface.
#
# Once a parcel of gas from the synestia condenses into a liquid droplet, that parcel is no longer coupled to its fellow gas parcels. The surrounding vapor cannot support the condensate due to its higher density. The droplet loses pressure support from the surrounding gas and begins to fall towards the midplane (z = 0), pulled by gravitational forces. As the droplet falls, it passes through gas which imparts a drag on the droplet and causes the droplet to slow down and lose angular velocity. The orbit of the droplets will decay over time in accord with the strength of the gas drag acceleration.
#
# Will this droplet make it to the midplane and help form a moonlet? Will this droplet re-vaporize? Or will this droplet be lost to the planet-like central region of the synestia? To answer this, we need to know how strong the gas drag force is. What affects the strength of gas drag?
#
# ## Gas Drag Equation

# ![FigureGasDrag.png](Figures/FigureGasDrag.png)

# Caption. Gas drag acting on a sphere with a cross sectional area A of $\pi$r$_{body}^2$. The sphere has a mass m$_{body}$ and is moving with a velocity u$_{body}$ through a gas with density $\rho_{gas}$ moving at a velocity u$_{gas}$. As the sphere moves through the gas, gas particles will pile up in front of the sphere, creating a high pressure region. The gas drag that arises as a result of the velocity differential between the gas and sphere imparts an acceleration a$_{drag}$ on the sphere. Credit: G. O. Hollyday.

# The magnitude of the gas drag acceleration can be estimated by assuming that the moonlets can experience simple gas drag due to laminar flow around a solid sphere. The equation for simple gas drag acceleration is as follows:
#
# $$a_{drag}=\frac{C_D\ A \ \rho_{gas} \ (u_{gas} - u_{body})^2}{2 \ m_{body}}$$
#
# Simple gas drag acceleration a$_{drag}$ is dependent on:
# * the density $\rho_{gas}$ of the gas that drags a body of mass m$_{body}$,
# * the cross sectional area A of the body that the gas is interacting with (typically we assume solar system bodies are spherical, so they would have a circular cross section),
# * the gas drag coefficient C$_D$ (which depends on what thermal regime we are in), and
# * the relative velocity between the gas and the body being dragged u$_{gas}$-u$_{body}$.
#
# The relative velocity between the gas and a body is more important than how fast the body moves. If the body were traveling at the same speed and in the same direction as the gas surrounding it, it wouldn't matter how fast the body was traveling, there would be no gas drag acting on the body. The relative speed is what affects the magnitude of the gas drag the most. Think of the difference in velocity between a synestia's gas field and a moonlet as the headwind a moonlet must battle as it orbits within a synestia. If the gas and moonlet are moving in the same direction, but the gas is moving faster, the gas will help pull the moonlet along, much like how wind might push you from behind and help you run or move faster. If the gas and moonlet are moving in opposite directions, then the gas will slow the moonlet, much like how you might struggle to step or bicycle against a strong gust of a wind coming at you.
#
# Test this as you interact with the figure below, which plots the gas drag acceleration for a number of relative velocities. Notice how each of the parameters ("Size of moonlet", "Gas drag coefficient", and "Gas density") affects the magnitude of the gas drag.

# ```{margin} Running a code cell
# Access interactive features by 'Launch CoLab' or 'Launch Binder' from the rocket logo at the top of the page. When the interactive environment is ready, place your cursor in the code cell and press shift-return to execute the code. CoLab launches more quickly.
# ```
# Click the + symbol to see the code that generates the next interactive feature.

# In[1]:

import numpy as np
import math
from ipywidgets import *
import matplotlib.pyplot as plt

vel_rel = np.linspace(0.0, 1000., 100) #relative velocity between gas and body in m/s
dens = 3000. #lunar density in kg/m3; 3.0 g/cc

def dragfunc(rad_body, C_D, dens_gas):
    #mass of body in kg; assume spherical body (and spherical cross section)
    mass_body = dens*(4./3.)*math.pi*(rad_body**3)
    #acceleration imparted on body due to gas drag in m/s^2
    a_drag = 0.5*C_D*math.pi*(rad_body**2)*dens_gas*(vel_rel**2)/mass_body
    plt.plot(vel_rel, a_drag)
    plt.xlabel('velocity of the gas relative to the velocity of the moonlet (m/s)')
    plt.ylabel('gas drag acceleration (m/s$^2$)')
    plt.show()
    plt.close()

#RELATIVE SCALE IN GENERAL: how big the moonlets are (marble, golf ball, grain of sand,
#yoga ball, soccer field, etc.), how fast the relative and absolute velocity is
style = {'description_width': 'initial'}
layout = {'width': '400px'}
get_ipython().run_line_magic('matplotlib', 'inline')
interact(dragfunc,
         rad_body=FloatLogSlider(value=1, base=10, min=-3, max=6, step=1,
                                 description='Size of moonlet (m)', style=style, layout=layout),
         C_D=FloatSlider(value=2, min=0, max=2.5, step=0.5, description='Gas drag coefficient C$_D$ (-)',
                         continuous_update=False, readout=True, readout_format='.1f', style=style, layout=layout),
         dens_gas=FloatLogSlider(value=1, base=10, min=-3, max=3, step=1,
                                 description='Gas density (kg/m$^3$)', style=style, layout=layout)
        )

#acceleration in terms of what we experience on Earth
accel_label=['on board space shuttle ISS','elevator headed up','gravitational pull on surface of Earth',
             'roller coaster launch','space shuttle launch, gravitron amusement ride','top of roller coaster loop',
             'maximum heavy braking on Formula One Car','maximum acceleration permitted on piloted planes',
             'hard slap to the face, car crash','shock capacity of mechanical wrist watch',
             'Mantis Shrimp claw strike','proton in Large Hadron Collider']
accel_g=['1e-6 g','0.2 g','1 g','1.5 g','3 g','5 g','6.3 g','10 g','100 g','5000 g','10400 g','1.9e8 g']
left_box0=widgets.VBox([widgets.Label(str(i)) for i in accel_label])
right_box0=widgets.VBox([widgets.Label(str(i)) for i in accel_g])

#velocity in terms of what we experience on Earth
vel_label=['Walking','Biking','Usain Bolt sprinting','Train','Car moving 100 mph (160 km/h)',
           'Helicopter','Commerical airplane','Speed of sound in air (sound barrier)']
vel_ms=['1.5 m/s','4.5 m/s','12.5 m/s','27 m/s','45 m/s','70 m/s','220 m/s','343 m/s']
left_box1=widgets.VBox([widgets.Label(str(i)) for i in vel_label])
right_box1=widgets.VBox([widgets.Label(str(i)) for i in vel_ms])

widgets.HBox([left_box0,right_box0,left_box1,right_box1])

# Caption. Comparison of gas drag acceleration on various sized bodies for same surrounding gas. The gas drag acceleration is plotted in units of m/s$^2$ against the relative velocity in m/s.  For a more intuitive idea of how strong these accelerations are, or how fast these speeds are, there is also a list of different accelerations and speeds that you may encounter just below the plot above. The accelerations are given in terms of 'g' where 1 g = 9.8 m/s$^2$ (i.e. 1 g is the gravitational acceleration here on the surface of Earth). To the right of the accelerations is a second list of various speeds in m/s to calibrate against the relative velocity between the droplet at the gas.

# For the same gas field, gas drag acceleration is greater when the mass of the moonlet is smaller and the gas drag coefficient and relative velocity are larger. Gas drag will be stronger in a denser gas field.
#
# ## Decaying Orbits: Moonlets Dragged by Gas in a Synestia
# We would like to understand how gas drag affects how a single moonlet travels through the vapor in a synestia.
#
# In order to form a moon, moonlets must survive long enough to collide and grow into a body that is large enough to not be heavily affected by gas drag and will remain in orbit within the synestia without breaking apart from gravitational forces. Otherwise, the masses of the moonlets will be added to the center of a synestia (extends to about 10,000 km for the example Earth-mass synestia in the plot below) and lost from the moon-forming region of a synestia (about $>$ 20,000 km for the example Earth-mass synestia in the plot below). Moons are formed in the midplane, because they grow from colliding with other moonlets that reach the midplane. In the interactive below, you will see how the gas drag significantly decays and circularizes the orbit of a moonlet over time.
#
# Which scenarios result in the shortest decay times?
#
# In potential Moon-forming synestias, the gas drag coefficient is 2 and this has been set for you; changing the coefficient has a relatively minimal effect on the gas drag magnitude anyways. The velocity of the moonlet and the surrounding gas and the gas density gradient is calculated according to what is expected in a typical Moon-forming synestia.
#
# You will be able to choose where (in xyz) to initialize the moonlet. Where you initialize a moonlet indicates where that moonlet condensed from the vapor in a synestia. You can control the size of the moonlets, though it is more likely for smaller mass moonlets (roughly cm-sized) to form at the photosphere (z $>$ 10,000,000 m) and for moonlets in the midplane (z = 0) to be more massive (due to the addition of moonlets raining from the photosphere). Note that the moonlet's starting cylindrical radius will default to at least 20,000,000 m (starting x and y can't both be zero).
#
# Note: The interactive below is not run in real time (i.e. it doesn't update as you change the slider bars / drop-down menu). To start the interactive, click on the "Run Interact" button. You may have to be patient while the plots load. There is a blue progress bar located under the "Run Interact" button to indicate your place in the gas drag calculation. The gas drag calculation can be lengthy. If the notebook does not respond after 5 minutes, restart the notebook. If you are struggling to produce a plot in a reasonable amount of time, try these four scenarios:
# 1. 230 story building, x = 20,000,000 m, y = 20,000,000 m, and z = 0 m
# 2. 230 story building, x = 20,000,000 m, y = 20,000,000 m, and z = 0 m
# 3. The Moon, x = 20,000,000 m, y = 20,000,000 m, and z = 15,000,000 m
# 4. Cow, x = 20,000,000 m, y = 20,000,000 m, and z = 0 m

# ```{margin} Running a code cell
# Access interactive features by 'Launch Binder' or 'Launch CoLab' from the rocket logo at the top of the page. When the interactive environment is ready, place your cursor in the code cell and press shift-return to execute the code. If using CoLab (loads faster), you need to edit the code cell as directed to install rebound.
# ```
# Click the + symbol to see the code that generates the next interactive feature.

# In[2]:

#Dear Reader: If using CoLab, remove the # from the line below and shift-return to install rebound
#!pip install rebound

# In[3]:

# interactive feature
import math
import numpy as np
import rebound
import matplotlib.pyplot as plt
from IPython.display import clear_output

Mass_syn=5e24 #Earth mass synestia in kg
dens=3000. #density of moonlet kg/m3
G_hr=(6.674e-11)*(3600**2) #gravitational constant converted to m^3 kg^-1 hr^-2
const1=10.55#10.55#10 to 11
const2=0.86#0.86#0.85 to 0.9
const3=0.9e35#0.9e35 #1e33 to 1e37
const4=-4.7#-4.7 #-4.5 to -5

def pltdrag(r0,x0,y0,z0):
    clear_output(wait=True)
    f = FloatProgress(min=0, max=100) # instantiate the bar
    display(f) # display the bar
    m0=4.*math.pi*dens*(r0**3)/3. #mass of moonlet in kg
    s0=(x0**2 + y0**2)**0.5 #initial moonlet distance from rotational axis in m
    v_gas=0.95*(G_hr*Mass_syn/s0)**0.5 #sub-Keplerian gas velocity m/hr #0.01*(10.**(0.15599001870969695*np.log10(s0) - 1.6401471197424731))*((G_hr*Mass_syn/s0)**0.5)
    vx0=-v_gas*y0/s0 #x-component of sub-Keplerian gas velocity m/hr
    vy0=v_gas*x0/s0 #y-component of sub-Keplerian gas velocity m/hr

    #intialize orbit simulation
    sim_drag = rebound.Simulation() #start simulation
    sim_drag.units = ('Hr', 'M', 'Kg') #use SI units, seconds not available
    sim_drag.add(m=Mass_syn) #add central body: Earth mass synestia
    sim_drag.add(m=m0,r=r0) #add moonlet to simulation
    ps_drag = sim_drag.particles #initialize particle class (holds xyz position information)
    ps_drag[1].x=x0 #set initial x position of moonlet in m
    ps_drag[1].y=y0 #set initial y position of moonlet in m
    ps_drag[1].z=z0 #set initial z position of moonlet in m
    ps_drag[1].vx=vx0 #set initial x velocity of moonlet in m/hr #intial velocities will be same as gas
    ps_drag[1].vy=vy0 #set initial y velocity of moonlet in m/hr
    ps_drag[1].vz=0 #set initial z velocity of moonlet in m/hr #assume gas is not moving in z
    sim_drag.integrator = "ias15"#"whfast"#"leapfrog"#symplectic integrator #"ias15" #default integrator; adaptive time-stepping
    sim_drag.min_dt=0.001 #make sure doesn't pick too small of a time step
    sim_drag.dt=0.001 #start with this time step
    f.value += 5

    #define gas drag acceleration
    def drag(sim_drag):
        #drag coefficient is 2 in potential Moon-forming synestias
        ps_x = ps_drag[1].x
        ps_y = ps_drag[1].y
        ps_z = ps_drag[1].z
        ps_vx = ps_drag[1].vx
        ps_vy = ps_drag[1].vy
        ps_vz = ps_drag[1].vz
        s=np.sqrt(ps_x**2 + ps_y**2) #cylindrical radius of moonlet in m
        z_s=const1*(s**const2) #scale height fit in m
        rho_g=const3*(s**const4)*np.exp(-(ps_z/z_s)**2) #gas density fit in kg/m^3
        v_g=0.95*np.sqrt(G_hr*Mass_syn/s) #sub-Keplerian gas velocity in m/hr #(10.**(0.15599001870969695*np.log10(s) - 1.6401471197424731))*((reb_sim.G*Mass_syn/s)**0.5)
        theta=np.arctan(ps_y/ps_x)
        v_gx=-np.abs(v_g*np.sin(theta))*np.sign(ps_y) #x component of gas velocity in m/hr
        v_gy=np.abs(v_g*np.cos(theta))*np.sign(ps_x) #y component of gas velocity in m/hr
        #assumption: speed of gas in z, v_gz, will be zero for synestia in hydrostatic equilibrium
        v_rel_x=ps_vx - v_gx
        v_rel_y=ps_vy - v_gy
        v_rel=np.sqrt(v_rel_x**2 + v_rel_y**2 + ps_vz**2)
        factor = math.pi*(ps_drag[1].r**2)*rho_g*v_rel/ps_drag[1].m #common factor for gas drag in components
        ps_drag[1].ax += -factor*v_rel_x
        ps_drag[1].ay += -factor*v_rel_y
        ps_drag[1].az += -factor*ps_vz

    #tell simulation to add gas drag and update velocities with each time step
    sim_drag.additional_forces = drag
    sim_drag.force_is_velocity_dependent = 1
    sim_drag.move_to_com() #move to center of mass (central body is at origin)

    Noutputs = 2000 #increase this to get a smoother curve for orbit
    if r0 < 10:
        times = np.linspace(0.,600., Noutputs)
    elif r0 > 1e4:
        times = np.linspace(0.,600., Noutputs)
    else:
        times = np.linspace(0.,200., Noutputs) #times to output position information of particle

    #intialize holder arrays for positional output from simulation
    xholder = np.zeros(Noutputs)
    yholder = np.zeros(Noutputs)
    zholder = np.zeros(Noutputs)
    sholder = np.zeros(Noutputs)
    f.value += 5
    fstep = 81./len(times)
    for k,time in enumerate(times):
        sim_drag.integrate(time)
        xtemp = ps_drag[1].x
        ytemp = ps_drag[1].y
        stemp = np.sqrt(xtemp**2 + ytemp**2)
        f.value += fstep
        if stemp <= 10e6:
            #exit out of simulation when reach the inner region (10,000,000 m)
            #shorten arrays to exclude zeros at end
            xholder = xholder[:k]
            yholder = yholder[:k]
            zholder = zholder[:k]
            sholder = sholder[:k]
            times = times[:k]
            f.value = 91
            print('The length of time the moonlet spends in the disk-like region is {0:.1f} hr.'.format(times[-1]))
            break
        else:
            xholder[k] = xtemp  # This stores the data which allows us to plot it later
            yholder[k] = ytemp
            zholder[k] = ps_drag[1].z
            sholder[k] = stemp
    if sholder[-1] > 11e6:
        print('The moonlet remains in the disk-like region at {0:.1f} hr.'.format(times[-1]))
    xholder = xholder*(1e-3)
    yholder = yholder*(1e-3)
    zholder = zholder*(1e-3)
    sholder = sholder*(1e-3)
    rholder = np.sqrt(sholder**2 + zholder**2)
    ymax=np.amax(rholder)

    #plot
    fig = plt.figure(figsize=(18,4))
    ax = plt.subplot(131)
    plt.plot(xholder, yholder)
    plt.plot(x0/1e3,y0/1e3,'x',markersize=10)
    plt.plot(0,0,'k*')
    plt.plot(xholder[-1], yholder[-1],'o')
    plt.xlabel('x (km)')
    plt.ylabel('y (km)')
    plt.ylim([-ymax,ymax])
    plt.xlim([-ymax,ymax])
    ax.set_aspect(1)

    f.value += 3

    ax2 = plt.subplot(132)
    plt.plot(sholder, zholder)
    plt.plot(s0/1e3,z0/1e3,'x',markersize=10)
    plt.plot(0,0,'k*')
    plt.plot(sholder[-1], zholder[-1],'o')
    plt.xlabel('distance to rotational axis r$_{xy}$ (km)')
    plt.ylabel('vertical distance off midplane z (km)')
    plt.ylim([-ymax,ymax])
    plt.xlim(xmin=0)
    ax2.axis('equal')

    f.value += 3

    ax3 = plt.subplot(133)
    plt.plot(times,rholder)
    plt.plot(0,rholder[0],'x',markersize=10)
    plt.plot(times[-1],rholder[-1],'o')
    plt.xlabel('time (hr)')
    plt.ylabel('radius r (km)')
    plt.show()
    plt.close()

    f.value += 3

x0slider = FloatSlider(value=20e6, min=0, max=60e6, step=20e6, description='x starting position of moon (m)',
                         continuous_update=False, readout=True, readout_format='.0f', style=style, layout=layout)
y0slider = FloatSlider(value=20e6, min=0, max=60e6, step=20e6, description='y starting position of moon (m)',
                         continuous_update=False, readout=True, readout_format='.0f', style=style, layout=layout)

def rxy_min_x0(*args):
    rxy0 = np.sqrt(x0slider.value**2 + y0slider.value**2)
    if x0slider.value < 20e6:
        y0slider.value = 20e6

def rxy_min_y0(*args):
    rxy0 = np.sqrt(x0slider.value**2 + y0slider.value**2)
    if y0slider.value < 20e6:
        x0slider.value = 20e6

x0slider.observe(rxy_min_x0, 'value')
y0slider.observe(rxy_min_y0, 'value')
interact_manual(pltdrag,
         r0=widgets.Dropdown(
             options=[('Coarse sand (mm)', 1.e-3),
                      ('Diameter of a belly button (cm)', 1.e-2),
                      ('Coconut (dm)', 0.1),
                      ('Cow (m)', 1.),
                      ('230 story building (km)', 1.e3),
                      ('Width of the United States (Mm)', 1.e6),
                      ('The Moon (3.5 Mm)', 3.5e6)],
             value=1.e3, description='Size (diameter) of moonlet:', continuous_update=False, style=style),
         x0=x0slider,
         y0=y0slider,
         z0=FloatSlider(value=0, min=0, max=30e6, step=5e6, description='z starting position of moon (m)',
                         continuous_update=False, readout=True, readout_format='.1f', style=style, layout=layout),
    )

# Caption. The effect of gas drag acceleration on the orbit of various sized spherical bodies for the same surrounding varying-density gas. The simulation ends when, during any point in the condensate's orbit, the condensate's distance from the rotational axis reaches the planet-like region, which is r = 10,000 km in this case. The time it takes the moonlet to reach this distance is called the decay time. For bodies initialized off the midplane, the particle oscillates a bit around the midplane as it settles into gravitational equilibrium. The particle overshoots the midplane a little initially because it cannot slow down immediately at the midplane, loses speed as it turns back around but gains a little speed as it gets accelerated back up towards the midplane, and so on, until there is no movement in z.

# You may have noticed in the interactive above that a condensate decays more quickly when:
# 1. the gas field of a synestia is denser,
# 2. a condensate forms in the midplane (z = 0) as opposed to in the photosphere (high z),
# 3. and a condensate accretes closer to the center of a synestia.
#
# Why is this?
#
# ### 1. Gas Drag Acceleration Increases with Gas Field Density
#
# The balance of forces changes for a condensate:
#
# $$\omega^2r = \frac{GM}{r^2}\left(1 + \frac{3J_2a_{eq}^2}{2r^2}\left(\frac{5z^2}{r^2} - 1\right)\right) + a_{drag}$$
#
# where the gas slows down the condensate instead of supporting it. The denser the gas field surrounding the condensate, the stronger the gas drag acceleration that pulls on the condensate and the quicker the angular velocity of the condensate decreases. As a consequence, the angular velocity deficit between the condensate and the equivalent Keplerian orbit is greater, and the condensate's orbit decays more rapidly.
#
# ### 2. Gravity Off the Midplane Tempers Gas Drag Acceleration
#
# Net accelerations imparted on a droplet by gas drag and gravity continuously change as the droplet travels in a synestia.
#
# When a droplet condenses out of a synestia near the photosphere, it loses pressure support, but will initially have the same velocity as the gas around it. The gas is moving at just below Keplerian velocities in r$_{xy}$ = (x, y) with zero velocity in z. At the start, the moonlet experiences very little gas drag and is being pulled towards the midplane (z = 0) by gravity, no longer balanced by gas pressure in z.
#
# As gravity accelerates the moonlet towards the midplane, this creates a large velocity differential between the gas (zero velocity in z) and the moonlet. As gas drag begins to act in z, the particle slows in its descent to the midplane. Simultaneously, gravity in r$_{xy}$ will move the droplet to a smaller r$_{xy}$ on a Keplerian orbit. However, gas drag will slow this inward motion.
#
# Gas drag perturbs the moonlet's orbit, resulting in an elliptical, or eccentric, orbit. This temporarily prevents inward spiraling of the particle since the elliptical orbit moves the moonlet away from the center of its synestia. The higher off the midplane, the more elliptical the orbit. The farther the elliptical orbit moves the droplet, the slower the inward decay of the droplet.
#
# ### 3. Greater Distance from Rotational Axis Allows More Orbital Travel Time
#
# The farther away from the rotational axis the particle condenses and begins its orbital decay, the slower the particle decays.
#
# As a condensate decays due to gas drag, the condensate will approach the central body closely on one side of the orbit (nearest approach). The location of the condensate's nearest approach will travel closer to the central body with time. If a moonlet travels close enough to the central body, it will become a part of the interior (planet-like center) of a synestia.
#
# The greater the distance between the condensate and the planet-like region, the greater the distance the condensate can travel before entering the planet-like region of a synestia. This is a given, but it takes a longer time to travel a larger distance. If the condensate initializes at a greater distance from the center of a synestia, then the longer it will take for the condensate to spiral inwards and decay.
#
# ### Orbital Decay Time Dependence on Mass
#
# Gas drag has a smaller effect on large particles (it is difficult to change their motion due to their large mass) and small particles (very little surface area for the gas to exert drag on). Very small droplets and very large moonlets survive the longest in the interior of a synestia, whereas medium-sized particles decay the quickest. Droplets and moonlets that decay more rapidly than they can grow will spiral into the planet-like center of a synestia.
#
# Based on the fastest accretion time for the Moon in cooling simulations of Earth-mass synestias, the minimum time in which a Moon-mass could grow is about 10 years (87,600 hrs) ([Lock et al., 2018](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2017JE005333)). Compared to this, for orbits (calculated by the interactive above) that begin in the midplane close to the central body, a single 1-km-sized moonlet could survive 10's of hours, which is the fastest decay time of all the moonlet and droplet sizes. In contrast, the quickest decay times (calculated by the interactive plot above) for a meter-sized moonlet and a Moon-sized moonlet are at least several orders of magnitude longer.
#
# The quick decay of km-sized moonlets is a barrier for lunar growth within a synestia. As smaller moonlets reach km-scale, they must grow quickly to enter a size regime where the gas won't decay the moonlet's orbit as quickly.
#
# #### Continually Growing Masses
#
# It is important to note that a single moonlet within a synestia is not alone. There are many droplets condensing out of the photosphere at once, which creates a nearly continuous torrential downpour of rock rain as the synestia cools. These droplets will collide with one another and grow, increasing their mass. If the droplets grow large enough, they become moonlets. As the moonlets fall to the midplane, they also accrete surrounding droplets and moonlets. Once droplets and moonlets reach the midplane, they can accrete to other material in the midplane, and therefore slow their decay. How quickly the droplets and moonlets combine into larger masses relative to the decay of their orbits via gas drag is what will dictate whether, ultimately, the moon could survive.
#
# #### Survival of the Moon Depends on the Growth of the Lunar Seed
#
# In the stages just after the giant impact that forms a synestia, material in the midplane either quickly binds together via gravitational clumping or was part of the impacting bodies and didn't fully melt or vaporize, and acts as a lunar seed for droplets and moonlets to hold onto as they spiral past. The lunar seed acts as the core of the moon and scoops up droplets and moonlets within reach (under the moon's gravitational influence). The lunar seed must survive until the gas in a synestia disappears. The larger the mass of the lunar seed, the slower its orbit will decay due to gas drag. Orbital decay of droplets and mooonlets can help grow the moon (lunar seed), depending on where a given droplet or moonlet is in relation to the lunar seed.
#
# #### Comparison with Canonical Model
#
# In a canonical (planet-disk) Moon-forming case, the characteristic moonlet size that forms via spontaneous gravitational clumping of material in the disk is 10 km ([Salmon & Canup, 2012](https://iopscience.iop.org/article/10.1088/0004-637X/760/1/83/meta)). In a gas regime such as the one seen in the interactive above, these 10-km clumps would not last very long in a synestia (10's of hours). However, according to moon-forming simulations of the canonical case, these same sized moonlets would survive despite their short gas drag decay timescales, because they accrete quickly onto the lunar seed as they spiral in towards the planet ([Thompson & Stevenson, 1988](https://authors.library.caltech.edu/39035/)). Without a lunar seed, moonlets would not survive very long.
#
# ### Vaporization of Droplets May Slow Their Decay But Decrease Mass that Feeds Lunar Seed
#
# There is another factor as well. As a droplet descends into a synestia, the temperature and pressure of the gas surrounding the droplet will increase, with maximum temperature and pressure values at the midplane. A droplet must combat vaporization in order to provide liquid material for the moon to grow from. The good news is that if a droplet vaporizes, its decay ceases as it buoys up towards the photosphere and remains in the moon-forming region of a synestia. This is a bit more complicated considering there are interactions with other droplets occurring, which may involve thermal buffering (some droplets at lower z act as shields and provide temperature buffers for the droplets behind them at greater z). However, if the mass of condensed material reaching the midplane is less due to vaporization of droplets, then less material accretes to the lunar seed at any given time. The lunar seed may decay more quickly if this is true. I will not explore the effects of vaporization here, but keep it in mind as you evaluate the decay times in the interactive figure above.

# ## References
#
# Lock, S. J., Stewart, S. T., Petaev, M. I., Leinhardt, Z. M., Mace, M. T., Jacobsen, S. B., & Ćuk, M. (2018). The origin of the Moon within a terrestrial synestia. Journal of Geophysical Research: Planets (American Geophysical Union), 123 (4), 910-951.
#
# Salmon, J., & Canup, R. M. (2012). Lunar accretion from a Roche-interior fluid disk. The Astrophysical Journal, 760 (83), 1-18.
#
# Thompson, C., & Stevenson, D. J. (1988). Gravitational instability in two-phase disks and the origin of the Moon. The Astrophysical Journal, 333 (1), 452-481.


text:
# Under MIT licence, see LICENCE.txt

import time
from typing import Optional, List

import numpy as np

from Util import Pose
from Util.ai_command import Idle, CmdBuilder, MoveTo
from Util.constant import BALL_RADIUS, ROBOT_RADIUS, KickForce
from ai.GameDomainObjects import Player
from ai.STA.Action.GoBehind import GoBehind
from ai.STA.Tactic.tactic import Tactic
from ai.STA.Tactic.tactic_constants import Flags
from ai.states.game_state import GameState

ORIENTATION_DEADZONE = 0.2
DISTANCE_TO_KICK_REAL = ROBOT_RADIUS * 3.4
DISTANCE_TO_KICK_SIM = ROBOT_RADIUS + BALL_RADIUS
COMMAND_DELAY = 1.5

class PassToPlayer(Tactic):
    def __init__(self, game_state: GameState, player: Player, target: Pose=Pose(), args: Optional[List[str]]=None):
        super().__init__(game_state, player, target, args)
        self.current_state = self.kick_charge
        self.next_state = self.kick_charge
        self.last_time = time.time()
        if args:
            self.target_id = int(args[0])
        else:
            self.target_id = 1

    def kick_charge(self):
        if time.time() - self.last_time > COMMAND_DELAY:
            self.next_state = self.get_behind_ball
            self.last_time = time.time()

        return CmdBuilder().addKick().addForceDribbler().build()

    def get_behind_ball(self):
        self.status_flag = Flags.WIP
        player = self.player.pose.position
        ball = self.game_state.ball_position
        target = self.game_state.get_player_position(self.target_id)

        vector_player_2_ball = ball - player
        vector_player_2_ball /= vector_player_2_ball.norm

        if self._is_player_towards_ball_and_target():
            self.next_state = self.grab_ball
        else:
            self.next_state = self.get_behind_ball
        return GoBehind(self.player, self.game_state.ball_position, target, 120)

    def grab_ball(self):
        if self._get_distance_from_ball() < 120:
            self.next_state = self.kick
            self.last_time = time.time()
        elif self._is_player_towards_ball_and_target(-0.9):
            self.next_state = self.grab_ball
        else:
            self.next_state = self.get_behind_ball
        return MoveTo(Pose(self.game_state.ball.position, self.player.pose.orientation))

    def kick(self):
        if self._get_distance_from_ball() > 300:
            self.next_state = self.halt
            self.last_time = time.time()
        elif time.time() - self.last_time < COMMAND_DELAY:
            self.next_state = self.kick
        else:
            self.next_state = self.kick_charge
        return CmdBuilder().addKick(KickForce.HIGH).addMoveTo(self.target).build()

    def halt(self):
        if self.status_flag == Flags.INIT:
            self.next_state = self.kick_charge
        else:
            self.status_flag = Flags.SUCCESS
        return Idle

    def _get_distance_from_ball(self):
        return (self.player.pose.position - self.game_state.ball_position).norm

    def _is_player_towards_ball_and_target(self, fact=-0.99):
        player = self.player.pose.position
        ball = self.game_state.ball_position
        target = self.game_state.get_player_position(self.target_id)

        vector_player_2_ball = ball - player
        vector_target_2_ball = ball - target
        vector_player_2_ball /= vector_player_2_ball.norm
        vector_target_2_ball /= vector_target_2_ball.norm
        vector_player_dir = np.array([np.cos(self.player.pose.orientation),
                                      np.sin(self.player.pose.orientation)])
        if np.dot(vector_player_2_ball.array, vector_target_2_ball.array) < fact:
            if np.dot(vector_player_dir, vector_target_2_ball.array) < fact:
                return True
        return False

max_stars_repo_name: RoboCupULaval/StrategyAI

max_issues_repo_name: RoboCupULaval/StrategyIA

max_forks_repo_name: RoboCupULaval/StrategyIA


text:
import tensorflow as tf
import numpy as np
import cv2
import math
from dataclasses import dataclass
import time

@dataclass
class identified_object:
    x: int
    y: int
    radius: int
    area:int

@dataclass
class circle_to_draw:
    x: int
    y: int
    radius: int

# Lets load the model
model = tf.keras.models.load_model('spikeball-nn.h5')

# Lets load the video

cap = cv2.VideoCapture("spikeball.mov")

fg = cv2.createBackgroundSubtractorMOG2(history=25, detectShadows=False)

fourcc = cv2.VideoWriter_fourcc('m', 'p', '4', 'v')
writer = cv2.VideoWriter("spikeball-with-ball.mp4", fourcc, cap.get(cv2.CAP_PROP_FPS), (1920, 1080))

ret, frame = cap.read()

fno = 1

ball_list: list[circle_to_draw] = [None] * 30

while ret:
    ret, frame = cap.read()
    if ret:
        prev = time.perf_counter()
        frame_copy = frame.copy()
        # Lets do some processing
        fgmask = fg.apply(frame)
        fgmask = cv2.erode(fgmask, None, iterations=1)
        fgmask = cv2.boxFilter(fgmask, -1, (5, 5))
        # Now lets take that fg mask and grab its contours
        contours, hierarchy = cv2.findContours(fgmask, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
        current = time.perf_counter()
        print(f"Took {current - prev} to apply mask and find contours")
        prev = time.perf_counter()
        evaluated_contours = 0
        added_ball = False
        for c in contours:
            circle = cv2.minEnclosingCircle(c)
            area = cv2.contourArea(c)
            center = (int(circle[0][0]), int(circle[0][1]))
            x = center[0]
            y = center[1]
            radius = int(circle[1])
            circle_area = math.pi * radius * radius
            # lets make sure the image can be grabed due to range limitaitons:
            if 25 < x and x < 1895 and 25 < y and y < 1055:
                # Now lets check to see if the radius meets spec
                # if not 2 < radius < 20:
                #     cv2.circle(frame, center=center, radius=radius, color=(255, 0, 0), thickness=3)
                # if not area * .5 < circle_area < area * 1.5:
                #     cv2.circle(frame, center=center, radius=radius, color=(255, 255, 255), thickness=3)
                # if not 150 < area < 1000 and not circle_area * .5 < area < circle_area * 1.5:
                #     cv2.circle(frame, center=center, radius=radius, color=(255, 255, 255), thickness=3)
                if 150 < area < 1000 and circle_area * .5 < area < circle_area * 1.5:
                    evaluated_contours += 1
                    #item = identified_object(x, y, radius, area)
                    # Now lets grab the image
                    image = frame_copy[np.ix_(range(y - 25, y + 25),range(x - 25, x + 25))][:,:,::-1]
                    # Now lets pipe it into the network
                    img_array = tf.expand_dims(image, axis=0) # Create a batch
                    nn_start = time.perf_counter()
                    y_proba = model(img_array)
                    # y_proba = model.predict(img_array)
                    nn_end = time.perf_counter()
                    #print(f"Took {nn_end - nn_start} to run the network")
                    score = tf.nn.softmax(y_proba[0])
                    #print(score)
                    if score[1] > 0.8:
                        #cv2.circle(frame, center=center, radius=radius, color=(0, 255, 0), thickness=-1)
                        ball_list.append(circle_to_draw(x, y, radius))
                        added_ball = True
                        break
                    # else:
                    #     cv2.circle(frame, center=center, radius=radius, color=(0, 0, 255), thickness=3)
        if not added_ball:
            ball_list.append(None)
        for ball in ball_list[:-1]:
            if ball is not None:
                #print(ball.x, ball.y, ball.radius)
                cv2.circle(frame, center=(ball.x, ball.y), radius=(ball.radius  // 2), color=(0, 255, 0), thickness=-1)
        current = time.perf_counter()
        print(f"Took {current - prev} to iterate over {len(contours)} contours and eval {evaluated_contours} of them")
        writer.write(frame)
        print(f"Wrote frame {fno}")
        fno += 1
        ball_list.pop(0)

writer.release()


text:
import numpy as np

from anomalydetection_app.inspect_element import is_selected_from_n_clicks

rng = np.random.default_rng()

def sin_data(x):
    new_xs = x / len(x) * 2 * np.pi * 2
    data = np.sin(new_xs)
    return data

def sin_cos_data(x):
    new_xs = x / len(x) * 2 * np.pi * 8
    sin_cos = np.sin(new_xs) + np.cos(new_xs / 2)
    return sin_cos

def linear_data(x):
    linear = x / 1000
    return linear

def normal_noise_per_time_point(noise_factor, this_x, time_point_noise_probability):
    noise = add_normal_noise_to_series(np.zeros(len(this_x)), noise_factor)
    noise_locations = rng.uniform(size=len(this_x)) <= time_point_noise_probability
    noise = noise * noise_locations
    return noise

def exponentially_distributed_noise(this_x, time_point_noise_probability):
    noise = np.zeros(len(this_x))
    if time_point_noise_probability > 0:
        anomaly_locations = get_anomaly_locations(this_x, time_point_noise_probability)
        noise[anomaly_locations] = 1
    return noise

def exponentially_distributed_cluster_noise(this_x, time_point_noise_probability):
    cluster_size = 5
    time_point_noise_probability = time_point_noise_probability / cluster_size
    noise = np.zeros(len(this_x))
    if time_point_noise_probability > 0:
        anomaly_locations = get_anomaly_locations(this_x, time_point_noise_probability)
        final_locations = construct_clustered_anomalies(this_x, anomaly_locations,
                                                        half_cluster=int(np.floor(cluster_size / 2)))
        noise[final_locations] = 1
    return noise

def get_anomaly_locations(x, time_point_noise_probability):
    n_expected_events = int(len(x) * time_point_noise_probability)
    time_to_anomaly = rng.exponential(1 / time_point_noise_probability, size=n_expected_events)
    anomaly_locations = [int(item) for item in np.cumsum(time_to_anomaly) if int(item) < len(x)]
    return anomaly_locations

def construct_clustered_anomalies(x, anomaly_locations, half_cluster):
    final_locations = []
    for location in anomaly_locations:
        append_locations_in_cluster(final_locations, half_cluster, location)
    final_locations = np.array(final_locations, dtype=int)
    final_locations = final_locations[(final_locations >= 0) & (final_locations < len(x))]
    return final_locations

def append_locations_in_cluster(final_locations, half_cluster, location):
    for i in range(-half_cluster, half_cluster, 1):
        final_locations.append(int(location) + i)

def add_normal_noise_to_series(series, noise_factor):
    len_series = len(series)
    return series + rng.normal(scale=noise_factor, size=len_series)

def add_noise_to_data(data, exp_cluster_noise_clicks, exp_noise_clicks, noise_factor, normal_noise_clicks,
                      time_point_noise_probability):
    if is_selected_from_n_clicks(exp_noise_clicks):
        data = data + exponentially_distributed_noise(data, time_point_noise_probability)
    if is_selected_from_n_clicks(exp_cluster_noise_clicks):
        data = data + exponentially_distributed_cluster_noise(data, time_point_noise_probability)
    if is_selected_from_n_clicks(normal_noise_clicks):
        data = data + normal_noise_per_time_point(noise_factor, data, time_point_noise_probability)
    return data

def simulate_data_pattern(xs, linear_clicks, sin_clicks, sin_cos_clicks):
    data = np.zeros(shape=len(xs))
    if is_selected_from_n_clicks(sin_clicks):
        data = data + sin_data(xs)
    if is_selected_from_n_clicks(sin_cos_clicks):
        data = data + sin_cos_data(xs)
    if is_selected_from_n_clicks(linear_clicks):
        data = data + linear_data(xs)
    return data