text:
#encoding: utf-8
"""
pathint -- Path integration figure showing time-series data (without cues)

Created by Joe Monaco on 2010-10-07.

Copyright (c) 2009-2011 Johns Hopkins University. All rights reserved.

This software is provided AS IS under the terms of the Open Source MIT License.
See http://www.opensource.org/licenses/mit-license.php.
"""

# Library imports
import os
import numpy as np
from numpy import pi
import matplotlib as mpl
import matplotlib.pylab as plt

# Package imports
from ..core.analysis import BaseAnalysis
from ..vmo import VMOModel
from ..session import VMOSession
from ..placemap import CirclePlaceMap
from ..compare import correlation_matrix, correlation_diagonals
from ..tools.stats import smooth_pdf
from ..tools.filters import halfwave
from ..tools.images import array_to_image

class PathIntFigure(BaseAnalysis):

    """
    Create the example path integration dataset for producing the basic
    introductory figure to the path integration mechanism of the
    VMOModel model
    """

    label = 'path integration'

    def collect_data(self, **kwargs):
        """Run basic VMOModel path integration simulation and save relevant
        data for the intro figure of the paper.

        Additional keyword arguments are passed on to VMOModel.
        """
        # Set up and run the path integration model
        self.out('Running simulation...')
        model = VMOModel(N_outputs=500, N_theta=1000, N_cues=1, C_W=0.05,
            gamma_local=0, gamma_distal=0, **kwargs)
        model.advance()
        pm = model.post_mortem()

        # Save session data object for visualization
        session_fn = 'session.tar.gz'
        self.out('Computing and saving responses to %s...'%repr(session_fn))
        SD = VMOSession(model, save_I=True, save_E=True)
        SD.tofile(os.path.join(self.datadir, session_fn))
        self.results['session_file'] = session_fn

        # Save other simulation data
        data_fn = 'data.npz'
        self.out('Saving other data to %s...'%repr(data_fn))
        np.savez(os.path.join(self.datadir, data_fn), vel=pm.vel)
        self.results['data_file'] = data_fn

        # Good-bye!
        self.out('All done!')

    def create_plots(self, examples=[0,1], lap_inset=[3,4]):
        """Create figure with basic data panels
        """
        self.figure = {}
        figsize = 10, 12
        plt.rcParams['figure.figsize'] = figsize
        self.figure['pathint'] = f = plt.figure(figsize=figsize)
        f.suptitle(self.label.title())

        # Load the simulation data if necessary
        if 'SD' in self.results:
            SD = self.results['SD']
        else:
            self.results['SD'] = SD = VMOSession.fromfile(
                os.path.join(self.datadir, self.results['session_file']))
        if 'data' in self.results:
            data = self.results['data']
        else:
            self.results['data'] = data = np.load(
                os.path.join(self.datadir, self.results['data_file']))
        os.chdir(self.datadir)
        self.out.outfd = file('figure.log', 'w')

        # Track angle
        tlim = (SD.t[0], SD.t[-1])
        alpha_star = SD.alpha - SD.alpha[0]
        alpha_star[alpha_star<=0] += 2*pi
        ax = plt.subplot(12,1,1)
        ax.plot(SD.t, alpha_star, 'k-', lw=1.5)
        ax.set_xlim(tlim)
        ax.set_axis_off()

        # V_x
        vlim = (-40,40)
        bar_x = 100
        bar_y = vlim[0]+4
        ax = plt.subplot(12,1,2)
        ax.plot(SD.t, data['vel'][:,0,0], 'k-', lw=0.5, zorder=5)
        ax.plot(tlim, [0,0], '-', c='0.6', lw=1, zorder=2)
        ax.plot([bar_x,bar_x+30], [bar_y]*2, 'k-', lw=1.5, zorder=3) # scale bar: 30 s
        ax.plot([bar_x+30]*2, [bar_y,bar_y+10], 'k-', lw=1.5, zorder=3) # 10 cm/s
        self.out('Velocity trace scale bar = 30s x 10cm/s')
        ax.set_xlim(tlim)
        ax.set_ylim(vlim)
        ax.set_axis_off()

        # V_y
        ax = plt.subplot(12,1,3)
        ax.plot(SD.t, data['vel'][:,0,1], 'k-', lw=0.5, zorder=5)
        ax.plot(tlim, [0,0], '-', c='0.6', lw=1, zorder=2)
        ax.set_xlim(tlim)
        ax.set_ylim(vlim)
        ax.set_axis_off()

        # Ex 1
        ax1 = plt.subplot(12,1,4)
        ax1.plot(SD.t, SD.I_cache[examples[0]], 'b-', lw=0.5, zorder=5)
        # ax1.plot(tlim, [0,0], 'k-', lw=1, zorder=7)
        ax1.axis('tight')
        ax1.set_xlim(tlim)
        ax1.set_axis_off()

        # Ex 2
        ax2 = plt.subplot(12,1,5)
        ax2.plot(SD.t, SD.I_cache[examples[1]], 'b-', lw=0.5, zorder=5)
        # ax2.plot(tlim, [0,0], 'k-', lw=1, zorder=7)
        ax2.axis('tight')
        ax2.set_xlim(tlim)
        ax2.set_axis_off()

        # Set y limits for examples
        max_I = 1.03*np.max(np.absolute(SD.I_cache[examples]))
        ax1.set_ylim(-max_I, max_I)
        ax2.set_ylim(-max_I, max_I)

        # Lap indicators
        rect_kw = dict(fill=True, linewidth=0, zorder=0, alpha=0.2)
        cols = ['b', 'g']
        for sp in 2, 3, 4, 5:
            ax = plt.subplot(12,1,sp)
            for lap in xrange(len(SD.laps) - 1):
                lap_rect = mpl.patches.Rectangle((SD.laps[lap], vlim[0]),
                    SD.laps[lap+1]-SD.laps[lap], vlim[1]-vlim[0],
                    fc=cols[lap%2], **rect_kw)
                ax.add_artist(lap_rect)

        # Example insets
        crop = slice((SD.t >= SD.laps[lap_inset[0]]).nonzero()[0][0],
            (SD.t >= SD.laps[lap_inset[1]+1]).nonzero()[0][0])
        examples = np.asarray(examples)
        E_max = 1.05*np.absolute(SD.E_cache[examples,crop]).max()
        for subp in ((12,1,6), (12,1,7)), ((12,2,15), (12,2,17)):
            ax = plt.subplot(*subp[0])
            ax.plot(SD.t[crop], SD.I_cache[examples[0],crop], 'b-', lw=0.5)
            ax.plot(SD.t[crop], SD.E_cache[examples[0],crop], 'k-', lw=1.5, aa=True)
            ax.plot([SD.t[crop.start], SD.t[crop.stop]], [SD.thresh]*2, 'k--')
            ax.axis('tight')
            ax.set_ylim(-E_max, E_max)
            ax.set_axis_off()
            ax = plt.subplot(*subp[1])
            ax.plot(SD.t[crop], SD.I_cache[examples[1],crop], 'b-', lw=0.5)
            ax.plot(SD.t[crop], SD.E_cache[examples[1],crop], 'k-', lw=1.5, aa=True)
            ax.plot([SD.t[crop.start], SD.t[crop.stop]], [SD.thresh]*2, 'k--')
            ax.axis('tight')
            ax.set_ylim(-E_max, E_max)
            ax.set_axis_off()
            bar_x = SD.laps[lap_inset[0]]
            bar_y = -E_max+2
            ax.plot([bar_x,bar_x+5], [bar_y]*2, 'k-', lw=1.5, zorder=3) # scale bar: 5 s
            ax.plot([bar_x+5]*2, [bar_y,bar_y+5], 'k-', lw=1.5, zorder=3) # 5 Hz
        self.out('Input scale bar = 5s x 5Hz')

        # Start new figure
        self.figure['responses'] = f = plt.figure(figsize=figsize)
        f.suptitle('Output Responses')

        # Input peak distro
        E_max = SD.E_cache.max()
        ax = plt.subplot(421)
        ax.plot(*smooth_pdf(SD.E_cache.max(axis=1)), c='k', lw=1.5, aa=True)
        ax.plot([SD.thresh]*2, [0, 1.05*ax.get_ylim()[1]], 'r-', zorder=5)
        ax.axis('tight')
        ax.set_xlim(0, 1.05*E_max)

        # Input/output trajectories
        # ax = plt.subplot(8,6,4)
        # ax.scatter(x=SD.x[::4], y=SD.y[::4], c=SD.E_cache[examples[0],::4],
        #     s=1, lw=0, marker='o')
        # ax.axis('equal')
        # ax.set_axis_off()

        ax = plt.subplot(843)
        ax.scatter(x=SD.x[::4], y=SD.y[::4],
            c=halfwave(SD.E_cache[examples[0],::4]-SD.thresh),
            s=1, lw=0, marker='o')
        ax.axis('equal')
        ax.set_axis_off()

        ax = plt.subplot(844)
        ax.pcolor(SD.R_laps[examples[0]].T, cmap=mpl.cm.jet)
        self.out('Example 1, lap pcolor max = %.2f'%SD.R_laps[examples[0]].max())
        ax.axis('tight')
        ax.set_axis_off()
        array_to_image(np.flipud(SD.R_laps[examples[0]].T), 'ex1_pcolor.png',
            cmap=mpl.cm.jet)

        # ax = plt.subplot(8,6,10)
        # ax.scatter(x=SD.x[::4], y=SD.y[::4], c=SD.E_cache[examples[1],::4],
        #     s=1, lw=0, marker='o')
        # ax.axis('equal')
        # ax.set_axis_off()

        ax = plt.subplot(847)
        ax.scatter(x=SD.x[::4], y=SD.y[::4],
            c=halfwave(SD.E_cache[examples[1],::4]-SD.thresh),
            s=1, lw=0, marker='o')
        ax.axis('equal')
        ax.set_axis_off()

        ax = plt.subplot(848)
        ax.pcolor(SD.R_laps[examples[1]].T, cmap=mpl.cm.jet)
        self.out('Example 2, lap pcolor max = %.2f'%SD.R_laps[examples[1]].max())
        ax.axis('tight')
        ax.set_axis_off()
        array_to_image(np.flipud(SD.R_laps[examples[1]].T), 'ex2_pcolor.png',
            cmap=mpl.cm.jet)

        # Population responses
        # ax = plt.subplot(434)
        R = SD.get_population_matrix()
        # ax.imshow(R, aspect='auto', interpolation='nearest')
        # ax.set_axis_off()
        array_to_image(R, 'pop_response.png',
            cmap=mpl.cm.jet)
        self.out('Pop matrix max = %.2f'%SD.R.max())

        # ax = plt.subplot(435)
        PM = CirclePlaceMap(SD)
        # ax.imshow(PM.coverage_maps, cmap=mpl.cm.gray_r, aspect='auto',
        #     interpolation='nearest')
        # ax.set_axis_off()
        array_to_image(PM.coverage_maps, 'field_extents.png',
            cmap=mpl.cm.gray_r)
        self.out('Num active units = %d'%PM.num_active)

        # ax = plt.subplot(8, 3, 9)
        C = correlation_matrix(SD)
        # ax.imshow(np.flipud(C))
        # ax.axis('equal')
        # ax.set_axis_off()
        array_to_image(np.flipud(C), 'corr_matrix.png',
            cmap=mpl.cm.jet)

        # ax = plt.subplot(8, 3, 12)
        D = correlation_diagonals(C, centered=True)
        # ax.plot(D[0], D[1], 'k-', lw=1.5, aa=True)
        # ax.axis('tight')
        # ax.set_ylim(0, 1.05)
        width = (D[0,1]-D[0,0]) * (D[1]>=D[1].max()/2).sum()
        self.out('Corr width at half max = %.1f degrees'%width)

        # Population distros
        # ax = plt.subplot(8, 3, 13)
        # ax.plot(*smooth_pdf(PM.field_size, sd=10), c='k', lw=1.5)
        # ax.axis('tight')
        # newy = 1.1*ax.get_ylim()[1]
        # ax.set_xlim(-10, 200)
        # ax.set_ylim(ymax=newy)
        med_size = np.median(PM.field_size)
        # ax.plot([med_size]*2, [0, newy], 'r-', zorder=1)
        self.out('Number of fields = %d'%PM.num_fields.sum())
        self.out('Median field size = %.1f degrees'%med_size)

        # ax = plt.subplot(8, 3, 14)
        # ax.plot(*smooth_pdf(PM.maxima[:,1]), c='k', lw=1.5)
        # ax.axis('tight')
        # newy = 1.1*ax.get_ylim()[1]
        # ax.set_xlim(-1, 12)
        # ax.set_ylim(ymax=newy)
        med_rate = np.median(PM.maxima[:,1])
        # ax.plot([med_rate]*2, [0, newy], 'r-', zorder=1)
        self.out('Median firing rate = %.2f Hz'%med_rate)

        # ax = plt.subplot(8, 3, 15)
        # ax.plot(*smooth_pdf(SD.I_rate[SD.sortix], sd=0.25), c='k', lw=1.5)
        # ax.axis('tight')
        # newy = 1.1*ax.get_ylim()[1]
        # ax.set_xlim(0, 12)
        # ax.set_ylim(ymax=newy)
        med_info = np.median(SD.I_rate[SD.sortix])
        # ax.plot([med_info]*2, [0, newy], 'r-', zorder=1)
        self.out('Median spatial info = %.1f bits/spike'%med_info)

        plt.draw()
        plt.rcParams['figure.figsize'] = plt.rcParamsDefault['figure.figsize']
        self.out.outfd.close()

