text:
"""
REDS dataset
support reading images from lmdb, image folder and memcached
"""
import os.path as osp
import random
import pickle
import logging
import numpy as np
import cv2
import lmdb
import torch
import torch.utils.data as data
import data.util as util

try:
    import mc  # import memcached
except ImportError:
    pass

logger = logging.getLogger("base")
imgShape = (3, 128, 128)
lrShape = (3, 32, 32)

class LicensePlateDataset(data.Dataset):
    """
    Reading the training LicensePlateDataset dataset
    key example: XXXX_X_Y_ZZZZZZ,
    X: Directory
    Y: Total number of images
    Z: Image name
    GT: Ground-Truth;
    LQ: Low-Quality, e.g., low-resolution/blurry/noisy/compressed frames
    support reading N LQ frames, N = 1, 3, 5, 7
    """

    def __init__(self, opt, isVal=False):
        super(LicensePlateDataset, self).__init__()
        self.opt = opt
        # temporal augmentation
        self.interval_list = opt["interval_list"]
        self.random_reverse = opt["random_reverse"]
        logger.info(
            "Temporal augmentation interval list: [{}], with random reverse is {}.".format(
                ",".join(str(x) for x in opt["interval_list"]), self.random_reverse
            )
        )

        self.half_N_frames = opt["N_frames"] // 2
        self.GT_root, self.LQ_root = opt["dataroot_GT"], opt["dataroot_LQ"]
        self.data_type = self.opt["data_type"]
        self.LR_input = (
            False if opt["GT_size"] == opt["LQ_size"] else True
        )  # low resolution inputs
        #### directly load image keys

        # Check if this is validation dataset
        if isVal:
            logger.info("Using cache keys for validating: {}".format(opt["cache_keys"]))
        else:
            logger.info("Using cache keys for training: {}".format(opt["cache_keys"]))

        cache_keys_path = opt["cache_keys"]

        if cache_keys_path:
            logger.info("Using cache keys - {}.".format(cache_keys_path))
            self.paths_GT = pickle.load(open(cache_keys_path, "rb"))
        # remove the REDS4 for testing
        #         self.paths_GT = [
        #             v
        #             for v in self.paths_GT
        #             if v.split("_")[0] not in ["000", "011", "015", "020"]
        #         ]
        #         assert self.paths_GT, "Error: GT path is empty."

        if self.data_type == "lmdb":
            self.GT_env, self.LQ_env = None, None
        elif self.data_type == "mc":  # memcached
            self.mclient = None
        elif self.data_type == "img":
            pass
        else:
            raise ValueError("Wrong data type: {}".format(self.data_type))

    def _init_lmdb(self):
        # https://github.com/chainer/chainermn/issues/129
        self.GT_env = lmdb.open(
            self.opt["dataroot_GT"],
            readonly=True,
            lock=False,
            readahead=False,
            meminit=False,
        )
        self.LQ_env = lmdb.open(
            self.opt["dataroot_LQ"],
            readonly=True,
            lock=False,
            readahead=False,
            meminit=False,
        )

    def _ensure_memcached(self):
        if self.mclient is None:
            # specify the config files
            server_list_config_file = None
            client_config_file = None
            self.mclient = mc.MemcachedClient.GetInstance(
                server_list_config_file, client_config_file
            )

    def _read_img_mc(self, path):
        """ Return BGR, HWC, [0, 255], uint8"""
        value = mc.pyvector()
        self.mclient.Get(path, value)
        value_buf = mc.ConvertBuffer(value)
        img_array = np.frombuffer(value_buf, np.uint8)
        img = cv2.imdecode(img_array, cv2.IMREAD_UNCHANGED)
        return img

    def _read_img_mc_BGR(self, path, name_a, name_b):
        """ Read BGR channels separately and then combine for 1M limits in cluster"""
        img_B = self._read_img_mc(osp.join(path + "_B", name_a, name_b + ".png"))
        img_G = self._read_img_mc(osp.join(path + "_G", name_a, name_b + ".png"))
        img_R = self._read_img_mc(osp.join(path + "_R", name_a, name_b + ".png"))
        img = cv2.merge((img_B, img_G, img_R))
        return img

    def __getitem__(self, index):
        if self.data_type == "mc":
            self._ensure_memcached()
        elif self.data_type == "lmdb":
            if (self.GT_env is None) or (self.LQ_env is None):
                self._init_lmdb()

        scale = self.opt["scale"]
        GT_size = self.opt["GT_size"]
        key = self.paths_GT[index]
        name_a, totalKeyPerDir, name_b = key.rsplit("_", 2)

        # Restore the origin key to match file name
        key = name_a + "_" + name_b
        totalKeyPerDir = int(totalKeyPerDir)
        center_frame_idx = int(name_b)

        #### determine the neighbor frames
        interval = random.choice(self.interval_list)
        if self.opt["border_mode"]:
            direction = 1  # 1: forward; 0: backward
            N_frames = self.opt["N_frames"]
            if self.random_reverse and random.random() < 0.5:
                direction = random.choice([0, 1])
            if center_frame_idx + interval * (N_frames - 1) > totalKeyPerDir:
                direction = 0
            elif center_frame_idx - interval * (N_frames - 1) < 0:
                direction = 1
            # get the neighbor list
            if direction == 1:
                neighbor_list = list(
                    range(
                        center_frame_idx,
                        center_frame_idx + interval * N_frames,
                        interval,
                    )
                )
            else:
                neighbor_list = list(
                    range(
                        center_frame_idx,
                        center_frame_idx - interval * N_frames,
                        -interval,
                    )
                )
            name_b = "{:06d}".format(neighbor_list[0])
        else:
            # ensure not exceeding the borders
            while (
                center_frame_idx + self.half_N_frames * interval > totalKeyPerDir
            ) or (center_frame_idx - self.half_N_frames * interval < 1):
                center_frame_idx = random.randint(1, totalKeyPerDir)
            # get the neighbor list
            neighbor_list = list(
                range(
                    center_frame_idx - self.half_N_frames * interval,
                    center_frame_idx + self.half_N_frames * interval + 1,
                    interval,
                )
            )
            if self.random_reverse and random.random() < 0.5:
                neighbor_list.reverse()
            name_b = "{:06d}".format(neighbor_list[self.half_N_frames])

        assert (
            len(neighbor_list) == self.opt["N_frames"]
        ), "Wrong length of neighbor list: {}".format(len(neighbor_list))

        #### get the GT image (as the center frame)
        if self.data_type == "mc":
            img_GT = self._read_img_mc_BGR(self.GT_root, name_a, name_b)
            img_GT = img_GT.astype(np.float32) / 255.0
        elif self.data_type == "lmdb":
            img_GT = util.read_img(self.GT_env, key, imgShape)
        else:
            img_GT = util.read_img(
                None, osp.join(self.GT_root, name_a, name_b + ".png")
            )

        #### get LQ images
        # LQ_size_tuple = (3, 32, 32) if self.LR_input else (3, 128, 128)
        LQ_size_tuple = lrShape
        img_LQ_l = []
        for v in neighbor_list:
            img_LQ_path = osp.join(self.LQ_root, name_a, "{:06d}.png".format(v))
            if self.data_type == "mc":
                if self.LR_input:
                    img_LQ = self._read_img_mc(img_LQ_path)
                else:
                    img_LQ = self._read_img_mc_BGR(
                        self.LQ_root, name_a, "{:06d}".format(v)
                    )
                img_LQ = img_LQ.astype(np.float32) / 255.0
            elif self.data_type == "lmdb":
                # Log info if the size is not matched
                if LQ_size_tuple != lrShape:
                    print(LQ_size_tuple, name_a, v)

                img_LQ = util.read_img(
                    self.LQ_env, "{}_{:06d}".format(name_a, v), LQ_size_tuple
                )
            else:
                img_LQ = util.read_img(None, img_LQ_path)
            img_LQ_l.append(img_LQ)

        if self.opt["phase"] == "train":
            C, H, W = LQ_size_tuple  # LQ size
            # randomly crop
            if self.LR_input:
                LQ_size = GT_size // scale
                rnd_h = random.randint(0, max(0, H - LQ_size))
                rnd_w = random.randint(0, max(0, W - LQ_size))
                img_LQ_l = [
                    v[rnd_h : rnd_h + LQ_size, rnd_w : rnd_w + LQ_size, :]
                    for v in img_LQ_l
                ]
                rnd_h_HR, rnd_w_HR = int(rnd_h * scale), int(rnd_w * scale)
                img_GT = img_GT[
                    rnd_h_HR : rnd_h_HR + GT_size, rnd_w_HR : rnd_w_HR + GT_size, :
                ]
            else:
                rnd_h = random.randint(0, max(0, H - GT_size))
                rnd_w = random.randint(0, max(0, W - GT_size))
                img_LQ_l = [
                    v[rnd_h : rnd_h + GT_size, rnd_w : rnd_w + GT_size, :]
                    for v in img_LQ_l
                ]
                img_GT = img_GT[rnd_h : rnd_h + GT_size, rnd_w : rnd_w + GT_size, :]

            # augmentation - flip, rotate
            img_LQ_l.append(img_GT)
            rlt = util.augment(img_LQ_l, self.opt["use_flip"], self.opt["use_rot"])
            img_LQ_l = rlt[0:-1]
            img_GT = rlt[-1]

        # stack LQ images to NHWC, N is the frame number
        img_LQs = np.stack(img_LQ_l, axis=0)
        # BGR to RGB, HWC to CHW, numpy to tensor
        img_GT = img_GT[:, :, [2, 1, 0]]
        img_LQs = img_LQs[:, :, :, [2, 1, 0]]
        img_GT = torch.from_numpy(
            np.ascontiguousarray(np.transpose(img_GT, (2, 0, 1)))
        ).float()
        img_LQs = torch.from_numpy(
            np.ascontiguousarray(np.transpose(img_LQs, (0, 3, 1, 2)))
        ).float()
        return {"LQs": img_LQs, "GT": img_GT, "key": key}

    def __len__(self):
        return len(self.paths_GT)


text:
"""
Validate calculation bfield calculation from line segments
=================================================================
"""

import numpy as np
import matplotlib.pyplot as plt
import time as t
from bfieldtools.line_magnetics import magnetic_field

""" Bfield calculation from circular current loops using elliptic integrals
"""

def double_factorial(n):
    if n <= 0:
        return 1
    else:
        return n * double_factorial(n - 2)

