text:
"""
Created by Constantin Philippenko, 10th January 2022.
"""

import copy
import math
from abc import abstractmethod, ABC

from numpy.linalg import inv
from scipy.linalg import sqrtm
from tqdm import tqdm

from src.CompressionModel import CompressionModel
from src.JITProduct import *
from src.PickleHandler import pickle_saver
from src.SyntheticDataset import MAX_SIZE_DATASET
from src.Utilities import print_mem_usage

ONLY_ADDITIVE_NOISE = False
USE_MOMENTUM = False
BETA = 0
REGULARIZATION = 0

DISABLE = False
CORRECTION_SQUARE_COV = False
CORRECTION_DIAG = False

def log_sampling_xaxix(size_dataset):
    log_len = np.int(math.log10(size_dataset))
    residual_len = math.log10(size_dataset) - log_len
    log_xaxis = [[math.pow(10, a) * math.pow(10, i / 100) for i in range(100)] for a in range(log_len)]
    log_xaxis.append([math.pow(10, log_len) * math.pow(10, i / 100) for i in range(int(100 * residual_len))])
    log_xaxis = np.concatenate(log_xaxis, axis=None)
    log_xaxis = np.unique(log_xaxis.astype(int))
    return log_xaxis

class SeriesOfSGD:

    def __init__(self, list_of_sgd) -> None:
        super().__init__()
        self.series = []
        for serie in list_of_sgd:
            assert isinstance(serie, SGDRun), "The object added to the series is not of type SGDRun."
            self.series.append(serie)

    def append(self, list_of_sgd):
        for serie in list_of_sgd:
            assert isinstance(serie, SGDRun), "The object added to the series is not of type SGDRun."
            self.series.append(serie)

    def save(self, filename: str):
        pickle_saver(self, filename)

class SGDRun:

    def __init__(self, size_dataset, all_avg_w, last_w, losses, avg_losses, diag_cov_gradients, label=None) -> None:
        super().__init__()
        self.size_dataset = size_dataset
        self.all_avg_w = all_avg_w
        self.last_w = last_w
        self.losses = np.array(losses)
        self.avg_losses = np.array(avg_losses)
        self.log_xaxis = log_sampling_xaxix(size_dataset)
        self.diag_cov_gradients = diag_cov_gradients
        self.label = label

class SGD(ABC):
    NB_EPOCH = 1

    def __init__(self, synthetic_dataset, reg: int = REGULARIZATION) -> None:
        super().__init__()
        self.do_logistic_regression = synthetic_dataset.do_logistic_regression
        self.synthetic_dataset = synthetic_dataset
        self.w_star = self.synthetic_dataset.w_star
        self.SIZE_DATASET, self.DIM = self.synthetic_dataset.size_dataset, self.synthetic_dataset.dim
        self.additive_stochastic_gradient = ONLY_ADDITIVE_NOISE
        self.root_square_upper_sigma = sqrtm(self.synthetic_dataset.upper_sigma)
        self.inv_root_square_upper_sigma = inv(self.root_square_upper_sigma)

        self.Q, self.D = np.identity(self.DIM), np.identity(self.DIM)
        self.approx_hessian = np.identity(self.DIM)
        self.debiased_hessian = np.identity(self.DIM)
        self.reg = reg
        self.compressor = None

    def compute_empirical_risk(self, w, data, labels):
        if CORRECTION_SQUARE_COV:
            data = data @ self.inv_root_square_upper_sigma.T
        if CORRECTION_DIAG:
            data = data @ self.Q
        return 0.5 * np.linalg.norm(data @ w - labels) ** 2 / len(labels)

    def compute_true_risk(self, w, data, labels):
        if data is None:
            return 0
        if CORRECTION_SQUARE_COV:
            w_star = self.inv_root_square_upper_sigma @ self.w_star
            return constant_product(0.5, vectorial_norm(minus(w, w_star)))
        elif CORRECTION_DIAG:
            w_star = matrix_vector_product(self.Q.T, self.w_star)
            return wAw_product(0.5, minus(w, w_star), self.D)
        return wAw_product(0.5, minus(w, self.w_star), self.synthetic_dataset.upper_sigma)

    def compute_stochastic_gradient(self, w, data, labels, index):
        x, y = data[index], labels[index]
        if CORRECTION_SQUARE_COV:
            x = matrix_vector_product(self.inv_root_square_upper_sigma, x)
        elif CORRECTION_DIAG:
            x = matrix_vector_product(self.Q.T, x)
        return (x @ w - y) * x

    def compute_full_gradient(self, w, data, labels):
        return (data @ w - labels) @ data / len(labels)

    def compute_additive_stochastic_gradient(self, w, data, labels, index):
        x, y = data[index], labels[index]
        if self.additive_stochastic_gradient and self.do_logistic_regression:
            raise ValueError("Compute only additive stochastic gradient is not possible in the logistic setting")
        if CORRECTION_SQUARE_COV:
            x = matrix_vector_product(self.inv_root_square_upper_sigma, x)
        elif CORRECTION_DIAG:
            x = matrix_vector_product(self.Q.T, x)
        return self.D.dot(w) - y * x

    def sgd_update(self, w, gradient, gamma):
        return w - gamma * gradient - self.reg * (w - self.synthetic_dataset.w0)

    def get_step_size(self, it: int, gamma: int, deacreasing_step_size: bool = False):
        if deacreasing_step_size:
            return 1 / (np.sqrt(it) * max(self.synthetic_dataset.eigenvalues))
        return gamma

    def gradient_descent(self, label: str = None, deacreasing_step_size: bool = False) -> SGDRun:
        log_xaxis = log_sampling_xaxix(self.synthetic_dataset.size_dataset)
        current_w = self.synthetic_dataset.w0
        all_avg_w = [current_w]
        avg_w = copy.deepcopy(current_w)
        it = 1
        losses = [self.compute_true_risk(current_w, self.synthetic_dataset.X_complete, self.synthetic_dataset.Y)]
        avg_losses = [losses[-1]]
        for epoch in range(self.NB_EPOCH):
            indices = np.arange(self.SIZE_DATASET)
            for idx in tqdm(indices, disable=DISABLE):
                if idx % MAX_SIZE_DATASET == 0 and idx != 0:
                    print("Regenerating ...")
                    self.synthetic_dataset.regenerate_dataset()
                gamma = self.get_step_size(it, self.synthetic_dataset.gamma, deacreasing_step_size)
                it += 1

                if CORRECTION_DIAG:
                    self.Q, self.D = self.synthetic_dataset.Q, self.synthetic_dataset.D #diagonalization(self.debiased_hessian) #

                if self.additive_stochastic_gradient:
                    grad = self.compute_additive_stochastic_gradient(current_w, self.synthetic_dataset.X_complete,
                                                                     self.synthetic_dataset.Y, idx % MAX_SIZE_DATASET)
                else:
                    grad = self.compute_stochastic_gradient(current_w, self.synthetic_dataset.X_complete,
                                                            self.synthetic_dataset.Y, idx % MAX_SIZE_DATASET)
                grad = self.gradient_processing(grad)

                if idx == 0:
                    self.approx_hessian = np.kron(grad, grad).reshape((self.DIM, self.DIM))
                else:
                    self.approx_hessian = np.kron(grad, grad).reshape((self.DIM, self.DIM)) / it + self.approx_hessian * (it - 1)/ it

                current_w = self.sgd_update(current_w, grad, gamma)
                avg_w = current_w / it + avg_w * (it - 1) / it
                all_avg_w.append(avg_w)
                if idx in log_xaxis[1:]:
                    losses.append(self.compute_true_risk(current_w, self.synthetic_dataset.X_complete, self.synthetic_dataset.Y))
                    avg_losses.append(self.compute_true_risk(avg_w, self.synthetic_dataset.X_complete, self.synthetic_dataset.Y))

        print_mem_usage("End of sgd descent ...")

        if self.synthetic_dataset.use_ortho_matrix:
            cov_matrix = self.synthetic_dataset.ortho_matrix.T.dot(self.approx_hessian).dot(self.synthetic_dataset.ortho_matrix)
        else:
            cov_matrix = self.approx_hessian

        return SGDRun(self.synthetic_dataset.size_dataset, all_avg_w, current_w, losses, avg_losses, np.diag(cov_matrix),
                      label=label)
        # return SGDRun(current_w, losses, avg_losses, np.diag(self.synthetic_dataset.Q.T @ self.approx_hessian @ self.synthetic_dataset.Q), label=label)

    @abstractmethod
    def gradient_processing(self, grad):
        pass

class SGDVanilla(SGD):

    def gradient_processing(self, grad):
        return grad

class SGDNoised(SGD):

    def gradient_processing(self, grad):
        return grad + np.random.normal(0, 1, size=self.DIM)

class SGDCompressed(SGD):

    def __init__(self, synthetic_dataset, compressor: CompressionModel) -> None:
        super().__init__(synthetic_dataset)
        self.compressor = compressor

    def gradient_processing(self, grad):
        return self.compressor.compress(grad)

class SGDSportisse(SGD):

    def __init__(self, synthetic_dataset, compressor: CompressionModel) -> None:
        super().__init__(synthetic_dataset)
        self.compressor = compressor

    def gradient_processing(self, grad):
        return grad

    def compute_stochastic_gradient(self, w, data, labels, index):
        """Can be used only in the MISSING VALUE MODE."""
        x, y = self.synthetic_dataset.X[index], self.synthetic_dataset.Y[index]
        p = self.synthetic_dataset.estimated_p
        g = x * (w @ x - y) - (1 - p) * np.diag(x ** 2) @ w
        return g

class SGDNaiveSparsification(SGDCompressed):

    def gradient_processing(self, grad):
        return self.compressor.compress(grad)

    def gradient_processing(self, grad):
        return grad

    def compute_stochastic_gradient(self, w, data, labels, index):
        x, y = data[index], labels[index]
        x = self.synthetic_dataset.sparsificator.compress(x)
        g = x * (w @ x - y)
        return g

reason: from numpy,from scipy


text:
"""Natural Policy Gradient Optimization."""
from dowel import logger, tabular
import numpy as np
import copy

import tensorflow as tf

from metarl.misc import tensor_utils as np_tensor_utils
from metarl.tf.algos import NPO