max_issues_repo_name: jdmonaco/vmo-feedback-model

max_forks_repo_name: jdmonaco/vmo-feedback-model

reason: import numpy,from numpy


text:
"""
(CRNN) models training

Training strategy:
------------------
1. the following pairs of classes will be treated the same:

2. the following classes will be determined by the special detectors:
    PR, LAD, RAD, LQRSV, Brady,
    (potentially) SB, STach

3. models will be trained for each tranche separatly:

4. one model will be trained using the whole dataset (consider excluding tranche C? good news is that tranche C mainly consists of "Brady" and "STach", which can be classified using the special detectors)

References: (mainly tips for faster and better training)
-----------
1. https://efficientdl.com/faster-deep-learning-in-pytorch-a-guide/
2. (optim) https://www.fast.ai/2018/07/02/adam-weight-decay/
3. (lr) https://spell.ml/blog/lr-schedulers-and-adaptive-optimizers-YHmwMhAAACYADm6F
4. more....
"""

import os, sys, time, textwrap, argparse, logging
from copy import deepcopy
from typing import Tuple, Dict, Any, List, NoReturn, Optional

import numpy as np
np.set_printoptions(precision=5, suppress=True)
import torch
from torch import nn
from torch import Tensor
import torch.nn.functional as F
from torch.nn.parallel import DistributedDataParallel as DDP, DataParallel as DP
from torch.utils.data.dataset import Dataset
from torch.utils.data import DataLoader

try:
    import torch_ecg
except ModuleNotFoundError:
    import sys
    from os.path import dirname, abspath
    sys.path.insert(0, dirname(dirname(dirname(abspath(__file__)))))

from torch_ecg.cfg import CFG, DEFAULTS
from torch_ecg.utils.utils_nn import default_collate_fn as collate_fn
from torch_ecg.utils.trainer import BaseTrainer

from model import ECG_CRNN_CINC2021
from dataset import CINC2021
from scoring_metrics import evaluate_scores
from cfg import BaseCfg, TrainCfg, ModelCfg

CINC2021.__DEBUG__ = False
ECG_CRNN_CINC2021.__DEBUG__ = False

if BaseCfg.torch_dtype == torch.float64:
    torch.set_default_tensor_type(torch.DoubleTensor)

class CINC2021Trainer(BaseTrainer):
    """
    """
    __name__ = "CINC2021Trainer"

    def __init__(self,
                 model:nn.Module,
                 model_config:dict,
                 train_config:dict,
                 device:Optional[torch.device]=None,
                 lazy:bool=True,
                 **kwargs:Any,) -> NoReturn:
        """ finished, checked,

        Parameters
        ----------
        model: Module,
            the model to be trained
        model_config: dict,
            the configuration of the model,
            used to keep a record in the checkpoints
        train_config: dict,
            the configuration of the training,
            including configurations for the data loader, for the optimization, etc.
            will also be recorded in the checkpoints.
            `train_config` should at least contain the following keys:
                "monitor": str,
                "loss": str,
                "n_epochs": int,
                "batch_size": int,
                "learning_rate": float,
                "lr_scheduler": str,
                    "lr_step_size": int, optional, depending on the scheduler
                    "lr_gamma": float, optional, depending on the scheduler
                    "max_lr": float, optional, depending on the scheduler
                "optimizer": str,
                    "decay": float, optional, depending on the optimizer
                    "momentum": float, optional, depending on the optimizer
        device: torch.device, optional,
            the device to be used for training
        lazy: bool, default True,
            whether to initialize the data loader lazily
        """
        super().__init__(model, CINC2021, model_config, train_config, device, lazy)

    def _setup_dataloaders(self, train_dataset:Optional[Dataset]=None, val_dataset:Optional[Dataset]=None) -> NoReturn:
        """ finished, checked,

        setup the dataloaders for training and validation

        Parameters
        ----------
        train_dataset: Dataset, optional,
            the training dataset
        val_dataset: Dataset, optional,
            the validation dataset
        """
        if train_dataset is None:
            train_dataset = self.dataset_cls(config=self.train_config, training=True, lazy=False)

        if self.train_config.debug:
            val_train_dataset = train_dataset
        else:
            val_train_dataset = None
        if val_dataset is None:
            val_dataset = self.dataset_cls(config=self.train_config, training=False, lazy=False)

         # https://discuss.pytorch.org/t/guidelines-for-assigning-num-workers-to-dataloader/813/4
        num_workers = 4

        self.train_loader = DataLoader(
            dataset=train_dataset,
            batch_size=self.batch_size,
            shuffle=True,
            num_workers=num_workers,
            pin_memory=True,
            drop_last=False,
            collate_fn=collate_fn,
        )

        if self.train_config.debug:
            self.val_train_loader = DataLoader(
                dataset=val_train_dataset,
                batch_size=self.batch_size,
                shuffle=True,
                num_workers=num_workers,
                pin_memory=True,
                drop_last=False,
                collate_fn=collate_fn,
            )
        else:
            self.val_train_loader = None
        self.val_loader = DataLoader(
            dataset=val_dataset,
            batch_size=self.batch_size,
            shuffle=True,
            num_workers=num_workers,
            pin_memory=True,
            drop_last=False,
            collate_fn=collate_fn,
        )

    def run_one_step(self, *data:Tuple[torch.Tensor, torch.Tensor]) -> Tuple[torch.Tensor, torch.Tensor]:
        """

        Parameters
        ----------
        data: tuple of Tensors,
            the data to be processed for training one step (batch),
            should be of the following order:
            signals, labels, *extra_tensors

        Returns
        -------
        preds: Tensor,
            the predictions of the model for the given data
        labels: Tensor,
            the labels of the given data
        """
        signals, labels = data
        signals = signals.to(self.device)
        labels = labels.to(self.device)
        preds = self.model(signals)
        return preds, labels

    @torch.no_grad()
    def evaluate(self, data_loader:DataLoader) -> Dict[str, float]:
        """
        """
        self.model.eval()

        all_scalar_preds = []
        all_bin_preds = []
        all_labels = []

        for signals, labels in data_loader:
            signals = signals.to(device=self.device, dtype=self.dtype)
            labels = labels.numpy()
            all_labels.append(labels)

            if torch.cuda.is_available():
                torch.cuda.synchronize()
            preds, bin_preds = self._model.inference(signals)
            all_scalar_preds.append(preds)
            all_bin_preds.append(bin_preds)

        all_scalar_preds = np.concatenate(all_scalar_preds, axis=0)
        all_bin_preds = np.concatenate(all_bin_preds, axis=0)
        all_labels = np.concatenate(all_labels, axis=0)
        classes = data_loader.dataset.all_classes

        if self.val_train_loader is not None:
            msg = f"all_scalar_preds.shape = {all_scalar_preds.shape}, all_labels.shape = {all_labels.shape}"
            self.log_manager.log_message(msg, level=logging.DEBUG)
            head_num = 5
            head_scalar_preds = all_scalar_preds[:head_num,...]
            head_bin_preds = all_bin_preds[:head_num,...]
            head_preds_classes = [np.array(classes)[np.where(row)] for row in head_bin_preds]
            head_labels = all_labels[:head_num,...]
            head_labels_classes = [np.array(classes)[np.where(row)] for row in head_labels]
            for n in range(head_num):
                msg = textwrap.dedent(f"""
                ----------------------------------------------
                scalar prediction:    {[round(n, 3) for n in head_scalar_preds[n].tolist()]}
                binary prediction:    {head_bin_preds[n].tolist()}
                labels:               {head_labels[n].astype(int).tolist()}
                predicted classes:    {head_preds_classes[n].tolist()}
                label classes:        {head_labels_classes[n].tolist()}
                ----------------------------------------------
                """)
                self.log_manager.log_message(msg)

        auroc, auprc, accuracy, f_measure, f_beta_measure, g_beta_measure, challenge_metric = \
            evaluate_scores(
                classes=classes,
                truth=all_labels,
                scalar_pred=all_scalar_preds,
                binary_pred=all_bin_preds,
            )
        eval_res = dict(
            auroc=auroc,
            auprc=auprc,
            accuracy=accuracy,
            f_measure=f_measure,
            f_beta_measure=f_beta_measure,
            g_beta_measure=g_beta_measure,
            challenge_metric=challenge_metric,
        )

        # in case possible memeory leakage?
        del all_scalar_preds, all_bin_preds, all_labels

        self.model.train()

        return eval_res

    @property
    def batch_dim(self) -> int:
        """
        batch dimension, for CinC2021, it is 0,
        """
        return 0

    @property
    def extra_required_train_config_fields(self) -> List[str]:
        """
        """
        return []

    @property
    def save_prefix(self) -> str:
        return f"{self._model.__name__}_{self.model_config.cnn.name}_epoch"

    def extra_log_suffix(self) -> str:
        return super().extra_log_suffix() + f"_{self.model_config.cnn.name}"

