text:
import imutils
import cv2
import numpy as np
import os
import pandas as pd
from PIL import Image, ImageStat, ImageEnhance
import util
from datetime import datetime

out_dir_path = './out'
small_plaques = False
debug_mode = False

def debug_info(file_path, img):
    if debug_mode:
        cv2.imwrite(file_path, img)

def adjust_gamma(image, gamma=1.0):
    # build a lookup table mapping the pixel values [0, 255] to
    # their adjusted gamma values
    inv_gamma = 1.0 / gamma
    table = np.array([((i / 255.0) ** inv_gamma) * 255 for i in np.arange(0, 256)]).astype("uint8")

    return cv2.LUT(image, table)

def process_image(image, contrast):
    gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
    debug_info("./test_pic_grey.jpg", gray)

    # different values for different plaque sizes
    h = 6
    if small_plaques:
        h = 3
    gray = cv2.fastNlMeansDenoising(gray, h=h)
    debug_info("./test_pic_grey_thresh_denoise.jpg", gray)

    # gray = unsharp_mask(gray)
    # cv2.imwrite("./test_pic_grey_unsharp.jpg", gray)

    blurred = cv2.GaussianBlur(gray, (7, 7), 0)
    # blurred = cv2.medianBlur(gray, 9)
    debug_info("./test_pic_blur.jpg", blurred)

    # ET point where it depends on a pic
    high_contrast = cv2.convertScaleAbs(blurred, alpha=contrast, beta=0)
    debug_info("./test_pic_high.jpg", high_contrast)

    gamma_test = adjust_gamma(high_contrast, 7.1)
    debug_info("./test_pic_green_gamma_0.jpg", gamma_test)

    high_contrast = adjust_gamma(high_contrast, 1.0)
    debug_info("./test_pic_green_gamma.jpg", high_contrast)

    # binary = cv2.adaptiveThreshold(high_contrast, 500, cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY, 55, 2)
    # binary = cv2.adaptiveThreshold(high_contrast, 500, cv2.ADAPTIVE_THRESH_GAUSSIAN_C , cv2.THRESH_BINARY, 55, 2)

    plate_only = cv2.adaptiveThreshold(high_contrast, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY_INV, 9, 1)
    # plate_only = cv2.threshold(high_contrast, 65, 255, cv2.THRESH_BINARY_INV)
    debug_info("./test_pic_plate_only.jpg", plate_only)

    # ret, thresh = cv2.threshold(high_contrast, 162, 255, cv2.THRESH_BINARY_INV)

    # ET added
    # blockSize affects large/small plaques (circles in circles)
    # change blockSize based on the AREA
    #block_size = 265
    block_size = 231
    if small_plaques:
        block_size = 49

    thresh = cv2.adaptiveThreshold(high_contrast, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY_INV,
                                   block_size, 2)
    # thresh = cv2.adaptiveThreshold(high_contrast, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY_INV, 49, 2)
    debug_info("./test_pic_grey_adapt_thresh.jpg", thresh)

    # laplacian = cv2.Laplacian(blur, -1, ksize=17, delta=-50)
    # laplacian = cv2.Laplacian(thresh, cv2.CV_64F)
    laplacian = cv2.Laplacian(thresh, cv2.CV_8UC1)
    debug_info("./test_pic_laplacian.jpg", laplacian)
    # gray_lapl = cv2.cvtColor(laplacian, cv2.COLOR_BGR2GRAY)

    # binary = cv2.threshold(laplacian, 165, 255, cv2.THRESH_BINARY)
    # binary = cv2.adaptiveThreshold(laplacian, 500, cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY, 55, 2)
    # cv2.imwrite("./test_pic_green_laplacian_binary.jpg", binary)

    clr_high_contrast = cv2.cvtColor(high_contrast, cv2.COLOR_GRAY2BGR)
    return laplacian, high_contrast, clr_high_contrast
    # return binary, high_contrast, clr_high_contrast

def get_contours(binary_image):
    contours = cv2.findContours(binary_image, cv2.RETR_CCOMP, cv2.CHAIN_APPROX_SIMPLE)
    #contours = cv2.findContours(binary_image.copy(), cv2.RETR_CCOMP, cv2.CHAIN_APPROX_SIMPLE)
    #contours = cv2.findContours(binary_image.copy(), cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    #contours = cv2.findContours(binary_image.copy(), cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)

    return imutils.grab_contours(contours)

def draw_contours(image, green_df, red_df, other_df, plate_df):
    image_copy = image.copy()
    for index, green in green_df.iterrows():
        draw_one_contour(image_copy, green, (0, 255, 0))
        if debug_mode:
            for index, red in red_df.iterrows():
                draw_one_contour(image_copy, red, (0, 0, 255))
            for index, other in other_df.iterrows():
                # draw_one_contour(image_copy, other, (200, 150, 150))
                draw_one_contour(image_copy, other, (100, 100, 100))
            for index, plate in plate_df.iterrows():
                draw_one_contour(image_copy, plate, (0, 128, 255))

    return image_copy

def draw_one_contour(image, c_df, color):
    m = cv2.moments(c_df['CONTOURS'])
    if m["m00"] != 0:
        cx = int((m["m10"] / m["m00"]))
        cy = int((m["m01"] / m["m00"]))
    else:
        cx = 0
        cy = 0

    pd.set_option('display.precision', 2)
    image_w_contours = cv2.drawContours(image, [c_df['HULL']], -1, color, 1)
    # cv2.putText(image, f"#{c_df['INDEX_COL']}:{c_df['ENCL_DIAMETER_MM']}", (cx, cy), cv2.FONT_HERSHEY_SIMPLEX, 0.5,
    #             color, 1)

    if c_df['DIAMETER_MM'] == 0:
        cv2.putText(image, f"#{c_df['INDEX_COL']}:{c_df['DIAMETER_PXL']}", (cx, cy), cv2.FONT_HERSHEY_SIMPLEX, 0.5,
                    (25, 51, 0), 1)
    else:
        cv2.putText(image, f"#{c_df['INDEX_COL']}:{c_df['DIAMETER_MM']}", (cx, cy), cv2.FONT_HERSHEY_SIMPLEX, 0.5,
                    (25,51,0), 1)
    return image_w_contours

def get_image_paths(image, directory):
    images = []
    if image:
        images.append(image)
    else:
        for r, d, f in os.walk(directory):
            for file in f:
                if os.path.splitext(file)[1] in ['.tif', '.tiff', '.jpg', '.jpeg', '.png']:
                    images.append(os.path.join(r, file))
    return images

def write_images(out_dir, output_image, binary_image, high_contrast_image, image_path):
    now = datetime.now() # current date and time
    fileTime = now.strftime("%Y-%m-%d_%H-%M-%S")

    syncPath = f'{out_dir}/out_{fileTime}_{os.path.split(image_path)[1]}'

    cv2.imwrite(syncPath, output_image)
    # cv2.imwrite(f'{out_dir}/out_red-{os.path.split(image_path)[1]}', output_image_red)
    #cv2.imwrite(f'{out_dir}/contrast-{os.path.split(image_path)[1]}', high_contrast_image)
    #cv2.imwrite(f'{out_dir}/binary-{os.path.split(image_path)[1]}', binary_image)
    return syncPath

def write_data(out_dir, image_path, green_df, red_df, other_df, syncPath = 'out/out_0_0.jpg'):

    dataPath = write_one_data(out_dir, image_path, syncPath, green_df)
    return dataPath
    #write_one_data(out_dir, image_path, 'green', green_df)
    #write_one_data(out_dir, image_path, 'red', red_df)
    #write_one_data(out_dir, image_path, 'other', other_df)

def write_one_data(out_dir, image_path, syncPath, df):

    #image_file_name = os.path.split(image_path)[1]
    image_file_name = os.path.split(syncPath)[1]

    image_name = os.path.splitext(image_file_name)[0]

    path = ''
    if df.empty != True:
        #path = f'{out_dir}/data-{prefix}-{image_name}.csv'
        path = f'{out_dir}/{image_name}.csv'

        df.to_csv(path_or_buf=path, index = False,
                  #columns=['INDEX_COL', 'AREA_PXL', 'DIAMETER_PXL', 'AREA_MM2', 'DIAMETER_MM'])
                columns=['INDEX_COL', 'AREA_PXL', 'DIAMETER_PXL', 'AREA_MM2', 'DIAMETER_MM','ENCL_CENTER'])
                  #columns=['INDEX_COL', 'AREA_PXL', 'DIAMETER_PXL', 'PERIMETER_PXL', 'ENCL_CENTER', 'ENCL_DIAMETER_PXL', 'AREA_MM2', 'DIAMETER_MM'])
    return path

def calc_AREA_PXL_diff(contour_df):
    encl_AREA_PXL = 3.1415 * (contour_df['ENCL_DIAMETER_PXL'] ** 2) / 4
    return abs(1 - contour_df['AREA_PXL'] / encl_AREA_PXL)

def prepare_df(contours):
    pd.options.display.float_format = '{:,.2f}'.format
    df = pd.DataFrame(contours, columns=['CONTOURS'])
    # df = pd.DataFrame(imutils.grab_contours(contours), columns=['CONTOURS'])
    df['HULL'] = df.apply(lambda x: cv2.convexHull(x['CONTOURS']), axis=1)
    df['AREA_PXL'] = df.apply(lambda x: cv2.contourArea(x['HULL']), axis=1)
    encl_circle = df.apply(lambda x: cv2.minEnclosingCircle(x['HULL']), axis=1)
    df['ENCL_CENTER'] = encl_circle.str[0]
    df['ENCL_DIAMETER_PXL'] = encl_circle.str[1] * 2
    df['PERIMETER_PXL'] = df.apply(lambda x: f"{cv2.arcLength(x['HULL'], True):.2f}", axis=1)

    return df

def filter_contours(contours, image_size):
    df = prepare_df(contours)
    # df = prepare_df(imutils.grab_contours(contours))

    # filter_other = df.apply(lambda x: x['AREA_PXL'] < 100 or x['AREA_PXL'] > 100000, axis=1)
    #filter_other = df.apply(lambda x: x['AREA_PXL'] < 100 or (x['AREA_PXL'] > 100 and calc_AREA_PXL_diff(x)) > 0.21, axis=1)
    #other_df = df[filter_other]
    #wo_other_df = df[~filter_other]

    plaque_minimum_size = image_size[1]/15
    if small_plaques:
        filter_green = df.apply(lambda x: x['AREA_PXL'] < 100000 and x['AREA_PXL'] > plaque_minimum_size and calc_AREA_PXL_diff(x) < 0.21 and x['ENCL_DIAMETER_PXL'] > 12, axis=1)
    else:
        filter_green = df.apply(lambda x: x['AREA_PXL'] < 100000 and x['AREA_PXL'] > plaque_minimum_size and calc_AREA_PXL_diff(x) < 0.21 and x['ENCL_DIAMETER_PXL'] > 30, axis=1)

    #filter_green = df.apply(lambda x: x['AREA_PXL'] < 100000 and x['AREA_PXL'] > 100 and calc_AREA_PXL_diff(x) < 0.21, axis=1)
    filter_plate = df.apply(lambda x: x['AREA_PXL'] > 100000 and calc_AREA_PXL_diff(x) < 0.21, axis=1)
    filter_red = df.apply(lambda x: x['AREA_PXL'] > plaque_minimum_size and calc_AREA_PXL_diff(x) > 0.21 and x['ENCL_DIAMETER_PXL'] < 30, axis=1)
    filter_other = df.apply(lambda x: x['AREA_PXL'] < plaque_minimum_size or calc_AREA_PXL_diff(x) > 0.21, axis=1)
    #filter_other = df.apply(lambda x: calc_AREA_PXL_diff(x) > 0.23, axis=1)
    #filter_other = df.copy()
    # filter_plate = filter_plate.apply(lambda x: calc_AREA_PXL_diff(x) < 0.21, axis=1)

    green_df = df[filter_green]
    red_df = df[filter_red]
    plate_df = df[filter_plate]
    other_df = df[filter_other]

    green_df.reset_index()
    green_df = reindex(green_df)

    red_df.reset_index()
    red_df = reindex(red_df)

    other_df.reset_index()
    other_df = reindex(other_df)

    plate_df.reset_index()
    plate_df = reindex(plate_df)

    return green_df, red_df, other_df, plate_df

def reindex(df):
    dfNew = df.copy()
    dfNew['INDEX_COL'] = df.index
    return dfNew

def unsharp_mask(image, kernel_size=(3, 3), sigma=1.0, amount=1.0, threshold=0):
    """Return a sharpened version of the image, using an unsharp mask."""
    # ET
    blurred = cv2.GaussianBlur(image, kernel_size, sigma)
    sharpened = float(amount + 1) * image - float(amount) * blurred
    sharpened = float(amount + 1) * image - float(amount) * image
    sharpened = np.maximum(sharpened, np.zeros(sharpened.shape))
    sharpened = np.minimum(sharpened, 255 * np.ones(sharpened.shape))
    sharpened = sharpened.round().astype(np.uint8)
    if threshold > 0:
        low_contrast_mask = np.absolute(image - blurred) < threshold
        np.copyto(sharpened, image, where=low_contrast_mask)
    return sharpened

