text:
from __future__ import print_function
import os, sys, glob
import numpy as np
import cv2
import scipy.spatial
import matplotlib
import wavepy
import matplotlib.pyplot as plt
# import pyvista as pv
# import plotly.graph_objects as go

# Set Parameters
DEFLECTOMETRY_GRADMAP_FILTERSIZE = 49

# Render parameters
MESH_Ka = [0.0, 0.0, 0.0]
MESH_Kd = [0.9, 0.9, 0.9]
MESH_Ks = [0.1, 0.1, 0.1]
MESH_Ns = 5.000000
MESH_d = 0.010000
MESH_illum = 1

# Scaling factors
MESH_DEPTH_SCALINGFACTOR = 5
MESH_SIGMA_GLOBCORR = 5

class Mesh:
    def __init__(self, name, height, width, cropMask=0):

        self.name = name
        self.height = height
        self.width = width
        self.vertex = np.zeros((self.height * self.width, 3))
        # set 2d plane
        self.setVertex2D()
        self.setFace()
        if cropMask == 0:
            self.mask = np.array([0, 0, self.height, self.width])
        else:
            self.mask = cropMask

    def setVertex2D(self):
        [i, j] = np.meshgrid(np.r_[1:self.height + 1], np.r_[1:self.width + 1])
        self.vertex2D = np.concatenate((i.flatten('F')[:, None], j.flatten('F')[:, None]), axis=1).astype(np.int16)
        self.vertex[:, 0:2] = self.vertex2D

    # options = 'Qt QbB Qc'
    # self.faces = scipy.spatial.Delaunay(v, qhull_options=options).simplices + 1

    def setFace(self):
        # N = 2
        options = 'Qt QbB Qc'
        self.faces = scipy.spatial.Delaunay(self.vertex2D, qhull_options=options).simplices + 1

    def setDepth(self):
        # Set x and y direction
        p = self.normals[..., 0] / self.normals[..., 2]
        q = self.normals[..., 1] / self.normals[..., 2]
        p = np.nan_to_num(p)
        q = np.nan_to_num(q)
        # TODO: Low pass filter to filter out the noise (BilateralFilter)
        # Compute depth as a imaginary number (real: depth, img: noise)
        # Using Frankot Chellappa method with FFT
        #x = sfg.frankotchellappa(p, q, False)
        x = wavepy.surface_from_grad.frankotchellappa(p, q, False)
        depth = x.real
        # depth_smooth = cv2.GaussianBlur(depth, (sigma_GlobCorr, sigma_GlobCorr), ((sigma_GlobCorr + 1) / 2))
        depth *= MESH_DEPTH_SCALINGFACTOR
        # cv.bilateralFilter(src, dst, 9, 75, 75, cv.BORDER_DEFAULT)
        depth_corr = cv2.bilateralFilter(depth.astype(np.float32), MESH_SIGMA_GLOBCORR, 90, 90)
        # depth_corr = depth
        # cv2.imwrite('./depth_before.png', (depth*255).astype(np.uint8))
        # depth = cv2.blur(depth.astype(np.float32), (5,5))
        # cv2.imwrite('./depth_after.png', (depth*255).astype(np.uint8))
        depth_crop = depth_corr[self.mask[0]:self.mask[2], self.mask[1]:self.mask[3]]
        self.vertex[:, 2] = depth_crop.reshape(-1)
        return depth

    def setNormal(self, normals):
        self.normals = normals

    def setTexture(self, texture):
        #self.texture = (cv2.flip(texture, -1) * 255).astype(np.uint8)
        self.texture = cv2.flip(texture, -1)

    # @nb.jit(nopython=True)
    def exportOBJ(self, path, withTexture=True):
        # totalCount = self.vertex.shape[0] + self.height * self.width + self.faces.shape[0]
        # progress = ProgressBar(totalCount, fmt=ProgressBar.FULL)
        print('Exporting OBJ...')
        N = self.normals[self.mask[0]:self.mask[2], self.mask[1]:self.mask[3]].reshape(-1, 3)
        filename = os.path.join(os.path.normpath(path), "mesh_" + self.name + ".obj")
        f = open(filename, 'w')
        f.write('mtllib material_' + self.name + '.mtl\n')
        f.write('usemtl Textured\n')
        for i in range(self.vertex.shape[0]):
            # progress.current += 1
            # progress()
            f.write('v %f %f %f\n' % (self.vertex[i, 0], self.vertex[i, 1], self.vertex[i, 2]))
            f.write('vn %f %f %f\n' % (N[i, 0], N[i, 1], N[i, 2]))

        if withTexture:
            for i in range(self.height):
                for j in range(self.width):
                    # progress.current += 1
                    # progress()
                    u = (i) / self.height
                    v = 1 - (j) / self.width
                    f.write('vt %f %f\n' % (v, u))

        for i in range(self.faces.shape[0]):
            # progress.current += 1
            # progress()
            # f.write('f ' + str(self.faces[i, 0]) + ' ' + str(self.faces[i, 1]) + ' ' + str(self.faces[i, 2]) + '\n')
            f.write(
                'f ' + str(self.faces[i, 0]) + '/' + str(self.faces[i, 0]) + '/' + str(self.faces[i, 0]) + ' ' + str(
                    self.faces[i, 1]) + '/' + str(self.faces[i, 1]) + '/' + str(self.faces[i, 1]) + ' ' + str(
                    self.faces[i, 2]) + '/' + str(self.faces[i, 2]) + '/' + str(self.faces[i, 2]) + '\n')

        f.close()

        if withTexture:
            # export mtl
            filename = os.path.join(os.path.normpath(path), "material_" + self.name + ".mtl")
            f = open(filename, 'w')
            f.write('newmtl Textured\n')
            f.write('Ka %f %f %f\n' % (MESH_Ka[0], MESH_Ka[1], MESH_Ka[2]))
            f.write('Kd %f %f %f\n' % (MESH_Kd[0], MESH_Kd[1], MESH_Kd[2]))
            f.write('Ks %f %f %f\n' % (MESH_Ks[0], MESH_Ks[1], MESH_Ks[2]))
            f.write('Ns %f\n' % (MESH_Ns))
            f.write('d %f\n' % (MESH_d))
            f.write('illum %d\n' % (MESH_illum))
            f.write("map_Kd texture_" + self.name + ".jpg\n")
            f.close()

            # export texture
            filename = os.path.join(os.path.normpath(path), "texture_" + self.name + ".jpg")
            # cv2.imwrite(filename, self.texture)
            # cv2.imwrite(filename, np.array((utilities.imgBrighten(self.texture, 10) * 255), dtype=np.uint8))
            cv2.imwrite(filename, self.texture)

max_issues_repo_name: clkimsdu/GhostScan

max_forks_repo_name: clkimsdu/GhostScan


text:
#!/usr/bin/env python

"""
SYNOPSIS

    process_brightness.py [-h] [--verbose] [--plot] [-v, --version] <filename> <outputfile_base_name>

DESCRIPTION

    TODO This describes how to use this script.
    This docstring will be printed by the script if there is an error or
    if the user requests help (-h or --help).

EXAMPLES

    process_brightness.py IMG0092.m4v

AUTHOR

    Ethan Gutmann - gutmann@ucar.edu

LICENSE

    This script is in the public domain.

REQUIREMENTS

    numpy, scipy, interception
    if plotting : matplotlib
    if saving to netcdf : xarray

VERSION
    1.0

"""
from __future__ import absolute_import, print_function, division

import sys
import os
import traceback
import argparse
import time

import numpy as np
from scipy.signal import tukey
from scipy.ndimage.filters import median_filter as medfilt
from matplotlib import pyplot as plt

from interception import video_reader as vr

global verbose
verbose = False

def filter_best(data, ampl):
    ny, nx = data.shape
    output = np.zeros(data.shape, dtype=np.int)
    for i in range(1,ny-2):
        for j in range(1,nx-2):
            pos = np.argmax(ampl[i-1:i+2, j-1:j+2])
            output[i,j] = data[i-1:i+2, j-1:j+2].flat[pos]
    return medfilt(output, (3,3)).astype('i')