class RL2NPO(NPO):
    """Natural Policy Gradient Optimization.

    This is specific for RL^2
    (https://arxiv.org/pdf/1611.02779.pdf).

    NPO is used as an inner algorithm for RL^2. When sampling for RL^2,
    there are more than one environments to be sampled from. In the original
    implementation, within each task/environment, all rollouts sampled will be
    concatenated into one single rollout, and fed to the inner algorithm.
    Thus, returns and advantages are calculated across the rollout.

    """

    def __init__(self,
                 env_spec,
                 policy,
                 baseline,
                 scope=None,
                 max_path_length=500,
                 discount=0.99,
                 gae_lambda=1,
                 center_adv=True,
                 positive_adv=False,
                 fixed_horizon=False,
                 pg_loss='surrogate',
                 lr_clip_range=0.01,
                 max_kl_step=0.01,
                 optimizer=None,
                 optimizer_args=None,
                 policy_ent_coeff=0.0,
                 use_softplus_entropy=False,
                 use_neg_logli_entropy=False,
                 stop_entropy_gradient=False,
                 entropy_method='no_entropy',
                 flatten_input=True,
                 fit_baseline='before',
                 name='NPO'):
        self._fit_baseline = fit_baseline
        if fit_baseline != 'before':
            raise ValueError('RL2NPO only supports fit baseline before')
        super().__init__(env_spec=env_spec,
                         policy=policy,
                         baseline=baseline,
                         scope=scope,
                         max_path_length=max_path_length,
                         discount=discount,
                         gae_lambda=gae_lambda,
                         center_adv=center_adv,
                         positive_adv=positive_adv,
                         fixed_horizon=fixed_horizon,
                         pg_loss=pg_loss,
                         lr_clip_range=lr_clip_range,
                         max_kl_step=max_kl_step,
                         optimizer=optimizer,
                         optimizer_args=optimizer_args,
                         policy_ent_coeff=policy_ent_coeff,
                         use_softplus_entropy=use_softplus_entropy,
                         use_neg_logli_entropy=use_neg_logli_entropy,
                         stop_entropy_gradient=stop_entropy_gradient,
                         entropy_method=entropy_method,
                         flatten_input=flatten_input)

    def optimize_policy(self, itr, samples_data):
        """Optimize policy.

        Args:
            itr (int): Iteration number.
            samples_data (dict): Processed sample data.
                See process_samples() for details.

        """
        if self._fit_baseline == 'before':
            self._fit_baseline_with_data(samples_data)
            samples_data['baselines'] = self._get_baseline_prediction(
                samples_data)

        policy_opt_input_values = self._policy_opt_input_values(samples_data)
        # Train policy network
        logger.log('Computing loss before')
        loss_before = self._optimizer.loss(policy_opt_input_values)
        logger.log('Computing KL before')
        policy_kl_before = self._f_policy_kl(*policy_opt_input_values)
        logger.log('Optimizing')
        self._optimizer.optimize(policy_opt_input_values)
        logger.log('Computing KL after')
        policy_kl = self._f_policy_kl(*policy_opt_input_values)
        logger.log('Computing loss after')
        loss_after = self._optimizer.loss(policy_opt_input_values)
        tabular.record('{}/LossBefore'.format(self.policy.name), loss_before)
        tabular.record('{}/LossAfter'.format(self.policy.name), loss_after)
        tabular.record('{}/dLoss'.format(self.policy.name),
                       loss_before - loss_after)
        tabular.record('{}/KLBefore'.format(self.policy.name),
                       policy_kl_before)
        tabular.record('{}/KL'.format(self.policy.name), policy_kl)
        pol_ent = self._f_policy_entropy(*policy_opt_input_values)
        tabular.record('{}/Entropy'.format(self.policy.name), np.mean(pol_ent))

        if self._fit_baseline == 'after':
            self._fit_baseline_with_data(samples_data)

        ev = np_tensor_utils.explained_variance_1d(samples_data['baselines'],
                                                   samples_data['returns'],
                                                   samples_data['valids'])
        tabular.record('{}/ExplainedVariance'.format(self.baseline.name), ev)

    def _get_baseline_prediction(self, samples_data):
        paths = samples_data['paths']
        baselines = [self.baseline.predict(path) for path in paths]
        return np_tensor_utils.pad_tensor_n(baselines, self.max_path_length)

ext: py

reason: import numpy


text:
# -*- coding: utf-8 -*-
"""
Created on Thu Oct 30 12:31:07 2014

@author: wangronin

"""

import pdb
import sys
import numpy as np
from copy import deepcopy

from pyDOE import lhs

from owck import OWCK as GaussianProcess

from scipy.stats import norm
from numpy.random import rand
from numpy import ones, array, sqrt, nonzero, min, max, sum, inf

# optimizers
from cma_es_class import cma_es
from scipy.optimize import fmin_l_bfgs_b

import random
import warnings

normcdf, normpdf = norm.cdf, norm.pdf

def ei(model, plugin=None):

    def __ei(X):

        X = np.atleast_2d(X)

        X = X.T if X.shape[1] != model.X.shape[1] else X

        n_sample = X.shape[0]

        if True:
            #here you did de-standardization? Not needed with OWCK
            y = model.y# * np.std(model.y) + np.mean(model.y)
        else:
            y = model.y * model.y_std + model.y_mean

        fmin = min(y) if plugin is None else plugin

        y_pre = []
        mse = []
        for sample in X:
            y_sample, mse_sample = model.predict(sample, eval_MSE=True)
            y_pre.append(y_sample[0])
            mse.append(mse_sample[0])

        y_pre = np.array(y_pre)
        mse = np.array(mse)
        y_pre = y_pre.reshape(n_sample)
        mse = mse.reshape(n_sample)

        sigma = sqrt(mse)

        value = (fmin - y_pre) * normcdf((fmin - y_pre) / sigma) + \
            sigma * normpdf((fmin - y_pre) / sigma)

        return value

    return __ei

def ei_dx(model, plugin=None):

    def __ei_dx(X):

        X = np.atleast_2d(X)

        X = X.T if X.shape[1] != model.X.shape[1] else X

        if True:
            #here you did de-standardization? Not needed with OWCK
            y = model.y# * np.std(model.y) + np.mean(model.y)
        else:
            y = model.y * model.y_std + model.y_mean

        fmin = min(y) if plugin is None else plugin

        y, sd2 = model.predict(X, eval_MSE=True)
        sd = np.sqrt(sd2)

        y_dx, sd2_dx = model.gradient(X)
        sd_dx = sd2_dx / (2. * sd)

        xcr = (fmin - y) / sd
        xcr_prob, xcr_dens = normcdf(xcr), normpdf(xcr)

        grad = -y_dx * xcr_prob + sd_dx * xcr_dens

        return grad

    return __ei_dx

def pi_dx(self, x):
    pass

def minei(model, dim, plugin=None):

    def __minei(X):

        X = np.atleast_2d(X)

        X = X.T if X.shape[0] == dim else X

        n_sample = X.shape[0]

        if True:
            #here you did de-standardization? Not needed with OWCK
            y = model.y# * np.std(model.y) + np.mean(model.y)
        else:
            y = model.y * model.y_std + model.y_mean

        fmin = min(y) if plugin is None else plugin

        y_pre = []
        mse = []
        for sample in X:
            y_sample, mse_sample = model.predict(sample, eval_MSE=True)
            y_pre.append(y_sample[0])
            mse.append(mse_sample[0])

        #print "y",y_pre
        y_pre = np.array(y_pre)
        mse = np.array(mse)
        y_pre = y_pre.reshape(n_sample)
        mse = mse.reshape(n_sample)
        #print "y",y_pre
        #print "mse", mse

        sigma = sqrt(mse)

        #print (fmin - y_pre) / sigma

        value = (fmin - y_pre) * normcdf((fmin - y_pre) / sigma) + \
            sigma * normpdf((fmin - y_pre) / sigma)
        #print value

        return -1*value

    return __minei

def pi(model, dim, plugin=None):

    def __pi(X):

        X = np.atleast_2d(X)

        X = X.T if X.shape[0] == dim else X

        y = model.y * model.y_std + model.y_mean

        fmin = min(y) if plugin is None else plugin

        y_pre, mse = model.predict(X, eval_MSE=True)
        sigma = sqrt(mse)

        value = (fmin - y_pre) * normcdf((fmin - y_pre) / sigma) + \
            sigma * normpdf((fmin - y_pre) / sigma)

        return value

    return __pi