def check_duplicate_centers(obj, obj_df):
    for x in obj_df.iterrows():

        if obj['INDEX_COL'] != x[0]:
            close_centers_x = abs(obj['ENCL_CENTER'][0] - x[1]['ENCL_CENTER'][0]) <= 5
            close_centers_y = abs(obj['ENCL_CENTER'][1] - x[1]['ENCL_CENTER'][1]) <= 5

            #print("obj encl center x = " + str(obj['ENCL_CENTER'][0]) + ", x encl center x = " + str( x[1]['ENCL_CENTER'][0]))
            #print("obj encl center y = " + str(obj['ENCL_CENTER'][1]) + ", x encl center y = " + str(x[1]['ENCL_CENTER'][1]))

            #print(str(x[0]) + ", difference 0 " + str(abs(obj['ENCL_CENTER'][0] - x[1]['ENCL_CENTER'][0])) + ", difference 1 " + str(abs(obj['ENCL_CENTER'][1] - x[1]['ENCL_CENTER'][1])))

            if close_centers_x is True and close_centers_y is True:
                return True

    return False

def check_duplicate_diameters(obj, obj_df):
    for x in obj_df.iterrows():
        if obj['INDEX_COL'] != x[0]:
            close_diameter = abs(obj['ENCL_DIAMETER_PXL'] - x[1]['ENCL_DIAMETER_PXL']) <= 25

            if close_diameter is True:
                return True

def getPlateSize(image):
    gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
    high_contrast = cv2.convertScaleAbs(gray, alpha=None, beta=0)
    # high_contrast = adjust_gamma(high_contrast, 1.0)

    gamma_test = adjust_gamma(gray, 7.1)
    # cv2.imwrite("./test_pic_green_gamma_0.jpg", gamma_test)

    binary_image = cv2.adaptiveThreshold(gamma_test, 500, cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY, 55, 2)

    contours = get_contours(binary_image)
    image_size = binary_image.shape
    plate_df = filter_contours(contours, image_size)
    output = draw_contours(gamma_test, plate_df, None, None)
    # write_images(out_dir_path, output, binary_image, high_contrast, "./out")

def check_duplicate_plaques(obj_df):
    for green in obj_df.iterrows():
        if check_duplicate_centers(green[1], obj_df):
            obj_df = obj_df.drop(obj_df[obj_df.index == green[0]].index)

    return obj_df

def check_duplicate_plates(obj_df):
    for green in obj_df.iterrows():
        if check_duplicate_diameters(green[1], obj_df):
            obj_df = obj_df.drop(obj_df[obj_df.index == green[0]].index)

    return obj_df

def calculate_size_mm(plate_size, obj_df, plate_df):

    if obj_df.size > 0:
        obj_df['DIAMETER_PXL'] = obj_df.apply(lambda x: f"{np.sqrt(float(x['AREA_PXL'])/np.pi)*2:.2f}",
                                                                axis=1)
    if plate_size and obj_df.size > 0:
        max_plate_diameter = plate_df['ENCL_DIAMETER_PXL'].max()
        pxl_per_mm = float(plate_size) / float(max_plate_diameter)
        obj_df['ENCL_DIAMETER_MM'] = obj_df.apply(lambda x: f"{x['ENCL_DIAMETER_PXL'] * pxl_per_mm:.2f}",
                                                                axis=1)
        obj_df['AREA_MM2'] = obj_df.apply(lambda x: f"{x['AREA_PXL'] * pxl_per_mm * pxl_per_mm:.2f}",
                                                        axis=1)
        obj_df['DIAMETER_MM'] = obj_df.apply(lambda x: f"{np.sqrt(float(x['AREA_MM2'])/np.pi)*2:.2f}",
                                                                axis=1)
    else:
        obj_df['ENCL_DIAMETER_MM'] = 0
        obj_df['AREA_MM2'] = 0
        obj_df['DIAMETER_MM'] = 0

    return obj_df

def get_brightness( im_file ):
   im = Image.open(im_file)
   im.convert('L')
   stat = ImageStat.Stat(im)
   return stat.mean[0]

def main(args):
    np.warnings.filterwarnings('ignore', category=np.VisibleDeprecationWarning)

    image_paths = get_image_paths(args['image'], args['directory'])
    plate_size = args['plate_size']

    global small_plaques
    global debug_mode
    small_plaques = args['small_plaque']
    debug_mode = args['debug']

    resList = []

    for image_path in image_paths:
        print("Processing " + image_path)

        # adjust too bright image
        #average image brightness
        image_brightness = get_brightness(image_path)
        im = Image.open(image_path)

        image = cv2.imread(image_path)
        if image_brightness > 70:
            #im = Image.open(image_tmp)
            enhancer = ImageEnhance.Brightness(im)
            factor = 0.5
            im_output = enhancer.enhance(factor)

            im_output.save(image_path + '.tmp' + os.path.splitext(image_path)[1])
            image = cv2.imread(image_path + '.tmp' + os.path.splitext(image_path)[1])

        binary_image, high_contrast, clr_high_contrast = process_image(image, 2.5)
        #       cv2.imshow("Binary image", binary_image)
        contours = get_contours(binary_image)
        image_size = binary_image.shape

        green_df, red_df, other_df, plate_df = filter_contours(contours, image_size)

        # Filter plaques duplicates (circle in circle) in valid plaques and plates
        green_df_copy = green_df.copy()
        plate_df_copy = plate_df.copy()
        green_df_copy = check_duplicate_plaques(green_df_copy)
        plate_df_copy = check_duplicate_plaques(plate_df_copy)

        red_df_copy = red_df.copy()
        other_df_copy = other_df.copy()

        green_df_copy = calculate_size_mm(plate_size, green_df_copy, plate_df)
        red_df_copy = calculate_size_mm(plate_size, red_df_copy, plate_df)
        other_df_copy = calculate_size_mm(plate_size, other_df_copy, plate_df)
        plate_df_copy = calculate_size_mm(plate_size, plate_df_copy, plate_df)

        if(green_df_copy.size > 0):
            # get Petri dish size and adjust plaques sizes
            green_df_copy['MEAN_COLOUR'] = green_df_copy.apply(lambda x: get_mean_grey_colour(high_contrast, x['CONTOURS']),
                                                               axis=1)

            green_df_copy['MEAN_COLOUR'] = green_df_copy.apply(
                lambda x: get_mean_grey_colour(clr_high_contrast, x['CONTOURS']),
                axis=1)

            #Remove extra black contours
            #TODO dark plaques are being filtered out (online_image.png
            filter_dev_colour = green_df_copy.apply(lambda x: abs(x['MEAN_COLOUR'] ) < 40, axis=1)
            green_df_copy = green_df_copy[~filter_dev_colour]
        else:
            green_df_copy['MEAN_COLOUR'] = 0

        green_df_copy = renumerate_df(green_df_copy)
        renumerate_df(red_df)
        renumerate_df(other_df)

        output = draw_contours(clr_high_contrast, green_df_copy, red_df_copy, other_df_copy, plate_df_copy)

        if not os.path.exists(out_dir_path):
            os.makedirs(out_dir_path)
        syncPath = write_images(out_dir_path, output, binary_image, high_contrast, image_path)

        #Delete the temporary file
        if os.path.exists(image_path + '.tmp' + os.path.splitext(image_path)[1]):
            os.remove(image_path + '.tmp' + os.path.splitext(image_path)[1])

        # format float values
        # green_df_copy['ENCL_DIAMETER_MM'] = green_df_copy['ENCL_DIAMETER_MM'].apply(lambda x: f"{x:.2f}")
        dataPath = write_data(out_dir_path, image_path, green_df_copy, red_df_copy, other_df, syncPath=syncPath)
        print("Process completed successfully")
        print(str(len(green_df_copy)) + ' plaques were found\n')

        resList.append((len(green_df_copy),syncPath,dataPath))

    return resList

def renumerate_df(df):
    df_new = df.copy()
    df = df_new.reset_index()

    df_new['INDEX_COL'] = df.index + 1
    return df_new

def get_mean_grey_colour(img, contour):
    contour_mask = np.zeros(img.shape[:2], dtype="uint8")
    contour_mask = cv2.drawContours(contour_mask, [contour], -1, 255, -1)
    mean = cv2.mean(img, contour_mask)
    return mean[0]

if __name__ == '__main__':
    args = util.parse_args()
    main(args)


text:
import os
import sys
import pytest
import inspect
import numpy as np
import atomsci.ddm.pipeline.featurization as feat
import atomsci.ddm.pipeline.parameter_parser as parse
import deepchem as dc
import atomsci.ddm.pipeline.model_datasets as model_dataset
import atomsci.ddm.pipeline.model_wrapper as model_wrapper
import atomsci.ddm.pipeline.model_pipeline as MP
import atomsci.ddm.pipeline.perf_data as perf_data

import utils_testing as utils
import copy

"""This testing script assumes that /ds/data/public/delaney/delaney-processed.csv is still on the same path on twintron. Assumes that the dataset_key: /ds/projdata/gsk_data/GSK_derived/PK_parameters/gsk_blood_plasma_partition_rat_crit_res_data.csv under the bucket gskdata and with the object_oid: 5af0e6368003ff018de33db5 still exists.
"""

#The dataset object from file is a delaney dataset using an ecfp featurizer with a default scaffold split.

datastore_is_down = utils.datastore_status()
MP_delaney_ecfp_train_valid_test_random = utils.delaney_pipeline()
(delaney_params, mdl_dataset_delaney, delaney_df) = utils.delaney_objects()

general_params = {'dataset_key' : './delaney-processed.csv',
'featurizer': 'ecfp',
'response_cols': 'measured log solubility in mols per litre',
'id_col': 'Compound ID',
'smiles_col': 'smiles',
'output_dir': 'pytest',
'model_type' : 'NN',
'splitter' : 'scaffold',
'prediction_type' : 'regression',
'baseline_epoch' : '7',
'max_epochs' : '10',
'datastore': 'False',
'save_results': 'False'}

DD = dc.data.datasets.DiskDataset

#***********************************************************************************
def test_create_model_wrapper():
    """
        Args:
        params (Namespace) : Parameters passed to the model pipeline
                     featurizer (Featurization): Object managing the featurization of compounds
                                 ds_client (DatastoreClient): Interface to the file datastore

                                                  Returns:
                                                  model (pipeline.Model): Wrapper for DeepChem, sklearn or other model.

                                                              Raises:
ValueError: Only params.model_type = 'NN' or 'RF' is supported.

Dependencies:
None

Calls:
DCNNModelWrapper, DCRFModelWrapper
    """
    inp_params = parse.wrapper(general_params)
    featurization = feat.create_featurization(inp_params)
    mdl = model_wrapper.create_model_wrapper(inp_params, featurization)
    mdl.setup_model_dirs()
    # testing for correct attribute initialization with model_type == "NN"
    test = []
    test.append(mdl.params.model_type == 'NN')
    test.append(isinstance(mdl.featurization,feat.DynamicFeaturization))
    test.append(mdl.output_dir == inp_params.output_dir)
    test.append(mdl.model_dir == inp_params.output_dir + '/' + 'model')
    test.append(mdl.best_model_dir == inp_params.output_dir + '/' + 'best_model')
    test.append(mdl.baseline_model_dir == inp_params.output_dir + '/' + 'baseline_epoch_model')
    test.append(mdl.transformers == [])
    test.append(mdl.transformers_x == [])
    test.append(isinstance(mdl, model_wrapper.DCNNModelWrapper))

    # testing for correct attribute initialization with model_type == "RF"
    temp_params = copy.deepcopy(inp_params)
    temp_params.model_type = 'RF'
    featurization = feat.create_featurization(temp_params)
    mdl_RF = model_wrapper.create_model_wrapper(temp_params, featurization)
    test.append(isinstance(mdl_RF, MP.model_wrapper.DCRFModelWrapper))
    test.append(mdl_RF.params.model_type == 'RF')

    # assertion for all tests
    assert all(test)

    #testing for Exception with model_type not in ['NN','RF']
    with pytest.raises(ValueError):
        temp_params.model_type = 'wrong'
        mdl_wrong = model_wrapper.create_model_wrapper(temp_params, featurization)

#***********************************************************************************
def test_super_create_transformers():
    """
    Args:
    model_dataset: The ModelDataset object that handles the current dataset

    Returns:
    self.transformers
    self.transformers_x
    self.params.transformer_key
    self.params.transformer_oid (if datastore)

    Raises:
    Exception when failing to save to the datastore

    Dependencies:
    create_featurization
    create_model_dataset
    model_dataset.load_full_dataset
    model_dataset.get_dataset_tasks
    model_dataset.check_task_columns
    model_dataset.get_featurized_data
    Requires (self.params.prediction_type == 'regression' and self.params.transformers == True) or len(self.transformers) > 0

    Calls:
    self.featurization.create_feature_transformer
    dsf.upload_pickle_to_DS

    """
    #set up for a model wrapper with regression and NN.

    inp_params = parse.wrapper(general_params)
    featurization = feat.create_featurization(inp_params)
    data_obj_ecfp = model_dataset.create_model_dataset(inp_params, featurization, ds_client = None)
    df_delaney = data_obj_ecfp.load_full_dataset()
    data_obj_ecfp.get_dataset_tasks(df_delaney)
    data_obj_ecfp.check_task_columns(df_delaney)
    data_obj_ecfp.get_featurized_data()
    mdl = model_wrapper.create_model_wrapper(inp_params, data_obj_ecfp.featurization)
    mdl.setup_model_dirs()

    #testing correct model_wrapper build with regression and NN
    test = []
    test.append(mdl.params.prediction_type == 'regression')
    test.append(mdl.params.model_type == 'NN')
    mdl.create_transformers(data_obj_ecfp)
    test.append(isinstance(mdl.transformers[0], dc.trans.transformers.NormalizationTransformer))
    test.append(mdl.transformers_x == [])
    #testing saving of transformer to correct location:
    transformer_path = os.path.join(mdl.output_dir, 'transformers.pkl')
    test.append(os.path.isfile(transformer_path))

    # TODO: test proper saving of the transformer to the datastore

    # TODO: test when transformers is False:
    inp_params.prediction_type = 'classification'
    mdl = model_wrapper.create_model_wrapper(inp_params, featurization)
    test.append(mdl.transformers == [])
    test.append(mdl.transformers_x == [])
    assert all(test)