def main (filename, output_file,
          plot_data = False,
          save_results = False,
          normalize = True,
          min_amplitude = 0.001,
          min_freq = 0.2,
          max_freq = 2.5,
          reduction = 2,
          max_times = 50,
          start_time = 2):

    file_base = output_file # filename.split(".")[0]
    prefix = "t{}_{}_".format(max_times, start_time)

    if os.path.isfile("{prefix}amplitudes_{}.png".format(file_base, prefix=prefix)):
        if verbose: print("Already processed: "+filename)
        sys.exit()

    video = vr.Video_Reader(filename)

    nx = int(video.shape[2]/reduction)
    ny = int(video.shape[1]/reduction)

    if verbose: print("loading data")
    t0 = time.time()

    frame_rate = video.metadata["r_frame_rate"]
    fps = float(frame_rate.split("/")[0]) / float(frame_rate.split("/")[1])
    if verbose: print("FPS:{}".format(fps))
    max_times *= fps
    start_time *= fps
    max_times = int(max_times)
    start_time = int(start_time)

    data = np.zeros((min(video.shape[0],max_times), ny, nx), dtype=np.float)

    if start_time > 0:
        for i in range(start_time):
            _ = video.next()

    for i,v in enumerate(video):
        if i==0:
            initial_image = v[:,:,:]

        if i<data.shape[0]:
            if reduction > 1:
                data[i,:,:] = v[:ny*reduction, :nx*reduction, 2].reshape((ny,reduction, nx, reduction)).mean(axis=1).mean(axis=2)
            else:
                data[i,:,:] = v[:,:,2]
        else:
            break
    print("finished: {:5.3} seconds\n".format(time.time()-t0))

    if plot_data:
        plt.figure(figsize=(20,10))
        if verbose: print("plotting initial image")
        t0 = time.time()
        plt.imshow(initial_image, origin="upper")
        plt.savefig("{prefix}image_{}.png".format(file_base, prefix=prefix))
        if verbose: print("finished: {:5.3} seconds\n".format(time.time()-t0))

    if verbose: print("computing fft")
    t0 = time.time()
    nt = min(max_times, min(data.shape[0], i))
    if verbose: print(data.shape[0], nt, max_times)
    window = tukey(nt,0.1)
    data = data[:nt,:,:]# * window[:,np.newaxis,np.newaxis]
    f_data = np.zeros(data.shape, dtype=np.float)

    if normalize:
        data_mean = data.mean(axis=0)
        for i in range(nt):
            data[i] -= data_mean

    # data = np.transpose(data, (1,2,0))

    # for memory efficiency purposes
    for i in range(data.shape[1]):
        f_data[:,i,:] = np.abs(np.fft.fft(data[:,i,:], axis=0)) / nt

    if verbose: print("finished: {:5.3} seconds\n".format(time.time()-t0))

    if verbose: print("finding best frequencies")
    t0 = time.time()
    freqs = np.fft.fftfreq(nt, 1/fps)
    if verbose: print(freqs.min(), freqs.max())
    bottom = np.where(freqs > min_freq)[0][0]
    top = np.where(freqs > max_freq)[0][0]
    if verbose: print(bottom, top)

    if save_results:
        import xarray as xr
        xr.DataArray(f_data[bottom:top]).to_netcdf(file_base+".nc")

    best = np.argmax(f_data[bottom:top], axis=0)
    ampl = np.max(f_data[bottom:top], axis=0) / data.mean(axis=0)

    # best = filter_best(best, ampl)

    freq_data = freqs[best+bottom]

    masked_freq_data = np.ma.array(freq_data, mask = ampl < min_amplitude)
    if verbose: print("finished: {:5.3} seconds\n".format(time.time()-t0))

    if plot_data:
        if verbose: print("plotting frequency image")
        t0 = time.time()
        plt.clf();
        plt.imshow(masked_freq_data, vmax=2, vmin=0.1, origin="upper", cmap=plt.cm.jet)
        plt.colorbar()
        plt.savefig("{prefix}fft_{}.png".format(file_base, prefix=prefix))
        if verbose: print("finished: {:5.3} seconds\n".format(time.time()-t0))

        if verbose: print("plotting amplitude image")
        t0 = time.time()
        plt.clf()
        plt.imshow(ampl, origin="upper", vmax=2, vmin=0)
        plt.colorbar()
        plt.savefig("{prefix}amplitudes_{}.png".format(file_base, prefix=prefix))
        if verbose: print("finished: {:5} seconds\n".format(time.time()-t0))

if __name__ == '__main__':
    try:
        parser= argparse.ArgumentParser(description='Compute the dominant brightness frequencies in a video. ',
                                        formatter_class=argparse.ArgumentDefaultsHelpFormatter)
        parser.add_argument('filename', action='store', default="movie.mp4", help="Name of movie file to process")
        parser.add_argument('-o',       dest="outputfile",                  action='store',  default="out_", help="prefix of output files")
        parser.add_argument('--mina',   dest="min_amplitude",   nargs="?",  action='store',  default=0.001, type=float, help="minimum amplitude to plot")
        parser.add_argument('--minf',   dest="min_freq",        nargs="?",  action='store',  default=0.2,   type=float, help="minimum frequency to examine")
        parser.add_argument('--maxf',   dest="max_freq",        nargs="?",  action='store',  default=2.5,   type=float, help="maximum frequency to examine")
        parser.add_argument('--reduce', dest="reduction",       nargs="?",  action='store',  default=2,     type=float, help="factor to reduce x/y dimensions by")
        parser.add_argument('--nt',     dest="max_times",       nargs="?",  action='store',  default=50,    type=float, help="maximum time period to examine [s]")
        parser.add_argument('--st',     dest="start_time",      nargs="?",  action='store',  default=2,     type=float, help="start time [s]")
        parser.add_argument('--plot_data',   dest='plot_data',     action='store_true',  default=False, help='verbose output')
        parser.add_argument('--save_results',dest='save_results',  action='store_true',  default=False, help='verbose output')
        parser.add_argument('--normalize',   dest='normalize',     action='store_true',  default=False, help='verbose output')
        parser.add_argument('-v', '--version', action='version', version='Process Brightness 1.0')
        parser.add_argument('--verbose',dest='verbose', action='store_true',  default=False, help='verbose output')
        args = parser.parse_args()

        verbose = args.verbose

        exit_code = main(args.filename, args.outputfile,
                          plot_data = args.plot_data,
                          save_results = args.save_results,
                          normalize = args.normalize,
                          min_amplitude = args.min_amplitude,
                          min_freq = args.min_freq,
                          max_freq = args.max_freq,
                          reduction = args.reduction,
                          max_times = args.max_times,
                          start_time = args.start_time)

        if exit_code is None:
            exit_code = 0
        sys.exit(exit_code)
    except KeyboardInterrupt as e: # Ctrl-C
        raise e
    except SystemExit as e: # sys.exit()
        raise e
    except Exception as e:
        print('ERROR, UNEXPECTED EXCEPTION')
        print(str(e))
        traceback.print_exc()
        os._exit(1)

ext: py

reason: import numpy,from scipy


text:
from typing import Dict

import numpy as np

def hypothesis_test(result, col_coef_dict: Dict[str, float]):
    """
    Perform a hypothesis test with one or multiple columns

    :param result: regression result
    :param col_coef_dict: keys are names of columns in DataFrame used for regression. Values are
        coefficients to apply to those columns.
    :return:
    """
    r_matrix = _create_r_matrix_for_hypothesis_test(result, col_coef_dict)
    return result.t_test(r_matrix, use_t=True)

def _create_r_matrix_for_hypothesis_test(result, col_coef_dict):
    """
    col_coef_dict: dict, keys are names of columns in DataFrame used for regression. Values are
                   coefficients to apply to those columns.
    """
    index_coef_dict = _replace_dict_keys_with_key_index_in_list(col_coef_dict, result.model.exog_names)
    r_matrix = _numpy_array_from_index_coef_dict(result, index_coef_dict)

    return r_matrix

def _numpy_array_from_index_coef_dict(result, index_coef_dict):
    arr = np.zeros(len(result.model.exog_names))
    _modify_array_from_index_coef_dict(arr, index_coef_dict)

    return arr

def _modify_array_from_index_coef_dict(arr, index_coef_dict):
    """
    Note: inplace
    """
    for index, coef in index_coef_dict.items():
        arr[index] = coef

def _replace_dict_keys_with_key_index_in_list(dict_, list_):
    return {list_.index(key): value for key, value in dict_.items()}

max_stars_repo_name: nickderobertis/regtools

max_issues_repo_name: nickderobertis/regtools

max_forks_repo_name: nickderobertis/regtools

reason: import numpy


text:
# coding: utf-8
'''
xgboost log model
'''

import xgboost as xgb
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split

from util import Path
from models.losses import evalmape

np.random.seed(42)

def run():
    dataset_store = pd.HDFStore(Path.dataset_path + '/dataset_0608.h5')
    train = dataset_store['train']
    train_X = train.drop(['total', 'persentage'], axis=1)
    train_y = train['persentage']
    test_X = dataset_store['test']
    test_X_base = test_X['month_avg']

    dtrain = xgb.DMatrix(data=train_X, label=train_y, feature_names=train_X.columns)
    dtest = xgb.DMatrix(data=test_X, feature_names=test_X.columns)

    param = {'booster': 'dart', #gbtree',
             'max_depth': 6,
             'min_child_weight': 2,
             'learning_rate': 0.05, # 0.03
             'subsample': 0.8,
             'colsample_bytree': 0.9,
             'objective': 'reg:linear',
             'silent': 1,
             'lambda': 0.4,
             'alpha': 0.4
             }

    evallist = [(dtrain, 'train')]
    num_round = 250
    bst = xgb.train(param, dtrain, num_round, evallist,  early_stopping_rounds=20) # feval=evalmape,

    fi_weight = bst.get_score(importance_type='weight')
    print("feature importance:")
    print(fi_weight)

    pred = bst.predict(dtest)

    #result = np.exp(pred).round().astype(int)
    result = (pred+1)*test_X_base.values
    result = result.round().astype(int)
    dataset_store.close()

    return result


text:
#
# author: Jungtaek Kim (jtkim@postech.ac.kr)
# last updated: October 8, 2021
#
"""It defines an abstract class of Bayesian optimization."""

import numpy as np
import abc
import time
import qmcpy

from bayeso import constants
from bayeso.utils import utils_bo
from bayeso.utils import utils_common
from bayeso.utils import utils_logger