class ego:

    def __init__(self,
                 dim,
                 n_init_sample,
                 fitness,
                 n_step,
                 lb=None,
                 ub=None,
                 gp_model=None,
                 criterion='EI',
                 strategy='niching',
                 sampling='LHS',
                 parallel_eval=False,
                 is_minimize=True,
                 useOWCK=True,
                 method="tree",
                 n_cluster=5,
                 solver="CMA",
                 verbose=False,
                 random_seed=999
                 ):

        self.n_init_sample = n_init_sample
        self.fitness = fitness if is_minimize else lambda x: -fitness(x)
        self.dim = dim
        self.sampling = sampling
        self.is_minimize = is_minimize
        self.parallel_eval = parallel_eval
        self.useOWCK = useOWCK
        self.method = method
        self.n_cluster = n_cluster
        self.solver = solver
        self.random_seed = random_seed
        self.criterion = criterion

        # restart upper bound for L-BFGS-B algorithm
        # TODO: this number should be related to the dimensionality
        self.n_restarts_optimizer = 50

        self.n_step = n_step     # num
        self.mode = 0            # '0' mode stands for sequential evaluation
                                 # '1' mode stands for parallel evaluation
        if q is not None:
            self.q = q
            self.mode = 1

        x_lb = np.atleast_2d(lb)
        x_ub = np.atleast_2d(ub)

        self.x_lb = x_lb if x_lb.shape[0] == 1 else x_lb.T
        self.x_ub = x_ub if x_ub.shape[0] == 1 else x_ub.T

        if gp_model is None: # if no GPR model provided
        	self.DOE()
        	self.__model_creation(X, y)

        # TODO
        # create the infill-criterion

        self.strategy = strategy if self.mode == 1 else None

        if (self.mode == 1) and (not strategy in ['CLmin', 'CLmax', 'niching']):
            raise Exception('Not supported strategy for multi-point EGO!')

        self.verbose = verbose

    def DOE(self):
    	# get old random state
        old_state1 = np.random.get_state()
        old_state2 = random.getstate()

        # set the new seed
        random.seed(self.random_seed)
        np.random.seed(self.random_seed)

        if self.sampling == 'LHS':
            X = lhs(dim, samples=n_init_sample, criterion='cm') * (self.x_ub - \
                self.x_lb) + self.x_lb
        elif self.sampling == 'uniform':
            X = rand(n_init_sample, dim) * (self.x_ub - self.x_lb) + \
                self.x_lb

        # restore the random state
        np.random.set_state(old_state1)
        random.setstate(old_state2)

        y = self.__parallel_evaluation(X) if self.parallel_eval else \
            array([self.fitness(x) for x in X])

        # create the meta model
        self.X, self.y = X, y

    def __model_creation(self, X, y):

        # bounds of hyper-parameters for optimization (not needed for normalization)
        thetaL = 1e-3 * ones(self.dim) * (self.x_ub - self.x_lb)
        thetaU = 10 * ones(self.dim) * (self.x_ub - self.x_lb)

        # initial search point for hyper-parameters
        theta0 = rand(1, self.dim) * (thetaU - thetaL) + thetaL
        if self.useOWCK:
            self.model = OWCK(regr='constant', corr='matern',
                 n_cluster=self.n_cluster, min_leaf_size=0,
                 cluster_method=self.method, overlap=0.0, beta0=None,
                 storage_mode='full', verbose=False, theta0=theta0, thetaL=thetaL,
                 thetaU=thetaU, optimizer='fmin_cobyla', random_start=1, nugget=0.000001,
                 normalize=False,
                 is_parallel=False)
        else:
            self.model = GaussianProcess(regr='constant',
                                         corr='matern',
                                         theta0=theta0,
                                         thetaL=thetaL,
                                         normalize=False,
                                         nugget=0.000001,
                                         thetaU=thetaU)
        # model fitting
        self.model.fit(X, y)

    def __parallel_evaluation(self, X):

        from mpi4py import MPI
        from pandas import DataFrame

        comm = MPI.COMM_WORLD

        n_procs, _ = X.shape

        # Spawning processes to test kriging mixture
        comm = MPI.COMM_SELF.Spawn(sys.executable, args=['ego_evaluation.py'],
                                   maxprocs=n_procs)

        # scatter the models and data
        comm.bcast(self.fitness, root=MPI.ROOT)
        comm.scatter([(k, X[k, :]) \
            for k in range(n_procs)], root=MPI.ROOT)

        # Synchronization while the children process are performing
        # heavy computations...
        comm.Barrier()

        # Gether the fitted model from the childrenn process
        # Note that 'None' is only valid in master-slave working mode
        result = comm.gather(None, root=MPI.ROOT)

        comm.Disconnect()

        # register the measures
        data = DataFrame([[d['index'], d['y']] \
            for d in result], columns=['index', 'y'])
        data.sort('index', inplace=True)

        return array(data['y'])

    def optimize(self):
        y_hist = []
        if self.mode == 0:
            best_x, best_y, new_X, new_y, y_hist = self.nstep()
        elif self.mode == 1:
            best_x, best_y, new_X, new_y = self.nstep_parallel()

        return best_x, best_y, new_X, new_y, y_hist

    def nstep(self):
        y_hist = []
        for i in range(self.n_step):
            new_x = self.max_ei(self.solver, min(self.y))
            new_x = self.__remove_duplicate(new_x)

            if len(new_x) == 0:
                y_hist.append( min(self.y))
                continue

            new_y = []
            for x in new_x:
                new_y.append(self.fitness(x)[0])
            new_y = np.array(new_y)
            #print "found new x", new_x
            #print "with value:",new_y

            self.X = np.r_[self.X, new_x]
            self.y = np.r_[self.y, new_y]

            y_hist.append(min(self.y))

            if self.useOWCK and self.method=="tree":
                self.model.updateModel(np.atleast_2d(new_x),new_y)
                #print(len(self.model.y))
            else:
                self.__model_creation(self.X, self.y)

            #print ('step {}, best fitness {}'.format(i+1, min(self.y)) )

        # set up the best solution found
        idx = nonzero(self.y == min(self.y))[0]
        best_x = self.X[idx, :]
        best_y = self.y[idx]
        best_y = best_y if self.is_minimize else -best_y

        return best_x, best_y, self.X[self.n_init_sample:, :], \
            self.y[self.n_init_sample:], y_hist

    def nstep_parallel(self):

        X_per_step = []
        y_per_step = []
        for i in range(self.n_step):

            if self.strategy == 'CL':
                new_x, _ = self.__parallel_CL()

            elif self.strategy == 'niching':
                new_x, _ = self.__parallel_niching()

            new_x = self.__remove_duplicate(new_x)
            if new_x is not []:
                new_y = self.__parallel_evaluation(new_x) if self.parallel_eval \
                    else array([self.fitness(x) for x in new_x])
            else:
                new_y = []

            if len(new_y) == 0:
                continue

            X_per_step.append(new_x)
            y_per_step.append(new_y)

            self.X = np.r_[self.X, new_x]
            self.y = np.r_[self.y, new_y]

            self.__model_creation(self.X, self.y)

            print('step {}, best fitness {}'.format(i+1, min(new_y)))

        # set up the best solution found
        idx = nonzero(self.y == min(self.y))
        best_x = self.X[idx, :]
        best_y = self.y[idx]
        best_y = best_y if self.is_minimize else -best_y

        return best_x, best_y, X_per_step, y_per_step

    def __remove_duplicate(self, new_x):

        # TODO: show a warning here
        new_x = np.atleast_2d(new_x)
        samples = []
        for x in new_x:
            if not any(sum(np.isclose(self.X, x), axis=1)):
                samples.append(x)

        return array(samples)

    def __parallel_CL():
        pass

    def __parallel_niching(self):

        from niching_dr1 import niching_dr1

        x_lb, x_ub = self.x_lb, self.x_ub

        q = self.q
        q_eff = q + 3
        rho = sqrt(sum((x_ub - x_lb) ** 2, axis=1)) / 2. / q**(1./self.dim)
        kappa = 1e2
        sigma0 = 0.1 * max(x_ub - x_lb)
        eval_budget = q_eff*1e3
        fit_ei = ei(self.model, self.dim)

        peaks, fit = niching_dr1(fit_ei, self.dim, x_lb.T, x_ub.T, q, q_eff, rho,
                            kappa, sigma0, eval_budget, is_minimize=False)

        return peaks.T, fit

    def __get_tree_boundary(self):
        """
        Get the variable boundary for each leaf node
        """

        def recurse(tree, node_id, decision_path, bounds):
            """
            Traverse the decision tree in preoder to retrieve the
            decision path for leaf nodes
            """
            decision_path_ = deepcopy(decision_path)
            left_child = tree.children_left[node_id]
            right_child = tree.children_right[node_id]
            split_feature = tree.feature[node_id]
            threshold = tree.threshold[node_id]

            if split_feature == -2: # leaf node
                leaf_node_path[node_id] = decision_path_
                leaf_bounds[node_id] = bounds
                return
            else:
                decision_path_.append(node_id) # pre-order visit
                __ = bounds[split_feature]

                # proceed to the left child
                bounds_left = deepcopy(bounds)
                bounds_left[split_feature] = update_bounds(__, leq=threshold)
                recurse(tree, left_child, decision_path_, bounds_left)

                # proceed to the right child
                bounds_right = deepcopy(bounds)
                bounds_right[split_feature] = update_bounds(__, geq=threshold)
                recurse(tree, right_child, decision_path_, bounds_right)

        def update_bounds(bounds, leq=None, geq=None):
            lb, ub = bounds
            if leq is not None:
                if leq <= lb:
                    raise Exception
                ub = leq if leq < ub else ub

            if geq is not None:
                if geq >= ub:
                    raise Exception
                lb = geq if geq > lb else lb

            return [lb, ub]

        leaf_node_path = {}
        leaf_bounds = {}
        tree = self.model.clusterer.tree_

        recurse(tree, 0, [],
                {i: [self.x_lb[0, i], self.x_ub[0, i]] for i in range(self.dim)})

        return leaf_bounds

    def max_ei(self, solver='BB', plugin=None):

        if plugin is None:
            plugin = min(self.y)

        # CMA-ES with Decision tree search
        if self.solver == 'CMA-tree':
            leaf_bounds = self.__get_tree_boundary()
            total_budget = 1e5*self.dim
            total_measure = np.prod(self.x_ub - self.x_lb)

            # TODO: this part should be parallelizable
            fopt = -inf
            for leaf_node, bounds in leaf_bounds.iteritems():

                x_lb, x_ub = zip(*[bounds[i] for i in bounds.keys()])
                x_lb = np.atleast_2d(x_lb)
                x_ub = np.atleast_2d(x_ub)

                measure = np.prod(x_ub - x_lb)
                eval_budget = int(measure / total_measure * total_budget)

                # Algorithm parameters
                opt = {
                       'sigma_init': 0.25 * np.max(x_lb - x_ub),
                       'eval_budget': eval_budget,
                       'f_target': -inf,
                       'lb': x_lb.T,
                       'ub': x_ub.T
                       }

                init_wcm = rand(1, self.dim) * (x_ub - x_lb) + x_lb
                fitnessfunc = ei(self.model)

                optimizer = cma_es(self.dim, init_wcm, fitnessfunc, opt,
                                   is_minimize=False)
                xopt_, fopt_, evalcount, _ = optimizer.optimize()

                if fopt_ > fopt:
                    fopt = fopt_
                    xopt = xopt_.T

            return xopt

        #standard scipy optimizers:
        if solver == "diff":
            from scipy.optimize import differential_evolution
            fitnessfunc = minei(self.model, self.dim, plugin)
            bounds = []
            for i in range(len(self.x_lb[0])):
                bounds.append([self.x_lb[0][i],self.x_ub[0][i]])
            #print bounds
            testfopt = fitnessfunc(np.ones(self.dim)*0)
            result = differential_evolution(fitnessfunc, bounds, maxiter=10000)
            #print "x:",result.x
            print "fopt", result.fun
            print "testfopt", testfopt
            print result.x
            return [result.x]

        elif solver == "basinhopping":

            from scipy.optimize import basinhopping
            fitnessfunc = minei(self.model, self.dim, plugin)
            bounds = []
            for i in range(len(self.x_lb[0])):
                bounds.append([self.x_lb[0][i],self.x_ub[0][i]])
            #print bounds
            testfopt = fitnessfunc(np.ones(self.dim)*0)
            minimizer_kwargs = {"method": "BFGS"}
            ri = random.randint(0, len(self.X))
            result = basinhopping(fitnessfunc,self.X[ri],  minimizer_kwargs=minimizer_kwargs, niter=5000)
            #print "x:",result.x
            print "fopt", result.fun
            print "testfopt", testfopt
            print result.x
            return [result.x]

        elif solver =="pso":
            from pyswarm import pso
            fitnessfunc = minei(self.model, self.dim, plugin)
            testfopt = fitnessfunc(np.ones(self.dim)*0)
            xopt, fopt = pso(fitnessfunc, self.x_lb, self.x_ub)
            print "fopt", fopt
            print "testfopt",testfopt
            print xopt
            return [xopt]

        elif solver == 'BB':   # Branch and Bound solver
            # TODO: implement this for Gaussian kernel
            pass
        elif solver == 'genoud':    # Genetic optimization with derivatives
            # TODO: somehow get it working from R package
            pass
        elif solver == 'BFGS':      # quasi-newton method with restarts

            EI, EI_dx = ei(self.model), ei_dx(self.model)
            def obj_func(x): return -EI(x), -EI_dx(x)

            fopt = np.inf
            eval_budget = 1000*self.dim
            bounds = np.c_[self.x_lb.T, self.x_ub.T]

            if not np.isfinite(bounds).all() and self.n_restarts_optimizer > 1:
                raise ValueError(
                        "Multiple optimizer restarts (n_restarts_optimizer>0) "
                        "requires that all bounds are finite.")

            # L-BFGS-B algorithm with restarts
            for iteration in range(self.n_restarts_optimizer):

                x0 = np.random.uniform(bounds[:, 0], bounds[:, 1])

                xopt_, fopt_, convergence_dict = \
                    fmin_l_bfgs_b(obj_func, x0, pgtol=1e-30, factr=1e5,
                                  bounds=bounds, maxfun=eval_budget)

                if convergence_dict["warnflag"] != 0:
                    warnings.warn("fmin_l_bfgs_b terminated abnormally with the "
                          " state: %s" % convergence_dict)

                xopt_list.append(xopt_)

                if fopt_ < fopt:
                    if self.verbose:
                        print 'iteration: ', iteration+1, convergence_dict['funcalls']

                    xopt, fopt = xopt_, fopt_

                eval_budget -= convergence_dict['funcalls']

                if eval_budget <= 0:
                    break

        elif solver == 'CMA':        # CMA-ES solver

            # TODO: use BI-POP/ IPOP-CMA-ES here
            # Algorithm parameters
            opt = { \
                'sigma_init': 0.25 * np.max(self.x_lb - self.x_ub),
                'eval_budget': 1e5*self.dim,
                'f_target': -inf,
                'lb': self.x_lb.T,
                'ub': self.x_ub.T \
                }

            init_wcm = rand(1, self.dim) * (self.x_ub - self.x_lb) + \
                self.x_lb

            fitnessfunc = ei(self.model)
            testfopt = fitnessfunc(np.ones(self.dim)*0)
            optimizer = cma_es(self.dim, init_wcm, fitnessfunc, opt,
                               is_minimize=False) #maximize ei
            xopt, fopt, evalcount, _ = optimizer.optimize()
            xopt = xopt.T

        return xopt

xopt_list = []
if __name__ == '__main__':

    from utils import plot_contour_gradient
    import matplotlib.pyplot as plt

