text:
# -*- encoding: utf-8 -*-

from cogdata.utils.register import register
import os
import PIL
import torch
import json
from torch.utils.data import DataLoader

from .base_task import BaseTask
from cogdata.utils.logger import set_logger, get_logger
import numpy as np

from icetk import icetk as tokenizer

def txt_collate_fn(x):
    return x

@register
class IcetkTextTask(BaseTask):
    '''handle tokenization
    '''

    def __init__(self, saver, **kwargs) -> None:
        """config saver
        """
        self.saver = saver

    def get_transform_fn(self, transform=None, **kwargs):

        return transform

    def process(self, sub_datasets, progress_record=None, dataset_dir='', **kwargs):
        """ Use cuda to process batch data from dataloader,
            save via Saver,
            report progress every 1/5000 ?
            final commit saver

        Parameters
        ----------
        sub_datasets:[Dataset]
            All datasets in processing list
        progress_record:ProgressBar
            The progress bar for this task
        dataset_dir:str
            The path of the dataset folder

        Returns
        -------
        int
            0 - Process successfully
        """
        text_format = kwargs.get('text_format', 'raw_txt')
        batch_size = kwargs.get('batch_size', 128)
        num_workers = kwargs.get('dataloader_num_workers', 2)
        ratio = kwargs.get('ratio', 1)

        cnt, total_cnt = 0, sum([len(dataset) for dataset in sub_datasets])
        for dataset in sub_datasets:
            loader = DataLoader(dataset, batch_size=batch_size, shuffle=False,
                                num_workers=num_workers, collate_fn=txt_collate_fn, pin_memory=True)
            for batch_txts in loader:

                if dataset.use_bytes_as_length:
                    cnt += sum(len(x) for x in batch_txts)
                else:
                    cnt += len(batch_txts)  # may not full batch
                if cnt > total_cnt * ratio:
                    break
                ret = []
                for txt in batch_txts:
                    if len(txt) == 0:
                        continue
                    if text_format == 'raw_txt':
                        ret.extend(tokenizer.encode(txt))
                        ret.append(tokenizer['</s>']) # </s> TODO
                    else:
                        tokenized_line = process_special_line(txt, text_format)
                        if tokenized_line is not None:
                            ret.extend(tokenized_line)
                            ret.append(tokenizer['</s>']) # </s> TODO

                if len(ret) == 0:
                    continue
                data = torch.tensor(ret)

                self.saver.save(data)

                if dataset.use_bytes_as_length and (cnt // (total_cnt // 1000)) % 5 == 0:
                    if progress_record is not None:
                        progress_record.update(cnt, total_cnt)
                elif cnt // batch_size % 50 == 0:
                    # get_logger().info("progress {}/{}".format(cnt, total_cnt))
                    if progress_record is not None:
                        progress_record.update(cnt, total_cnt)
            self.saver.commit()
        return 0

def process_special_line(data, txt_format,
                         pile_valid_set=['Pile-CC', "OpenWebText2", "Wikipedia (en)"]):
    if txt_format == 'jsonl': # default ["text"]
        data = json.loads(data)
        text = data.get("text", None)
        return tokenizer.encode(text)
    elif txt_format == 'ccnews':
        data = json.loads(data)
        text = ""
        title = data.get("title", None)
        description = data.get("description", None)
        maintext = data.get("maintext", None)
        if title:
            text += title.strip() + " "
        if description and (not maintext or not maintext.startswith(description)):
            text += description.strip() + " "
        if maintext:
            text += maintext
        if len(text) > 100:
            return tokenizer.encode(text)
        else:
            return None
    elif txt_format == 'baike':
        data = json.loads(data)
        text = data.get("title", "") + data.get("abstract", "") + data.get("content", "")
        if text:
            return tokenizer.encode(text)
        else:
            return None
    elif txt_format == 'zhihu':
        data = json.loads(data)
        ans_length = len(data.get("ans-content", ""))
        ans_up = data.get("ans-up-num", "")
        ans_up = int(ans_up) if ans_up else 0
        if ans_length > 100 or ans_up > 1000:
            qtitle = data["q_title"]
            qcontent = data["q-content"]
            if qcontent is None:
                qcontent = ""
            else:
                qcontent = "问题描述：" + qcontent
            # user = data.get("user-signature", "")
            # if len(user) > 0:
            #     user = "回答用户：" + user
            text = "问题：" + qtitle + qcontent + "回答：" + data["ans-content"]
            return tokenizer.encode(text)
        else:
            return None
    elif txt_format == 'zhidao':
        data = json.loads(data)
        if "title" not in data:
            return None
        qtitle = data["title"]
        qcontent = data.get("content", "")
        if len(qcontent) > 0:
            qcontent = "问题描述：" + qcontent
        qtitle = "问题：" + qtitle
        text_all = qtitle +  qcontent + "回答："
        if "best_answer" in data:
            text = data["best_answer"]["content"]
            if len(text) > 10:
                text_all += text
        for answer in data.get("other_answers", []):
            text = answer["content"]
            if len(text) > 100:
                text_all += text
        return tokenizer.encode(text_all)
    elif txt_format == 'pile':
        data = json.loads(data)
        source = data["meta"].get("pile_set_name", None)
        if source not in pile_valid_set:
            return None
        text = data.get('text', '')
        return tokenizer.encode(text)
    else:
        raise ValueError(f'Unknown text format: \"{txt_format}\"')

max_stars_repo_name: Sleepychord/cogdata

max_issues_repo_name: Sleepychord/cogdata

max_forks_repo_name: Sleepychord/cogdata


text:
import numpy as np
import matplotlib.pyplot as plt
import os
from utilities import load_pickle

def draw_curve(a, save_folder, patient_name):
    plt.figure()
    plt.subplot(311)
    plt.plot(a["training_loss"], label='train_loss',color='blue')
    plt.plot(a["validation_loss"], label='val_loss',color='orange')
    plt.title("loss")
    plt.legend()
    plt.subplot(312)
    plt.plot(a["training_acc"], label='train_acc',color='blue')
    plt.plot(a["validation_acc"], label='vali_acc',color='orange')
    plt.plot(a["test_acc"], label='test_acc', color='red')
    plt.title("accuracy")
    plt.ylim([0.5, 1.0])
    plt.legend()
    plt.subplot(313)
    #plt.plot(a["training_f1"], label='train_f1')
    plt.plot(a["validation_f1"], label='vali_f1',color='orange')
    plt.plot(a["test_f1"], label='test_f1', color='red')
    plt.title("F1-score")
    plt.legend()
    plt.suptitle(patient_name)
    plt.savefig(os.path.join(save_folder, patient_name+".jpg"))
    plt.close()
if __name__ == "__main__":

    '''
    This file is to plot the training curve,

    path: in result of training folder,  there is a file called  training_curve_a.pkl and training_curve_s.pkl
          which contains train stats.
    out_folder : the output folder to dump the plot
    '''

    path = ""
    out_folder = ""
    for p_name in sorted(os.listdir(path)):
        stats = load_pickle(os.path.join(path, p_name,"training_curve_a.pkl"))
        draw_curve(stats,out_folder,p_name+"_artifacts")
        stats = load_pickle(os.path.join(path, p_name,"training_curve_s.pkl"))
        draw_curve(stats,out_folder,p_name+"_spike")

ext: py

max_stars_repo_name: roychowdhuryresearch/HFO-Classification

max_issues_repo_name: roychowdhuryresearch/HFO-Classification

max_forks_repo_name: roychowdhuryresearch/HFO-Classification

reason: import numpy


text:
import numpy as np
from q_learn_linear_func_approx import execute

config = {
  'MAX_EPISODES': 2500,
  # 'MAX_STEPS_PER_EPISODE': 1000,
  'LEARNING_RATE': 0.0005,
  # 'DISCOUNT_RATE': 0.99,
  'EXPLORATION_DECAY': 0.001,
  'MIN_EXPLORATION_RATE': 0.01,
  'DO_RENDER': True,
  # 'SHOW_STATS_PLOT': True,
}

def reward_function_cartpole(env_reward, state, next_state, step):
  cart_position, cart_speed, pole_angle, pole_speed = [np.abs(s) for s in state]
  next_cart_position, next_cart_speed, next_pole_angle, next_pole_speed = [np.abs(s) for s in next_state]
  decrease_in_angle = pole_angle - next_pole_angle
  decrease_in_pole_speed = pole_speed - next_pole_speed
  reward = 1 + 100 * decrease_in_angle - next_pole_speed  # - pole_angle #  - (next_pole_speed - pole_speed) + (cart_speed - next_cart_speed)
  return reward # env_reward - next_pole_speed - next_pole_angle

def run():
  rewards = execute('CartPole-v1', reward_function_cartpole, **config)

if __name__ == '__main__':
  run()

max_stars_repo_name: johanlahti/reinforcement-algorithms

max_issues_repo_name: johanlahti/reinforcement-algorithms

max_forks_repo_name: johanlahti/reinforcement-algorithms


text:
import numpy as np

def soda_strategy_discount(n_energy, n_nosugar):
    if n_energy <= n_nosugar:
        discount = -0.2;
    else:
        discount = 0.2;

    return discount

def soda_strategy_nodiscount(n_energy, n_nosugar):
    discount = 0.0;
    return discount

def soda_strategy_param(n_energy, n_nosugar, T, V):
    assert len(T) == len(V)
    diff = abs(n_energy - n_nosugar)
    sg = np.sign(n_energy - n_nosugar)
    if diff < T[0]:
        return 0
    if T[0] < diff <= T[1]:
        return sg * V[0]
    if T[1] < diff <= T[2]:
        return sg * V[1]
    if T[2] < diff:
        return sg * V[2]
    return 0

max_stars_repo_name: GeoffNN/BanditsGameTheoryRL

max_issues_repo_name: GeoffNN/BanditsGameTheoryRL

max_forks_repo_name: GeoffNN/BanditsGameTheoryRL


text:
import GPy
import GPyOpt
import numpy as np
import copy
from gym.spaces import *

from GPy.kern import Matern52, RBF
from GPyOpt.models.gpmodel import GPModel, GPModel_MCMC
from GPyOpt.optimization.acquisition_optimizer import ContextManager, AcquisitionOptimizer
from GPyOpt.methods.bayesian_optimization import BayesianOptimization

__all__ = ['CGPUCBPolicy', 'dict2ker', 'space2domain',
           'DEFAULT_BO_CONFIG', 'DEFAULT_SURROGATE_CONFIG']

DEFAULT_SURROGATE_CONFIG = {
    # Kernel over observations
    'ob_ker': {'type': RBF,
               'kwargs': {
                     'lengthscale': {'v': 1.,
                                     'prior': GPy.priors.Gamma,
                                     'prior_params': {'E': 1.,
                                                      'V': 0.5}},
                     'variance': {'v': 1.,
                                     'prior': GPy.priors.Gamma,
                                     'prior_params': {'E': 1.,
                                                      'V': 0.5}}
                 }},
    # Kernel over acitons
    'ac_ker': {'type': RBF,
                 'kwargs': {
                     'lengthscale': {'v': 1.,
                                     'prior': GPy.priors.Gamma,
                                     'prior_params': {'E': 1.,
                                                      'V': 0.5}},
                     'variance': {'v': 1.,
                                     'prior': GPy.priors.Gamma,
                                     'prior_params': {'E': 1.,
                                                      'V': 0.5}}
                 }},
    # Likelihood
    'lik': {'type': GPy.likelihoods.Gaussian,
                         'kwargs': {
                             'variance': {'v': 0.0001,
                                          'prior': GPy.priors.Gamma,
                                          'prior_params': {'E': 0.01,
                                                           'V': 0.01}}
                         }},
    # How to compose the observation kernel and the action one
    'composition': 'prod',

    # Wheter to use a sparse GP model (does not work with MCMC)
    'sparse': False,
    'num_inducing': 10,

    # Whether to use MCMC to marginalize the hyperparameters
    'MCMC': False,

    # Number of restarts or iterations per restart for MLE estimate of
    # hyperparamers if MCMC is False
    'optimize_restarts': 5,
    'max_iters': 0,

    # Parameters for Markov chain and Monte Carlo estimate if MCMC True
    'n_samples': 10,
    'subsample_interval': 10,
    'n_burnin': 100,
    'step_size': 0.1
}

DEFAULT_BO_CONFIG = {
    'beta': lambda t: 2,
    'maximize': True,
    'model_update_interval': 1,
    'normalize_Y': False
}

class CGPUCBPolicy(object):
    """
    Class for Contextual GPUCB policy.

    This class implements a myopic policy using contextual GPUCB. In this
    case, the state of the system is used as a context and the actions as
    arms in a multi-armed bandit problem. It can be effective when long term
    planning is not relevant to solve the task.

    ob_space: gym.spaces
        Observation space of the MDP
    ac_space: gym.spaces
        Action space of the MDP
    new_surrogate_config: dict
        Dictionary to update the default config of the surrogate model. See
        src.contextual_bandits.policies for the default.
    new_BO_config: dict
        Dictionary to update the default config of the BO problem. See
        src.contextual_bandits.policies for the default.
    """
    def __init__(self, ob_space, ac_space, new_surrogate_config={},
                 new_BO_config={}, n_env=1, n_steps=1, n_batch=1):
        self.ob_space = ob_space
        self.ac_space = ac_space

        # Define BO domain from gym spaces
        ob_dom = space2domain(ob_space, 'obs')
        ac_dom = space2domain(ac_space, 'act')
        self.domain = ob_dom + ac_dom

        # Update the configs
        gp = copy.deepcopy(DEFAULT_SURROGATE_CONFIG)
        gp.update(new_surrogate_config)
        BO_config = copy.deepcopy(DEFAULT_BO_CONFIG)
        BO_config.update(new_BO_config)

        # Initialize kernels and GP model
        gp['ob_ker']['kwargs']['input_dim'] = {'v': len(ob_dom)}
        gp['ob_ker']['kwargs']['active_dims'] = {'v': np.arange(len(ob_dom))}
        gp['ac_ker']['kwargs']['input_dim'] = {'v': len(ac_dom)}
        gp['ac_ker']['kwargs']['active_dims'] = {'v': np.arange(
            len(ob_dom), len(ob_dom) + len(ac_dom))}

        ob_ker = dict2ker(gp['ob_ker'])
        ac_ker = dict2ker(gp['ac_ker'])
        lik = dict2ker(gp['lik'])

        if gp['composition'] == 'prod':
            ker = ob_ker * ac_ker
        elif gp['composition'] == 'sum':
            ker = ob_ker + ac_ker
        else:
            raise ValueError('Only product and sum composition implemented')

        if gp['MCMC']:
            self.model = CustomGPModel_MCMC(
                kernel=ker, lik=lik, verbose=False, n_samples=gp['n_samples'],
                subsample_interval=gp['subsample_interval'],
                n_burnin=gp['n_burnin'], step_size=gp['step_size'])
            self.acquisition_type = 'LCB_MCMC'
            self.sparse = False
        else:
            self.model = CustomGPModel(
                kernel=ker, lik=lik, sparse=gp['sparse'], num_inducing=gp[
                'num_inducing'], optimize_restarts=gp['optimize_restarts'],
                max_iters=gp['max_iters'], verbose=False)
            self.sparse = gp['sparse']
            self.acquisition_type = 'LCB'

        # Bayesian optimization
        self.BO = None
        self.beta = BO_config['beta']
        self.maximize = BO_config['maximize']
        self.model_update_interval = BO_config['model_update_interval']
        self.normalize_Y = BO_config['normalize_Y']

        # Optimization problem when we are performing pure exploitation
        self.exploit_func = None
        self.exploit_optimizer = None

    def step(self, obs, explore=True):
        """
        Returns the policy for a single step.

        Parameters
        ----------
        obs:
            Observation from env.
        explore: bool
             If True, it suggests also exploratory actions.
        """
        # We need at least one data point before creating the BO model.
        # Until then, we sample from the action space
        if self.BO is None:
            return self.ac_space.sample()
        else:
            # Add new data to gp and optimize hyperparameters if max_iters > 0
            self.BO._update_model(self.BO.normalization_type)

            # Create context dict
            if not hasattr(obs, '__iter__'):
                obs = np.array([obs])
            obs = np.asarray(obs)
            obs_copy = obs.copy()
            context = {'obsvar_{}'.format(i): v for i, v in
                       enumerate(obs_copy.flatten())}

            # Maximizing the UCB criterion trades-off exploration/exploitation
            if explore:
                self.BO.num_acquisitions += 1
                self.BO.acquisition.exploration_weight = self.beta(self.BO.num_acquisitions)
                self.BO.context = context
                x = self.BO._compute_next_evaluations()

            # In case we don't explore, we maximize the posterior mean
            else:
                # Optimize the exploit function (posterior mean)
                self.exploit_optimizer.context_manager = ContextManager(
                    space=self.BO.space, context=context)
                x, fx = self.exploit_optimizer.optimize(f=self.exploit_func)

            # Remove context (i.e. state) and cast to correct type
            x = np.squeeze(x)
            x_no_context = x[len(obs):]
            if isinstance(self.ac_space, Discrete):
                x_no_context = int(x_no_context[0])
            elif isinstance(self.ac_space, MultiDiscrete):
                x_no_context.astype(np.int64)

            return x_no_context

    def update(self, obs, act, r):
        """
        Add data point(s) and, if necessary, create the BO model.

        Parameters
        ---------
        obs: np.array
            Observation (n_data x obs_space_dim)
        act: np.array
            Action (n_data x obs_action_dim)
        r: float or np.array
            Reward
        """
        newX = np.atleast_2d(np.hstack((obs, act)))
        newY = np.atleast_2d(r)

        if self.BO is None:
            self.set_XY(newX, newY)  # Automatically creates BO
        else:
            oldX, oldY = self.get_XY()
            self.set_XY(np.vstack((oldX, newX)),
                        np.vstack((oldY, newY)))

    def set_XY(self, X, Y):
        """
        Set the values of the data points in the BO model.

        Since BO minimizes by default, we store -Y in case the optimize flag is
        set to True.

        Parameters
        ---------
        X: np.ndarray
            Array of inputs (n_data x input_dim)
        X: np.ndarray
            Array of observations (n_data x 1)
        """
        if X is not None and Y is not None:
            if self.maximize:
                Y = np.copy(Y)
                Y *= -1
            if self.BO is None:
                self._create_BO(X, Y)
            else:
                self.BO.X = X
                self.BO.Y = Y
        else:
            print('Skipping set XY step as the values provided are None')

    def _create_BO(self, X, Y):
        """
        Create the internal BO object with the values provided in the config.

        Together with the BO object, we create the exploit function and the
        exploit function optimizer that depend on it and are used to provide
        purely exploitative actions.

        Parameters
        ---------
        X: np.ndarray
            Array of inputs (n_data x input_dim)
        X: np.ndarray
            Array of observations (n_data x 1)
        """
        # Create BO
        self.BO = BayesianOptimization(f=None,
                                       domain=self.domain,
                                       model=self.model,
                                       acquisition_type=self.acquisition_type,
                                       model_update_interval=self.model_update_interval,
                                       normalize_Y=self.normalize_Y,
                                       X=X,
                                       Y=Y)
        # Create the internal GP model
        self.BO._update_model(self.BO.normalization_type)
        # The posterior mean is the function we optimize when we want to
        # exploit
        def exploit_f(x):
            x = np.atleast_2d(x)
            m, _ = self.predict(x, with_noise=False)
            if self.maximize:
                return -m
            else:
                return m

        self.exploit_func = exploit_f
        self.exploit_optimizer = AcquisitionOptimizer(self.BO.space, 'lbfgs')

    def predict(self, X, with_noise=False):
        """
        Predict mean and standard deviation of the reward at specified inputs.

        Parameters
        ----------
        X: np.ndarray
            Array of inputs where to make predictions (n_data x input_dim)
        with_noise: bool
            If true, predicts the function value + noise
        """
        if self.BO is not None:
            self.BO._update_model(self.BO.normalization_type)
            mu, sigma = self.model.predict(X, with_noise=with_noise)
            if self.maximize:
                mu *= -1
            return mu, sigma
        else:
            raise NotImplementedError('Not implemented the predict function '
                                      'when the BO model and its GP are not '
                                      'istanciated')

    def get_XY(self):
        """
        Get the input and observations.
        """
        if self.BO is not None:
            X = np.copy(self.BO.X)
            # We flip the sign in maximization case as BO stores -Y internally
            Y = -np.copy(self.BO.Y) if self.maximize else np.copy(self.BO.Y)
        else:
            X = Y = None
        return X, Y

class CustomGPModel(GPModel):
    """
    GP model that works with BO class with custom kern/lik.
    """
    def __init__(self, kernel, lik, num_inducing=10, sparse=False, **kwargs):
        super(CustomGPModel, self).__init__(**kwargs)
        # Overwrite kernel initialization of GPModel
        self.kernel = kernel
        self.lik = lik
        self.num_inducing = num_inducing
        self.kwargs = kwargs
        self.sparse = sparse

    def _create_model(self, X, Y):
        if not self.sparse:
            self.model = GPy.core.GP(X, Y, kernel=self.kernel,
                                     likelihood=self.lik)
        else:
            # Sample points from unit [0, 1]^d to reach the desired number
            # of points
            Z = np.random.rand(self.num_inducing, X.shape[1])

            self.model = GPy.core.SparseGP(X, Y, Z, kernel=self.kernel,
                                           likelihood=self.lik)

    def copy(self):
        copied_model = CustomGPModel(kernel=self.kernel,
                                     lik=self.lik,
                                     **self.kwargs)
        if self.model is not None:
            copied_model._create_model(self.model.X, self.model.Y)
            copied_model.updateModel(self.model.X, self.model.Y, None, None)
        return copied_model

class CustomGPModel_MCMC(GPModel_MCMC):
    """
    GP_MCMC model that works with BO class with custom kern/lik.
    """
    def __init__(self, kernel, lik, **kwargs):
        super(CustomGPModel_MCMC, self).__init__(**kwargs)
        # Overwrite kernel initialization of GPModel_MCMC
        self.kernel = kernel
        self.lik = lik
        self.kwargs = kwargs

    def _create_model(self, X, Y):
        self.model = GPy.core.GP(X, Y, kernel=self.kernel,
                                     likelihood=self.lik)

    def copy(self):
        copied_model = CustomGPModel_MCMC(kernel=self.kernel,
                                          lik=self.lik,
                                          **self.kwargs)
        if self.model is not None:
            copied_model._create_model(self.model.X, self.model.Y)
            copied_model.updateModel(self.model.X, self.model.Y, None, None)
        return copied_model

def space2domain(space, base_name=''):
    """Helper function to create GPyOpt domain from gym.space"""
    domain = []

    if isinstance(space, Box):
        if space.dtype in [np.int64, np.int32, np.int16]:
            vtype = 'discrete'
        elif space.dtype in [np.float64, np.float32, np.float16]:
            vtype = 'continuous'
        else:
            raise ValueError('Not supported dtype for conversion: {'
                             '}'.format(space.dtype))

        # GPyOpt cannot handle unbounded domains
        low = np.clip(space.low, -1e5, 1e5).flatten()
        high = np.clip(space.high, -1e5, 1e5).flatten()
        for i, (l, u) in enumerate(zip(low, high)):
            domain.append(
                {'name': '{}var_{}'.format(base_name, i),
                 'type': vtype,
                 'domain': (l, u),
                 'dimensionality': 1})
    elif isinstance(space, Discrete):
        domain.append({'name': '{}var_{}'.format(base_name, 0),
                       'type': 'discrete',
                       'domain': np.arange(space.n),
                       'dimensionality': 1})
    elif isinstance(space, MultiDiscrete):
        for i, n in enumerate(space.nvec):
            domain.append({'name': '{}var_{}'.format(base_name, i),
                           'type': 'discrete',
                           'domain': np.arange(n),
                           'dimensionality': 1})
    else:
        raise NotImplementedError('gym space to GPyOpt domain conversion '
                                  'implemented only for Box, Discrete and '
                                  'MultiDiscrete')
    return domain

def dict2ker(d):
    """
    Helper function to create kernel/likelihood from dictionary.

    Example dict structure
        d = {'type': RBF,
             'kwargs': {
                        'input_dim': {'v': 1},
                        'lengthscale': {'v': 0.1, 'prior': LogGaussian,
                                        'prior_params': {'mu': 0.1,
                                                         'sigma': 0.2},
                                        'constraint': 'bounded',
                                        'constraint_params': {'lower':1,
                                                             'upper':2}}
                        }}
    Should return
    k = GPy.kern.RBF(input_dim=1, lengthscale=0.1)
    k.lengthscale.set_prior(GPy.priors.LogGaussian(0.1, 0.2))
    k.lengthscale.constrain_bounded(1, 2)
    """
    # Initialize kernel
    tmp_dict = {param_name: param['v'] for param_name, param in d[
        'kwargs'].items()}
    k = d['type'](**tmp_dict)

    # Set priors and constraints
    for param_name, param in d['kwargs'].items():
        prior = param.get('prior')
        constraint = param.get('constraint')
        if prior is not None:
            prior_params = param.get('prior_params', {})
            if prior == GPy.priors.Gamma:
                k[param_name].set_prior(prior.from_EV(**prior_params))
            else:
                k[param_name].set_prior(prior(**prior_params))
        if constraint is not None:
            constraint_params = param.get('constraint_params')
            if constraint == 'bounded':
                k[param_name].constrain_bounded(**constraint_params)
            elif constraint == 'fixed':
                k[param_name].constrain_fixed(**constraint_params)
            elif constraint == 'negative':
                k[param_name].constrain_negative()
            elif constraint == 'positive':
                k[param_name].constrain_positive()
    return k


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

def get_data():
    exc = input("Enter name of exchange the data was pulled from: ")
    asset = input("Enter name of the asset you wish to work with: ")
    return({"exc":exc,"asset":asset})

def gen():
    excas = get_data()
    total = [["STATS","Mean","Variance","Volatility"]]
    ohlc = ["open","high","low","close","volume"]
    try:
        exas = pd.read_csv(excas["exc"] + excas["asset"] + "1.csv")
    except:
        print("File Not Found")
        get_data()
    for i in range(len(ohlc)):
        dub = []
        dub.append(excas["asset"]+ohlc[i])
        dub.append(exas[excas["asset"]+ohlc[i]].mean())
        dub.append(exas[excas["asset"]+ohlc[i]].var())
        dub.append(exas[excas["asset"]+ohlc[i]].std())
        total.append(dub)
    total = pd.DataFrame(total).transpose()
    file = (excas["exc"] + "stats.csv")
    path = os.getcwd() + "\\" + file
    if(os.path.isfile(path)):
        print("Existing file")
        with open(file, 'a') as f:
            full.to_csv(f,index=False,header=False)
    else:
        with open(file, 'a') as f:
            full.to_csv(f,index=False,header=ohlc)
    print("")
    print("Success!")
    print("")

def corr():
    exc = input("Enter Exchange to work with: ")
    asA = input("Enter First Asset/Pair:  ")
    asB = input("Enter Second Asset/Pair:  ")
    full = []
    ohlc = ["","open","high","low","close","volume"]
    try:
        csvA = pd.read_csv(exc + asA + "1.csv")
        csvB = pd.read_csv(exc + asB + "1.csv")
    except:
        print("File Not Found")
        exc = input("Enter Exchange to work with: ")
        asA = input("Enter First Asset/Pair:  ")
        asB = input("Enter Second Asset/Pair:  ")
    full.append(asA + "/" + asB)
    for i in range(len(ohlc)):
        try:
            full.append(csvA[asA + ohlc[i]].corr(csvB[asB+ohlc[i]]))
        except:
            print(ohlc[i])
    full = pd.DataFrame(full).transpose()
    file = (exc + "corr.csv")
    path = os.getcwd() + "\\" + file
    if(os.path.isfile(path)):
        print("Existing file")
        with open(file, 'a') as f:
            full.to_csv(f,index=False,header=False)
    else:
        with open(file, 'a') as f:
            full.to_csv(f,index=False,header=ohlc)
    print("")
    print("Success!")
    print("")

def rollvol():
    asA = get_data()
    full = []
    ohlc = ["","open","high","low","close","volume"]
    try:
        csvA = pd.read_csv(asA["exc"] + asA["asset"] + "1.csv")
    except:
        print("File Not Found")
        get_data()
    print(" N / A temporary")
    '''for i in range(len(ohlc)):
        dub = []
        dub.append(pd.rolling_std('''

def rollhr():
    print("Coming Soon....")

def hrat():
    print("Coming Soon....")

max_stars_repo_name: nlysov/old-trading-project

max_issues_repo_name: nlysov/old-trading-project

max_forks_repo_name: nlysov/old-trading-project


text:
#%% Load libraries.
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns

import pandas as pd
import numpy as np

#%% Load datasets from CSVs.
facebook_data = pd.read_csv('data/facebook.csv')
microsoft_data = pd.read_csv('data/microsoft.csv')

facebook_data = facebook_data.set_index('Date')
microsoft_data = microsoft_data.set_index('Date')

#%% Describe loaded datasets.
display(facebook_data.head(5))
display(microsoft_data.head(5))

display(facebook_data.index.min(), facebook_data.index.max())
display(microsoft_data.index.min(), microsoft_data.index.max())

display(facebook_data.shape)
display(microsoft_data.shape)

display(facebook_data.columns.values)
display(microsoft_data.columns.values)

display(facebook_data.describe())
display(microsoft_data.describe())

#%% Visualize loaded datasets.
fb_2017_2018 = facebook_data.loc['2017-01-01':'2018-12-31'].drop('Volume', axis=1)
ms_2016 = microsoft_data.loc['2016-01-01':'2016-12-31'].drop('Volume', axis=1)

fig, [axis_0, axis_1] = plt.subplots(1, 2, figsize=(15, 5))
axis_0.set_title('Facebook 2017-2018')
fb_2017_2018.plot(ax=axis_0)
axis_1.set_title('Microsoft 2016')
ms_2016.plot(ax=axis_1)
plt.show()

fb_2017_closeyesterday = facebook_data.loc['2017-01-01':'2017-12-31']
fb_2017_closeyesterday = fb_2017_closeyesterday.drop('Volume', axis=1)
fb_2017_closeyesterday['close_yesterday'] = fb_2017_closeyesterday['Close'].shift(1)
fb_2017_closeyesterday.plot()
plt.show()

#%% Add new features to dataset
def fin_df_add_features(df):
    """Add new important features to financial dataset.
    Adding is done inplace.
    """

    # Close price for tomorrow.
    df['close_tomorrow'] = df['Close'].shift(-1)
    # Next day Close price difference.
    df['diff'] = df['close_tomorrow'] - df['Close']
    # Daily return.
    df['return'] = df['diff'] / df['Close']
    # Direction of the Close price.
    df['direction'] = [
        1 if df.loc[el_idx, 'diff'] > 0
        else 0
        for el_idx in df.index
    ]
    # Moving average over 3 days.
    # moving_average_days_cnt = 3
    # df['mov_average_3'] = (df['Close'] + df['Close'].shift(1) + df['Close'].shift(2)) / moving_average_days_cnt
    df['ma_10'] = df['Close'].rolling(10).mean()
    df['ma_40'] = df['Close'].rolling(40).mean()
    df['ma_50'] = df['Close'].rolling(50).mean()
    df['ma_200'] = df['Close'].rolling(200).mean()

    # todo: dropna for the dataframe?

def plot_custom_close_ma_40_200(df, title=""):
    df['Close'].plot(label='Close price')
    df['ma_10'].plot(label='Moving average 10 days')
    df['ma_40'].plot(label='Moving average 40 days')
    df['ma_50'].plot(label='Moving average 50 days')
    df['ma_200'].plot(label='Moving average 200 days')
    plt.title(title)
    plt.legend()
    plt.show()

facebook_data_ext = facebook_data.copy()
fin_df_add_features(facebook_data_ext)
display(facebook_data_ext.head(3))

microsoft_data_ext = microsoft_data.copy()
fin_df_add_features(microsoft_data_ext)
display(microsoft_data_ext.head(3))

# Plot moving averages
# ToNote: both moving averages SMOOTH the original Close price.
# ma_40 - FAST SIGNAL (MORE CLOSELY associated with Close price) (reflects price over the SHORT history)
# ma_200 - SLOW SIGNAL (reflects price over the LONG history)
plot_custom_close_ma_40_200(facebook_data_ext, 'facebook')
plot_custom_close_ma_40_200(microsoft_data_ext, 'microsoft')

# If ma40 > ma200 -> some traders (trend-following traders) believe the stock price will move upwards for a while.

#%% Build a simple trading strategy.
# Strategy: MA10 > MA50 -> buy and hold one share of stock (long 1 share of stock).
def fin_df_add_simplest_strategy_features(df):
    """Add important features for our simplest strategy"""

    # Denote, whether we long or not.
    df['shares'] = [
        1 if df.loc[el_idx, 'ma_10'] > df.loc[el_idx, 'ma_50']
        else 0
        for el_idx in df.index
    ]
    # Daily profit for our simplest strategy
    df['profit'] = [
        df.loc[el_idx, 'diff'] if df.loc[el_idx, 'shares'] == 1
        else 0
        for el_idx in df.index
    ]
    # Cumulative wealth
    df['wealth'] = df['profit'].cumsum()

facebook_data_trading = facebook_data_ext.drop('Volume', axis=1).copy()
microsoft_data_trading = microsoft_data_ext.drop('Volume', axis=1).copy()

fin_df_add_simplest_strategy_features(facebook_data_trading)
fin_df_add_simplest_strategy_features(microsoft_data_trading)

# Plot all the calculated data
# Note: does not describe 'Volume' values.
fig, [axis_0, axis_1] = plt.subplots(1, 2, figsize=(20, 5))
axis_0.set_title('facebook')
facebook_data_trading.plot(ax=axis_0)
axis_1.set_title('microsoft')
microsoft_data_trading.plot(ax=axis_1)
plt.show()

# NOTE: sometimes we tend to be bankrupt!

# Plot the profit data
fig, [axis_0, axis_1] = plt.subplots(1, 2, figsize=(20, 5))
axis_0.set_title('facebook')
facebook_data_trading['profit'].plot(ax=axis_0)
axis_1.set_title('microsoft')
microsoft_data_trading['profit'].plot(ax=axis_1)
plt.show()

# Plot cumulative wealth
fig, [axis_0, axis_1] = plt.subplots(1, 2, figsize=(20, 5))
axis_0.set_title('facebook')
facebook_data_trading['wealth'].plot(ax=axis_0)
axis_1.set_title('microsoft')
microsoft_data_trading['wealth'].plot(ax=axis_1)
plt.show()

#%% Print results if we would really buy & sell daily
def print_money_earned_spent(df, name):
    print(name)
    print('Total money you win: {0}'.format(df.loc[df.index[-2], 'wealth']))
    print('Total money you spent: {0}'.format(df.loc[df.index[0], 'Close']))
    print('Final close price: {0}\n'.format(df.loc[df.index[-1], 'Close']))

print_money_earned_spent(facebook_data_trading, 'facebook')
print_money_earned_spent(microsoft_data_trading, 'microsoft')

#%% [markdown]
# # Questions
# Can we find a better signal for trading? <br>
# How do you evaluate your performance of shared strategy correctly?

max_stars_repo_name: honchardev/Fun

max_issues_repo_name: honchardev/Fun

max_forks_repo_name: honchardev/Fun


text:
"""
Run the uncoverml pipeline for clustering, supervised learning and prediction.

.. program-output:: uncoverml --help
"""

import logging
import os

import joblib
import resource
import json
from os.path import isfile, splitext, exists
from pathlib import Path
import warnings

import click
import numpy as np
import matplotlib
matplotlib.use('Agg')

import uncoverml as ls
import uncoverml.cluster
import uncoverml.config
import uncoverml.features
import uncoverml.geoio
import uncoverml.learn
import uncoverml.mllog
import uncoverml.mpiops
import uncoverml.predict
import uncoverml.validate
import uncoverml.targets
from uncoverml.transforms.linear import WhitenTransform
from uncoverml.transforms import StandardiseTransform
from uncoverml import optimisation, hyopt
# from uncoverml.mllog import warn_with_traceback

log = logging.getLogger(__name__)
# warnings.showwarning = warn_with_traceback
warnings.filterwarnings(action='ignore', category=FutureWarning)
warnings.filterwarnings(action='ignore', category=DeprecationWarning)
warnings.filterwarnings(action='ignore', category=Warning)

@click.group()
@click.option('-v', '--verbosity',
              type=click.Choice(['DEBUG', 'INFO', 'WARNING', 'ERROR']),
              default='INFO', help='Level of logging')
def cli(verbosity):
    ls.mllog.configure(verbosity)

def run_crossval(x_all, targets_all, config):
    crossval_results = ls.validate.local_crossval(x_all, targets_all, config)
    ls.mpiops.run_once(ls.geoio.export_crossval, crossval_results, config)

@cli.command()
@click.argument('pipeline_file')
@click.option('-j', '--param_json', type=click.Path(exists=True), multiple=True,
              help='algorithm parameters json, possibly from a previous optimise job')
@click.option('-p', '--partitions', type=int, default=1,
              help='divide each node\'s data into this many partitions')
def learn(pipeline_file, param_json, partitions):
    config = ls.config.Config(pipeline_file)
    if len(param_json) > 0:
        # log.info('the following jsons are used as params \n', click.echo('\n'.join(param_json)))
        param_dicts = []  # list of dicts
        for pj in param_json:
            with open(pj, 'r') as f:
                log.info(f"{config.algorithm} params were updated using {param_json}")
                param_dicts.append(json.load(f))
        config.algorithm_args = {k: v for D in param_dicts for k, v in D.items()}
    param_str = f'Learning {config.algorithm} model with the following params:\n'
    for param, value in config.algorithm_args.items():
        param_str += "{}\t= {}\n".format(param, value)

    log.info(f"{param_str}")

    targets_all, x_all = _load_data(config, partitions)

    if config.cross_validate:
        run_crossval(x_all, targets_all, config)

    log.info("Learning full {} model".format(config.algorithm))
    model = ls.learn.local_learn_model(x_all, targets_all, config)

    ls.mpiops.run_once(ls.geoio.export_model, model, config)

    # use trained model
    if config.permutation_importance:
        ls.mpiops.run_once(ls.validate.permutation_importance, model, x_all,
                           targets_all, config)

    # if config.feature_importance:
    #     # model, x_all, targets_all, config: Config
    #     ls.mpiops.run_once(ls.validate.plot_feature_importance, model, x_all,
    #                        targets_all, config)

    log.info("Finished! Total mem = {:.1f} GB".format(_total_gb()))

def _load_data(config, partitions):
    if config.pickle_load:
        x_all = joblib.load(open(config.pickled_covariates, 'rb'))
        targets_all = joblib.load(open(config.pickled_targets, 'rb'))
        if config.cubist or config.multicubist:
            config.algorithm_args['feature_type'] = \
                joblib.load(open(config.featurevec, 'rb'))
        log.warning('Using  pickled targets and covariates. Make sure you have'
                    ' not changed targets file and/or covariates.')
    else:
        log.info('One or both pickled files were not '
                 'found. All targets will be intersected.')
        config.n_subchunks = partitions
        if config.n_subchunks > 1:
            log.info("Memory constraint forcing {} iterations "
                     "through data".format(config.n_subchunks))
        else:
            log.info("Using memory aggressively: "
                     "dividing all data between nodes")

        # Make the targets
        if config.train_data_pk and exists(config.train_data_pk):
            log.info('Reusing pickled training data')
            image_chunk_sets, transform_sets, targets = \
                joblib.load(open(config.train_data_pk, 'rb'))
        else:
            log.info('Intersecting targets as pickled train data was not '
                     'available')
            targets = ls.geoio.load_targets(shapefile=config.target_file,
                                            targetfield=config.target_property,
                                            conf=config)
            # Get the image chunks and their associated transforms
            image_chunk_sets = ls.geoio.image_feature_sets(targets, config)
            transform_sets = [k.transform_set for k in config.feature_sets]

        if config.rawcovariates:
            log.info('Saving raw data before any processing')
            ls.features.save_intersected_features_and_targets(image_chunk_sets,
                                                              transform_sets,
                                                              targets, config)

        if config.train_data_pk:
            joblib.dump([image_chunk_sets, transform_sets, targets],
                        open(config.train_data_pk, 'wb'))

        if config.rank_features:
            measures, features, scores = ls.validate.local_rank_features(
                image_chunk_sets,
                transform_sets,
                targets,
                config)
            ls.mpiops.run_once(ls.geoio.export_feature_ranks, measures,
                               features, scores, config)

        # need to add cubist cols to config.algorithm_args
        # keep: bool array corresponding to rows that are retained
        features, keep = ls.features.transform_features(image_chunk_sets,
                                                        transform_sets,
                                                        config.final_transform,
                                                        config)
        # learn the model
        # local models need all data
        x_all = ls.features.gather_features(features[keep], node=0)

        # We're doing local models at the moment
        targets_all = ls.targets.gather_targets(targets, keep, config, node=0)

        if config.pickle and ls.mpiops.chunk_index == 0:
            if hasattr(config, 'pickled_covariates'):
                joblib.dump(x_all, open(config.pickled_covariates, 'wb'))
            if hasattr(config, 'pickled_targets'):
                joblib.dump(targets_all, open(config.pickled_targets, 'wb'))

    return targets_all, x_all

@cli.command()
@click.argument('pipeline_file', type=click.Path(exists=True))
@click.option('-j', '--param_json', type=click.Path(exists=True), multiple=True,
              help='algorithm parameters json, possibly from a previous optimise job')
@click.option('-p', '--partitions', type=int, default=1,
              help='divide each node\'s data into this many partitions')
def optimise(pipeline_file: str, param_json: str, partitions: int) -> None:
    """Optimise model parameters using Bayesian regression."""
    if uncoverml.mpiops.chunks > 1 and ('PBS_NNODES' in os.environ and int(os.environ['PBS_NNODES']) > 1):
        raise NotImplementedError("Currently optimiser does not work with mpi. \n"
                                  "However it can utilise a whole NCI node with many CPUs!")
    config = ls.config.Config(pipeline_file)
    if len(param_json) > 0:
        # log.info('the following jsons are used as params \n', click.echo('\n'.join(param_json)))
        param_dicts = []  # list of dicts
        for pj in param_json:
            with open(pj, 'r') as f:
                log.info(f"{config.algorithm} params were updated using {param_json}")
                param_dicts.append(json.load(f))
        config.algorithm_args = {k: v for D in param_dicts for k, v in D.items()}

    param_str = f'Optimising {config.algorithm} model with the following base params:\n'
    for param, value in config.algorithm_args.items():
        param_str += "{}\t= {}\n".format(param, value)
    log.info(param_str)
    targets_all, x_all = _load_data(config, partitions)
    if ls.mpiops.chunk_index == 0:
        if config.hpopt:
            log.info("Using hyperopt package to optimise model params")
            hyopt.optimise_model(x_all, targets_all, config)
        else:
            log.info("Using scikit-optimise package to optimise model params")
            optimisation.bayesian_optimisation(x_all, targets_all, config)

    ls.mpiops.comm.barrier()
    log.info("Finished! Total mem = {:.1f} GB".format(_total_gb()))

    log.info("Finished optimisation of model parameters!")

@cli.command()
@click.argument('pipeline_file')
@click.option('-s', '--subsample_fraction', type=float, default=1.0,
              help='only use this fraction of the data for learning classes')
def cluster(pipeline_file, subsample_fraction):
    config = ls.config.Config(pipeline_file)

    for f in config.feature_sets:
        if not f.transform_set.global_transforms:
            raise ValueError('Standardise transform must be used for kmeans')
        for t in f.transform_set.global_transforms:
            if not isinstance(t, StandardiseTransform):
                raise ValueError('Only standardise transform is '
                                 'allowed for kmeans')

    config.subsample_fraction = subsample_fraction
    if config.subsample_fraction < 1:
        log.info("Memory contstraint: using {:2.2f}%"
                 " of pixels".format(config.subsample_fraction * 100))
    else:
        log.info("Using memory aggressively: dividing all data between nodes")

    if config.semi_supervised:
        semisupervised(config)
    else:
        unsupervised(config)
    log.info("Finished! Total mem = {:.1f} GB".format(_total_gb()))

def semisupervised(config):

    # make sure we're clear that we're clustering
    config.algorithm = config.clustering_algorithm
    config.cubist = False
    # Get the taregts
    targets = ls.geoio.load_targets(shapefile=config.class_file,
                                    targetfield=config.class_property,
                                    conf=config)

    # Get the image chunks and their associated transforms
    image_chunk_sets = ls.geoio.semisupervised_feature_sets(targets, config)
    transform_sets = [k.transform_set for k in config.feature_sets]

    features, _ = ls.features.transform_features(image_chunk_sets,
                                                 transform_sets,
                                                 config.final_transform,
                                                 config)
    features, classes = ls.features.remove_missing(features, targets)
    indices = np.arange(classes.shape[0], dtype=int)

    config.n_classes = ls.cluster.compute_n_classes(classes, config)
    model = ls.cluster.KMeans(config.n_classes, config.oversample_factor)
    log.info("Clustering image")
    model.learn(features, indices, classes)
    ls.mpiops.run_once(ls.geoio.export_cluster_model, model, config)

def unsupervised(config):
    # make sure we're clear that we're clustering
    config.algorithm = config.clustering_algorithm
    config.cubist = False
    # Get the image chunks and their associated transforms
    image_chunk_sets = ls.geoio.unsupervised_feature_sets(config)
    transform_sets = [k.transform_set for k in config.feature_sets]

    features, _ = ls.features.transform_features(image_chunk_sets,
                                                 transform_sets,
                                                 config.final_transform,
                                                 config)

    features, _ = ls.features.remove_missing(features)
    model = ls.cluster.KMeans(config.n_classes, config.oversample_factor)
    log.info("Clustering image")
    model.learn(features)
    ls.mpiops.run_once(ls.geoio.export_cluster_model, model, config)

@cli.command()
@click.argument('pipeline_file')
@click.argument('model_or_cluster_file')
@click.option('-p', '--partitions', type=int, default=1,
              help='divide each node\'s data into this many partitions')
def validate(pipeline_file, model_or_cluster_file, partitions):
    """Validate a model with out of sample shapefile."""
    with open(model_or_cluster_file, 'rb') as f:
        state_dict = joblib.load(f)

    model = state_dict["model"]

    config = ls.config.Config(pipeline_file)
    config.pickle_load = False
    config.target_file = config.oos_validation_file
    config.target_property = config.oos_validation_property

    targets_all, x_all = _load_data(config, partitions)

    ls.validate.oos_validate(targets_all, x_all, model, config)

    log.info("Finished OOS validation job! Total mem = {:.1f} GB".format(_total_gb()))

@cli.command()
@click.argument('model_or_cluster_file')
@click.option('-p', '--partitions', type=int, default=1,
              help='divide each node\'s data into this many partitions')
@click.option('-m', '--mask', type=str, default='',
              help='mask file used to limit prediction area')
@click.option('-r', '--retain', type=int, default=None,
              help='mask values where to predict')
def predict(model_or_cluster_file, partitions, mask, retain):

    with open(model_or_cluster_file, 'rb') as f:
        state_dict = joblib.load(f)

    model = state_dict["model"]
    config = state_dict["config"]
    config.cluster = True if splitext(model_or_cluster_file)[1] == '.cluster' \
        else False
    config.mask = mask if mask else config.mask
    if config.mask:
        config.retain = retain if retain else config.retain

        if not isfile(config.mask):
            config.mask = ''
            log.info('A mask was provided, but the file does not exist on '
                     'disc or is not a file.')

    config.n_subchunks = partitions
    if config.n_subchunks > 1:
        log.info("Memory contstraint forcing {} iterations "
                 "through data".format(config.n_subchunks))
    else:
        log.info("Using memory aggressively: dividing all data between nodes")

    image_shape, image_bbox, image_crs = ls.geoio.get_image_spec(model, config)

    outfile_tif = config.name + "_" + config.algorithm
    predict_tags = model.get_predict_tags()
    if not config.outbands:
        config.outbands = len(predict_tags)

    image_out = ls.geoio.ImageWriter(image_shape, image_bbox, image_crs,
                                     outfile_tif,
                                     config.n_subchunks, config.output_dir,
                                     band_tags=predict_tags[
                                               0: min(len(predict_tags),
                                                      config.outbands)],
                                     **config.geotif_options)

    for i in range(config.n_subchunks):
        log.info("starting to render partition {}".format(i+1))
        ls.predict.render_partition(model, i, image_out, config)

    # explicitly close output rasters
    image_out.close()

    if config.cluster and config.cluster_analysis:
        if ls.mpiops.chunk_index == 0:
            ls.predict.final_cluster_analysis(config.n_classes,
                                              config.n_subchunks)

    # ls.predict.final_cluster_analysis(config.n_classes,
    #                                   config.n_subchunks)

    if config.thumbnails:
        image_out.output_thumbnails(config.thumbnails)
    log.info("Finished! Total mem = {:.1f} GB".format(_total_gb()))

@cli.command()
@click.argument('pipeline_file')
@click.option('-p', '--partitions', type=int, default=1,
              help='divide each node\'s data into this many partitions')
@click.option('-m', '--mask', type=str, default='',
              help='mask file used to limit prediction area')
@click.option('-r', '--retain', type=int, default=None,
              help='mask values where to predict')
def pca(pipeline_file, partitions, mask, retain):
    config = ls.config.Config(pipeline_file)
    assert config.pca, "Not a pca analysis. Please include the pca block in your yaml!!"
    config.mask = mask if mask else config.mask
    if config.mask:
        config.retain = retain if retain else config.retain

        if not isfile(config.mask):
            config.mask = ''
            log.info('A mask was provided, but the file does not exist on '
                     'disc or is not a file.')

    config.n_subchunks = partitions
    if config.n_subchunks > 1:
        log.info("Memory contstraint forcing {} iterations "
                 "through data".format(config.n_subchunks))
    else:
        log.info("Using memory aggressively: dividing all data between nodes")

    image_shape, image_bbox, image_crs = ls.geoio.get_image_spec_from_nchannels(config.n_components, config)
    tr_whiten = __validate_pca_config(config)
    outfile_tif = config.name + "_pca"

    image_out = ls.geoio.ImageWriter(image_shape, image_bbox, image_crs,
                                     outfile_tif,
                                     config.n_subchunks, config.output_dir,
                                     band_tags=[f'_pc_{n}' for n in range(1, config.n_components+1)],
                                     **config.geotif_options)

    for i in range(config.n_subchunks):
        log.info("starting to render partition {}".format(i+1))
        # TODO: ideally want to take a random sample of each covariate and compute whiten stats
        ls.predict.export_pca(i, image_out, config)

    # explicitly close output rasters
    image_out.close()

    log.info("Finished! Total mem = {:.1f} GB".format(_total_gb()))

def __validate_pca_config(config):
    transform_sets = [k.transform_set for k in config.feature_sets]
    for t in transform_sets:
        assert len(t.image_transforms) == 0  # validation that there are no image or global transforms
        assert len(t.global_transforms) == 0
    assert len(config.final_transform.global_transforms) == 1
    assert isinstance(config.final_transform.global_transforms[0], WhitenTransform)
    return config.final_transform.global_transforms[0]

def _total_gb():
    # given in KB so convert
    my_usage = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / (1024**2)
    # total_usage = mpiops.comm.reduce(my_usage, root=0)
    total_usage = ls.mpiops.comm.allreduce(my_usage)
    return total_usage

max_stars_repo_name: basaks/uncover-ml

max_issues_repo_name: basaks/uncover-ml

max_forks_repo_name: basaks/uncover-ml


text:
# Copyright (c) 2021 PaddlePaddle 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.

from __future__ import print_function

import unittest
import numpy as np
import paddle.fluid.core as core
from paddle.fluid.tests.unittests.op_test import OpTest, convert_float_to_uint16

from paddle.fluid.tests.unittests.test_conv2d_transpose_op import conv2dtranspose_forward_naive
from paddle import enable_static

def conv2d_bias_naive(out, bias):
    _, out_c, _, _ = out.shape

    for l in range(out_c):
        out[:, l, :, :] = out[:, l, :, :] + bias[l]
    return out

@unittest.skipIf(not core.supports_bfloat16(),
                 "place does not support BF16 evaluation")
class TestConv2DTransposeBF16MKLDNNOp(OpTest):
    def test_check_output(self):
        self.check_output_with_place(core.CPUPlace())

    def test_check_grad(self):
        pass

    def test_check_grad_no_input(self):
        pass

    def test_check_grad_no_filter(self):
        pass

    def init_op_type(self):
        self.data_format = "NCHW"
        self.op_type = 'conv2d_transpose'
        self._cpu_only = True

    def init_test_case(self):
        self.pad = [0, 0]
        self.fuse_bias = False
        self.use_mkldnn = True
        self.is_test = True
        self.bias_size = None
        self.fuse_activation = ""
        self.fuse_alpha = 0.0
        self.fuse_beta = 0.0
        self.stride = [1, 1]
        self.dilations = [1, 1]
        self.input_size = [2, 3, 5, 5]  # NCHW
        f_c = self.input_size[1]
        self.filter_size = [f_c, 6, 3, 3]
        self.groups = 1
        self.output_size = None
        self.output_padding = []
        self.data_format = "NCHW"
        self.pad = [0, 0]
        self.padding_algorithm = "EXPLICIT"
        self.force_fp32_output = False

    def setUp(self):
        self.input_type = np.uint16
        self.dtype = np.uint16
        self.mkldnn_data_type = "bfloat16"
        self.init_op_type()
        self.init_test_case()

        input = np.random.random(self.input_size).astype(np.float32)
        filter = np.random.random(self.filter_size).astype(np.float32)

        self.attrs = {
            'strides': self.stride,
            'paddings': self.pad,
            'padding_algorithm': self.padding_algorithm,
            'groups': self.groups,
            'dilations': self.dilations,
            'is_test': self.is_test,
            'use_mkldnn': self.use_mkldnn,
            'mkldnn_data_type': self.mkldnn_data_type,
            'force_fp32_output': self.force_fp32_output,
            'data_format': self.data_format,
            'fuse_activation': self.fuse_activation,
            'fuse_alpha': self.fuse_alpha,
            'fuse_beta': self.fuse_beta
        }
        if self.output_size is not None:
            self.attrs['output_size'] = self.output_size

        if len(self.output_padding) > 0:
            self.attrs['output_padding'] = self.output_padding

        output = conv2dtranspose_forward_naive(input, filter,
                                               self.attrs).astype(np.float32)

        if self.input_type is not np.float32:
            input = convert_float_to_uint16(input)

        self.inputs = {
            'Input': input.view(self.input_type),
            'Filter': OpTest.np_dtype_to_fluid_dtype(filter)
        }

        if self.fuse_bias and self.bias_size is not None:
            bias = np.random.random(self.bias_size).astype(np.float32)
            output = conv2d_bias_naive(output, bias)
            output = output.astype(np.float32)
            self.attrs['fuse_bias'] = self.fuse_bias
            self.inputs['Bias'] = OpTest.np_dtype_to_fluid_dtype(bias)

        if self.fuse_activation == "relu":
            output = np.maximum(output, 0).astype(np.float32)
        output = output.astype(np.float32)

        if not self.force_fp32_output:
            output = convert_float_to_uint16(output, self.attrs['data_format'])

        self.outputs['Output'] = output

class TestMKLDNNFuseBias(TestConv2DTransposeBF16MKLDNNOp):
    def init_test_case(self):
        super(TestMKLDNNFuseBias, self).init_test_case()
        self.pad = [1, 1]
        self.fuse_bias = True
        self.bias_size = [6]

class TestMKLDNNWithPad(TestConv2DTransposeBF16MKLDNNOp):
    def init_test_case(self):
        super(TestMKLDNNWithPad, self).init_test_case()
        self.pad = [1, 1]
        self.input_size = [2, 3, 10, 10]

class TestMKLDNNWithStride(TestConv2DTransposeBF16MKLDNNOp):
    def init_test_case(self):
        super(TestMKLDNNWithStride, self).init_test_case()
        self.pad = [1, 1]
        self.stride = [2, 2]
        self.input_size = [2, 3, 6, 6]  # NCHW

class TestMKLDNNWithAsymPad(TestConv2DTransposeBF16MKLDNNOp):
    def init_test_case(self):
        super(TestMKLDNNWithAsymPad, self).init_test_case()
        self.pad = [0, 0, 1, 2]
        self.padding_algorithm = "EXPLICIT"

class TestMKLDNNWithSamePad(TestConv2DTransposeBF16MKLDNNOp):
    def init_test_case(self):
        super(TestMKLDNNWithSamePad, self).init_test_case()
        self.pad = [0, 0]
        self.padding_algorithm = "SAME"

class TestMKLDNNWithValidPad(TestConv2DTransposeBF16MKLDNNOp):
    def init_test_case(self):
        super(TestMKLDNNWithValidPad, self).init_test_case()
        self.pad = [1, 1]
        self.padding_algorithm = "VALID"

class TestMKLDNNWithValidPad_NHWC(TestMKLDNNWithValidPad):
    def init_test_case(self):
        super(TestMKLDNNWithValidPad_NHWC, self).init_test_case()
        self.data_format = 'NHWC'
        N, C, H, W = self.input_size
        self.input_size = [N, H, W, C]

class TestConv2DTransposeMKLDNNWithDilationsExplicitPad(
        TestConv2DTransposeBF16MKLDNNOp):
    def init_test_case(self):
        super(TestConv2DTransposeMKLDNNWithDilationsExplicitPad,
              self).init_test_case()
        self.stride = [2, 1]
        self.dilations = [1, 2]
        self.groups = 1
        self.input_size = [4, 3, 8, 7]  # NCHW
        f_c = self.input_size[1]
        self.filter_size = [f_c, 6, 4, 3]
        self.pad = [1, 3, 2, 1]
        self.padding_algorithm = "EXPLICIT"

if __name__ == '__main__':
    enable_static()
    unittest.main()

max_stars_repo_name: zmxdream/Paddle

max_issues_repo_name: zmxdream/Paddle

max_forks_repo_name: zmxdream/Paddle


text:
import gym.envs.mujoco.humanoid as humanoid
import numpy as np

class HumanoidEnv(humanoid.HumanoidEnv):
    def _get_obs(self):
        data = self.sim.data
        return np.concatenate([data.qpos.flat[2:],
                               data.qvel.flat,
                               ])


text:
import numpy as np
from loader import *
from scipy import sparse
import math

import os, sys
utils_path = os.path.dirname(os.path.abspath(__file__)) + os.path.sep + "/../utils"
sys.path.append(utils_path)
import sparse_matrix

def loadAllCutSongciLines(path):
    f = open(path, "r")
    lines = f.readlines()
    cs = "".join(lines).split("$$$34567$$$")
    f.close()
    ret = []
    for line in lines:
        line = line.strip().strip('/').replace(',', '\n').replace('.', '\n')
        ls = line.split('\n')
        for l in ls:
            if len(l) == 0 or l == '$$$34567$$$':
                continue
            ret.append(l)
    return ret

dict_words = loadPZYDict('./output/words_pzy.txt')
word_idx = {k: i for i, k in enumerate(dict_words)}
dict_word_itmes = list(dict_words.items())
print(len(dict_word_itmes))

lines = loadAllCutSongciLines("./res/all_songci_cut.txt")
print("lines: ", len(lines))

cs_len = len(lines)
word_len = len(dict_words)
mat = np.zeros((word_len, cs_len), dtype='float32')

for i, line in enumerate(lines):
    words = line.split('/')
    for w in words:
        w = w.strip()
        if w in [',', '.'] or len(w) == 0:
            continue
        if word_idx.get(w) is None:
            #print('************************', w)
            continue
        else:
            mat[word_idx[w], i] += 1

coo = sparse.csr_matrix(mat)
sparse_matrix.save_sparse_csr("./output/w_c3.npy", coo)

reason: import numpy,from scipy


text:
import numpy as np
import tensorflow as tf
import matplotlib.pyplot as plt
from util.sprint import sfill, sfloat

NUM_OBS = 100
NUM_CPTS = 3
TRUE_CPTS = np.array([20, 30, 75])

# from the changepoints, we can compute segments
TRUE_SEGMENT_LENGTHS = np.diff([0] + TRUE_CPTS.tolist() + [NUM_OBS])
TRUE_SEGMENT_LABELS = [j for j in range(NUM_CPTS + 1)
                       for _ in range(TRUE_SEGMENT_LENGTHS[j])]
TRUE_SIGMA = np.array([2, 5, 1, 3])

TRUE_INTERCEPT = 1
TRUE_SLOPE = 5
TRUE_DELTAS = np.array([-9, 3, 1])
TRUE_PARAMS = np.array([TRUE_INTERCEPT, TRUE_SLOPE] + TRUE_DELTAS.tolist())

# seasonality represented by Fourier sinusoidal wave
TRUE_SEASONAL_PERIOD = 7
TRUE_SEASONAL_ORDER = 1
TRUE_SEASONAL_BETA = np.array([2, 2])

INIT_BETA_PARAMS = {'mean': 0.0, 'stddev': 1.0}
SMOOTHING_PARAM = 1.0
LEARNING_RATE = 0.001
MAX_ITER = 10000
TOL_PARAM, TOL_LOSS, TOL_GRAD = 1e-8, 1e-8, 1e-8
RANDOM_SEED = 0

MAX_CHARS = 15

def hinge(z):
    return np.max([0, z])

# generate time grid
t_grid = range(1, NUM_OBS + 1)

def evaluate_spline(t):
    """Evaluates spline with known knots at time t"""
    return TRUE_INTERCEPT + TRUE_SLOPE * t + \
           np.sum([delta * hinge(t - cpt)
                   for delta, cpt in zip(TRUE_DELTAS, TRUE_CPTS)])

# generate mean trend using spline formulation
np.random.seed(0)
mu_grid = [evaluate_spline(t) for t in t_grid]

def create_design_matrix(t_grid):
    """Creates design matrix for spline with known notes over grid of time"""
    columns = [np.repeat(1, NUM_OBS), t_grid]
    columns.extend([[hinge(t - cpt) for t in t_grid] for cpt in TRUE_CPTS])
    return np.stack(columns, axis=1)

# generate mean trend using design matrix formulation (equivalent to above)
np.random.seed(0)
mu_grid = np.matmul(create_design_matrix(t_grid), TRUE_PARAMS).reshape(-1, 1)

def fourier_expansion(t_grid, period, order):
    """Creates design matrix of Fourier terms evaluated at each t"""
    fourier_term_matrix = [
        trig_fun((2.0 * np.pi * (k + 1) * np.array(t_grid) / period))
        for k in range(order)
        for trig_fun in (np.sin, np.cos)
        ]
    return np.stack(fourier_term_matrix, axis=1)

# generate weekly seasonal trend using design matrix formulation
s_grid = np.matmul(fourier_expansion(t_grid,
                           period=TRUE_SEASONAL_PERIOD,
                           order=TRUE_SEASONAL_ORDER),
                   TRUE_SEASONAL_BETA).reshape(-1, 1)

# generate noise
e_grid = np.array([np.random.normal(loc=0.0,
                                    scale=TRUE_SIGMA[TRUE_SEGMENT_LABELS[t]])
                   for t in range(NUM_OBS)]).reshape(-1, 1)

# compute observations
y_obs = mu_grid + s_grid.sum(axis=1).reshape(-1, 1) + e_grid

# center and scale the data
CENTER = y_obs.mean()
SCALE = y_obs.std()
y_obs = (y_obs - CENTER) / SCALE

# tensor for data
y = tf.placeholder(dtype=tf.float32, shape=[None, 1])
x = tf.placeholder(dtype=tf.float32, shape=[None, NUM_OBS])

# tensors representing parameters
beta = tf.Variable(initial_value=tf.random_normal(shape=[NUM_OBS, 1],
                                                  seed=RANDOM_SEED,
                                                  **INIT_BETA_PARAMS),
                   dtype=tf.float32)

# differnece matrix is constant
D = np.zeros(shape=[NUM_OBS - 2, NUM_OBS])
for index_row in range(NUM_OBS - 2):
    D[index_row, index_row] = 1.0
    D[index_row, index_row + 1] = -2.0
    D[index_row, index_row + 2] = 1.0
D = tf.constant(value=D, dtype=tf.float32)

# smoothing param
smoothing_param = tf.constant(SMOOTHING_PARAM, dtype=tf.float32)

# loss function
squared_error = tf.reduce_sum(tf.square(y - beta))
# squared_error = tf.reduce_sum(tf.square(y - tf.matmul(x, beta)))
regularizer = smoothing_param * tf.norm(tf.matmul(D, beta), ord=1)
fused_lasso_loss = squared_error + regularizer

# optimizer
optimizer = tf.train.AdamOptimizer(learning_rate=LEARNING_RATE)
train_op = optimizer.minimize(loss=fused_lasso_loss)

# gradient
grad_beta = tf.gradients(fused_lasso_loss, beta)

with tf.Session() as sess:
    sess.run(fetches=tf.global_variables_initializer())

    i = 1
    obs_beta = sess.run(fetches=[beta])
    obs_loss = sess.run(fetches=[fused_lasso_loss], feed_dict={y: y_obs})
    obs_grad_beta = sess.run(fetches=[grad_beta], feed_dict={y: y_obs})

    print(' {} | {} | {} | {}'
          .format(sfill('iter', len(str(MAX_ITER)), '>'),
                  sfill('|beta|', MAX_CHARS, '^'),
                  sfill('loss', MAX_CHARS, '^'),
                  sfill('grad', MAX_CHARS, '^')))

    while True:
        # gradient step
        sess.run(fetches=train_op, feed_dict={y: y_obs})

        # update parameters
        new_beta = sess.run(fetches=beta)
        diff_norm = np.linalg.norm(np.subtract(new_beta, obs_beta))

        # update loss
        new_loss = sess.run(fetches=fused_lasso_loss, feed_dict={y: y_obs})
        loss_diff = np.abs(new_loss - obs_loss[-1])

        # update gradient
        new_grad_beta = sess.run(fetches=grad_beta, feed_dict={y: y_obs})
        grad_norm = np.linalg.norm(new_grad_beta)

        obs_beta.append(new_beta)
        obs_loss.append(new_loss)
        obs_grad_beta.append(new_grad_beta)

        if (i - 1) % 100 == 0:
            print(' {} | {} | {} | {}'
                  .format(sfill(i, len(str(MAX_ITER))),
                          sfloat(np.linalg.norm(new_beta), MAX_CHARS),
                          sfloat(new_loss, MAX_CHARS),
                          sfloat(grad_norm, MAX_CHARS)))

        if diff_norm < TOL_PARAM:
            print('Parameter convergence in {} iterations!'.format(i))
            break

        if loss_diff < TOL_LOSS:
            print('Loss function convergence in {} iterations!'.format(i))
            break

        if grad_norm < TOL_GRAD:
            print('Gradient convergence in {} iterations!'.format(i))
            break

        if i >= MAX_ITER:
            print('Max number of iterations reached without convergence.')
            break

        i += 1

# visualize results
# plt.plot(t_grid, mu_grid, color='blue', label='trend')
plt.plot(t_grid, mu_grid + s_grid, color='green', label='trend + seasonal')
plt.scatter(x=t_grid, y=SCALE * y_obs + CENTER, s=7, color='black', marker='o',
            label='obs')
plt.plot(t_grid, SCALE * obs_beta[-1] + CENTER, color='red', label='est')
plt.legend(loc='upper right')
plt.show()

max_stars_repo_name: kyleclo/tensorflow-mle

max_issues_repo_name: kyleclo/tensorflow-mle

max_forks_repo_name: kyleclo/tensorflow-mle


text:
# Creator: Newell Moser
# Python 3 (compatible with Python 2)
# 7/16/2015
#
# EASTER EGG
#
# Simple text-based game that randomly places a "plane" in space. Based on
# user-specified launch angle and velocity, a projectile is launched with
# normal Earth gravity. The streamline is visualized, and the game is over when
# the projectile collides with the plane.
#
from __future__ import print_function
import random
import math
#import numpy as np

# Screen class
class Screen:
    def __init__(self):
        self.__width = 71 # INTEGERS; Screen is initially set to 51 characters wide
        self.__height = 36 # and 26 characters tall
        self.__xUnit = 2 # 2 meters/character in the x-direction
        self.__yUnit = 4 # 4 meters/character in the y-direction
        self.__win = [[' ' for m in range(0,self.__width)] for n in range(0,self.__height)]
        self.clear() # Clear the just initialized array representing the to-be printed "window"

    # Function methods to get or set data members/attributes
    def getWidth(self):
        return self.__width
    def getHeight(self):
        return self.__height
    def getXUnit(self):
        return self.__xUnit
    def getYUnit(self):
        return self.__yUnit
    def setWidth(self,width):
        self.__width = int(width)
        self.__win = [[' ' for m in range(0,self.__width)] for n in range(0,self.__height)]
        self.clear()
    def setHeight(self,height):
        self.__height = int(height)
        self.__win = [[' ' for m in range(0,self.__width)] for n in range(0,self.__height)]
        self.clear()
    def setXUnit(self,xUnit):
        self.__xUnit = xUnit
    def setYUnit(self,yUnit):
        self.__yUnit = yUnit

    # Clears the screen array to spaces and creates the border
    def clear(self):
        for m in range(0,self.__height):
            for n in range(0,self.__width):
                if m >= self.__height-1 and n <= 0:
                    self.__win[m][n] = '#'
                elif m <= 0 or m >= self.__height-1:
                    self.__win[m][n] = '-'
                elif n <= 0 or n >= self.__width-1:
                    self.__win[m][n] = '|'
                else:
                    self.__win[m][n] = ' '

    # Prints the screen array to the standard output (i.e. the screen)
    def show(self):
        print("\n"*80)
        print("*** PLANE DESTROYER ***     (Type -1 to exit)")
        for m in range(0,self.__height):
            for n in range(0,self.__width):
                print(self.__win[m][n], end="")
            print("")

    # Writes a text-character to a location in the screen array.
    # It's assumed that the location is given in meters, thus will
    # be automatically converted into character spacing.
    def addText(self, text, xLoc, yLoc):
        try:
            if len(text) > 1:
                print("Cannot add text to screen, must be only one character.")
        except TypeError:
            print("Oops, can only print text characters to the screen!")

        #print(xArr, math.floor(xArr))
        xArr = float(xLoc) / self.__xUnit
        yArr = (self.__height - 1.0) - float(yLoc) / self.__yUnit

        if xArr < 0 or xArr >= self.__width:
            print("X-location is out of the domain. Current max width is ", (self.__width-1)*self.__xUnit, " meters.")
        elif yArr < 0 or yArr >= self.__height:
            print("Y-location is out of the domain. Current max height is ", (self.__height-1)*self.__yUnit, " meters.")
        else:
            self.__win[self.regRound(yArr)][self.regRound(xArr)] = text
            return [self.regRound(xArr),self.regRound(yArr)]

    # Python 3.x uses Banker's rounding, this is a function to utilize
    # traditional rounding techniques
    def regRound(self,value):
        fracPart = math.modf(value)
        fracPart = fracPart[0]

        if fracPart < 0.5:
            return int(math.floor(value))
        else:
            return int(math.ceil(value))

# Python 3.x uses Banker's rounding, this is a function to utilize
# traditional rounding techniques
def regRound(value):
    fracPart = math.modf(value)
    fracPart = fracPart[0]

    if fracPart < 0.5:
        return int(math.floor(value))
    else:
        return int(math.ceil(value))

def checkCollision(shipXY,bulletXY):
    if bulletXY[0] >= shipXY[0]-1 and bulletXY[0] <= shipXY[0]+1:
        if bulletXY[1] == shipXY[1]:
            return True
    return False

def frange(start,stop, step=1.0):
    while start < stop:
        yield start
        start +=step

def runEasterEgg():
    PI = 355.0/113.0
    GRAV = 9.81

    # Initialize a Screen object
    scn = Screen()

    # Store maximum values for width and height [meters] that can be plotted
    xUnit = scn.getXUnit() #Unit conversions in meters/character
    yUnit = scn.getYUnit()
    maxWidth = (scn.getWidth() - 1) * xUnit
    maxHeight = (scn.getHeight() - 1) * yUnit

    # Initialize the ship's location
    shipX = random.randint(round(maxWidth/5),round(maxWidth - xUnit))
    shipY = random.randint(round(maxHeight/5),round(maxHeight - yUnit))

    scn.addText('<',shipX - xUnit,shipY)
    shipWin = scn.addText('=',shipX,shipY)
    scn.addText('>',shipX + xUnit,shipY)
    scn.show()

    # This is the main game loop
    collFlag = False
    exFlag = 1
    while exFlag != -1:

        print("Current location of enemy plane: (X,Y) = (",shipX,",",shipY,") meters")

        flag = False
        while flag == False:
            try:
                ang = float(input("Launch Angle (0 - 90deg): "))
                if ang == -1:
                    break
                spd = float(input("Projectile Speed (30 - 70 [m/s]): "))
                if spd == -1:
                    break
                elif ang >= 0 and ang <= 90 and spd > 0:
                    flag = True
                else:
                    print("Not valid numbers...")
            except ValueError:
                print("Not valid numbers...")

        if ang == -1 or spd == -1:
            exFlag = -1
            continue

        scn.clear()
        scn.addText('<',shipX - xUnit,shipY)
        scn.addText('=',shipX,shipY)
        scn.addText('>',shipX + xUnit,shipY)

        ang = ang*PI/180 #Convert to radians
        vx = spd*math.cos(ang)
        vy = spd*math.sin(ang)
        totT = 2*spd*math.sin(ang)/GRAV

        for x in frange(0.0,maxWidth,xUnit):
            tcur = x/vx
            ycur = vy*tcur-0.5*GRAV*tcur**2

            if ycur >= 0 and ycur <= maxHeight:
                bulletWin = scn.addText('*',x,ycur)

            collFlag = checkCollision(shipWin,bulletWin)
            if collFlag == True:
                break

        if collFlag == False:
            scn.show()
        else:
            scn.show()
            print("We got 'em Captain!")
            break

easterEgg = False


text:
"""
This module provides a useful way to sample from a stationary Gaussian process.

By assuming stationary, Gaussian, mean-zero increments, the process is uniquely
defined by its MSD. Here we use that observation to generate sample traces from
such processes, given the MSD.
"""

import os,sys
from copy import deepcopy

import numpy as np
from scipy.linalg import cholesky, toeplitz

from tracklib import Trajectory, TaggedSet

def sampleMSD(msd, n=1, isCorr=False, subtractMean=True):
    r"""
    Sample traces from a stationary Gaussian process with given MSD.

    Parameters
    ----------
    msd : (N, ) np.array
        the MSD from which to sample. We use the convention that ``msd[t]``
        should be the MSD at time lag ``t``, i.e. ``msd[0] = 0``.
    n : integer, optional
        Number of traces to generate
    isCorr : bool, optional
        whether the input array is MSD (default) or autocorrelation.
    subtractMean : bool
        whether to subtract the mean of each trace. See Notes below.

    Returns
    -------
    traj : (N, n) array
        The generated traces.

    See also
    --------
    dataset

    Notes
    -----
    Strictly speaking, only the ensemble of displacements is well-defined,
    because this is the one we assume steady state for. For the actual traces,
    we can add an arbitrary offset, and the reasonable thing to do here depends
    on the process we are sampling. For a purely diffusive process for example
    it makes sense to have them all start from zero, since there is no steady
    state anyways. If however, we are sampling from an MSD that plateaus at
    long times, we basically also have a steady state for the traces and thus
    the long-run mean converges.  In that case of course it makes sense to have
    this mean be identical for all traces. The option `!subtractMean` does
    exactly this: if ``True``, we subtract the mean from all generated traces,
    i.e. fix the ensemble mean to zero. Note that this is not precisely correct
    for finite-length traces, but probably the best we can do.

    **Algorithm:** from the MSD :math:`\mu(k)` we calculate the increment
    autocorrelation function :math:`\gamma(k)` as

    .. math:: \gamma(k) = \frac{1}{2} \left( \mu(|k+1|) + \mu(|k-1|) - 2\mu(|k|) \right) \,.

    The covariance matrix of the increment process we want to sample from is
    then the Toeplitz matrix with diagonals :math:`\gamma(k)`.
    """
    if not isCorr:
        msd[0] = 0
        msd = np.insert(msd, 0, msd[1])
        corr = 0.5 * (msd[2:] + msd[:-2] - 2*msd[1:-1])
    else: # pragma: no cover
        corr = msd

    try:
        L = cholesky(toeplitz(corr), lower=True)
    except np.linalg.LinAlgError: # pragma: no cover
        vals = np.linalg.eigvalsh(toeplitz(corr))
        raise RuntimeError("Correlation not positive definite. First 5 eigenvalues: {}".format(vals[:5]))
    steps = L @ np.random.normal(size=(len(corr), n))

    trajs = np.insert(np.cumsum(steps, axis=0), 0, 0, axis=0)
    if subtractMean:
        trajs = trajs - np.mean(trajs, axis=0)

    return trajs

def dataset(msd, N=1, Ts=None, d=3, **kwargs):
    """
    Generate a dataset of MSD sampled trajectories.

    All keyword arguments not mentioned below will be forwarded to `sampleMSD`.

    Parameters
    ----------
    msd : (N,) np.ndarray
        the MSD to sample from.
    N : int, optional
        number of particles per trajectory
    Ts : list of int, optional
        list of trajectory lengths, i.e. this determines number and length of
        trajectories. Any ``None`` entry will be replaced by the maximum possible
        value, ``len(msd)``. If there are values bigger than that, raises a
        `!ValueError`. If not specified, will default to 100 trajectories of
        maximum length.
    d : int, optional
        spatial dimension of the trajectories to sample

    Returns
    -------
    `TaggedSet` of `Trajectory`
        the generated data set.

    See also
    --------
    sampleMSD, tracklib.trajectory.Trajectory, tracklib.taggedset.TaggedSet

    Notes
    -----
    The input MSD is assumed to be the goal for the generated trajectories,
    i.e. incorporate the prefactor ``N*d``.
    """
    if Ts is None: # pragma: no cover
        Ts = 100*[None]

    Tmax = len(msd)
    for i, T in enumerate(Ts):
        if T is None:
            Ts[i] = Tmax
        elif T > Tmax: # pragma: no cover
            raise ValueError("Cannot sample trajectory of length {} from MSD of length {}".format(T, Tmax))
        elif T % 1 != 0: # pragma: no cover
            raise ValueError("Found non-integer trajectory length: {} (at index {})".format(T, i))

    # Timing execution of sampleMSD indicates that as long as we have a
    # reasonable distribution of lengths (e.g. exponential with mean 10% of
    # len(msd)), sampling all traces at the same time is faster, even though we
    # might generate a bunch of unused data.
    kwargs['n'] = len(Ts)*N*d
    traces = sampleMSD(msd/(N*d), **kwargs)

    def gen():
        for iT, T in enumerate(Ts):
            mytraces = [traces[:T, ((iT*N + n)*d):((iT*N + n+1)*d)] for n in range(N)]
            yield Trajectory.fromArray(mytraces)

    return TaggedSet(gen(), hasTags=False)

def control(dataset, msd=None):
    """
    Generate a stationary control data set to the one given.

    The control will look exactly like the original in all "meta"-aspects
    (number of trajectories, their length, any meta data, etc.), but the
    trajectories will be sampled from a stationary Gaussian process with the
    given MSD (or the empirical ensemble MSD of the original).

    The mean of each trajectory will be set to coincide with the mean of the
    sister trajectory it is generated from.

    Parameters
    ----------
    dataset : `TaggedSet` of `Trajectory`
        the dataset to generate a control for
    msd : (T,) np.ndarray, callable, or None
        the MSD to use for sampling, either directly (i.e. as array) or as
        function that will be evaluated on the needed lag times (where the unit
        of time is one frame). Note that this will be divided by
        (#loci)x(#dimensions) before sampling scalar traces, matching the usual
        notion of MSD of (e.g.) multidimensional trajectories.

    Returns
    -------
    `TaggedSet` of `Trajectory`
        the generated control data set.

    Notes
    -----
    Generation from empirical MSDs does not always work, since they might be
    noisy. In that case, provide a smoothed version as `!msd`.
    """
    if callable(msd):
        maxlen = max([len(traj) for traj in dataset])
        msd = msd(np.arange(maxlen))
    elif msd is None:
        from tracklib.analysis import MSD # bad style :( We cannot do the import at
                                          # import time (i.e. on top), because
                                          # the analysis module depends on the
                                          # models one.
        msd = MSD(dataset)
        del MSD

    N = dataset.map_unique(lambda traj : traj.N)
    d = dataset.map_unique(lambda traj : traj.d)
    msd = msd / (N*d) # reassign!

    def gen():
        for (traj, mytags) in dataset(giveTags=True):
            try:
                traces = sampleMSD(msd[:len(traj)], n=traj.N*traj.d, subtractMean=True)
            except np.linalg.LinAlgError: # pragma: no cover
                raise RuntimeError("Could not generate trajectories from provided (or ensemble) MSD. Try using something cleaner.")
            newdata = np.array([traces[:, (i*traj.d):((i+1)*traj.d)] for i in range(traj.N)])
            newdata += np.nanmean(traj.data, axis=1, keepdims=True)
            newdata[np.isnan(traj.data)] = np.nan

            newmeta = deepcopy(traj.meta)
            # Remove those meta entries that explicitly depend on the data
            for key in ['MSD', 'MSDmeta', 'chi2scores']:
                try:
                    del newmeta[key]
                except KeyError:
                    pass

            yield (Trajectory.fromArray(newdata, **newmeta), deepcopy(mytags))

    return TaggedSet(gen())

max_stars_repo_name: SGrosse-Holz/tracklib

max_issues_repo_name: SGrosse-Holz/tracklib

max_forks_repo_name: SGrosse-Holz/tracklib