class BaseBO(abc.ABC):
    """
    It is a Bayesian optimization class.

    :param range_X: a search space. Shape: (d, 2).
    :type range_X: numpy.ndarray
    :param str_surrogate: the name of surrogate model.
    :type str_surrogate: str.
    :param str_acq: the name of acquisition function.
    :type str_acq: str.
    :param str_optimizer_method_bo: the name of optimization method for
        Bayesian optimization.
    :type str_optimizer_method_bo: str.
    :param normalize_Y: flag for normalizing outputs.
    :type normalize_Y: bool.
    :param debug: flag for printing log messages.
    :type debug: bool.

    """

    def __init__(self,
        range_X: np.ndarray,
        str_surrogate: str,
        str_acq: str,
        str_optimizer_method_bo: str,
        normalize_Y: bool,
        debug: bool
    ):
        """
        Constructor method

        """

        assert isinstance(range_X, np.ndarray)
        assert isinstance(str_surrogate, str)
        assert isinstance(str_acq, str)
        assert isinstance(str_optimizer_method_bo, str)
        assert isinstance(normalize_Y, bool)
        assert isinstance(debug, bool)
        assert len(range_X.shape) == 2
        assert range_X.shape[1] == 2
        assert (range_X[:, 0] <= range_X[:, 1]).all()
        assert str_surrogate in constants.ALLOWED_SURROGATE \
            + constants.ALLOWED_SURROGATE_TREES
        assert str_acq in constants.ALLOWED_BO_ACQ
        assert str_optimizer_method_bo in constants.ALLOWED_OPTIMIZER_METHOD_BO \
            + constants.ALLOWED_OPTIMIZER_METHOD_BO_TREES

        self.range_X = range_X
        self.num_dim = range_X.shape[0]
        self.str_surrogate = str_surrogate
        self.str_acq = str_acq
        self.str_optimizer_method_bo = utils_bo.check_optimizer_method_bo(
            str_optimizer_method_bo, range_X.shape[0], debug)
        self.normalize_Y = normalize_Y
        self.logger = utils_logger.get_logger('bo_w_{}'.format(str_surrogate))
        self.debug = debug

    def _get_random_state(self, seed: constants.TYPING_UNION_INT_NONE):
        """
        It returns a random state, defined by `seed`.

        :param seed: None, or a random seed.
        :type seed: NoneType or int.

        :returns: a random state.
        :rtype: numpy.random.RandomState

        :raises: AssertionError

        """

        assert isinstance(seed, (int, constants.TYPE_NONE))

        if seed is not None:
            state_random = np.random.RandomState(seed)
        else:
            state_random = np.random.RandomState()

        return state_random

    def _get_bounds(self) -> constants.TYPING_LIST:
        """
        It returns list of range tuples, obtained from `self.range_X`.

        :returns: list of range tuples.
        :rtype: list

        """

        list_bounds = []
        for elem in self.range_X:
            list_bounds.append(tuple(elem))
        return list_bounds

    def _get_samples_grid(self, num_grids: int=constants.NUM_GRIDS_AO) -> np.ndarray:
        """
        It returns grids of `self.range_X`.

        :param num_grids: the number of grids.
        :type num_grids: int., optional

        :returns: grids of `self.range_X`. Shape: (`num_grids`:math:`^{\\text{d}}`, d).
        :rtype: numpy.ndarray

        :raises: AssertionError

        """

        assert isinstance(num_grids, int)

        initials = utils_common.get_grids(self.range_X, num_grids)
        return initials

    def _get_samples_uniform(self, num_samples: int,
        seed: constants.TYPING_UNION_INT_NONE=None
    ) -> np.ndarray:
        """
        It returns `num_samples` examples uniformly sampled.

        :param num_samples: the number of samples.
        :type num_samples: int.
        :param seed: None, or random seed.
        :type seed: NoneType or int., optional

        :returns: random examples. Shape: (`num_samples`, d).
        :rtype: numpy.ndarray

        :raises: AssertionError

        """

        assert isinstance(num_samples, int)
        assert isinstance(seed, (int, constants.TYPE_NONE))

        state_random = self._get_random_state(seed)

        list_initials = []
        for _ in range(0, num_samples):
            list_initial = []
            for elem in self.range_X:
                list_initial.append(state_random.uniform(elem[0], elem[1]))
            list_initials.append(np.array(list_initial))
        initials = np.array(list_initials)
        return initials

    def _get_samples_gaussian(self, num_samples: int,
        seed: constants.TYPING_UNION_INT_NONE=None
    ) -> np.ndarray:
        """
        It returns `num_samples` examples sampled from Gaussian distribution.

        :param num_samples: the number of samples.
        :type num_samples: int.
        :param seed: None, or random seed.
        :type seed: NoneType or int., optional

        :returns: random examples. Shape: (`num_samples`, d).
        :rtype: numpy.ndarray

        :raises: AssertionError

        """

        assert isinstance(num_samples, int)
        assert isinstance(seed, (int, constants.TYPE_NONE))

        state_random = self._get_random_state(seed)

        list_initials = []

        for _ in range(0, num_samples):
            list_initial = []
            for elem in self.range_X:
                new_mean = (elem[1] + elem[0]) / 2.0
                new_std = (elem[1] - elem[0]) / 4.0

                cur_sample = state_random.randn() * new_std + new_mean
                cur_sample = np.clip(cur_sample, elem[0], elem[1])

                list_initial.append(cur_sample)
            list_initials.append(np.array(list_initial))
        initials = np.array(list_initials)
        return initials

    def _get_samples_sobol(self, num_samples: int,
        seed: constants.TYPING_UNION_INT_NONE=None
    ) -> np.ndarray:
        """
        It returns `num_samples` examples sampled from Sobol' sequence.

        :param num_samples: the number of samples.
        :type num_samples: int.
        :param seed: None, or random seed.
        :type seed: NoneType or int., optional

        :returns: examples sampled from Sobol' sequence. Shape: (`num_samples`, d).
        :rtype: numpy.ndarray

        :raises: AssertionError

        """

        assert isinstance(num_samples, int)
        assert isinstance(seed, (int, constants.TYPE_NONE))

        sampler = qmcpy.Sobol(self.num_dim, seed=seed, graycode=True)
        samples = sampler.gen_samples(num_samples)

        samples = samples * (self.range_X[:, 1].flatten() - self.range_X[:, 0].flatten()) \
            + self.range_X[:, 0].flatten()
        return samples

    def _get_samples_halton(self, num_samples: int,
        seed: constants.TYPING_UNION_INT_NONE=None
    ) -> np.ndarray:
        """
        It returns `num_samples` examples sampled by Halton algorithm.

        :param num_samples: the number of samples.
        :type num_samples: int.
        :param seed: None, or random seed.
        :type seed: NoneType or int., optional

        :returns: examples sampled by Halton algorithm. Shape: (`num_samples`, d).
        :rtype: numpy.ndarray

        :raises: AssertionError

        """

        assert isinstance(num_samples, int)
        assert isinstance(seed, (int, constants.TYPE_NONE))

        sampler = qmcpy.Halton(self.num_dim, randomize='OWEN', generalize=False, seed=seed)
        samples = sampler.gen_samples(num_samples)

        samples = samples * (self.range_X[:, 1].flatten() - self.range_X[:, 0].flatten()) \
            + self.range_X[:, 0].flatten()
        return samples

    def get_samples(self, str_sampling_method: str,
        fun_objective: constants.TYPING_UNION_CALLABLE_NONE=None,
        num_samples: int=constants.NUM_SAMPLES_AO,
        seed: constants.TYPING_UNION_INT_NONE=None,
    ) -> np.ndarray:
        """
        It returns `num_samples` examples, sampled by a sampling method `str_sampling_method`.

        :param str_sampling_method: the name of sampling method.
        :type str_sampling_method: str.
        :param fun_objective: None, or objective function.
        :type fun_objective: NoneType or callable, optional
        :param num_samples: the number of samples.
        :type num_samples: int., optional
        :param seed: None, or random seed.
        :type seed: NoneType or int., optional

        :returns: sampled examples. Shape: (`num_samples`, d).
        :rtype: numpy.ndarray

        :raises: AssertionError

        """

        assert isinstance(str_sampling_method, str)
        assert callable(fun_objective) or fun_objective is None
        assert isinstance(num_samples, int)
        assert isinstance(seed, (int, constants.TYPE_NONE))
        assert str_sampling_method in constants.ALLOWED_SAMPLING_METHOD

        if str_sampling_method == 'grid':
            assert fun_objective is not None
            if self.debug:
                self.logger.debug('For this option, num_samples is used as num_grids.')
            samples = self._get_samples_grid(num_grids=num_samples)
            samples = utils_bo.get_best_acquisition_by_evaluation(samples, fun_objective)
        elif str_sampling_method == 'uniform':
            samples = self._get_samples_uniform(num_samples, seed=seed)
        elif str_sampling_method == 'gaussian':
            samples = self._get_samples_gaussian(num_samples, seed=seed)
        elif str_sampling_method == 'sobol':
            samples = self._get_samples_sobol(num_samples, seed=seed)
        elif str_sampling_method == 'halton':
            samples = self._get_samples_halton(num_samples, seed=seed)
        else:
            raise NotImplementedError('get_samples: allowed str_sampling_method,\
                but it is not implemented.')

        if self.debug:
            self.logger.debug('samples:\n%s', utils_logger.get_str_array(samples))

        return samples

    def get_initials(self, str_initial_method: str, num_initials: int,
        seed: constants.TYPING_UNION_INT_NONE=None,
    ) -> np.ndarray:
        """
        It returns `num_initials` examples, sampled by a sampling method `str_initial_method`.

        :param str_initial_method: the name of sampling method.
        :type str_initial_method: str.
        :param num_initials: the number of samples.
        :type num_initials: int.
        :param seed: None, or random seed.
        :type seed: NoneType or int., optional

        :returns: sampled examples. Shape: (`num_samples`, d).
        :rtype: numpy.ndarray

        :raises: AssertionError

        """

        assert isinstance(str_initial_method, str)
        assert isinstance(num_initials, int)
        assert isinstance(seed, (int, constants.TYPE_NONE))
        assert str_initial_method in constants.ALLOWED_INITIALIZING_METHOD_BO

        return self.get_samples(str_initial_method, num_samples=num_initials, seed=seed)

    @abc.abstractmethod
    def compute_posteriors(self): # pragma: no cover
        """
        It is an abstract method.

        """

        pass

    @abc.abstractmethod
    def compute_acquisitions(self): # pragma: no cover
        """
        It is an abstract method.

        """

        pass

    @abc.abstractmethod
    def optimize(self): # pragma: no cover
        """
        It is an abstract method.

        """

        pass

