text:
import os
import argparse
import skimage.io
import numpy as np

def describe(im):
    print(f"\tdtype: {im.dtype}")
    print(f"\tshape: {im.shape}")
    print(f"\tmin: {im.min()}")
    print(f"\tmean: {im.mean()}")
    print(f"\tmax: {im.max()}")

def get_extention(filename):
    return filename.split('.')[-1]

def getFileType(filename):
    return filename.split('.')[0].split('_')[-1]

def isLabel(filename):
    return getFileType(filename) == 'label'

def isRaw(filename):
    return getFileType(filename) == 'raw'

def getFileName(filename):
    return filename.split('_')[0]

def find_matching_files(raw_dir, label_dir, file_ext = None):
    files =[]
    for raw_file in os.listdir(raw_dir):
        raw_name = getFileName(raw_file)
        if file_ext and file_ext != get_extention(raw_file):
            continue
        for label_file in os.listdir(label_dir):
            if file_ext and file_ext != get_extention(label_file):
                continue
            label_name = getFileName(label_file)
            if raw_name == label_name:
                files.append((raw_file, label_file))
    return files

if __name__ == '__main__':
    # parse arguments
    parser = argparse.ArgumentParser()
    parser.add_argument('--label_dir', type=str)
    parser.add_argument('--raw_dir', type=str)
    parser.add_argument('--dst_dir', type=str)
    parser.add_argument('--max_slice', type=int, default= -1 )
    config = parser.parse_args()

    # Assert all exist
    assert os.path.isdir(config.label_dir), "label dir must exist"
    assert os.path.isdir(config.raw_dir), "Raw dir must exist"
    assert os.path.isdir(config.dst_dir), "Dst dir must exist"

    # find matching names in label and raw path
    matched_files = find_matching_files(config.raw_dir, config.label_dir, file_ext='tif')

    # Go through files in label
    for raw_file,label_file in matched_files:
        ## Open the file
        raw = skimage.io.imread(os.path.join(config.raw_dir, raw_file))
        label = skimage.io.imread(os.path.join(config.label_dir, label_file))

        ## Continue if all label is empty (bg == 1)
        if np.all(label==1): continue
        ## Convert background to 0
        if label.min() == 1:
            label[label==1] = 0
            label[label==label.max()] = 1
        elif label.min() > 1:
            assert False

        ## Loop through slices
        for z_slice in range(len(label)):
            if z_slice == config.max_slice: break
            #if np.all(label[z_slice] == 0): continue
            ### Name slice
            new_raw_file = getFileName(raw_file) + '_slice'+str(z_slice).zfill(3) + '_raw.tif'
            new_label_file = getFileName(label_file) + '_slice'+str(z_slice).zfill(3) + '_label.tif'
            ### Save slice in dst
            skimage.io.imsave(os.path.join(config.dst_dir, new_raw_file), raw[z_slice: z_slice+1])
            skimage.io.imsave(os.path.join(config.dst_dir, new_label_file), label[z_slice: z_slice+1])

max_stars_repo_name: qin-yu/Plantseg-Cellpose-Augmentation

max_issues_repo_name: qin-yu/Plantseg-Cellpose-Augmentation

max_forks_repo_name: qin-yu/Plantseg-Cellpose-Augmentation


text:
# %%
# 除法 取的小数结果
3/2
# %%
# 除法 整除 最后结果去掉小数部分
3//2
# %%
# 分数表示

from fractions import Fraction
f = Fraction(38/32)
g = Fraction(38,32)
h = Fraction(3832.323)
print(f)
print(g)
print(h)
# %%
# 复数
i= 3+4j
# %%
from sympy import Symbol, expand, factor, solve
from sympy.plotting import ploi
# %%
x = Symbol('x')
# %%
x
# %%

ext: py

max_stars_repo_name: andyliumathematics/mlnotes

max_issues_repo_name: andyliumathematics/mlnotes

max_forks_repo_name: andyliumathematics/mlnotes

reason: from sympy


text:
from sklearn.svm import SVC
import numpy as np
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import confusion_matrix
from sklearn.model_selection import KFold

LAD_recordings = np.load('../data/LAD_train.npy')
print(len(LAD_recordings))
y = np.load('../data/labels_500.npy')
X_train, X_test, y_train, y_test = train_test_split(LAD_recordings, y, test_size=0.3, random_state=0)

clf = make_pipeline(StandardScaler(), SVC(gamma='auto'))
# clf.fit(X_train, y_train)
clf.fit(LAD_recordings, y)

# y_pred = clf.predict(X_test)
y_pred = clf.predict(LAD_recordings)

tn, fp, fn, tp = confusion_matrix(y, y_pred).ravel()
acc = (tn + tp) / (fp + fn + tn + tp)
print(acc)

scores = cross_val_score(clf, LAD_recordings, y, cv=5)

print(scores)

max_stars_repo_name: skywalkerzhang/ECG-classifier

max_issues_repo_name: skywalkerzhang/ECG-classifier

max_forks_repo_name: skywalkerzhang/ECG-classifier

reason: import numpy


text:
# -*- coding: utf-8 -*-
# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright Holders: Rene Milk, Stephan Rave, Felix Schindler
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)

from __future__ import absolute_import, division, print_function

import numpy as np
from types import FunctionType, MethodType

from pymor.core import BasicInterface
from pymor.core.pickle import dumps, loads, dumps_function, PicklingError
from pymor.grids.subgrid import SubGrid
from pymor.operators.basic import NumpyMatrixBasedOperator

is_equal_ignored_attributes = \
    ((SubGrid, {'_uid', '_CacheableInterface__cache_region', '_SubGrid__parent_grid'}),
     (NumpyMatrixBasedOperator, {'_uid', '_CacheableInterface__cache_region', '_assembled_operator'}),
     (BasicInterface, {'_uid', '_CacheableInterface__cache_region'}))

is_equal_dispatch_table = {}

def func_with_closure_generator():
    x = 42
    def func_with_closure():
        print(x)
    return func_with_closure

cell_type = type(func_with_closure_generator().func_closure[0])

def assert_is_equal(first, second):

    seen = []

    def _assert_is_equal(first, second):

        for x in seen:
            if x is first:
                return

        seen.append(first)

        for c, m in is_equal_dispatch_table.iteritems():
            if type(first) == c:
                assert m(first, second)
                return

        ignored_attributes = set()
        for c, v in is_equal_ignored_attributes:
            if isinstance(first, c):
                ignored_attributes = v
                break

        assert type(first) == type(second)

        if isinstance(first, np.ndarray):
            assert np.all(first == second)
        elif isinstance(first, (list, tuple)):
            assert len(first) == len(second)
            for u, v in zip(first, second):
                _assert_is_equal(u, v)
        elif isinstance(first, dict):
            assert set(first.keys()) == set(second.keys())
            for k, u in first.iteritems():
                _assert_is_equal(u, second.get(k))
        elif isinstance(first, FunctionType):
            for k in ['func_closure', 'func_code', 'func_dict', 'func_doc', 'func_name']:
                _assert_is_equal(getattr(first, k), getattr(second, k))
        elif isinstance(first, MethodType):
            _assert_is_equal(first.im_func, second.im_func)
            _assert_is_equal(first.im_self, second.im_self)
        elif isinstance(first, cell_type):
            _assert_is_equal(first.cell_contents, second.cell_contents)
        elif not isinstance(first, BasicInterface):
            assert first == second
        else:
            assert (set(first.__dict__.keys()) - ignored_attributes) == (set(second.__dict__.keys()) - ignored_attributes)
            for k, v in first.__dict__.iteritems():
                if not k in ignored_attributes:
                    _assert_is_equal(v, second.__dict__.get(k))

    _assert_is_equal(first, second)

def assert_picklable(o):
    s = dumps(o)
    o2 = loads(s)
    assert_is_equal(o, o2)

def assert_picklable_without_dumps_function(o):

    def dumps_function_raise(function):
        raise PicklingError('Cannot pickle function {}'.format(function))

    old_code = dumps_function.func_code
    dumps_function.func_code = dumps_function_raise.func_code
    try:
        s = dumps(o)
        o2 = loads(s)
        assert_is_equal(o, o2)
    finally:
        dumps_function.func_code = old_code

max_stars_repo_name: ftalbrecht/pymor

max_issues_repo_name: ftalbrecht/pymor

max_forks_repo_name: ftalbrecht/pymor


text:
from __future__ import division, print_function, absolute_import

import mock
import pytest
import importlib
import numpy as np
import scipy.spatial.distance as sci_dist

from numpy.testing import assert_almost_equal, assert_raises

from flithic.distance import distance

@pytest.fixture
def u():
    return np.random.rand(100,1)

@pytest.fixture
def v():
    return np.random.rand(100,1)

@pytest.fixture
def v_nan(v):
    v_nan = v.copy()
    v_nan[:20] = np.nan
    return v_nan

@pytest.fixture
def test_bivariate(funcname, u, v, v_nan, *args, **kwargs):
    # get functions
    f_scipy = getattr(sci_dist, funcname)
    f_flithic = getattr(distance, funcname)

    # test need for function
    assert( np.isnan(f_scipy(u, v_nan, *args)).any() )

    # test funcs are equivalent
    scipy_ = f_scipy(u, v, *args, **kwargs)
    flithic_ = f_flithic.__wrapped__(u, v, *args, **kwargs)

    assert_almost_equal(scipy_, flithic_)

    # test nan version works like we think
    scipy_nan = f_flithic(u[20:], v[20:], *args, **kwargs)
    flithic_nan = f_flithic(u, v_nan, *args, **kwargs)

    assert_almost_equal( scipy_nan, flithic_nan)

@pytest.fixture
def check_bivariate(funcname, u, v, v_nan, **kwargs):
    # get functions
    f_scipy = getattr(sci_dist, funcname)

    assert_almost_equal( f_scipy(u,v), f_scipy(u, v_nan))

def get_vnan_bool(v_nan):
    v_nan_bool = v_nan.copy()
    v_nan_bool[~np.isnan(v_nan)] = v_nan_bool[~np.isnan(v_nan)] > 0.5

    return v_nan_bool

# =============================================================