def get_args(**kwargs:Any):
    """
    """
    cfg = deepcopy(kwargs)
    parser = argparse.ArgumentParser(
        description="Train the Model on CINC2021",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    parser.add_argument(
        "-l", "--leads",
        type=int, default=12,
        help="number of leads",
        dest="n_leads")
    parser.add_argument(
        "-t", "--tranches",
        type=str, default="",
        help="the tranches for training",
        dest="tranches_for_training")
    parser.add_argument(
        "-b", "--batch-size",
        type=int, default=128,
        help="the batch size for training",
        dest="batch_size")
    parser.add_argument(
        "-c", "--cnn-name",
        type=str, default="multi_scopic_leadwise",
        help="choice of cnn feature extractor",
        dest="cnn_name")
    parser.add_argument(
        "-r", "--rnn-name",
        type=str, default="none",
        help="choice of rnn structures",
        dest="rnn_name")
    parser.add_argument(
        "-a", "--attn-name",
        type=str, default="se",
        help="choice of attention structures",
        dest="attn_name")
    parser.add_argument(
        "--keep-checkpoint-max", type=int, default=20,
        help="maximum number of checkpoints to keep. If set 0, all checkpoints will be kept",
        dest="keep_checkpoint_max")
    # parser.add_argument(
    #     "--optimizer", type=str, default="adam",
    #     help="training optimizer",
    #     dest="train_optimizer")
    parser.add_argument(
        "-d", "--debug", action="store_true",
        help="train with more debugging information",
        dest="debug")

    args = vars(parser.parse_args())

    cfg.update(args)

    return CFG(cfg)

if __name__ == "__main__":
    train_config = get_args(**TrainCfg)

    tranches = train_config.tranches_for_training
    if tranches:
        classes = train_config.tranche_classes[tranches]
    else:
        classes = train_config.classes

    if train_config.n_leads == 12:
        model_config = deepcopy(ModelCfg.twelve_leads)
    elif train_config.n_leads == 6:
        model_config = deepcopy(ModelCfg.six_leads)
    elif train_config.n_leads == 4:
        model_config = deepcopy(ModelCfg.four_leads)
    elif train_config.n_leads == 3:
        model_config = deepcopy(ModelCfg.three_leads)
    elif train_config.n_leads == 2:
        model_config = deepcopy(ModelCfg.two_leads)
    model_config.cnn.name = train_config.cnn_name
    model_config.rnn.name = train_config.rnn_name
    model_config.attn.name = train_config.attn_name

    ECG_CRNN_CINC2021.__DEBUG__ = False
    model = ECG_CRNN_CINC2021(
        classes=classes,
        n_leads=train_config.n_leads,
        config=model_config,
    )

    if torch.cuda.device_count() > 1:
        model = DP(model)
        # model = DDP(model)
    model.to(device=DEFAULTS.device)

    trainer = CINC2021Trainer(
        model=model,
        model_config=model_config,
        train_config=train_config,
        device=DEFAULTS.device,
        lazy=False,
    )

    try:
        best_model_state_dict = trainer.train()
    except KeyboardInterrupt:
        try:
            sys.exit(0)
        except SystemExit:
            os._exit(0)

ext: py

reason: import numpy


text:
import numpy as np
from sklearn.datasets import load_iris

from pypadre.core.model.pipeline.pipeline import DefaultPythonExperimentPipeline
from pypadre.examples.base_example import example_app

app = example_app()

@app.dataset(name="iris", columns=['sepal length (cm)', 'sepal width (cm)', 'petal length (cm)',
                                   'petal width (cm)', 'class'], target_features='class')
def dataset():
    data = load_iris().data
    target = load_iris().target.reshape(-1, 1)
    return np.append(data, target, axis=1)

@app.parameter_map()
def config():
    return {"n_neighbors": 20}

@app.custom_splitter(name="Custom splitter", reference_git=__file__)
def custom_splitter(dataset, **kwargs):
    idx = np.arange(dataset.size[0])
    cutoff = int(len(idx) / 2)
    return idx[:cutoff], idx[cutoff:], None

@app.estimator(config=config, reference_git=__file__)
def estimator(X_train, y_train, *args, **kwargs):
    from sklearn.neighbors import KNeighborsClassifier
    knn = KNeighborsClassifier(*args)
    knn.fit(X_train, y_train)
    return knn

@app.evaluator(task_type="Classification", reference_git=__file__)
def evaluator(model, X_test, *args, **kwargs):
    y_predicted = model.predict(X_test)
    probabilities = None
    if hasattr(model, "predict_proba"):
        probabilities = model.predict_proba(X_test)
    return y_predicted, probabilities

@app.experiment(dataset=dataset, reference_git=__file__, splitting=custom_splitter,
                estimator=estimator, evaluator=evaluator,
                experiment_name="Iris KNN", project_name="Examples", ptype=DefaultPythonExperimentPipeline)
def experiment():
    pass

max_stars_repo_name: padre-lab-eu/pypadre

max_issues_repo_name: padre-lab-eu/pypadre

max_forks_repo_name: padre-lab-eu/pypadre


text:
import os
import unittest
import numpy as np

from numpy.testing import assert_array_equal

from amset.interpolation.interpolate import BoltzTraP2Interpolater
from pymatgen.io.vasp import Vasprun

test_dir = os.path.dirname(os.path.abspath(__file__))
amset_files = os.path.join(test_dir, '..', '..', '..', 'test_files',
                           'AlCuS2_mp-4979')

class TestBoltzTraP2Interpolater(unittest.TestCase):
    """Tests for interpolating a band structure using BoltzTraP2."""

    def setUp(self):
        vr = Vasprun(os.path.join(amset_files, 'vasprun.xml'))
        bs = vr.get_band_structure()
        num_electrons = vr.parameters['NELECT']

        self.kpoints = np.array(vr.actual_kpoints)
        self.interpolater = BoltzTraP2Interpolater(bs, num_electrons)

    def test_initialisation(self):
        """Test coefficients and parameters are calculated correctly."""
        self.interpolater.initialize()
        params = self.interpolater._parameters
        assert_array_equal(params[0][0], [[0, 0, 0]])
        self.assertAlmostEqual(params[1][0][0], 9.14055614)
        self.assertAlmostEqual(params[2][0][0].real, -2.33144546)

    def test_get_energies(self):
        """Test getting the interpolated energy, velocity and effective mass."""

        # test just getting energy
        energies = self.interpolater.get_energies(
            self.kpoints, 25, return_velocity=False,
            return_effective_mass=False)
        self.assertEqual(energies.shape, (138,))
        self.assertAlmostEqual(energies[0], 3.852399483908641)

        # test energy + velocity
        energies, velocities = self.interpolater.get_energies(
            self.kpoints, 25, return_velocity=True,
            return_effective_mass=False)
        self.assertEqual(energies.shape, (138,))
        self.assertAlmostEqual(energies[0], 3.852399483908641)
        self.assertEqual(velocities.shape, (138, 3))
        self.assertAlmostEqual(velocities[10][0], 5.401166156893334e+6,
                               places=0)

        # test energy + effective_mass
        energies, effective_masses = self.interpolater.get_energies(
            self.kpoints, 25, return_velocity=False,
            return_effective_mass=True)
        self.assertEqual(energies.shape, (138,))
        self.assertAlmostEqual(energies[0], 3.852399483908641)
        self.assertEqual(effective_masses.shape, (138, 3, 3))
        self.assertAlmostEqual(effective_masses[10][0][0], 0.2242569613626494)

        # test energy + velocity + effective_mass
        energies, velocities, effective_masses = self.interpolater.get_energies(
            self.kpoints, 25, return_velocity=True,
            return_effective_mass=True)
        self.assertEqual(energies.shape, (138,))
        self.assertAlmostEqual(energies[0], 3.852399483908641)
        self.assertEqual(velocities.shape, (138, 3))
        self.assertAlmostEqual(velocities[10][0], 5.401166156893334e+6,
                               places=0)
        self.assertEqual(effective_masses.shape, (138, 3, 3))
        self.assertAlmostEqual(effective_masses[10][0][0], 0.2242569613626494)

    def test_get_energies_scissor(self):
        """Test scissoring of band energies."""

        # test valence band
        energies = self.interpolater.get_energies(
            self.kpoints, 25, return_velocity=False,
            return_effective_mass=False, scissor=1.)
        self.assertEqual(energies.shape, (138,))
        self.assertAlmostEqual(energies[0], 3.852399483908641 - 0.5)

        # test conduction band
        energies = self.interpolater.get_energies(
            self.kpoints, 33, return_velocity=False,
            return_effective_mass=False, scissor=1.)
        self.assertEqual(energies.shape, (138,))
        self.assertAlmostEqual(energies[0], 7.301700765 + 0.5)

    def test_get_energies_multiple_bands(self):
        """Test getting the interpolated data for multiple bands."""

        # test just getting energy
        energies = self.interpolater.get_energies(
            self.kpoints, [25, 35], return_velocity=False,
            return_effective_mass=False)
        self.assertEqual(energies.shape, (2, 138,))
        self.assertAlmostEqual(energies[0][0], 3.852399483908641)
        self.assertAlmostEqual(energies[1][0], 9.594401616456384)

        # test energy + velocity
        energies, velocities = self.interpolater.get_energies(
            self.kpoints, [25, 35], return_velocity=True,
            return_effective_mass=False)
        self.assertEqual(energies.shape, (2, 138,))
        self.assertAlmostEqual(energies[0][0], 3.852399483908641)
        self.assertAlmostEqual(energies[1][0], 9.594401616456384)
        self.assertEqual(velocities.shape, (2, 138, 3))
        self.assertAlmostEqual(velocities[0][10][0], 5.401166156893334e+6,
                               places=0)
        self.assertAlmostEqual(velocities[1][10][0], 1.1686720831758736e+8,
                               places=0)

        # test energy + effective_mass
        energies, effective_masses = self.interpolater.get_energies(
            self.kpoints, [25, 35], return_velocity=False,
            return_effective_mass=True)
        self.assertEqual(energies.shape, (2, 138,))
        self.assertAlmostEqual(energies[0][0], 3.852399483908641)
        self.assertAlmostEqual(energies[1][0], 9.594401616456384)
        self.assertEqual(effective_masses.shape, (2, 138, 3, 3))
        self.assertAlmostEqual(effective_masses[0][10][0][0], 0.224256961362649)
        self.assertAlmostEqual(effective_masses[1][10][0][0], 0.009057103344700)

        # test energy + velocity + effective_mass
        energies, velocities, effective_masses = self.interpolater.get_energies(
            self.kpoints, [25, 35], return_velocity=True,
            return_effective_mass=True)
        self.assertEqual(energies.shape, (2, 138,))
        self.assertAlmostEqual(energies[0][0], 3.852399483908641)
        self.assertAlmostEqual(energies[1][0], 9.594401616456384)
        self.assertEqual(velocities.shape, (2, 138, 3))
        self.assertAlmostEqual(velocities[0][10][0], 5.401166156893334e+6,
                               places=0)
        self.assertAlmostEqual(velocities[1][10][0], 1.1686720831758736e+8,
                               places=0)
        self.assertEqual(effective_masses.shape, (2, 138, 3, 3))
        self.assertAlmostEqual(effective_masses[0][10][0][0], 0.224256961362649)
        self.assertAlmostEqual(effective_masses[1][10][0][0], 0.009057103344700)

    def test_get_energies_all_bands(self):
        # test all bands
        energies, velocities, effective_masses = self.interpolater.get_energies(
            self.kpoints, None, return_velocity=True,
            return_effective_mass=True)
        self.assertEqual(energies.shape, (96, 138,))
        self.assertAlmostEqual(energies[25][0], 3.852399483908641)
        self.assertAlmostEqual(energies[35][0], 9.594401616456384)
        self.assertEqual(velocities.shape, (96, 138, 3))
        self.assertAlmostEqual(velocities[25][10][0], 5.401166156893334e+6,
                               places=0)
        self.assertAlmostEqual(velocities[35][10][0], 1.1686720831758736e+8,
                               places=0)
        self.assertEqual(effective_masses.shape, (96, 138, 3, 3))
        self.assertAlmostEqual(effective_masses[25][10][0][0], 0.22425696136264)
        self.assertAlmostEqual(effective_masses[35][10][0][0], 0.00905710334470)

    def test_get_dos(self):
        """Test generating the interpolated DOS."""
        dos = self.interpolater.get_dos([10, 10, 10], emin=-10, emax=10,
                                        width=0.075)
        self.assertEqual(dos.shape, (20000, 2))
        self.assertEqual(dos[0][0], -10)
        self.assertAlmostEqual(dos[15000][1], 3.5362612128412807)

    def test_get_extrema(self):
        """Test getting the band structure extrema."""

        # test VBM
        extrema = self.interpolater.get_extrema(31, e_cut=1.)
        np.testing.assert_array_almost_equal(extrema[0], [0., 0.0, 0.0], 10)
        np.testing.assert_array_almost_equal(extrema[1],
                                             [0.0972, 0.4028, 0.0972], 4)
        np.testing.assert_array_almost_equal(extrema[2],
                                             [0.5, 0.5, -0.5], 10)
        np.testing.assert_array_almost_equal(extrema[3],
                                             [-0.0365, 0.0365, 0.5], 4)

        # test CBM
        extrema = self.interpolater.get_extrema(32, e_cut=1.)
        np.testing.assert_array_almost_equal(extrema[0], [0., 0.0, 0.0], 10)
        np.testing.assert_array_almost_equal(extrema[1], [0., 0.0, 0.5], 10)
        np.testing.assert_array_almost_equal(extrema[2],
                                             [0.5, 0.5, -0.5], 4)
        np.testing.assert_array_almost_equal(extrema[3],
                                             [0.4167, 0.4167, -0.0019], 4)

max_stars_repo_name: Navolo/amset

max_issues_repo_name: Navolo/amset

max_forks_repo_name: Navolo/amset


text:
#
import tensorflow as tf
import numpy as np
import auxiliary as aux
from termcolor import colored

from tensorflow.python.util.tf_export import tf_export

from tensorflow.python.ops import array_ops
from tensorflow.python.ops import init_ops
from tensorflow.python.ops import math_ops
from tensorflow.python.ops import nn_ops
from tensorflow.python.ops import variable_scope as vs

from tensorflow.python.ops.rnn_cell_impl import RNNCell
from tensorflow.contrib.layers import fully_connected

sigmoid = math_ops.sigmoid
tanh = math_ops.tanh
matm = math_ops.matmul
mul = math_ops.multiply
relu = nn_ops.relu
sign = math_ops.sign

def rotation_operator(x, y, n_hidden, eps = 1e-12):
	"""Rotation between two tensors: R(x,y) is unitary and takes x to y.

	Args:
		x: a tensor from where we want to start
		y: a tensor at which we want to finish
		eps: the cutoff for the normalizations (avoiding division by zero)
	Returns:
		a tensor, which is the orthogonal rotation operator R(x,y)
	"""

	size_batch = tf.shape(x)[0]
	hidden_size = n_hidden

	#construct the 2x2 rotation
	u = tf.nn.l2_normalize(x, 1, epsilon = eps)
	costh = tf.reduce_sum(u * tf.nn.l2_normalize(y, 1, epsilon = eps), 1)
	sinth = tf.sqrt(1 - costh ** 2)
	step1 = tf.reshape(costh, [size_batch, 1])
	step2 = tf.reshape(sinth, [size_batch, 1])
	Rth = tf.reshape(tf.concat([step1, -step2, step2, step1], axis = 1), [size_batch, 2, 2])

	#get v and concatenate u and v
	v = tf.nn.l2_normalize(y - tf.reshape(tf.reduce_sum(u * y, 1), [size_batch,1]) * u, 1, epsilon = eps)
	step3 = tf.concat([tf.reshape(u, [size_batch, 1, hidden_size]),
					  tf.reshape(v, [size_batch, 1, hidden_size])],
					  axis = 1)

	#do the batch matmul
	step4 = tf.reshape(u, [size_batch, hidden_size, 1])
	step5 = tf.reshape(v, [size_batch, hidden_size, 1])

	return (tf.eye(hidden_size, batch_shape = [size_batch]) -
		   tf.matmul(step4, tf.transpose(step4, [0,2,1])) -
		   tf.matmul(step5, tf.transpose(step5, [0,2,1])) +
		   tf.matmul(tf.matmul(tf.transpose(step3, [0,2,1]), Rth), step3))

def rotation_components(x, y, eps = 1e-12):
	"""Components for the operator R(x,y)
	   Together with `rotate` achieves best memory complexity: O(N_batch * N_hidden)

	Args:
		x: a tensor from where we want to start
		y: a tensor at which we want to finish
		eps: the cutoff for the normalizations (avoiding division by zero)
	Returns:
		Four components: u, v, [u,v] and R'(theta)
	"""

	size_batch = tf.shape(x)[0]
	hidden_size = tf.shape(x)[1]

	#construct the 2x2 rotation
	u = tf.nn.l2_normalize(x, 1, epsilon = eps)
	costh = tf.reduce_sum(u * tf.nn.l2_normalize(y, 1, epsilon = eps), 1)
	sinth = tf.sqrt(1 - costh ** 2)
	step1 = tf.reshape(costh, [size_batch, 1])
	step2 = tf.reshape(sinth, [size_batch, 1])
	Rth = tf.reshape(tf.concat([step1, -step2, step2, step1], axis = 1), [size_batch, 2, 2])

	#get v and concatenate u and v
	v = tf.nn.l2_normalize(y - tf.reshape(tf.reduce_sum(u * y, 1), [size_batch,1]) * u, 1, epsilon = eps)
	step3 = tf.concat([tf.reshape(u, [size_batch, 1, hidden_size]),
					  tf.reshape(v, [size_batch, 1, hidden_size])],
					  axis = 1)

	#do the batch matmul
	step4 = tf.reshape(u, [size_batch, hidden_size, 1])
	step5 = tf.reshape(v, [size_batch, hidden_size, 1])
	return step4, step5, step3, Rth

def rotate(v1, v2, v):
	"""Rotates v via the rotation R(v1,v2)

	Args:
		v: a tensor, which is the vector we want to rotate
		== to define R(v1,v2) ==
		v1: a tensor from where we want to start
		v2: a tensor at which we want to finish

	Returns:
		A tensor: the vector R(v1,v2)[v]
	"""
	size_batch = tf.shape(v1)[0]
	hidden_size = tf.shape(v1)[1]

	U  = rotation_components(v1, v2)
	h = tf.reshape(v, [size_batch, hidden_size, 1])

	return	(v + tf.reshape(
							- tf.matmul(U[0], tf.matmul(tf.transpose(U[0], [0,2,1]), h))
							- tf.matmul(U[1], tf.matmul(tf.transpose(U[1], [0,2,1]), h))
							+ tf.matmul(tf.transpose(U[2], [0,2,1]), tf.matmul(U[3], tf.matmul(U[2], h))),
							[size_batch, hidden_size]
						))

class RUMCell(RNNCell):
	"""Rotational Unit of Memory

	lambda = 0;
	uses `rotate` to implement the `Rotation` efficiently.
	"""

	def __init__(self,
				 hidden_size,
				 activation = None,
				 reuse = None,
				 kernel_initializer = None,
				 bias_initializer = None,
				 T_norm = None,
				 eps = 1e-12,
				 use_zoneout = False,
				 zoneout_keep_h = 0.9,
				 use_layer_norm = False,
				 is_training = False
				 ):
		"""Initialization of the RUM cell.

		Args:
			hidden_size: number of neurons in hidden state
			acitvation_tmp: activation of the temporary new state
			activation_tar: activation of the target
			activation_emb: activation of the embedded input
			T_norm: norm for time normalization, `eta` in the paper
			eps: the cutoff for the normalizations
			use_zoneout: zoneout, True or False
			use_layer_norm: batch normalization, True or False
			is_training: marker for the zoneout
		"""
		super(RUMCell, self).__init__(_reuse = reuse)
		self._hidden_size = hidden_size
		self._activation = activation or relu
		self._T_norm = T_norm
		self._kernel_initializer = kernel_initializer or aux.orthogonal_initializer(1.0)
		self._bias_initializer = bias_initializer
		self._eps = eps
		self._use_zoneout = use_zoneout
		self._zoneout_keep_h = zoneout_keep_h
		self._use_layer_norm = use_layer_norm
		self._is_training = is_training

	@property
	def state_size(self):
		return self._hidden_size

	@property
	def output_size(self):
		return self._hidden_size

	def call(self, inputs, state):
		with vs.variable_scope("gates"):
			bias_ones = self._bias_initializer
			if self._bias_initializer is None:
				dtype = [a.dtype for a in [inputs, state]][0]
				bias_ones = init_ops.constant_initializer(1.0, dtype = dtype)
			value = fully_connected(inputs=tf.concat([inputs, state],axis=1),
									num_outputs=2*self._hidden_size,
									activation_fn=None,
									biases_initializer=bias_ones,
									weights_initializer=aux.rum_ortho_initializer())
			r, u = array_ops.split(value = value, num_or_size_splits = 2, axis = 1)
			u = sigmoid(u)
			if self._use_layer_norm:
				concat = tf.concat([r, u], 1)
				concat = aux.layer_norm_all(concat, 2, self._hidden_size, "LN_r_u")
				r, u = tf.split(concat, 2, 1)
		with vs.variable_scope("candidate"):
			x_emb = fully_connected(inputs=inputs,
									num_outputs=self._hidden_size,
									activation_fn=None,
									biases_initializer=self._bias_initializer,
									weights_initializer=self._kernel_initializer)
			state_new = rotate(x_emb, r, state)
			if self._use_layer_norm:
				c = self._activation(aux.layer_norm(x_emb + state_new, "LN_c"))
			else:
				c = self._activation(x_emb + state_new)
		new_h = u * state + (1 - u) * c
		if self._T_norm != None:
			new_h = tf.nn.l2_normalize(new_h, 1, epsilon = self._eps) * self._T_norm
		if self._use_zoneout:
			new_h = aux.rum_zoneout(new_h, state, self._zoneout_keep_h, self._is_training)
		return new_h, new_h

	def zero_state(self, batch_size, dtype):
		h = tf.zeros([batch_size, self._hidden_size], dtype=dtype)
		return h

class ARUMCell(RNNCell):
	"""(Associative) Rotational Unit of Memory

	lambda = 1 (tuning in the associative memory);
	uses `rotation_operator` to update the rotation recursively.
	"""

	def __init__(self,
				 hidden_size,
				 activation = None,
				 reuse = None,
				 kernel_initializer = None,
				 bias_initializer = None,
				 T_norm = None,
				 eps = 1e-12,
				 use_zoneout = False,
				 zoneout_keep_h = 0.9,
				 use_layer_norm = False,
				 is_training = False,
				 lambda_pow = 0
				 ):
		"""Initialization of the Associative RUM cell.

		Args:
			hidden_size: number of neurons in hidden state
			acitvation_tmp: activation of the temporary new state
			activation_tar: activation of the target
			activation_emb: activation of the embedded input
			T_norm: norm for time normalization, `eta` in the paper
			eps: the cutoff for the normalizations
			use_zoneout: zoneout, True or False
			use_layer_norm: batch normalization, True or False
			is_training: marker for the zoneout
			lambda_pow: the power for the associative memory (an integer)
		"""
		super(ARUMCell, self).__init__(_reuse = reuse)
		self._hidden_size = hidden_size
		self._activation = activation or relu
		self._T_norm = T_norm
		self._kernel_initializer = kernel_initializer or aux.orthogonal_initializer(1.0)
		self._bias_initializer = bias_initializer
		self._eps = eps
		self._use_zoneout = use_zoneout
		self._zoneout_keep_h = zoneout_keep_h
		self._use_layer_norm = use_layer_norm
		self._is_training = is_training
		self._lambda_pow = lambda_pow

	@property
	def state_size(self):
		return self._hidden_size * (self._hidden_size + 1)

	@property
	def output_size(self):
		return self._hidden_size

	def call(self, inputs, state):
		#extract the associative memory and the state
		size_batch = tf.shape(state)[0]
		assoc_mem, state = tf.split(state, [self._hidden_size * self._hidden_size, self._hidden_size], 1)
		assoc_mem = tf.reshape(assoc_mem, [size_batch, self._hidden_size, self._hidden_size])

		with vs.variable_scope("gates"):
			bias_ones = self._bias_initializer
			if self._bias_initializer is None:
				dtype = [a.dtype for a in [inputs, state]][0]
				bias_ones = init_ops.constant_initializer(1.0, dtype = dtype)
			value = fully_connected(inputs=tf.concat([inputs, state],axis=1),
									num_outputs=2*self._hidden_size,
									activation_fn=None,
									biases_initializer=bias_ones,
									weights_initializer=aux.rum_ortho_initializer())
			r, u = array_ops.split(value = value, num_or_size_splits = 2, axis = 1)
			u = sigmoid(u)
			if self._use_layer_norm:
				concat = tf.concat([r, u], 1)
				concat = aux.layer_norm_all(concat, 2, self._hidden_size, "LN_r_u")
				r, u = tf.split(concat, 2, 1)
		with vs.variable_scope("candidate"):
			x_emb = fully_connected(inputs=inputs,
									num_outputs=self._hidden_size,
									activation_fn=None,
									biases_initializer=self._bias_initializer,
									weights_initializer=self._kernel_initializer)
			tmp_rotation = rotation_operator(x_emb, r, self._hidden_size)
			Rt = tf.matmul(assoc_mem, tmp_rotation)

			state_new = tf.reshape(tf.matmul(Rt, tf.reshape(state, [size_batch, self._hidden_size, 1])), [size_batch, self._hidden_size])
			if self._use_layer_norm:
				c = self._activation(aux.layer_norm(x_emb + state_new, "LN_c"))
			else:
				c = self._activation(x_emb + state_new)
		new_h = u * state + (1 - u) * c
		if self._T_norm != None:
			new_h = tf.nn.l2_normalize(new_h, 1, epsilon = self._eps) * self._T_norm
		if self._use_zoneout:
			new_h = aux.rum_zoneout(new_h, state, self._zoneout_keep_h, self._is_training)

		Rt = tf.reshape(Rt, [size_batch, self._hidden_size * self._hidden_size])
		new_state = tf.concat([Rt, new_h], 1)
		return new_h, new_state

	def zero_state(self, batch_size, dtype):
		e = tf.eye(self._hidden_size, batch_shape = [batch_size])
		e = tf.reshape(e, [batch_size, self._hidden_size * self._hidden_size])
		c = tf.zeros([batch_size, self._hidden_size], dtype=dtype)
		h = tf.concat([e, c], 1)
		return h

max_forks_repo_name: ekcole/RUM-Tensorflow


text:
# coding: utf-8

# ### Trabalho de SAD
# ### Aluno: Diego Dias Dos Santos Martins.  matrícula: 0050005813.
# ### Professor: Alex Salgado.
# * Preencha sua resposta no próprio arquivo do Jupyter e depois me envie no link que vou colocar no basecamp.

# # Tarefa simples de Classificação
#
# Referência: https://github.com/justmarkham/scikit-learn-videos/blob/master/03_getting_started_with_iris.ipynb
#
# 1 - importar os modulos python para machine learn e carregar o arquivo fruit_data_with_colors2.xlsx usando o método read_excel do pandas

# In[30]:

get_ipython().magic(u'matplotlib notebook')
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.model_selection import train_test_split
fruits = pd.read_excel('fruit_data_with_colors2.xlsx')

# 2 - Exibir os primeiros registros desta tabela

# In[31]:

fruits.head()

# ## Machine learning terminology
# Each row is an observation (also known as: sample, example, instance, record)
# Each column is a feature (also known as: predictor, attribute, independent variable, input, regressor, covariate)
#
# 3.1 - Quantas observações têm nessa base de dados?
# 3.2 - Quantas "features" têm nessa base de dados?

# In[32]:

##############################################################################################
#3.1 - 59 Observações\n"
#3.2 - 7 features"

# Each value we are predicting is the response (also known as: target, outcome, label, dependent variable)
# Classification is supervised learning in which the response is categorical
# Regression is supervised learning in which the response is ordered and continuous
#
# 4.1 - Quantas respostas temos nessa base?
# 4.2 - Para fazer nossas previsões em cima desta base, devemos usar algoritmo de Classificação ou Regressão?

# In[33]:

##############################################################################################

#4.1 - 59 respostas\n",
#4.2 - Algoritmo de Classificação para efetuar previsões exatas a respeito das frutas."

# # store feature matrix in "X"
# X = ???
#
# # store response vector in "y"
# y = ???
#
# 5 - Como você pode gerar a matriz X de feature e o vetor y de respostas? Dica: use essas 3 features para X ['mass', 'width', 'height'] e a resposta deve ser ['fruit_label'].

# In[34]:

X = fruits[['mass','width','height']]
y = fruits['fruit_label']

# ## Training a machine learning model with scikit-learn
#
#
# K-nearest neighbors (KNN) classification
# 6.1 - Usando o algoritmo de KNN (com 1 vizinho, k=1), qual a previsão? Adivinhe qual é a fruta com massa 20g, largura 4.3 cm, altura 5.5 cm, ou seja, com as seguintes features (mass = 20,	width=4.3,	height=5.5).
#
#
# Referência: https://github.com/justmarkham/scikit-learn-videos/blob/master/04_model_training.ipynb
#
#

# In[41]:

from sklearn.neighbors import KNeighborsClassifier
knn = KNeighborsClassifier(n_neighbors=1)
knn.fit(X, y)
knn.predict([[20, 4.3, 5.5]])

#

# In[44]:

from sklearn.neighbors import KNeighborsClassifier
knn = KNeighborsClassifier(n_neighbors=5)
knn.fit(X, y)
knn.predict([[20, 4.3, 5.5]])

# e 6.3 - LogisticRegression
#
#

# ## Evaluation procedure #1: Train and test on the entire dataset
# 7 - Usando o método de avaliação de acurácia (Treinar e testar na base de dados inteira), qual dos 3 métodos do item 6 é mais eficiente?
#
# Referência: https://github.com/justmarkham/scikit-learn-videos/blob/master/05_model_evaluation.ipynb
#

# 7.1 - Acurácia usando o algoritmo de KNN (com 1 vizinho, k=1)

# In[45]:

from sklearn import metrics
from sklearn.neighbors import KNeighborsClassifier
knn = KNeighborsClassifier(n_neighbors=1)
knn.fit(X, y)
y_pred = knn.predict(X)
print(metrics.accuracy_score(y, y_pred))

# 7.2 - Acurácia usando o algoritmo de KNN (com 5 vizinho, k=5)

# In[46]:

from sklearn import metrics
from sklearn.neighbors import KNeighborsClassifier
knn = KNeighborsClassifier(n_neighbors=5)
knn.fit(X, y)
y_pred = knn.predict(X)
print(metrics.accuracy_score(y, y_pred))

# 7.3 - Acurácia usando o algoritmo de LogisticRegression
#

# ## Evaluation procedure #2: Train/test split
# Usando o método de avaliação de acurácia (Treinar e testar SEPARADAMENTE), qual dos 3 métodos do item 6 é mais eficiente?
# Referência: https://github.com/justmarkham/scikit-learn-videos/blob/master/05_model_evaluation.ipynb

# 8.1 - Acurácia usando o algoritmo de KNN (com 1 vizinho, k=1)

# In[51]:

from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4, random_state=4)
knn = KNeighborsClassifier(n_neighbors=1)
knn.fit(X_train, y_train)
y_pred = knn.predict(X_test)
print(metrics.accuracy_score(y_test, y_pred))