max_stars_repo_name: jungtaekkim/bayeso

max_issues_repo_name: POSTECH-CVLab/bayeso

max_forks_repo_name: POSTECH-CVLab/bayeso


text:
#!/usr/bin/env python3
"""
Load OWL as NetworkX (NX) graph, for NX analytics.
For each class, associate (cluster) with an ancestral group of defined
size or other criteria.
"""
###
import sys,os,json,logging,argparse,time,tqdm
import pandas as pd
import networkx as nx
from networkx.convert_matrix import from_pandas_edgelist

###
def SubclassCount(G, n):
   """Recursive all-subclass count."""
   N_sub = 0
   for nn in G.successors(n):
     N_sub += 1
     N_sub += SubclassCount(G, nn)
   return N_sub
###
def SubclassCount_InSet(G, n, nodeset):
   """Recursive all-subclass-in-set count. Slow."""
   N_sub = 0
   for nn in G.successors(n):
     in_set = bool(n in nodeset)
     if in_set:
       N_sub += 1
     N_sub += SubclassCount_InSet(G, nn, nodeset)
   return N_sub
###
def Level(G, n):
  """Level relative to ancestral root[s] at level 0"""
  roots = [n for n,d in G.in_degree() if d==0]
  min_pathlen = min([nx.shortest_path_length(G, root, n) if nx.has_path(G, root, n) else sys.maxsize for root in roots])
  return min_pathlen
###
def Graph2CYJS(G, ofile):
  cyjs = nx.cytoscape_data(G)
  fout.write(json.dumps(cyjs, indent=2))
###
def Groups2TSV(groups, fout):
  if groups is None: return
  logging.info(f"Output rows: {groups.shape[0]}; cols: {groups.shape[1]}")
  groups.to_csv(fout, '\t', index=False)
###
def GraphSummary(G):
  logging.info(nx.info(G))
  logging.info(f"nodes: {G.number_of_nodes()}; edges: {G.number_of_edges()}; directed: {G.is_directed()}")
  logging.info(f"connected: {nx.is_weakly_connected(G)}")
  logging.info(f"connected components: {nx.number_weakly_connected_components(G)}")
  logging.info(f"DAG (mono-hierarchy): {nx.is_directed_acyclic_graph(G)}")
  logging.info(f"Tree: {nx.is_tree(G)}")
  logging.info(f"Forest: {nx.is_forest(G)}")
  #
  roots = [n for n,d in G.in_degree() if d==0]
  logging.info(f"Roots: {len(roots)}")
  for root in roots:
    logging.info("root: {}: {}".format(root, nx.get_node_attributes(G, 'label')[root] if root in nx.get_node_attributes(G, 'label') else ''))
  #
  leafs = [n for n,d in G.out_degree() if d==0]
  logging.info(f"Leafs: {len(leafs)}")
  #
  singles = set(roots) & set(leafs)
  logging.info(f"Singles: {len(singles)}")
###
def GraphCount_InSet(G, nodeset):
  N_in = 0
  for n in G.nodes:
    if n in nodeset:
      N_in+=1
  return N_in
###
def Cluster(G, nodeSet, min_groupsize, max_level, setname):
  """Find ancestors with N_subclasses >= min_groupsize."""
  i_node=0; grouplist=[]; tq=None;
  roots = [n for n,d in G.in_degree() if d==0]
  logging.debug(f"Roots: {len(roots)}")
  N_inset = GraphCount_InSet(G, nodeSet)
  logging.info(f"nodeSet count: {len(nodeSet)}")
  logging.info(f"Graph nodes in nodeSet: {N_inset}")
  if N_inset==0:
    logging.error("Zero graph nodes in nodeSet.")
    return None
  for n in G.nodes:
    if tq is None: tq = tqdm.tqdm(total=len(G.nodes), unit="nodes")
    tq.update()
    i_node+=1
    level = Level(G, n)
    if level>max_level: continue
    N_sub = SubclassCount(G, n)
    label = nx.get_node_attributes(G, 'label')[n]
    in_set = bool(n in nodeSet)
    #is_root = bool(n in roots)
    N_sub_set = SubclassCount_InSet(G, n, nodeSet)
    logging.debug(f"{i_node}/{G.number_of_nodes()}. {n}: {label}; level={level}; N_sub={N_sub}; N_sub_{setname}={N_sub_set}")
    if N_sub_set >= min_groupsize:
      logging.debug(f"{i_node}/{G.number_of_nodes()}. {n}: {label}; level={level}; N_sub={N_sub}; N_sub_{setname}={N_sub_set} (> min_groupsize={min_groupsize}")
      grouplist.append((n,label,level,in_set,N_sub,N_sub_set))
  groups = pd.DataFrame({
	'Id': [Id for Id,label,level,in_set,N_sub,N_sub_set in grouplist],
	'label': [label for Id,label,level,in_set,N_sub,N_sub_set in grouplist],
	'level': [level for Id,label,level,in_set,N_sub,N_sub_set in grouplist],
	f'in_{setname}': [in_set for Id,label,level,in_set,N_sub,N_sub_set in grouplist],
	'N_sub': [N_sub for Id,label,level,in_set,N_sub,N_sub_set in grouplist],
	f'N_sub_{setname}': [N_sub_set for Id,label,level,in_set,N_sub,N_sub_set in grouplist]
	}).sort_values(by=[f'N_sub_{setname}', 'N_sub'], ascending=False)
  #print(groups.head(10)) #DEBUG
  print(groups[groups[f'in_{setname}']].head(18)) #DEBUG
  return groups

#############################################################################
if __name__=="__main__":
  parser = argparse.ArgumentParser(description="NetworkX analytics, from edgelist, optional node list, attributes.", epilog="Clustering into groups by common ancestor. Designed for EFO.")
  OPS = ['summary', 'graph2cyjs', 'cluster']
  parser.add_argument("op", choices=OPS, help='OPERATION')
  parser.add_argument("--i_edge", required=True, dest="ifile_edge", help="Input edgelist (TSV)")
  parser.add_argument("--i_node_set", dest="ifile_node_set", help="Input node set (IDs)")
  parser.add_argument("--i_node_attr", dest="ifile_node_attr", help="Input node attributes (TSV)")
  parser.add_argument("--o", dest="ofile", help="output (TSV|CYJS|etc.)")
  parser.add_argument("--setname", default="myset", help="Node set name")
  parser.add_argument("--graphname", help="Assign this name to graph.")
  parser.add_argument("--min_groupsize", type=int, default=100)
  parser.add_argument("--max_level", type=int, default=10)
  parser.add_argument("-v", "--verbose", action="count", default=0)
  args = parser.parse_args()

  logging.basicConfig(format='%(levelname)s:%(message)s', level = logging.DEBUG if args.verbose>0 else logging.INFO)

  t0 = time.time()

  logging.info(f"Reading {args.ifile_edge}")
  efo_edges = pd.read_csv(args.ifile_edge, "\t", dtype=str)
  G = from_pandas_edgelist(efo_edges, source="source", target="target", edge_attr="edge_attr", create_using=nx.DiGraph)
  G.graph['name'] = args.graphname

  fout = open(args.ofile, "w") if args.ofile else sys.stdout

  ###
  if args.ifile_node_attr:
    logging.info(f"Reading {args.ifile_node_attr}")
    efo_nodes = pd.read_csv(args.ifile_node_attr, "\t", index_col="id")
    node_attr = efo_nodes.to_dict(orient='index')
    nx.set_node_attributes(G, node_attr)
  #
  ###
  if args.ifile_node_set:
    nodeSetIds = set()
    logging.info(f"Reading {args.ifile_node_set}")
    with open(args.ifile_node_set) as fin:
      for line in fin:
        nodeSetIds.add(line.strip())
    logging.info(f"nodeSetIds in {args.setname}: {len(nodeSetIds)}")
    nx.set_node_attributes(G, {nodeId:{'in_%s'%args.setname:True} for nodeId in nodeSetIds})

  ###
  if args.op == 'summary':
    GraphSummary(G)
  #
  elif args.op == 'graph2cyjs':
    logging.info(f"Writing {args.ofile}")
    Graph2CYJS(G, fout)
  #
  ###
  elif args.op == 'cluster':
    if not args.ifile_node_set:
      parser.error(f"--i_node_set required for: {args.op}")
    GraphSummary(G)
    groups = Cluster(G, nodeSetIds, args.min_groupsize, args.max_level, args.setname)
    Groups2TSV(groups, fout)

  logging.info(('Elapsed time: %s'%(time.strftime('%Hh:%Mm:%Ss', time.gmtime(time.time()-t0)))))