#***********************************************************************************
def test_super_transform_dataset():
    """
    Args:
    dataset: The DeepChem DiskDataset that contains a dataset

    Returns:
    transformed_dataset

    Raises:
    None

    Dependencies:
    model_dataset.create_transformers

    Calls:
    None

    """
    #set up for a model wrapper with regression and NN.
    inp_params = parse.wrapper(general_params)
    featurization = feat.create_featurization(inp_params)
    data_obj_ecfp = model_dataset.create_model_dataset(inp_params, featurization, ds_client = None)
    df_delaney = data_obj_ecfp.load_full_dataset()
    data_obj_ecfp.get_dataset_tasks(df_delaney)
    data_obj_ecfp.check_task_columns(df_delaney)
    data_obj_ecfp.get_featurized_data()
    mdl = model_wrapper.create_model_wrapper(inp_params, data_obj_ecfp.featurization)
    mdl.setup_model_dirs()
    mdl.create_transformers(data_obj_ecfp)
    dataset = mdl.transform_dataset(data_obj_ecfp.dataset)

    test = []
    # checking that the dataset is the correct type
    test.append(isinstance(dataset, DD))
    # since this is not descriptor featurization, the X values for the datasets should be the same
    test.append((dataset.X == data_obj_ecfp.dataset.X).all())
    # and the response values should be the same length:
    test.append(len(dataset.y) == len(data_obj_ecfp.dataset.y))
    test.append(len(dataset.y) == len(dataset.ids))
    assert all(test)

#***********************************************************************************
#def test_train_NN_graphconv_scaffold_defaults():
#    """
#
#    Args:
#    pipeline (ModelPipeline): The ModelPipeline instance for this model run.
#
#    Dependencies:
#    ModelPipeline creation
#    featurization creation
#    creation of model_wrapper
#    mp.load_featurize_data
#
#    Calls:
#    create_perf_data
#    perf_data.accumulate_preds
#    perf_data.comput_perf_metrics
#    data.combined_training-data()
#    self._copy_model
#    """
#    # set up for the default graphconv model.
#    general_params['featurizer'] = 'graphconv'
#    inp_params = parse.wrapper(general_params)
#    mp = MP.ModelPipeline(inp_params)
#    mp.featurization = feat.create_featurization(inp_params)
#    mp.model_wrapper = model_wrapper.create_model_wrapper(inp_params, mp.featurization, mp.ds_client)
#    # asserting that the default model is created with the correct layer sizes, dropouts, model_dir, and mode
#    test2 = []
#    test2.append(mp.model_wrapper.params.layer_sizes == [64, 64, 128])
#    print(mp.model_wrapper.params.dropouts)
#    test2.append(mp.model_wrapper.params.dropouts == [0.25,0.25,0.25])
#    test2.append(isinstance(mp.model_wrapper.model, dc.models.tensorgraph.models.graph_models.GraphConvModel))
#    test2.append(mp.model_wrapper.model.model_dir == mp.model_wrapper.model_dir)
#    test2.append(mp.model_wrapper.model.graph_conv_layers == [64,64])
#    test2.append(mp.model_wrapper.model.dropout == [0.25,0.25,0.25])
#    test2.append(mp.model_wrapper.model.mode == 'regression')
#    test2.append(mp.model_wrapper.model.dense_layer_size == 128)
#    assert all(test2)
#
#    # calling dependencies for model_wrapper.train()
#    mp.model_wrapper.setup_model_dirs()
#    mp.load_featurize_data()
#    mp.model_wrapper.train(mp)
#
#    mdl = mp.model_wrapper
#    test3 = []
#    #asserting that attributes are properly created during train
#    test3.append(isinstance(mdl.data, model_dataset.FileDataset))
#    test3.append(isinstance(mdl.best_epoch, int))
#    test3.append(isinstance(mdl.train_perf_data, list))
#    test3.append(isinstance(mdl.train_perf_data[-1],perf_data.SimpleRegressionPerfData))
#    test3.append(isinstance(mdl.valid_perf_data, list))
#    test3.append(isinstance(mdl.valid_perf_data[-1],perf_data.SimpleRegressionPerfData))
#    test3.append(isinstance(mdl.train_epoch_perfs, np.ndarray))
#    test3.append(isinstance(mdl.valid_epoch_perfs, np.ndarray))
#    assert all(test3)
#
#    # For this particular dataset with graphconv and scaffold splitting, the training r2 tends to reach 0.4 by 10 epochs, setting the threshold to 0.3. If the training r2 is less than 0.4, there is most likely something wrong with the model training.
#    assert max(mdl.train_epoch_perfs) > 0.3
#    # For this particular dataset with graphconv, the valid r2 tends to reach 0.1 by 10 epochs
#    assert max(mdl.valid_epoch_perfs) > 0.1
#
#    # checking that the real_vals within train_perf_data and valid_perf_data match the datasets they originated from.
#    test4 = []
#    test4.append((mp.data.train_valid_dsets[0][0].y == mp.model_wrapper.train_perf_data[-1].real_vals).all())
#    test4.append((mp.data.train_valid_dsets[0][1].y == mp.model_wrapper.valid_perf_data[-1].real_vals).all())
#    # checking that the ids match as well
#    test4.append((mp.data.train_valid_dsets[0][0].ids == mp.model_wrapper.train_perf_data[-1].ids).all())
#    test4.append((mp.data.train_valid_dsets[0][1].ids == mp.model_wrapper.valid_perf_data[-1].ids).all())
#    # checking that the predicted values have the correct length and type
#    test4.append(len(mp.model_wrapper.train_perf_data[-1].pred_vals) == len(mp.data.train_valid_dsets[0][0].y))
#    test4.append(isinstance(mp.model_wrapper.train_perf_data[-1].pred_vals[0][0], np.float32))
#    test4.append(len(mp.model_wrapper.valid_perf_data[-1].pred_vals) == len(mp.data.train_valid_dsets[0][1].y))
#    test4.append(isinstance(mp.model_wrapper.valid_perf_data[-1].pred_vals[0][0], np.float32))
#    assert all(test4)

#***********************************************************************************
def test_train_NN_graphconv_scaffold_inputs():
    """

    Args:
    pipeline (ModelPipeline): The ModelPipeline instance for this model run.

    Dependencies:
    ModelPipeline creation
    featurization creation
    creation of model_wrapper
    mp.load_featurize_data

    Calls:
    create_perf_data
    perf_data.accumulate_preds
    perf_data.comput_perf_metrics
    data.combined_training-data()
    self._copy_model
    """
    # checking that the layers, dropouts, and learning rate are properly added to the deepchem graphconv model
    general_params['featurizer'] = 'graphconv'
    general_params['layer_sizes'] = '100,100,10'
    general_params['dropouts'] = '0.3,0.3,0.1'
    general_params['uncertainty'] = False
    inp_params = parse.wrapper(general_params)
    mp = MP.ModelPipeline(inp_params)
    mp.featurization = feat.create_featurization(inp_params)
    mp.model_wrapper = model_wrapper.create_model_wrapper(inp_params, mp.featurization, mp.ds_client)
    # asserting that the correct model is created with the correct layer sizes, dropouts, model_dir, and mode by default
    test1 = []
    test1.append(mp.model_wrapper.params.layer_sizes == [100, 100, 10])
    test1.append(mp.model_wrapper.params.dropouts == [0.3,0.3,0.1])
    # checking that parameters are properly passed to the deepchem model object
    test1.append(isinstance(mp.model_wrapper.model, dc.models.tensorgraph.models.graph_models.GraphConvModel))
    test1.append(mp.model_wrapper.model.model_dir == mp.model_wrapper.model_dir)
    test1.append(mp.model_wrapper.model.graph_conv_layers == [100,100])
    test1.append(mp.model_wrapper.model.dropout == [0.3,0.3,0.1])
    test1.append(mp.model_wrapper.model.mode == 'regression')
    test1.append(mp.model_wrapper.model.dense_layer_size == 10)
    assert all(test1)

    #***********************************************************************************
    def test_super_get_train_valid_pred_results():
        """
        Args:
        perf_data: A PerfData object that stores the predicted values and metrics
        Returns:
        dict: A dictionary of the prediction results

            Raises:
        None

        Dependencies:
        create_perf_data

        Calls:
        perf_data.get_prediction_results()

        """
        pass
    # should be tested in perf_data.get_prediction_results()
    # should still be called to make sure that the function is callable

    #***********************************************************************************
    def test_super_get_test_perf_data():
        """
        Args:
        model_dir (str): Directory where the saved model is stored
        model_dataset (DiskDataset): Stores the current dataset and related methods

        Returns:
        perf_data: PerfData object containing the predicted values and metrics for the current test dataset

            Raises:
        None

        Dependencies:
        A model must be in model_dir
        model_dataset.test_dset must exist

        Calls:
        create_perf_data
        self.generate_predictions
        perf_data.accumulate_preds
        """
        pass
        # mostly tested in accumulate_preds, but should be tested to ensure taht the predictions are properly being called

    #***********************************************************************************
    def test_super_get_test_pred_results():
        """
        Args:
        model_dir (str): Directory where the saved model is stored
        model_dataset (DiskDataset): Stores the current dataset and related methods

        Returns:
        dict: A dictionary containing the prediction values and metrics for the current dataset.

            Raises:
        None

        Dependencies:
        A model must be in model_dir
        model_dataset.test_dset must exist

        Calls:
        self.get_test_perf_data
        perf_data.get_prediction_results
        """
        pass
        #mostly tested in perf_data.get_prediction_results

    #***********************************************************************************
    def test_super_get_full_dataset_perf_data():
        """
        Args:
        model_dataset (DiskDataset): Stores the current dataset and related methods

        Returns:
        perf_data: PerfData object containing the predicted values and metrics for the current full dataset

            Raises:
        None

        Dependencies:
        A model must already be trained

        Calls:
        create_perf_data
        self.generate_predictions
        self.accumulate_preds
        """
        pass

    #***********************************************************************************
    def test_super_get_full_dataset_pred_results():
        """
        Args:
        model_dataset (DiskDataset): Stores the current dataset and related methods
        Returns:
        dict: A dictionary containing predicted values and metrics for the current full dataset

            Raises:
        None

        Dependencies:
        A model was already be trained.

        Calls:
        get_full_dataset_perf_data
        self.get_prediction_results()
        """
        pass

ext: py

max_stars_repo_name: vgutta/AMPL

max_issues_repo_name: vgutta/AMPL

max_forks_repo_name: vgutta/AMPL

reason: import numpy


text:
import requests
import sys
import h5py
import numpy as np

def get(path, params=None, savedir=None):
    # make HTTP GET request to path
    headers = {"api-key":"44140abd45569b745a4e3e7cf1284ce0"}
    r = requests.get(path, params=params, headers=headers)

    # raise exception if response code is not HTTP SUCCESS (200)
    r.raise_for_status()

    if r.headers['content-type'] == 'application/json':
        return r.json() # parse json responses automatically

    if 'content-disposition' in r.headers:
        filename = r.headers['content-disposition'].split("filename=")[1]
        if savedir != None:
            filename = savedir + filename
        with open(filename, 'wb') as f:
            f.write(r.content)
        return filename # return the filename string

    return r

def HaloProgenitors(haloID):
    '''
    haloID is the subhalo's ID in snap_099
    return a dict = {'SnapNum' : SubfindID}
    '''
    url = 'http://www.tng-project.org/api/TNG100-1/snapshots/99/subhalos/%d'%haloID
    try:
        info = get(url)
        sublink = get(info['trees']['sublink_mpb'])
    except:
        print(sys.exc_info()[0])
        return -1

    with h5py.File(sublink, 'r') as f:
        #Find halo's Subfind ID with redshift(ie:SnapNum), and save the dict in '/Raid0/zhouzb/diskHalo_Sublink/'
        snap_num = np.array(f['SnapNum'])
        subfind_ID = np.array(f['SubfindID'])

    Progenitors_dict = {}
    for i in range(len(snap_num)):
        Progenitors_dict['%d'%snap_num[i]] = subfind_ID[i]

    return Progenitors_dict