# 8.2 - Acurácia usando o algoritmo de KNN (com 5 vizinho, k=5)

# In[52]:

from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4, random_state=4)
knn = KNeighborsClassifier(n_neighbors=5)
knn.fit(X_train, y_train)
y_pred = knn.predict(X_test)
print(metrics.accuracy_score(y_test, y_pred))

# 8.3 - Acurácia usando o algoritmo de LogisticRegression

# In[55]:

from sklearn.linear_model import LogisticRegression
logreg = LogisticRegression()
logreg.fit(X_train, y_train)
y_pred = logreg.predict(X_test)
print(metrics.accuracy_score(y_test, y_pred))

# ## Can we locate an even better value for K?
# Faça um programa python para testar diferentes valores de K entre 1 e 25 e diga através de gráfico, qual foi o melhor valor de K que obteve a melhor performance?
#

# In[58]:

k_range = list(range(1, 26))
scores = []
for k in k_range:
    knn = KNeighborsClassifier(n_neighbors=k)
    knn.fit(X_train, y_train)
    y_pred = knn.predict(X_test)
    scores.append(metrics.accuracy_score(y_test, y_pred))

import matplotlib.pyplot as plt

get_ipython().magic(u'matplotlib inline')

plt.plot(k_range, scores)
plt.xlabel('Value of K for KNN')
plt.ylabel('Testing Accuracy')