max_stars_repo_name: unmtransinfo/gwas-explorer

max_issues_repo_name: unmtransinfo/gwas-explorer

max_forks_repo_name: unmtransinfo/gwas-explorer

reason: import networkx,from networkx


text:
#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt

# see weakstudy.slurm for run script and slurm.10883 for results

# the base run
basic = "./ice -ice_verif 2 -ts_type beuler -ice_tf 10.0 -ice_dtinit 10.0 -snes_monitor -snes_converged_reason -ksp_converged_reason -pc_type mg -pc_mg_levels $(( $LEV - 1 )) -da_refine $LEV"

# on mx=1536:
#   max D_SIA 1.550 m2s-1
#   dtD 7.017e-03 a   <-- time step limit for Forward Euler

dtratios = 10.0 / (0.007017 / np.array([1.0,2.0,4.0,8.0])**2)

# mx, rank, Newtons, time
raw  = [[1536,  1,  9,  146.907], #2m16.907s
        [3072,  4,  14, 339.587], #5m39.587s
        [6144,  16, 16, 630.244], #10m30.244s
        [12288, 64, 17, 960.936]] #16m0.936s
data = np.array(raw)

lev = [9,10,11,12]
L = 1800.0e3

mx = data[:,0]
m = mx * mx
h = L / mx
rank = data[:,1]
snes = data[:,2]
time = data[:,3]

def writeit(fname):
    print 'saving figure %s' % fname
    plt.savefig(fname,bbox_inches='tight')

plt.figure(1)
plt.plot(lev,time / snes,'o',ms=12)
plt.grid('on')
dof = [r'$2.4\times 10^6$', r'$9.4\times 10^6$', r'$3.8\times 10^7$', r'$1.5\times 10^8$']
xlabeldof = []
for k in range(4):
    xlabeldof += ['%d\n%s\n%d m' % (rank[k],dof[k],h[k])]
plt.xticks(lev,xlabeldof)
plt.ylabel('time (s) per Newton step')
plt.axis([8.8,12.2,10.0,70.0])
plt.yticks([20, 40, 60],['20','40','60'])
writeit('timepernewtonweak.pdf')

plt.text(lev[0]+0.1,60.0,r'the good news: $\frac{\Delta t}{\Delta t_{FE}}$',color='r',fontsize=30)
for k in range(4):
    plt.text(lev[k]-0.2,30.0,'%.1e' % dtratios[k],color='r',fontsize=12)
writeit('withdtratioweak.pdf')

#plt.show()


text:
from scipy import mat, sin, zeros
K = mat('1 0 0;0 2 0;0 0 3')
M = mat('4 1 0;1 4 1;0 1 2')/6.0
r = mat('0;-1;1')
def load(t): return r*sin(7.0*t) #
h = 0.005; duration = 6.0
# linear acceleration coefficients
A = 3.0*M;   V = 6.0*M/h
Flex  = (K + 6.0*M/(h*h)).I
MI = M.I
# initial state
t = 0
x, v, p = mat(zeros((3,1))), mat(zeros((3,1))), r*sin(omega*t)
a = MI*(p - K*x)
print "%12.9f" % t, ' '.join(["%12.9f" % x_i for x_i in x])
# iterate with linear acceleration algorithm
while (t + h/2) < duration:
    t  = t + h
    dp = r*sin(omega*t) - p
    dx = Flex*(dp + A*a + V*v)
    dv = 3.0*dx/h - 3*v - a*h/2.0
    x, v, p = x + dx, v + dv, p + dp
    a  = MI*(p - K*x)
    print "%12.9f" % t, ' '.join(["%12.9f" % x_i for x_i in x])

reason: from scipy


text:
from __future__ import annotations

from dataclasses import dataclass
from dataclasses import field
from numpy import integer
from numpy.typing import NBitBase
from pentagram.parse.word import WordComment
from pentagram.parse.word import WordIdentifier
from pentagram.parse.word import WordLine
from pentagram.parse.word import WordNumber
from pentagram.parse.word import WordTerm
from typing import Generic
from typing import TypeVar

@dataclass
class GroupTerm:
    pass

TBit = TypeVar("TBit", bound=NBitBase)

@dataclass
class GroupNumber(GroupTerm, Generic[TBit]):
    value: integer[TBit]
    value_type: type[integer[TBit]] = field(init=False)

    def __post_init__(self) -> None:
        self.value_type = type(self.value)
        assert issubclass(
            self.value_type, integer
        ), self.value_type

@dataclass
class GroupIdentifier(GroupTerm):
    name: str

@dataclass
class GroupComment(GroupTerm):
    text: str

@dataclass
class Group(GroupTerm):
    lines: list[GroupLine]

@dataclass
class GroupLine:
    terms: list[GroupTerm]

def parse_group(lines: list[WordLine]) -> Group:
    groups = [Group([])]
    substantial_indent_level = -1

    def group_end(indent_level: int) -> None:
        while indent_level < len(groups) - 1:
            groups[-2].lines[-1].terms.append(groups.pop())

    for line in lines:
        if not line.terms:
            indent_level = substantial_indent_level + 1
        elif len(line.terms) == 1 and isinstance(
            line.terms[0], WordComment
        ):
            indent_level = min(
                substantial_indent_level + 1,
                line.indent // 2,
            )
        else:
            assert line.indent % 2 == 0, line
            indent_level = line.indent // 2
            assert indent_level <= len(groups), line
            substantial_indent_level = indent_level
        if indent_level == len(groups):
            groups.append(Group([]))
        group_end(indent_level)
        groups[-1].lines.append(
            GroupLine(parse_group_terms(line.terms))
        )

    group_end(0)
    return groups[0]

def parse_group_terms(
    terms: list[WordTerm],
) -> list[GroupTerm]:
    return [parse_one_group_term(term) for term in terms]

def parse_one_group_term(term: WordTerm) -> GroupTerm:
    if isinstance(term, WordNumber):
        return GroupNumber(term.value)
    elif isinstance(term, WordIdentifier):
        return GroupIdentifier(term.name)
    elif isinstance(term, WordComment):
        return GroupComment(term.text)
    else:
        raise AssertionError(term)

max_stars_repo_name: pentagram-lang/pentagram

max_issues_repo_name: tacit-lang/tacit

max_forks_repo_name: tacit-lang/tacit

reason: from numpy


text:
# Last modified date: 12 Apr 2021
import os
import sys
import numpy as np
import pandas as pd
sys.path.insert(1, os.getcwd())

from datetime import timedelta
from tqdm.auto import tqdm
from warnings import warn
from datetime import date
from decimal import Decimal
from sklearn.linear_model import LinearRegression

from helper import Notification

class PriceAction():

    def __init__(self, price, o='open', h='high', l='low', c='close', v='volume', max_periods=None):
        """
            Derive Price Action Data from raw OHLCV Price Data.
            Parameters:
            -----------
            price (Pandas DataFrame): column: open, high, low, close, volume, index: int range
            o (str): column label for "open" price, default open
            h (str): column label for "high" price, default high
            l (str): column label for "low" price, default low
            c (str): column label for "close" price, default close
            v (str): column label for "volume" price, default volume
            max_period (int): Most recent number of price data for calculation, default to using full data
        """

        remaining_cols = list(set([o, h, l, c, v, 'date']) - set(price.columns))
        assert not remaining_cols, f"Missing columns: {remaining_cols}"

        self.price = price.copy()
        self.price = self.price.sort_values(by='date', ascending=True).reset_index(drop=True)

        data_size = self.price.shape[0]

        if max_periods is not None:
            if data_size > max_periods:
                self.price = self.price.iloc[-max_periods:]