#    plt.style.use('ggplot')
    np.random.seed(370)

    fig_width = 22
    fig_height = fig_width * 9 / 16

    # To verify the corr gradient computation
    x_lb = [-5, -5]
    x_ub = [5, 5]

    n_sample = 10

    def fitness(X):
        X = np.atleast_2d(X)
        return np.sin(np.sum(X, axis=1) ** 2.0)

    X = np.random.rand(n_sample, 2) * n_sample - 5
    y = fitness(X)

    model = GaussianProcess(corr='matern', cluster_method='tree',
                            n_cluster=3, theta0=[1e-1, 5*1e-1], normalize=False)
    model.fit(X, y)

    # Plot EI and PI landscape
    fig0, ax0 = plt.subplots(1, 1, figsize=(fig_width, fig_height),
                                    subplot_kw={'aspect':'equal'}, dpi=100)

    ax0.set_xlim([x_lb[0], x_ub[0]])
    ax0.set_ylim([x_lb[1], x_ub[1]])

    f = ei(model)
    grad = ei_dx(model)

    f_data = np.load('fitness.npy')
    grad_data = np.load('grad.npy')

    plot_contour_gradient(ax0, f, grad, x_lb, x_ub, title='Expected Improvement',
                          n_per_axis=100)

    # diagnostic comparisons of EI maximization methods
    color = ['k', 'g']
    for i, method in enumerate(['CMA-tree', 'BFGS']):
        optimizer = ego(2, 10,
                        fitness, 1,
                        lb=x_lb,
                        ub=x_ub,
                        useOWCK=False,
                        n_cluster=1,
                        solver=method,
                        verbose=True,
                        random_seed=0)

        optimizer.model = model
        optimizer.X = model.X
        optimizer.y = model.y

        _, __, new_X, new_y, ___ = optimizer.optimize()

        ax0.plot(new_X[-1, 0], new_X[-1, 1], ls='none', marker='.', ms=15,
                 mfc=color[i], mec='none', alpha=0.7)

        if len(xopt_list) != 0:
            __ = np.array(xopt_list)
            ax0.plot(__[:, 0], __[:, 1], ls='none', marker='.', ms=15,
                 mfc=color[i], mec='none', alpha=0.7, clip_on=False)

    plt.show()

max_stars_repo_name: Basvanstein/knn-variance

max_issues_repo_name: Basvanstein/knn-variance

max_forks_repo_name: Basvanstein/knn-variance

reason: import numpy,from numpy,from scipy


text:
"""An augmentation that adds random noise to an image."""
import random
from math import floor
from enum import Enum
import numpy as np
from discolight.params.params import Params
from .augmentation.types import ColorAugmentation, BoundedNumber
from .decorators.accepts_probs import accepts_probs

class NoiseType(Enum):

    """The noise type to be applied in the SaltAndPepperNoise augmentation."""

    RGB = "RGB"
    SnP = "SnP"

@accepts_probs
class SaltAndPepperNoise(ColorAugmentation):

    """Add salt and pepper or RGB noise to the given image."""

    def __init__(self, replace_probs, pepper, salt, noise_type):
        """Construct a SaltAndPepperNoise augmenation.

        You should probably use the augmentation factory or Discolight
        library interface to construct augmentations. Only invoke
        this constructor directly if you know what you are doing.
        """
        super().__init__()

        self.replace_probs = replace_probs
        self.pepper = pepper
        self.salt = salt
        self.noise_type = noise_type

    @staticmethod
    def params():
        """Return a Params object describing constructor parameters."""
        return Params().add("replace_probs", "", float,
                            0.1).add("pepper", "The color of the pepper",
                                     BoundedNumber(int, 0, 255),
                                     0).add("salt", "The color of the salt",
                                            BoundedNumber(int, 0, 255),
                                            255).add("noise_type",
                                                     "The type of noise",
                                                     NoiseType, "RGB")

    def augment_img(self, img, _bboxes):
        """Augment an image."""
        np.random.seed(floor(random.random() * 1000000))

        if self.noise_type == "SnP":
            random_matrix = np.random.rand(img.shape[0], img.shape[1])
            img[random_matrix >= (1 - self.replace_probs)] = self.salt
            img[random_matrix <= self.replace_probs] = self.pepper
        elif self.noise_type == "RGB":
            random_matrix = np.random.rand(img.shape[0], img.shape[1],
                                           img.shape[2])
            img[random_matrix >= (1 - self.replace_probs)] = self.salt
            img[random_matrix <= self.replace_probs] = self.pepper
        return img

max_stars_repo_name: tingalinga/discolight

max_issues_repo_name: tingalinga/discolight

max_forks_repo_name: tingalinga/discolight


text:
import numpy as np
import tensorflow as tf
from algos.tf1.offPolicy.baseOffPolicy import OffPolicy
from algos.tf1.sac_auto import core
from algos.tf1.sac_auto.core import get_vars

class SAC_AUTO(OffPolicy):
    def __init__(self,
                 act_dim, obs_dim, a_bound,
                 mlp_actor_critic=core.mlp_actor_critic,
                 ac_kwargs=dict(),
                 seed=0,
                 replay_size=int(1e6),
                 gamma=0.9,
                 polyak=0.99,
                 pi_lr=1e-3, q_lr=1e-3,
                 act_noise=0.1, target_noise=0.2,
                 noise_clip=0.5, policy_delay=2,
                 sess_opt=None,
                 sess=None,
                 batch_size=256,
                 buffer=None,
                 per_flag=True,
                 her_flag=True,
                 goal_selection_strategy="future",
                 n_sampled_goal=4,
                 action_l2=0.0,
                 clip_return=None,
                 state_norm=True,
                 alpha='auto',
                 ):
        super(SAC_AUTO, self).__init__(act_dim, obs_dim, a_bound,
                                       mlp_actor_critic=core.mlp_actor_critic,
                                       ac_kwargs=ac_kwargs, seed=seed,
                                       replay_size=replay_size, gamma=gamma, polyak=polyak,
                                       pi_lr=pi_lr, q_lr=q_lr, batch_size=batch_size, act_noise=act_noise,
                                       target_noise=target_noise, noise_clip=noise_clip,
                                       policy_delay=policy_delay, sess_opt=sess_opt,
                                       per_flag=per_flag, her_flag=her_flag,
                                       goal_selection_strategy=goal_selection_strategy,
                                       n_sampled_goal=n_sampled_goal, action_l2=action_l2,
                                       clip_return=clip_return, state_norm=state_norm)

        # Main outputs from computation graph
        with tf.variable_scope('main'):
            self.mu, self.pi, logp_pi, logp_pi2, q1, q2, q1_pi, q2_pi, = mlp_actor_critic(self.x_ph,
                                                                                self.x2_ph,
                                                                                self.a_ph,
                                                                                **ac_kwargs)

        # Target value network
        with tf.variable_scope('target'):
            _, _, logp_pi_, _, _, _, q1_pi_, q2_pi_ = mlp_actor_critic(self.x2_ph,
                                                                       self.x2_ph,
                                                                       self.a_ph,
                                                                       **ac_kwargs)

        # Count variables
        var_counts = tuple(core.count_vars(scope) for scope in
                           ['main/pi', 'main/q1', 'main/q2', 'main/v', 'main'])
        print(('\nNumber of parameters: \t pi: %d, \t' + \
               'q1: %d, \t q2: %d, \t v: %d, \t total: %d\n') % var_counts)
        # 重新修改下面这段!
        target_entropy = (-np.prod(act_dim))

        log_alpha = tf.get_variable('log_alpha', dtype=tf.float32, initializer=0.0)
        alpha = tf.exp(log_alpha)

        alpha_loss = tf.reduce_mean(-log_alpha * tf.stop_gradient(logp_pi + target_entropy))

        alpha_optimizer = tf.train.AdamOptimizer(learning_rate=1e-4,
                                                 name='alpha_optimizer')
        train_alpha_op = alpha_optimizer.minimize(loss=alpha_loss, var_list=[log_alpha])

        # Min Double-Q:
        min_q_pi = tf.minimum(q1_pi_, q2_pi_)

        # Targets for Q and V regression
        v_backup = tf.stop_gradient(min_q_pi - alpha * logp_pi2)
        q_backup = self.r_ph + gamma * (1 - self.d_ph) * v_backup

        # Soft actor-critic losses
        pi_loss = tf.reduce_mean(alpha * logp_pi - q1_pi)
        q1_loss = 0.5 * tf.reduce_mean((q_backup - q1) ** 2)
        q2_loss = 0.5 * tf.reduce_mean((q_backup - q2) ** 2)
        # value_loss = q1_loss + q2_loss
        if self.per_flag:
            # 也许可以选q2,但是一般来说q1和q2值相差不大。
            self.abs_errors = tf.abs(q_backup - q1)
            value_loss = self.ISWeights * (q1_loss + q2_loss)
        else:
            # 正常的！
            value_loss = q1_loss + q2_loss

        # Policy train op
        # (has to be separate from value train op, because q1_pi appears in pi_loss)
        pi_optimizer = tf.train.AdamOptimizer(learning_rate=self.actor_lr)
        train_pi_op = pi_optimizer.minimize(pi_loss, var_list=get_vars('main/pi'))

        # Value train op
        # (control dep of train_pi_op because sess.run otherwise evaluates in nondeterministic order)
        value_optimizer = tf.train.AdamOptimizer(learning_rate=self.critic_lr)
        value_params = get_vars('main/q')
        with tf.control_dependencies([train_pi_op]):
            train_value_op = value_optimizer.minimize(value_loss,
                                                      var_list=value_params)

        # Polyak averaging for target variables
        # (control flow because sess.run otherwise evaluates in nondeterministic order)
        with tf.control_dependencies([train_value_op]):
            target_update = tf.group([tf.assign(v_targ, polyak * v_targ + (1 - polyak) * v_main)
                                      for v_main, v_targ in zip(get_vars('main'), get_vars('target'))])

        # All ops to call during one training step
        if self.per_flag:
            self.step_ops = [pi_loss,
                             q1_loss, q2_loss,
                             q1, q2,
                             logp_pi, alpha,
                             train_pi_op,
                             train_value_op,
                             target_update,
                             train_alpha_op,
                             self.abs_errors]
        else:
            self.step_ops = [pi_loss,
                             q1_loss, q2_loss,
                             q1, q2,
                             logp_pi, alpha,
                             train_pi_op,
                             train_value_op,
                             target_update,
                             train_alpha_op]

        # Initializing targets to match main variables
        target_init = tf.group([tf.assign(v_targ, v_main)
                                for v_main, v_targ in zip(get_vars('main'), get_vars('target'))])

        if sess_opt:
            gpu_options = tf.GPUOptions(per_process_gpu_memory_fraction=sess_opt)
            self.sess = tf.Session(config=tf.ConfigProto(gpu_options=gpu_options))
        else:
            self.sess = tf.Session()
        self.sess.run(tf.global_variables_initializer())
        self.sess.run(target_init)

    def get_action(self, s, noise_scale=0):
        if self.norm is not None:
            s = self.norm.normalize(v=s)
        if not noise_scale:
            noise_scale = self.action_noise

        act_op = self.mu if noise_scale == 0 else self.pi
        a = self.sess.run(act_op, feed_dict={x_ph: o.reshape(1, -1)})[0]
        # a = self.sess.run(self.pi, feed_dict={self.x_ph: s.reshape(1, -1)})[0]
        a += noise_scale * np.random.randn(self.act_dim)
        return np.clip(a, -self.a_bound, self.a_bound)

    def learn(self, batch_size=100,
              actor_lr_input=0.001,
              critic_lr_input=0.001,
              ):

        if self.per_flag:
            tree_idx, batch_memory, ISWeights = self.replay_buffer.sample(batch_size=batch_size)
            batch_states, batch_actions, batch_rewards, batch_states_, batch_dones = [], [], [], [], []
            for i in range(batch_size):
                batch_states.append(batch_memory[i][0])
                batch_actions.append(batch_memory[i][1])
                batch_rewards.append(batch_memory[i][2])
                batch_states_.append(batch_memory[i][3])
                batch_dones.append(batch_memory[i][4])

            feed_dict = {self.x_ph: np.array(batch_states),
                         self.x2_ph: np.array(batch_states_),
                         self.a_ph: np.array(batch_actions),
                         self.r_ph: np.array(batch_rewards),
                         self.d_ph: np.array(batch_dones),
                         self.actor_lr: actor_lr_input,
                         self.critic_lr: critic_lr_input,
                         self.ISWeights: ISWeights
                         }

            outs = self.sess.run(self.step_ops,
                                 feed_dict)
            abs_errors = outs[-1]
            pi_loss, q1_loss, q2_loss, q1, q2, logp_pi, alpha, _, _, _,  train_alpha_op, abs_errors = outs
            self.replay_buffer.batch_update(tree_idx,
                                            abs_errors)  # update priority
            self.learn_step += 1
            return outs
        else:
            batch = self.replay_buffer.sample_batch(batch_size)
            feed_dict = {self.x_ph: batch['obs1'],
                         self.x2_ph: batch['obs2'],
                         self.a_ph: batch['acts'],
                         self.r_ph: batch['rews'],
                         self.d_ph: batch['done'],
                         self.actor_lr: actor_lr_input,
                         self.critic_lr: critic_lr_input,
                         }
            outs = self.sess.run(self.step_ops,
                                 feed_dict)
            self.learn_step += 1
            return outs

    def get_action(self, s, noise_scale=0):
        if self.norm is not None:
            s = self.norm.normalize(v=s)
        if not noise_scale:
            act_op = self.mu
        else:
            act_op = self.pi
        a = self.sess.run(act_op,
                          feed_dict={self.x_ph: s.reshape(1, -1)})[0]
        return np.clip(a, -self.a_bound, self.a_bound)