# Utilizando este seu último ajuste de K, faça novamente a seguinte previsão. Adivinhe qual é a fruta com massa 20g, largura 4.3 cm, altura 5.5 cm, ou seja, com as seguintes features (mass = 20, width=4.3, height=5.5

# In[60]:

knn = KNeighborsClassifier(n_neighbors=11)

knn.fit(X, y)

knn.predict([[20,4.3,5.5]])

# In[ ]:


text:
import numpy as np
from numpy.testing import assert_allclose
import unittest
from pb_bss_eval.distribution import GaussianTrainer

class TestGaussian(unittest.TestCase):
    def test_gaussian(self):
        samples = 10000
        mean = np.ones((3,))
        covariance = 2 * np.eye(3)
        x = np.random.multivariate_normal(mean, covariance, size=(samples,))
        model = GaussianTrainer().fit(x)
        assert_allclose(model.mean, mean, atol=0.1)
        assert_allclose(model.covariance, covariance, atol=0.1)

    def test_diagonal_gaussian(self):
        samples = 10000
        mean = np.ones((3,))
        covariance = 2 * np.eye(3)
        x = np.random.multivariate_normal(mean, covariance, size=(samples,))
        model = GaussianTrainer().fit(x, covariance_type="diagonal")
        assert_allclose(model.mean, mean, atol=0.1)
        assert_allclose(model.covariance, np.diag(covariance), atol=0.1)

    def test_spherical_gaussian(self):
        samples = 10000
        mean = np.ones((3,))
        covariance = 2 * np.eye(3)
        x = np.random.multivariate_normal(mean, covariance, size=(samples,))
        model = GaussianTrainer().fit(x, covariance_type="spherical")
        assert_allclose(model.mean, mean, atol=0.1)
        assert_allclose(
            model.covariance, np.mean(np.diag(covariance)), atol=0.1
        )