#             else:
#                 warn(f"The data size {data_size} is less than the specified period {max_periods}. Skip slicing.")

        self.o = self.price[o]
        self.h = self.price[h]
        self.l = self.price[l]
        self.c = self.price[c]
        self.v = self.price[v]

    def get_difference(self, n_period=1, normalize=True, inplace=True):
        """
            To calculate the difference in close price between N days
        """
        close_diff = self.c.diff(n_period)
        if normalize:
            close_diff_pct = close_diff.abs() / (close_diff.abs() + self.c)

        if inplace:
            self.price[f'close_diff_{n_period}'] = close_diff
            if normalize:
                self.price[f'close_diff_pct_{n_period}'] = close_diff_pct

        return close_diff

    def get_body_size(self, inplace=True):
        """
            To calculate the size of the candlestick body (without +/- sign)
        """
        body_size = (self.c - self.o).abs()
        if inplace: self.price['body_size'] = body_size
        return body_size

    def get_shadow_to_range_ratio(self, inplace=True):
        """
            To calculate the ratio of non-body part in candletick to range of candlestick
        """
        shadow = np.where(self.c > self.o, self.h - self.c, self.c - self.l)
        uncertainty = shadow / (self.h - self.l)
        uncertainty[uncertainty == np.inf] = 0
        uncertainty = uncertainty.fillna(0)
        if inplace: self.price['uncertainty'] = uncertainty
        return uncertainty

    def get_overlapped_size(self, inplace=True):
        """
            To calculate the percentage of intersected price range between candlesticks
        """
        def overlapping(x):
            a = np.arange(x[0], x[1], 0.005)
            b = np.arange(x[2], x[3], 0.005)
            o = len(set(b) & set(a))
            if len(a) == 0: return 0
            return o / len(a)

        comb = pd.concat([self.c, self.o], axis=1)
        top = comb.max(axis=1)
        bottom = comb.min(axis=1)
        overlapped = (pd.concat([bottom.shift(1), top.shift(1), bottom, top], axis=1)
                     .fillna(0).apply(overlapping, axis=1))

        if inplace: self.price['overlapped'] = overlapped
        return overlapped

    def get_volume_total_ratio(self, n_period=10, inplace=True):
        """
            To calculate the ratio of specified date volume to the total volume of last N period
        """
        volume_sum = self.v.rolling(n_period).sum()
        volume_pct = self.v / volume_sum
        volume_pct[volume_pct == np.inf] = 0
        if inplace: self.price['volume_total_ratio'] = volume_pct
        return volume_pct

    def get_close_to_high_ratio(self, n_period=10, inplace=True):
        """
            To calculate the ratio of closing price to previous high in last N period
        """
        close_to_high_ratio = self.c / self.h.rolling(n_period).max().fillna(self.c)
        if inplace: self.price['close_to_high_ratio'] = close_to_high_ratio
        return close_to_high_ratio

    def get_accumulated_gain(self, n_period=10, inplace=True):
        """
            To calculate the accumulated gain/loss of price between open and close for last N period
        """
        accumulated_gain = (self.c - self.o).rolling(n_period).sum().fillna(0)
        self.price['accumulated_gain'] = accumulated_gain
        return accumulated_gain

    def get_volume_distribution(self, n_period=10, inplace=True):
        """
            To calculate the skewness of volume for last N period
        """
        volume_distribution = self.v.rolling(n_period).skew().fillna(0)
        if inplace: self.price['volume_distribution'] = volume_distribution
        return volume_distribution

    def get_target(self, n_period=10, inplace=True):
        """
            To calculate the percentage of gain over next N period
        """
        gained_count = ((self.c - self.o) > 0).rolling(n_period).mean()
        if inplace: self.price['target'] = gained_count
        return gained_count

    def get_range(self, n_period=10, inplace=True):
        """
            To calculate the position of latest price based on last n period range
        """
        high = self.c.iloc[-n_period:].max()
        low = self.c.iloc[-n_period:].min()
        current = self.c.iloc[-1]
        if abs(high - low) < 0.005:
            n_range = None
        else:
            n_range = (current - low) / (high - low)

        if inplace: self.price[f'{n_period} days range'] = n_range
        return n_range

def get_feature(price, max_periods=None):
    if price.shape[0] <= 1: return price
    pa = PriceAction(price, max_periods=max_periods)

    pa.get_accumulated_gain(5)
    pa.get_body_size()
    try:
        pa.get_overlapped_size()
    except Exception as e:
        print(e)
    pa.get_shadow_to_range_ratio()

    pa.get_close_to_high_ratio(5)
    pa.get_difference(1)
    pa.get_difference(5)
    pa.get_difference(20)

    pa.get_volume_distribution(5)
    pa.get_volume_total_ratio(5)
    pa.get_target(5)

    pa.get_range(50)
    pa.get_range(100)
    pa.get_range(200)
    return pa.price

FOLDER = "C:/Users/ycwong/Desktop/Bursa Quickview/data"
BURSA_METADATA_FILE = 'Bursa Metadata.xlsx'
PRICE_DATA_FILE = "prices.pkl"
today = date.today().strftime("%Y-%m-%d")
OUTPUT_FOLDER = 'C:/Users/ycwong/Desktop/Bursa Quickview/output/'
OUTPUT_FILE = f"result ({today}).xlsx"
MAX_LAG = 200

bursa = pd.read_excel(os.path.join(FOLDER, BURSA_METADATA_FILE))
prices = pd.read_pickle(os.path.join(FOLDER, PRICE_DATA_FILE))
feature_list = [pd.Series(dtype=float)] * len(prices.keys())

for i, code in enumerate(tqdm(list(prices.keys()))):
    if len(code) != 4: continue
    price = prices[code]
    if price is None:
        print(f"No data for {code}, skip process.")
        continue
    feature = get_feature(price, MAX_LAG)
    feature_list[i] = feature.iloc[-1]

features = pd.concat(feature_list, axis=1)
features = features.T.reset_index(drop=True)
features.sort_values(by=['close_to_high_ratio', 'volume_total_ratio'], inplace=True)

meta_cols = ["bursacode", "name", "alias", "description", "economicsectorcode", "industrygroupcode"]
bursa_meta = bursa[meta_cols].dropna(subset=['bursacode']).rename(columns={"bursacode": "code"})
results = bursa_meta.merge(features, on='code', how='right')
results['date'] = results['date'].dt.date
results['close'] = results['close'].astype(float)
results = results[results['code'].str.len() <= 4]
results = results.drop(['open', 'high', 'low', 'volume'], axis=1)
results = results.rename(columns={'economicsectorcode': 'Sector', 'industrygroupcode': 'Industry'})
results = results[results['date'] >= results['date'].max() - timedelta(days=3)]
results = results.dropna(subset=['name'])
results.sort_values(by=['Sector', 'Industry', '50 days range', '100 days range', '200 days range'])
os.chdir(OUTPUT_FOLDER)
results.to_excel(OUTPUT_FILE, index=False)

tables = []
noti = Notification(TO="ycfkjc@Hotmail.com")

for lag in [1, 5, 20]:
    col = f'close_diff_{lag}'
    pct_col = f'close_diff_pct_{lag}'
    table = results[['code', 'alias', 'Sector', 'Industry', 'close', col, pct_col]].copy()
#     table = table.sort_values(by=[pct_col, col], ascending=False)
    table['Industry'] = table['Industry'].fillna('-')
    table = table[table['close'] > 0.2]
    table[col] = table[col].apply(noti.rounding)
    table[pct_col] = table[pct_col].apply(noti.rounding)
    table = table.dropna(axis=0, subset=[col, pct_col])
    top_gain = noti.create_html(table.sort_values(by=[pct_col, col], ascending=False).head(10), gain=True)
    top_loss = noti.create_html(table.sort_values(by=[pct_col, col], ascending=True).head(10), gain=False)
    tables.append(top_gain)
    tables.append(top_loss)

HTML = noti.get_style() + '

'.join(tables)
noti.send_email(f"Daily KLSE Screener ({today})", HTML=HTML, FILENAMES=[OUTPUT_FILE])


text:
__author__ = 'julia'
# !/usr/bin/python

try:
    from mayavi import mlab
except ImportError:
    print 'mayavi not installed, visualization methods will not be available'
from gurobipy import *
from gurobipy import Model
import numpy as np
from skeleton import networkx_utils
from skeleton import knossos_utils
import utils

import math

class ILPParameters(object):
    # As proposed in python Cookbook. (usage a bit like a c struct)
    def __init__(self, dummy_edge_cost=100, distance_cost=1, comb_angle_cost=5, angle_cost_factor=0,
                 selection_cost=-50, pairwise=True, exclusive_distance=False, exclusive_distance_threshold=25, **kwds):
        # First set default values
        default_dic = {}
        default_dic['dummy_edge_cost'] = dummy_edge_cost
        default_dic['distance_cost'] = distance_cost
        default_dic['comb_angle_cost'] = comb_angle_cost
        default_dic['angle_cost_factor'] = angle_cost_factor
        default_dic['selection_cost'] = selection_cost
        default_dic['pairwise'] = pairwise
        default_dic['exclusive_distance'] = exclusive_distance
        default_dic['exclusive_distance_threshold'] = exclusive_distance_threshold

        self.__dict__.update(default_dic)
        self.__dict__.update(kwds)

def get_angle(g, i, cur_dir_vec, noise_sigma=False):
    # Get ground truth direction vectors
    # #Get ground truth direction vectors
    # pos = g.nx_graph.node[node_id]['position']
    # if tuple(pos) in all_skeletons.coords_to_id:
    #     node_id = all_skeletons.coords_to_id[tuple(pos)]
    # # if 'dir_vector' in all_skeletons.nx_graph.node[node_id]:
    #     dir_vec_of_node = all_skeletons.nx_graph.node[node_id]['dir_vector']
    # else:
    #     dir_vec_of_node = np.random.random((1, 3))
    if not 'dir_vector' in g.nx_graph.node[i]:
        dir_vec_of_node = [0, 0, 1]
    else:
        dir_vec_of_node = g.nx_graph.node[i]['dir_vector']

    # dir_vec_of_node = np.array([0, 0, 1])
    angle = utils.angle_between_vecs(dir_vec_of_node, cur_dir_vec)
    if noise_sigma:
        mu = 0
        sigma = noise_sigma
        angle_noise = np.abs(sigma * np.random.randn() + mu)
        # angle_noise = np.random.uniform(0, np.pi/2.)
        angle += angle_noise

    if angle > np.pi:
        angle = angle - np.pi

    if angle > np.pi / 2.:
        angle = np.abs(angle - np.pi)

    if math.isnan(angle):
        print "angle is not defined, check why"
        angle = 0
    return angle