max_stars_repo_name: kaixindelele/DRLib


text:
# -*- coding: utf-8 -*-
"""
Created on Fri Sep 11 15:50:18 2020

@author: bbate
"""
#
#%% libraries
#
import pandas as pd
import numpy as np
from datetime import datetime
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde
from joblib import Parallel, delayed
#
#%% function def
#
def slow_function():
    for rep in range(5):
        x = pd.Series(np.random.normal(loc = 0, scale = 1, size = 1000))
        for i in x.index:
            z = x[i]**2
            p = np.sqrt(np.abs(x[i]))
            x[i] = z + p
        x = np.random.normal(loc = 5 * np.random.randint(1, 3), scale = 1) + x
    return x
#
#%% parallel function
#
def parallel_runs(n_trials):
#
    results = []
    for trial in range(n_trials):
           results = results + [slow_function()]
    return results
#
#%% parallel slow function demo
#
start = datetime.now()
parallel_processes = 8
trials_per_process = 3
parallel_results = (Parallel(n_jobs = -1)([delayed(parallel_runs)(trials_per_process)
                                           for _ in range(parallel_processes)]))
print('parallel time required: ', (datetime.now() - start).total_seconds())
#
serial_results = []
start = datetime.now()
for i in range(parallel_processes * trials_per_process):
    serial_results = serial_results + [slow_function()]
print('serial time required: ', (datetime.now() - start).total_seconds())
#
#%% vusualize
#
x_values = np.linspace(0, 25, 1000)
fig, ax = plt.subplots()
for job in parallel_results:
    for run in job:
        kde = gaussian_kde(run)
        density = kde(x_values)
        ax.fill(x_values, density, alpha = 0.3)
plt.show()
#
fig, ax = plt.subplots()
for job in serial_results:
    kde = gaussian_kde(job)
    density = kde(x_values)
    ax.fill(x_values, density, alpha = 0.3)
plt.show()
#

reason: import numpy,from scipy


text:
import numpy as np

def solve(expectancy, cutoff1, cutoff2):
    return 71

max_stars_repo_name: kylebegovich/MathPlayground

max_issues_repo_name: kylebegovich/MathPlayground

max_forks_repo_name: kylebegovich/MathPlayground


text:
# Copyright 2016 The TensorFlow Authors. All Rights Reserved.
#
# 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 imageio
import tensorflow.compat.v1 as tf
from spatial_transformer import ElasticTransformer
import numpy as np

# Input image retrieved from:
# https://raw.githubusercontent.com/skaae/transformer_network/master/cat.jpg
im = imageio.imread("data/cat.jpg")
im = im / 255.0
im = im.astype("float32")

# input batch
batch_size = 4
batch = np.expand_dims(im, axis=0)
batch = np.tile(batch, [batch_size, 1, 1, 1])

# Let the output size of the elastic transformer be quarter of the image size.
outsize = (int(im.shape[0] / 4), int(im.shape[1] / 4))

# Elastic Transformation Layer
stl = ElasticTransformer(outsize)
stlc = ElasticTransformer(outsize, interp_method="bicubic")
stlcm = ElasticTransformer(outsize, interp_method="bicubic", masked=True, cval=0.5)

# Run session
def main(x):
    # Random jitter of identity parameters
    theta = 0.1 * tf.random_normal([batch_size, stl.param_dim])
    result = stl.transform(x, theta)
    resultc = stlc.transform(x, theta)
    resultcm = stlcm.transform(x, theta)
    return result, resultc, resultcm

result_, resultc_, resultcm_ = main(batch)

# save our result
for i in range(result_.shape[0]):
    imageio.imsave("elastic" + str(i) + ".png", result_[i])
    imageio.imsave("elastic" + str(i) + "_bicubic.png", resultc_[i])
    imageio.imsave("elastic" + str(i) + "_bicubic_masked_grey.png", resultcm_[i])

max_stars_repo_name: bamford/spatial-transformer-tensorflow

max_issues_repo_name: bamford/spatial-transformer-tensorflow

max_forks_repo_name: bamford/spatial-transformer-tensorflow


text:
import xarray as xr
import numpy as np
import pytest

from xrspatial import aspect