text:
import numpy as np
import tensorflow as tf
from tensorflow.contrib import rnn

#Tensorboard summaries
def variable_summaries(var, var_name):
    """Attach a lot of summaries to a Tensor (for TensorBoard visualization)."""
    with tf.name_scope('summaries'):
        mean = tf.reduce_mean(var)
        tf.summary.scalar('mean_'+ var_name, mean)
        with tf.name_scope('stddev'):
            stddev = tf.sqrt(tf.reduce_mean(tf.square(var - mean)))
        tf.summary.scalar('stddev_'+ var_name, stddev)
        tf.summary.histogram('histogram_' + var_name, var)

class LSTMCompSeq2Seq(object):
    def __init__(self, encoder_steps, feature_size, decoder_steps, label_size,
                 hidden_size, sess, auto_optimizer, pred_optimizer, reverse=False, unconditional_decode=False):

        self.reverse = reverse
        self.unconditional_decode = unconditional_decode
        self.sess = sess
        self.encoder_steps = encoder_steps
        self.feature_size = feature_size
        self.decoder_steps = decoder_steps
        self.label_size = label_size
        self.hidden_size = hidden_size

        self.auto_optimizer = auto_optimizer
        self.auto_lr = self.auto_optimizer._lr
        self.pred_optimizer = pred_optimizer
        self.pred_lr = self.pred_optimizer._lr

        #initializing placeholder
        inputs = tf.placeholder(tf.float32, [None, encoder_steps, feature_size], name='inputs_placeholder')
        labels = tf.placeholder(tf.float32, [None, decoder_steps, label_size], name='labels_placeholder')
        self.inputs = inputs
        self.labels = labels

        #reshape to [[batch_size, timestep],...,[b,t]]
        inputs = [tf.squeeze(step, [1]) for step in tf.split(inputs, encoder_steps, 1)]
        labels = [tf.squeeze(step, [1]) for step in tf.split(labels, decoder_steps, 1)]

        #init. LSTM Cells
        self._encoder_cell = rnn.MultiRNNCell([rnn.BasicLSTMCell(size) for size in hidden_size])
        self._decoder_cell = rnn.MultiRNNCell([rnn.BasicLSTMCell(size) for size in hidden_size])
        self._predictor_cell = rnn.MultiRNNCell([rnn.BasicLSTMCell(size) for size in hidden_size])

        #building graph
        encoder_outputs, encoder_state = self.encode(inputs)

        with tf.variable_scope('autoencoder'):
            decoder_outputs = self.decode(self._decoder_cell, inputs, self.feature_size, encoder_state, self.encoder_steps, self.unconditional_decode,
                                          self.reverse, is_guided=False)
        #auto loss
        self.autoencoder_loss = tf.losses.mean_squared_error(self.inputs, decoder_outputs)
        self.auto_grad = self.auto_optimizer.compute_gradients(self.autoencoder_loss)
        self.auto_train_op = self.auto_optimizer.apply_gradients(self.auto_grad)
        with tf.name_scope('auto-train'):
            self.auto_grad_summ_op = tf.summary.merge([tf.summary.histogram("%s-grad" % v.name, g) for g, v in self.auto_grad])
        tf.summary.scalar('auto_loss', self.autoencoder_loss)

        ###############
        ## predictor ##
        ###############
        with tf.variable_scope('predictor'):
            predictor_outputs = self.decode(self._predictor_cell, labels, self.label_size, encoder_state, self.decoder_steps, self.unconditional_decode,
                                           reverse=False, is_guided=True)

        #pred loss
        self.predictor_loss = tf.losses.mean_squared_error(self.labels, predictor_outputs)
        self.pred_grad = self.pred_optimizer.compute_gradients(self.predictor_loss)
        self.pred_grad = [(op, var) for op, var in self.pred_grad if op is not None]
        self.pred_train_op = self.pred_optimizer.apply_gradients(self.pred_grad)
        with tf.name_scope('pred-train'):
            self.pred_grad_summ_op = tf.summary.merge([tf.summary.histogram("%s-grad" % v.name, g) for g, v in self.pred_grad])
        tf.summary.scalar('pred_loss', self.predictor_loss)

    def train_autoencoder(self, inputs, learn_rate, merged):
        labels_holder = self.sess.run(tf.zeros([tf.shape(inputs)[0], self.decoder_steps, self.label_size], dtype=tf.float32))
        return self.sess.run([self.autoencoder_loss, merged, self.auto_grad_summ_op, self.auto_grad, self.auto_train_op],
                             feed_dict={self.inputs:inputs, self.labels:labels_holder, self.auto_lr:learn_rate})

    def train_predictor(self, inputs, labels, learn_rate, merged):
        return self.sess.run([self.predictor_loss, merged, self.pred_grad_summ_op, self.pred_grad, self.pred_train_op],
                            feed_dict={self.inputs:inputs, self.labels:labels, self.pred_lr:learn_rate})

    def encode(self, inputs):
        with tf.variable_scope('encoder', reuse=tf.AUTO_REUSE):
            encoder_outputs, encoder_state = rnn.static_rnn(self._encoder_cell, inputs, dtype=tf.float32)
            #tensorboard dist. on encoder_state
            for i in range(len(encoder_state)):
                variable_summaries(encoder_state[i][0], 'enc_hidden_'+str(i))
                variable_summaries(encoder_state[i][1], 'enc_output_'+str(i))
            return encoder_outputs, encoder_state

    def decode(self, lstmcell, labels, label_size, encoder_state, future_steps, unconditional_decode, reverse, is_guided):
        with tf.variable_scope('decoder', reuse=tf.AUTO_REUSE) as var_scope:
            #initializing variables
            d_weights = tf.get_variable('d_weights', shape=[self.hidden_size[-1], label_size],
                                        initializer=tf.truncated_normal_initializer, dtype=tf.float32)
            d_biases = tf.get_variable('d_biases', shape=[label_size],
                                       initializer=tf.ones_initializer, dtype=tf.float32)

            #tensorboard summary
            variable_summaries(d_weights, 'dec_weights')
            variable_summaries(d_biases, 'dec_biases')

            if is_guided:
                #guided training, feeding labels at each timestep t
                decoder_state = encoder_state #using last generated input
                init_input = [tf.zeros([tf.shape(decoder_state[-1][1])[0], label_size], dtype=tf.float32)]

                decoder_inputs_ = init_input + labels[:-1]
                decoder_outputs_, decoder_state = rnn.static_rnn(lstmcell, decoder_inputs_,
                                                                initial_state=decoder_state, dtype=tf.float32)

                decoder_outputs = []
                for i in range(len(decoder_outputs_)):
                    if i > 0:
                        var_scope.reuse_variables()
                    decoder_outputs.append(tf.matmul(decoder_outputs_[i], d_weights) + d_biases)
                outputs = tf.transpose(tf.stack(decoder_outputs), [1, 0, 2])
                return outputs

            else:
                decoder_state = encoder_state #using encoder state
                decoder_outputs = []
                decoder_input_ = tf.zeros([tf.shape(decoder_state[-1][1])[0], label_size],
                                          dtype=tf.float32) #initialize decoder input

                for step in range(future_steps):
                    if step > 0:
                        var_scope.reuse_variables()
                    decoder_input_, decoder_state = lstmcell(decoder_input_, decoder_state)
                    decoder_input_ = tf.matmul(decoder_input_, d_weights) + d_biases
                    decoder_outputs.append(decoder_input_)

                if reverse and (future_steps==self.encoder_steps):
                    decoder_outputs = decoder_outputs[::-1] #reverse outputs' order
                outputs = tf.transpose(tf.stack(decoder_outputs), [1, 0, 2])
                return outputs

    def predict(self, inputs, predict_price=True):
        """ input as [batch_size, steps, num_features] """
        cell = self._predictor_cell if predict_price else self._decoder_cell
        steps = self.decoder_steps if predict_price else self.encoder_steps
        output_size = self.label_size if predict_price else self.feature_size
        scope = 'predictor' if predict_price else 'autoencoder'

        inputs = [tf.squeeze(step, [1]) for step in tf.split(inputs, self.encoder_steps, 1)]
        _, encoder_state = self.encode(inputs)

        with tf.variable_scope(scope):
            decoder_outputs = self.decode(cell, None, output_size, encoder_state, steps, self.unconditional_decode,
                                          reverse=self.reverse, is_guided=False)
        return self.sess.run(decoder_outputs)