def create_graph_from_solution(g_old, solution):
    print " creating new graph"
    nx_skeleton = networkx_utils.NxSkeleton()
    position_dic = g_old.node_feature_to_dic('position')
    nx_skeleton.initialize_from_edgelist(position_dic, solution, only_from_edgelist=True)
    return nx_skeleton

def check_consistency(comb_solution, solution):
    print "checking consistency"
    for comb_edge, sol in comb_solution.iteritems():
        node_id = comb_edge[0]
        neighbour1 = comb_edge[1]
        neighbour2 = comb_edge[2]
        if sol == 1.0:
            # Check for those combi edges that are on, that also the corresponding edge variables are on
            assert solution[node_id, neighbour1] == 1
            assert solution[node_id, neighbour2] == 1
        else:
            assert solution[node_id, neighbour1] + solution[node_id, neighbour2] <= 1

def add_pairwaisefactors_to_model(model, g, params, variables):
    # Add additional variables that represent the combination of two edges respectively
    edge_combinations = {}
    for node_id, node_attr in g.nx_graph.nodes_iter(data=True):
        neighbours = g.nx_graph.neighbors(node_id)
        neighbours.sort()
        for ii in range(len(neighbours) - 1):
            neighbour_id1 = neighbours[ii]
            for jj in range(len(neighbours) - ii - 1):
                neighbour_id2 = neighbours[jj + ii + 1]
                spanning_angle = node_attr['spanning_angles'][(neighbour_id1, neighbour_id2)]
                edge_combinations[node_id, neighbour_id1, neighbour_id2] = \
                    model.addVar(obj=(np.abs(spanning_angle - np.pi) * params.comb_angle_cost) ** 2, vtype=GRB.BINARY,
                                 name='node_%i_edge_%i_%i' % (node_id, neighbour_id1, neighbour_id2))

    model.update()
    # Add corresponding constraints to new introduced variables (specifically,
    # combination of egdes should only be switched on if both corresponding edges are switched on)

    for node_id, node_attr in g.nx_graph.nodes_iter(data=True):
        neighbours = g.nx_graph.neighbors(node_id)
        neighbours.sort()
        for ii in range(len(neighbours) - 1):
            neighbour_id1 = neighbours[ii]
            for jj in range(len(neighbours) - ii - 1):
                neighbour_id2 = neighbours[jj + ii + 1]
                model.addConstr(
                    edge_combinations[node_id, neighbour_id1, neighbour_id2] <= variables[node_id, neighbour_id1])
                model.addConstr(
                    edge_combinations[node_id, neighbour_id1, neighbour_id2] <= variables[node_id, neighbour_id2])
                model.addConstr(edge_combinations[node_id, neighbour_id1, neighbour_id2] + 1 >=
                                variables[node_id, neighbour_id2] + variables[node_id, neighbour_id1])
    model.update()
    return model, edge_combinations

def add_proximity_exclusivity(model, g, selection_cost_variables, distance_threshold):
    for node_id, node_attr in g.nx_graph.nodes_iter(data=True):
        neighbours = g.nx_graph.neighbors(node_id)
        neighbours.sort()
        neighbours_in_proximity = []
        for neighbour in neighbours:
            distance = g.nx_graph[node_id][neighbour]['distance']
            if distance < distance_threshold:
                neighbours_in_proximity.append(neighbour)
        if neighbours_in_proximity:
            model.addConstr(selection_cost_variables[node_id] +
                            1. / len(neighbours_in_proximity) * quicksum(
                                selection_cost_variables[neighbour] for neighbour in neighbours_in_proximity) <= 1)
    model.update()
    return model

def calculate_ILP_linear(g, params,
                         verbose=True, visualize=False, timelimit=120):
    # This is a reimplementaiton of the original version of setting up the ILP for microtubule tracking. (as described
    # in paper)
    if verbose:
        print "starting constructing ILP"
        for key, value in params.__dict__.iteritems():
            print key, value
    model = Model()
    if timelimit > 0:
        model.setParam(GRB.Param.TimeLimit, timelimit)
    variables = {}
    # Add variables to the model, per edge in the graph one variable
    for edge in g.nx_graph.edges_iter(data=True):
        i = edge[0]
        j = edge[1]
        dist = edge[2]['distance']
        cur_dir_vector = edge[2]['dir_vector']
        angle_cost = 0
        if params.angle_cost_factor != 0:
            angle1 = get_angle(g, i, cur_dir_vector)
            angle2 = get_angle(g, j, cur_dir_vector)
            angle_cost = (angle1 + angle2) * params.angle_cost_factor

        variables[i, j] = model.addVar(obj=(dist * params.distance_cost + angle_cost) ** 1, vtype=GRB.BINARY,
                                       name='e' + str(i) + '_' + str(j))
        variables[j, i] = variables[i, j]

    sink_variables = {}
    selection_cost_variables = {}
    for node_id in g.nx_graph.nodes_iter():
        # Add dummy node (sink), here the costs decide about how expensive it is to open / close a chain
        sink_variables[node_id] = model.addVar(obj=params.dummy_edge_cost, vtype=GRB.BINARY,
                                               name='dummy_%s' % str(node_id))
        # Add cost for selecting a candidate (this is uniformly set,
        # but can potentially exchanged with adaptive costs from underlying pipeline)
        selection_cost_variables[node_id] = model.addVar(obj=params.selection_cost, vtype=GRB.BINARY,
                                                         name='selectioncost_%s' % str(node_id))
    model.update()

    # Add consistency constraints for guaranteeing coherent topology
    # If a candidate is selected, two edge variables should also be selected
    for node_id in g.nx_graph.nodes_iter():
        neighbours = g.nx_graph.neighbors(node_id)
        model.addConstr(2 * selection_cost_variables[node_id] - (quicksum(variables[node_id, j] for j in neighbours) +
                                                                 sink_variables[node_id]) <= 0,
                        name='nodeconstraint1_%i' % node_id)
        model.addConstr(2 * selection_cost_variables[node_id] - (quicksum(variables[node_id, j] for j in neighbours) +
                                                                 sink_variables[node_id]) >= 0,
                        name='nodeconstraint2_%i' % node_id)

    # model.addConstr(selection_cost_variables[0] == 1,
    #                 name='nodeconstraint3_%i' %node_id)
    model.update()
    # If a link is selected, the two incident candidates on this link also need to be selected
    # The other way round is not true, if a candidate is selected
    for edge in g.nx_graph.edges_iter():
        node_id1 = edge[0]
        node_id2 = edge[1]
        model.addConstr(
            2 * variables[edge] - selection_cost_variables[node_id1] - selection_cost_variables[node_id2] <= 0,
            name='edgeconstraint_%i_%i' % edge)

    model.update()
    print "Intermediate number of variables before adding combination aspect", model.NumVars

    if params.pairwise:
        # Microtubules are rigid. Add cost for bending
        model, edge_combinations = add_pairwaisefactors_to_model(model, g, params, variables)
    if params.exclusive_distance:
        # Too many candidates for the same microtubules cause parallel occurence
        model = add_proximity_exclusivity(model, g, selection_cost_variables, params.exclusive_distance_threshold)
    if verbose:
        print "starting Optimization"
        print "Number of variables %i" % model.NumVars
        print "Number of linear constraints %i" % model.NumConstrs
    # model.display()
    model.optimize()

    print 'Obj:', model.ObjVal
    # Contains a label for each edge
    solution = model.getAttr("X", variables)
    if params.pairwise:
        comb_solution = model.getAttr('X', edge_combinations)
        check_consistency(comb_solution, solution)
    new_edgelist = []
    for edge, sol in solution.iteritems():
        if sol == 1:
            new_edgelist.append(edge)

    if visualize:
        g.visualize(edgelist=new_edgelist, only_edges_from_edgelist=True)
        mlab.show()

    # Create new graph from solution
    nx_skeleton_solution = create_graph_from_solution(g, new_edgelist)
    nx_skeleton_solution.print_statistics()

    return nx_skeleton_solution, model