def _do_sparse_array(data_array):
    import random
    indx = list(zip(*np.where(data_array)))
    pos = random.sample(range(data_array.size), data_array.size//2)
    indx = np.asarray(indx)[pos]
    r = indx[:, 0]
    c = indx[:, 1]
    data_half = data_array.copy()
    data_half[r, c] = 0
    return data_half

def _do_gaussian_array():
    _x = np.linspace(0, 50, 101)
    _y = _x.copy()
    _mean = 25
    _sdev = 5
    X, Y = np.meshgrid(_x, _y, sparse=True)
    x_fac = -np.power(X-_mean, 2)
    y_fac = -np.power(Y-_mean, 2)
    gaussian = np.exp((x_fac+y_fac)/(2*_sdev**2)) / (2.5*_sdev)
    return gaussian
#
# -----

data_random = np.random.random_sample((100, 100))
data_random_sparse = _do_sparse_array(data_random)
data_gaussian = _do_gaussian_array()

def test_aspect_transfer_function():
    """
    Assert aspect transfer function
    """
    da = xr.DataArray(data_gaussian, dims=['y', 'x'], attrs={'res': 1})
    da_aspect = aspect(da)
    assert da_aspect.dims == da.dims
    assert da_aspect.coords == da.coords
    assert da_aspect.attrs == da.attrs
    assert da.shape == da_aspect.shape
    assert pytest.approx(da_aspect.data.max(), .1) == 360.
    assert pytest.approx(da_aspect.data.min(), .1) == 0.

max_stars_repo_name: ocefpaf/xarray-spatial

max_issues_repo_name: ocefpaf/xarray-spatial

max_forks_repo_name: ocefpaf/xarray-spatial


text:
#!/usr/bin/env python
# -*-coding:utf-8-*-
# Author: nomalocaris <nomalocaris.top>
"""Geolife 轨迹提取相关函数"""
import random
import os
from utilis import cal_time_interval
from utilis import ProgressBar
import numpy as np

def extract_dataset(lat_lon,
                    good_traj_file,
                    init_path,
                    base_path,
                    rand_num=14650):
    """
    对数据进行抽取
    :param lat_lon:经纬度的界限值
    :param good_traj_file:用于存储Trajectory里面的文件名
    :param init_path:抽取后数据的存储路径
    :param base_path:原始数据
    :param rand_num:所有数据的数量,14650
    :return:
    """
    max_time = 0
    D = []
    file_list = []
    traj_fs = os.listdir(base_path)
    p1 = ProgressBar(len(traj_fs), '提取轨迹')
    for t in range(len(traj_fs)):  # 0-181名用户的轨迹数据
        p1.update(t)
        base_name = traj_fs[t]
        file_path = base_path + "/" + base_name + '/Trajectory/'
        files = os.listdir(file_path)  # 打开文件夹
        for f in files:  # 遍历文件
            line_count = 0  # 轨迹段计数
            # print(file_path + f)
            with open(file_path + f) as fr:
                # for line in fr.readlines():
                #     print(line)
                # print(type(fr.readlines()))
                T = []
                prev_time = 0
                lines = fr.readlines()
                for id in range(6, len(lines)):
                    now_line = lines[id].strip().split(',')
                    now_line_latitude = float(now_line[0])
                    now_line_longitude = float(now_line[1])
                    # 检验是否在划定范围内
                    if lat_lon[0] < now_line_latitude < lat_lon[1] \
                            and lat_lon[2] < now_line_longitude < lat_lon[3]:
                        now_line_time = now_line[6]
                        line_count += 1
                        if line_count > 1:
                            now_time = now_line_time
                            t_interval = cal_time_interval(prev_time, now_time)
                            if t_interval > max_time:
                                # 轨迹间隔超过2mins, 退出读取
                                max_time = t_interval
                            if t_interval > 6950:
                                # 轨迹间隔超过2mins, 退出读取
                                break
                        prev_time = now_line_time
                        T.append((now_line_latitude, now_line_longitude))
                    else:
                        # 经纬度超过范围
                        if line_count:
                            break
                        else:
                            continue
                if line_count < 3:
                    # 轨迹长度小于3，跳过
                    continue
                D.append(T)
                file_list.append(base_name + '_' + f)
    print("最大用时：", max_time)
    rand_ind = random.sample([i for i in range(len(file_list))], rand_num)  # 随机抽取轨迹
    with open(good_traj_file, 'w') as f2:
        p2 = ProgressBar(rand_num, '写入文件')
        for i in range(rand_num):
            p2.update(i)
            f2.writelines(file_list[rand_ind[i]] + '\n')
            with open(init_path + file_list[rand_ind[i]], 'w') as f3:
                for step in D[rand_ind[i]]:
                    f3.writelines(str(step[0]) + ',' + str(step[1]) + '\n')

def EDistant(m, n):
    """
    计算欧式距离
    :param m:数据1
    :param n:数据2
    :return:
    """
    return np.sqrt((m-n)[0]**2 + (m-n)[1]**2)

def check_data():
    length_ = []
    base_path_list = os.listdir('../data/Geolife Trajectories 1.3/Trajectories7000/')
    for path in base_path_list:
        file_object = open('../data/Geolife Trajectories 1.3/Trajectories7000/' + path, 'r')
        length_.append(len(file_object.readlines()))
    print(length_)
    print(np.mean(length_))
    print(np.std(length_, ddof=1))

if __name__ == '__main__':
    # extract_dataset([39.6, 40.8, 115.8, 117.4], '../data/Geolife Trajectories 1.3/good_traj_file.txt',
    #                 '../data/Geolife Trajectories 1.3/Trajectories/',
    #                 '../data/Geolife Trajectories 1.3/Geolife Trajectories 1.3/Data', 14650)
    check_data()

max_stars_repo_name: hahahahahzt/DP-Star

max_issues_repo_name: hahahahahzt/DP-Star

max_forks_repo_name: hahahahahzt/DP-Star


text:
import numpy as np
from grid_map import grid_map
from ir_sim.env import env_base
from pathlib import Path
from reinforcement_learning import reinforcement_learning
import sys
from pathlib import Path

animation = False

# path set
cur_path = sys.path[0]
map_path = cur_path + '/map_matrix.npy'
reward_path = cur_path + '/reward_matrix.npy'
image_path = Path(__file__).parent / 'image'
gif_path = Path(__file__).parent / 'gif'

# load map and reward matrix
map_matrix = np.load(map_path)
reward_matrix = np.load(reward_path)

# set environment and grid_map
env = env_base(world_width = 20, world_height = 20)
grid_map = grid_map(map_matrix=map_matrix, reward_matrix=reward_matrix, start_index=(2, 2))

# run Q learning algorithm
rl = reinforcement_learning(grid_map.state_space, grid_map.action_space, grid_map)
state_action_value = rl.Q_learning()

# policy performance test
cur_index = grid_map.start_index

for i in range(300):
    action_index = np.argmax( state_action_value[cur_index[0], cur_index[1], :])
    cur_index, _, _, done = grid_map.step(cur_index, action_index)
    grid_map.set_path(cur_index)
    grid_map.draw_map()

    if animation: env.save_fig(image_path, i)

    if done:
        print('done')
        break

if animation: env.save_ani(image_path, gif_path, ani_name='Q_learning', keep_len=10)

grid_map.show_map()

# print(policy_value)

max_stars_repo_name: Intelligent-Robot-Course/Homework

max_issues_repo_name: Intelligent-Robot-Course/Homework

max_forks_repo_name: Intelligent-Robot-Course/Homework


text:
from statsmodels.tsa.vector_ar.vecm import coint_johansen
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.api import VAR
import statsmodels.api as sm
import copy
import pingouin
from scipy.stats import pearsonr
p_val = .05

def adfuller_df(df, maxlag, regression = "c"):
    dct = {}
    for key, val in df.items():
        dct[key] = adfuller(val, maxlag=maxlag, autolag = "aic", regression = regression)[1]#.pvalue

    return dct

def cointegration_df(df):
    dct = {}
    for key1 in df:
        dct[key1] = {}
        for key2 in df:
            if key1 == key2:
                dct[key1][key2] = np.nan
            else:
                dct[key1][key2] = coint_johansen(df[[key1,key2]], det_order=0, k_ar_diff=1).lr1
    return pd.DataFrame(dct)

max_stars_repo_name: jlcatonjr/ProjectRepositories

max_issues_repo_name: jlcatonjr/ProjectRepositories

max_forks_repo_name: jlcatonjr/ProjectRepositories

reason: from scipy,import statsmodels,from statsmodels


text:
#!/usr/bin/env python

from collections import namedtuple
import networkx as nx
from networkx.algorithms.simple_paths import all_simple_paths
from networkx.drawing.nx_pydot import to_pydot

Route = namedtuple('Route', ['start', 'end', 'via', 'distance'])

graph = nx.MultiDiGraph()

class RouteTypes:
    ROAD = ('Road', 'red')
    MULTI_USE = ('Multi-use Trail', 'green')
    SHARROW = ('Sharrow', 'orange')
    GREENWAY = ('Greenway', 'orange')
    MAJOR_SEP = ('Road, Major Separation', 'yellow')
    MINOR_SEP = ('Road, Minor Separation', 'yellow')

def node(name, color):
    graph.add_node(name, color=color)

def point_node(name):
    graph.add_node()

def route(start, end, via, distance, route_type):
    distance = float(distance)
    graph.add_edge(start, end,
        via=via,
        weight=distance,
        route_type=route_type[0],
        color=route_type[1],
        label='{0}\n{1:0.1f}'.format(via, distance))

graph.add_node('Home', color='green')
for bridge in ['Ballard Locks', 'Ballard Bridge', 'Fremont Bridge', 'University Bridge', 'Montlake Bridge']:
    graph.add_node(bridge, color='blue')
for point in ["Fisherman's Terminal", 'Ballmer Yard', 'Diamond Marina', 'Fred Meyer']:
    graph.add_node(point, shape='point')
graph.add_node('Work', color='red')

route('Home', 'Golden Gardens', 'Golden Gardens Dr', 1.7, RouteTypes.ROAD)
route('Golden Gardens', 'Ballard Locks', 'B-G Trail', 1.9, RouteTypes.MULTI_USE)
route('Home', 'Ballard Locks', '32nd Ave', 2, RouteTypes.SHARROW)
route('Home', 'Ballard Locks', '28th Ave', 1.7, RouteTypes.ROAD)
route('Ballard Locks', "Fisherman's Terminal", 'Commodore Way', 1.1, RouteTypes.ROAD)
route("Fisherman's Terminal", 'Ballmer Yard', 'Gilman Ave', 1.1, RouteTypes.MAJOR_SEP)
route('Ballmer Yard', 'Myrtle Edwards', 'Elliot Bay Trail', 2.1, RouteTypes.MULTI_USE)
route('Myrtle Edwards', 'Work', 'Mercer St', 1.6, RouteTypes.ROAD)
route('Myrtle Edwards', 'Work', 'Broad St', 1.5, RouteTypes.ROAD)
route("Fisherman's Terminal", 'Diamond Marina', 'Ship Canal Trail', 2.3, RouteTypes.MULTI_USE)
route('Diamond Marina', 'Lake Union Park', 'Lake Union Loop', 1.5, RouteTypes.MULTI_USE)
route('Lake Union Park', 'Work', 'Terry Ave', 0.4, RouteTypes.ROAD)

route('Home', 'Ballard Bridge', '24th Ave', 2.5, RouteTypes.MINOR_SEP)
route('Home', 'Ballard Bridge', '17th Ave', 2.6, RouteTypes.SHARROW)
route('Ballard Bridge', 'Diamond Marina', 'Ship Canal Trail', 2.6, RouteTypes.MULTI_USE)
route('Ballard Bridge', 'Myrtle Edwards', 'Elliot Bay Trail', 3.5, RouteTypes.MULTI_USE)

route('Home', 'Fred Meyer', '24th Ave', 2.9, RouteTypes.MINOR_SEP)
route('Home', 'Fred Meyer', '17th Ave', 3.2, RouteTypes.SHARROW)
route('Home', 'Fred Meyer', '8th Ave', 3.2, RouteTypes.MINOR_SEP)

route('Fred Meyer', 'Fremont Bridge', 'B-G Trail', 1.2, RouteTypes.MULTI_USE)
route('Fremont Bridge', 'Diamond Marina', '-', 0.3, RouteTypes.MULTI_USE)
route('Fremont Bridge', 'Work', 'Dexter Ave', 2.4, RouteTypes.MINOR_SEP)

route('Fred Meyer', 'University Bridge', 'B-G Trail', 3.2, RouteTypes.MULTI_USE)
route('University Bridge', 'Lake Union Park', 'Eastlake Ave', 2.2, RouteTypes.MINOR_SEP)

route('Fred Meyer', 'Montlake Bridge', 'B-G Trail', 4.0, RouteTypes.MULTI_USE)
route('Montlake Bridge', 'Lake Union Park', 'Delmar Dr', 3.3, RouteTypes.MINOR_SEP)

route('Home', 'Greenlake', '77th St', 2.8, RouteTypes.ROAD)
route('Greenlake', 'University Bridge', 'Roosevelt Way', 2.9, RouteTypes.MINOR_SEP)
route('Greenlake', 'Fremont Bridge', 'Stone Way', 3.3, RouteTypes.MINOR_SEP)

num_paths = len(list(all_simple_paths(graph, 'Home', 'Work')))
print("Calculated {0} paths".format(num_paths))

def dfs(node, target, path=[]):
    if node == target:
        yield path
    for s, e, data in graph.edges(node, data=True):
        new_path = list(path)
        new_path.append((s, e, data))
        yield from dfs(e, target, new_path)

for path in dfs('Home', 'Work'):
     label = 'Home'
     distance = 0.
     for s, e, data in path:
         label += ' -> ' + data['via'] + ' -> ' + e
         distance += data['weight']
     print('{0:0.2f}\t{1}'.format(distance, label))

gv = to_pydot(graph)
gv.write('bike-to-work.dot')
gv.write_png('bike-to-work.png', prog='dot')

max_stars_repo_name: ungood/bike-to-work

max_issues_repo_name: ungood/bike-to-work

max_forks_repo_name: ungood/bike-to-work

reason: import networkx,from networkx


text:
import pandas as pd
import os
import numpy as np

def merge(enhancer_features_path, window_features_path,
		  patch_features_path,promoter2k_features_path,
		  complete_feature_save_path):
	enhancer_features_path_list = os.listdir(enhancer_features_path)
	window_features_path_list = os.listdir(window_features_path)
	patch_features_path_list = os.listdir(patch_features_path)
	promoter2k_features_path_list = os.listdir(promoter2k_features_path)

	abs_paths = []
	for enhancer_path in enhancer_features_path_list:
		abs_enhancer_path = enhancer_features_path+enhancer_path
		abs_paths.append(abs_enhancer_path)

	for window_path in window_features_path_list:
		abs_window_path = window_features_path+window_path
		abs_paths.append(abs_window_path)

	for patch_path in patch_features_path_list:
		abs_patch_path = patch_features_path+patch_path
		abs_paths.append(abs_patch_path)

	for promoter2k_path in promoter2k_features_path_list:
		abs_promoter2k_path =  promoter2k_features_path+promoter2k_path
		abs_paths.append(abs_promoter2k_path)

	complete_data = pd.read_csv(abs_paths[0])
	print(complete_data.shape)
	for abs_path in abs_paths[1:]:
		current_df = pd.read_csv(abs_path)
		#current_df = current_df.drop(['chrom-Enh','chromStart','chromEnd','TSS','label','id'],axis=1)
		#'''
		complete_data = complete_data.merge(current_df, how='left', on=['chrom-Enh',
																		'chromStart',
																		'chromEnd',
																		'TSS',
																		'label',
																		'id'])
		#'''
		#complete_data = pd.concat([complete_data,current_df],axis=1)
		print(complete_data.shape,abs_path)

	complete_data.to_csv(complete_feature_save_path, index=False)
	print(complete_data.shape)

def merge_enhancer(use_real_ep_folder, use_fake_ep_folder,
				  cellline, region):
	use_real_ep_files = os.listdir(use_real_ep_folder)
	use_fake_ep_files = os.listdir(use_fake_ep_folder)
	intersect = np.intersect1d(use_real_ep_files, use_fake_ep_files)

	all_merged = pd.DataFrame()
	for file in intersect:
		histon_name = file.split('_')[0]
		if region == 'enhancerfeatures':
			use_real_ep_columns = ['chrom-Enh','chromStart','chromEnd','TSS','label','id']
			enhancer_features_cols = [histon_name+'_'+x for x in ['enhancer_avg_sigval','enhancer_max_sigval','enhancer_min_sigval',
															 'enhancer_std_sigval','enhancer_var_sigval','enhancer_avg_peak',
															 'enhancer_max_peak','enhancer_min_peak','enhancer_std_peak','enhancer_var_peak',
													 		  'enhancer_peak_density']]
			use_real_ep_columns.extend(enhancer_features_cols)

			real_data = pd.read_csv(use_real_ep_folder+file,names=use_real_ep_columns)
		else:
			real_data = pd.read_csv(use_real_ep_folder+file)
		fake_data = pd.read_csv(use_fake_ep_folder+file)
		merged_tmp = fake_data.append(real_data)
		merged_tmp.to_csv('../merged-ef/'+cellline+'/'+region+'/'+file,index=False)
		print(merged_tmp.shape)

if __name__ == '__main__':
	'''
	merge_enhancer('../../data/GM12878/midfile/enhancerfeatures/','../data/GM12878/muti-midfile/enhancerfeatures/',
				'GM12878','enhancerfeatures')

	merge_enhancer('../../data/GM12878/midfile/windowfeatures/','../data/GM12878/muti-midfile/windowfeatures/',
				'GM12878','windowfeatures')
	merge_enhancer('../../data/GM12878/midfile/patchfeatures/','../data/GM12878/muti-midfile/patchfeatures/',
				'GM12878','patchfeatures')

	merge_enhancer('../../data/GM12878/midfile/promoter2k/','../data/GM12878/muti-midfile/promoter2k/',
				'GM12878','promoter2k')

	'''

	'''
	merge('../data/NHEK/muti-midfile/enhancerfeatures/',
		'../data/NHEK/muti-midfile/windowfeatures/',
		'../data/NHEK/muti-midfile/patchfeatures/',
		'../data/NHEK/muti-midfile/promoter2k/',
		'../data/NHEK/muti-midfile/train.csv')
	'''
	for item in ['NHEK']:#['k562','GM12878','IMR90','HUVEC','Hela-S3','NHEK']:
		merge('../data/{}/muti-midfile/enhancerfeatures/'.format(item),
			'../data/{}/muti-midfile/windowfeatures/'.format(item),
			'../data/{}/muti-midfile/patchfeatures/'.format(item),
			'../data/{}/muti-midfile/promoter2k/'.format(item),
			'../data/{}/muti-midfile/train.csv'.format(item))

max_stars_repo_name: zhouweiyg/EpPredictor

max_issues_repo_name: zhouweiyg/EpPredictor

max_forks_repo_name: zhouweiyg/EpPredictor


text:
# $Id: __init__.py,v 1.1 2010-10-07 15:57:29 wirawan Exp $
#
# wpylib.math.stats.array_stats module
# Created: 20140404
# Wirawan Purwanto
#

"""
wpylib.math.stats.array_stats module
Tools for studying the statistics of an array, the difference of
two or more arrays, etc.

"""

import os
import sys
import numpy

# Functions below were imported from V2b_inspect.py, 20140404

"""
CONVENTIONS

- Usually, M1 is the matrix under examination, M2 is the reference matrix.

PRINCIPLES

- Make sure functions are dimension-independent as much as possible
  (i.e. not limited to 1D or 2D only).
- Text output should not hardwired only to sys.stdout via vanilla print
  statement.
"""

def rms(m):
  return numpy.sqrt(numpy.sum(m**2) / m.size)

def report_diff_stat(M1, M2, out=sys.stdout):
  # Original function name: statdiff
  """Studies the difference of two arrays (or matrices).
  Usually, M1 is the matrix under examination, M2 is the reference
  matrix.
  Prints out standard report to a given text file, or returns
  the report as a string.
  """
  dM = M1 - M2
  if len(M1.shape) == 2:
    dM_diag = numpy.diagonal(dM)
    stats = (
      dM.max(),
      dM.min(),
      dM.mean(),
      numpy.abs(dM).mean(),
      rms(dM),
      dM_diag.max(),
      dM_diag.min(),
      dM_diag.mean(),
      numpy.abs(dM_diag).mean(),
      rms(dM_diag),
    )
    if out == tuple:
      return stats

    rslt = """\
- max difference  = %13.6e
- min difference  = %13.6e
- mean difference = %13.6e
- mean abs diff   = %13.6e
- rms difference  = %13.6e
- max difference  = %13.6e   on diagonal
- min difference  = %13.6e   on diagonal
- mean difference = %13.6e   on diagonal
- mean abs diff   = %13.6e   on diagonal
- rms difference  = %13.6e   on diagonal
""" % stats
  else:
    stats = (
      dM.max(),
      dM.min(),
      dM.mean(),
      numpy.abs(dM).mean(),
      rms(dM),
    )
    if out == tuple:
      return stats

    rslt = """\
- max difference  = %13.6e
- min difference  = %13.6e
- mean difference = %13.6e
- mean abs diff   = %13.6e
- rms difference  = %13.6e
""" % stats

  if out == str:
    return rslt
  else:
    out.write(rslt)
    out.flush()

def maxadiff(M1, M2):
  # Original function name: maxdiff
  """Returns the maximum absolute difference between two matrices.
  Used for checking whether two matrices are identical."""
  return numpy.abs(M1 - M2).max()

def maxadiff_by_col(M1, M2):
  # Original function name: maxdiff_by_col
  """Assuming two-dimensional array (where first index is the row
  index), computes the maximum absolute difference.
  """
  return numpy.abs(M1 - M2).max(axis=0)

def maxradiff(M1, M2, ztol=1e-15):
  # Original function name: maxrdiff
  """Returns the maximum *relative* absolute difference between two matrices.
  Used for checking whether two matrices are identical.

  CAVEATS:
  - Zero elements in M2 (reference matrix) are converted to unity for
    calculating the relative difference.
  """
  (mapprox, mref) = (M1, M2)
  diff = numpy.asarray(numpy.abs(mapprox - mref))
  mref0 = numpy.asarray(numpy.abs(mref))
  numpy.putmask(mref0, mref0 < ztol, [1.0])

  return (diff / mref0).max()

def rdiff(M1, M2, ztol=1e-15, Abs=False):
  """Returns the *relative* [absolute] difference between two matrices.
  Used for checking whether two matrices are identical.
  Here, M1 is the array being compared, and M2 is the reference array.

  CAVEATS:
  - Zero elements in M2 (reference matrix) are converted to unity for
    calculating the relative difference.
  """
  from numpy import abs, array, asarray
  (mapprox, mref) = (M1, M2)
  if Abs:
    diff = abs(asarray(mapprox) - asarray(mref))
    mref0 = abs(asarray(mref))
    numpy.putmask(mref0, mref0 < ztol, [1.0])
  else:
    diff = asarray(mapprox) - asarray(mref)
    mref0 = array(mref, copy=True)  # must make a copy
    numpy.putmask(mref0, abs(mref0) < ztol, [1.0])

  return (diff / mref0)

class ArrayStat(object):
  """A class to compute the statistics of an array.
  """
  def __init__(self, mat, save=False):
    mat = numpy.asarray(mat)
    amat = numpy.abs(mat)
    if save:
      self.mat = mat
    self.min = mat.min()
    self.max = mat.max()
    self.mean = mat.mean()
    self.amin = amat.min()
    self.amax = amat.max()
    self.amean = amat.mean()
    self.rms = numpy.sqrt(numpy.sum(amat**2) / amat.size)

  def report(self, out=sys.stdout):
    """Prints out a standard report.
    """
    #print_histogram(study_sparsity(self.diff, (log_delta-5.5,log_delta+4.5+0.1)), xaxis='left_edge')
    rslt = """\
. min       = %12.6g
. max       = %12.6g
. mean      = %12.6g
. absmin    = %12.6g
. absmax    = %12.6g
. absmean   = %12.6g
. rms       = %12.6g
""" % (self.min, self.max, self.mean,
       self.amin, self.amax, self.amean,
       self.rms)

    if out == str:
      return rslt
    else:
      out.write(rslt)
      out.flush()

def study_sparsity(M, binrange=(-16.5,6.5,1.0), kind='log10x', minzero=1e-30, args={}):
  # Adapted from V2b_inspect.py research module.
  """Study the distribution of values or sparsity of a given array, M,
  using histogram.

  Returns a tuple of length N; the first (N-2) elements are the
  histogram() output, while the last two are the number of elements
  whose values are smaller (greater) than the lower (upper) bound of
  the histogram bind range.

  Argument 'kind': By default, the 'x' values of the histogram is on
  logarithmic-10 scale.
  Other common log scales can be chosen using 'log2x' and 'logx'.
  Choose 'linx' if you want the linear scale instead.

  Argument 'binrange' can be one of the following:
  - 3-tuple: the left edge of the first bin, the right edge of the last bin,
    and the width of each bin
  - 2-tuple: the left MIDPOINT of the first bin, the right MIDPOINT
    of the last bin; the width is assumed to be 1.

  The defaults are a 3-tuple, intended for log10x scale, where
  we check order-of-magnitudes between 1e-16 and 1e+5 (roughly).
  """
  from numpy import abs, count_nonzero, arange, histogram, log10, log2, log
  if kind == 'log10x': # usual way, for broad view
    log_val = log10(abs(M.flatten()) + minzero)
  elif kind == 'log2x':
    log_val = log2(abs(M.flatten()) + minzero)
  elif kind == 'logx':
    log_val = log(abs(M.flatten()) + minzero)
  elif kind == 'linx': # linear x scale, usually for more detailed view
    log_val = abs(M.flatten())
  else:
    raise ValueError, "Invalid kind=%s" % (str(kind))
  if len(binrange) == 3:
    binedges = numpy.arange(*binrange)
  elif len(binrange) == 2:
    l = binrange[0]-0.5
    r = binrange[1]+0.5
    binedges = numpy.arange(l,r)
  else:
    raise ValueError, "Invalid binrange parameter value"
  #print binedges
  hist = histogram(log_val, bins=binedges, **args)
  # Count values outside the range being considered:
  l = 0
  r = 0
  llim = binedges[0]
  rlim = binedges[-1]
  l = count_nonzero(log_val < llim)
  r = count_nonzero(log_val > rlim)
  return hist + (l,r)

def print_histogram(hist, xaxis='midpt',
                    # output formatting
                    xwidth=3, xprec=0,
                    out=sys.stdout):
  # Adapted from V2b_inspect.py research module.
  """Prints histogram in an easier way to visualize in text fashion.
  """
  if len(hist) >= 4:
    # special case: for study_sparsity output.
    (bar,edges,lexcess,rexcess) = hist
  else:
    # for all other cases--
    (bar,edges) = hist[:2]
    lexcess = 0
    rexcess = 0
  if xaxis == 'midpt':
    midpt = (edges[1:] + edges[:-1]) * 0.5
  elif xaxis == 'leftpt':
    midpt = edges[:-1]
  elif xaxis == 'rightpt':
    midpt = edges[1:]
  else:
    raise ValueError, "Invalid xaxis parameter value"
  width = max(int(numpy.log10(max(bar.max(), lexcess, rexcess))+1), xwidth)
  #print "width = %d" % width
  barfmt = " %" + str(width) + "d"
  midfmt = " %" + str(width) + "." + str(xprec) + "f"
  if out == str:
    from StringIO import StringIO
    out = StringIO()
    def retval():
      return out.getvalue()
  else:
    def retval():
      return None

  def Print(x):
    out.write(x)

  out.write((barfmt % lexcess) \
            + "".join([ barfmt % i for i in bar ]) \
            + (barfmt % rexcess) \
            + "\n")
  out.write(" " * width + "<" \
            + "".join([ midfmt % i for i in midpt ]) \
            + " " * width + ">" \
            + "\n")

  return retval()

reason: import numpy,from numpy


text:
"""
Transfer Learning exp1
1) groups - only CLL and normal
   model : base - 9F all CLL and normal samples
           target - 5F - increasing sample size (start with very few samples)

2) groups - only rare subtypes ( no CLL, MBl, normal)
   model : base - 9F all rare subtypes samples
           target - 5F - increasing sample size (start with very few samples)
"""

from flowcat import classifier, utils, io_functions
from flowcat.constants import DEFAULT_CLASSIFIER_CONFIG, GROUPS, DEFAULT_CLASSIFIER_ARGS
from flowcat import flowcat_api as fc_api
from flowcat.classifier import som_dataset
from flowcat.classifier.models import create_model_multi_input

import numpy as np
import math

def create_class_weight(labels_dict, mu=0.15):
    total = np.sum(list(labels_dict.values()))
    keys = labels_dict.keys()
    class_weight = dict()

    for key in keys:
        score = math.log(mu*total/float(labels_dict[key]))
        class_weight[key] = score if score > 1.0 else 1.0

    return class_weight

SOM_DATASET = utils.URLPath("/data/flowcat-data/2020-04_merged_train/MLL5F")
OUTPUT = utils.URLPath("/data/flowcat-data/2020-04_merged_train/TL/class_weights/mu_30/model_30")

def create_class_weight(labels_dict, mu=0.30):
    total = np.sum(list(labels_dict.values()))
    keys = labels_dict.keys()
    class_weight = dict()

    for key in keys:
        score = math.log(mu*total/float(labels_dict[key]))
        class_weight[key] = score if score > 1.0 else 1.0

    return class_weight

LOGGER = utils.logs.setup_logging(None, "classify")

groups = GROUPS
tubes = ("1")
#tubes = ("1", "2")

mapping = None
dataset = som_dataset.SOMDataset.from_path(SOM_DATASET)

train_dataset, validate_dataset = fc_api.prepare_classifier_train_dataset(dataset, split_ratio=0.3, groups= groups, mapping= mapping, balance=None)
labels_dict = train_dataset.group_count

config = classifier.SOMClassifierConfig(**{
"tubes": {tube: dataset.config[tube] for tube in tubes},
"groups": groups,
"pad_width": 2,
"mapping": mapping,
"cost_matrix": None,
})

class_weight = create_class_weight(labels_dict)
#class_weight = utils.classification_utils.calculate_group_weights(labels_dict)
class_weight = {
    i: class_weight.get(g, 1.0) for i, g in enumerate(groups)
    }
print(class_weight)
model = fc_api.train_som_classifier(train_dataset, validate_dataset, config, class_weights= class_weight)

model.save(OUTPUT)
model.save_information(OUTPUT)

max_stars_repo_name: NandithaMallesh/mergeTL

max_issues_repo_name: NandithaMallesh/mergeTL

max_forks_repo_name: NandithaMallesh/mergeTL


text:
import cv2
import numpy as np
import copy
import matplotlib.pyplot as plt
from scipy.ndimage import label

def threshold(im, value):
    new_image = np.zeros(im.shape)

    for i in range(im.shape[0]):
        for j in range(im.shape[1]):
                if im[i,j] > value:
                    new_image[i,j] = 0
                else:
                    new_image[i,j] = 1

    return new_image

def histogram(im):
    hist, be = np.histogram(im, 256, (0, 256))
    return hist

def add_neighbor(collection, tup):
    if (tup[0] > 0) and (tup[1] > 0):
        if tup[0] > tup[1]:
            collection.add((tup[1], tup[0]))
        else:
            collection.add(tup)

def identify_neighbors(im_binary):

    colored = copy.deepcopy(im_binary)
    ctr = 1
    neighbors = set()
    for i in range(1, colored.shape[0]):
        for j in range(1, colored.shape[1]-1):
            if colored[i,j] == 1:
                mask = [colored[i-1, j-1],
                        colored[i - 1, j],
                        colored[i - 1, j + 1],
                        colored[i, j - 1]]
                indices = np.nonzero(mask)[0]

                if len(indices) == 0: ## vsude v okoli jsou hodnoty pozadi
                    ctr += 1
                    colored[i,j] = ctr
                    continue
                else: ## v okoli jsou i jine hodnoty
                    no_zeros = []
                    no_zeros[:] = (value for value in mask if value != 0)

                    if np.all(no_zeros == mask[indices[0]]): # vsechny hodnoty v okoli jsou stejne
                        colored[i,j] = mask[indices[0]]
                    else:
                        colored[i,j] = np.min(no_zeros)
                        add_neighbor(neighbors, (colored[i-1, j-1], colored[i-1, j]))
                        add_neighbor(neighbors, (colored[i-1, j-1], colored[i - 1, j + 1]))
                        add_neighbor(neighbors, (colored[i - 1, j + 1], colored[i, j - 1]))

    return colored, neighbors

def join_neighbors(im_colored, neighbors_set):
    im_joined = im_colored
    for i in range(im_joined.shape[0]):
        for j in range(im_joined.shape[1]):
            if im_joined[i,j] > 0:
                matches = []
                matches[:] = (value for value in neighbors_set if value[1] == im_joined[i,j])

                if len(matches) > 0: # nalezli se sousede
                    im_joined[i,j] = min(matches, key = lambda t: t[0])[0]
    return im_joined

## returns format (dictionary) {area_number: [y_t, x_t, sum]}
def area_information(im_colored):
    # {area_number: [y_sum, x_sum, sum]}
    info = dict()
    for i in range(im_colored.shape[0]):
        for j in range(im_colored.shape[1]):
                area_n = im_colored[i,j]
                if area_n > 0:
                    if info.get(area_n, 0) == 0:
                        info[area_n] = [i, j, 1]
                    else:
                        info[area_n] = [sum(x) for x in zip(info[area_n], [i,j,1])]

    for key in info:
        #devide by sum
        info[key][0] = int(np.floor(info[key][0] / info[key][2]))
        info[key][1] = int(np.floor(info[key][1] / info[key][2]))
    return info

def draw_centers(original_image, segments_info, color):

    for key in segments_info:

        x = segments_info[key][1]
        y = segments_info[key][0]

        for m in range(y-2, y+3):
            for n in range (x-2, x+3):
                original_image[m,n] = color

def invert(binary_im):
    for i in range(binary_im.shape[0]):
        for j in range(binary_im.shape[1]):
            if binary_im[i,j] == 0:
                binary_im[i,j] = 1
            else:
                binary_im[i,j] = 0

image = cv2.imread(r'images/cv08_im1.bmp')
im_copy = copy.deepcopy(image)
thresholded = threshold(image[:,:,0],100)
B = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(7,7))
bw1 = cv2.morphologyEx(thresholded, cv2.MORPH_OPEN, B)