def SelectDisk(snap_num):
    '''
    Input Snapshot number, like snap_num = 99 (z=0)
    Select disk galaxies, return haloID of them.
    '''
    #select halo stellar particles > 40000
    with h5py.File('/Raid0/zhouzb/TNGdata/offsets_0%d.hdf5'%snap_num,'r') as offset:
        haloSBT = (np.array(offset['Subhalo']['SnapByType']))[:,4]
    #Total halo number, it also be an index of haloID

    halolen = haloSBT[1:]
    halolen = np.append(halolen, halolen.max()) - haloSBT

    with h5py.File('/Raid0/zhouzb/TNGdata/stellar_circs.hdf5','r') as cir:
        haloID = np.array(cir['Snapshot_%d'%snap_num]['SubfindID'])
        cir07frac = np.array(cir['Snapshot_%d'%snap_num]['CircAbove07Frac'])
        MassTensor = np.array(cir['Snapshot_%d'%snap_num]['MassTensorEigenVals'])
    #circularity parameter ϵ > 0.2
    cir_mask = cir07frac > 0.2

    #flatness of the galaxy is defined as the ratio M1=(M2M3)**0.5 , disk galaxy's flatness < 0.7
    MassTensor = MassTensor[cir_mask]
    haloID = haloID[cir_mask]

    flat = MassTensor[:,0]/(MassTensor[:,1]*MassTensor[:,2])**0.5
    flat_mask = flat < 0.7

    haloID = haloID[flat_mask]

    mas_mask = halolen[haloID] >= 40000
    haloID = haloID[mas_mask]
    return haloID

'''
snap_91 z=0.1
snap_84 z=0.2
snap_78 z=0.3
snap_72 z=0.4
snap_67 z=0.5
snap_59 z=0.7
snap_50 z=1.0
snap_40 z=1.5
snap_33 z=2.0
'''

barred = np.load('F:/Linux/data/099fig/barredID.npy')
snap = [67, 59, 50, 40, 33]
errorHalo = []

for z in snap:
    print('Downloading snap_%d'%z)
    disk = SelectDisk(z)

    for haloID in disk:
        cutoff_url = 'http://www.tng-project.org/api/TNG100-1/snapshots/%d/subhalos/%d/cutout.hdf5?stars=Masses,Coordinates,Velocities,GFM_StellarFormationTime'%(z, haloID)
        try:
            get(cutoff_url, savedir='F:\\Linux\\data\\TNG\\snap_%d\\'%z)
        except:
            print("halo %d in snap_%d Fail:"%(haloID, z), sys.exc_info()[0])

            errorHalo.append(haloID)
            break
        else:
            print('halo %d in snap_%d downloaded'%(haloID, z))

    print('halo %d in all snapshot download Completed'%haloID)

if len(errorHalo) == 0:
    print('All done.')
else:
    print('%d halo download faild'%len(errorHalo))
    print("Error halo's ID were saved in 'F:\\Linux\\data\\TNG\\downError.log.npy'.")
    np.save('F:\\Linux\\data\\TNG\\downError.log.npy', errorHalo)


text:
def PieChart1(filepath):
    import pandas as pd
    import numpy as np

    data1 = pd.read_csv(filepath)
    data1_clean = (
    data1
        .dropna(axis=0)
        .sort_values(by = 'yearID', ascending = True)
    )
    AwardPlayer = pd.read_csv(filepath)
    AwardPlayer_clean = (
    AwardPlayer
        .drop(['lgID','tie','notes'],axis=1)
        .dropna(axis = 0)
    )
    merged = data1_clean.merge(AwardPlayer_clean, on="playerID", how="outer")
    merged_clean = (
    merged
        .dropna(axis=0)
    )

    new_df = merged_clean.loc[merged_clean['awardID'].isin(['TSN All-Star','Gold Glove','Silver Slugger','Most Valuable Player ','TSN Pitcher of the Year','Cy Young Award','Rookie of the Year','TSN Fireman of the Year','TSN Player of the Year']
    )]

    return new_df

def PieChart2(filepath):
    import pandas as pd
    import numpy as np

    data1 = pd.read_csv(filepath)
    data1_clean = (
    data1
        .dropna(axis=0)
        .sort_values(by = 'yearID', ascending = True)
    )
    HallofFame = pd.read_csv(filepath)
    HallofFame_clean = (
    HallofFame
        .drop(['votedBy','needed_note','category'],axis=1)
        .dropna(axis = 0)
    )
    merged = data1_clean.merge(HallofFame_clean, on="playerID", how="outer")
    merged_clean = (
    merged
        .dropna(axis=0)
    )
    new_df = merged_clean.loc[merged_clean['teamID'].isin(['BRO','SLN','PIT','BSN','CHN','CIN','WS1','TEX','NYA'])]

    return new_df2

def PieChart3(filepath):
    import pandas as pd
    import numpy as np

    AwardPlayer = pd.read_csv(filepath)
    AwardPlayer_clean = (
    AwardPlayer
        .drop(['lgID','tie','notes'],axis=1)
        .dropna(axis = 0)
    )
    Salaries = pd.read_csv(filepath)
    Salaries_clean = (
    Salaries
        .drop(['lgID'],axis=1)
        .dropna(axis = 0)
    )
    merged = AwardPlayer_clean.merge(Salaries_clean, on="playerID", how="outer")
    merged_clean = (
    merged
        .dropna(axis=0)
    )
    new_df = merged_clean.loc[merged_clean['awardID'].isin(['TSN Player of the Year','Pitching Triple Crown','Triple Crown','Rookie of the Year','TSN All-Star','Hutch Award','Gold Glove','World Series MVP']
    )]

    return new_df3

def Plotgraph1(filepath):
    import pandas as pd
    import numpy as np

    data1 = pd.read_csv(filepath)
    old_games = data1[data1['yearID'] < 1891]
    new_games = data1[data1['yearID'] > 2000]
    new_games = new_games.dropna().reset_index(drop=True)
    new_games.head()

    return new_df4


text:
import numpy as np

array=np.array([2,4,56,8,98])

def MaxOfArray(array):
  return np.max(array)
  print(MaxOfArray)

a=np.array([2,4,5])
b=["asf","asfh","hkjhsdf"]
c="casa"
d={"name":5}
l=(3,5)
print()

#Differenza tra array e lista è che la lista è un contenitore di oggetti

listval=[a,b,c,d,l]

#L'array è numerico e si possono fare delle operazioni negli oggetti contenuti all'interno

#La lista ha una dimensione, l'array ha una shape, cioè può avere più dimensioni, come una matrice, mentre una lista è vettoriale

a2=np.array([[1,2,3],[4,5,6],[7,8,9]])

print(a2.shape)

#Come vediamo qua ci da la shape di a2 che ha 3 righe e 3 colonne

max_stars_repo_name: therealsapo/machine-learning-fav

max_issues_repo_name: therealsapo/machine-learning-fav

max_forks_repo_name: therealsapo/machine-learning-fav


text:
import numpy as np

DIFF_CASE = 32  # distance between lower and upper case in ASCII

def p(*args, verbose=True):
    if verbose:
        print(*args)

def read_input(filepath):
    with open(filepath) as f:
        data = f.read().rstrip()
    return data

def encode_polymer(string):
    chars = [ord(c) for c in string]
    return np.array(chars)

def decode_polymer(numbers):
    chars = [chr(n) for n in numbers]
    return ''.join(chars)

def react(polymer):
    while True:
        diffs = np.abs(np.diff(polymer))
        reactions = diffs == DIFF_CASE
        if not reactions.any():
            break
        polymer = reduce_polymer(polymer, reactions)
    return polymer

# Recursive
def react_recursive(polymer):
    """
    Needs something like sys.setrecursionlimit(len(data)) in main first
    """
    diffs = np.abs(np.diff(polymer))
    reactions = diffs == DIFF_CASE
    if reactions.any():
        polymer = reduce_polymer(polymer, reactions)
        return react_recursive(polymer)
    else:
        return polymer

def reduce_polymer(polymer, reactions):
    # Where two letters are going to disappear
    indices = np.flatnonzero(reactions)
    # But not if there are three together
    indices = np.delete(indices, np.flatnonzero(np.diff(indices) == 1) + 1)
    short_polymer = np.delete(polymer, [indices, indices + 1])
    return short_polymer

def remove_monomers(string, char):
    lower = char
    upper = char.upper()
    stripped = string.replace(lower, '').replace(upper, '')
    return stripped

def get_length_shortest_polymer(string):
    all_chars = set(string.lower())
    lengths = np.empty(len(all_chars), np.uint16)
    for i, char in enumerate(all_chars):
        substring = remove_monomers(string, char)
        polymer = encode_polymer(substring)
        resulting_polymer = decode_polymer(react(polymer))
        lengths[i] = len(resulting_polymer)
    return lengths.min()

def part_1(string):
    polymer = encode_polymer(string)
    short_polymer = react(polymer)
    result_string = decode_polymer(short_polymer)
    return result_string

def part_2(string):
    answer = get_length_shortest_polymer(string)
    return answer

if __name__ == "__main__":
    verbose = False

    example_data = read_input('example_1.txt')
    data = read_input('input.txt')

    p('Part 1', verbose=verbose)
    example_result_string = part_1(example_data)
    example_1 = len(example_result_string)
    p('Example:', example_1, verbose=verbose)
    result_string = part_1(data)
    answer_1 = len(result_string)
    p('Answer:', answer_1, verbose=verbose)

    p(verbose=verbose)

    p('Part 2', verbose=verbose)
    example_2 = part_2(example_result_string)
    p('Example:', example_2, verbose=verbose)
    answer_2 = part_2(result_string)
    p('Answer:', answer_2, verbose=verbose)


text:
# Author: Mathieu Blondel
#         Derek Lim
# License: BSD 3 clause

import argparse
import sys
from pathlib import Path
import os

import numpy as np
import matplotlib.pylab as plt
from matplotlib import patches as patches
from sklearn.metrics.pairwise import euclidean_distances
try:
    import joblib
except ImportError:
    from sklearn.externals import joblib
import ot

from smoothot.dual_solvers import solve_semi_dual, get_plan_from_semi_dual
from smoothot.dual_solvers import NegEntropy, SquaredL2
import dataset

# make needed directories if they do not already exist
root_dir = os.path.dirname(os.path.abspath(__file__))
Path(os.path.join(root_dir, 'images')).mkdir(exist_ok=True)
Path(os.path.join(root_dir, 'res')).mkdir(exist_ok=True)

parser = argparse.ArgumentParser()
parser.add_argument('--n_colors', type=int, default=256, help='number of color clusters')
parser.add_argument('--method', type=str, default='l2_sd', help='OT method')
parser.add_argument('--gamma', type=float, default=1.0, help='regularization parameter')
parser.add_argument('--max_iter', type=int, default=1000)
parser.add_argument('--img1', type=str, default='comunion')
parser.add_argument('--img2', type=str, default='autumn')
args = parser.parse_args()

def map_img(T, img1, img2, weights1):
    """Transfer colors from img2 to img1"""
    return np.dot(T / weights1[:, np.newaxis], img2)

n_colors = args.n_colors
method = args.method
gamma = args.gamma
max_iter = args.max_iter
img1 = args.img1
img2 = args.img2
pair = img1+'-'+img2

print("images:", img1, '|', img2)
print("n_colors:", n_colors)
print("gamma:", gamma)
print("max_iter:", max_iter)
print()

hist1, hist2, C, centers1, centers2, labels1, labels2, shape1, shape2 = \
dataset.load_color_transfer(img1=img1, img2=img2, n_colors=n_colors,
                            transpose=False)
m = len(hist1)
n = len(hist2)

# Obtain transportation plan.
if method == "l2_sd":
    regul = SquaredL2(gamma=gamma)
    alpha = solve_semi_dual(hist1, hist2, C, regul, max_iter=max_iter, tol=1e-6)
    T = get_plan_from_semi_dual(alpha, hist2, C, regul)
    name = "Squared 2-norm"
elif method == "ent_sd":
    regul = NegEntropy(gamma=gamma)
    alpha = solve_semi_dual(hist1, hist2, C, regul, max_iter=max_iter, tol=1e-6)
    T = get_plan_from_semi_dual(alpha, hist2, C, regul)
    name = "Entropy"
elif method == "lp":
    T = ot.emd(hist1, hist2, C)
    name = "Unregularized"
else:
    raise ValueError("Invalid method")

sparsity = np.sum(T > 1e-10) / T.size
print("Sparsity:", sparsity)

T1 = np.sum(T, axis=1)
Tt1 = np.sum(T, axis=0)
err_a = hist1 - np.sum(T, axis=1)
err_b = hist2 - np.sum(T, axis=0)
print("Marginal a", np.dot(err_a, err_a))
print("Marginal b", np.dot(err_b, err_b))

#print(np.sum(T * C) + 0.5 / gamma * np.dot(err_a, err_a) + 0.5 / gamma * np.dot(err_b, err_b))
print('Objective value:', np.sum(T * C))

T_ = ot.emd(hist1, hist2, C)
print('Unregularized objective value:', np.sum(T_ * C))

img1 = centers1[labels1]
img2 = centers2[labels2]

centers1_mapped = map_img(T, centers1, centers2, T1)
img1_mapped = centers1_mapped[labels1]

centers2_mapped = map_img(T.T, centers2, centers1, Tt1)
img2_mapped = centers2_mapped[labels2]

fig = plt.figure(figsize=(6,6))

ax = fig.add_subplot(221)
ax.imshow(img1.reshape(shape1))
ax.axis("off")

ax = fig.add_subplot(222)
ax.imshow(img2.reshape(shape2))
ax.axis("off")

ax = fig.add_subplot(223)
ax.imshow(img1_mapped.reshape(shape1))
ax.axis("off")

ax = fig.add_subplot(224)
ax.imshow(img2_mapped.reshape(shape2))
ax.axis("off")

plt.tight_layout()

# plot original and transformed images
out = "%s/images/%s_%d_%s_%0.3e.jpg" % (root_dir, method, n_colors, pair, gamma)
plt.savefig(out)
print()
print('Saved image to:', out)

out = "%s/res/img_%s_%d_%s_%0.3e.pkl" % (root_dir, method, n_colors, pair, gamma)
tup = (img1_mapped.reshape(shape1), img2_mapped.reshape(shape2), sparsity)