def calculate_ILP(g, quadratic=False, dummy_edge_cost=100, distance_cost=1, comb_angle_cost=5,
                  verbose=True, visualize=False, linear=False, angle_cost_factor=0):
    print "starting constructing ILP"
    model = Model()
    # model.setParam(GRB.Param.TimeLimit, 7200.0)
    variables = {}
    # Add variables to the model, per edge in the graph one variable
    for edge in g.nx_graph.edges_iter(data=True):
        i = edge[0]
        j = edge[1]
        dist = edge[2]['distance']

        cur_dir_vector = edge[2]['dir_vector']
        angle_cost = 0
        if angle_cost_factor != 0:
            angle1 = get_angle(g, i, cur_dir_vector)
            angle2 = get_angle(g, j, cur_dir_vector)
            angle_cost = (angle1 + angle2) * angle_cost_factor
            # angle = utils.angle_between_vecs(dir_vec_of_node, cur_dir_vec)
        # print dist
        # print angle1, angle2, cur_dir_vector
        variables[i, j] = model.addVar(obj=(dist * distance_cost + angle_cost) ** 1, vtype=GRB.BINARY,
                                       name='e' + str(i) + '_' + str(j))
        variables[j, i] = variables[i, j]
    # This is the number for the dummy node. This guarantees that the id is not already occcupied by a node id
    n = np.max(g.nx_graph.nodes()) + 1

    g.add_geom_features_to_edges()
    g.add_geom_features_to_nodes()
    for node_id in g.nx_graph.nodes_iter():
        # Add dummy node, here the costs decide about how expensive it is to open / close a chain
        variables[node_id, n] = model.addVar(obj=dummy_edge_cost, vtype=GRB.BINARY, name='dummy1_%s' % str(node_id))
        variables[node_id, n + 1] = model.addVar(obj=dummy_edge_cost, vtype=GRB.BINARY, name='dummy2_%s' % str(node_id))
        variables[n, node_id] = variables[node_id, n]
        variables[n + 1, node_id] = variables[node_id, n + 1]
    model.update()

    for node_id in g.nx_graph.nodes_iter():
        neighbours = g.nx_graph.neighbors(node_id)
        neighbours.extend([n, n + 1])
        model.addConstr(quicksum(variables[node_id, j] for j in neighbours) == 2, name='node_%i' % node_id)

    model.update()
    print "Intermediate number of variables before adding combination aspect", model.NumVars
    if quadratic:
        if verbose:
            print "adding quadratic terms"
        quad = QuadExpr()
        for node_id, node_attr in g.nx_graph.nodes_iter(data=True):
            neighbours = g.nx_graph.neighbors(node_id)
            neighbours.sort()

            # Add the quadratic terms for the combination of each edge pair incident to the same node

            for ii in range(len(neighbours) - 1):
                neighbour_id1 = neighbours[ii]
                for jj in range(len(neighbours) - ii - 1):
                    neighbour_id2 = neighbours[jj + ii + 1]
                    spanning_angle = node_attr['spanning_angles'][(neighbour_id1, neighbour_id2)]
                    quad.add(variables[node_id, neighbour_id1] * variables[node_id, neighbour_id2] * (
                        np.abs(spanning_angle - np.pi) * comb_angle_cost) ** 2)

        # Add linear costs to quadratic expression
        for edge in g.nx_graph.edges_iter(data=True):
            i = edge[0]
            j = edge[1]
            dist = edge[2]['distance']
            cur_dir_vector = edge[2]['dir_vector']
            angle1 = get_angle(g, i, cur_dir_vector)
            angle2 = get_angle(g, j, cur_dir_vector)
            angle_cost = (angle1 + angle2) * angle_cost_factor
            total_cost = angle_cost + dist * distance_cost
            quad.add(variables[i, j] * (total_cost))

        for node_id in g.nx_graph.nodes_iter():
            # Add dummy node, here the costs decide about how expensive it is to open / close a chain
            quad.add(variables[node_id, n] * dummy_edge_cost)
            quad.add(variables[node_id, n + 1] * dummy_edge_cost)
        print "model collection costs finished"
        model.setObjective(quad)
        model.update()
        print "Number of variables after adding quadratic component", model.NumVars

    elif linear:
        # Add additional variables that represent the combination of two edges respectively
        edge_combinations = {}
        for node_id, node_attr in g.nx_graph.nodes_iter(data=True):
            neighbours = g.nx_graph.neighbors(node_id)
            neighbours.sort()
            for ii in range(len(neighbours) - 1):
                neighbour_id1 = neighbours[ii]
                for jj in range(len(neighbours) - ii - 1):
                    neighbour_id2 = neighbours[jj + ii + 1]
                    spanning_angle = node_attr['spanning_angles'][(neighbour_id1, neighbour_id2)]
                    edge_combinations[node_id, neighbour_id1, neighbour_id2] = \
                        model.addVar(obj=(np.abs(spanning_angle - np.pi) * comb_angle_cost) ** 2, vtype=GRB.BINARY,
                                     name='node_%i_edge_%i_%i' % (node_id, neighbour_id1, neighbour_id2))

        model.update()
        # Add corresponding constraints to new introduced variables (specifically,
        # combination of egdes should only be switched on if both corresponding edges are switched on)

        for node_id, node_attr in g.nx_graph.nodes_iter(data=True):
            neighbours = g.nx_graph.neighbors(node_id)
            neighbours.sort()
            for ii in range(len(neighbours) - 1):
                neighbour_id1 = neighbours[ii]
                for jj in range(len(neighbours) - ii - 1):
                    neighbour_id2 = neighbours[jj + ii + 1]
                    model.addConstr(
                        edge_combinations[node_id, neighbour_id1, neighbour_id2] <= variables[node_id, neighbour_id1])
                    model.addConstr(
                        edge_combinations[node_id, neighbour_id1, neighbour_id2] <= variables[node_id, neighbour_id2])
                    model.addConstr(edge_combinations[node_id, neighbour_id1, neighbour_id2] + 1 >=
                                    variables[node_id, neighbour_id2] + variables[node_id, neighbour_id1])
        model.update()

    if verbose:
        print "starting Optimization"
        print "Number of variables %i" % model.NumVars
        print "Number of linear constraints %i" % model.NumConstrs
    model.optimize()
    # Contains a label for each edge
    solution = model.getAttr("X", variables)
    if linear:
        comb_solution = model.getAttr("X", edge_combinations)

    new_edgelist = []

    count2 = 0
    for edge, sol in solution.iteritems():
        if n + 1 in edge or n in edge:
            if sol == 1:
                count2 += 1
        elif sol == 1:
            if not edge in new_edgelist or not (edge[1], edge[0]) in new_edgelist:
                new_edgelist.append(edge)

    if visualize:
        g.visualize(edgelist=new_edgelist, only_edges_from_edgelist=True)
        mlab.show()

    # Create new graph from solution
    graph_solution = create_graph_from_solution(g, new_edgelist)
    graph_solution.print_statistics()
    if linear:
        check_consistency(comb_solution, solution)

    if linear:
        return graph_solution, [model, solution, comb_solution]
    else:
        return graph_solution, [model, solution]

def validate_ILP(ground_truth_filename, test_filename, validation_parameter,
                 results_outputfile, prior, scaling=[4.6, 4.6, 50], gt_additive=[0, 0, 0]):
    print 'starting'
    scaling = np.array(scaling)
    gt_nx_skeleton = networkx_utils.from_filename_to_nx_skeleton(ground_truth_filename, scaling)
    gt_nx_skeleton.scale_positions([1, 1, 1], additive=gt_additive)
    test_nx_skeleton = networkx_utils.from_filename_to_nx_skeleton(test_filename, scaling)
    print test_nx_skeleton.nx_graph.number_of_nodes(), 'number of nodes'
    if test_nx_skeleton.nx_graph.number_of_nodes() < 3:
        print " no nodes in skeleton"
        return False
    # Remove nodes from skeleton that are above a specific distance (mask)
    masked_skeleton = networkx_utils.filter_nodes_based_on_distance(gt_nx_skeleton, test_nx_skeleton,
                                                                    max_radius=validation_parameter['masked_radius'])

    test_skeleton = networkx_utils.match_components_of_two_nx_skeletons_hungarian(gt_nx_skeleton, masked_skeleton,
                                                                                max_radius=validation_parameter[
                                                                                'allowed_max_distance'])
    true_positive_skeleton, relabeled_test_skeleton, false_positive_nx_skeleton = \
        networkx_utils.align_two_skeletons_based_on_partner(gt_nx_skeleton, test_skeleton)
    print 'starting validation'
    test_edgelist = relabeled_test_skeleton.nx_graph.edges()
    ground_truth_edgelist = gt_nx_skeleton.nx_graph.edges()
    edge_errors = utils.from_lists_to_precision_recall(ground_truth_edgelist, test_edgelist)
    if not edge_errors:
        return False
    for key, value in edge_errors.iteritems():
        utils.write_data_to_h5(results_outputfile, value, key + "_edgebased", overwrite=True, verbose=True)
        utils.write_data_to_h5(results_outputfile, test_filename, 'validation_filename', overwrite=True)
        utils.write_data_to_h5(results_outputfile, prior, 'prior', overwrite=True)
        utils.print_dic(edge_errors)
    masked_skeleton.scale_positions(scale_factor=1. / scaling)
    outputfilename = test_filename.replace('.nml', '_masked_skeleton.nml')
    knossos_utils.from_nx_graphs_to_knossos([masked_skeleton.nx_graph], outputfilename)

    true_positive_skeleton.scale_positions(scale_factor=1. / scaling)
    outputfilename = test_filename.replace('.nml', '_true_positive_skeleton.nml')
    knossos_utils.from_nx_graphs_to_knossos([true_positive_skeleton.nx_graph], outputfilename)

    false_positive_nx_skeleton.scale_positions(scale_factor=1. / scaling)
    outputfilename = test_filename.replace('.nml', '_false_positive_skeleton.nml')
    knossos_utils.from_nx_graphs_to_knossos([false_positive_nx_skeleton.nx_graph], outputfilename)
    return True

max_stars_repo_name: juliabuhmann/microtubules

max_issues_repo_name: juliabuhmann/microtubules

max_forks_repo_name: juliabuhmann/microtubules

reason: import numpy,import networkx