# TODO
#def test_directed_hausdorff(u, v, v_nan):
    #test_bivariate("directed_hausdorff", u, v, v_nan)

# =============================================================

# TODO
#def test_minkowski(u, v, v_nan):
#

def test_sqeuclidean(u, v, v_nan):
    test_bivariate("sqeuclidean", u, v, v_nan)

def test_correlation(u, v, v_nan):
    test_bivariate("correlation", u, v, v_nan)

def test_cosine(u, v, v_nan):
    test_bivariate("cosine", u, v, v_nan)

def test_hamming(u, v, v_nan):
    check_bivariate("hamming", u, v, v_nan)

def test_jaccard(u, v, v_nan):
    check_bivariate("jaccard", u, v, v_nan)

def test_kulsinski(u, v, v_nan):
    test_bivariate("kulsinski", u, v, v_nan)

def test_seuclidean(u, v, v_nan):
    V = np.random.rand(*v.shape)
    test_bivariate("seuclidean", u, v, v_nan, V)

def test_citybock(u, v, v_nan):
    test_bivariate("cityblock", u, v, v_nan)

def test_mahalanobis(u, v, v_nan):
    VI = np.random.rand(u.shape[0], v.shape[0])
    test_bivariate("mahalanobis", u, v, v_nan, VI)

def test_chebyshev(u, v, v_nan):
    test_bivariate("chebyshev", u, v, v_nan)

def test_braycurtis(u, v, v_nan):
    test_bivariate("braycurtis", u, v, v_nan)

def test_canberra(u, v, v_nan):
    check_bivariate("canberra", u, v, v_nan)
    #test_bivariate("canberra", u, v, v_nan)

def test_yule(u, v, v_nan):
    test_bivariate("yule", u, v, v_nan)

def test_rogerstanimoto(u, v, v_nan):
    test_bivariate("rogerstanimoto", u, v, v_nan)

def test_russellrao(u, v, v_nan):
    test_bivariate("russellrao", u, v, v_nan)

def test_sokalmichener(u, v, v_nan):
    test_bivariate("sokalmichener", u, v, v_nan)

def test_sokalsneath(u, v, v_nan):
    test_bivariate("sokalsneath", u, v, v_nan)
    test_bivariate("sokalsneath", u > 0.5, v > 0.5, get_vnan_bool(v_nan))

# =============================================================

reason: import numpy,from numpy,import scipy


text:
"""
A module which contains the routines needed for computing the spillage values.
"""

#: Conversion between Atomic Units and Bohr
AU_to_A = 0.52917721092

class MatrixMetadata:
    """
    This class contains the information stored in the sparse matrix
    metadata file.

    Args:
      filename (str): the name of the metadata file to initialize from.

    Attributes:
      atoms (list): a list of atoms with their positions and associated
        basis functions.
      matdim (int): the dimension of the matrix.
    """
    def __init__(self, filename):
        from BigDFT.Atom import Atom
        self.atoms = []
        symbol_lookup = []
        with open(filename, "r") as ifile:
            # Read the first line
            matinfo = next(ifile).split()
            self.matdim, natoms, ntypes = [int(x) for x in matinfo[:3]]
            # Units
            next(ifile)
            # skip geocode
            line = next(ifile)
            # skip shift?
            line = next(ifile)
            # get the symbol lookup information
            for i in range(0, ntypes):
                line = next(ifile).split()
                symbol_lookup.append(line[2])
            # Get the atom positions
            for i in range(0, natoms):
                line = next(ifile).split()
                adict = {}
                adict["sym"] = symbol_lookup[int(line[0]) - 1]
                adict["r"] = [float(x) for x in line[1:4]]
                # THIS IS BECAUSE THE METADATA FILE CURRENTLY PRINTS THE
                # WRONG UNITS!
                adict["units"] = "bohr"
                adict["indices"] = []
                self.atoms.append(Atom(adict))
            # Get the indices
            for i in range(0, self.matdim):
                line = next(ifile).split()
                self.atoms[int(line[0]) - 1]["indices"].append(i)

    def get_frag_indices(self, system):
        """
        Retrive the indices of the matrix associated with each fragment
        of a system.

        Args:
          system (Fragments.System): the set of fragments to get the indices
            of.

        Return:
          (dict): a mapping from fragment to indices.

        Return:
          (dict): a mapping from fragment to indices.
        """
        matching = system.compute_matching(self.atoms)

        frag_indices = {}
        for fragid in system:
            frag_indices[fragid] = []
            for at in matching[fragid]:
                frag_indices[fragid] += self.atoms[at]["indices"]

        return frag_indices

def compute_spillage_values(sinvxh, sinvxh2, frag_indices, target):
    """
    Computes the actual spillage values.

    Args:
      sinvxh (scipy.sparse.csc): S^-1 * H
      sinvxh2 (scipy.sparse.csc): (S^-1 * H)^2
      frag_indices (dict): list of indices associated with each fragment.
      target (int): which fragment is the target fragment.

    Returns:
      (tuple): ``H2F``,``HGHF``, where the first contribution is the target
      value that the sum of the other contribution should reach.
      Such second term is the interaction spillage between the target and
      each fragment.
    """
    from numpy import trace
    indices_f = frag_indices[target]
    spillage_values = {}

    # Compute the denominator tr(HRfHRf)
    denom = sinvxh[:, indices_f]
    denom = denom[indices_f, :]
    denom = denom.dot(denom)
    denom_t = trace(denom.todense())

    # Compute the left side tr(HS-1HRf)
    H2T = sinvxh2[:, indices_f]
    H2T = H2T[indices_f, :]
    left_t = trace(H2T.todense())

    # Compute the right side values tr(HRgHRf)
    for id_g in frag_indices:
        indices_g = frag_indices[id_g]
        TFH = sinvxh[indices_f, :]
        TFHTG = TFH[:, indices_g]

        TGH = sinvxh[indices_g, :]
        TGHTF = TGH[:, indices_f]

        right_mat = (TFHTG.dot(TGHTF))
        right_t = trace(right_mat.todense())
        spillage_values[id_g] = right_t / denom_t

    return left_t / denom_t, spillage_values

def serial_compute_spillbase(sfile, hfile):
    """
    This routine computes the matrix (S^-1 * H)^2 using python.

    You will need this for computing the spillage values.

    Args:
      sfile (str): the file name of the overlap matrix. Must be in ccs format
      hfile (str): the file name of the hamiltonian. Must be in ccs format

     Returns:
       (scipy.sparse.csc, scipy.sparse.csc): (S^-1 * H), (S^-1 * H)^2
    """
    from scipy.sparse.linalg import inv

    # Read from file
    smat = _read_ccs(sfile)
    hmat = _read_ccs(hfile)

    # Compute the matrix (S^-1H)^2
    sinv = inv(smat)
    sinvxh = sinv.dot(hmat)
    sinvxh2 = sinvxh.dot(sinvxh)

    return sinvxh, sinvxh2

def serial_compute_puritybase(sfile, dfile):
    """
    This routine computes the matrix K*S using python.

    You will need this for computing the purity values.

    Args:
     sfile (str): the file name of the overlap matrix. Must be in ccs format
     dfile (str): the file name of the density kernel. Must be in ccs format

    Returns:
      (scipy.sparse.csc, scipy.sparse.csc): K*S
    """
    # Read from file
    smat = _read_ccs(sfile)
    dmat = _read_ccs(dfile)

    return dmat.dot(smat)

def _read_ccs(fname):
    """
    Read a CCS matrix from file into a scipy csc matrix.

    This routine uses the CCS matrix files that CheSS generates.
    Args:
      fname (str): name of the file.

    Result:
      scipy.sparse.csc_matrix: the matrix in the file.
    """
    from scipy.sparse import csc_matrix

    data = []
    indices = []
    indptr = []

    with open(fname, "r") as ifile:
        # Read the meta data
        line = next(ifile)
        split = line.split()
        matdim = int(split[0])
        nel = int(split[2])
        split = next(ifile).split()
        indptr = [int(x) - 1 for x in split]

        # Read the indices
        added = 0
        while(added < nel):
            split = next(ifile).split()
            values = [int(x) - 1 for x in split]
            indices.extend(values)
            added += len(values)

        # Read the data
        for line in ifile:
            data.append(float(line))

    matrix = csc_matrix((data, indices, indptr), shape=(matdim, matdim))
    return matrix

reason: from numpy,from scipy


text:
"""
Generate a video of an agent playing the game.
"""

import random

from anyrl.utils.ffmpeg import export_video
import numpy as np
import torch

from obs_tower2.model import ACModel
from obs_tower2.states import StateEnv
from obs_tower2.util import big_obs, create_single_env

def main():
    def image_fn():
        env = StateEnv(create_single_env(random.randrange(15, 20), clear=False))
        model = ACModel()
        tail_model = ACModel()
        model.load_state_dict(torch.load('save.pkl', map_location='cpu'))
        tail_model.load_state_dict(torch.load('save_tail.pkl', map_location='cpu'))
        model.to(torch.device('cuda'))
        tail_model.to(torch.device('cuda'))
        state, obs = env.reset()
        floor = 0
        while True:
            output = (model if floor < 10 else tail_model).step(np.array([state]), np.array([obs]))
            (state, obs), rew, done, info = env.step(output['actions'][0])
            floor = info['current_floor']
            yield big_obs(obs[..., -3:], info)
            if done:
                break
        env.close()

    export_video('export.mp4', 168, 168, 10, image_fn())

if __name__ == '__main__':
    main()


text:
import numpy as np

np.random.seed(0)
a = "computerscience"
b = "philosophy"