joblib.dump(tup, out)
print('Saved pickle to:', out)

plt.show()

if n_colors <= 32:
    # plot transport plan and color histogram

    def draw_blocks(ax, T):
        # find contiguous chunks between coefficients
        for k, attn_row in enumerate(T):
            brk = np.diff(attn_row)
            brk = np.where(brk != 0)[0]
            brk = np.append(0, brk + 1)
            brk = np.append(brk, T.shape[0])

            right_border = True
            for s, t in zip(brk[:-1], brk[1:]):
                if attn_row[s:t].sum() == 0:
                    right_border = False
                    continue
                lines = [(s, k), (t, k), (t, k + 1), (s, k + 1)]
                lines = np.array(lines, dtype=np.float) - 0.5
                path = patches.Polygon(lines, facecolor='none', linewidth=1.5,
                                       alpha=1, joinstyle='round',
                                       closed=not right_border,
                                       edgecolor='#999999')
                ax.add_patch(path)
                right_border = True

    def draw_border(ax, T):
        lines = [(0, 0), (0, T.shape[1]), (T.shape[0], T.shape[1]), (T.shape[0], 0)]
        lines = np.array(lines) - 0.5
        path = patches.Polygon(lines, facecolor='none', linewidth=1.5,
                                   alpha=1, joinstyle='round',
                                   closed=True, edgecolor='#999999')
        scatter_axes.add_patch(path)

    plt.figure(figsize=(5,5))

    scatter_axes = plt.subplot2grid((3, 3), (1, 0), rowspan=2, colspan=2)
    x_hist_axes = plt.subplot2grid((3, 3), (0, 0), colspan=2,
                                   sharex=scatter_axes)

    scatter_axes.imshow(T, cmap=plt.cm.Reds, interpolation="nearest")
    scatter_axes.axis("off")

    scatter_axes.set_xlabel(name)

    draw_blocks(scatter_axes, T)
    draw_border(scatter_axes, T)

    bar_list = x_hist_axes.bar(np.arange(m), height=hist1,
                               width=(0.80 * n_colors) / m, edgecolor="#999999",
                               linewidth=1.0, color=centers1)
    x_hist_axes.axis("off")

    #if method in ("l2_primal", "l2_sp"):
    if True:
        y_hist_axes = plt.subplot2grid((3, 3), (1, 2), rowspan=2,
                                   sharey=scatter_axes)
        bar_list = y_hist_axes.barh(np.arange(n), width=hist2,
                                    height=(0.80 * n_colors) / n,
                                    edgecolor="#999999", linewidth=1.0, color=centers2)
        y_hist_axes.axis("off")

    plt.tight_layout()
    plt.subplots_adjust(wspace=0, hspace=0, bottom=0.12)

    plt.show()

max_stars_repo_name: mblondel/smooth-ot

max_issues_repo_name: mblondel/smooth-ot

max_forks_repo_name: mblondel/smooth-ot


text:
import numpy as np
import re
from .basic import BasicProducer

def _extract_index(name):
    token = '\[\d+\]'
    match = re.search(token, name)
    if match:
        substr = match.group(0)
        index = int(substr.strip('[]'))
        newName = name.replace(substr, '')

        return newName, index
    return name, None

def _extract_values(obj, namesAndIndices):
    for (name, index) in namesAndIndices:
        if index is not None:
            yield getattr(obj, name)[index]
        else:
            yield getattr(obj, name)

class InvariantMass(BasicProducer):

    def __init__(self, outputName, fourMomentum1, fourMomentum2):
        '''
            fourMomentum = [energy, px, py, pz]
        '''
        super(InvariantMass, self).__init__(outputName)
        self.input1 = [_extract_index(x) for x in fourMomentum1]
        self.input2 = [_extract_index(x) for x in fourMomentum2]

        self.v1 = []
        self.v2 = []
        self.v = []

    def _value(self, obj):
        self.v1[:] = list(_extract_values(obj, self.input1))
        self.v2[:] = list(_extract_values(obj, self.input2))
        self.v[:] = np.add(self.v1, self.v2)

        energySquared = pow(self.v[0], 2)
        momentumSquared = np.dot(self.v[1:], self.v[1:])
        return [np.sqrt(energySquared - momentumSquared)]

max_stars_repo_name: benkrikler/alphatwirl-interface

max_issues_repo_name: benkrikler/alphatwirl-interface

max_forks_repo_name: benkrikler/alphatwirl-interface


text:
# encoding=utf8

"""Evolution strategy algorithm module."""

from math import ceil

import numpy as np
from numpy import random as rand
from numpy.linalg import (
	norm,
	cholesky as chol,
	eig,
	solve,
	lstsq
)

from WeOptPy.algorithms.interfaces import (
	Algorithm,
	Individual,
	default_individual_init
)
from WeOptPy.util.utility import objects2array

__all__ = [
	'EvolutionStrategy1p1',
	'EvolutionStrategyMp1',
	'EvolutionStrategyMpL',
	'EvolutionStrategyML',
	'CovarianceMatrixAdaptionEvolutionStrategy',
	'IndividualES'
]

class IndividualES(Individual):
	r"""Individual for Evolution Strategies.

	See Also:
		* :class:`WeOptPy.algorithms.Individual`
	"""
	def __init__(self, **kwargs):
		r"""Initialize individual.

		Args:
			kwargs (Dict[str, Any]): Additional arguments.

		See Also:
			* :func:`WeOptPy.algorithms.Individual.__init__`
		"""
		Individual.__init__(self, **kwargs)
		self.rho = kwargs.get('rho', 1)

class EvolutionStrategy1p1(Algorithm):
	r"""Implementation of (1 + 1) evolution strategy algorithm. Uses just one individual.

	Algorithm:
		(1 + 1) Evolution Strategy Algorithm

	Date:
		2018

	Authors:
		Klemen Berkovič

	License:
		MIT

	Reference URL:

	Reference paper:

	Attributes:
		Name (List[str]): List of strings representing algorithm names.
		mu (int): Number of parents.
		k (int): Number of iterations before checking and fixing rho.
		c_a (float): Search range amplification factor.
		c_r (float): Search range reduction factor.

	See Also:
		* :class:`WeOptPy.algorithms.Algorithm`
	"""
	Name = ['EvolutionStrategy1p1', 'EvolutionStrategy(1+1)', 'ES(1+1)']

	@staticmethod
	def type_parameters():
		r"""Get dictionary with functions for checking values of parameters.

		Returns:
			Dict[str, Callable]:
				* mu (Callable[[int], bool])
				* k (Callable[[int], bool])
				* c_a (Callable[[Union[float, int]], bool])
				* c_r (Callable[[Union[float, int]], bool])
				* epsilon (Callable[[float], bool])
		"""
		return {
			'mu': lambda x: isinstance(x, int) and x > 0,
			'k': lambda x: isinstance(x, int) and x > 0,
			'c_a': lambda x: isinstance(x, (float, int)) and x > 1,
			'c_r': lambda x: isinstance(x, (float, int)) and 0 < x < 1,
			'epsilon': lambda x: isinstance(x, float) and 0 < x < 1
		}

	def set_parameters(self, mu=1, k=10, c_a=1.1, c_r=0.5, epsilon=1e-20, **ukwargs):
		r"""Set the arguments of an algorithm.

		Args:
			mu (Optional[int]): Number of parents
			k (Optional[int]): Number of iterations before checking and fixing rho
			c_a (Optional[float]): Search range amplification factor
			c_r (Optional[float]): Search range reduction factor

		See Also:
			* :func:`WeOptPy.algorithms.Algorithm.setParameters`
		"""
		Algorithm.set_parameters(self, n=mu, itype=ukwargs.pop('itype', IndividualES), **ukwargs)
		self.mu, self.k, self.c_a, self.c_r, self.epsilon = mu, k, c_a, c_r, epsilon

	def mutate(self, x, rho):
		r"""Mutate individual.

		Args:
			x (Individual): Current individual.
			rho (float): Current standard deviation.

		Returns:
			Individual: Mutated individual.
		"""
		return x + self.normal(0, rho, len(x))

	def update_rho(self, rho, k):
		r"""Update standard deviation.

		Args:
			rho (float): Current standard deviation.
			k (int): Number of succesfull mutations.

		Returns:
			float: New standard deviation.
		"""
		phi = k / self.k
		if phi < 0.2: return self.c_r * rho if rho > self.epsilon else 1
		elif phi > 0.2: return self.c_a * rho if rho > self.epsilon else 1
		else: return rho

	def init_population(self, task):
		r"""Initialize starting individual.

		Args:
			task (Task): Optimization task.

		Returns:
			Tuple[Individual, float, list, dict]:
				1, Initialized individual.
				2, Initialized individual fitness/function value.
				3. Additional arguments.
				4. Additional keyword arguments:
					* ki (int): Number of successful rho update.
		"""
		c, ki = IndividualES(task=task, rnd=self.Rand), 0
		return c, c.f, [], {'ki': ki}

	def run_iteration(self, task, c, fpop, xb, fxb, ki, *args, **dparams):
		r"""Core function of EvolutionStrategy(1+1) algorithm.

		Args:
			task (Task): Optimization task.
			pop (Individual): Current position.
			fpop (float): Current position function/fitness value.
			xb (Individual): Global best position.
			fxb (float): Global best function/fitness value.
			ki (int): Number of successful updates before rho update.
			args (list): Additional arguments.
			dparams (dict): Additional keyword arguments.

		Returns:
			Tuple[Individual, float, Individual, float, list, dict]:
				1, Initialized individual.
				2, Initialized individual fitness/function value.
				3. New global best solution.
				4. New global best solutions fitness/objective value.
				5. Additional arguments.
				6. Additional keyword arguments:
					* ki (int): Number of successful rho update.
		"""
		if task.Iters % self.k == 0: c.rho, ki = self.update_rho(c.rho, ki), 0
		cn = objects2array([task.repair(self.mutate(c.x, c.rho), self.Rand) for _i in range(self.mu)])
		cn_f = np.asarray([task.eval(cn[i]) for i in range(len(cn))])
		ib = np.argmin(cn_f)
		if cn_f[ib] < c.f:
			c.x, c.f, ki = cn[ib], cn_f[ib], ki + 1
			if cn_f[ib] < fxb: xb, fxb = self.get_best(cn[ib], cn_f[ib], xb, fxb)
		return c, c.f, xb, fxb, args, {'ki': ki}

class EvolutionStrategyMp1(EvolutionStrategy1p1):
	r"""Implementation of (mu + 1) evolution strategy algorithm. Algorithm creates mu mutants but into new generation goes only one individual.

	Algorithm:
		(:math:`\mu + 1`) Evolution Strategy Algorithm

	Date:
		2018

	Authors:
		Klemen Berkovič

	License:
		MIT

	Reference URL:

	Reference paper:

	Attributes:
		Name (List[str]): List of strings representing algorithm names.

	See Also:
		* :class:`WeOptPy.algorithms.EvolutionStrategy1p1`
	"""
	Name = ['EvolutionStrategyMp1', 'EvolutionStrategy(mu+1)', 'ES(m+1)']

	def set_parameters(self, **kwargs):
		r"""Set core parameters of EvolutionStrategy(mu+1) algorithm.

		Args:
			kwargs (Dict[str, Any]):

		See Also:
			* :func:`WeOptPy.algorithms.EvolutionStrategy1p1.setParameters`
		"""
		mu = kwargs.pop('mu', 40)
		EvolutionStrategy1p1.set_parameters(self, mu=mu, **kwargs)