def bfield_iterative(x, y, z, x_c, y_c, z_c, r, I, n):
    """Compute b field of a current loop using an iterative method to estimate
    elliptic integrals.

    Parameters:
        x, y, z: Evaluation points in 3D space. Accepts matrices and integers.
        x_c, y_c, z_c: Coordinates of the center point of the current loop.
        r: Radius of the current loop.
        I: Current of the current loop.
        n: Number of terms in the serie expansion.

    Returns:
        bfiels (N_points, 3) at evaluation points.

    This calculation is based on paper by Robert A. Schill, Jr (General
    Relation for the Vector Magnetic Field of aCircular Current Loop:
    A Closer Look). DOI: 10.1109/TMAG.2003.808597
    """

    st2 = t.time()
    np.seterr(divide="ignore", invalid="ignore")
    u0 = 4 * np.pi * 1e-7
    Y = u0 * I / (2 * np.pi)

    # Change to cylideric coordinates
    rc = np.sqrt(np.power(x - x_c, 2) + np.power(y - y_c, 2))

    # Coefficients for estimating elliptic integrals with nth degree series
    # expansion using Legendre polynomials
    m = 4 * r * rc / (np.power((rc + r), 2) + np.power((z - z_c), 2))
    K = 1
    E = 1
    for i in range(1, n + 1):
        K = K + np.square(
            double_factorial(2 * i - 1) / double_factorial(2 * i)
        ) * np.power(m, i)
        E = E - np.square(
            double_factorial(2 * i - 1) / double_factorial(2 * i)
        ) * np.power(m, i) / (2 * i - 1)

    K = K * np.pi / 2
    E = E * np.pi / 2

    # Calculation of radial and axial components of B-field
    Brc = (
        Y
        * (z - z_c)
        / (rc * np.sqrt(np.power((rc + r), 2) + np.power((z - z_c), 2)))
        * (
            -K
            + E
            * (np.power(rc, 2) + np.power(r, 2) + np.power((z - z_c), 2))
            / (np.power((rc - r), 2) + np.power((z - z_c), 2))
        )
    )

    Bz = (
        Y
        / (np.sqrt(np.power((rc + r), 2) + np.power((z - z_c), 2)))
        * (
            K
            - E
            * (np.power(rc, 2) - np.power(r, 2) + np.power((z - z_c), 2))
            / (np.power((rc - r), 2) + np.power((z - z_c), 2))
        )
    )

    # Set nan and inf values to 0
    Brc[np.isinf(Brc)] = 0
    Bz[np.isnan(Bz)] = 0
    Brc[np.isnan(Brc)] = 0
    Bz[np.isinf(Bz)] = 0

    # Change back to cartesian coordinates
    Bx = Brc * (x - x_c) / rc
    By = Brc * (y - y_c) / rc

    # Change nan values from coordinate transfer to 0
    Bx[np.isnan(Bx)] = 0
    By[np.isnan(By)] = 0

    B = np.zeros((3, X.size), dtype=np.float64)
    B[0] = Bx.flatten()
    B[1] = By.flatten()
    B[2] = Bz.flatten()

    et2 = t.time()
    print("Execution time for iterative method is:", et2 - st2)

    return B.T

""" Plot field of a circular current path
"""
x = np.linspace(-1, 1, 100)
Ntheta = 10000
theta = np.linspace(0, 2 * np.pi, Ntheta)
vertices = np.zeros((Ntheta, 3), dtype=np.float64)
vertices[:, 0] = np.cos(theta) * 0.1
vertices[:, 1] = np.sin(theta) * 0.1
vertices[:, 2] = 0.2

X, Y = np.meshgrid(x, x, indexing="ij")
Z = np.zeros((x.size, x.size))

points = np.zeros((3, X.size), dtype=np.float64)
points[0] = X.flatten()
points[1] = Y.flatten()

b1 = magnetic_field(vertices, points.T)  # Calculates discretised bfield
b2 = bfield_iterative(X, Y, Z, 0, 0, 0.2, 0.1, 1, 25)  # Calculates bfield iteratively

# Error between two calculation methods.
berr = (b2 - b1) / b1 * 100
BE = berr.T[2]  # By changing the index, errors in different components can be obtained
ind = np.where(np.abs(BE) > 0.1)  # The limit for significant error is set to 0.1%
bpoints = points.T[ind]

from mayavi import mlab

mlab.figure(1)
q = mlab.quiver3d(*points, *b1.T)
q.glyph.glyph_source.glyph_position = "center"
mlab.plot3d(*vertices.T)

mlab.figure(2)
q = mlab.quiver3d(*points, *b2.T)
q.glyph.glyph_source.glyph_position = "center"
mlab.plot3d(*vertices.T)

plt.figure(3)
plt.hist(berr.T[2], bins=50, density=True, histtype="bar")
plt.title("Histogram of error between calculation methods.")
plt.xlabel("%")

#%% Plot the b-field vectors exceeding the error limit
if len(bpoints > 0):

    from mayavi import mlab

    mlab.figure(3)
    q = mlab.quiver3d(*bpoints.T, *b1[ind].T)
    q.glyph.glyph_source.glyph_position = "center"
    mlab.plot3d(*vertices.T)

    q = mlab.quiver3d(*bpoints.T, *b2[ind].T)
    q.glyph.glyph_source.glyph_position = "center"

ext: py

max_stars_repo_name: SAFedorov/bfieldtools

max_issues_repo_name: SAFedorov/bfieldtools

max_forks_repo_name: SAFedorov/bfieldtools

reason: import numpy


text:
from contextlib import contextmanager
import math
import os
import random
import psutil
import sys
import time

import numpy as np
import torch

def seed_torch(seed=None, random_seed=True):
    if random_seed or seed is None:
        seed = np.random.randint(0, 1000000)
    random.seed(seed)
    os.environ["PYTHONHASHSEED"] = str(seed)
    np.random.seed(seed)
    torch.manual_seed(seed)
    torch.cuda.manual_seed_all(seed)
    torch.backends.cudnn.deterministic = True

@contextmanager
def trace(title, logger=None):
    t0 = time.time()
    p = psutil.Process(os.getpid())
    m0 = p.memory_info()[0] / 2. ** 30
    yield
    m1 = p.memory_info()[0] / 2. ** 30
    delta = m1 - m0
    sign = '+' if delta >= 0 else '-'
    delta = math.fabs(delta)
    message = f'[{m1:.1f}GB({sign}{delta:.1f}GB):{time.time() - t0:.1f}sec] {title} '
    print(message)
    if logger is not None:
        logger.info(message)


text:
import numpy as np
from PIL import Image
from usuelles import *
from TestPerlin import *

def SquareDiamond(taille, ecart, poids, nb):
    """Retourne une matrice de taille * taille en utilisant l'algorithme Square-Diamond avec un écart-type ecart et un poids poids sur le coeff rando. nb est le nombre d'itérations avec coeff aléatoire, le reste étant juste de l'interpolation"""
    matrice = np.zeros((taille, taille))
    matrice[0, 0] = poids * np.random.normal(scale = ecart)
    matrice[0, taille - 1] = poids * np.random.normal(scale = ecart)
    matrice[taille - 1, 0] = poids * np.random.normal(scale = ecart)
    matrice[taille - 1, taille - 1] = poids * np.random.normal(scale = ecart)
    padding = taille - 1
    while padding > 1:
        pas = padding // 2
        print(padding, pas)
        if padding % 1 != 0:
            print("Problème de dimension")
            return 0
        #Diamond
        for x in range(pas, taille, padding):
            for y in range(pas, taille, padding):
                if nb > 0 :
                    matrice[x, y] = (matrice[x - pas, y - pas] + matrice[x - pas, y + pas] + matrice[x + pas, y - pas] + matrice[x + pas, y + pas]) / 4 + poids * np.random.normal(scale = ecart)
                else :
                    matrice[x, y] = (matrice[x - pas, y - pas] + matrice[x - pas, y + pas] + matrice[x + pas, y - pas] + matrice[x + pas, y + pas]) / 4
        decalage = 0
        for x in range(0, taille, pas):
            if decalage == 0:
                decalage = pas
            else:
                decalage = 0
            for y in range(decalage, taille, padding):
                tot = 0
                n = 0
                if x >= pas:
                    tot += matrice[x - pas, y]
                    n += 1
                if x + pas < taille:
                    tot += matrice[x + pas, y]
                    n += 1
                if y >= pas:
                    tot += matrice[x, y - pas]
                    n += 1
                if y + pas < taille:
                    tot += matrice[x, y + pas]
                if nb > 0 :
                    matrice[x, y] = tot / n + poids * np.random.normal(scale = ecart)
                else :
                    matrice[x, y] = tot / n
        padding = pas
        nb -= 1
    return matrice

def printMatGrad(mat, couleurs = "gradient.png"):
    """Retourne la matrice de couleurs avec le gradient en entrée."""
    (matrice, mini, maxi) = mat, np.min(mat), np.max(mat)
    print(mini, maxi)
    n = len(matrice)
    image = Image.new('RGB', (n, n))
    gradient = Image.open(couleurs)
    lGrad = gradient.size[0]
    for x in range(n):
        print(x)
        for y in range(n):
            image.putpixel((x, y), gradient.getpixel( (int( map(matrice[x, y], (mini, maxi), (0, lGrad - 1)) ), 1)) )
    nom = str(np.random.randint(0, 99999))
    print(nom)
    image.save(nom + '.png')
    return True