n_examples = 5
print("
")
for i in range(n_examples):
    s = ""
    # Just like merging in merge sort
    i1 = 0
    i2 = 0
    while i1 < len(a) and i2 < len(b):
        i12 = min(len(a), i1+np.random.randint(1, 4))
        i22 = min(len(b), i2+np.random.randint(1, 4))
        s += "" + a[i1:i12] + ""
        s += "" + b[i2:i22] + ""
        i1 = i12
        i2 = i22
    while i1 < len(a):
        i12 = min(len(a), i1+np.random.randint(4))
        s += "" + a[i1:i12] + ""
        i1 = i12
    while i2 < len(b):
        i22 = min(len(b), i2+np.random.randint(4))
        s += "" + b[i2:i22] + ""
        i2 = i22
    print("
", s, "
")
print("
")


text:
import os
import gym
import cv2
import numpy as np
from collections import deque

import game.wrapped_flappy_bird as game

from parameters import ENV, FRAME_SKIP, FRAME_BUFFER_SIZE

from PIL import Image
import imageio

def onehot(action):
    action_onehot = [0] * 2
    action_onehot[action] = 1
    return action_onehot

class Environment:

    def __init__(self, worker_index=0):

        self.frame_buffer = deque(maxlen=FRAME_BUFFER_SIZE)
        self.reset()

        self.render = False
        self.images = []
        self.tmp_name = 'tmp_{}.png'.format(worker_index)

    def get_state_size(self):
        return (84, 84, FRAME_BUFFER_SIZE)

    def get_action_size(self):
        return 2

    def set_render(self, render):
        self.render = render

    def _convert_process_buffer(self):
        # Convert RGB images into grayscale images
        state = map(lambda x: cv2.resize(cv2.cvtColor(x, cv2.COLOR_RGB2GRAY),
                                         (84, 90)),
                    self.frame_buffer)

        # Cut the images to get a 84x84 image
        state = map(lambda x: x[1:85, :, np.newaxis], state)

        # Concatenate the frames to get a single state
        return np.concatenate(list(state), axis=2)

    def reset(self):
        self.env = game.GameState(1, False)

        # Reset the frame buffer with FRAME_BUFFER_SIZE frames
        for _ in range(FRAME_BUFFER_SIZE):
            frame, r, done = self.env.frame_step(onehot(0))
            self.frame_buffer.append(frame)

        return self._convert_process_buffer()

    def act(self, action, gif=False):
        if not gif:
            return self._act(action)
        else:
            return self._act_gif(action)

    def _act(self, action):

        s, r, done = self.env.frame_step(onehot(action))
        self.frame_buffer.append(s)
        return self._convert_process_buffer(), r, done, ""

    def _act_gif(self, action):

        r = 0
        i, done = 0, False
        while i < (FRAME_SKIP + 1) and not done:

            s_, r_tmp, done = self.env.frame_step(onehot(action))
            r += r_tmp
            i += 1

            # Save image
            img = Image.fromarray(s_)
            img.save(self.tmp_name)
            self.images.append(imageio.imread(self.tmp_name))

        self.frame_buffer.append(s_)
        return self._convert_process_buffer(), r, done, ""

    def save_gif(self, path, i):
        path = path + "_{}.gif".format(i)
        os.makedirs(os.path.dirname(path), exist_ok=True)
        imageio.mimsave(path, self.images, duration=1)
        self.images = []

    def close(self):
        pass

max_stars_repo_name: SuReLI/Deep-RL-agents

max_issues_repo_name: SuReLI/Deep-RL-agents

max_forks_repo_name: SuReLI/Deep-RL-agents


text:
# -*- coding: iso-8859-15 -*-
"""
CALFEM Geometry module

Contains functions and classes for describing geometry.
"""

import numpy as np

class Geometry:
    '''
    Instances of GeoData can hold geometric data and be passed to
    GmshMesher in pycalfem_Mesh to mesh the geometry.
    '''

    def __init__(self):
        self.points = {}  # dict of [[x, y, z], elSize, marker]
        # dict of [curvTypestring, [p1, p2, ... pn], marker, elementsOnCurve, distributionString, distributionVal]
        self.curves = {}
        # dict of [SurfaceTypeString, [c1, c2 ... cn], [[c1, c2 ... cm], ... [c1, ... ck]], ID, marker, is_structured]. c means curve-ID.
        self.surfaces = {}
        # dict of [[s1, s2 ..], [[s1,s2...],[s1,s2..],..], ID, marker, is_structured]
        self.volumes = {}
        # This is automatically set to True if a 3D point is added.
        self.is3D = False
        self._pointIDspecified = False
        self._nextPointID = 0
        self._curveIDspecified = False
        self._nextcurveID = 0
        self._surfaceIDspecified = False
        self._nextsurfaceID = 0
        self._volumeIDspecified = False
        self._nextvolumeID = 0

    def removePoint(self, ID):
        '''Removes the point with this ID'''
        self.points.pop(ID)

    def removeCurve(self, ID):
        '''Removes the curve with this ID'''
        self.curves.pop(ID)

    def removeSurface(self, ID):
        '''Removes the surface with this ID'''
        self.surfaces.pop(ID)

    def removeVolume(self, ID):
        '''Removes the volume with this ID'''
        self.volumes.pop(ID)

    def getPointCoords(self, IDs=None):
        '''
        Returns an N-by-3 list of point coordinates if the parameter is
        a list of IDs. If the parameter is just a single integer then
        a single coordinate (simple 3-element list) is returned.
        If the parameter is undefined (or None) all point coords will be returned
        '''
        if IDs == None:
            return [p[0] for p in self.points.values()]
        try:
            pointCoords = [self.points[pID][0] for pID in IDs]
        except TypeError:  # IDs was not iterable. Probably just a single ID.
            pointCoords = self.points[IDs][0]
        return pointCoords

    def pointsOnCurves(self, IDs):
        '''
        Returns a list of all geometric points (not nodes) on the curves
        specified in IDs. IDs may be an integer or a list of integers.
        '''
        return self._subentitiesOnEntities(IDs, self.curves, 1)

    def stuffOnSurfaces(self, IDs):
        '''
        Returns lists of all geometric points and curves on the surfaces
        specified in IDs. IDs may be an integer or a list of integers
        '''
        curveSet = self._subentitiesOnEntities(
            IDs, self.surfaces, 1)  # Curves on the outer edges
        curveSet.update(self._subentityHolesOnEntities(
            IDs, self.surfaces, 2))  # Curves on the holes
        # Points on the curves of these surfaces.
        pointList = self.pointsOnCurves(curveSet)
        return pointList, list(curveSet)

    def stuffOnVolumes(self, IDs):
        '''
        Returns lists of all geometric points, curves, and surfaces on the volumes
        specified in IDs. IDs may be an integer or a list of integers
        '''
        surfaceSet = self._subentitiesOnEntities(IDs, self.surfaces, 0)
        surfaceSet.update(self._subentitiesOnEntities(IDs, self.surfaces, 1))
        pointList, curveList = self.stuffOnSurfaces(surfaceSet)
        return pointList, curveList, list(surfaceSet)

    def _subentitiesOnEntities(self, IDs, entityDict, index):
        '''
        Duplicate code. Gets the IDs of the subentities that
        make up an entity, i.e. the points that define a curve or
        the curves that define a surface. Note that only the outer
        subentities of surfaces and volumes can be extracted with
        this function. For holes use _subentityHolesOnEntities().
        '''
        theSet = set()
        try:
            for ID in IDs:
                theSet.update(entityDict[ID][index])
        except TypeError:  # IDs is not iterable, so it is probably a single ID
            theSet.update(entityDict[IDs][index])
        return theSet

    def _subentityHolesOnEntities(self, IDs, entityDict, index):
        '''Duplicate code. Does the same thing as _subentitiesOnEntities(), but for holes'''
        theSet = set()
        try:
            for ID in IDs:
                for hole in entityDict[ID][index]:
                    theSet.update(hole)
        except TypeError:  # IDs is not iterable, so it is probably a single ID
            for hole in entityDict[IDs][index]:
                theSet.update(hole)
        return theSet

    def addPoints(self, points, markers=None, ids=None, elSizes=None):
        '''
        Add points from a numpy-array
        '''
        nPoints, dims = points.shape

        if dims == 2:
            for row in points:
                self.addPoint(row.tolist())
        elif dims == 3:
            for row in points:
                self.addPoint(row.tolist())

    def addPoint(self, coord, ID=None, marker=0, el_size=1):
        '''
        Adds a point.

        Parameters:
        coord     - [x, y] or [x, y, z].
                    List, not array.

        ID        - Positive integer ID of this point. If left unspecified the
                    point will be assigned the smallest unused point-ID.
                    It is recommended to specify all point-IDs or none.

        marker    - Marker applied to this point. Default 0.
                    It is not a good idea to apply non-zero markers to points
                    that are control points on B-splines or center points on
                    circles/ellipses, since this can lead to "loose" nodes
                    that are not part of any elements.

        elSize    - The size of elements at this point. Default 1. Use to make
                    a mesh denser or sparser here. Only affects unstructured
                    meshes
        '''
        if len(coord) == 3:  # A 3D point is inserted.
            self.is3D = True
        else:  # The point is in 2D (we assume)
            # Pad with a 0. (we store points as 3D points for consistency's sake)
            coord = coord+[0]

        if ID == None:  # ID is not specified. Get an ID for this point:
            ID = self._getNewPointID()
        else:
            self._pointIDspecified = True

        self.points[ID] = [coord, el_size, marker]

    def bounding_box_2d(self):
        """Calculate bounding box geometry"""

        min_x = 1e300
        max_x = -1e300
        min_y = 1e300
        max_y = -1e300

        for point_id in self.points.keys():
            p = self.points[point_id]
            if p[0][0] > max_x:
                max_x = p[0][0]
            if p[0][0] < min_x:
                min_x = p[0][0]
            if p[0][1] > max_y:
                max_y = p[0][1]
            if p[0][1] < min_y:
                min_y = p[0][1]

        return min_x, max_x, min_y, max_y

    def addSplines(self, points):
        '''
        Add splines from numpy array
        '''
        nPoints, dims = points.shape

        if dims == 2:
            for row in points:
                self.addSpline(row.tolist())
        elif dims == 3:
            for row in points:
                splineDef = row.tolist()
                self.addSpline(splineDef[:-1], marker=splineDef[2])

    def addSpline(self, points, ID=None, marker=0, el_on_curve=None, el_distrib_type=None, el_distrib_val=None):
        '''
        Adds a Spline curve

        points    - List of indices of control points that make a Spline
                    [p1, p2, ... , pn]

        ID        - Positive integer ID of this curve. If left unspecified the
                    curve will be assigned the smallest unused curve-ID.
                    It is recommended to specify all curve-IDs or none.

        marker    - Integer. Marker applied to this curve. Default 0.

        elOnCurv  - Positive integer. Elements on curve.
                    The number of element edges that will be distributed
                    along this curve. Only works for structured meshes.

        el_distrib_type -
                    String. Either "bump" or "progression".
                    Determines how the density of elements vary along the curve
                    for structured meshes. Only works for structured meshes.
                    elOnCurv and el_distrib_val must be be defined if this param
                    is used.

        el_distrib_val -
                    Float. Determines how severe the element distribution is.
                    Only works for structured meshes. elOnCurv and
                    el_distrib_type must be be defined if this param is used.

                        bump:
                    Smaller value means elements are bunched up at the edges
                    of the curve, larger means bunched in the middle.

                        progression:
                    The edge of each element along this curve (from starting
                    point to end) will be larger than the preceding one by
                    this factor.
                    el_distrib_val = 2 meaning for example that each line element
                    in the series will be twice as long as the preceding one.
                    el_distrib_val < 1 makes each element smaller than the
                    preceeding one.
        '''
        self._addCurve("Spline", points, ID, marker, el_on_curve,
                       el_distrib_type, el_distrib_val)

    def addBSpline(self, points, ID=None, marker=0, el_on_curve=None,  el_distrib_type=None, el_distrib_val=None):
        '''
        Adds a B-Spline curve

        points    - List of indices of control points that make a B-spline
                    [p1, p2, ... , pn]

        ID        - Positive integer ID of this curve. If left unspecified the
                    curve will be assigned the smallest unused curve-ID.
                    It is recommended to specify all curve-IDs or none.

        marker    - Integer. Marker applied to this curve. Default 0.

        elOnCurv  - Positive integer. Elements on curve.
                    The number of element edges that will be distributed
                    along this curve. Only works for structured meshes.

        el_distrib_type -
                    String. Either "bump" or "progression".
                    Determines how the density of elements vary along the curve
                    for structured meshes. Only works for structured meshes.
                    elOnCurv and el_distrib_val must be be defined if this param
                    is used.

        el_distrib_val -
                    Float. Determines how severe the element distribution is.
                    Only works for structured meshes. elOnCurv and
                    el_distrib_type must be be defined if this param is used.

                        bump:
                    Smaller value means elements are bunched up at the edges
                    of the curve, larger means bunched in the middle.

                        progression:
                    The edge of each element along this curve (from starting
                    point to end) will be larger than the preceding one by
                    this factor.
                    el_distrib_val = 2 meaning for example that each line element
                    in the series will be twice as long as the preceding one.
                    el_distrib_val < 1 makes each element smaller than the
                    preceeding one.
        '''
        self._addCurve("BSpline", points, ID, marker,
                       el_on_curve, el_distrib_type, el_distrib_val)

    def addCircle(self, points, ID=None, marker=0, el_on_curve=None, el_distrib_type=None, el_distrib_val=None):
        '''
        Adds a Circle arc curve.

        points    - list of 3 indices of point that make a circle arc smaller
                    than Pi.
                    [startpoint, centerpoint, endpoint]

        ID        - Positive integer ID of this curve. If left unspecified the
                    curve will be assigned the smallest unused curve-ID.
                    It is recommended to specify all curve-IDs or none.

        marker    - Marker applied to this curve. Default 0.

        elOnCurv  - Elements on curve.
                    The number of element edges that will be distributed
                    along this curve. Only works for structured meshes.

        el_distrib_type -
                    String. Either "bump" or "progression".
                    Determines how the density of elements vary along the curve
                    for structured meshes. Only works for structured meshes.
                    elOnCurv and el_distrib_val must be be defined if this param
                    is used.

        el_distrib_val -
                    Float. Determines how severe the element distribution is.
                    Only works for structured meshes. elOnCurv and
                    el_distrib_type must be be defined if this param is used.

                        bump:
                    Smaller value means elements are bunched up at the edges
                    of the curve, larger means bunched in the middle.

                        progression:
                    The edge of each element along this curve (from starting
                    point to end) will be larger than the preceding one by
                    this factor.
                    el_distrib_val = 2 meaning for example that each line element
                    in the series will be twice as long as the preceding one.
                    el_distrib_val < 1 makes each element smaller than the
                    preceeding one.
        '''
        if len(points) != 3:
            raise IndexError(
                "Circle: points must be a list of 3 positive integers denoting point indices")
        self._addCurve("Circle", points, ID, marker, el_on_curve,
                       el_distrib_type, el_distrib_type)

    def addEllipse(self, points, ID=None, marker=0, el_on_curve=None, el_distrib_type=None, el_distrib_val=None):
        '''
        Adds a Ellipse arc curve.

        points    - List of 4 indices of point that make a ellipse arc smaller
                    than Pi.
                    [startpoint, centerpoint, mAxisPoint, endpoint]
                    Startpoint is the starting point of the arc.
                    Centerpoint is the point at the center of the ellipse.
                    MAxisPoint is any point on the major axis of the ellipse.
                    Endpoint is the end point of the arc.

        ID        - Positive integer ID of this curve. If left unspecified the
                    curve will be assigned the smallest unused curve-ID.
                    It is recommended to specify all curve-IDs or none.

        marker    - Integer. Marker applied to this curve. Default 0.

        elOnCurv  - Positive integer. Elements on curve.
                    The number of element edges that will be distributed
                    along this curve. Only works for structured meshes.

        el_distrib_type -
                    String. Either "bump" or "progression".
                    Determines how the density of elements vary along the curve
                    for structured meshes. Only works for structured meshes.
                    elOnCurv and el_distrib_val must be be defined if this param
                    is used.

        el_distrib_val -
                    Float. Determines how severe the element distribution is.
                    Only works for structured meshes. elOnCurv and
                    el_distrib_type must be be defined if this param is used.

                        bump:
                    Smaller value means elements are bunched up at the edges
                    of the curve, larger means bunched in the middle.

                        progression:
                    The edge of each element along this curve (from starting
                    point to end) will be larger than the preceding one by
                    this factor.
                    el_distrib_val = 2 meaning for example that each line element
                    in the series will be twice as long as the preceding one.
                    el_distrib_val < 1 makes each element smaller than the
                    preceeding one.
        '''
        if len(points) != 4:
            raise IndexError(
                "Ellipse: points must be a list of 4 positive integers denoting point indices")
        self._addCurve("Ellipse", points, ID, marker,
                       el_on_curve, el_distrib_type, el_distrib_val)

    def _addCurve(self, name, points, ID, marker, el_on_curve, el_distrib_type, el_distrib_val):
        '''Duplicate code goes here!'''
        if ID == None:
            ID = self._getNewCurveID()
        else:
            self._curveIDspecified = True

        if el_distrib_type != None:
            # transform into lowercase except the first letter upper case.
            el_distrib_type = el_distrib_type.lower().title()
            if el_distrib_type not in ["Bump", "Progression"]:
                raise ValueError(
                    "el_distrib_type must be a string, either \"bump\" or \"progression\". Curve with ID=%i was incorrect." % ID)
            if el_distrib_val == None:
                raise ValueError(
                    "If el_distrib_type is defined then el_distrib_val must be given a (float) value")

        self.curves[ID] = [name, points, marker,
                           el_on_curve, el_distrib_type, el_distrib_val]

    def addSurface(self, outer_loop, holes=[], ID=None, marker=0):
        '''
        Adds a plane surface (flat).
        Parameters:
        outer_loop - List of curve IDs that make up the outer boundary of
                    the surface. The curves must lie in the same plane.

        holes     - List of lists of curve IDs that make up the inner
                    boundaries of the surface. The curves must lie in the
                    same plane.

        ID        - Positive integer ID of this surface. If left unspecified
                    the surface will be assigned the smallest unused surface-ID.
                    It is recommended to specify all surface-IDs or none.

        marker    - Integer. Marker applied to this surface. Default 0.
        '''
        # TODO: Possibly check if outer_loop is an actual loop and if the holes are correct.
        self._addSurf("Plane Surface", outer_loop, holes,
                      ID, marker, is_structured=False)

    def addRuledSurface(self, outer_loop, ID=None, marker=0):
        '''
        Adds a Ruled Surface (bent surface).
        Parameters:
        outer_loop - List of 3 or 4 curve IDs that make up the boundary of
                    the surface.

        ID        - Positive integer ID of this surface. If left unspecified
                    the surface will be assigned the smallest unused surface-ID.
                    It is recommended to specify all surface-IDs or none.

        marker    - Integer. Marker applied to this surface. Default 0.
        '''
        if len(outer_loop) not in [3, 4]:
            raise IndexError(
                "Ruled Surface: outer_loop must be a list of 3 or 4 positive integers denoting curve indices")
        self._addSurf("Ruled Surface", outer_loop, [],
                      ID, marker, is_structured=False)

    def addStructuredSurface(self, outer_loop, ID=None, marker=0):
        '''
        Adds a Structured Surface.
        Parameters:
        outer_loop - List of 4 curve IDs that make up the boundary of
                    the surface. The curves must be structured, i.e. their
                    parameter 'elOnCurv' must be defined.

        ID        - Positive integer ID of this surface. If left unspecified
                    the surface will be assigned the smallest unused surface-ID.
                    It is recommended to specify all surface-IDs or none.

        marker    - Integer. Marker applied to this surface. Default 0.
        '''
        self._checkIfProperStructuredQuadBoundary(outer_loop, ID)
        self._addSurf("Ruled Surface", outer_loop, [],
                      ID, marker, is_structured=True)

    def _addSurf(self, name, outer_loop, holes, ID, marker, is_structured):
        '''For duplicate code'''
        # TODO: check if the curves in outer_loop actually exist. Maybe print a warning.
        if ID == None:
            ID = self._getNewSurfaceID()
        else:
            self._surfaceIDspecified = True

        # Catch the easy mistake of making holes a list of ints rather than a list of lists of ints.
        for hole in holes:
            try:
                [h for h in hole]
            except TypeError:
                raise TypeError(
                    "Hole " + str(hole) + " is not iterable. Parameter holes must be a list of lists of integers")

        self.surfaces[ID] = [name, outer_loop,
                             holes, ID, marker, is_structured]

    def addVolume(self, outer_surfaces, holes=[], ID=None, marker=0):
        '''Adds a Volume
        Parameters:
        outer_surfaces - List of surface IDs that make up the outer boundary of
                        the volume.

        holes         - List of lists of surface IDs that make up the inner
                        boundaries of the volume.

        ID            - Positive integer ID of this volume. If left unspecified
                        the volume will be assigned the smallest unused volume-ID.
                        It is recommended to specify all volume-IDs or none.

        marker        - Integer. Marker applied to this volume. Default 0.'''
        self._addVolume(outer_surfaces, holes, ID, marker, is_structured=False)

    def addStructuredVolume(self, outer_surfaces, ID=None, marker=0):
        '''Adds a Structured Volume
        Parameters:
        outer_surfaces - List of surface IDs that make up the outer boundary of
                        the volume. The surfaces must be Structured Surfaces.

        ID            - Positive integer ID of this volume. If left unspecified
                        the volume will be assigned the smallest unused volume-ID.
                        It is recommended to specify all volume-IDs or none.

        marker        - Integer. Marker applied to this volume. Default 0.'''
        # TODO: Check input. (see if surfaces are structured)
        self._addVolume(outer_surfaces, [], ID, marker, is_structured=True)

    def _addVolume(self,  outer_surfaces, holes, ID, marker, is_structured):
        '''Duplicate code'''
        # TODO: Check input (outer_surfaces and holes[i] must be closed surfaces)
        if ID == None:
            ID = self._getNewVolumeID()
        else:
            self._volumeIDspecified = True
        self.volumes[ID] = [outer_surfaces, holes, ID, marker, is_structured]

    def pointMarker(self, ID, marker):
        self.setPointMarker(ID, marker)

    def setPointMarker(self, ID, marker):
        '''Sets the marker of the point with the ID'''
        self.points[ID][2] = marker

    def curveMarker(self, ID, marker):
        self.setCurveMarker(ID, marker)

    def setCurveMarker(self, ID, marker):
        '''Sets the marker of the curve with the ID'''
        self.curves[ID][2] = marker

    def surfaceMarker(self, ID, marker):
        self.setSurfaceMarker(ID, marker)

    def setSurfaceMarker(self, ID, marker):
        '''Sets the marker of the surface with the ID'''
        self.surfaces[ID][4] = marker

    def setVolumeMarker(self, ID, marker):
        '''Sets the marker of the volume with the ID'''
        self.volumes[ID][3] = marker

    def _checkIfProperStructuredQuadBoundary(self, outer_loop, ID):
        '''Checks if the four edges of a quad-shaped superelement exist and
        are correct, i.e el_on_curve of opposite curves are equal.'''
        if len(outer_loop) != 4:
            raise IndexError(
                "Structured Surface: outer_loop must be a list of 4 positive integers denoting curve indices")

        try:
            c0 = self.curves[outer_loop[0]]
            c1 = self.curves[outer_loop[1]]
            c2 = self.curves[outer_loop[2]]
            c3 = self.curves[outer_loop[3]]
        except KeyError:
            raise KeyError(
                "Structured Surface: Attempted construction of StructuredSurface with ID=%s from a curve that does not exist" % ID)

        if None in [c0, c1, c2, c3]:
            raise Exception(
                "Attempted to create structured surface from non-structured boundary curves.")

        # Check if the number of elements on opposite curves match.
        if(c0[-3] != c2[-3] or c1[-3] != c3[-3]):
            raise Exception("Structured Surface: The outer_loop of StructuredSurface %i is not properly " +
                            "constructed. The reason could be that the number of elements (elOnCurv) on " +
                            "opposite pairs of curves are different")

    def _getNewPointID(self):
        if not self._pointIDspecified:
            self._nextPointID += 1
            return self._nextPointID - 1
        else:
            return self._smallestFreeKey(self.points)

    def _getNewCurveID(self):
        if not self._curveIDspecified:
            self._nextcurveID += 1
            return self._nextcurveID - 1
        else:
            return self._smallestFreeKey(self.curves)

    def _getNewSurfaceID(self):
        if not self._surfaceIDspecified:
            self._nextsurfaceID += 1
            return self._nextsurfaceID - 1
        else:
            return self._smallestFreeKey(self.surfaces)

    def _getNewVolumeID(self):
        if not self._volumeIDspecified:
            self._nextvolumeID += 1
            return self._nextvolumeID - 1
        else:
            return self._smallestFreeKey(self.volumes)

    def _smallestFreeKey(self, dictionary):
        '''Finds the smallest unused key in the dict.'''
        sortedkeys = sorted(dictionary)
        for i in range(len(dictionary)):
            if sortedkeys[i] != i:
                return i

    point = addPoint
    spline = addSpline
    line = addSpline
    bspline = addBSpline
    circle = addCircle
    ellipse = addEllipse
    surface = addSurface
    surf = addSurface
    ruledSurface = addRuledSurface
    ruled_surface = addRuledSurface
    ruled_surf = addRuledSurface
    structuredSurface = addStructuredSurface
    structured_surface = addStructuredSurface
    struct_surf = addStructuredSurface
    volume = addVolume
    structuredVolume = addStructuredVolume
    structured_volume = addStructuredVolume
    get_point_coords = getPointCoords
    curve_marker = curveMarker

def geometry():
    return Geometry()

max_stars_repo_name: CALFEM/calfem-py

max_issues_repo_name: CALFEM/calfem-py

max_forks_repo_name: CALFEM/calfem-py

reason: import numpy,from numpy


text:
import matplotlib.pyplot as plt
import numpy
import os
import sys
sys.path.append(r'/home/jorge/ParaView-5.8.0-RC3-MPI-Linux-Python3.7-64bit/lib/python3.7/site-packages/')
from paraview.simple import *
from subprocess import check_output

cwd_base = os.getcwd().rsplit(sep = os.sep, maxsplit = 2)
cwd = os.path.join(cwd_base[0], 'results','')
stdout = check_output('ls' +' {}'.format(cwd), shell=True)
files = stdout.decode().split(sep='\n')
reader = OpenDataFile([cwd + a for a in files])
print(reader.PointData.keys())
print(reader.PointData.values())
p_vel = reader.PointData["Proton - Solar wind-velocity"]
print(type(p_vel))
print(type(p_vel.FieldData.GetArray(7)))
#help(reader)

max_stars_repo_name: jagarciap/SCSI

max_issues_repo_name: jagarciap/SCSI

max_forks_repo_name: jagarciap/SCSI


text:
import numpy as np

from meegkit import dss, sns, tspca
from meegkit.utils import demean, fold, unfold

import matplotlib.pyplot as plt

def test_tspca_sns_dss():  # TODO
    """Test TSPCA, SNS, DSS.

    Requires data stored in a time X channels X trials matrix.

    Remove environmental noise with TSPCA (shifts=-50:50).
    Remove sensor noise with SNS.
    Remove non-repeatable components with DSS.
    """
    # Random data (time*chans*trials)
    data = np.random.random((800, 102, 200))
    ref = np.random.random((800, 3, 200))

    # remove means
    noisy_data = demean(data)
    noisy_ref = demean(ref)

    # Apply TSPCA
    # -------------------------------------------------------------------------
    # shifts = np.arange(-50, 51)
    # print('TSPCA...')
    # y_tspca, idx = tspca.tsr(noisy_data, noisy_ref, shifts)[0:2]
    # print('\b OK!')
    y_tspca = noisy_data

    # Apply SNS
    # -------------------------------------------------------------------------
    nneighbors = 10
    print('SNS...')
    y_tspca_sns, r = sns.sns(y_tspca, nneighbors)
    print('\b OK!')

    # apply DSS
    # -------------------------------------------------------------------------
    print('DSS...')
    # Keep all PC components
    y_tspca_sns = demean(y_tspca_sns)
    print(y_tspca_sns.shape)
    todss, fromdss, _, _ = dss.dss1(y_tspca_sns)
    print('\b OK!')

    # c3 = DSS components
    y_tspca_sns_dss = fold(
        np.dot(unfold(y_tspca_sns), todss), y_tspca_sns.shape[0])

    return y_tspca, y_tspca_sns, y_tspca_sns_dss

def test_tsr(show=True):
    """Test time-shift regression."""
    sr = 200
    nsamples = 10000
    nchans = 10
    x = np.random.randn(nsamples, nchans)

    # artifact + harmonics
    artifact = np.sin(np.arange(nsamples) / sr * 2 * np.pi * 10)[:, None]
    artifact[artifact < 0] = 0
    artifact = 5 * artifact ** 3
    signal = x + artifact

    signal -= np.mean(signal, keepdims=True)
    artifact -= np.mean(artifact, keepdims=True)

    # Without shifts
    y, idx, mean_total, weights = tspca.tsr(
        signal,
        artifact,
        shifts=[0])

    if show:
        f, ax = plt.subplots(2, 1, num='without shifts')
        ax[0].plot(y[:500, 0], 'grey', label='recovered signal')
        ax[0].plot(x[:500, 0], ':', label='real signal')
        ax[1].plot((y - x)[:500], label='residual')
        ax[0].legend()
        ax[1].legend()
        # plt.show()

    # Test residual almost 0.0
    np.testing.assert_almost_equal(y - x, np.zeros_like(y), decimal=1)

    # With shifts. We slide the input array by one sample, and check that the
    # artifact is successfully regressed.
    y, idx, mean_total, weights = tspca.tsr(
        signal,
        np.roll(artifact, 1, axis=0),
        shifts=[-1, 0, 1])

    if show:
        f, ax = plt.subplots(3, 1, num='with shifts')
        ax[0].plot(signal[:500], label='signal + noise')
        ax[1].plot(x[:500, 0], 'grey', label='real signal')
        ax[1].plot(y[:500, 0], ':', label='recovered signal')
        ax[2].plot((signal - y)[:500, 0], label='before - after')
        ax[0].legend()
        ax[1].legend()
        ax[2].legend()
        plt.show()

if __name__ == '__main__':
    # import pytest
    # pytest.main([__file__])
    test_tsr()


text:
#!/usr/bin/env python

import keras
import numpy as np
from imgaug import augmenters as iaa
from keras.applications.mobilenet import preprocess_input
from keras.preprocessing.image import image

class CustomDataGenerator(keras.utils.Sequence):
    """
    Takes input of image paths and corresponding labels and generates batches of tensor image data with real-time
    data augmentation. The data will be looped over (in batches).
    """

    def __init__(self, images_paths, labels, batch_size=64, image_dimensions=(224, 224, 3), shuffle=False,
                 augment=False):
        self.images_paths = images_paths  # array of image paths.
        self.labels = labels  # array of labels
        self.batch_size = batch_size  # batch size.
        self.dim = image_dimensions  # image dimensions.
        self.shuffle = shuffle  # shuffle Boolean (default False).
        self.augment = augment  # augment Boolean (default False).
        self.on_epoch_end()

    def __len__(self):
        """Denotes the number of batches per epoch"""
        return int(np.floor(len(self.images_paths) / self.batch_size))

    def on_epoch_end(self):
        """Updates indexes after each epoch"""
        self.indexes = np.arange(len(self.images_paths))
        if self.shuffle:
            np.random.shuffle(self.indexes)

    def __getitem__(self, index):
        """Generate one batch of data"""
        # selects indices of data for next batch
        indexes = self.indexes[index * self.batch_size: (index + 1) * self.batch_size]

        # select data and load images
        labels = np.array([self.labels[k] for k in indexes])
        images = [image.img_to_array(image.load_img(self.images_paths[k], target_size=(224, 224))) for k in indexes]
        # preprocess and augment data
        if self.augment:
            images = self.augmentor(images)
        images = np.array([preprocess_input(img) for img in images])
        return images, labels

    def augmentor(self, images):
        """Apply data augmentation"""
        sometimes = lambda aug: iaa.Sometimes(1, aug)
        seq = iaa.Sequential(
            [
                iaa.Fliplr(0.5),  # horizontally flip 50% of images
                sometimes(iaa.Affine(
                    translate_percent={"x": (-0.1, 0.1), "y": (-0.1, 0.1)},
                    # translate by -20 to +20 percent (per axis)
                    rotate=(-10, 10),  # rotate by -10 to +10 degrees
                ))
            ]
        )
        return seq.augment_images(images)

max_stars_repo_name: djordjebatic/AffectNet

max_issues_repo_name: djordjebatic/AffectNet

max_forks_repo_name: djordjebatic/AffectNet


text:
import base64
from django.http import JsonResponse
from django.shortcuts import render
import io
import matplotlib
matplotlib.use("Agg") # Required to redirect locally
import matplotlib.patches as patches
import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
import tensorflow as tf
from urllib.parse import urlparse
import urllib.request
from .models import URLImage

#   Helper functions for the view functions

def draw_box(box, image_np):
    '''
    This function a NumPy array for a box, and a NumPy array for an image and
    returns a matplotlib pyplot figure and matplotlib pyplot axes indicating a
    box drawn around those coordinates.
    '''

    #   Expand the box by 50%
    box += np.array([-(box[2] - box[0]) / 2,
                     -(box[3] - box[1]) / 2,
                     (box[2] - box[0]) / 2,
                     (box[3] - box[1]) / 2])

    #   Create a large figure to capture all the data
    fig = plt.figure(figsize = (8. * image_np.shape[1] / max(image_np.shape[0], image_np.shape[1]),
                                8. * image_np.shape[0] / max(image_np.shape[0], image_np.shape[1])))

    ax = plt.Axes(fig, [0., 0., 1., 1.])
    fig.add_axes(ax)

    #   Draw blurred boxes around the box
    ax.add_patch(patches.Rectangle((0, 0),
                                   box[1] * image_np.shape[1],
                                   image_np.shape[0],
                                   linewidth = 0,
                                   edgecolor = "none",
                                   facecolor = "w",
                                   alpha=2. / 3.))
    ax.add_patch(patches.Rectangle((box[3] * image_np.shape[1], 0),
                                   image_np.shape[1],
                                   image_np.shape[0],
                                   linewidth = 0,
                                   edgecolor = "none",
                                   facecolor = "w",
                                   alpha = 2. / 3.))
    ax.add_patch(patches.Rectangle((box[1] * image_np.shape[1], 0),
                                   (box[3] - box[1]) * image_np.shape[1],
                                   box[0] * image_np.shape[0],
                                   linewidth = 0,
                                   edgecolor = "none",
                                   facecolor = "w",
                                   alpha = 2. / 3.))
    ax.add_patch(patches.Rectangle((box[1] * image_np.shape[1], box[2] * image_np.shape[0]),
                                   (box[3] - box[1]) * image_np.shape[1],
                                   image_np.shape[0],
                                   linewidth = 0,
                                   edgecolor = "none",
                                   facecolor = "w",
                                   alpha = 2. / 3.))
    return fig, ax

def load_image_into_numpy_array(image):
    '''
    This function takes an image, and returns a Numpy array representation of
    that image.
    '''

    #   Convert 4 channel images into a RGB image for consistency
    if (image.mode == "RGBA") or (image.mode == "CMYK"):
        image = image.convert("RGB")

    (im_width, im_height) = image.size
    return np.array(image.getdata()).reshape(
        (im_height, im_width, 3)).astype(np.uint8)

def open_url_http(url):
    '''
    This function takes a url string and returns an open request from that url.
    '''

    #   Include http prefix if the url omits it
    parse_object = urlparse(url)
    if not parse_object.scheme:
        url = "http://" + url

    #   Request and open the url
    url_request = urllib.request.Request(url)
    return urllib.request.urlopen(url_request)

def go_find_wally(working_graph, url):
    '''
    This function takes a TensorFlow graph and a url string, and returns a Boolean
    indicating if Wally is found and shows an image with his location.
    '''

    tf_config = tf.ConfigProto(intra_op_parallelism_threads = 1,
                               inter_op_parallelism_threads = 1,
                               device_count={"CPU": 1})

    with tf.Session(config = tf_config, graph = working_graph) as sess:

        found = False

        #   Loading variables
        response = open_url_http(url)
        image_np = load_image_into_numpy_array(Image.open(response))
        image_tensor = working_graph.get_tensor_by_name("image_tensor:0")
        boxes = working_graph.get_tensor_by_name("detection_boxes:0")
        scores = working_graph.get_tensor_by_name("detection_scores:0")
        classes = working_graph.get_tensor_by_name("detection_classes:0")
        num_detections = working_graph.get_tensor_by_name("num_detections:0")

        #   Actual detection
        (boxes, scores, classes, num_detections) = sess.run(
            [boxes, scores, classes, num_detections],
            feed_dict = {image_tensor: np.expand_dims(image_np, axis = 0)})

        #   If there is a high enough largest score, draw a box around the prediction
        #   and indicate that Wally is found; otherwise, return an opaque image
        if scores[0][0] > 0.05:
            fig, ax = draw_box(boxes[0][0], image_np)
            found = True
        else:
            fig, ax = draw_box(boxes[0][0] * 0, image_np)

        #   Show image and return if found or not
        ax.imshow(image_np)
        plt.axis("off")
        plt.show()
        return found

#   View functions for the different Webpages on Finding Wally

def error_500(request):
    '''
    This function renders a custom 500 error page telling the user about the
    limitations of the free PythonAnywhere site.
    '''
    data = {}
    return render(request, "error_500.html", data)

def more_information(request):
    '''
    This function renders a request for the More Information page.
    '''

    return render(request, "more_information.html")

def home(request):
    '''
    This function renders a request for the Home page with a variable for all
    urls of the images in the Find Wally carousel.
    '''

    images = None

    #   If there is a POST request to add an image, get the submitted url,
    #   create a URLImage object, and add it to the collection
    if request.POST.get("url_image"):

        #   Get and open the url
        url_image = request.POST.get("url_image")
        response = open_url_http(url_image)

        #   Create the object if the url points to an image
        r = dict(response.getheaders())
        if "image" in r["Content-Type"]:
            URLImage.objects.get_or_create(url = url_image)

    images = URLImage.objects.all()

    #   Render the Home page with a list of all the image urls as a variable
    return render(request, "home.html", {
        "images": images,
    })

def find_wally(request):
    '''
    This function renders a request for a caption and image indicating the
    location of Wally from a search. (This reponse is found within a modal from
    clicking the FIND WALLY button.)
    '''

    #   Location of image and neural network
    #   IMPORTANT: locations to load the relevant image, and the neural network to examine the image for Wally
    image_url = request.GET.get("current_url")
    model_path = "/home/findingwally/mysite/findingwally/trained_model/frozen_inference_graph.pb"

    detection_graph = tf.Graph()
    with detection_graph.as_default():
        od_graph_def = tf.GraphDef()

        #   Read the frozen inference graph
        with tf.gfile.GFile(model_path, "rb") as fid:
            serialized_graph = fid.read()
            od_graph_def.ParseFromString(serialized_graph)
            tf.import_graph_def(od_graph_def, name = "")

        #   Find Wally, show image, and return Boolean
        is_found = go_find_wally(detection_graph, image_url)

    # Redirect the figure to a BytesIO object
    f = io.BytesIO()
    plt.savefig(f, format="png", facecolor = (0.95, 0.95, 0.95), dpi = 400) # High DPI to keep details, around 4K DPI
    plt.clf()
    encoded_img = base64.b64encode(f.getvalue()).decode("utf-8").replace("\n", "")
    f.close()

    # Add the contents of the BytesIO object to the response as a PNG
    if is_found:
        return JsonResponse('

 ' % encoded_img, safe=False)
    else:
        return JsonResponse('

 ' % encoded_img, safe=False)

max_stars_repo_name: pauldnguyen/findingwally

max_issues_repo_name: pauldnguyen/findingwally

max_forks_repo_name: pauldnguyen/findingwally


text:
from __future__ import absolute_import

import ctypes
from .._base import _LIB
from .. import ndarray as _nd

def log_link(input, output, stream=None):
    assert isinstance(input, _nd.NDArray)
    assert isinstance(output, _nd.NDArray)
#    import numpy as np
#    import hetu as ht
#    if input.ctx==ht.gpu(0):
#        np.savetxt("softmax_result2.txt", input.asnumpy())
    _LIB.DLGpuLog(input.handle, output.handle, stream.handle if stream else None)

def log_grad_link(output_grad, input, input_grad, stream=None):
    assert isinstance(output_grad, _nd.NDArray)
    assert isinstance(input, _nd.NDArray)
    assert isinstance(input_grad, _nd.NDArray)

#    import numpy as np
#    import hetu as ht
#    if output_grad.ctx==ht.gpu(0):
#        np.savetxt("log_output_grad2.txt", output_grad.asnumpy())

    _LIB.DLGpuLogGrad(output_grad.handle, input.handle, input_grad.handle, stream.handle if stream else None)

max_stars_repo_name: codecaution/Hetu

max_issues_repo_name: codecaution/Hetu

max_forks_repo_name: codecaution/Hetu


text:
#!/usr/bin/env python
# Licensed under a 3-clause BSD style license - see LICENSE.rst

import os
import sys
import time
import re
import json
import shutil
import logging
import glob

# Matplotlib setup: use Agg backend for command-line (non-interactive) operation
import matplotlib
if __name__ == '__main__':
    matplotlib.use('Agg')

import numpy as np
import sherpa.ui as ui
from mpi4py import MPI
from Chandra.Time import DateTime
import pyyaks.context as pyc

import clogging   # get rid of this or something
import xija

freeze_pars = ('.*',
               '.*__T_e',
               '.*__pi_.*',
               '.*__pf_.*',
               '.*__tau_.*',
               '.*__tau_sc',
               '.*__p_ampl',
               )

#thaw_pars = (#'tcylaft6__pf_045',
             # 'tmzp_my__pf_.*',
             #'tcylfmzm__pf_.*',
             # '.*__p_ampl',
             #'.*__tau_sc',
             #'tcylfmzm__tau_t.*',
             #'.*__tau_ext',
             #'tau_tmzp_cnt',
             #)

class CalcModel(object):
    def __init__(self, model, comm=None, outdir='', pardir=''):
        self.model = model
        self.comm = comm
        if self.comm:
            comm.bcast(dict(cmd='init', model=model.name, tstart=model.tstart, tstop=model.tstop,
                            outdir=model.outdir, pardir=model.pardir),
                       root=MPI.ROOT)

    def __call__(self, parvals, x):
        fit_logger.info('Calculating params:')
        for parname, parval, newparval in zip(self.model.parnames, self.model.parvals, parvals):
            if parval != newparval:
                fit_logger.info('  {0}: {1}'.format(parname, newparval))
        self.model.parvals = parvals

        if self.comm:
            self.comm.bcast(dict(cmd='calc_model', parvals=parvals), root=MPI.ROOT)

        return np.ones_like(x)

class CalcStat(object):
    def __init__(self, model, comm=None):
        self.model = model
        self.comm = comm
        self.cache_fit_stat = {}
        self.min_fit_stat = None

    def __call__(self, data, model, staterror=None, syserror=None, weight=None):
        parvals_key = tuple('%.4e' % x for x in self.model.parvals)
        try:
            fit_stat = self.cache_fit_stat[parvals_key]
            logger.debug('nmass_model: Cache hit %s' % str(parvals_key))
        except KeyError:
            if self.comm:
                msg = dict(cmd='calc_stat', data=data, model=model, staterror=staterror,
                           syserror=syserror, weight=weight)
                comm.bcast(msg, root=MPI.ROOT)
                fit_stat = np.array(0.0, 'd')
                comm.Reduce(None, [fit_stat, MPI.DOUBLE], op=MPI.SUM, root=MPI.ROOT)
                fit_stat = fit_stat.tolist() / comm.Get_size()
            else:
                fit_stat = self.model.calc_stat()

        fit_logger.info('Fit statistic: %.4f' % fit_stat)
        self.cache_fit_stat[parvals_key] = fit_stat

        if self.min_fit_stat is None or fit_stat < self.min_fit_stat:
            self.min_fit_stat = fit_stat
            self.min_parvals = self.model.parvals.copy()

        return fit_stat, np.ones(1)

def fit_model(model,
             comm=None,
             method='simplex',
             config=None,
             nofit=None,
             freeze_pars=freeze_pars,
             thaw_pars=[],
             ):

    dummy_data = np.zeros(1)
    dummy_times = np.arange(1)
    ui.load_arrays(1, dummy_times, dummy_data)

    ui.set_method(method)
    ui.get_method().config.update(config or sherpa_configs.get(method, {}))

    ui.load_user_model(CalcModel(model, comm), 'xijamod')
    ui.add_user_pars('xijamod', model.parnames)
    ui.set_model(1, 'xijamod')

    fit_parnames = set()
    for parname, parval in zip(model.parnames, model.parvals):
        getattr(xijamod, parname).val = parval
        fit_parnames.add(parname)
        if any([re.match(x + '$', parname) for x in freeze_pars]):
            fit_logger.info('Freezing ' + parname)
            ui.freeze(getattr(xijamod, parname))
            fit_parnames.remove(parname)
        if any([re.match(x + '$', parname) for x in thaw_pars]):
            fit_logger.info('Thawing ' + parname)
            ui.thaw(getattr(xijamod, parname))
            fit_parnames.add(parname)
            if 'tau' in parname:
                getattr(xijamod, parname).min = 0.1

    calc_stat = CalcStat(model, comm)
    ui.load_user_stat('xijastat', calc_stat, lambda x: np.ones_like(x))
    ui.set_stat(xijastat)

    if fit_parnames and not nofit:
        ui.fit(1)
    else:
        model.calc()

def make_out_dir():
    out_dir = files['out_dir'].abs
    out_dirs = glob.glob(out_dir + '.*')
    out_dirs_search = re.compile(r'{0}\.(\d+)'.format(out_dir)).search
    out_indexes = [int(out_dirs_search(x).group(1)) for x in out_dirs if out_dirs_search(x)]
    if out_indexes:
        n_outs = max(out_indexes)
    else:
        n_outs = 0
    if os.path.exists(out_dir):
        os.rename(out_dir, out_dir + '.%d' % (n_outs + 1))
    os.makedirs(out_dir)

def get_options():
    from optparse import OptionParser
    parser = OptionParser()
    parser.add_option("--days",
                      type='float',
                      default=90,
                      help="Number of days in fit interval (default=90")
    parser.add_option("--stop",
                      help="Stop time of fit interval (default=NOW - 7 days)")
    parser.add_option("--method",
                      default='simplex',
                      help="Fit method (default=simple)")
    parser.add_option("--model",
                      default='minusz',
                      help="Model to predict (default=minusz)")
    parser.add_option("--thaw-pars",
                      default='solarheat__tephin__P_.* solarheat__tephin__dP_.*',
                      help="List of parameters (space-separated) to thaw (default=none)")
    parser.add_option("--nproc",
                      default=0,
                      type='int',
                      help="Number of processors (default=1)")
    parser.add_option("--ftol",
                      default=1e-4,
                      type='float',
                      help="ftol convergence parameter (default=1e-4")
    parser.add_option("--outdir",
                      default='fit',
                      help="Output directory within <msid> root (default=fit)")
    parser.add_option("--pardir",
                      help="Directory containing model params file (default=<model>)")
    parser.add_option("--quiet",
                      default=False,
                      action='store_true',
                      help="Suppress screen output")
    parser.add_option("--nofit",
                      action="store_true",
                      help="Do not fit (default=False)")
    (opt, args) = parser.parse_args()
    return (opt, args)

# Default configurations for fit methods
sherpa_configs = dict(
    simplex = dict(ftol=1e-3,
                   finalsimplex=0,   # converge based only on length of simplex
                   maxfev=1000),
    )

opt, args = get_options()

src = pyc.CONTEXT['src'] if 'src' in pyc.CONTEXT else pyc.ContextDict('src')
files = (pyc.CONTEXT['file'] if 'file' in pyc.CONTEXT else
         pyc.ContextDict('files', basedir=os.getcwd()))
files.update(xija.files)

fit_logger = clogging.config_logger('fit', level=clogging.INFO)
sherpa_logger = logging.getLogger("sherpa")
loggers = (fit_logger, sherpa_logger)
if opt.quiet:
    for logger in loggers:
        for h in logger.handlers:
            logger.removeHandler(h)

# Use supplied stop time or NOW - 7 days (truncated to nearest day)
stop = opt.stop or DateTime(DateTime().secs - 7 * 86400).date[:8]
start = DateTime(DateTime(stop).secs - opt.days * 86400).date[:8]

src['model'] = opt.model
src['outdir'] = opt.outdir
src['pardir'] = opt.pardir or opt.model

model_spec = json.load(open(files['model_spec.json'].abs, 'r'))
model = xija.ThermalModel(opt.model, start, stop, model_spec=model_spec)

### UGHH, need to clean this up, probably have model elements accumulate dynamically
### instead of a make step.
if not opt.nproc:
    model.make()

model.outdir = src['outdir'].val
model.pardir = src['pardir'].val

thaw_pars = opt.thaw_pars.split()

make_out_dir()

hdlr = logging.FileHandler(files['fit_log'].abs, 'w')
hdlr.setLevel(logging.INFO)
hdlr.setFormatter(logging.Formatter('%(message)s'))
for logger in loggers:
    logger.addHandler(hdlr)

fit_logger.info('Running: ' + ' '.join(sys.argv))
fit_logger.info('Start: {0}'.format(start))
fit_logger.info('Stop: {0}'.format(stop))
fit_logger.info('Options: {0}'.format(opt))
fit_logger.info('')
fit_logger.info('Start time: {0}'.format(time.ctime()))

if opt.nproc:
    fit_logger.info('fit_nmass: Spawning processes')
    comm = MPI.COMM_SELF.Spawn(sys.executable, args=['worker.py'], maxprocs=opt.nproc)
    fit_logger.info('fit_nmass: Finished spawning processes')
    fit_logger.info('fit_nmass: comm.Get_size() = {0}'.format(comm.Get_size()))
else:
    comm = None

try:
    fit_model(model, comm=comm, method=opt.method, nofit=opt.nofit, thaw_pars=thaw_pars)
    model.write(files['fit_spec.json'].abs)

    kwargs = {'savefig': True}
    if comm:
        comm.bcast(dict(cmd='model', func='plot_fit_resids', kwargs=kwargs), root=MPI.ROOT)
    else:
        model.plot_fit_resids(**kwargs)

    if comm:
        pass
        #comm.bcast(dict(cmd='model', func='plot_hist', kwargs=kwargs), root=MPI.ROOT)
    else:
        if 0:
            model.plot_hist(select='all', **kwargs)
            model.plot_hist(select='hot', **kwargs)
            model.write_fit_data()
            model.write_html_report()

finally:
    if comm:
        comm.bcast({'cmd': 'stop'}, root=MPI.ROOT)
        semaphore = np.array(0.0, 'd')
        comm.Reduce(None, [semaphore, MPI.DOUBLE], op=MPI.SUM, root=MPI.ROOT)

        time.sleep(1)
        comm.Disconnect()

    fit_logger.info('Stop time: {0}'.format(time.ctime()))

max_stars_repo_name: jzuhone/xija

max_issues_repo_name: jzuhone/xija

max_forks_repo_name: jzuhone/xija


text:
import time
import random
import os
import logging
from shutil import rmtree
import numpy as np
import torch
import schnetpack as spk

__all__ = ["ScriptError", "setup_run", "get_environment_provider"]

class ScriptError(Exception):
    pass

def setup_run(args):
    if args.mode == "train":

        # build modeldir
        if args.overwrite and os.path.exists(args.modelpath):
            logging.info("existing model will be overwritten...")
            rmtree(args.modelpath)
        if not os.path.exists(args.modelpath):
            os.makedirs(args.modelpath)

        _set_random_seed(args.seed)
        train_args = args
    else:
        # check if modelpath is valid
        if not os.path.exists(args.modelpath):
            raise ScriptError(
                "The selected modeldir does not exist " "at {}!".format(args.modelpath)
            )

        # load training arguments
        train_args = args

    # apply alias definitions
    train_args = _apply_aliases(train_args)
    return train_args

def _set_random_seed(seed: int):
    """
    This function sets the random seed (if given) or creates one for torch and numpy random state initialization

    Parameters
    ----------
    seed : int or None, default=None
        if seed not present, it is generated based on time
    """
    if seed is None:
        seed = int(time.time() * 1000.0)
        # Reshuffle current time to get more different seeds within shorter time intervals
        # Taken from https://stackoverflow.com/questions/27276135/python-random-system-time-seed
        # & Gets overlapping bits, << and >> are binary right and left shifts
        seed = (
            ((seed & 0xFF000000) >> 24)
            + ((seed & 0x00FF0000) >> 8)
            + ((seed & 0x0000FF00) << 8)
            + ((seed & 0x000000FF) << 24)
        )
    # random
    random.seed(seed)
    # Numpy
    np.random.seed(seed)
    # Pytorch
    torch.manual_seed(seed)
    torch.cuda.manual_seed_all(seed)
    torch.backends.cudnn.deterministic = True
    logging.info("Random state initialized with seed {:<10d}".format(seed))

def _apply_aliases(args):
    # force alias for custom dataset
    if args.dataset == "custom":
        if args.force is not None:
            if args.derivative is not None:
                raise ScriptError(
                    "Force and derivative define the same property. Please don`t use "
                    "both."
                )
            args.derivative = args.force
            args.negative_dr = True

            # add rho value if selected
            if "force" in args.rho.keys():
                args.rho["derivative"] = args.rho.pop("force")

    return args

def get_environment_provider(args, device):
    if args.environment_provider == "simple":
        return spk.environment.SimpleEnvironmentProvider()
    elif args.environment_provider == "ase":
        return spk.environment.AseEnvironmentProvider(cutoff=args.cutoff)
    elif args.environment_provider == "torch":
        return spk.environment.TorchEnvironmentProvider(
            cutoff=args.cutoff, device="cpu"
        )
    else:
        raise NotImplementedError


text:
import femm
import numpy as np
from os.path import basename, splitext

from numpy import zeros, pi, roll, mean, max as np_max, min as np_min
from ....Functions.FEMM.update_FEMM_simulation import update_FEMM_simulation
from ....Functions.FEMM.comp_FEMM_torque import comp_FEMM_torque
from ....Functions.FEMM.comp_FEMM_Phi_wind import comp_FEMM_Phi_wind
from ....Classes.MeshSolution import MeshSolution
from ....Classes.Mesh import Mesh
from ....Classes.Solution import Solution
from SciDataTool import DataLinspace, DataTime
from os.path import join

def solve_FEMM(self, output, sym, FEMM_dict):

    # Loading parameters for readibilitys
    angle = output.mag.angle

    L1 = output.simu.machine.stator.comp_length()
    Nt_tot = output.mag.Nt_tot  # Number of time step
    Na_tot = output.mag.Na_tot  # Number of angular step
    save_path = self.get_path_save(output)

    if (
        hasattr(output.simu.machine.stator, "winding")
        and output.simu.machine.stator.winding is not None
    ):
        qs = output.simu.machine.stator.winding.qs  # Winding phase number
        Npcpp = output.simu.machine.stator.winding.Npcpp
        Phi_wind_stator = zeros((Nt_tot, qs))
    else:
        Phi_wind_stator = None

    # Create the mesh
    femm.mi_createmesh()

    # Initialize results matrix
    Br = zeros((Nt_tot, Na_tot))
    Bt = zeros((Nt_tot, Na_tot))
    Tem = zeros((Nt_tot, 1))

    lam_int = output.simu.machine.get_lamination(True)
    lam_ext = output.simu.machine.get_lamination(False)
    Rgap_mec_int = lam_int.comp_radius_mec()
    Rgap_mec_ext = lam_ext.comp_radius_mec()

    if self.is_get_mesh or self.is_save_FEA:
        meshFEMM = [Mesh() for ii in range(Nt_tot)]
        solutionFEMM = [Solution() for ii in range(Nt_tot)]
    else:
        meshFEMM = [Mesh()]
        solutionFEMM = [Solution()]

    # Compute the data for each time step
    for ii in range(Nt_tot):
        # Update rotor position and currents
        update_FEMM_simulation(
            output=output,
            materials=FEMM_dict["materials"],
            circuits=FEMM_dict["circuits"],
            is_mmfs=self.is_mmfs,
            is_mmfr=self.is_mmfr,
            j_t0=ii,
            is_sliding_band=self.is_sliding_band,
        )
        # try "previous solution" for speed up of FEMM calculation
        if self.is_sliding_band:
            try:
                base = basename(self.get_path_save_fem(output))
                ans_file = splitext(base)[0] + ".ans"
                femm.mi_setprevious(ans_file, 0)
            except:
                pass

        # Run the computation
        femm.mi_analyze()
        femm.mi_loadsolution()

        # Get the flux result
        if self.is_sliding_band:
            for jj in range(Na_tot):
                Br[ii, jj], Bt[ii, jj] = femm.mo_getgapb("bc_ag2", angle[jj] * 180 / pi)
        else:
            Rag = (Rgap_mec_ext + Rgap_mec_int) / 2
            for jj in range(Na_tot):
                B = femm.mo_getb(Rag * np.cos(angle[jj]), Rag * np.sin(angle[jj]))
                Br[ii, jj] = B[0] * np.cos(angle[jj]) + B[1] * np.sin(angle[jj])
                Bt[ii, jj] = -B[0] * np.sin(angle[jj]) + B[1] * np.cos(angle[jj])

        # Compute the torque
        Tem[ii] = comp_FEMM_torque(FEMM_dict, sym=sym)

        if (
            hasattr(output.simu.machine.stator, "winding")
            and output.simu.machine.stator.winding is not None
        ):
            # Phi_wind computation
            Phi_wind_stator[ii, :] = comp_FEMM_Phi_wind(
                qs, Npcpp, is_stator=True, Lfemm=FEMM_dict["Lfemm"], L1=L1, sym=sym
            )

        # Load mesh data & solution
        if self.is_get_mesh or self.is_save_FEA:
            meshFEMM[ii], solutionFEMM[ii] = self.get_meshsolution(
                self.is_get_mesh, self.is_save_FEA, save_path, ii
            )

    # Shift to take into account stator position
    roll_id = int(self.angle_stator * Na_tot / (2 * pi))
    Br = roll(Br, roll_id, axis=1)
    Bt = roll(Bt, roll_id, axis=1)

    # Store the results
    Time = DataLinspace(
        name="time",
        unit="s",
        symmetries={},
        initial=output.mag.time[0],
        final=output.mag.time[-1],
        number=Nt_tot,
    )
    Angle = DataLinspace(
        name="angle",
        unit="rad",
        symmetries={},
        initial=angle[0],
        final=angle[-1],
        number=Na_tot,
    )
    output.mag.Br = DataTime(
        name="Airgap radial flux density",
        unit="T",
        symbol="B_r",
        axes=[Time, Angle],
        values=Br,
    )
    output.mag.Bt = DataTime(
        name="Airgap tangential flux density",
        unit="T",
        symbol="B_r",
        axes=[Time, Angle],
        values=Bt,
    )
    output.mag.Tem = Tem
    output.mag.Tem_av = mean(Tem)
    if output.mag.Tem_av != 0:
        output.mag.Tem_rip = abs((np_max(Tem) - np_min(Tem)) / output.mag.Tem_av)
    output.mag.Phi_wind_stator = Phi_wind_stator
    output.mag.FEMM_dict = FEMM_dict

    if self.is_get_mesh:
        cond = (not self.is_sliding_band) or (Nt_tot == 1)
        output.mag.meshsolution = MeshSolution(
            name="FEMM_magnetic_mesh",
            mesh=meshFEMM,
            solution=solutionFEMM,
            is_same_mesh=cond,
        )

    if self.is_save_FEA:
        save_path_fea = join(save_path, "MeshSolutionFEMM.json")
        output.mag.meshsolution.save(save_path_fea)

    if (
        hasattr(output.simu.machine.stator, "winding")
        and output.simu.machine.stator.winding is not None
    ):
        # Electromotive forces computation (update output)
        self.comp_emf()
    else:
        output.mag.emf = None

max_stars_repo_name: magnetron/pyleecan

max_issues_repo_name: magnetron/pyleecan

max_forks_repo_name: magnetron/pyleecan