class EvolutionStrategyMpL(EvolutionStrategy1p1):
	r"""Implementation of (mu + lambda) evolution strategy algorithm. Mulation creates lambda individual. Lambda individual compete with mu individuals for survival, so only mu individual go to new generation.

	Algorithm:
		(:math:`\mu + \lambda`) Evolution Strategy Algorithm

	Date:
		2018

	Authors:
		Klemen Berkovič

	License:
		MIT

	Reference URL:

	Reference paper:

	Attributes:
		Name (List[str]): List of strings representing algorithm names
		lam (int): TODO

	See Also:
		* :class:`WeOptPy.algorithms.EvolutionStrategy1p1`
	"""
	Name = ['EvolutionStrategyMpL', 'EvolutionStrategy(mu+lambda)', 'ES(m+l)']

	@staticmethod
	def type_parameters():
		r"""TODO.

		Returns:
			Dict[str, Any]:
				* lam (Callable[[int], bool]): TODO.

		See Also:
			* :func:`WeOptPy.algorithms.EvolutionStrategy1p1`
		"""
		d = EvolutionStrategy1p1.type_parameters()
		d['lam'] = lambda x: isinstance(x, int) and x > 0
		return d

	def set_parameters(self, lam=45, **ukwargs):
		r"""Set the arguments of an algorithm.

		Arguments:
			lam (int): Number of new individual generated by mutation.

		See Also:
			* :func:`WeOptPy.algorithms.EvolutionStrategy1p1.setParameters`
		"""
		EvolutionStrategy1p1.set_parameters(self, init_pop_func=default_individual_init, **ukwargs)
		self.lam = lam

	def update_rho(self, pop, k):
		r"""Update standard deviation for population.

		Args:
			pop (numpy.ndarray): Current population.
			k (int): Number of successful mutations.
		"""
		phi = k / self.k
		if phi < 0.2:
			for i in pop: i.rho = self.c_r * i.rho
		elif phi > 0.2:
			for i in pop: i.rho = self.c_a * i.rho

	def change_count(self, c, cn):
		r"""Update number of successful mutations for population.

		Args:
			c (numpy.ndarray): Current population.
			cn (numpy.ndarray): New population.

		Returns:
			int: Number of successful mutations.
		"""
		k = 0
		for e in cn:
			if e not in c: k += 1
		return k

	def mutate_rand(self, pop, task):
		r"""Mutate random individual form population.

		Args:
			pop (numpy.ndarray): Current population.
			task (Task): Optimization task.

		Returns:
			numpy.ndarray: Random individual from population that was mutated.
		"""
		i = self.randint(self.mu)
		return task.repair(self.mutate(pop[i].x, pop[i].rho), rnd=self.Rand)

	def init_population(self, task):
		r"""Initialize starting population.

		Args:
			task (Task): Optimization task.

		Returns:
			Tuple[numpy.ndarray, numpy.ndarray, list, dict]:
				1. Initialized populaiton.
				2. Initialized populations function/fitness values.
				3. Additional arguments.
				4. Additional keyword arguments:
					* ki (int): Number of successful mutations.

		See Also:
			* :func:`WeOptPy.algorithms.Algorithm.initPopulation`
		"""
		c, fc, args, d = Algorithm.init_population(self, task)
		d.update({'ki': 0})
		return c, fc, args, d

	def run_iteration(self, task, c, fpop, xb, fxb, ki, *args, **dparams):
		r"""Core function of EvolutionStrategyMpL algorithm.

		Args:
			task (Task): Optimization task.
			c (numpy.ndarray): Current population.
			fpop (numpy.ndarray): Current populations fitness/function values.
			xb (numpy.ndarray): Global best individual.
			fxb (float): Global best individuals fitness/function value.
			ki (int): Number of successful mutations.
			args (list): Additional arguments.
			dparams (dict): Additional keyword arguments.

		Returns:
			Tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, float, list, dict]:
				1. New population.
				2. New populations function/fitness values.
				3. New global best solution.
				4. New global best solutions fitness/objective value.
				5. Additional arguments.
				6. Additional keyword arguments:
					* ki (int): Number of successful mutations.
		"""
		if task.Iters % self.k == 0: _, ki = self.update_rho(c, ki), 0
		cn = objects2array([IndividualES(x=self.mutate_rand(c, task), task=task, rnd=self.Rand) for _ in range(self.lam)])
		cn = np.append(cn, c)
		cn = objects2array([cn[i] for i in np.argsort([i.f for i in cn])[:self.mu]])
		ki += self.change_count(c, cn)
		fcn = np.asarray([x.f for x in cn])
		xb, fxb = self.get_best(cn, fcn, xb, fxb)
		return cn, fcn, xb, fxb, args, {'ki': ki}

class EvolutionStrategyML(EvolutionStrategyMpL):
	r"""Implementation of (mu, lambda) evolution strategy algorithm. Algorithm is good for dynamic environments. Mu individual create lambda chields. Only best mu chields go to new generation. Mu parents are discarded.

	Algorithm:
		(:math:`\mu + \lambda`) Evolution Strategy Algorithm

	Date:
		2018

	Authors:
		Klemen Berkovič

	License:
		MIT

	Reference URL:

	Reference paper:

	Attributes:
		Name (List[str]): List of strings representing algorithm names

	See Also:
		* :class:`WeOptPy.algorithm.EvolutionStrategyMpL`
	"""
	Name = ['EvolutionStrategyML', 'EvolutionStrategy(mu,lambda)', 'ES(m,l)']

	def new_pop(self, pop):
		r"""Return new population.

		Args:
			pop (numpy.ndarray): Current population.

		Returns:
			numpy.ndarray: New population.
		"""
		pop_s = np.argsort([i.f for i in pop])
		if self.mu < self.lam: return objects2array([pop[i] for i in pop_s[:self.mu]])
		npop = list()
		for i in range(int(ceil(float(self.mu) / self.lam))): npop.extend(pop[:self.lam if (self.mu - i * self.lam) >= self.lam else self.mu - i * self.lam])
		return objects2array(npop)

	def init_population(self, task):
		r"""Initialize starting population.

		Args:
			task (Task): Optimization task.

		Returns:
			Tuple[numpy.ndarray, numpy.ndarray, list, dict]:
				1. Initialized population.
				2. Initialized populations fitness/function values.
				3. Additional arguments.
				4. Additional keyword arguments.

		See Also:
			* :func:`WeOptPy.algorithm.EvolutionStrategyMpL.initPopulation`
		"""
		return EvolutionStrategyMpL.init_population(self, task)

	def run_iteration(self, task, c, fpop, xb, fxb, *args, **dparams):
		r"""Core function of EvolutionStrategyML algorithm.

		Args:
			task (Task): Optimization task.
			c (numpy.ndarray): Current population.
			fpop (numpy.ndarray): Current population fitness/function values.
			xb (numpy.ndarray): Global best individual.
			fxb (float): Global best individuals fitness/function value.
			args (list): Additional arguments.
			dparams (dict): Additional keyword arguments.

		Returns:
			Tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, float, list, dict]:
				1. New population.
				2. New populations fitness/function values.
				3. New global best solution.
				4. New global best solutions fitness/objective value.
				5. Additional arguments.
				6. Additional keyword arguments.
		"""
		cn = objects2array([IndividualES(x=self.mutate_rand(c, task), task=task, rand=self.Rand) for _ in range(self.lam)])
		c = self.new_pop(cn)
		fc = np.asarray([x.f for x in c])
		xb, fxb = self.get_best(c, fc, xb, fxb)
		return c, fc, xb, fxb, args, dparams

def CovarianceMaatrixAdaptionEvolutionStrategyF(task, epsilon=1e-20, rnd=rand):
	lam, alpha_mu, hs, sigma0 = (4 + np.round(3 * np.log(task.D))) * 10, 2, 0, 0.3 * task.range()
	mu = int(np.round(lam / 2))
	w = np.log(mu + 0.5) - np.log(range(1, mu + 1))
	w = w / np.sum(w)
	mueff = 1 / np.sum(w ** 2)
	cs = (mueff + 2) / (task.D + mueff + 5)
	ds = 1 + cs + 2 * max(np.sqrt((mueff - 1) / (task.D + 1)) - 1, 0)
	ENN = np.sqrt(task.D) * (1 - 1 / (4 * task.D) + 1 / (21 * task.D ** 2))
	cc, c1 = (4 + mueff / task.D) / (4 + task.D + 2 * mueff / task.D), 2 / ((task.D + 1.3) ** 2 + mueff)
	cmu, hth = min(1 - c1, alpha_mu * (mueff - 2 + 1 / mueff) / ((task.D + 2) ** 2 + alpha_mu * mueff / 2)), (1.4 + 2 / (task.D + 1)) * ENN
	ps, pc, C, sigma, M = np.full(task.D, 0.0), np.full(task.D, 0.0), np.eye(task.D), sigma0, np.full(task.D, 0.0)
	x = rnd.uniform(task.lower(), task.upper())
	x_f = task.eval(x)
	while not task.stop_cond_i():
		pop_step = np.asarray([rnd.multivariate_normal(np.full(task.D, 0.0), C) for _ in range(int(lam))])
		pop = np.asarray([task.repair(x + sigma * ps, rnd) for ps in pop_step])
		pop_f = np.apply_along_axis(task.eval, 1, pop)
		isort = np.argsort(pop_f)
		pop, pop_f, pop_step = pop[isort[:mu]], pop_f[isort[:mu]], pop_step[isort[:mu]]
		if pop_f[0] < x_f: x, x_f = pop[0], pop_f[0]
		M = np.sum(w * pop_step.T, axis=1)
		ps = solve(chol(C).conj() + epsilon, ((1 - cs) * ps + np.sqrt(cs * (2 - cs) * mueff) * M + epsilon).T)[0].T
		sigma *= np.exp(cs / ds * (norm(ps) / ENN - 1)) ** 0.3
		ifix = np.where(sigma == np.inf)
		if np.any(ifix): sigma[ifix] = sigma0
		if norm(ps) / np.sqrt(1 - (1 - cs) ** (2 * (task.Iters + 1))) < hth: hs = 1
		else: hs = 0
		delta = (1 - hs) * cc * (2 - cc)
		pc = (1 - cc) * pc + hs * np.sqrt(cc * (2 - cc) * mueff) * M
		C = (1 - c1 - cmu) * C + c1 * (np.tile(pc, [len(pc), 1]) * np.tile(pc.reshape([len(pc), 1]), [1, len(pc)]) + delta * C)
		for i in range(mu): C += cmu * w[i] * np.tile(pop_step[i], [len(pop_step[i]), 1]) * np.tile(pop_step[i].reshape([len(pop_step[i]), 1]), [1, len(pop_step[i])])
		E, V = eig(C)
		if np.any(E < epsilon):
			E = np.fmax(E, 0)
			C = lstsq(V.T, np.dot(V, np.diag(E)).T, rcond=None)[0].T
	return x, x_f

class CovarianceMatrixAdaptionEvolutionStrategy(Algorithm):
	r"""Implementation of (mu, lambda) evolution strategy algorithm. Algorithm is good for dynamic environments. Mu individual create lambda chields. Only best mu chields go to new generation. Mu parents are discarded.

	Algorithm:
		(:math:`\mu + \lambda`) Evolution Strategy Algorithm

	Date:
		2018

	Authors:
		Klemen Berkovič

	License:
		MIT

	Reference URL:
		https://arxiv.org/abs/1604.00772

	Reference paper:
		Hansen, Nikolaus. "The CMA evolution strategy: a tutorial." arXiv preprint arXiv:1604.00772 (2016).

	Attributes:
		Name (List[str]): List of names representing algorithm names
		epsilon (float): TODO
	"""
	Name = ['CovarianceMatrixAdaptionEvolutionStrategy', 'CMA-ES', 'CMAES']
	epsilon = 1e-20

	@staticmethod
	def type_parameters():
		r"""Return functions for checking values of parameters.

		Return:
			Dict[str, Callable[[Any], bool]]:
				* n: Check if number of individuals is :math:`\in [0, \infty]`.
		"""
		return {
			'epsilon': lambda x: isinstance(x, (float, int)) and 0 < x < 1
		}

	def set_parameters(self, epsilon=1e-20, **ukwargs):
		r"""Set core parameters of CovarianceMatrixAdaptionEvolutionStrategy algorithm.

		Args:
			epsilon (float): Small number.
			ukwargs (Dict[str, Any]): Additional arguments.
		"""
		Algorithm.set_parameters(self, **ukwargs)
		self.epsilon = epsilon

	def run_task(self, task):
		r"""TODO.

		Args:
			task (Task): Optimization task.

		Returns:
			TODO.
		"""
		return CovarianceMaatrixAdaptionEvolutionStrategyF(task, self.epsilon, rnd=self.Rand)

# vim: tabstop=3 noexpandtab shiftwidth=3 softtabstop=3

reason: import numpy,from numpy


text:
import math
import numpy as np

def fill_hist(h_bins, h_low, h_up, values, dt_offset):
  bin_width = (float(h_up) - float(h_low)) / float(h_bins)
  hist = [0.0] * (h_bins + 2) #add 2 extra bins for underflow and overflow

  for i in range(0, len(values)):
    v = values[i] + dt_offset
    filled = False

    for ii in range(0, h_bins):
      if v >= (ii * bin_width + h_low) and v < ((ii+1) * bin_width + h_low):
        hist[ii+1] += 1
        filled = True
        break

    if filled == True:
      continue
    elif v < h_low:
      hist[0] += 1
    elif v >= h_up:
      hist[len(hist)-1] += 1

  hist = [float(x)/float(float(len(values))-hist[0]-hist[-1]) for x in hist] #normalise excluding flow bins

  return hist

def remove_flow(hist): #remove the flow bins
  hist.remove(hist[-1])
  hist.remove(hist[0])

  return hist

def diff_hist(hist1, hist2, scale):
  if len(hist1) != len(hist2):
    print("histograms with different number of bins cannot be merged!")
    exit()

  h1max = max(hist1) * scale
  h2max = max(hist2) * scale
  sum_hist = [0.0] * len(hist1)
  for i in range(len(sum_hist)):
    if hist1[i] != 0 and hist2[i] != 0:
      sum_hist[i] = hist1[i] + hist2[i]
    else:
      sum_hist[i] = h1max + h2max
  diff_hist = [abs(hist2[ii] - hist1[ii]) for ii in range(len(hist1))]

  metric = 0
  for iii in range(len(diff_hist)):
    metric += sum_hist[iii] * diff_hist[iii]

  return metric

def minimise(mlist, dt0, dt_step, dt_N, polyN, fit_range):
  dt_list = [(dt0 + i*dt_step) for i in range(dt_N)]

  #setting fit range
  i_fit_range = int(fit_range / dt_step) + 1
  i_min_metric = np.argmin(mlist)
  i_low = 0
  i_up = dt_N - 1
  if (i_min_metric - i_fit_range) > 0:
    i_low = i_min_metric - i_fit_range
  if (i_min_metric + i_fit_range) < (dt_N - 1):
    i_up = i_min_metric + i_fit_range
  coeff = np.polyfit(dt_list[i_low:i_up], mlist[i_low:i_up], polyN) #polynomial fit, coeff[0]*x^polyN + coeff[1]*x^(polyN-1) + ...

  min_dt = 0
  func_prev = 0
  for ii in range(len(dt_list)):
    func = 0
    for iii in range(polyN):
      func += coeff[iii] * (polyN - iii) * pow(dt_list[ii], polyN-1-iii)

    if ii == 0:
      func_prev = func
      continue
    if func * func_prev <= 0:
      min_dt = (dt_list[ii] + dt_list[ii-1])/2.0
      break
    func_prev = func

  return min_dt