def interpMatrix(matrice, Taillei, Taillef, methode = 'linear'):
    """Faire l'interpolation de la matrice."""
    if methode == 'linear' :
        interp = lambda t : t
    elif methode == 'polynome':
        interp = lambda t : 6 * (t ** 5) - 12 * (t ** 4) + 10 * (t ** 3)
        #interp = lambda t : t ** 2
    else:
        print('Erreur de méthode')
        return 0
    nombre = (Taillef - Taillei) / (Taillei - 1)
    if int(nombre) != nombre:
        print('Erreur de format')
        return 0
    nombre = int(nombre)
    facteur = 1 + nombre
    print(nombre, facteur)
    res = np.zeros((Taillef, Taillef))
    for x in range(Taillef):
        print('interp' + str(x))
        for y in range(Taillef):
            dx = x % facteur
            dy = y % facteur
            if dx == 0 and dy == 0:
                res[x, y] = matrice[x // facteur, y // facteur]
            else:
                qx = x // facteur
                qy = y // facteur
                cx = dx / facteur
                cy = dy / facteur
                ex = 1 - cx
                ey = 1 - cy
                if dx == 0:
                    res[x, y] = matrice[qx, qy] * interp(ey) + matrice[qx, qy + 1] * interp(cy)
                elif dy == 0:
                    res[x, y] = matrice[qx, qy] * interp(ex) + matrice[qx + 1, qy] * interp(cx)
                else:
                    res[x, y] = matrice[qx, qy] * (interp(ex * ey)) + matrice[qx + 1, qy + 1] * (interp(cx * cy)) + matrice[qx, qy + 1] * (interp(ex * cy)) + matrice[qx + 1, qy] * (interp(cx * ey))
    #Remap sur [0, 1]
    mini = np.min(res)
    maxi = np.max(res)
    res = (res - mini) * ( 1 - 0) / (maxi - mini) + 0
    return res

def taillemini(nombre):
    return int(np.log2(nombre)) + 1

def ajouterPerlin(matrice, octaves, pers, poids):
    """Ajoute à la matrice un bruit de Perlin."""
    Taille = 2 ** taillemini(matrice.shape[0])
    bruit = generate_octaves((Taille, Taille), 2, pers, octaves)
    res = bruit[0:matrice.shape[0], 0:matrice.shape[0]]
    return matrice + (poids * res)


text:
import pyfx
import numpy as np
import pytest

def test_position_by_string():

    test_array_size = 11

    # Test array check
    test_array = np.zeros((test_array_size))
    with pytest.raises(ValueError):
        pyfx.util.helper.position_by_string(test_array.shape)

    # Test default (middle center)
    test_array = np.zeros((test_array_size,test_array_size))
    r, c = pyfx.util.helper.position_by_string(test_array.shape)
    assert r == test_array_size//2
    assert c == test_array_size//2

    # Test bad entry (middle middle)
    test_array = np.zeros((test_array_size,test_array_size))
    with pytest.raises(ValueError):
        pyfx.util.helper.position_by_string(test_array.shape,"middlemiddle")

    r = ["","top","middle","bottom"]
    r_value = [test_array_size//2,0,test_array_size//2,test_array_size-1]
    y_value = [test_array_size//2,test_array_size-1,test_array_size//2,0]

    c = ["","left","center","right"]
    c_value = [test_array_size//2,0,test_array_size//2,test_array_size-1]
    x_value = [test_array_size//2,0,test_array_size//2,test_array_size-1]

    for i in range(len(r)):
        for j in range(len(c)):
            pos_str = "{}{}".format(r[i],c[j])
            r_out, c_out = pyfx.util.helper.position_by_string(test_array.shape,
                                                               pos_str)
            assert r_out == r_value[i]
            assert c_out == c_value[j]

            x_out, y_out = pyfx.util.helper.position_by_string(test_array.shape,
                                                               pos_str,
                                                               return_xy=True)
            assert x_out == x_value[j]
            assert y_out == y_value[i]

max_stars_repo_name: harmsm/pyfx

max_issues_repo_name: harmsm/pyfx

max_forks_repo_name: harmsm/pyfx


text:
import numpy as np
import scipy.io
from scipy.io import loadmat
from sklearn.cross_validation import train_test_split
from sklearn.metrics import accuracy_score
from sklearn.metrics import confusion_matrix
import warnings
import time

from sklearn.neighbors import KNeighborsClassifier
from sklearn.tree import DecisionTreeClassifier
from sklearn.svm import SVC
from sklearn.ensemble import RandomForestClassifier, AdaBoostClassifier
from sklearn.naive_bayes import GaussianNB
from sklearn.linear_model import LogisticRegression

names = ["Nearest Neighbors", "Linear SVM", "Gaussian SVM", "Polynomial SVM", "Decision Tree",
		 "Random Forest", "AdaBoost Classifier", "Logistic Regression", "Naive Bayes"]

classifiers = [
	KNeighborsClassifier(n_neighbors=15, weights='distance'),
	SVC(kernel="linear", C=3.4),
	SVC(kernel="rbf", C=3.4, gamma=0.1),
	SVC(kernel="poly", C=3.4, degree=2, gamma=0.1),
	DecisionTreeClassifier(),
	RandomForestClassifier(n_estimators=300, n_jobs=-1),
	AdaBoostClassifier(n_estimators=70),
	LogisticRegression(random_state=1, C=0.4),
	GaussianNB()]

def main():
	#ignore all warnings
	#warnings.filterwarnings("ignore")
	start = time.time()

	#load the data
	trainX = np.load('trainX_feat.npy')
	testX = np.load('testX_feat.npy')
	trainY = np.load('trainY_feat.npy')
	testY = np.load('testY_feat.npy')
	print('\n!!! Data Loading Completed !!!\n')

	for name, clf in zip(names, classifiers):
		clf.fit(trainX, trainY)
		pred = clf.predict(testX)
		print('classifier : %s, Accuracy :%f%% ' %(name, accuracy_score(testY, pred)*100))
		print('total : %d, correct : %d, incorrect : %d\n' %(len(pred), np.sum(pred == testY), np.sum(pred != testY)))

		#print(confusion_matrix(testY, pred)) #uncomment to get confusion matrix
	print('Test Time : %f Minutes\n' %((time.time()-start)/60))

if __name__ == '__main__':
	main()

max_stars_repo_name: akhilpm/Masters-Project

max_issues_repo_name: akhilpm/Masters-Project

max_forks_repo_name: akhilpm/Masters-Project

reason: import numpy,import scipy,from scipy


text:
import shutil
from collections import defaultdict
import numpy as np
import pandas as pd
import torch
import torch.nn.functional as F
from sklearn.model_selection import train_test_split
from torch import nn
from torch.utils.data import Dataset, DataLoader
from transformers import (
    BertModel,
    BertTokenizer,
    AdamW
)
import argparse
import os
from tqdm import tqdm
import requests
from torchtext.utils import download_from_url, extract_archive
from torchtext.datasets.text_classification import URLS
import sys
import argparse
import logging
import pyarrow as pa
import pyarrow.csv as pv
import pyarrow.parquet as pq
from pathlib import Path
import subprocess

if __name__ == "__main__":

    import sys
    import json

    data_set = json.loads(sys.argv[1])[0]['dataset_url']
    output_path = json.loads(sys.argv[2])[0]

    Path(output_path).mkdir(parents=True, exist_ok=True)

    dataset_tar = download_from_url(
            data_set, root="./")
    extracted_files = extract_archive(dataset_tar)

    ag_news_csv = pv.read_csv("ag_news_csv/train.csv")

    pq.write_table(ag_news_csv, os.path.join(output_path, "ag_news_data.parquet"))

    entry_point=["ls", "-R", output_path]
    run_code = subprocess.run(entry_point, stdout=subprocess.PIPE)
    print(run_code.stdout)

    # for root, dirs, files in os.walk(".", topdown = False):
    #     for name in files:
    #         print(os.path.join(root, name))
    #     for name in dirs:
    #         print(os.path.join(root, name))


text:
# python 3.6
"""diffuses target images to context images with In-domain GAN Inversion.

Basically, this script first copies the central region from the target image to
the context image, and then performs in-domain GAN inversion on the stitched
image. Different from `intert.py`, masked reconstruction loss is used in the
optimization stage.

NOTE: This script will diffuse every image from `target_image_list` to every
image from `context_image_list`.
"""

import os
import argparse
from tqdm import tqdm
import numpy as np

from utils.inverter import StyleGANInverter
from utils.logger import setup_logger
from utils.visualizer import HtmlPageVisualizer
from utils.visualizer import load_image, resize_image

def parse_args():
  """Parses arguments."""
  parser = argparse.ArgumentParser()
  parser.add_argument('model_name', type=str, help='Name of the GAN model.')
  parser.add_argument('target_list', type=str,
                      help='List of target images to diffuse from.')
  parser.add_argument('context_list', type=str,
                      help='List of context images to diffuse to.')
  parser.add_argument('-o', '--output_dir', type=str, default='',
                      help='Directory to save the results. If not specified, '
                           '`./results/diffusion` will be used by default.')
  parser.add_argument('-s', '--crop_size', type=int, default=110,
                      help='Crop size. (default: 110)')
  parser.add_argument('-x', '--center_x', type=int, default=125,
                      help='X-coordinate (column) of the center of the cropped '
                           'patch. This field should be adjusted according to '
                           'dataset and image size. (default: 125)')
  parser.add_argument('-y', '--center_y', type=int, default=145,
                      help='Y-coordinate (row) of the center of the cropped '
                           'patch. This field should be adjusted according to '
                           'dataset and image size. (default: 145)')
  parser.add_argument('--learning_rate', type=float, default=0.01,
                      help='Learning rate for optimization. (default: 0.01)')
  parser.add_argument('--num_iterations', type=int, default=100,
                      help='Number of optimization iterations. (default: 100)')
  parser.add_argument('--num_results', type=int, default=5,
                      help='Number of intermediate optimization results to '
                           'save for each sample. (default: 5)')
  parser.add_argument('--loss_weight_feat', type=float, default=5e-5,
                      help='The perceptual loss scale for optimization. '
                           '(default: 5e-5)')
  parser.add_argument('--viz_size', type=int, default=256,
                      help='Image size for visualization. (default: 256)')
  parser.add_argument('--gpu_id', type=str, default='0',
                      help='Which GPU(s) to use. (default: `0`)')
  parser.add_argument('--batch_size', type=int, default=4,
                      help='Batch size. (default: 4)')
  return parser.parse_args()

def main():
  """Main function."""
  args = parse_args()
  os.environ["CUDA_VISIBLE_DEVICES"] = args.gpu_id
  assert os.path.exists(args.target_list)
  target_list_name = os.path.splitext(os.path.basename(args.target_list))[0]
  assert os.path.exists(args.context_list)
  context_list_name = os.path.splitext(os.path.basename(args.context_list))[0]
  output_dir = args.output_dir or f'results/diffusion'
  job_name = f'{target_list_name}_TO_{context_list_name}'
  logger = setup_logger(output_dir, f'{job_name}.log', f'{job_name}_logger')

  logger.info(f'Loading model.')
  inverter = StyleGANInverter(
      args.model_name,
      learning_rate=args.learning_rate,
      iteration=args.num_iterations,
      reconstruction_loss_weight=1.0,
      perceptual_loss_weight=args.loss_weight_feat,
      regularization_loss_weight=0.0,
      logger=logger)
  image_size = inverter.G.resolution

  # Load image list.
  logger.info(f'Loading target images and context images.')
  target_list = []
  with open(args.target_list, 'r') as f:
    for line in f:
      target_list.append(line.strip())
  num_targets = len(target_list)
  context_list = []
  with open(args.context_list, 'r') as f:
    for line in f:
      context_list.append(line.strip())
  num_contexts = len(context_list)
  num_pairs = num_targets * num_contexts

  # Initialize visualizer.
  save_interval = args.num_iterations // args.num_results
  headers = ['Target Image', 'Context Image', 'Stitched Image',
             'Encoder Output']
  for step in range(1, args.num_iterations + 1):
    if step == args.num_iterations or step % save_interval == 0:
      headers.append(f'Step {step:06d}')
  viz_size = None if args.viz_size == 0 else args.viz_size
  visualizer = HtmlPageVisualizer(
      num_rows=num_pairs, num_cols=len(headers), viz_size=viz_size)
  visualizer.set_headers(headers)

  # Diffuse images.
  logger.info(f'Start diffusion.')
  latent_codes = []
  for target_idx in tqdm(range(num_targets), desc='Target ID', leave=False):
    # Load target.
    target_image = resize_image(load_image(target_list[target_idx]),
                                (image_size, image_size))
    visualizer.set_cell(target_idx * num_contexts, 0, image=target_image)
    for context_batch_idx in tqdm(range(0, num_contexts, args.batch_size),
                            desc='Context ID', leave=False):
      context_images = []
      for it in range(args.batch_size):
        context_idx = context_batch_idx + it
        if context_idx >= num_contexts:
          continue
        row_idx = target_idx * num_contexts + context_idx
        context_image = resize_image(load_image(context_list[context_idx]),
                                     (image_size, image_size))
        context_images.append(context_image)
        visualizer.set_cell(row_idx, 1, image=context_image)
      code, viz_results = inverter.easy_diffuse(target=target_image,
                                                context=np.asarray(context_images),
                                                center_x=args.center_x,
                                                center_y=args.center_y,
                                                crop_x=args.crop_size,
                                                crop_y=args.crop_size,
                                                num_viz=args.num_results)
      for key, values in viz_results.items():
        context_idx = context_batch_idx + key
        row_idx = target_idx * num_contexts + context_idx
        for viz_idx, viz_img in enumerate(values):
          visualizer.set_cell(row_idx, viz_idx + 2, image=viz_img)
      latent_codes.append(code)

  # Save results.
  os.system(f'cp {args.target_list} {output_dir}/target_list.txt')
  os.system(f'cp {args.context_list} {output_dir}/context_list.txt')
  np.save(f'{output_dir}/{job_name}_inverted_codes.npy',
          np.concatenate(latent_codes, axis=0))
  visualizer.save(f'{output_dir}/{job_name}.html')

if __name__ == '__main__':
  main()


text:
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
#
#   See COPYING file distributed along with the PyMVPA package for the
#   copyright and license terms.
#
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
"""Unit tests for PyMVPA pattern handling"""

import unittest
import numpy as np

from mvpa2.testing.datasets import datasets
from mvpa2.base.node import ChainNode
from mvpa2.datasets.base import dataset_wizard, Dataset
from mvpa2.generators.partition import NFoldPartitioner
from mvpa2.mappers.slicing import StripBoundariesSamples
from mvpa2.generators.splitters import Splitter
from mvpa2.generators.partition import *

from mvpa2.testing.tools import ok_, assert_array_equal, assert_true, \
        assert_false, assert_equal, assert_not_equal, reseed_rng

class SplitterTests(unittest.TestCase):

    @reseed_rng()
    def setUp(self):
        self.data = dataset_wizard(np.random.normal(size=(100,10)),
                            targets=[ i%4 for i in range(100) ],
                            chunks=[ i//10 for i in range(100)])

    def test_splitattr_deprecation(self):
        # Just a smoke test -- remove for 2.1 release
        nfs = NFoldPartitioner()
        _ = nfs.splitattr

    def test_reprs(self):
        # very very basic test to see that there is no errors in reprs
        # of partitioners
        import mvpa2.generators.partition as mgp
        for sclass in (x for x in dir(mgp) if x.endswith('Partitioner')):
            args = (1,)
            if sclass == 'ExcludeTargetsCombinationsPartitioner':
                args += (1,1)
            pclass = getattr(mgp, sclass)
            r = repr(pclass(*args))
            assert_false('ERROR' in r)

    def test_simplest_cv_pat_gen(self):
        # create the generator
        nfs = NFoldPartitioner(cvtype=1)
        spl = Splitter(attr='partitions')
        # now get the xval pattern sets One-Fold CV)
        xvpat = [ list(spl.generate(p)) for p in nfs.generate(self.data) ]

        self.assertTrue( len(xvpat) == 10 )

        for i,p in enumerate(xvpat):
            self.assertTrue( len(p) == 2 )
            self.assertTrue( p[0].nsamples == 90 )
            self.assertTrue( p[1].nsamples == 10 )
            self.assertTrue( p[1].chunks[0] == i )

    def test_odd_even_split(self):
        oes = OddEvenPartitioner()
        spl = Splitter(attr='partitions')

        splits = [ list(spl.generate(p)) for p in oes.generate(self.data) ]

        self.assertTrue(len(splits) == 2)

        for i,p in enumerate(splits):
            self.assertTrue( len(p) == 2 )
            self.assertTrue( p[0].nsamples == 50 )
            self.assertTrue( p[1].nsamples == 50 )

        assert_array_equal(splits[0][1].sa['chunks'].unique, [1, 3, 5, 7, 9])
        assert_array_equal(splits[0][0].sa['chunks'].unique, [0, 2, 4, 6, 8])
        assert_array_equal(splits[1][0].sa['chunks'].unique, [1, 3, 5, 7, 9])
        assert_array_equal(splits[1][1].sa['chunks'].unique, [0, 2, 4, 6, 8])

        # check if it works on pure odd and even chunk ids
        moresplits = [ list(spl.generate(p)) for p in oes.generate(splits[0][0])]

        for split in moresplits:
            self.assertTrue(split[0] != None)
            self.assertTrue(split[1] != None)

    def test_half_split(self):
        hs = HalfPartitioner()
        spl = Splitter(attr='partitions')

        splits = [ list(spl.generate(p)) for p in hs.generate(self.data) ]

        self.assertTrue(len(splits) == 2)

        for i,p in enumerate(splits):
            self.assertTrue( len(p) == 2 )
            self.assertTrue( p[0].nsamples == 50 )
            self.assertTrue( p[1].nsamples == 50 )

        assert_array_equal(splits[0][1].sa['chunks'].unique, [0, 1, 2, 3, 4])
        assert_array_equal(splits[0][0].sa['chunks'].unique, [5, 6, 7, 8, 9])
        assert_array_equal(splits[1][1].sa['chunks'].unique, [5, 6, 7, 8, 9])
        assert_array_equal(splits[1][0].sa['chunks'].unique, [0, 1, 2, 3, 4])

        # check if it works on pure odd and even chunk ids
        moresplits = [ list(spl.generate(p)) for p in hs.generate(splits[0][0])]

        for split in moresplits:
            self.assertTrue(split[0] != None)
            self.assertTrue(split[1] != None)

    def test_n_group_split(self):
        """Test NGroupSplitter alongside with the reversal of the
        order of spit out datasets
        """
        # Test 2 groups like HalfSplitter first
        hs = NGroupPartitioner(2)

        for isreversed, splitter in enumerate((hs, hs)):
            if isreversed:
                spl = Splitter(attr='partitions', reverse=True)
            else:
                spl = Splitter(attr='partitions')
            splits = [ list(spl.generate(p)) for p in hs.generate(self.data) ]
            self.assertTrue(len(splits) == 2)

            for i, p in enumerate(splits):
                self.assertTrue( len(p) == 2 )
                self.assertTrue( p[0].nsamples == 50 )
                self.assertTrue( p[1].nsamples == 50 )

            assert_array_equal(splits[0][1-isreversed].sa['chunks'].unique,
                               [0, 1, 2, 3, 4])
            assert_array_equal(splits[0][isreversed].sa['chunks'].unique,
                               [5, 6, 7, 8, 9])
            assert_array_equal(splits[1][1-isreversed].sa['chunks'].unique,
                               [5, 6, 7, 8, 9])
            assert_array_equal(splits[1][isreversed].sa['chunks'].unique,
                               [0, 1, 2, 3, 4])

        # check if it works on pure odd and even chunk ids
        moresplits = [ list(spl.generate(p)) for p in hs.generate(splits[0][0])]

        for split in moresplits:
            self.assertTrue(split[0] != None)
            self.assertTrue(split[1] != None)

        # now test more groups
        s5 = NGroupPartitioner(5)

        # get the splits
        for isreversed, s5splitter in enumerate((s5, s5)):
            if isreversed:
                spl = Splitter(attr='partitions', reverse=True)
            else:
                spl = Splitter(attr='partitions')
            splits = [ list(spl.generate(p)) for p in s5splitter.generate(self.data) ]

            # must have 10 splits
            self.assertTrue(len(splits) == 5)

            # check split content
            assert_array_equal(splits[0][1-isreversed].sa['chunks'].unique,
                               [0, 1])
            assert_array_equal(splits[0][isreversed].sa['chunks'].unique,
                               [2, 3, 4, 5, 6, 7, 8, 9])
            assert_array_equal(splits[1][1-isreversed].sa['chunks'].unique,
                               [2, 3])
            assert_array_equal(splits[1][isreversed].sa['chunks'].unique,
                               [0, 1, 4, 5, 6, 7, 8, 9])
            # ...
            assert_array_equal(splits[4][1-isreversed].sa['chunks'].unique,
                               [8, 9])
            assert_array_equal(splits[4][isreversed].sa['chunks'].unique,
                               [0, 1, 2, 3, 4, 5, 6, 7])

        # Test for too many groups
        def splitcall(spl, dat):
            return list(spl.generate(dat))
        s20 = NGroupPartitioner(20)
        self.assertRaises(ValueError,splitcall,s20,self.data)

    @reseed_rng()
    def test_nfold_random_counted_selection_partitioner(self):
        # Lets get somewhat extensive but complete one and see if
        # everything is legit. 0.5 must correspond to 50%, in our case
        # 5 out of 10 unique chunks
        split_partitions = [
            tuple(x.sa.partitions)
            for x in NFoldPartitioner(0.5).generate(self.data)]
        # 252 is # of combinations of 5 from 10
        assert_equal(len(split_partitions), 252)

        # verify that all of them are unique
        assert_equal(len(set(split_partitions)), 252)

        # now let's limit our query
        kwargs = dict(count=10, selection_strategy='random')
        split10_partitions = [
            tuple(x.sa.partitions)
            for x in NFoldPartitioner(5, **kwargs).generate(self.data)]
        split10_partitions_ = [
            tuple(x.sa.partitions)
            for x in NFoldPartitioner(0.5, **kwargs).generate(self.data)]
        # to make sure that I deal with sets of tuples correctly:
        assert_equal(len(set(split10_partitions)), 10)
        assert_equal(len(split10_partitions), 10)
        assert_equal(len(split10_partitions_), 10)
        # and they must differ (same ones are possible but very very unlikely)
        assert_not_equal(split10_partitions, split10_partitions_)
        # but every one of them must be within known exhaustive set
        assert_equal(set(split_partitions).intersection(split10_partitions),
                     set(split10_partitions))
        assert_equal(set(split_partitions).intersection(split10_partitions_),
                     set(split10_partitions_))

    @reseed_rng()
    def test_nfold_random_counted_selection_partitioner_huge(self):
        # Just test that it completes in a reasonable time and does
        # not blow up as if would do if it was not limited by count
        kwargs = dict(count=10)
        ds = dataset_wizard(np.arange(1000).reshape((-1, 1)),
                            targets=range(1000),
                            chunks=range(500)*2)
        split_partitions_random = [
            tuple(x.sa.partitions)
            for x in NFoldPartitioner(100,  selection_strategy='random',
                                      **kwargs).generate(ds)]
        assert_equal(len(split_partitions_random), 10) # we get just 10

    def test_custom_split(self):
        #simulate half splitter
        hs = CustomPartitioner([(None,[0,1,2,3,4]),(None,[5,6,7,8,9])])
        spl = Splitter(attr='partitions')
        splits = [ list(spl.generate(p)) for p in hs.generate(self.data) ]
        self.assertTrue(len(splits) == 2)

        for i,p in enumerate(splits):
            self.assertTrue( len(p) == 2 )
            self.assertTrue( p[0].nsamples == 50 )
            self.assertTrue( p[1].nsamples == 50 )

        assert_array_equal(splits[0][1].sa['chunks'].unique, [0, 1, 2, 3, 4])
        assert_array_equal(splits[0][0].sa['chunks'].unique, [5, 6, 7, 8, 9])
        assert_array_equal(splits[1][1].sa['chunks'].unique, [5, 6, 7, 8, 9])
        assert_array_equal(splits[1][0].sa['chunks'].unique, [0, 1, 2, 3, 4])

        # check fully customized split with working and validation set specified
        cs = CustomPartitioner([([0,3,4],[5,9])])
        # we want to discared the unselected partition of the data, hence attr_value
        # these two splitters should do exactly the same thing
        splitters = (Splitter(attr='partitions', attr_values=[1,2]),
                     Splitter(attr='partitions', ignore_values=(0,)))
        for spl in splitters:
            splits = [ list(spl.generate(p)) for p in cs.generate(self.data) ]
            self.assertTrue(len(splits) == 1)

            for i,p in enumerate(splits):
                self.assertTrue( len(p) == 2 )
                self.assertTrue( p[0].nsamples == 30 )
                self.assertTrue( p[1].nsamples == 20 )

            self.assertTrue((splits[0][1].sa['chunks'].unique == [5, 9]).all())
            self.assertTrue((splits[0][0].sa['chunks'].unique == [0, 3, 4]).all())

    def test_label_splitter(self):
        oes = OddEvenPartitioner(attr='targets')
        spl = Splitter(attr='partitions')

        splits = [ list(spl.generate(p)) for p in oes.generate(self.data) ]

        assert_array_equal(splits[0][0].sa['targets'].unique, [0,2])
        assert_array_equal(splits[0][1].sa['targets'].unique, [1,3])
        assert_array_equal(splits[1][0].sa['targets'].unique, [1,3])
        assert_array_equal(splits[1][1].sa['targets'].unique, [0,2])

    def test_counted_splitting(self):
        spl = Splitter(attr='partitions')
        # count > #chunks, should result in 10 splits
        nchunks = len(self.data.sa['chunks'].unique)
        for strategy in Partitioner._STRATEGIES:
            for count, target in [ (nchunks*2, nchunks),
                                   (nchunks, nchunks),
                                   (nchunks-1, nchunks-1),
                                   (3, 3),
                                   (0, 0),
                                   (1, 1)
                                   ]:
                nfs = NFoldPartitioner(cvtype=1, count=count,
                                       selection_strategy=strategy)
                splits = [ list(spl.generate(p)) for p in nfs.generate(self.data) ]
                self.assertTrue(len(splits) == target)
                chosenchunks = [int(s[1].uniquechunks) for s in splits]

                # Test if configuration matches as well
                nsplits_cfg = len(nfs.get_partition_specs(self.data))
                self.assertEqual(nsplits_cfg, target)

                # Check if "lastsplit" dsattr was assigned appropriately
                nsplits = len(splits)
                if nsplits > 0:
                    # dummy-proof testing of last split
                    for ds_ in splits[-1]:
                        self.assertTrue(ds_.a.lastpartitionset)
                    # test all now
                    for isplit,split in enumerate(splits):
                        for ds_ in split:
                            ds_.a.lastpartitionset == isplit==nsplits-1

                # Check results of different strategies
                if strategy == 'first':
                    self.assertEqual(chosenchunks, range(target))
                elif strategy == 'equidistant':
                    if target == 3:
                        self.assertEqual(chosenchunks, [0, 3, 7])
                elif strategy == 'random':
                    # none is selected twice
                    self.assertTrue(len(set(chosenchunks)) == len(chosenchunks))
                    self.assertTrue(target == len(chosenchunks))
                else:
                    raise RuntimeError, "Add unittest for strategy %s" \
                          % strategy

    def test_discarded_boundaries(self):
        ds = datasets['hollow']
        # four runs
        ds.sa['chunks'] = np.repeat(np.arange(4), 10)
        # do odd even splitting for lots of boundaries in few splits
        part = ChainNode([OddEvenPartitioner(),
                          StripBoundariesSamples('chunks', 1, 2)])

        parts = [d.samples.sid for d in part.generate(ds)]

        # both dataset should have the same samples, because the boundaries are
        # identical and the same sample should be stripped
        assert_array_equal(parts[0], parts[1])

        # we strip 3 samples per boundary
        assert_equal(len(parts[0]), len(ds) - (3 * 3))

        for i in [9, 10, 11, 19, 20, 21, 29, 30, 31]:
            assert_false(i in parts[0])

    def test_slicing(self):
        hs = HalfPartitioner()
        spl = Splitter(attr='partitions')
        splits = list(hs.generate(self.data))
        for s in splits:
            # partitioned dataset shared the data
            assert_true(s.samples.base is self.data.samples)
        splits = [ list(spl.generate(p)) for p in hs.generate(self.data) ]

        # with numpy 1.7.0b1 "chaining" was deprecated so let's create
        # check function appropriate for the given numpy version
        _a = np.arange(5)
        __a = _a[:4][:3]
        if __a.base is _a:
            # 1.7.0b1
            def is_the_same_base(x, base=self.data.samples):
                return x.base is base
        elif __a.base.base is _a:
            # prior 1.7.0b1
            def is_the_same_base(x, base=self.data.samples):
                return x.base.base is base
        else:
            raise RuntimeError("Uknown handling of .base by numpy")

        for s in splits:
            # we get slicing all the time
            assert_true(is_the_same_base(s[0].samples))
            assert_true(is_the_same_base(s[1].samples))
        spl = Splitter(attr='partitions', noslicing=True)
        splits = [ list(spl.generate(p)) for p in hs.generate(self.data) ]
        for s in splits:
            # we no slicing at all
            assert_false(s[0].samples.base is self.data.samples)
            assert_false(s[1].samples.base is self.data.samples)
        nfs = NFoldPartitioner()
        spl = Splitter(attr='partitions')
        splits = [ list(spl.generate(p)) for p in nfs.generate(self.data) ]
        for i, s in enumerate(splits):
            # training only first and last split
            if i == 0 or i == len(splits) - 1:
                assert_true(is_the_same_base(s[0].samples))
            else:
                assert_true(s[0].samples.base is None)
            # we get slicing all the time
            assert_true(is_the_same_base(s[1].samples))
        step_ds = Dataset(np.random.randn(20,2),
                          sa={'chunks': np.tile([0,1], 10)})
        oes = OddEvenPartitioner()
        spl = Splitter(attr='partitions')
        splits = list(oes.generate(step_ds))
        for s in splits:
            # partitioned dataset shared the data
            assert_true(s.samples.base is step_ds.samples)
        splits = [ list(spl.generate(p)) for p in oes.generate(step_ds) ]
        assert_equal(len(splits), 2)
        for s in splits:
            # we get slicing all the time
            assert_true(is_the_same_base(s[0].samples, step_ds.samples))
            assert_true(is_the_same_base(s[1].samples, step_ds.samples))

def suite():  # pragma: no cover
    return unittest.makeSuite(SplitterTests)

if __name__ == '__main__':  # pragma: no cover
    import runner
    runner.run()

max_issues_repo_name: mih/PyMVPA

max_forks_repo_name: mih/PyMVPA


text:
import logging
from unittest.mock import Mock, patch

import numpy as np
import pytest
from ophyd.sim import make_fake_device

from pcdsdevices.epics_motor import OffsetMotor
from pcdsdevices.lodcm import (CHI1, CHI2, H1N, H2N, LODCM, Y1, Y2, Dectris,
                               Diode, Foil, LODCMEnergyC, LODCMEnergySi,
                               SimFirstTower, SimLODCM, SimSecondTower, YagLom)

logger = logging.getLogger(__name__)

def motor_setup(mot, pos=0):
    mot.user_readback.sim_put(0)
    mot.user_setpoint.sim_put(0)
    mot.user_setpoint.sim_set_limits((0, 0))
    mot.motor_spg.sim_put(2)

def make_fake_offset_ims(prefix, motor_pos=0, user_offset=0):
    class MyOffsetIMS(OffsetMotor):
        def set_current_position(self, position):
            # TODO: is this supposed to be:
            # position - motor.position
            # or the way i have it here?
            logger.debug('Set current position to %d', position)
            new_offset = self.motor.position - position
            self.user_offset.sim_put(new_offset)

    fake_ims = make_fake_device(MyOffsetIMS)('PREFIX',
                                             motor_prefix='MOTOR:PREFIX',
                                             name='fake_ims')

    motor_setup(fake_ims.motor, pos=motor_pos)
    fake_ims.user_offset.sim_put(user_offset)
    return fake_ims

@pytest.fixture(scope='function')
def fake_lodcm():
    FakeLODCM = make_fake_device(SimLODCM)

    # After the fake_lodcm setup:
    # fake_lodcm will have get_material() = 'C' by default
    # and get_reflection() = (1, 1, 1) by default

    lodcm = FakeLODCM('FAKE:LOM', name='fake_lom')
    lodcm.h1n_state.state.sim_put(1)
    lodcm.h1n_state.state.sim_set_enum_strs(['Unknown'] + H1N.states_list)
    lodcm.yag.state.sim_put(1)
    lodcm.yag.state.sim_set_enum_strs(['Unknown'] + YagLom.states_list)
    lodcm.dectris.state.sim_put(1)
    lodcm.dectris.state.sim_set_enum_strs(['Unknown'] + Dectris.states_list)
    lodcm.diode.state.sim_put(1)
    lodcm.diode.state.sim_set_enum_strs(['Unknown'] + Diode.states_list)
    lodcm.foil.state.sim_put(1)
    lodcm.foil.state.sim_set_enum_strs(['Unknown'] + Foil.states_list)
    lodcm.h1n_state.state.sim_put(1)

    # additional states for Crystal Tower 1
    # set y1_state and chi1_state to 'C'
    lodcm.y1_state.state.sim_put(1)
    lodcm.y1_state.state.sim_set_enum_strs(
        ['Unknown'] + Y1.states_list)
    lodcm.chi1_state.state.sim_put(1)
    lodcm.chi1_state.state.sim_set_enum_strs(
        ['Unknown'] + CHI1.states_list)

    # additional states for Crystal Tower 2
    # set y2_state, chi2_state and h2_state to 'C'
    lodcm.y2_state.state.sim_put(1)
    lodcm.y2_state.state.sim_set_enum_strs(
        ['Unknown'] + Y2.states_list)
    lodcm.chi2_state.state.sim_put(1)
    lodcm.chi2_state.state.sim_set_enum_strs(
        ['Unknown'] + CHI2.states_list)
    lodcm.h2n_state.state.sim_put(1)
    lodcm.h2n_state.state.sim_set_enum_strs(
        ['Unknown'] + H2N.states_list)

    # set the reflection to default to (1, 1, 1) for both towers and both
    # materials
    lodcm.tower1.diamond_reflection.sim_put((1, 1, 1))
    lodcm.tower1.silicon_reflection.sim_put((1, 1, 1))
    lodcm.tower2.diamond_reflection.sim_put((1, 1, 1))
    lodcm.tower2.silicon_reflection.sim_put((1, 1, 1))

    lodcm.th1Si = make_fake_offset_ims('TH1SI:PREFIX', user_offset=-23)
    lodcm.th2Si = make_fake_offset_ims('TH2C:PREFIX', user_offset=-23)
    lodcm.th1C = make_fake_offset_ims('TH1C:PREFIX')
    lodcm.th2C = make_fake_offset_ims('TH2C:PREFIX')

    lodcm.z1Si = make_fake_offset_ims('TH1SI:PREFIX')
    lodcm.z2Si = make_fake_offset_ims('TH2C:PREFIX')
    lodcm.z1C = make_fake_offset_ims('TH1C:PREFIX')
    lodcm.z2C = make_fake_offset_ims('TH2C:PREFIX')

    lodcm.tower2.x2_retry_deadband.sim_put(-1)
    lodcm.tower2.z2_retry_deadband.sim_put(-1)

    return lodcm

# Crystal Tower 1 Setup - usefull for testing CrystalTower1 independently
@pytest.fixture(scope='function')
def fake_tower1():
    FakeLODCM = make_fake_device(SimFirstTower)
    tower1 = FakeLODCM('FAKE:TOWER1', name='fake_t1')

    tower1.diamond_reflection.sim_put((1, 1, 1))
    tower1.silicon_reflection.sim_put((1, 1, 1))
    # set h1n_state to 'OUT'
    tower1.h1n_state.state.sim_put(1)
    tower1.h1n_state.state.sim_set_enum_strs(
        ['Unknown'] + H1N.states_list)
    # set y1_state to 'C'
    tower1.y1_state.state.sim_put(1)
    tower1.y1_state.state.sim_set_enum_strs(
        ['Unknown'] + Y1.states_list)
    # set chi1_state to 'C'
    tower1.chi1_state.state.sim_put(1)
    tower1.chi1_state.state.sim_set_enum_strs(
        ['Unknown'] + CHI1.states_list)
    return tower1

# Crystal Tower 2 Setup - usefull for testing CrystalTower1 independently
@pytest.fixture(scope='function')
def fake_tower2():
    FakeLODCM = make_fake_device(SimSecondTower)
    tower2 = FakeLODCM('FAKE:TOWER2', name='fake_t2')

    tower2.diamond_reflection.sim_put((2, 2, 2))
    tower2.silicon_reflection.sim_put((2, 2, 2))
    # set y2_state to 'C'
    tower2.y2_state.state.sim_put(1)
    tower2.y2_state.state.sim_set_enum_strs(
        ['Unknown'] + Y2.states_list)
    # set chi2_state to 'C'
    tower2.chi2_state.state.sim_put(1)
    tower2.chi2_state.state.sim_set_enum_strs(
        ['Unknown'] + CHI2.states_list)
    # set h2n_state to 'C'
    tower2.h2n_state.state.sim_put(1)
    tower2.h2n_state.state.sim_set_enum_strs(
        ['Unknown'] + H2N.states_list)

    return tower2

# LODCM Energy Si - usefull for testing LODCMEnergySi independently
@pytest.fixture(scope='function')
def fake_energy_si():
    FakeLODCMEnergy = make_fake_device(LODCMEnergySi)
    energy = FakeLODCMEnergy('FAKE:ENERGY:SI', name='fake_energy_si')

    motor_setup(energy.th1Si.motor)
    motor_setup(energy.th2Si.motor)
    motor_setup(energy.z1Si.motor)
    motor_setup(energy.z2Si.motor)
    motor_setup(energy.dr)

    return energy

# LODCM Energy C - usefull for testing LODCMEnergyC independently
@pytest.fixture(scope='function')
def fake_energy_c():
    FakeLODCMEnergy = make_fake_device(LODCMEnergyC)
    energy = FakeLODCMEnergy('FAKE:ENERGY:C', name='fake_energy_c')

    motor_setup(energy.th1C.motor)
    motor_setup(energy.th2C.motor)
    motor_setup(energy.z1C.motor)
    motor_setup(energy.z2C.motor)
    motor_setup(energy.dr)

    return energy

def test_lodcm_destination(fake_lodcm):
    logger.debug('test_lodcm_destination')
    lodcm = fake_lodcm
    dest = lodcm.destination
    assert isinstance(dest, list)
    for d in dest:
        assert isinstance(d, str)

    lodcm.h1n_state.move('OUT')
    assert len(lodcm.destination) == 1
    assert lodcm.main_line in lodcm.destination

    lodcm.h1n_state.move('C')
    assert lodcm.main_line in lodcm.destination
    assert lodcm.mono_line in lodcm.destination
    assert len(lodcm.destination) == 2

    # Block the mono line
    lodcm.yag.move('IN')
    assert len(lodcm.destination) == 1
    assert lodcm.mono_line not in lodcm.destination
    assert lodcm.main_line in lodcm.destination
    lodcm.h1n_state.move('Si')
    assert len(lodcm.destination) == 0
    lodcm.yag.move('OUT')
    assert len(lodcm.destination) == 1
    assert lodcm.mono_line in lodcm.destination

    # Unknown state
    lodcm.h1n_state.state.sim_put('Unknown')
    assert len(lodcm.destination) == 0

def test_lodcm_branches(fake_lodcm):
    logger.debug('test_lodcm_branches')
    lodcm = fake_lodcm
    branches = lodcm.branches
    assert isinstance(branches, list)
    for b in branches:
        assert isinstance(b, str)

def test_lodcm_remove_dia(fake_lodcm):
    logger.debug('test_lodcm_remove_dia')
    lodcm = fake_lodcm
    lodcm.yag.insert(wait=True)
    cb = Mock()
    lodcm.remove_dia(moved_cb=cb, wait=True)
    assert cb.called
    assert lodcm.yag.position == 'OUT'

def test_hutch_foils():
    FakeFoil = make_fake_device(Foil)
    assert 'Zn' in FakeFoil('XPP', name='foil').in_states
    assert 'Ge' in FakeFoil('XCS', name='foil').in_states

def test_move_energy(fake_lodcm):
    lom = fake_lodcm
    # with material 'Si' and reflection as (1, 1, 1)
    with patch('pcdsdevices.lodcm.LODCM.get_material',
               return_value='Si'):
        lom.energy.move(10, wait=False)
        assert np.isclose(lom.energy.th1Si.wm(), 11.402710639982848)
        assert np.isclose(lom.energy.th2Si.wm(), 11.402710639982848)
        assert np.isclose(lom.energy.z1Si.wm(), -713.4828146545175)
        assert np.isclose(lom.energy.z2Si.wm(), 713.4828146545175)
        assert np.isclose(lom.energy.dr.wm(), 22.805421279965696)

def test_tweak_x(fake_lodcm):
    lom = fake_lodcm
    # with th2Si == 23 (init_pos) in fake_lodcm setup
    # with x = 3 => z = 2.897066324421222
    lom.tweak_x(3, 'Si', wait=False)
    assert lom.x2.wm() == 3
    assert lom.z2.wm() == 2.897066324421222
    # should move relative to current position
    # x2 = 3, z2 = 2.897066324421222
    lom.tweak_x(10, 'Si', wait=False)
    # x = 10 + current position
    # z = 9.656887748070739 + current position
    assert lom.x2.wm() == 13
    assert lom.z2.wm() == 9.656887748070739 + 2.897066324421222

def test_tweak_parallel(fake_lodcm):
    lom = fake_lodcm
    # with th2Si == 23 (init_pos) in fake_lodcm setup
    # with Si, x2 = 0, z2 = 0 (init_pos)
    # x2.mvr(p * np.sin(th * np.pi / 180)) => 0 + 1.1721933854678213
    # z2.mvr(p * np.cos(th * np.pi / 180)) => 24 + 2.761514560357321
    with patch('pcdsdevices.lodcm.LODCM.get_material',
               return_value='Si'):
        lom.tweak_parallel(3)
        assert lom.x2.wm() == 1.1721933854678213
        assert lom.z2.wm() == 2.761514560357321
        lom.tweak_parallel(3)
        # with now x2 = 1.1721933854678213, z2 = 2.761514560357321
        assert lom.x2.wm() == 1.1721933854678213 * 2
        assert lom.z2.wm() == 2.761514560357321 * 2

def test_set_energy(fake_lodcm, monkeypatch):
    lom = fake_lodcm
    lom.set_energy(10, material='Si', reflection=(1, 1, 1))

    assert lom.th1Si.wm() == 11.402710639982848
    assert lom.th2Si.wm() == 11.402710639982848
    assert lom.z1Si.wm() == -713.4828146545175
    assert lom.z2Si.wm() == 713.4828146545175

def test_tower1_crystal_type(fake_tower1):
    tower1 = fake_tower1
    # tower1 si: y1_state = Si, chi1_state = Si, h1n = Si or Out
    tower1.h1n_state.move('OUT')
    tower1.y1_state.move('Si')
    tower1.chi1_state.move('Si')
    assert tower1.is_silicon()
    tower1.h1n_state.move('Si')
    assert tower1.is_silicon()
    # tower1 c: y1_state = C, chi1_state = C, h1n = C or Out
    tower1.h1n_state.move('OUT')
    tower1.y1_state.move('C')
    tower1.chi1_state.move('C')
    assert tower1.is_diamond()
    tower1.h1n_state.move('C')
    assert tower1.is_diamond()
    # tower1 non matching states:
    tower1.h1n_state.move('OUT')
    tower1.y1_state.move('C')
    tower1.chi1_state.move('Si')
    with pytest.raises(ValueError):
        tower1.get_material()

def test_tower2_crystal_type(fake_tower2):
    tower2 = fake_tower2
    # tower2 si: y2_state = Si, chi2_state = Si, h2n_state = Si
    tower2.h2n_state.move('Si')
    tower2.y2_state.move('Si')
    tower2.chi2_state.move('Si')
    assert tower2.is_silicon()
    # tower2 c: y2_state = C, chi2_state = C, h2n_state = C
    tower2.h2n_state.move('C')
    tower2.y2_state.move('C')
    tower2.chi2_state.move('C')
    assert tower2.is_diamond()
    # tower2 non matching states:
    tower2.h2n_state.move('C')
    tower2.y2_state.move('Si')
    tower2.chi2_state.move('Si')
    with pytest.raises(ValueError):
        tower2.get_material()

def test_get_reflection_tower1(fake_tower1):
    tower1 = fake_tower1
    # defalts to (1, 1, 1), see fake_tower1 setup
    # we also have C as the material
    res = tower1.get_reflection()
    assert res == (1, 1, 1)
    tower1.diamond_reflection.sim_put((2, 2, 2))
    res = tower1.get_reflection()
    assert res == (2, 2, 2)
    with patch('pcdsdevices.lodcm.CrystalTower1.is_diamond',
               return_value=None):
        with patch('pcdsdevices.lodcm.CrystalTower1.is_silicon',
                   return_value=None):
            with pytest.raises(ValueError):
                tower1.get_reflection()

def test_get_reflection_tower2(fake_tower2):
    tower2 = fake_tower2
    # defalts to (2, 2, 2), see fake_tower2 setup
    # we also have C as the material
    res = tower2.get_reflection()
    assert res == (2, 2, 2)
    tower2.diamond_reflection.sim_put((2, 0, 0))
    res = tower2.get_reflection()
    assert res == (2, 0, 0)
    with patch('pcdsdevices.lodcm.CrystalTower2.is_diamond',
               return_value=None):
        with patch('pcdsdevices.lodcm.CrystalTower2.is_silicon',
                   return_value=None):
            with pytest.raises(ValueError):
                tower2.get_reflection()

def test_get_reflection_lodcm(fake_lodcm):
    lodcm = fake_lodcm
    # as of now tower 1: (1, 1, 1), tower 2: (1, 1, 1) - match
    res = lodcm.get_reflection()
    assert res == (1, 1, 1)
    # change tower 2: (2, 2, 2) - should not match
    lodcm.tower2.diamond_reflection.sim_put((2, 2, 2))
    lodcm.tower2.silicon_reflection.sim_put((2, 2, 2))
    with pytest.raises(ValueError):
        lodcm.get_reflection()

def test_get_material_lodcm(fake_lodcm):
    lom = fake_lodcm
    # defaults to C, in fake_lodcm stup
    assert lom.get_material() == 'C'

    # change the material for one of the towers so that it does not match
    # CrystalTower2 = 'Si', CrystalTower1 = 'C' - do not match
    lom.tower2.h2n_state.move('Si')
    lom.tower2.y2_state.move('Si')
    lom.tower2.chi2_state.move('Si')

    with pytest.raises(ValueError):
        lom.get_material()

def test_calc_geometry_c(fake_energy_c):
    energy = fake_energy_c
    logger.info('Testing with C, (1, 1, 1)')
    with patch('pcdsdevices.lodcm.LODCMEnergyC.get_reflection',
               return_value=(1, 1, 1)):
        th, z = energy.calc_geometry(energy=10)
        assert th == 17.51878596767417
        assert z == 427.8469911590626
        th, z = energy.calc_geometry(energy=6)
        assert th == 30.112367750143974
        assert z == 171.63966796710216

    logger.info('Testing with C, (2, 2, 0)')
    with patch('pcdsdevices.lodcm.LODCMEnergyC.get_reflection',
               return_value=(2, 2, 0)):
        th, z = energy.calc_geometry(energy=10)
        assert th == 29.443241721774093
        assert z == 181.06811018121837
        th, z = energy.calc_geometry(energy=6)
        assert th == 55.01163934185574
        assert z == -109.32912449140694

def test_calc_geometry_si(fake_energy_si):
    energy = fake_energy_si
    logger.info('Testing with Si, (1, 1, 1)')
    with patch("pcdsdevices.lodcm.LODCMEnergySi.get_reflection",
               return_value=(1, 1, 1)):
        th, z = energy.calc_geometry(energy=10)
        assert th == 11.402710639982848
        assert z == 713.4828146545175
        th, z = energy.calc_geometry(energy=6)
        assert th == 19.23880622548293
        assert z == 377.45432488131866

    logger.info('Testing with Si, (2, 2, 0)')
    with patch("pcdsdevices.lodcm.LODCMEnergySi.get_reflection",
               return_value=(2, 2, 0)):
        th, z = energy.calc_geometry(energy=10)
        assert th == 18.835297041786244
        assert z == 388.56663653387835
        th, z = energy.calc_geometry(energy=6)
        assert th == 32.55312408478254
        assert z == 139.21560118646275

def test_get_energy_c(fake_energy_c):
    energy = fake_energy_c
    energy.th1C.user_offset.sim_put(-23)

    with patch("pcdsdevices.lodcm.LODCMEnergyC.get_reflection",
               return_value=(1, 1, 1)):
        # offset of 23, motor position 0
        # current th1C.wm() should be 23
        assert energy.th1C.wm() == 23
        res = energy.get_energy()
        assert res == 7.7039801344046515

def test_get_energy_si(fake_energy_si):
    energy = fake_energy_si

    energy.th1Si.user_offset.sim_put(-23)
    with patch("pcdsdevices.lodcm.LODCMEnergySi.get_reflection",
               return_value=(1, 1, 1)):
        # offset of 23, motor position 0
        # current th1Si.wm() should be 23
        assert energy.th1Si.wm() == 23
        res = energy.get_energy()
        assert res == 5.059840436879476

def test_forward_si(fake_energy_si):
    energy = fake_energy_si

    energy.th1Si.user_offset.sim_put(-23)
    with patch("pcdsdevices.lodcm.LODCMEnergySi.get_reflection",
               return_value=(1, 1, 1)):
        res = energy.get_energy()
        assert res == 5.059840436879476
        res = energy.forward(5.059840436879476)
        assert np.isclose(res.dr, 46)
        assert np.isclose(res.th1Si, 23)
        assert np.isclose(res.th2Si, 23)
        assert np.isclose(res.z1Si, -289.70663244212216)
        assert np.isclose(res.z2Si, 289.70663244212216)

def test_forward_c(fake_energy_c):
    energy = fake_energy_c

    energy.th1C.user_offset.sim_put(-23)
    with patch("pcdsdevices.lodcm.LODCMEnergyC.get_reflection",
               return_value=(1, 1, 1)):
        res = energy.get_energy()
        assert res == 7.7039801344046515
        res = energy.forward(7.7039801344046515)
        assert np.isclose(res.dr, 46)
        assert np.isclose(res.th1C, 23)
        assert np.isclose(res.th2C, 23)
        assert np.isclose(res.z1C, -289.70663244212216)
        assert np.isclose(res.z2C, 289.70663244212216)

def test_inverse_si(fake_energy_si):
    energy = fake_energy_si
    energy.th1Si.user_offset.sim_put(-23)
    with patch("pcdsdevices.lodcm.LODCMEnergySi.get_reflection",
               return_value=(1, 1, 1)):
        res = energy.inverse(23)
        assert np.isclose(res[0], 5.05984044)

def test_inverse_c(fake_energy_c):
    energy = fake_energy_c
    energy.th1C.user_offset.sim_put(-23)
    with patch("pcdsdevices.lodcm.LODCMEnergyC.get_reflection",
               return_value=(1, 1, 1)):
        res = energy.inverse(23)
        assert np.isclose(res[0], 7.7039801344046515)

def test_lodcm_move_energy_si(fake_lodcm):
    lodcm = fake_lodcm
    # with material 'Si', defaulted in fake_lodcm setup
    with patch("pcdsdevices.lodcm.LODCMEnergySi.get_reflection",
               return_value=(1, 1, 1)):
        with patch('pcdsdevices.lodcm.LODCM.get_material',
                   return_value='Si'):
            lodcm.energy_si.move(10, wait=False)
            assert np.isclose(lodcm.energy.th1Si.wm(), 11.402710639982848)
            assert np.isclose(lodcm.energy.th2Si.wm(), 11.402710639982848)
            assert np.isclose(lodcm.energy.dr.wm(), 11.402710639982848*2)
            assert np.isclose(lodcm.energy.z1Si.wm(), -713.4828146545175)
            assert np.isclose(lodcm.energy.z2Si.wm(), 713.4828146545175)

@pytest.mark.timeout(5)
def test_lodcm_disconnected():
    LODCM('TST:LOM', name='test_lom')

max_stars_repo_name: ghalym/pcdsdevices

max_issues_repo_name: ghalym/pcdsdevices

max_forks_repo_name: ghalym/pcdsdevices


text:
import warnings
import nibabel as nib
import pathlib as pal
from scipy import io as sio

def niak_scrubbing(img_p, extra_p, out_p, clobber=False):
    """

    :param img_p: pathlib path to the functional image file
    :param extra_p: pathlib path to the .mat file that contains the scrubbing mask
    :param out_p: pathlib path to the output functional image file
    :param clobber: if true, overwrite existing output image file
    :return:
    """
    if not issubclass(type(out_p), pal.Path):
        out_p = pal.Path(out_p)

    if out_p.is_file() and not clobber:
        warnings.warn(f'{out_p.name} already exists and clobber = {clobber}. Not touching anything:\n {out_p}')
        return out_p.is_file()

    img = nib.load(str(img_p))
    extra = sio.loadmat(str(extra_p))

    scrub_mask = extra['mask_scrubbing'].squeeze()
    # Check that the length of the mask matches that of the temporal dimension
    if not len(scrub_mask) == img.shape[-1]:
        raise Exception(f'Shape mismatch between {img_p.name} and {extra_p.name}: {img.shape} vs {len(scrub_mask)}')

    masked_img = nib.Nifti1Image(img.get_data()[..., scrub_mask != 1], affine=img.affine, header=img.header)
    nib.save(masked_img, str(out_p))
    return out_p.is_file()

reason: from scipy


text:
'''
    Classes and functions for plotting the actual results of the simulations.
    Interpolating version
'''

# Dependencies: gitpython,cycler

import collections
import itertools
import json
import os
import os.path as op
import sys

#import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from cycler import cycler

from main_help import *

plt.close('all')

myColors = [[27, 158, 119],
 [217, 95, 2],
 [117, 112, 179],
 [231, 41, 138],
 [102, 166, 30],
 [230, 171, 2],
 [166, 118, 29],
 [102, 102, 102]]

sty = op.join(thispath, "swept.mplstyle")
plt.style.use(sty)

hxColors = ["#{:02x}{:02x}{:02x}".format(r,g,b) for r, g, b in myColors]

plt.rc('axes', prop_cycle=cycler('color', hxColors) + cycler('marker', ['D', 'o', 'h', '*', '^', 'x', 'v', '8']))

xlbl = "Grid Size"

dftype=pd.core.frame.DataFrame #Compare types

schemes = ["Classic", "Swept"]

schs = dict(zip([k[0] for k in schemes], schemes))

meas = {"time": "time per timestep (us)", "efficiency": "MGridPts/s", "spd": "Speedup", "spdg":"Speedup vs GPU", "btpb":"Best tpb comparison"}

fc = {'time':'min', 'efficiency': 'max'}

def rowSelect(df, n):
    return df.iloc[::len(df.index)//n, :]

def formatSubplot(f):
    nsub = len(f.axes)
    if nsub == 4:
        f.tight_layout(pad=0.2, w_pad=0.75, h_pad=0.75)
        f.subplots_adjust(top=0.9, bottom=0.08, right=0.85, hspace=0.3, wspace=0.3)
    if nsub > 4:
        f.tight_layout(pad=0.2, w_pad=0.75, h_pad=0.75)
        f.subplots_adjust(top=0.9, bottom=0.08, hspace=0.3, wspace=0.3)

    return f

def cartProd(y):
    return pd.DataFrame(list(itertools.product(*y))).values

def nxRange(mmax, n):
    rng=n//10
    return np.around(np.logspace(*(np.log2(mmax['nX'])), base=2, num=n), decimals=1)[rng:-rng]

def simpleNGrid(uniq, mmax):
    xint = [uniq['tpb'], uniq['gpuA']]
    xint.append(nxRange(mmax, 100))
    return xint

# At this point the function is independent
def getBestAll(df, respvar):
    f = fc[respvar]
    ifx = 'idx'+f
    ti = list(df.columns.names)
    ti.remove('metric')
    ixvar = list(ti)[0]
    to = list(df.columns.get_level_values('metric').unique())
    to.remove(respvar)
    oxvar = to[0]

    dfb = pd.DataFrame(df[respvar].apply(f, axis=1), columns=[respvar])
    dfb[ixvar] = pd.DataFrame(df[respvar].apply(ifx, axis=1))

    bCollect = []
    for i, k in zip(dfb.index, dfb[ixvar].values):
        bCollect.append(df[oxvar].loc[i, k])

    dfb[oxvar] = bCollect
    return dfb

class Perform(object):
    def __init__(self, df, name):
        self.oFrame = df
        self.title = name
        self.cols = list(df.columns.values)
        self.nobs = len(df)
        self.xo = self.cols[:3]
        self.uniques, self.minmaxes = {}, {}
        self.iFrame = pd.DataFrame()
        self.bFrame = pd.DataFrame()

        for k in self.xo:
            tmp = self.oFrame[k].unique()
            self.uniques[k] = tmp
            self.minmaxes[k] = [tmp.min() , tmp.max()]

        self.minmaxes['nX'] = [self.oFrame.groupby(self.xo[:2]).min()['nX'].max(),
                                self.oFrame.groupby(self.xo[:2]).max()['nX'].min()]

    def __str__(self):
        ms = "%s \n %s \n Unique Exog: \n" % (self.title, self.oFrame.head())

        for i, s in enumerate(self.uniques):
            ms = ms + self.cols[i] + ": \n " + str(s) + "\n"
        return ms

    def efficient(self, df=pd.DataFrame()):
        if not df.empty:
            df['efficiency'] = df['nX']/df['time']
            return df

        self.oFrame['efficiency'] = self.oFrame['nX']/self.oFrame['time']
        if not self.iFrame.empty:
            self.iFrame['efficiency'] = self.iFrame['nX']/self.iFrame['time']

    def plotRaw(self, subax, respvar, legstep=1):
        legax = self.xo[1] if subax==self.xo[0] else self.xo[0]

        drops = self.uniques[legax][1::legstep]
        saxVal = self.uniques[subax]
        ff = []
        ad = {}
        if respvar=='time':
            kwarg = {'loglog': True}
        else:
            kwarg = {'logx': True}

        for i in range(len(saxVal)//4):
            f, ai = plt.subplots(2,2)
            ap = ai.ravel()
            ff.append(f)
            for aa, t in zip(ap, saxVal[i::2]):
                ad[t] = aa

        for k, g in self.oFrame.groupby(subax):
            for kk, gg in g.groupby(legax):
                if kk in drops:
                    gg.plot(x='nX', y=respvar, ax=ad[k], grid=True, label=kk, **kwarg)

            ad[k].set_title(k)
            ad[k].set_ylabel(meas[respvar])
            ad[k].set_xlabel(xlbl)
            hd, lb = ad[k].get_legend_handles_labels()
            ad[k].legend().remove()

        for i, fi in enumerate(ff):
            fi = formatSubplot(fi)
            fi.legend(hd, lb, 'upper right', title=legax, fontsize="medium")
            fi.suptitle(self.title + " - Raw " +  respvar +  " by " + subax)

        return ff

    def getBest(self, subax, respvar, df=pd.DataFrame()):
        f = fc[respvar]
        ifx = 'idx'+f
        legax = self.xo[1] if subax==self.xo[0] else self.xo[0]
        blankidx = pd.MultiIndex(levels=[[],[]], labels=[[],[]], names=['metric', subax])
        fCollect = pd.DataFrame(columns=blankidx)

        if df.empty:
            df = self.iFrame

        for k, g in df.groupby(subax):
            gg = g.pivot(index='nX', columns=legax, values=respvar)
            fCollect[respvar, k] = gg.apply(f, axis=1)
            fCollect[legax, k] = gg.apply(ifx, axis=1)

        return fCollect

def addIntercept(perf):
    iprod = cartProd(perf.uniques['tpb'], perf.uniques['gpuA'], [0.0], [0.0])
    zros = pd.DataFrame(iprod, columns=perf.cols)
    return pd.concat(perf.oFrame, zros)


text:
import numpy as np
import sem
import matplotlib.pyplot as plt
#import cProfile

yv = np.linspace(0.,1.,7)
zv = np.linspace(0.,1.,7)
Tu = 0.1 * np.ones_like(yv)
uv = np.zeros_like(yv)
L = np.ones_like(yv)*0.1

theSEM = sem.AnisoSEM(yv, zv, Tu, Tu, Tu, uv, L, 1.)
# theSEM.plot_input()

zg, yg = np.meshgrid(zv, yv)

theSEM.set_grid(yg, zg)

nstep_cycle = 20
ncycle = 1000
dt = L[0]/nstep_cycle
tau = L[0]*ncycle
nt = int(tau/dt)

theSEM.loop(dt, nt)
#cProfile.run("theSEM.loop(dt, nt)",sort="cumtime")

theSEM.plot_output(yv, nstep_cycle)

plt.show()


text:
##~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~##
##                                                                                ##
##  This file forms part of the Underworld geophysics modelling application.         ##
##                                                                                ##
##  For full license and copyright information, please refer to the LICENSE.md file  ##
##  located at the project root, or contact the authors.                             ##
##                                                                                ##
##~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~##
import underworld as uw
import underworld._stgermain as _stgermain
import underworld.systems.sle as sle
import underworld.libUnderworld as libUnderworld

class SolveLinearSystem(_stgermain.StgCompoundComponent):
    """
    To solve for x in the equation Ax=b
    """
    _objectsDict = {  "_system" : "SystemLinearEquations" }
    _selfObjectName = "_system"

    def __init__(self, AMat, bVec, xVec, **kwargs):

        if not xVec:
            raise ValueError("You must specify a 'xVec' parameter.")
        if not AMat:
            raise ValueError("You must specify a 'AMat' parameter.")
        if not bVec:
            raise ValueError("You must specify a 'bVec' parameter.")

        if not isinstance( xVec, uw.systems.sle.SolutionVector):
            raise TypeError( "Provided 'xVec' must be of 'SolutionVector' class." )
        self._solutionVector = xVec
        if not isinstance( bVec, uw.systems.sle.SolutionVector):
            raise TypeError( "Provided 'bVec' must be of 'SolutionVector' class." )
        self._fvector = bVec
        if not isinstance( AMat, uw.systems.sle.AssembledMatrix):
            raise TypeError( "Provided 'AMat' must be of 'AssembledMatrix' class." )
        self._kmatrix = AMat

        self._solver = None
        self._swarm  = None

        super(SolveLinearSystem, self).__init__(**kwargs)

    def _setup(self):
        uw.libUnderworld.StGermain.Stg_ObjectList_Append( self._cself.forceVectors, self._fvector._cself )
        uw.libUnderworld.StGermain.Stg_ObjectList_Append( self._cself.stiffnessMatrices, self._kmatrix._cself )
        uw.libUnderworld.StGermain.Stg_ObjectList_Append( self._cself.solutionVectors, self._solutionVector._cself )

    def _add_to_stg_dict(self,componentDictionary):
        # call parents method
        super(SolveLinearSystem,self)._add_to_stg_dict(componentDictionary)

    def solve(self):
        """
        Solve system
        """
        if not self._solver:
            self._solver = uw.systems.Solver(self)
        self._solver.solve()

class MeshVariable_Projection(_stgermain.StgCompoundComponent):
    """
    This class provides functionality for projecting data
    from any underworld function onto a provided mesh variable.

    For the variable :math:`\\bf{U} =  \\bf{u}_a N_a` and function :math:`F`,
    we wish to determine appropriate values for :math:`\\bf{u}_a` such
    that :math:`\\bf{U} \\simeq F`.

    Two projection methods are supported; weighted averages and weighted
    residuals. Generally speaking, weighted averages provide robust low
    order results, while weighted residuals give higher accuracy but
    spurious results for *difficult* functions :math:`F`.

    The weighted average method is defined as:

    .. math::
         \\bf{u}_a = \\frac{\\int_{\\Omega} \\bf{F} N_a \\partial\\Omega }{\\int_{\\Omega} N_a \\partial\\Omega }

    The weighted residual method constructs an SLE which is then solved to
    determine the unknowns:

    .. math::
         \\bf{u}_a\\int_{\\Omega} N_a N_b \\partial\\Omega = \\int_{\\Omega} \\bf{F} N_b \\partial\\Omega

    Parameters
    ----------
    meshVariable : underworld.mesh.MeshVariable
        The variable you wish to project the function onto.
    fn : underworld.function.Function
        The function you wish to project.
    voronoi_swarm : underworld.swarm.Swarm
        Optional. If a voronoi_swarm is provided, voronoi type integration is
        utilised to integrate across elements. The provided swarm is used as the
        basis for the voronoi integration. If no swarm is provided, Gauss
        integration is used.
    gauss_swarm : underworld.swarm.GaussIntegrationSwarm
        If provided this gauss_swarm will be used for (gaussian) numerical
        integration rather than a default gauss integration swarm that is
        automatically generated and dependent on the element order of the mesh.
        NB: if a voronoi_swarm is defined it OVERRIDES this gauss_swarm as the
        preferred integration swarm (quadrature method).
    type : int, default=0
        Projection type.  0:`weighted average`, 1:`weighted residual`

    Notes
    -----
    Constructor must be called collectively by all processes.

    Examples
    --------
    >>> import underworld as uw
    >>> import numpy as np
    >>> mesh = uw.mesh.FeMesh_Cartesian()
    >>> U = uw.mesh.MeshVariable( mesh, 1 )

    Lets cast a constant value onto this mesh variable

    >>> const = 1.23456
    >>> projector = uw.utils.MeshVariable_Projection( U, const, type=0 )
    >>> np.allclose(U.data, const)
    False
    >>> projector.solve()
    >>> np.allclose(U.data, const)
    True

    Now cast mesh coordinates onto a vector variable

    >>> U_coord = uw.mesh.MeshVariable( mesh, 2 )
    >>> projector = uw.utils.MeshVariable_Projection( U_coord, uw.function.coord(), type=1 )
    >>> projector.solve()
    >>> np.allclose(U_coord.data, mesh.data)
    True

    Project one mesh variable onto another

    >>> U_copy = uw.mesh.MeshVariable( mesh, 2 )
    >>> projector = uw.utils.MeshVariable_Projection( U_copy, U_coord, type=1 )
    >>> projector.solve()
    >>> np.allclose(U_copy.data, U_coord.data)
    True

    Project the coords to the submesh (usually the constant mesh)

    >>> U_submesh = uw.mesh.MeshVariable( mesh.subMesh, 2 )
    >>> projector = uw.utils.MeshVariable_Projection( U_submesh, U_coord, type=1 )
    >>> projector.solve()
    >>> np.allclose(U_submesh.data, mesh.subMesh.data)
    True

    Create swarm, then project particle owning elements onto mesh

    >>> U_submesh = uw.mesh.MeshVariable( mesh.subMesh, 1 )
    >>> swarm = uw.swarm.Swarm(mesh)
    >>> swarm.populate_using_layout(uw.swarm.layouts.PerCellSpaceFillerLayout(swarm,4))
    >>> projector = uw.utils.MeshVariable_Projection( U_submesh, swarm.owningCell, type=1 )
    >>> projector.solve()
    >>> np.allclose(U_submesh.data[:,0], mesh.data_elgId)
    True

    """
    _objectsDict = {  "_system" : "SystemLinearEquations" }
    _selfObjectName = "_system"

    def __init__(self, meshVariable=None, fn=None, voronoi_swarm=None, type=0, gauss_swarm=None, **kwargs):

        if not meshVariable:
            raise ValueError("You must specify a mesh variable via the 'meshVariable' parameter.")
        if not isinstance( meshVariable, uw.mesh.MeshVariable):
            raise TypeError( "Provided 'meshVariable' must be of 'MeshVariable' class." )
        self._meshVariable = meshVariable

        if not fn:
            raise ValueError("You must specify a function via the 'fn' parameter.")
        try:
            _fn = uw.function.Function.convert(fn)
        except Exception as e:
            raise uw._prepend_message_to_exception(e, "Exception encountered. Note that provided 'fn' must be of or convertible to 'Function' class.\nEncountered exception message:\n")

        if voronoi_swarm and not isinstance(voronoi_swarm, uw.swarm.Swarm):
            raise TypeError( "Provided 'swarm' must be of 'Swarm' class." )
        self._swarm = voronoi_swarm
        if voronoi_swarm and meshVariable.mesh.elementType=='Q2':
            import warnings
            warnings.warn("Voronoi integration may yield unsatisfactory results for Q2 mesh.")

        if not type in [0,1]:
            raise ValueError( "Provided 'type' must take a value of 0 (weighted average) or 1 (weighted residual)." )
        self.type = type

        eqNum = sle.EqNumber(meshVariable)
        # create force vectors
        self._fvector = sle.AssembledVector(meshVariable, eqNum)

        # determine the required geometry mesh.  for submesh, use the parent mesh.
        geometryMesh = self._meshVariable.mesh
        if hasattr(self._meshVariable.mesh.generator, 'geometryMesh'):
            geometryMesh = self._meshVariable.mesh.generator.geometryMesh

        # setup the gauss integration swarm
        if gauss_swarm != None:
            if not isinstance(gauss_swarm,uw.swarm.GaussIntegrationSwarm):
                raise RuntimeError( "Provided 'gauss_swarm' must be a GaussIntegrationSwarm object" )
            intswarm = gauss_swarm
        else:
            intswarm = uw.swarm.GaussIntegrationSwarm(geometryMesh)

        # we will use voronoi if that has been requested by the user, else use
        # gauss integration.
        if self._swarm:
            intswarm = self._swarm._voronoi_swarm
            # need to ensure voronoi is populated now, as assembly terms will call
            # initial test functions which may require a valid voronoi swarm
            self._swarm._voronoi_swarm.repopulate()

        self._fn = _fn

        self._forceVecTerm = sle.VectorAssemblyTerm_NA__Fn(   integrationSwarm=intswarm,
                                                              assembledObject=self._fvector,
                                                              fn=_fn,
                                                              mesh=geometryMesh )

        if self.type == 0:
            # create copy guy
            self._copyMeshVariable = meshVariable.copy()
            # create unity array of required dimensionality
            self._unityArray = []
            for ii in range(self._meshVariable.nodeDofCount):
                self._unityArray.append(1.)
            self.solve = self._solve_average
        else:
            # create solutions vector
            self._solutionVector = sle.SolutionVector(meshVariable, eqNum)
            # and matrices
            self._kmatrix = sle.AssembledMatrix( self._solutionVector, self._solutionVector, rhs=self._fvector )
            # matrix term
            self._kMatTerm = sle.MatrixAssemblyTerm_NA__NB__Fn(  integrationSwarm=intswarm,
                                                                 assembledObject=self._kmatrix,
                                                                 fn = 1.0,
                                                                 mesh=geometryMesh )
            self._solver = None
            self.solve = self._solve_residual

        super(MeshVariable_Projection, self).__init__(**kwargs)

    @property
    def fn(self):
        return self._fn

    @fn.setter
    def fn(self, value):
        _fn = uw.function._function.Function.convert(value)
        if not isinstance( _fn, uw.function.Function):
            raise ValueError( "Provided 'fn' must be of, or convertible to, 'Function' class." )
        self._forceVecTerm.fn = _fn
        self._fn = _fn

    def _setup(self):
        uw.libUnderworld.StGermain.Stg_ObjectList_Append( self._cself.forceVectors, self._fvector._cself )
        if self.type == 1:
            uw.libUnderworld.StGermain.Stg_ObjectList_Append( self._cself.stiffnessMatrices, self._kmatrix._cself )
            uw.libUnderworld.StGermain.Stg_ObjectList_Append( self._cself.solutionVectors, self._solutionVector._cself )

    def _add_to_stg_dict(self,componentDictionary):
        # call parents method
        super(MeshVariable_Projection,self)._add_to_stg_dict(componentDictionary)

    def _solve_average(self):
        """
        Solve the projection for the current state of the provided function.
        """
        # first assemble \int{Fn.N}
        if self._swarm:
            self._swarm._voronoi_swarm.repopulate()
        libUnderworld.StgFEM.ForceVector_Zero( self._fvector._cself )
        libUnderworld.StgFEM.ForceVector_GlobalAssembly_General( self._fvector._cself )
        libUnderworld.StgFEM.SolutionVector_UpdateSolutionOntoNodes( self._fvector._cself );

        # now do again for \int{N}, but first create copy
        self._copyMeshVariable.data[:] = self._meshVariable.data[:]
        self._forceVecTerm.fn = self._unityArray
        libUnderworld.StgFEM.ForceVector_Zero( self._fvector._cself )
        libUnderworld.StgFEM.ForceVector_GlobalAssembly_General( self._fvector._cself )
        libUnderworld.StgFEM.SolutionVector_UpdateSolutionOntoNodes( self._fvector._cself );

        # right, now divide
        self._meshVariable.data[:] = self._copyMeshVariable.data[:] / self._meshVariable.data[:]
        # done! return to correct function
        self._forceVecTerm.fn = self._fn

    def _solve_residual(self):
        """
        Solve the projection given the current state of the provided function.
        """
        if not self._solver:
            self._solver = uw.systems.Solver(self)
        self._solver.solve()