def generate_train_sample(X, y, batch_size, X_time_step, y_time_step):
    if len(y.shape) != 2:
        y = np.reshape(y, (len(y),1))
    available_idxs = range(len(X) - X_time_step - y_time_step)
    start_batch_idxs = np.random.choice(available_idxs, batch_size, replace=False)
    X_batch_idxs = [list(range(i, i+X_time_step)) for i in start_batch_idxs]
    y_batch_idxs = [list(range(i+X_time_step, i+X_time_step+y_time_step)) for i in start_batch_idxs]
    return np.take(X, X_batch_idxs, axis=0), np.take(y, y_batch_idxs, axis=0)

def generate_test_sample(X, y, total_samples, start, X_time_step, y_time_step):
    X_batch_idxs = [list(range(start+i, start+i+X_time_step)) for i in range(total_samples)]
    y_batch_idxs = [list(range(start+i+X_time_step, start+i+X_time_step+y_time_step)) for i in range(total_samples)]
    return np.take(X, X_batch_idxs, axis=0), np.take(y, y_batch_idxs, axis=0)

def remove_overlap(data, steps):
    d = data[:steps-1]
    return np.append(d,(data[range(steps-1,len(data),steps)]))


text:
import numpy as np
from chaosencrypt import logistic
from smartplot import smartplot

bits = 16
N = 500
M = 1000

AV = np.linspace(0,2**bits-1,N)
#AV = np.linspace(0,2**(bits-4)-1,N)
#AV = np.arange(0,100)
T = np.empty((N,M+1),dtype='int64')
NV = np.repeat(AV,M+1)

for i,A in enumerate(AV):
	print('step',i+1)
	A = int(A)
	dl = logistic.discrete(A,bits)
	#dl.step(2**(bits-1),M)
	#T[i] = dl.steps(None,M)
	T[i] = dl.steps(6871,M)

smartplot([NV,T.reshape(-1)],{'linewidth':0,'markersize':0.2,'axform':['0x%X','0x%X']})

max_stars_repo_name: nfejes/chaotic-image-encryption

max_issues_repo_name: nfejes/chaotic-image-encryption

max_forks_repo_name: nfejes/chaotic-image-encryption


text:
# incrementally remove points from the redetection effort spline and calculate AIC for each

import pickle
import numpy as np

import sys
sys.path.insert(0,'../../..') # allows us to import undetected extinctions package

from undetected_extinctions.redetection_effort.redetection_effort import negLogLikelihoodSpline

# some parameters
# ---

# time range over which we'll do the redetection effort
t0 = 1822 # the Wallich collection
tf = 2015 # taken as our last date

results_dir = '../../../results/redetection_effort/backwards_stepping/' # where to put intermediate results

# where databases are located
# ---

# fname_startwith = '../../../results/redetection_effort/initial_fit/initial_fit.pkl'
fname_startwith = results_dir + 'backwards_stepping_150.pkl' # to restart backwards stepping

# obtain the initial fit
# ---

# read in data
f = open(fname_startwith, 'rb')
ss = pickle.load( f )   # explanatory string
cs = pickle.load( f )   # c_t points defining the spline
ts = pickle.load( f )   # t points defining the spline
spps = pickle.load( f ) # redetections info for each species used in the fitting
f.close()

# check the output
assert( ts[0] == t0 )  # t0 = 1822, the Wallich collection
assert( ts[-1] == tf ) # tf = 2015, taken as our last date

# calculate AIC
negLL = negLogLikelihoodSpline(cs, ts, spps)
no_params = len(cs) + len(ts)
AIC = 2*no_params + 2*negLL

# incrementally remove point that provides best AIC, store results
# ---

# storage of each of our steps
tsV = [ ts ]
csV = [ cs ]
AICV = [ AIC ]
no_paramsV = [ no_params ]

while no_params > 120:

    # find the best point to remove

    best_AIC = np.inf # initialise

    for idx in range( 1, len(ts)-1 ): # we never remove endpoints of timeseries

        # new are the same as old ...
        new_ts = [ t for t in ts]
        new_cs = [ c for c in cs]
        # .. but with idx removed
        del new_ts[idx]
        del new_cs[idx]

        # calculate new AIC
        negLL = negLogLikelihoodSpline(new_cs, new_ts, spps)
        no_params = len(new_cs) + len(new_ts)
        new_AIC = 2*no_params + 2*negLL

        if new_AIC < best_AIC:
            best_ts = new_ts
            best_cs = new_cs
            best_AIC = new_AIC

    # store the best found for the next update

    ts = best_ts; cs = best_cs; AIC = best_AIC
    tsV.append( ts ); csV.append( cs ); AICV.append( AIC ); no_paramsV.append( no_params )

    print('no params: ' + str(no_params) + ' AIC: ' + str(AIC) )

    # save it to a file

    fName = results_dir + 'backwards_stepping_' + str(no_params) + '.pkl'
    f = open(fName, 'wb')
    # a string explaining the pickle file
    ss  = 'Created by backwards_stepping.py.\n'
    ss += 'Contains the following:\n'
    ss += '0. ss, string: this string you are reading now.\n'
    ss += '1. cs, list of floats: c_t points defining the spline.\n'
    ss += '2. ts: list ints: t points defining the spline.\n'
    ss += '3. spps, dictionary: keys are names, and values are dictionary with frst, last, and redetns.\n'
    ss += '4. AIC, float: the AIC calculated for this spline.\n'
    pickle.dump( ss, f )
    pickle.dump( cs, f )
    pickle.dump( ts, f )
    pickle.dump( spps, f )
    pickle.dump( AIC, f )
    f.close()


text:
from __future__ import division, absolute_import, print_function

import sys
import itertools

import numpy as np
from numpy.testing import run_module_suite, assert_, assert_raises, assert_equal

from numpy.core.multiarray_tests import solve_diophantine, internal_overlap
from numpy.lib.stride_tricks import as_strided
from numpy.compat import long

if sys.version_info[0] >= 3:
    xrange = range

ndims = 2
size = 10
shape = tuple([size] * ndims)

MAY_SHARE_BOUNDS = 0
MAY_SHARE_EXACT = -1