colored, neighbors = identify_neighbors(bw1)

## spojeni sousedu (pruchod 2)
joined = join_neighbors(colored, neighbors)

## ziskani informaci (teziste, pocet pixelu, apod)
segments_info = area_information(joined)

## vykresleni stredu do objektu (cervena)
draw_centers(im_copy, segments_info, [255,0,0])

plt.figure()

plt.subplot(2,2,1)
plt.title("segmentation im1")
plt.imshow(bw1, 'gray')

plt.subplot(2,2,2)
plt.title("centers im1")
plt.imshow(im_copy)

image2 = cv2.imread(r'images/cv08_im2.bmp')
im2_copy = copy.deepcopy(image2)
image2 = cv2.cvtColor(image2, cv2.COLOR_BGR2YCR_CB)
im_cr = image2[:,:,1]
im_thresholded = threshold(im_cr, 120)
B2 = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(8,8))
bw2 = cv2.morphologyEx(im_thresholded, cv2.MORPH_CLOSE, B2)

invert(bw2)

colored, neighbors = identify_neighbors(bw2)

## spojeni sousedu (pruchod 2)
joined = join_neighbors(colored, neighbors)

## ziskani informaci (teziste, pocet pixelu, apod)
segments_info = area_information(joined)

## vykresleni stredu do objektu (cervena)
draw_centers(im2_copy, segments_info, [255,0,0])