max_stars_repo_name: calhewitt/snewpdag

max_issues_repo_name: calhewitt/snewpdag

max_forks_repo_name: calhewitt/snewpdag


text:
import numpy as np
import math
from gridmap2d import *
from bresenham_algorithm import *
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse, Circle
from matplotlib.animation import ArtistAnimation

def pose2transform(pose):
    """
    @brief convert (x, y , yaw) to transform matrix
    @param pose: (x, y, yaw)
    @return: 3*3 transform matrix
    """
    transform = np.identity(3)
    transform[0:2, 2] = pose[0:2]
    transform[0, 0] = math.cos(pose[2])
    transform[0, 1] = -math.sin(pose[2])
    transform[1, 0] = math.sin(pose[2])
    transform[1, 1] = math.cos(pose[2])
    return transform

def bresenham_line_by_world(world_map, A, B, prob):
    a = world_map.world2pixel(A)
    b = world_map.world2pixel(B)
    res = bresenham_line(a, b)
    for item in res:
        world_map.set_prob_by_pixel(item, prob)

def draw_circle_test(world_map):
    C = np.array([6, 6, 1]) / world_map.resolution
    res = bresenham_circle(C)
    for item in res:
        world_map.set_prob_by_pixel(item, world_map.probrange[1])

#    c(0,0)
#   /\
#  /  \
# /    \
# ------
#d(-0.6,-0.2) e(-0.6, 0.2)
def construct_robot_in_gridmap(pose):
    """
    @brief draw robot

    @param pose: (x, y, yaw) should be a numpy vector
    """
    transform = pose2transform(pose)
    c = np.array([0.0, 0.0, 1.0])
    d = np.array([-0.6, -0.2, 1.0])
    e = np.array([-0.6, 0.2, 1.0])
    c = transform @ c
    d = transform @ d
    e = transform @ e
    bresenham_line_by_world(world_map, e, d, world_map.probrange[1])
    bresenham_line_by_world(world_map, c, d, world_map.probrange[1])
    bresenham_line_by_world(world_map, c, e, world_map.probrange[1])

def construct_home(world_map):
    A = np.array([0, 0]) + 1
    B = np.array([0, 10]) + 1
    bresenham_line_by_world(world_map, A, B, world_map.probrange[1])
    A = np.array([0, 10]) + 1
    B = np.array([8, 10]) + 1
    bresenham_line_by_world(world_map, A, B, world_map.probrange[1])
    A = np.array([8, 10]) + 1
    B = np.array([8, 0]) + 1
    bresenham_line_by_world(world_map, A, B, world_map.probrange[1])
    A = np.array([8, 0]) + 1
    B = np.array([0, 0]) + 1
    bresenham_line_by_world(world_map, A, B, world_map.probrange[1])
    A = np.array([4, 0]) + 1
    B = np.array([4, 5]) + 1
    bresenham_line_by_world(world_map, A, B, world_map.probrange[1])
    A = np.array([4, 8]) + 1
    B = np.array([4, 10]) + 1
    bresenham_line_by_world(world_map, A, B, world_map.probrange[1])
    A = np.array([3, 6]) + 1
    B = np.array([3, 8]) + 1
    bresenham_line_by_world(world_map, A, B, world_map.probrange[1])
    A = np.array([0, 5]) + 1
    B = np.array([3, 5]) + 1
    bresenham_line_by_world(world_map, A, B, world_map.probrange[1])
    A = np.array([5, 5]) + 1
    B = np.array([8, 5]) + 1
    bresenham_line_by_world(world_map, A, B, world_map.probrange[1])
    A = np.array([0, 8]) + 1
    B = np.array([3, 8]) + 1
    bresenham_line_by_world(world_map, A, B, world_map.probrange[1])

ell1 = Ellipse(xy = (0.0, 0.0), width = 4, height = 8, angle = 30.0, facecolor= 'yellow', alpha=0.3)
print(ell1.get_gid())

world_map = gridmap2d((10, 12), 0.02, (-20, 120))
images = []
fig, ax = plt.subplots()
construct_home(world_map)

pose = np.array([3, 3, math.pi / 2])
construct_robot_in_gridmap(pose)

draw_circle_test(world_map)

im = plt.imshow(world_map.mapdata, cmap=plt.cm.get_cmap('hot'), animated=True,
    vmin=world_map.probrange[0], vmax=world_map.probrange[1])

fig.colorbar(im)
images.append([im])

ani = ArtistAnimation(fig, images, interval=1, blit=True,
                                    repeat_delay=10)
plt.show()

max_stars_repo_name: democheng/PythonRobotics

max_issues_repo_name: democheng/PythonRobotics

max_forks_repo_name: democheng/PythonRobotics


text:
#!/usr/bin/env python3
"""Finds latest versions of the CSVs for each subsystem, then plots the time
domain and X-Y data.

If provided, the first argument to this script is a filename regex that
restricts which CSVs are plotted to those that match the regex.
"""

import argparse
import matplotlib.pyplot as plt
import numpy as np
import os
import re

class UnitGroup:

    def __init__(self):
        # List of DataSeries objects (time and data column pairs)
        self.series = []

        # List of data labels
        self.labels = []

class NameGroup:

    def __init__(self, filename, series):
        self.filename = filename
        self.series = series

class DataSeries:

    def __init__(self, time, data):
        self.time = time
        self.data = data

def num_lines(filename):
    with open(filename) as f:
        i = 0
        for i, l in enumerate(f):
            pass
    return i + 1

def get_file_list(regex):
    # Get list of files in current directory
    files = [
        os.path.join(dp, f)
        for dp, dn, fn in os.walk(".")
        for f in fn
        if f.endswith(".csv")
    ]

    # Ignore files not matching optional pattern
    if regex:
        files = [f for f in files if re.search(regex, f)]

    # Maps subsystem name to tuple of csv_group and date (group is set of files
    # with states, inputs, or outputs suffixes and the same name stub like
    # "Flywheel")
    filtered = {}
    file_rgx = re.compile(
        r"(?P<name>^\./.*?[A-Za-z ]+)-(?P<date>\d{4}-\d{2}-\d{2}-\d{2}_\d{2}_\d{2})\.csv$"
    )
    for f in files:
        match = file_rgx.search(f)
        if not match:
            continue

        # If file is empty or only has header (that is, has no data), ignore it.
        # We ignore the case of one line of data because it might be truncated.
        if num_lines(f) <= 2:
            continue

        # If the file is a CSV with the correct name pattern, add it to the
        # filtered list. Files with newer dates override old ones in lexographic
        # ordering.
        name = match.group("name")
        date = match.group("date")
        if name not in filtered.keys() or filtered[name] < date:
            filtered[name] = date

    # Make filtered list of files
    files = []
    for name_stub in filtered.keys():
        files.append(name_stub + "-" + filtered[name_stub] + ".csv")

    return files

def make_groups(files):
    # Group files by category (sets of files with states, inputs, or outputs
    # suffixes and the same name stub like "Flywheel")
    category_rgx = re.compile(
        r"^\./.*?(?P<category>[A-Za-z ]+) (states|inputs|outputs)-")
    file_groups = {}
    if files:
        print("Loading CSVs...")
    else:
        print("No data to plot.")

    # Sorting the file list puts files into the order ["inputs", "outputs",
    # "states"]. This means data series will be loaded in the order of
    # ["inputs", "outputs", "references", "states"] (references are logged
    # before states). This produces the desired dataset layering on plots.
    for f in sorted(files):
        print(f"  {os.path.split(f)[1]}")

        match = category_rgx.search(f)

        if not match:
            # Couldn't find the file's category, so put it in its own category
            substr = os.path.split(f)[1]
            substr = substr[:substr.find("-")]
            file_groups[substr] = [f]
        else:
            category = match.group("category")

            # Create a new category if one doesn't exist, or add the file to the
            # existing category if it does
            if category not in file_groups.keys():
                file_groups[category] = [f]
            else:
                file_groups[category].append(f)
    return file_groups

def main():
    plt.rcParams.update({'figure.max_open_warning': 0})

    parser = argparse.ArgumentParser()
    parser.add_argument("-ymin",
                        dest="ymin",
                        type=float,
                        help="Y minimum for plots")
    parser.add_argument("-ymax",
                        dest="ymax",
                        type=float,
                        help="Y maximum for plots")
    parser.add_argument("regex", nargs="?")
    args = parser.parse_args()

    file_groups = make_groups(get_file_list(args.regex))
    if file_groups:
        print("Plotting...")

    # Within each group, make groups of datasets keyed on their unit, then plot
    # each group on their own figure
    unit_rgx = re.compile(r"^(?P<name>[\w\- ]+) \((?P<unit>.*?)\)$")
    for category, file_group in file_groups.items():
        unit_groups = {}
        name_groups = {}
        for filename in file_group:
            # Get labels from first row of file
            with open(filename) as f:
                labels = [
                    x.strip('"') for x in f.readline().rstrip().split(",")
                ]

            # Retrieve data from remaining rows of file. "skip_footer=1" skips
            # the last line because it may be incompletely written.
            data = np.genfromtxt(filename,
                                 delimiter=",",
                                 skip_header=1,
                                 skip_footer=1)

            times = data[:, 0:1]

            # Skips label in first column because that's always "Time (s)"
            for i, label in enumerate(labels[1:]):
                match = unit_rgx.search(label)
                name = match.group("name")
                unit = match.group("unit")

                if unit not in unit_groups.keys():
                    unit_groups[unit] = UnitGroup()
                # "i + 1" skips the time data column
                unit_groups[unit].series.append(
                    DataSeries(times, data[:, i + 1:i + 2]))
                unit_groups[unit].labels.append(name)

                # "i + 1" skips the time data column
                name_groups[name] = NameGroup(filename, data[:, i + 1:i + 2])

        # Plot time domain datasets
        print(f'  [vs time] {category} ({", ".join(unit_groups.keys())})')
        for unit, unit_group in unit_groups.items():
            fig, ax = plt.subplots(1, 1)
            ax.set_title(f"{category} ({unit})")

            for i in range(len(unit_group.series)):
                ax.plot(unit_group.series[i].time, unit_group.series[i].data)
                if args.ymin:
                    ax.set_ylim(bottom=args.ymin)
                if args.ymax:
                    ax.set_ylim(top=args.ymax)

            # First label is x axis label (time). The remainder are dataset
            # names.
            ax.set_xlabel("Time (s)")
            ax.set_ylabel(f"Data ({unit})")
            ax.legend(unit_group.labels)

        # Plot X-Y datasets. If the file doesn't have all the required keys,
        # skip it.
        if not (set(["X reference", "Y reference", "X estimate", "Y estimate"])
                - set(name_groups.keys())):
            print(f'  [y vs x] {category}')
            fig, ax = plt.subplots(1, 1)
            ax.set_title(f"{category} trajectory")

            ax.plot(name_groups["X reference"].series,
                    name_groups["Y reference"].series)
            ax.plot(name_groups["X estimate"].series,
                    name_groups["Y estimate"].series)

            ax.set_xlabel("X (m)")
            ax.set_ylabel("Y (m)")
            ax.legend(["Reference", "Estimate"])

            # This equalizes the X and Y axes so the trajectories aren't warped
            ax.axis("equal")

            if any([
                    "AutoNav" in group.filename or
                    "Galactic Search" in group.filename
                    for group in name_groups.values()
            ]):
                img = plt.imread("tools/images/2021-challenge-space.png")
                ax.imshow(img, extent=[0, 9.144, 0, 4.572])
            else:
                img = plt.imread("tools/images/2021-field.png")
                ax.imshow(img, extent=[0, 15.98295, 0, 8.21055])

    plt.show()

if __name__ == "__main__":
    main()


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#
#  IkaLog
#  ======
#  Copyright (C) 2015 Takeshi HASEGAWA
#
#  Licensed under the Apache License, Version 2.0 (the "License");
#  you may not use this file except in compliance with the License.
#  You may obtain a copy of the License at
#
#      http://www.apache.org/licenses/LICENSE-2.0
#
#  Unless required by applicable law or agreed to in writing, software
#  distributed under the License is distributed on an "AS IS" BASIS,
#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#  See the License for the specific language governing permissions and
#  limitations under the License.
#
import sys
import time

import cv2
import numpy as np

from ikalog.constants import special_weapons
from ikalog.scenes.stateful_scene import StatefulScene
from ikalog.utils import *
from ikalog.utils.character_recoginizer import *
from ikalog.utils.ikamatcher2.matcher import MultiClassIkaMatcher2 as MultiClassIkaMatcher