def _indices_for_nelems(nelems):
    """Returns slices of length nelems, from start onwards, in direction sign."""

    if nelems == 0:
        return [size // 2]  # int index

    res = []
    for step in (1, 2):
        for sign in (-1, 1):
            start = size // 2 - nelems * step * sign // 2
            stop = start + nelems * step * sign
            res.append(slice(start, stop, step * sign))

    return res

def _indices_for_axis():
    """Returns (src, dst) pairs of indices."""

    res = []
    for nelems in (0, 2, 3):
        ind = _indices_for_nelems(nelems)

        # no itertools.product available in Py2.4
        res.extend([(a, b) for a in ind for b in ind])  # all assignments of size "nelems"

    return res

def _indices(ndims):
    """Returns ((axis0_src, axis0_dst), (axis1_src, axis1_dst), ... ) index pairs."""

    ind = _indices_for_axis()

    # no itertools.product available in Py2.4

    res = [[]]
    for i in range(ndims):
        newres = []
        for elem in ind:
            for others in res:
                newres.append([elem] + others)
        res = newres

    return res

def _check_assignment(srcidx, dstidx):
    """Check assignment arr[dstidx] = arr[srcidx] works."""

    arr = np.arange(np.product(shape)).reshape(shape)

    cpy = arr.copy()

    cpy[dstidx] = arr[srcidx]
    arr[dstidx] = arr[srcidx]

    assert np.all(arr == cpy), 'assigning arr[%s] = arr[%s]' % (dstidx, srcidx)

def test_overlapping_assignments():
    """Test automatically generated assignments which overlap in memory."""

    inds = _indices(ndims)

    for ind in inds:
        srcidx = tuple([a[0] for a in ind])
        dstidx = tuple([a[1] for a in ind])

        yield _check_assignment, srcidx, dstidx

def test_diophantine_fuzz():
    # Fuzz test the diophantine solver
    rng = np.random.RandomState(1234)

    max_int = np.iinfo(np.intp).max

    for ndim in range(10):
        feasible_count = 0
        infeasible_count = 0

        min_count = 500//(ndim + 1)

        numbers = []
        while min(feasible_count, infeasible_count) < min_count:
            # Ensure big and small integer problems
            A_max = 1 + rng.randint(0, 11)**6
            U_max = rng.randint(0, 11)**6

            A_max = min(max_int, A_max)
            U_max = min(max_int-1, U_max)

            A = tuple(rng.randint(1, A_max+1) for j in range(ndim))
            U = tuple(rng.randint(0, U_max+2) for j in range(ndim))

            b_ub = min(max_int-2, sum(a*ub for a, ub in zip(A, U)))
            b = rng.randint(-1, b_ub+2)

            if ndim == 0 and feasible_count < min_count:
                b = 0

            X = solve_diophantine(A, U, b)

            if X is None:
                # Check the simplified decision problem agrees
                X_simplified = solve_diophantine(A, U, b, simplify=1)
                assert X_simplified is None, (A, U, b, X_simplified)

                # Check no solution exists (provided the problem is
                # small enough so that brute force checking doesn't
                # take too long)
                try:
                    ranges = tuple(xrange(0, a*ub+1, a) for a, ub in zip(A, U))
                except OverflowError:
                    # xrange on 32-bit Python 2 may overflow
                    continue

                size = 1
                for r in ranges:
                    size *= len(r)
                if size < 100000:
                    assert_(not any(sum(w) == b for w in itertools.product(*ranges)))
                    infeasible_count += 1
            else:
                # Check the simplified decision problem agrees
                X_simplified = solve_diophantine(A, U, b, simplify=1)
                assert X_simplified is not None, (A, U, b, X_simplified)

                # Check validity
                assert_(sum(a*x for a, x in zip(A, X)) == b)
                assert_(all(0 <= x <= ub for x, ub in zip(X, U)))
                feasible_count += 1

def test_diophantine_overflow():
    # Smoke test integer overflow detection
    max_intp = np.iinfo(np.intp).max
    max_int64 = np.iinfo(np.int64).max

    if max_int64 <= max_intp:
        # Check that the algorithm works internally in 128-bit;
        # solving this problem requires large intermediate numbers
        A = (max_int64//2, max_int64//2 - 10)
        U = (max_int64//2, max_int64//2 - 10)
        b = 2*(max_int64//2) - 10

        assert_equal(solve_diophantine(A, U, b), (1, 1))

def check_may_share_memory_exact(a, b):
    got = np.may_share_memory(a, b, max_work=MAY_SHARE_EXACT)

    assert_equal(np.may_share_memory(a, b),
                 np.may_share_memory(a, b, max_work=MAY_SHARE_BOUNDS))

    a.fill(0)
    b.fill(0)
    a.fill(1)
    exact = b.any()

    err_msg = ""
    if got != exact:
        err_msg = "    " + "\n    ".join([
            "base_a - base_b = %r" % (a.__array_interface__['data'][0] - b.__array_interface__['data'][0],),
            "shape_a = %r" % (a.shape,),
            "shape_b = %r" % (b.shape,),
            "strides_a = %r" % (a.strides,),
            "strides_b = %r" % (b.strides,),
            "size_a = %r" % (a.size,),
            "size_b = %r" % (b.size,)
        ])

    assert_equal(got, exact, err_msg=err_msg)

def test_may_share_memory_manual():
    # Manual test cases for may_share_memory

    # Base arrays
    xs0 = [
        np.zeros([13, 21, 23, 22], dtype=np.int8),
        np.zeros([13, 21, 23*2, 22], dtype=np.int8)[:,:,::2,:]
    ]

    # Generate all negative stride combinations
    xs = []
    for x in xs0:
        for ss in itertools.product(*(([slice(None), slice(None, None, -1)],)*4)):
            xp = x[ss]
            xs.append(xp)

    for x in xs:
        # The default is a simple extent check
        assert_(np.may_share_memory(x[:,0,:], x[:,1,:]))
        assert_(np.may_share_memory(x[:,0,:], x[:,1,:], max_work=None))

        # Exact checks
        check_may_share_memory_exact(x[:,0,:], x[:,1,:])
        check_may_share_memory_exact(x[:,::7], x[:,3::3])

        try:
            xp = x.ravel()
            if xp.flags.owndata:
                continue
            xp = xp.view(np.int16)
        except ValueError:
            continue

        # 0-size arrays cannot overlap
        check_may_share_memory_exact(x.ravel()[6:6],
                                     xp.reshape(13, 21, 23, 11)[:,::7])

        # Test itemsize is dealt with
        check_may_share_memory_exact(x[:,::7],
                                     xp.reshape(13, 21, 23, 11))
        check_may_share_memory_exact(x[:,::7],
                                     xp.reshape(13, 21, 23, 11)[:,3::3])
        check_may_share_memory_exact(x.ravel()[6:7],
                                     xp.reshape(13, 21, 23, 11)[:,::7])

    # Check unit size
    x = np.zeros([1], dtype=np.int8)
    check_may_share_memory_exact(x, x)
    check_may_share_memory_exact(x, x.copy())

def check_may_share_memory_easy_fuzz(get_max_work, same_steps, min_count):
    # Check that overlap problems with common strides are solved with
    # little work.
    x = np.zeros([17,34,71,97], dtype=np.int16)

    rng = np.random.RandomState(1234)

    def random_slice(n, step):
        start = rng.randint(0, n+1)
        stop = rng.randint(start, n+1)
        if rng.randint(0, 2) == 0:
            stop, start = start, stop
            step *= -1
        return slice(start, stop, step)

    feasible = 0
    infeasible = 0

    while min(feasible, infeasible) < min_count:
        steps = tuple(rng.randint(1, 11) if rng.randint(0, 5) == 0 else 1
                      for j in range(x.ndim))
        if same_steps:
            steps2 = steps
        else:
            steps2 = tuple(rng.randint(1, 11) if rng.randint(0, 5) == 0 else 1
                           for j in range(x.ndim))

        t1 = np.arange(x.ndim)
        rng.shuffle(t1)

        t2 = np.arange(x.ndim)
        rng.shuffle(t2)

        s1 = tuple(random_slice(p, s) for p, s in zip(x.shape, steps))
        s2 = tuple(random_slice(p, s) for p, s in zip(x.shape, steps2))
        a = x[s1].transpose(t1)
        b = x[s2].transpose(t2)

        bounds_overlap = np.may_share_memory(a, b)
        may_share_answer = np.may_share_memory(a, b)
        easy_answer = np.may_share_memory(a, b, max_work=get_max_work(a, b))
        exact_answer = np.may_share_memory(a, b, max_work=MAY_SHARE_EXACT)

        if easy_answer != exact_answer:
            # assert_equal is slow...
            assert_equal(easy_answer, exact_answer, err_msg=repr((s1, s2)))

        if may_share_answer != bounds_overlap:
            assert_equal(may_share_answer, bounds_overlap,
                         err_msg=repr((s1, s2)))

        if bounds_overlap:
            if exact_answer:
                feasible += 1
            else:
                infeasible += 1

def test_may_share_memory_easy_fuzz():
    # Check that overlap problems with common strides are always
    # solved with little work.

    check_may_share_memory_easy_fuzz(get_max_work=lambda a, b: 1,
                                     same_steps=True,
                                     min_count=2000)

def test_may_share_memory_harder_fuzz():
    # Overlap problems with not necessarily common strides take more
    # work.
    #
    # The work bound below can't be reduced much. Harder problems can
    # also exist but not be detected here, as the set of problems
    # comes from RNG.

    check_may_share_memory_easy_fuzz(get_max_work=lambda a, b: max(a.size, b.size)//2,
                                     same_steps=False,
                                     min_count=2000)

def test_shares_memory_api():
    x = np.zeros([4, 5, 6], dtype=np.int8)

    assert_equal(np.shares_memory(x, x), True)
    assert_equal(np.shares_memory(x, x.copy()), False)

    a = x[:,::2,::3]
    b = x[:,::3,::2]
    assert_equal(np.shares_memory(a, b), True)
    assert_equal(np.shares_memory(a, b, max_work=None), True)
    assert_raises(np.TooHardError, np.shares_memory, a, b, max_work=1)
    assert_raises(np.TooHardError, np.shares_memory, a, b, max_work=long(1))

def test_internal_overlap_diophantine():
    def check(A, U, exists=None):
        X = solve_diophantine(A, U, 0, require_ub_nontrivial=1)

        if exists is None:
            exists = (X is not None)

        if X is not None:
            assert_(sum(a*x for a, x in zip(A, X)) == sum(a*u//2 for a, u in zip(A, U)))
            assert_(all(0 <= x <= u for x, u in zip(X, U)))
            assert_(any(x != u//2 for x, u in zip(X, U)))

        if exists:
            assert_(X is not None, repr(X))
        else:
            assert_(X is None, repr(X))

    # Smoke tests
    check((3, 2), (2*2, 3*2), exists=True)
    check((3*2, 2), (15*2, (3-1)*2), exists=False)

def test_internal_overlap_slices():
    # Slicing an array never generates internal overlap

    x = np.zeros([17,34,71,97], dtype=np.int16)

    rng = np.random.RandomState(1234)

    def random_slice(n, step):
        start = rng.randint(0, n+1)
        stop = rng.randint(start, n+1)
        if rng.randint(0, 2) == 0:
            stop, start = start, stop
            step *= -1
        return slice(start, stop, step)

    cases = 0
    min_count = 5000

    while cases < min_count:
        steps = tuple(rng.randint(1, 11) if rng.randint(0, 5) == 0 else 1
                      for j in range(x.ndim))
        t1 = np.arange(x.ndim)
        rng.shuffle(t1)
        s1 = tuple(random_slice(p, s) for p, s in zip(x.shape, steps))
        a = x[s1].transpose(t1)

        assert not internal_overlap(a)
        cases += 1

def check_internal_overlap(a, manual_expected=None):
    got = internal_overlap(a)

    # Brute-force check
    m = set()
    ranges = tuple(xrange(n) for n in a.shape)
    for v in itertools.product(*ranges):
        offset = sum(s*w for s, w in zip(a.strides, v))
        if offset in m:
            expected = True
            break
        else:
            m.add(offset)
    else:
        expected = False

    # Compare
    if got != expected:
        assert_equal(got, expected, err_msg=repr((a.strides, a.shape)))
    if manual_expected is not None and expected != manual_expected:
        assert_equal(expected, manual_expected)
    return got

def test_internal_overlap_manual():
    # Stride tricks can construct arrays with internal overlap

    # We don't care about memory bounds, the array is not
    # read/write accessed
    x = np.arange(1).astype(np.int8)

    # Check low-dimensional special cases

    check_internal_overlap(x, False) # 1-dim
    check_internal_overlap(x.reshape([]), False) # 0-dim

    a = as_strided(x, strides=(3, 4), shape=(4, 4))
    check_internal_overlap(a, False)

    a = as_strided(x, strides=(3, 4), shape=(5, 4))
    check_internal_overlap(a, True)

    a = as_strided(x, strides=(0,), shape=(0,))
    check_internal_overlap(a, False)

    a = as_strided(x, strides=(0,), shape=(1,))
    check_internal_overlap(a, False)

    a = as_strided(x, strides=(0,), shape=(2,))
    check_internal_overlap(a, True)

    a = as_strided(x, strides=(0, -9993), shape=(87, 22))
    check_internal_overlap(a, True)

    a = as_strided(x, strides=(0, -9993), shape=(1, 22))
    check_internal_overlap(a, False)

    a = as_strided(x, strides=(0, -9993), shape=(0, 22))
    check_internal_overlap(a, False)

def test_internal_overlap_fuzz():
    # Fuzz check; the brute-force check is fairly slow

    x = np.arange(1).astype(np.int8)

    overlap = 0
    no_overlap = 0
    min_count = 100

    rng = np.random.RandomState(1234)

    while min(overlap, no_overlap) < min_count:
        ndim = rng.randint(1, 4)

        strides = tuple(rng.randint(-1000, 1000) for j in range(ndim))
        shape = tuple(rng.randint(1, 30) for j in range(ndim))

        a = as_strided(x, strides=strides, shape=shape)
        result = check_internal_overlap(a)

        if result:
            overlap += 1
        else:
            no_overlap += 1

if __name__ == "__main__":
    run_module_suite()

max_stars_repo_name: leifdenby/numpy

max_issues_repo_name: leifdenby/numpy

max_forks_repo_name: leifdenby/numpy