plt.subplot(2,2,3)
plt.title("segmentation im2")
plt.imshow(bw2, 'gray')

plt.subplot(2,2,4)
plt.title("centers im2")
plt.imshow(im2_copy)

plt.show()

max_stars_repo_name: david-salac/PythonMVSnippet

max_issues_repo_name: david-salac/PythonMVSnippet

max_forks_repo_name: david-salac/PythonMVSnippet


text:
import time
import numpy as np

try:
    from .options import TrainOptions
    from .datasets import create_dataset
    from .models import create_model
    from .util.visualizer import Visualizer
except ImportError:
    from options import TrainOptions
    from datasets import create_dataset
    from models import create_model
    from util.visualizer import Visualizer

if __name__ == '__main__':
    opt = TrainOptions().parse()
    dataset = create_dataset(opt)
    dataset_size = len(dataset)    # get the number of images in the dataset
    print('the number of training images = %d' % dataset_size)

    model = create_model(opt)      # create a model with model parameters in `opt`
    model.setup(opt)               # regular setup: load and print networks; create schedulers
    visualizer = Visualizer(opt)
    total_iters = 0

    for epoch in range(opt.start_epoch, opt.start_epoch + opt.n_epochs + 1):
        epoch_since = time.time()
        iter_data_since = time.time()
        epoch_iter = 0
        visualizer.reset()
        model.update_learning_rate()    # update learning rates in the beginning of every epoch
        for i, batch in enumerate(dataset):
            iter_since = time.time()
            if total_iters % opt.print_freq == 0:
                t_data = iter_since - iter_data_since

            total_iters += opt.batch_size
            epoch_iter += opt.batch_size
            model.set_input(batch)         # unpack data from dataset and apply preprocessing
            model.optimize_parameters()    # calculate loss functions, get gradients, update network weights

            if total_iters % opt.display_freq == 0:
                save_result = False
                model.compute_visuals()
                visuals = model.get_current_visuals()
                visualizer.display_current_results(visuals, epoch, save_result)

            if total_iters % opt.print_freq == 0:
                losses = model.get_current_losses()
                t_comp = (time.time() - iter_since) / opt.batch_size
                visualizer.print_current_losses(epoch, epoch_iter, losses, t_comp, t_data)
                if opt.display_id > 0:
                    visualizer.plot_current_losses(epoch, float(epoch_iter) / dataset_size, losses)

            if total_iters % opt.save_latest_freq == 0:
                print('saving the latest model (epoch %d, total_iters %d)' % (epoch, total_iters))
                save_suffix = 'iter_%d' % total_iters if opt.save_by_iter else 'latest'
                model.save_networks(save_suffix)

            iter_data_time = time.time()

        if epoch % opt.save_epoch_freq == 0:
            print('saving the model at the end of epoch %d, iters %d' % (epoch, total_iters))
            model.save_networks('latest')
            model.save_networks(epoch)

    print('end of epoch %d / %d \t Time Taken: %d sec' % (epoch, opt.start_epoch + opt.n_epochs, time.time() - epoch_since))

max_stars_repo_name: atomicoo/EnhanceIMG