##
# GameSpecialWeapon
#
# Detects if any special weapons are activated.
#
# EVENTS:
#
#     on_game_special_weapon()
#     Notifies special weapon activation in my team.
#         context['game']['special_weapon']
#             Id of special_weapon activated.
#
# HOW TO MANIPULATE THE MASKS:
#
#     This scene needs masks for each languages. To gather mask data,
#     Set "write_sample" in __init__() function to True.
#     It will write a sample image.
#
#     Filename of the samples should be:
#       masks/<IKALOG_LANG>/special_<SPECIAL_WEAPON_ID>.png
#
#     For SPECIAL_WEAPON_ID, check the array definition of
#      "special_weapon" array in ikalog/constants.py.
#
class GameSpecialWeapon(StatefulScene):

    # Called per Engine's reset.
    def reset(self):
        super(GameSpecialWeapon, self).reset()
        self.img_last_special = None

    def _match_phase1(self, context, img_special, img_last_special):
        #
        # Phase 1
        #
        # Crop the area special weapon message supposed to be appeared.
        # Compare with last frame, and check if it is (almost) same with
        # the last frame.
        #

        img_special_diff = abs(img_special - img_last_special)
        matched = bool(np.average(img_special_diff) < 90)
        return matched

    def _is_my_special_weapon(self, context, img_special_bgr):
        img = img_special_bgr[:, :150]

        img_s = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)[:, :, 2]
        img_s[matcher.MM_WHITE()(img) > 127] = 127

        img_s_hist = cv2.calcHist(img_s[:, :], [0], None, [5], [0, 256])
        img_s_hist_black = float(np.amax(img_s_hist[0:1]))
        img_s_hist_non_black = float(np.amax(img_s_hist[3:4]))
        return img_s_hist_black < img_s_hist_non_black

    def _state_default(self, context):
        if not self.is_another_scene_matched(context, 'GameTimerIcon'):
            return False

        frame = context['engine']['frame']
        if frame is None:
            return False

        # FIXME: this code works with the first special weapon only
        img_special_bgr = frame[260:260 + 24, 1006:1006 + 210, :]
        img_special = np.array(img_special_bgr)
        img_last_special = self.img_last_special
        self.img_last_special = img_special

        if img_last_special is None:
            return False

        if not self._match_phase1(context, img_special, img_last_special):
            return False

        #
        # Phase 2
        #
        # Check inkling image on right side
        img_sp_char = cv2.cvtColor(
            img_special_bgr[:, 150:210], cv2.COLOR_BGR2GRAY)

        laplacian_threshold = 60
        img_laplacian = cv2.Laplacian(img_sp_char, cv2.CV_64F)
        img_laplacian_abs = cv2.convertScaleAbs(img_laplacian)
        c_matched = bool(np.average(img_laplacian_abs) > 20)

        if not c_matched:
            return False

        # Phase 3
        # TODO: Background color

        # Phase 4
        # Forground text
        white_filter = matcher.MM_WHITE()
        img_sp_text = white_filter(img_special_bgr[:, 0:150])

        special = self._masks.match_best(img_special_bgr)[1]

        if self.write_samples:
            cv2.imwrite('training/_special_%s.png' %
                        time.time(), 255 - img_sp_text)

        if special is None:
            return False

        context['game']['special_weapon'] = special._id
        context['game']['special_weapon_is_mine'] = \
            self._is_my_special_weapon(context, img_special_bgr)
        self._call_plugins('on_game_special_weapon')

        self._switch_state(self._state_tracking)
        return True

    def _state_tracking(self, context):
        if not self.is_another_scene_matched(context, 'GameTimerIcon'):
            return False

        frame = context['engine']['frame']
        if frame is None:
            return False

        # FIXME
        img_special_bgr = frame[260:260 + 24, 1006:1006 + 210, :]
        img_special = np.array(img_special_bgr)
        img_last_special = self.img_last_special

        if img_last_special is None:
            return False

        special = self._masks.match_best(img_special_bgr)[1]
        if special is not None:
            self._call_plugins(
                'on_mark_rect_in_preview',
                [(1006, 260), (1006 + 210, 260 + 24)]
            )

            if context['game']['special_weapon'] == special._id:
                return True

        if self.matched_in(context, 150):
            return False

        self._switch_state(self._state_default)
        self.img_last_special = None
        return False

    def dump(self, context):
        # Not implemented :\
        pass

    def _analyze(self, context):
        pass

    # Called only once on initialization.
    def _init_scene(self, debug=False):
        #
        # To gather mask data, enable this.
        #
        self.write_samples = False

        # Load mask files.
        self._masks = MultiClassIkaMatcher()
        for special_weapon in special_weapons.keys():
            try:
                mask = IkaMatcher(
                    0, 0, 150, 24,
                    img_file='special_%s.png' % special_weapon,
                    threshold=0.90,
                    orig_threshold=0.20,
                    bg_method=matcher.MM_NOT_WHITE(),
                    fg_method=matcher.MM_WHITE(),
                    label='special/%s' % special_weapon,
                    call_plugins=self._call_plugins,
                    debug=debug,
                )
                mask._id = special_weapon
                self._masks.add_mask(mask)
            except:
                IkaUtils.dprint('%s: Failed to load mask for %s' %
                                (self, special_weapon))
                pass

if __name__ == "__main__":
    GameSpecialWeapon.main_func()

max_stars_repo_name: fetus-hina/IkaLog

max_issues_repo_name: fetus-hina/IkaLog

max_forks_repo_name: fetus-hina/IkaLog


text:
#! /usr/bin/env python
# -*- mode: python; coding: utf-8 -*-
# Copyright 2018 Peter Williams and collaborators
# Licensed under the MIT License.

"""More forward-time steady-state experimentation. Write:

   df/dt = V df/dx + Dxx d²f/dx²

for advection velocity V, diffusion constant Dxx. Set boundaries

   f(0) = 1,
   f(1) = 0.

Writing g = df/dx and solving for the steady-state, I find that after
defining Q = -V / Dxx, the steady-state solution is:

  f(x) = (exp(Q) - exp(Qx)) / (exp(Q) - 1)

Once again the generic time-forward SDE is:

  df/dt = -d(a f)/dx + 1/2 d²(b² f)/dx²

Since our constants do not vary spatially, it is easy to find that:

  a = -V    # note sign!
  b = sqrt(2 Dxx)

"""

import numpy as np
import omega as om

XMAX = 1.
Dxx = 0.25
V = 1.0
Q = -V / Dxx

def calculate(n_pseudo_particles=60000, delta_t=0.005):
    n_bins = 20
    x_edges = np.linspace(0., XMAX, n_bins + 1)
    x_centers = 0.5 * (x_edges[1:] + x_edges[:-1])
    grid = np.zeros(n_bins)
    a = -V
    b = np.sqrt(2 * Dxx)
    x0 = 0.

    #length_scale = 0.02
    #dt_b = 0.1 * length_scale**2 / b**2
    #dt_a = 0.1 * b**2 / a**2
    #print('delta_t:', dt_b, dt_a)

    for num in range(n_pseudo_particles):
        x = x0
        t_scale = 1.0

        while True:
            this_delta_t = delta_t * np.random.normal(loc=1, scale=0.2)
            if this_delta_t <= 0:
                this_delta_t = 0.1 * delta_t

            x += a * this_delta_t + b * np.sqrt(this_delta_t) * np.random.normal()
            if x <= 0 or x >= XMAX:
                break

            loc = min(int(x * n_bins / XMAX), n_bins - 1)
            grid[loc] += t_scale * this_delta_t
            t_scale = 1.

        #while x >= 0 and x <= XMAX:
        #    this_delta_t = delta_t * np.random.normal(loc=1, scale=0.2)
        #    if this_delta_t <= 0:
        #        this_delta_t = 0.1 * delta_t
        #
        #    loc = min(int(x * n_bins / XMAX), n_bins - 1)
        #    grid[loc] += t_scale * this_delta_t
        #    x += a * this_delta_t + b * np.sqrt(this_delta_t) * np.random.normal()
        #    t_scale = 1.

    grid /= grid[0] # close enough for jazz
    grid *= 0.9 # XXX hack to line up with analytic given binning
    return x_centers, grid

def main():
    import time
    from pwkit.ndshow_gtk3 import view

    t0 = time.time()
    x_centers1, grid1 = calculate(delta_t=0.0002)
    #x_centers2, grid2 = calculate(delta_t=0.00005)
    elapsed = time.time() - t0
    print('Calculated grid(s) in %.1f seconds' % elapsed)

    xex = np.linspace(0, XMAX, 50)
    yex = (np.exp(Q) - np.exp(Q * xex)) / (np.exp(Q) - 1)

    p = om.RectPlot()
    p.addXY(xex, yex, 'exact')
    p.addXY(x_centers1, grid1, '0.0002')
    #p.addXY(x_centers2, grid2, '0.00005')
    p.defaultKeyOverlay.hAlign = 0.95
    p.show()

if __name__ == '__main__':
    main()

max_stars_repo_name: pkgw/vernon

max_issues_repo_name: pkgw/vernon

max_forks_repo_name: pkgw/vernon


text:
import pandas as pd
import numpy as np
import re
from collections import defaultdict
import collections
from data.definition import LIVER_FINDINGS
from data.definition import TEMPLE_RADIO_1_2_SENTENCE_EMBEDDINGS, TEMPLE_RADIO_1_2_SENTENCES_WITH_TRANSLATION, DUPL_SENTENCES

from sentence_transformers import SentenceTransformer

def generate_clinical_bert_representation():
    # Use huggingface/transformers pre-trained model Bio_ClinicalBERT for mapping tokens to embeddings
    model = SentenceTransformer("emilyalsentzer/Bio_ClinicalBERT")

    first_sentences = pd.read_csv(LIVER_FINDINGS, header=0, usecols =["Sentence 1", "Sentence 1 Translation"])
    second_sentences = pd.read_csv(LIVER_FINDINGS, header=0, usecols =["Sentence 2", "Sentence 2 Translation"])
    first_sentences.columns = ["Sentence", "Translation"]
    second_sentences.columns = ["Sentence", "Translation"]

    df = pd.concat([first_sentences, second_sentences]).reset_index(drop=True)
    df = df[df['Sentence'].notna()]
    df['Translation'].replace(np.nan, 'empty', regex=True, inplace=True)
    df = df.drop_duplicates(subset=['Sentence'], keep='first')

    df['Sentence'] = [sent if sent.endswith('.') else sent + '.' for sent in df['Sentence'].values]
    df.Translation = [sent.strip() for sent in df.Translation.values]
    df['Sentence'] = [re.sub('(\d{1,2})/(\d{1,2})/(\d{4})', '', sent) for sent in df['Sentence'].values]
    df['Sentence'] = [re.sub('[0-9]+.[0-9]+', '', sent) for sent in df['Sentence'].values]
    df['Sentence'] = [re.sub('[0-9]+', '', sent) for sent in df['Sentence'].values]
    df['Sentence'] = [sent.replace('  ', ' ') for sent in df['Sentence'].values]
    df['Sentence'] = [sent.strip() for sent in df['Sentence'].values]

    df.to_csv(TEMPLE_RADIO_1_2_SENTENCES_WITH_TRANSLATION, index=False)

    df['Sentence'] = [sent.lower() for sent in df['Sentence'].values]
    sentence_embeddings = model.encode(df.Sentence.to_list())

    np.savetxt(TEMPLE_RADIO_1_2_SENTENCE_EMBEDDINGS, sentence_embeddings, delimiter=",")

if __name__ == '__main__':
    generate_clinical_bert_representation()


text:
import numpy as np
import libyay
a = np.array([1.0])
print libyay.yay(a)
print a

max_stars_repo_name: HTuennermann/pythonmodule

max_issues_repo_name: HTuennermann/pythonmodule

max_forks_repo_name: HTuennermann/pythonmodule


text:
import cv2
import numpy as np

net= cv2.dnn.readNet("yolov3_1000.weights", "yolov3.cfg" )

classes= []
with open("obj.names", "r") as f:
    classes = [line.strip() for line in f.readlines()]
# print(classes)

layer_names=net.getLayerNames()
output_layers = [layer_names[i[0]-1] for i in net.getUnconnectedOutLayers()]

img = cv2.imread("m148_2.jpg")
height, width, channels = img.shape

blob=cv2.dnn.blobFromImage(img, 0.00392, (416, 416), (0, 0, 0), True, crop=False)

# for b in blob:
# 	for n, img_blob in enumerate(b):
# 		cv2.imshow(str(n), img_blob)

net.setInput(blob)
outs = net.forward(output_layers)

class_ids = []
confidences = []
boxes = []
for out in outs:
	for detection in out:
		scores = detection[5:]
		class_id = np.argmax(scores)
		confidence = scores[class_id]
		if confidence > 0.5:
			center_x = int(detection[0] * width)
			center_y = int(detection[1] * height)
			w = int(detection[2] * width)
			h = int(detection[3] * height)

			# cv2.circle(img, (center_x, center_y), 10, (0, 255, 0), 2)
			x = int(center_x - w/2)
			y = int(center_y - h/2)

			boxes.append([x,y,w,h])
			confidences.append(float(confidence))
			class_ids.append(class_id)

			# cv2.rectangle(img, (x, y), (x + w, y + h), (0,255,0), 2)

indexes = cv2.dnn.NMSBoxes(boxes, confidences, 0.5, 0.4)
print(indexes)

font = cv2.FONT_HERSHEY_PLAIN
number_objects_detected = len(boxes)
for i in range(len(boxes)):
	if i in indexes:
		x, y, w, h = boxes[i]
		label = str(classes[class_ids[i]])
		cv2.rectangle(img, (x, y), (x + w, y + h), (0, 255, 0), 2)
		# cv2.putText(img, label, (x, y + 30), font, 1, (0, 0, 0), 3)

img = cv2.imshow("Image", img)
cv2.waitKey(0)

max_stars_repo_path: Detection using Yolo/detect.py

max_issues_repo_path: Detection using Yolo/detect.py

max_forks_repo_path: Detection using Yolo/detect.py