text:
# emacs: -*- mode: python-mode; 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 NiBabel package for the
#   copyright and license terms.
#
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
''' Very simple spatial image class

The image class maintains the association between a 3D (or greater)
array, and an affine transform that maps voxel coordinates to some real
world space.  It also has a ``header`` - some standard set of meta-data
that is specific to the image format - and ``extra`` - a dictionary
container for any other metadata.

It has attributes:

   * extra

methods:

   * .get_data()
   * .get_affine()
   * .get_header()
   * .set_shape(shape)
   * .to_filename(fname) - writes data to filename(s) derived from
     ``fname``, where the derivation may differ between formats.
   * to_file_map() - save image to files with which the image is already
     associated.
   * .get_shape() (Deprecated)

properties:

   * shape

classmethods:

   * from_filename(fname) - make instance by loading from filename
   * instance_to_filename(img, fname) - save ``img`` instance to
     filename ``fname``.

There are several ways of writing data.
=======================================

There is the usual way, which is the default::

    img.to_filename(fname)

and that is, to take the data encapsulated by the image and cast it to
the datatype the header expects, setting any available header scaling
into the header to help the data match.

You can load the data into an image from file with::

   img.from_filename(fname)

The image stores its associated files in its ``files`` attribute.  In
order to just save an image, for which you know there is an associated
filename, or other storage, you can do::

   img.to_file_map()

You can get the data out again with of::

    img.get_data()

Less commonly, for some image types that support it, you might want to
fetch out the unscaled array via the header::

    unscaled_data = img.get_unscaled_data()

Analyze-type images (including nifti) support this, but others may not
(MINC, for example).

Sometimes you might to avoid any loss of precision by making the
data type the same as the input::

    hdr = img.get_header()
    hdr.set_data_dtype(data.dtype)
    img.to_filename(fname)

Files interface
===============

The image has an attribute ``file_map``.  This is a mapping, that has keys
corresponding to the file types that an image needs for storage.  For
example, the Analyze data format needs an ``image`` and a ``header``
file type for storage:

   >>> import nibabel as nib
   >>> data = np.arange(24, dtype='f4').reshape((2,3,4))
   >>> img = nib.AnalyzeImage(data, np.eye(4))
   >>> sorted(img.file_map)
   ['header', 'image']

The values of ``file_map`` are not in fact files but objects with
attributes ``filename``, ``fileobj`` and ``pos``.

The reason for this interface, is that the contents of files has to
contain enough information so that an existing image instance can save
itself back to the files pointed to in ``file_map``.  When a file holder
holds active file-like objects, then these may be affected by the
initial file read; in this case, the contains file-like objects need to
carry the position at which a write (with ``to_files``) should place the
data.  The ``file_map`` contents should therefore be such, that this will
work:

   >>> # write an image to files
   >>> from io import BytesIO
   >>> file_map = nib.AnalyzeImage.make_file_map()
   >>> file_map['image'].fileobj = BytesIO()
   >>> file_map['header'].fileobj = BytesIO()
   >>> img = nib.AnalyzeImage(data, np.eye(4))
   >>> img.file_map = file_map
   >>> img.to_file_map()
   >>> # read it back again from the written files
   >>> img2 = nib.AnalyzeImage.from_file_map(file_map)
   >>> np.all(img2.get_data() == data)
   True
   >>> # write, read it again
   >>> img2.to_file_map()
   >>> img3 = nib.AnalyzeImage.from_file_map(file_map)
   >>> np.all(img3.get_data() == data)
   True

'''

try:
    basestring
except NameError:  # python 3
    basestring = str

import warnings

import numpy as np

from .filename_parser import types_filenames, TypesFilenamesError
from .fileholders import FileHolder
from .volumeutils import shape_zoom_affine

class HeaderDataError(Exception):
    ''' Class to indicate error in getting or setting header data '''
    pass

class HeaderTypeError(Exception):
    ''' Class to indicate error in parameters into header functions '''
    pass

class Header(object):
    ''' Template class to implement header protocol '''
    default_x_flip = True

    def __init__(self,
                 data_dtype=np.float32,
                 shape=(0,),
                 zooms=None):
        self.set_data_dtype(data_dtype)
        self._zooms = ()
        self.set_data_shape(shape)
        if not zooms is None:
            self.set_zooms(zooms)

    @classmethod
    def from_header(klass, header=None):
        if header is None:
            return klass()
        # I can't do isinstance here because it is not necessarily true
        # that a subclass has exactly the same interface as it's parent
        # - for example Nifti1Images inherit from Analyze, but have
        # different field names
        if type(header) == klass:
            return header.copy()
        return klass(header.get_data_dtype(),
                     header.get_data_shape(),
                     header.get_zooms())

    @classmethod
    def from_fileobj(klass, fileobj):
        raise NotImplementedError

    def write_to(self, fileobj):
        raise NotImplementedError

    def __eq__(self, other):
        return ((self.get_data_dtype(),
                 self.get_data_shape(),
                 self.get_zooms()) ==
                (other.get_data_dtype(),
                 other.get_data_shape(),
                 other.get_zooms()))

    def __ne__(self, other):
        return not self == other

    def copy(self):
        ''' Copy object to independent representation

        The copy should not be affected by any changes to the original
        object.
        '''
        return self.__class__(self._dtype, self._shape, self._zooms)

    def get_data_dtype(self):
        return self._dtype

    def set_data_dtype(self, dtype):
        self._dtype = np.dtype(dtype)

    def get_data_shape(self):
        return self._shape

    def set_data_shape(self, shape):
        ndim = len(shape)
        if ndim == 0:
            self._shape = (0,)
            self._zooms = (1.0,)
            return
        self._shape = tuple([int(s) for s in shape])
        # set any unset zooms to 1.0
        nzs = min(len(self._zooms), ndim)
        self._zooms = self._zooms[:nzs] + (1.0,) * (ndim-nzs)

    def get_zooms(self):
        return self._zooms

    def set_zooms(self, zooms):
        zooms = tuple([float(z) for z in zooms])
        shape = self.get_data_shape()
        ndim = len(shape)
        if len(zooms) != ndim:
            raise HeaderDataError('Expecting %d zoom values for ndim %d'
                                  % (ndim, ndim))
        if len([z for z in zooms if z < 0]):
            raise HeaderDataError('zooms must be positive')
        self._zooms = zooms

    def get_base_affine(self):
        shape = self.get_data_shape()
        zooms = self.get_zooms()
        return shape_zoom_affine(shape, zooms,
                                 self.default_x_flip)

    get_best_affine = get_base_affine

    def data_to_fileobj(self, data, fileobj):
        ''' Write image data to file in fortran order '''
        dtype = self.get_data_dtype()
        fileobj.write(data.astype(dtype).tostring(order='F'))

    def data_from_fileobj(self, fileobj):
        ''' Read data in fortran order '''
        dtype = self.get_data_dtype()
        shape = self.get_data_shape()
        data_size = int(np.prod(shape) * dtype.itemsize)
        data_bytes = fileobj.read(data_size)
        return np.ndarray(shape, dtype, data_bytes, order='F')

class ImageDataError(Exception):
    pass

class ImageFileError(Exception):
    pass

class SpatialImage(object):
    header_class = Header
    files_types = (('image', None),)
    _compressed_exts = ()

    ''' Template class for images '''
    def __init__(self, dataobj, affine, header=None,
                 extra=None, file_map=None):
        ''' Initialize image

        The image is a combination of (array, affine matrix, header), with
        optional metadata in `extra`, and filename / file-like objects contained
        in the `file_map` mapping.

        Parameters
        ----------
        dataobj : object
           Object containg image data.  It should be some object that retuns an
           array from ``np.asanyarray``.  It should have a ``shape`` attribute
           or property
        affine : None or (4,4) array-like
           homogenous affine giving relationship between voxel coordinates and
           world coordinates.  Affine can also be None.  In this case,
           ``obj.get_affine()`` also returns None, and the affine as written to
           disk will depend on the file format.
        header : None or mapping or header instance, optional
           metadata for this image format
        extra : None or mapping, optional
           metadata to associate with image that cannot be stored in the
           metadata of this image type
        file_map : mapping, optional
           mapping giving file information for this image format
        '''
        self._dataobj = dataobj
        if not affine is None:
            # Check that affine is array-like 4,4.  Maybe this is too strict at
            # this abstract level, but so far I think all image formats we know
            # do need 4,4.
            # Copy affine to isolate from environment.  Specify float type to
            # avoid surprising integer rounding when setting values into affine
            affine = np.array(affine, dtype=np.float64, copy=True)
            if not affine.shape == (4,4):
                raise ValueError('Affine should be shape 4,4')
        self._affine = affine
        if extra is None:
            extra = {}
        self.extra = extra
        self._header = self.header_class.from_header(header)
        # if header not specified, get data type from input array
        if header is None:
            if hasattr(dataobj, 'dtype'):
                self._header.set_data_dtype(dataobj.dtype)
        # make header correspond with image and affine
        self.update_header()
        if file_map is None:
            file_map = self.__class__.make_file_map()
        self.file_map = file_map
        self._load_cache = None
        self._data_cache = None

    @property
    def _data(self):
        warnings.warn('Please use ``dataobj`` instead of ``_data``; '
                      'We will remove this wrapper for ``_data`` soon',
                      FutureWarning,
                      stacklevel=2)
        return self._dataobj

    @property
    def dataobj(self):
        return self._dataobj

    @property
    def affine(self):
        return self._affine

    @property
    def header(self):
        return self._header

    def update_header(self):
        ''' Harmonize header with image data and affine

        >>> data = np.zeros((2,3,4))
        >>> affine = np.diag([1.0,2.0,3.0,1.0])
        >>> img = SpatialImage(data, affine)
        >>> hdr = img.get_header()
        >>> img.shape == (2, 3, 4)
        True
        >>> img.update_header()
        >>> hdr.get_data_shape() == (2, 3, 4)
        True
        >>> hdr.get_zooms()
        (1.0, 2.0, 3.0)
        '''
        hdr = self._header
        shape = self._dataobj.shape
        # We need to update the header if the data shape has changed.  It's a
        # bit difficult to change the data shape using the standard API, but
        # maybe it happened
        if hdr.get_data_shape() != shape:
            hdr.set_data_shape(shape)
        # If the affine is not None, and it is different from the main affine in
        # the header, update the heaader
        if self._affine is None:
            return
        if np.allclose(self._affine, hdr.get_best_affine()):
            return
        self._affine2header()

    def _affine2header(self):
        """ Unconditionally set affine into the header """
        RZS = self._affine[:3, :3]
        vox = np.sqrt(np.sum(RZS * RZS, axis=0))
        hdr = self._header
        zooms = list(hdr.get_zooms())
        n_to_set = min(len(zooms), 3)
        zooms[:n_to_set] = vox[:n_to_set]
        hdr.set_zooms(zooms)

    def __str__(self):
        shape = self.shape
        affine = self.get_affine()
        return '\n'.join((
                str(self.__class__),
                'data shape %s' % (shape,),
                'affine: ',
                '%s' % affine,
                'metadata:',
                '%s' % self._header))

    def get_data(self):
        """ Return image data from image with any necessary scalng applied

        If the image data is a array proxy (data not yet read from disk) then
        read the data, and store in an internal cache.  Future calls to
        ``get_data`` will return the cached copy.

        Returns
        -------
        data : array
            array of image data
        """
        if self._data_cache is None:
            self._data_cache = np.asanyarray(self._dataobj)
        return self._data_cache

    def uncache(self):
        """ Delete any cached read of data from proxied data

        Remember there are two types of images:

        * *array images* where the data ``img.dataobj`` is an array
        * *proxy images* where the data ``img.dataobj`` is a proxy object

        If you call ``img.get_data()`` on a proxy image, the result of reading
        from the proxy gets cached inside the image object, and this cache is
        what gets returned from the next call to ``img.get_data()``.  If you
        modify the returned data, as in::

            data = img.get_data()
            data[:] = 42

        then the next call to ``img.get_data()`` returns the modified array,
        whether the image is an array image or a proxy image::

            assert np.all(img.get_data() == 42)

        When you uncache an array image, this has no effect on the return of
        ``img.get_data()``, but when you uncache a proxy image, the result of
        ``img.get_data()`` returns to its original value.
        """
        self._data_cache = None

    @property
    def shape(self):
        return self._dataobj.shape

    def get_shape(self):
        """ Return shape for image

        This function deprecated; please use the ``shape`` property instead
        """
        warnings.warn('Please use the shape property instead of get_shape',
                      DeprecationWarning,
                      stacklevel=2)
        return self.shape

    def get_data_dtype(self):
        return self._header.get_data_dtype()

    def set_data_dtype(self, dtype):
        self._header.set_data_dtype(dtype)

    def get_affine(self):
        """ Get affine from image

        Please use the `affine` property instead of `get_affine`; we will
        deprecate this method in future versions of nibabel.
        """
        return self.affine

    def get_header(self):
        """ Get header from image

        Please use the `header` property instead of `get_header`; we will
        deprecate this method in future versions of nibabel.
        """
        return self.header

    def get_filename(self):
        ''' Fetch the image filename

        Parameters
        ----------
        None

        Returns
        -------
        fname : None or str
           Returns None if there is no filename, or a filename string.
           If an image may have several filenames assoctiated with it
           (e.g Analyze ``.img, .hdr`` pair) then we return the more
           characteristic filename (the ``.img`` filename in the case of
           Analyze')
        '''
        # which filename is returned depends on the ordering of the
        # 'files_types' class attribute - we return the name
        # corresponding to the first in that tuple
        characteristic_type = self.files_types[0][0]
        return self.file_map[characteristic_type].filename

    def set_filename(self, filename):
        ''' Sets the files in the object from a given filename

        The different image formats may check whether the filename has
        an extension characteristic of the format, and raise an error if
        not.

        Parameters
        ----------
        filename : str
           If the image format only has one file associated with it,
           this will be the only filename set into the image
           ``.file_map`` attribute. Otherwise, the image instance will
           try and guess the other filenames from this given filename.
        '''
        self.file_map = self.__class__.filespec_to_file_map(filename)

    @classmethod
    def from_filename(klass, filename):
        file_map = klass.filespec_to_file_map(filename)
        return klass.from_file_map(file_map)

    @classmethod
    def from_filespec(klass, filespec):
        warnings.warn('``from_filespec`` class method is deprecated\n'
                      'Please use the ``from_filename`` class method '
                      'instead',
                      DeprecationWarning, stacklevel=2)
        klass.from_filename(filespec)

    @classmethod
    def from_file_map(klass, file_map):
        raise NotImplementedError

    @classmethod
    def from_files(klass, file_map):
        warnings.warn('``from_files`` class method is deprecated\n'
                      'Please use the ``from_file_map`` class method '
                      'instead',
                      DeprecationWarning, stacklevel=2)
        return klass.from_file_map(file_map)

    @classmethod
    def filespec_to_file_map(klass, filespec):
        try:
            filenames = types_filenames(filespec,
                                        klass.files_types,
                                        trailing_suffixes=klass._compressed_exts)
        except TypesFilenamesError:
            raise ImageFileError('Filespec "%s" does not look right for '
                             'class %s ' % (filespec, klass))
        file_map = {}
        for key, fname in filenames.items():
            file_map[key] = FileHolder(filename=fname)
        return file_map

    @classmethod
    def filespec_to_files(klass, filespec):
        warnings.warn('``filespec_to_files`` class method is deprecated\n'
                      'Please use the ``filespec_to_file_map`` class method '
                      'instead',
                      DeprecationWarning, stacklevel=2)
        return klass.filespec_to_file_map(filespec)

    def to_filename(self, filename):
        ''' Write image to files implied by filename string

        Parameters
        ----------
        filename : str
           filename to which to save image.  We will parse `filename`
           with ``filespec_to_file_map`` to work out names for image,
           header etc.

        Returns
        -------
        None
        '''
        self.file_map = self.filespec_to_file_map(filename)
        self.to_file_map()

    def to_filespec(self, filename):
        warnings.warn('``to_filespec`` is deprecated, please '
                      'use ``to_filename`` instead',
                      DeprecationWarning, stacklevel=2)
        self.to_filename(filename)

    def to_file_map(self, file_map=None):
        raise NotImplementedError

    def to_files(self, file_map=None):
        warnings.warn('``to_files`` method is deprecated\n'
                      'Please use the ``to_file_map`` method '
                      'instead',
                      DeprecationWarning, stacklevel=2)
        self.to_file_map(file_map)

    @classmethod
    def make_file_map(klass, mapping=None):
        ''' Class method to make files holder for this image type

        Parameters
        ----------
        mapping : None or mapping, optional
           mapping with keys corresponding to image file types (such as
           'image', 'header' etc, depending on image class) and values
           that are filenames or file-like.  Default is None

        Returns
        -------
        file_map : dict
           dict with string keys given by first entry in tuples in
           sequence klass.files_types, and values of type FileHolder,
           where FileHolder objects have default values, other than
           those given by `mapping`
        '''
        if mapping is None:
            mapping = {}
        file_map = {}
        for key, ext in klass.files_types:
            file_map[key] = FileHolder()
            mapval = mapping.get(key, None)
            if isinstance(mapval, basestring):
                file_map[key].filename = mapval
            elif hasattr(mapval, 'tell'):
                file_map[key].fileobj = mapval
        return file_map

    @classmethod
    def load(klass, filename):
        return klass.from_filename(filename)

    @classmethod
    def instance_to_filename(klass, img, filename):
        ''' Save `img` in our own format, to name implied by `filename`

        This is a class method

        Parameters
        ----------
        img : ``spatialimage`` instance
           In fact, an object with the API of ``spatialimage`` -
           specifically ``get_data``, ``get_affine``, ``get_header`` and
           ``extra``.
        filename : str
           Filename, implying name to which to save image.
        '''
        img = klass.from_image(img)
        img.to_filename(filename)

    @classmethod
    def from_image(klass, img):
        ''' Class method to create new instance of own class from `img`

        Parameters
        ----------
        img : ``spatialimage`` instance
           In fact, an object with the API of ``spatialimage`` -
           specifically ``get_data``, ``get_affine``, ``get_header`` and
           ``extra``.

        Returns
        -------
        cimg : ``spatialimage`` instance
           Image, of our own class
        '''
        return klass(img._dataobj,
                     img._affine,
                     klass.header_class.from_header(img._header),
                     extra=img.extra.copy())

max_stars_repo_name: tobon/nibabel

max_issues_repo_name: tobon/nibabel

max_forks_repo_name: tobon/nibabel


text:
from __future__ import division
from numpy import linalg as LA

import numpy as np

np.set_printoptions(threshold=np.inf)
import sys
import librosa
import matplotlib.pyplot as plt

def buffer(signal, L, M):
    if M >= L:
        print ('Error: Overlapping windows cannot be larger than frame length!')
        sys.exit()

    len_signal = len(signal)

    print ('The signal length is %s: ' % (len_signal))

    K = np.ceil(len_signal / L).astype('int')  # num_frames

    print ('The number of frames \'K\' is %s: ' % (K))
    print ('The length of each frame \'L\' is %s: ' % (L))

    X_tmp = []
    k = 1
    while (True):
        start_ind = ((k - 1) * (L - M) + 1) - 1
        end_ind = ((k * L) - (k - 1) * M)
        if start_ind == len_signal:
            break
        elif (end_ind > len_signal):
            # print ('k=%s, [%s, %s] ' % (k, start_ind, end_ind - 1))
            val_in = len_signal - start_ind
            tmp_seg = np.zeros(L)
            tmp_seg[:val_in] = signal[start_ind:]
            X_tmp.append(tmp_seg)
            break
        else:
            # print ('k=%s, [%s, %s] ' % (k, start_ind, end_ind - 1))
            X_tmp.append(signal[start_ind:end_ind])
        k += 1

    return X_tmp

def unbuffer(X, hop):
    N, L = X.shape

    T = N + L * hop
    K = np.arange(0, N)
    x = np.zeros(T)
    H = np.hanning(N)
    for k in xrange(0, L):
        x[K] = x[K] + np.multiply(H, X[:, k])
        K = K + hop

    return x

class GAD():
    def __init__(self, X, params):

        self.X = X
        self.D = []
        self.params = params
        self.n_iter = self.params['rule_1']['n_iter']  # num_iterations
        self.verbose = self.params['verbose']

        self.K = self.X.shape[0]  # sample length
        self.L = self.X.shape[1]  # maximum atoms to be learned

        self.I = np.arange(0, self.L)
        self.set_ind = []

    def findResidualColumn(self):

        # Find residual column of R^l with lowest l1- to l2-norm ration
        tmp = []

        # COMPACT WAY TO DO IT
        # tmp2 = np.sum(np.abs(self.R),axis=0)/np.sqrt(np.sum(np.power(np.abs(self.R),2),axis=0))
        for k in self.I:
            r_k = self.R[:, k]
            tmp.append(LA.norm(r_k, 1) / LA.norm(r_k, 2))

        ind_k_min = np.nanargmin(tmp)  # nanargmin, nanmin
        k_min = tmp[ind_k_min]
        r_k_min = self.R[:, self.I[ind_k_min]]

        # Set the l-th atom to equal to normalized r_k
        psi = r_k_min / LA.norm(r_k_min, 2)

        # Add to the dictionary D and its index and shrinking set I
        self.D.append(psi)
        self.set_ind.append(self.I[ind_k_min])

        # COMPACT WAY TO DO IT
        # self.R = self.R - np.dot(np.outer(psi, psi), self.R)
        # self.R = np.delete(self.R, (ind_k_min), axis=1)

        # Compute the new residual for all columns k
        for k in self.I:
            r_k = self.R[:, k]
            alpha = np.dot(r_k, psi)
            self.R[:, k] = r_k - np.dot(psi, alpha)

        self.I = np.delete(self.I, ind_k_min)

    def iterative_GAD(self):

        # X columns w/ unit L2-norm
        # for k in xrange(0,self.X.shape[1]):
        #     self.X[:,k] = np.divide(self.X[:,k],LA.norm(self.X[:,k],2))

        if self.n_iter > self.L:
            print ('Cannot be learned more than %d atom!' % (self.L))
            sys.exit()

        # Initializating the residual matrix 'R' by using 'X'
        self.R = self.X.copy()

        print self.I.shape
        for l in xrange(0, self.n_iter):

            if self.verbose == True:
                print 'GAD iteration: ', l + 1

            self.findResidualColumn()

        self.D = np.vstack(self.D).T

        return self.D, self.set_ind

if __name__ == '__main__':
    L = 512  # frame length
    M = 500  # overlapping windows

    params = {

        'rule_1': {
            'n_iter': 10  # n_iter
        },

        'rule_2': {
            'error': 10 ** -7
        },

        'verbose': True
    }

    signal, fs = librosa.core.load('./dataset/source1.wav')

    X_tmp = buffer(signal, L, M)

    # new matrix LxK
    X = np.vstack(X_tmp).T.astype('float')

    # ???
    # if X.shape[1] < L:
    #     print 'The number of frames %s has to be greater than its own length %s'%(X.shape[1],X.shape[0])
    #     sys.exit()

    alg = GAD(X, params)

    D, I = alg.iterative_GAD()

    X_t = np.dot(np.dot(D, D.T), X)

    s_rec = unbuffer(X_t, L - M)

    plt.figure(1)
    plt.title('Original signal')
    plt.plot(signal)

    plt.figure(2)
    plt.title('Reconstructed signal')
    plt.plot(s_rec)

    plt.show()

ext: py

reason: import numpy,from numpy


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-

__author__ = "Evangelos A. Dimopoulos, Evan K. Irving-Pease"
__copyright__ = "Copyright 2020, University of Oxford"
__email__ = "antonisdim41@gmail.com"
__license__ = "MIT"

import numpy as np
import pandas as pd
import sys

from sklearn import metrics

def roc_per_read(taxonomy, lkhd_matrix, output):
    """Function to calculate the classification metrics"""

    print(
        "Calculating the TP, FP and FN rates ...", file=sys.stderr,
    )

    # load the list of reads and the species they came from

    taxa = pd.read_csv(taxonomy, names=["species", "reads"], sep="\t")
    taxa["species"] = taxa["species"].str.replace(" ", "_")
    taxa["reads"] = taxa["reads"] + "-1"

    # load ts and tv matrix and calculate dark and grey matter

    ts_tv_matrix = pd.read_csv(
        lkhd_matrix, sep=",", usecols=["Taxon", "Read_ID", "Dirichlet_Assignment"]
    )

    # Small df for the grey matter reads
    ts_tv_group = ts_tv_matrix.groupby("Read_ID").sum().squeeze()
    grey_matter = ts_tv_group.where(ts_tv_group == 0).replace(0, 1).fillna(0)
    grey_df = grey_matter[grey_matter == 1].to_frame().reset_index()
    grey_df["Taxon"] = "Grey_Matter"

    # Small df for the dark matter reads
    dark_df = taxa[~taxa["reads"].isin(ts_tv_matrix["Read_ID"])].rename(
        columns={"reads": "Read_ID", "species": "Taxon"}
    )
    dark_df["Dirichlet_Assignment"] = 1
    dark_df["Taxon"] = "Dark_Matter"

    # split them in assigned and unassigned reads
    tax_assigned_reads = ts_tv_matrix[ts_tv_matrix["Dirichlet_Assignment"] != 0]
    assigned_reads = pd.concat([tax_assigned_reads, grey_df, dark_df], axis=0)

    # merge all the reads no matter where they were assigned
    assigned_comp = pd.merge(
        assigned_reads, taxa, left_on=["Read_ID"], right_on=["reads"]
    )

    # calculate classification stats

    values_true = assigned_comp["species"]
    values_pred = assigned_comp["Taxon"]

    precision = metrics.precision_score(
        values_true, values_pred, average="weighted", zero_division=0
    )
    recall = metrics.recall_score(
        values_true, values_pred, average="weighted", zero_division=0
    )
    f1 = metrics.f1_score(values_true, values_pred, average="weighted", zero_division=0)

    with open(output, "w") as fout:
        print(round(precision, 3), round(recall, 3), round(f1, 3), sep="\t", file=fout)

    print(
        "Done with the classification metrics.", file=sys.stderr,
    )

if __name__ == "__main__":

    # noinspection PyUnresolvedReferences
    sys.stderr = open(snakemake.log[0], "w")

    # noinspection PyUnresolvedReferences
    roc_per_read(
        taxonomy=snakemake.input.taxonomy,
        lkhd_matrix=snakemake.input.lkhd_matrix,
        output=snakemake.output[0],
    )

reason: import numpy


text:
# -*- coding: utf-8 -*-
"""
Created on Tue Mar  5 10:48:44 2019

@author: Nate
"""

import numpy as np
import matplotlib.pyplot as plt
import pdb

#fig2, axes2 = plt.subplots()

##############################################################
##############################################################
#defining constants

#number of trials
N = 100000

#α=0.2 to 2 in increments of 0.2

alpha = [i*.2 for i in range(1,11)]
rmax_list = [(3.2/alph) for alph in alpha]
r0 = np.array([0.01,0,0])
energy = []
std_dev = []

#sigma = [4*i for i in alpha]
#rmax = 1/(2*alpha)

##############################################################
##############################################################
#defining functions

def get_delta_r(rmax):
    delta_r = np.array([rmax*(np.random.uniform()-.5),rmax*(np.random.uniform()-.5),rmax*(np.random.uniform()-.5)])
    return(delta_r)

#evaluate ratio
def evaluate_ratio(r, delta_r):
    global count
    A = np.exp(-2*alpha[i]*np.linalg.norm(r))
    B = np.exp(-2*alpha[i]*np.linalg.norm(r+delta_r))
    ratio = B/A

    if ratio > np.random.uniform():
        count+=1
        return(True)
    else:
        return(False)

def generate_point(r, delta_r):
    if evaluate_ratio(r, delta_r) == True:
        return(r+delta_r)
    else:
        return(r)

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

for i in range(len(alpha)):

    #generating points
    r = np.array([r0])
    count = 0
    en = 0

    for j in range(N):
        #r2=np.random.exponential(.5/alpha[i])
        #en +=-alpha[i]**2 + alpha[i]/np.linalg.norm(r2) - 1/np.linalg.norm(r2)

        #pdb.set_trace()
        r = np.append(r, [generate_point(r[len(r)-1], get_delta_r(rmax_list[i]))], axis = 0)
        en +=-(alpha[i]**2)/2 + alpha[i]/np.linalg.norm(r[j]) - 1/np.linalg.norm(r[j])

    print(rmax_list[i], 'accepted: ',count, 'rejected: ',N-count)
    print('Acceptance Ratio: ', count/N)
    energy.append(en/N)

    to_sum = 0
    for j in range(N):
        to_sum += (-(alpha[i]**2)/2 + alpha[i]/np.linalg.norm(r[j]) - 1/np.linalg.norm(r[j])-energy[i])**2

    std_dev.append(np.sqrt(to_sum)/np.sqrt(N))

    '''
    if i == 1:
        axes2.hist([np.linalg.norm(i) for i in r], bins = 100, normed = True)
        axes2.plot([n*.1 for n in range(0,3000)]    , [2*alpha[i]*np.e**(-2*alpha[i]*n*.1)  for n in range(0,3000)])
'''

print(energy)

##############################################################
##############################################################
#Plotting

fig1, axes1 = plt.subplots()

#axes1.scatter(alpha, [3/2*num-2**(3/2)/(np.pi**.5)*num**.5 for num in alpha])
#axes1.scatter(8/(9*np.pi), 3/2*(8/(9*np.pi))-2**(3/2)/(np.pi**.5)*(8/(9*np.pi))**.5)
#axes1.plot(alpha, energy)
axes1.errorbar(alpha, energy, yerr=std_dev, fmt='o', linestyle = '-')
axes1.plot(alpha, [(num**2)/2-num for num in alpha])

axes1.set_ylabel('Energy')
axes1.set_xlabel('alpha')
axes1.set_title("Energy vs alpha", va='bottom')

plt.show()

max_stars_repo_path: HW7/Problem 2.py

max_issues_repo_path: HW7/Problem 2.py

max_forks_repo_path: HW7/Problem 2.py


text:
import logging; log = logging.getLogger(__name__)
from importlib import resources
import struct
import math
import sys
import gl
import numpy as np

class BoxRenderer(gl.Pipeline):
    """Subprogram for rendering boxes."""
    MAX_BOXES = 8192
    vtxBufferFmt = "6f" # (x,y,z) A, (x,y,z) B
    colorBufferFmt = "4f" # (r,g,b,a) per box

    def __init__(self, parent, shader):
        self.parent = parent
        #self.shader = BoxShader(self.parent.ctx)
        self.shader = shader

        super().__init__(parent.ctx,
            vertex_shader   = self.shader,
            geometry_shader = self.shader,
            fragment_shader = parent.fragShader,
        )

        self._vtxSize   = struct.calcsize(self.vtxBufferFmt)
        self._colorSize = struct.calcsize(self.colorBufferFmt)
        self.bufVtxs   = self.ctx.Buffer(self._vtxSize * self.MAX_BOXES, True)
        self.bufColors = self.ctx.Buffer(self._colorSize * self.MAX_BOXES, True)
        self._nBoxes = 0

    def addBox(self, pointA, pointB, color):
        """Add a box with one corner at pointA and the other at pointB."""
        if self._nBoxes >= self.MAX_BOXES: return

        x1, y1, z1 = pointA
        x2, y2, z2 = pointB
        r, g, b, a = color

        vtxData = struct.pack(self.vtxBufferFmt, x1, y1, z1, x2, y2, z2)
        colData = struct.pack(self.colorBufferFmt, r, g, b, a)
        vtxOffs = self._nBoxes * self._vtxSize
        colOffs = self._nBoxes * self._colorSize
        self.bufVtxs.write(vtxData, vtxOffs)
        self.bufColors.write(colData, colOffs)
        self._nBoxes += 1

    def reset(self):
        self._nBoxes = 0

    def run(self):
        if self._nBoxes == 0: return
        self.ctx.glBlendFunc(self.ctx.GL_SRC_ALPHA, self.ctx.GL_ONE_MINUS_SRC_ALPHA)
        self.ctx.glEnable(self.ctx.GL_BLEND)
        with self:
            self._bindBuffers()
            self.shader.vao.render(self.ctx.GL_POINTS, count=self._nBoxes)

    def setMtxs(self, projection, modelView):
        self.shader.setMtxs(projection, modelView)

    def _bindBuffers(self):
        with self.shader as shader:
            # in vec3 vertA
            shader.attribs['vertA'].bindBuffer(
                self.bufVtxs, self.ctx.GL_FLOAT, 3,
                stride=6 * gl.Util.SIZEOF_FLOAT,
                offset=0 * gl.Util.SIZEOF_FLOAT)

            # in vec3 vertB
            shader.attribs['vertB'].bindBuffer(
                self.bufVtxs, self.ctx.GL_FLOAT, 3,
                stride=6 * gl.Util.SIZEOF_FLOAT,
                offset=3 * gl.Util.SIZEOF_FLOAT)

            # in vec4 vtxColor
            shader.attribs['vtxColor'].bindBuffer(
                self.bufColors, self.ctx.GL_FLOAT, 4,
                stride=4 * gl.Util.SIZEOF_FLOAT,
                offset=0 * gl.Util.SIZEOF_FLOAT)


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

def normalized_columns_initializer(std=1.0):
    def _initializer(shape, dtype=None, partition_info=None):
        out = np.random.randn(*shape).astype(np.float32)
        out *= std / np.sqrt(np.square(out).sum(axis=0, keepdims=True))
        return tf.constant(out)
    return _initializer

def flatten(x):
    return tf.reshape(x, [-1, np.prod(x.get_shape().as_list()[1:])])

def conv2d(x, num_filters, name, filter_size=(3, 3), stride=(1, 1), pad="SAME", dtype=tf.float32, collections=None):
    with tf.variable_scope(name):
        stride_shape = [1, stride[0], stride[1], 1]
        filter_shape = [filter_size[0], filter_size[1], int(x.get_shape()[3]), num_filters]

        # there are "num input feature maps * filter height * filter width"
        # inputs to each hidden unit
        fan_in = np.prod(filter_shape[:3])
        # each unit in the lower layer receives a gradient from:
        # "num output feature maps * filter height * filter width" /
        #   pooling size
        fan_out = np.prod(filter_shape[:2]) * num_filters
        # initialize weights with random weights
        w_bound = np.sqrt(6. / (fan_in + fan_out))

        w = tf.get_variable("W", filter_shape, dtype, tf.random_uniform_initializer(-w_bound, w_bound),
                            collections=collections)
        b = tf.get_variable("b", [1, 1, 1, num_filters], initializer=tf.constant_initializer(0.0),
                            collections=collections)
        return tf.nn.conv2d(x, w, stride_shape, pad) + b

def linear(x, size, name, initializer=None, bias_init=0):
    w = tf.get_variable(name + "/w", [x.get_shape()[1], size], initializer=initializer)
    b = tf.get_variable(name + "/b", [size], initializer=tf.constant_initializer(bias_init))
    return tf.matmul(x, w) + b

def categorical_sample(logits, d):
    value = tf.squeeze(tf.multinomial(logits - tf.reduce_max(logits, [1], keep_dims=True), 1), [1])
    return tf.one_hot(value, d)

class LSTMPolicy(object):
    def __init__(self, ob_space, ac_space):
        self.x = x = tf.placeholder(tf.float32, [None] + list(ob_space))

        for i in range(4):
            x = tf.nn.elu(conv2d(x, 32, "l{}".format(i + 1), [3, 3], [2, 2]))
        # introduce a "fake" batch dimension of 1 after flatten so that we can do LSTM over time dim
        x = tf.expand_dims(flatten(x), [0])

        size = 256
        lstm = rnn.rnn_cell.BasicLSTMCell(size, state_is_tuple=True)
        self.state_size = lstm.state_size
        step_size = tf.shape(self.x)[:1]

        c_init = np.zeros((1, lstm.state_size.c), np.float32)
        h_init = np.zeros((1, lstm.state_size.h), np.float32)
        self.state_init = [c_init, h_init]
        c_in = tf.placeholder(tf.float32, [1, lstm.state_size.c])
        h_in = tf.placeholder(tf.float32, [1, lstm.state_size.h])
        self.state_in = [c_in, h_in]

        state_in = rnn.rnn_cell.LSTMStateTuple(c_in, h_in)
        lstm_outputs, lstm_state = tf.nn.dynamic_rnn(
            lstm, x, initial_state=state_in, sequence_length=step_size,
            time_major=False)
        lstm_c, lstm_h = lstm_state
        x = tf.reshape(lstm_outputs, [-1, size])
        self.logits = linear(x, ac_space, "action", normalized_columns_initializer(0.01))
        self.vf = tf.reshape(linear(x, 1, "value", normalized_columns_initializer(1.0)), [-1])
        self.state_out = [lstm_c[:1, :], lstm_h[:1, :]]
        self.sample = categorical_sample(self.logits, ac_space)[0, :]
        self.var_list = tf.get_collection(tf.GraphKeys.TRAINABLE_VARIABLES, tf.get_variable_scope().name)

        self._last_state = self.get_initial_features()

    def get_initial_features(self):
        return self.state_init

    def act(self, ob, c, h, track_state=True):
        sess = tf.get_default_session()
        fetched = sess.run([self.sample, self.vf, self.logits] + self.state_out,
                           {self.x: [ob], self.state_in[0]: c, self.state_in[1]: h})

        if track_state:
            self._last_state = fetched[3:]

        return fetched

    def value(self, ob, c, h):
        sess = tf.get_default_session()
        return sess.run(self.vf, {self.x: [ob], self.state_in[0]: c, self.state_in[1]: h})[0]

max_stars_repo_name: sunblaze-ucb/rl-attack-vf

max_issues_repo_name: sunblaze-ucb/rl-attack-vf

max_forks_repo_name: sunblaze-ucb/rl-attack-vf


text:
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

plt.switch_backend('agg')

log1 = open('data-log/measure/atta-loss-1.log')
log2 = open('data-log/measure/atta-loss-2.log')
log3 = open('data-log/measure/atta-loss-3.log')
log4 = open('data-log/measure/atta-loss-4.log')
log5 = open('data-log/measure/atta-loss-5.log')
log6 = open('data-log/measure/atta-loss-6.log')
log7 = open('data-log/measure/atta-loss-7.log')
log8 = open('data-log/measure/atta-loss-8.log')
log9 = open('data-log/measure/atta-loss-9.log')
log10 = open('data-log/measure/atta-loss-10.log')

data1 = []
data2 = []
data3 = []
data4 = []
data5 = []
data6 = []
data7 = []
data8 = []
data9 = []
data10 = []

log_lines1 = log1.readlines()
log_lines2 = log2.readlines()
log_lines3 = log3.readlines()
log_lines4 = log4.readlines()
log_lines5 = log5.readlines()
log_lines6 = log6.readlines()
log_lines7 = log7.readlines()
log_lines8 = log8.readlines()
log_lines9 = log9.readlines()
log_lines10 = log10.readlines()
for i in range(len(log_lines1)):
  data1.append([eval(j) for j in log_lines1[i].split(' ')])
  data2.append([eval(j) for j in log_lines2[i].split(' ')])
  data3.append([eval(j) for j in log_lines3[i].split(' ')])
  data4.append([eval(j) for j in log_lines4[i].split(' ')])
  data5.append([eval(j) for j in log_lines5[i].split(' ')])
  data6.append([eval(j) for j in log_lines6[i].split(' ')])
  data7.append([eval(j) for j in log_lines7[i].split(' ')])
  data8.append([eval(j) for j in log_lines8[i].split(' ')])
  data9.append([eval(j) for j in log_lines9[i].split(' ')])
  data10.append([eval(j) for j in log_lines10[i].split(' ')])

print(len(data1))

x = np.array([i[1] for i in data1]) + 1

adv_loss1 = np.array([i[2] for i in data1])
adv_loss2 = np.array([i[2] for i in data2])
adv_loss3 = np.array([i[2] for i in data3])
adv_loss4 = np.array([i[2] for i in data4])
adv_loss5 = np.array([i[2] for i in data5])
adv_loss6 = np.array([i[2] for i in data6])
adv_loss7 = np.array([i[2] for i in data7])
adv_loss8 = np.array([i[2] for i in data8])
adv_loss9 = np.array([i[2] for i in data9])
adv_loss10 = np.array([i[2] for i in data10])

current_palette = sns.color_palette()

plt.plot(x, adv_loss1, color=current_palette[0], lw=2)
plt.plot(x, adv_loss2, color=current_palette[1], lw=2)
plt.plot(x, adv_loss3, color=current_palette[2], lw=2)
plt.plot(x, adv_loss4, color=current_palette[3], lw=2)
plt.plot(x, adv_loss5, color=current_palette[4], lw=2)
plt.plot(x, adv_loss6, color=current_palette[5], lw=2)
plt.plot(x, adv_loss7, color=current_palette[6], lw=2)
plt.plot(x, adv_loss8, color=current_palette[7], lw=2)
plt.plot(x, adv_loss9, color=current_palette[8], lw=2)
plt.plot(x, adv_loss10, color=current_palette[9], lw=2)
plt.yticks([0, 0.1, 0.2, 0.3, 0.4])

plt.xlabel("Attack iterations", fontsize=18)
plt.ylabel("Loss value", fontsize=18)
plt.tick_params(labelsize=13)

# plt.legend(fontsize='x-large')

plt.savefig('data-pic/multiple-checkpoints-attack-loss.pdf')


text:
import numpy
from ..Fragments import Fragments
from ..typing import SpectrumType

def select_by_relative_intensity(spectrum_in: SpectrumType, intensity_from: float = 0.0,
                                 intensity_to: float = 1.0) -> SpectrumType:
    """Keep only peaks within set relative intensity range (keep if
    intensity_from >= intensity >= intensity_to).

    Parameters
    ----------
    intensity_from:
        Set lower threshold for relative peak intensity. Default is 0.0.
    intensity_to:
        Set upper threshold for relative peak intensity. Default is 1.0.
    """
    if spectrum_in is None:
        return None

    spectrum = spectrum_in.clone()

    assert intensity_from >= 0.0, "'intensity_from' should be larger than or equal to 0."
    assert intensity_to <= 1.0, "'intensity_to' should be smaller than or equal to 1.0."
    assert intensity_from <= intensity_to, "'intensity_from' should be smaller than or equal to 'intensity_to'."

    if len(spectrum.peaks) > 0:
        scale_factor = numpy.max(spectrum.peaks.intensities)
        normalized_intensities = spectrum.peaks.intensities / scale_factor
        condition = numpy.logical_and(intensity_from <= normalized_intensities, normalized_intensities <= intensity_to)
        spectrum.peaks = Fragments(mz=spectrum.peaks.mz[condition],
                                   intensities=spectrum.peaks.intensities[condition])

    return spectrum

max_stars_repo_name: maximskorik/matchms

max_issues_repo_name: maximskorik/matchms

max_forks_repo_name: maximskorik/matchms


text:
import ecole
import numpy as np
import pyscipopt
from mllocalbranch_fromfiles import RegressionInitialK
from utility import instancetypes, instancesizes, incumbent_modes, lbconstraint_modes

# instance_type = instancetypes[1]
instance_size = instancesizes[0]
# incumbent_mode = 'firstsol'
lbconstraint_mode = 'symmetric'
samples_time_limit = 3
node_time_limit = 10

total_time_limit = 60
reset_k_at_2nditeration = True

for i in range(0, 1):
    instance_type = instancetypes[i]
    if instance_type == instancetypes[0]:
        lbconstraint_mode = 'asymmetric'
    else:
        lbconstraint_mode = 'symmetric'
    for j in range(1, 2):
        incumbent_mode = incumbent_modes[j]
        for k in range(0, 1):
            test_instance_size = instancesizes[k]

            print(instance_type + test_instance_size)
            print(incumbent_mode)
            print(lbconstraint_mode)

            regression_init_k = RegressionInitialK(instance_type, instance_size, lbconstraint_mode, incumbent_mode, seed=100)

            regression_init_k.generate_k_samples(t_limit=samples_time_limit, instance_size=instance_size)

            # regression_init_k.generate_regression_samples(t_limit=samples_time_limit, instance_size=instance_size)

            # regression_init_k.execute_regression(lr=0.00002, n_epochs=11) # setcovering small: lr=0.00002; capa-small: samne; independentset-small: first: lr=0.00002, root: lr=0.00003

            # regression_init_k.evaluate_localbranching(test_instance_size=test_instance_size,train_instance_size='-small', total_time_limit=total_time_limit, node_time_limit=node_time_limit, reset_k_at_2nditeration=reset_k_at_2nditeration)

            # regression_init_k.solve2opt_evaluation(test_instance_size='-small')

            # regression_init_k.primal_integral(test_instance_size=test_instance_size, total_time_limit=total_time_limit, node_time_limit=node_time_limit)


text:
import time
import grpc
import numpy as np
import platform
import sys
import rpc_client

from protos import jetsonrpc_pb2_grpc
from protos import jetsonrpc_pb2
from utils.protocol import encode_values
# from keyboard_driver import keyboard_val_gen

gamepad_running = False

detected_platform = None
if sys.platform == "win32" or sys.platform == "cygwin": detected_platform = "WINDOWS"
if sys.platform.startswith("linux"): detected_platform = "LINUX"

def gamepad_val_gen():
    if detected_platform == "WINDOWS":
        print("gamepad_encoder starting windows driver")
        from gamepad_driver_windows import get_gamepad_values, joyGetPosEx, p_info
    elif detected_platform == "LINUX":
        print("gamepad_encoder starting linux driver")
        import gamepad_driver_linux
        from gamepad_driver_linux import get_gamepad_values
        gamepad_driver_linux.start()
    else:
        print("Error: gamepad_encoder system platform not recognized! Aborting...")
        return

    prev_encoded_val = None
    while gamepad_running:
        if detected_platform == "WINDOWS" and joyGetPosEx(0, p_info) != 0:
            # Running this joyGetPosEx check is necessary for the windows code to update the gamepad values
            print("Error: gamepad_encoder gamepad disconnected")
            break

        values = get_gamepad_values()
        encoded_val = encode_values(*values) # python "spread"/"splat" operator
        if prev_encoded_val == encoded_val: # don't send messages if gamepad value hasn't changed
            time.sleep(0.05) # TODO adjust this rate to balance saving bandwidth and maintaining robot control. This could also be a parameter that gets set and passed in gui.py somehow
        else:
            print("gamepad_encoder sending motor values:", values)
            prev_encoded_val = encoded_val
            yield jetsonrpc_pb2.DDCommand(values=encoded_val)

    if detected_platform == "LINUX": # only linux driver needs a "stop()" command
        gamepad_driver_linux.stop()

def dummy_val_gen():
    import random
    while True:
        m1 = random.randint(0, 200)
        m2 = random.randint(0, 200)
        bucket_ladder_angle = random.randint(0, 200)
        bucket_ladder_extension = random.randint(0, 200)
        bucket_ladder_chain = random.randint(0, 200)
        deposit_angle = random.choice([0, 100, 200])
        conveyor = random.choice([0, 100, 200])

        yield jetsonrpc_pb2.DDCommand(values=
            encode_values(m1, m2, bucket_ladder_angle, bucket_ladder_extension, bucket_ladder_chain, deposit_angle, conveyor)
        )
        time.sleep(0.1)

def start(stub):
    print("gamepad_encoder starting...")
    global gamepad_running
    gamepad_running = True
    response = stub.SendDDCommand(gamepad_val_gen())
    # response = stub.SendDDCommand(dummy_val_gen()) # for sending fake data - will probably destroy robot if actually used
    print(response)

def stop():
    print("gamepad_encoder stopping...")
    global gamepad_running
    gamepad_running = False

if __name__ == '__main__':
    import network_config

    with grpc.insecure_channel("{}:{}".format(network_config.HOST, network_config.PORT)) as channel:
        try:
            stub = jetsonrpc_pb2_grpc.JetsonRPCStub(channel)
            print("Changing to drive state DIRECT_DRIVE...")
            rpc_client.change_drive_state(stub, jetsonrpc_pb2.DriveStateEnum.DIRECT_DRIVE)

            start(stub)
        except KeyboardInterrupt:
            stop()
        finally:
            print("Changing to drive state IDLE...")
            rpc_client.change_drive_state(stub, jetsonrpc_pb2.DriveStateEnum.IDLE)

max_stars_repo_name: MARS-UVA/mars-ui


text:
import numpy as np

def findCosineDistance(source_representation, test_representation):
    a = np.matmul(np.transpose(source_representation), test_representation)
    b = np.sum(np.multiply(source_representation, source_representation))
    c = np.sum(np.multiply(test_representation, test_representation))
    return 1 - (a / (np.sqrt(b) * np.sqrt(c)))

def findEuclideanDistance(source_representation, test_representation):
    euclidean_distance = source_representation - test_representation
    euclidean_distance = np.sum(np.multiply(euclidean_distance, euclidean_distance))
    euclidean_distance = np.sqrt(euclidean_distance)
    return euclidean_distance

def l2_normalize(x):
    return x / np.sqrt(np.sum(np.multiply(x, x)))

def findThreshold(model_name, distance_metric):

	base_threshold = {'cosine': 0.40, 'euclidean': 0.55, 'euclidean_l2': 0.75}

	thresholds = {
		'VGG-Face': {'cosine': 0.3828998554497957, 'euclidean': 0.5837138535827399, 'euclidean_l2': 0.8437118818678639},
		'OpenFace': {'cosine': 0.10, 'euclidean': 0.55, 'euclidean_l2': 0.55},
		'Facenet':  {'cosine': 0.6142496640305134, 'euclidean': 12.794783464670182, 'euclidean_l2': 1.0578221621432087},
		'DeepFace': {'cosine': 0.23, 'euclidean': 64, 'euclidean_l2': 0.64},
		'DeepID': 	{'cosine': 0.015, 'euclidean': 45, 'euclidean_l2': 0.17},
		'Dlib': 	{'cosine': 0.08838775356588767, 'euclidean': 0.5821402502891807, 'euclidean_l2': 0.4102203243676945},
		'ArcFace':  {'cosine': 0.620748013800026, 'euclidean': 5.081651272556998, 'euclidean_l2': 1.0885364669832316}
		}

	threshold = thresholds.get(model_name, base_threshold).get(distance_metric, 0.4)

	return threshold

max_stars_repo_name: thisKK/deepface

max_issues_repo_name: thisKK/deepface

max_forks_repo_name: thisKK/deepface


text:
import pandas as pd
import numpy as np
import asyncio
import os,sys,datetime

import scipy as sp
import statsmodels.tsa.stattools as sts
import matplotlib.pyplot as plt
import statsmodels.api as sm
import talib
import pickle
from sklearn.model_selection import train_test_split
from sklearn.base import BaseEstimator, TransformerMixin
import xgboost as xgb
from sklearn import metrics

from sklearn.impute import KNNImputer

reason: import numpy,import scipy,import statsmodels


text:
import numpy as np
import sys
# See https://github.com/YuyangL/SOWFA-PostProcess
sys.path.append('/home/yluan/Documents/SOWFA PostProcessing/SOWFA-Postprocess')
from FieldData import FieldData
from Postprocess.OutlierAndNoveltyDetection import InputOutlierDetection
from Preprocess.Tensor import processReynoldsStress, getBarycentricMapData, expandSymmetricTensor, contractSymmetricTensor
from Utility import interpolateGridData
from joblib import load, dump
import time as t
from PlottingTool import BaseFigure, Plot2D, Plot2D_Image, PlotContourSlices3D
from scipy import ndimage
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from copy import copy
import mpl_toolkits.mplot3d.art3d as art3d
import os

"""
User Inputs, Anything Can Be Changed Here
"""
# Name of the flow case in both RANS and LES
rans_case_name = 'RANS_Re10595'  # str
test_case_name = 'RANS_Re5600'  # str
les_case_name = 'LES_Breuer/Re_5600'  # str
# LES data name to read
les_data_name = 'Hill_Re_5600_Breuer.csv'  # str
# Absolute directory of this flow case
casedir = '/media/yluan/DNS/PeriodicHill'  # str
# Which time to extract input and output for ML
time = '5000'  # str/float/int or 'last'
# Interpolation method when interpolating mesh grids
interp_method = "nearest"  # "nearest", "linear", "cubic"

"""
Machine Learning Settings
"""
isoforest_name = 'IsoForest'
estimator_name = 'tbdt'
fs = 'grad(TKE)_grad(p)'
# Seed value for reproducibility
seed = 123
outlier_percent1, outlier_percent2, outlier_percent3, outlier_percent4, outlier_percent5 = \
    0.1, 0.08, 0.06, 0.04, 0.02
load_isoforest, save_estimator = True, False

"""
Plot Settings
"""
plot_tb, plot_g, plot_x = False, False, False
# When plotting, the mesh has to be uniform by interpolation, specify target size
uniform_mesh_size = 1e6  # int
# Limit for bij plot
bijlims = (-1/2., 2/3.)  # (float, float)
contour_lvl = 50  # int
alpha = 0.6  # int, float [0, 1]
gray = (80/255.,)*3
figheight_multiplier = 1.1
# Save anything when possible
save_fields = True  # bool
# Save figures and show figures
save_fig, show = True, False  # bool; bool
if save_fig:
    # Figure extension and DPI
    ext, dpi = 'png', 1000  # str; int

"""
Process User Inputs, No Need to Change
"""
# Average fields of interest for reading and processing
fields = ('U', 'k', 'p', 'omega',
          'grad_U', 'grad_k', 'grad_p')
# if fs == "grad(TKE)_grad(p)":
#     fields = ('U', 'k', 'p', 'omega',
#               'grad_U', 'grad_k', 'grad_p')
# elif fs == "grad(TKE)":
#     fields = ('k', 'omega',
#               'grad_U', 'grad_k')
# elif fs == "grad(p)":
#     fields = ('U', 'k', 'p', 'omega',
#               'grad_U', 'grad_p')
# else:
#     fields = ('k', 'omega',
#               'grad_U')

if estimator_name == 'tbdt': estimator_name = 'TBDT'
# Ensemble name of fields useful for Machine Learning
ml_field_ensemble_name = 'ML_Fields_' + rans_case_name
# Initialize case object
case = FieldData(casename=rans_case_name, casedir=casedir, times=time, fields=fields)
case_test = FieldData(casename=test_case_name, casedir=casedir, times=time, fields=fields)

"""
Machine Learning with Isolation Forest
"""
list_data_train = case.readPickleData(time, 'list_data_train_seed' + str(seed))
list_data_test = case_test.readPickleData(time, 'list_data_test_seed' + str(seed))
if '10595' not in test_case_name: list_data_train = list_data_test

cc_train, cc_test = list_data_train[0], list_data_test[0]
ccx_train, ccy_train, ccz_train = cc_train[:, 0], cc_train[:, 1], cc_train[:, 2]
ccx_test, ccy_test, ccz_test = cc_test[:, 0], cc_test[:, 1], cc_test[:, 2]
x_train, y_train, tb_train = list_data_train[1:4]
x_test, y_test, tb_test = list_data_test[1:4]

isoforest = load(case.result_paths[time] + isoforest_name + '.joblib') if load_isoforest else None
isoforest2 = load(case.result_paths[time] + isoforest_name + '2.joblib') if load_isoforest else None
isoforest3 = load(case.result_paths[time] + isoforest_name + '2.joblib') if load_isoforest else None
isoforest4 = load(case.result_paths[time] + isoforest_name + '2.joblib') if load_isoforest else None
isoforest5 = load(case.result_paths[time] + isoforest_name + '5.joblib') if load_isoforest else None

_, _, _, _, outliers, isoforest = InputOutlierDetection(x_train, x_test, y_train, y_test, outlier_percent=outlier_percent1, randstate=seed, isoforest=isoforest, n_estimators=1000)
_, _, _, _, outliers2, isoforest2 = InputOutlierDetection(x_train, x_test, y_train, y_test, outlier_percent=outlier_percent2, randstate=seed, isoforest=isoforest2, n_estimators=1000)
_, _, _, _, outliers3, isoforest3 = InputOutlierDetection(x_train, x_test, y_train, y_test, outlier_percent=outlier_percent3, randstate=seed, isoforest=isoforest3, n_estimators=1000)
_, _, _, _, outliers4, isoforest4 = InputOutlierDetection(x_train, x_test, y_train, y_test, outlier_percent=outlier_percent4, randstate=seed, isoforest=isoforest4, n_estimators=1000)
_, _, _, _, outliers5, isoforest5 = InputOutlierDetection(x_train, x_test, y_train, y_test, outlier_percent=outlier_percent5, randstate=seed, isoforest=isoforest5, n_estimators=1000)

if save_estimator and not load_isoforest:
    dump(isoforest, case.result_paths[time] + isoforest_name + '.joblib')
    dump(isoforest2, case.result_paths[time] + isoforest_name + '2.joblib')
    dump(isoforest3, case.result_paths[time] + isoforest_name + '3.joblib')
    dump(isoforest4, case.result_paths[time] + isoforest_name + '4.joblib')
    dump(isoforest5, case.result_paths[time] + isoforest_name + '5.joblib')

x_train_out = outliers['xtrain']
x_test_out = outliers['xtest']
y_train_out = outliers['ytrain']
y_test_out = outliers['ytest']
anomaly_idx_train = outliers['anomaly_idx_train']
anomaly_idx_test = outliers['anomaly_idx_test']
anomaly_idx_train2 = outliers2['anomaly_idx_train']
anomaly_idx_test2 = outliers2['anomaly_idx_test']
anomaly_idx_train3 = outliers3['anomaly_idx_train']
anomaly_idx_test3 = outliers3['anomaly_idx_test']
anomaly_idx_train4 = outliers4['anomaly_idx_train']
anomaly_idx_test4 = outliers4['anomaly_idx_test']
anomaly_idx_train5 = outliers5['anomaly_idx_train']
anomaly_idx_test5 = outliers5['anomaly_idx_test']

"""
Postprocess Machine Learning Predictions
"""
t0 = t.time()
_, eigval_test, _ = processReynoldsStress(y_test, make_anisotropic=False, realization_iter=0)
_, eigval_train, _ = processReynoldsStress(y_train, make_anisotropic=False, realization_iter=0)
t1 = t.time()
print('\nFinished processing Reynolds stress in {:.4f} s'.format(t1 - t0))

t0 = t.time()
xy_bary_test, rgb_bary_test = getBarycentricMapData(eigval_test)
xy_bary_train, rgb_bary_train = getBarycentricMapData(eigval_train)
rgb_bary_test_out = rgb_bary_test.copy()
rgb_bary_test_out[anomaly_idx_test], rgb_bary_test_out[anomaly_idx_test2], rgb_bary_test_out[anomaly_idx_test3], rgb_bary_test_out[anomaly_idx_test4], rgb_bary_test_out[anomaly_idx_test5] \
    = (0.4,)*3, (0.3,)*3, (0.2,)*3, (0.1,)*3, (0,)*3
rgb_bary_train_out = rgb_bary_train.copy()
rgb_bary_train_out[anomaly_idx_train], rgb_bary_train_out[anomaly_idx_train2], rgb_bary_train_out[anomaly_idx_train3], rgb_bary_train_out[anomaly_idx_train4], rgb_bary_train_out[anomaly_idx_train5] \
    = (0.4,)*3, (0.3,)*3, (0.2,)*3, (0.1,)*3, (0,)*3
t1 = t.time()
print('\nFinished getting Barycentric map data in {:.4f} s'.format(t1 - t0))

t0 = t.time()
ccx_test_mesh, ccy_test_mesh, _, rgb_bary_test_mesh = interpolateGridData(ccx_test, ccy_test, rgb_bary_test, mesh_target=uniform_mesh_size, interp=interp_method, fill_val=89/255.)
_, _, _, rgb_bary_test_out_mesh = interpolateGridData(ccx_test, ccy_test, rgb_bary_test_out, mesh_target=uniform_mesh_size, interp=interp_method, fill_val=89/255.)
ccx_train_mesh, ccy_train_mesh, _, rgb_bary_train_mesh = interpolateGridData(ccx_train, ccy_train, rgb_bary_train, mesh_target=uniform_mesh_size, interp=interp_method, fill_val=89/255.)
_, _, _, rgb_bary_train_out_mesh = interpolateGridData(ccx_train, ccy_train, rgb_bary_train_out, mesh_target=uniform_mesh_size, interp=interp_method, fill_val=89/255.)
t1 = t.time()
print('\nFinished interpolating mesh data for barycentric map in {:.4f} s'.format(t1 - t0))

"""
Visualize Outliers and Novelties in Barycentric Map
"""
rgb_bary_test_mesh = ndimage.rotate(rgb_bary_test_mesh, 90)
rgb_bary_train_mesh = ndimage.rotate(rgb_bary_train_mesh, 90)
rgb_bary_test_out_mesh = ndimage.rotate(rgb_bary_test_out_mesh, 90)
rgb_bary_train_out_mesh = ndimage.rotate(rgb_bary_train_out_mesh, 90)

xlabel, ylabel = (r'$x$ [m]', r'$y$ [m]')
geometry = np.genfromtxt(casedir + '/' + rans_case_name + '/'  + "geometry.csv", delimiter=",")[:, :2]
# Test barycentric map
figname = 'barycentric_periodichill_test_out'
bary_map = BaseFigure((None,), (None,), name=figname, xlabel=xlabel,
                      ylabel=ylabel, save=save_fig, show=show,
                      figdir=case.result_paths[time],
                      figheight_multiplier=0.7)
path = Path(geometry)
patch = PathPatch(path, linewidth=0., facecolor=bary_map.gray)
# patch is considered "a single artist" so have to make copy to use more than once
patches = []
for _ in range(236):
    patches.append(copy(patch))

patches = iter(patches)
extent_test = (ccx_test.min(), ccx_test.max(), ccy_test.min(), ccy_test.max())
extent_train = (ccx_train.min(), ccx_train.max(), ccy_train.min(), ccy_train.max())
bary_map.initializeFigure()
bary_map.axes.imshow(rgb_bary_test_mesh, origin='upper', aspect='equal', extent=extent_test)
bary_map.axes.imshow(rgb_bary_test_out_mesh, origin='upper', aspect='equal', extent=extent_test, alpha=alpha)
bary_map.axes.set_xlabel(bary_map.xlabel)
bary_map.axes.set_ylabel(bary_map.ylabel)
bary_map.axes.add_patch(next(patches))
if save_fig:
    plt.savefig(case.result_paths[time] + figname + '.' + ext, dpi=dpi)

plt.close()
# Train barycentric map
bary_map.name = 'barycentric_periodichill_train_out'
bary_map.initializeFigure()
bary_map.axes.imshow(rgb_bary_train_mesh, origin='upper', aspect='equal', extent=extent_train)
bary_map.axes.imshow(rgb_bary_train_out_mesh, origin='upper', aspect='equal', extent=extent_test, alpha=alpha)
bary_map.axes.set_xlabel(bary_map.xlabel)
bary_map.axes.set_ylabel(bary_map.ylabel)
bary_map.axes.add_patch(next(patches))
if save_fig:
    plt.savefig(case.result_paths[time] + bary_map.name + '.' + ext, dpi=dpi)

plt.close()

"""
Visualize Tij Along with X Outlier and Novelties
"""
if plot_tb:
    t0 = t.time()
    # Create Tij mesh ensemble, shape (n_x, n_y, n_bases, n_outputs)
    tb_test_mesh = np.empty((ccx_test_mesh.shape[0], ccy_test_mesh.shape[1], tb_test.shape[2], tb_test.shape[1]))
    tb_train_mesh = np.empty((ccx_train_mesh.shape[0], ccy_train_mesh.shape[1], tb_train.shape[2], tb_train.shape[1]))
    # Go through each Tij component, interpolate bases
    for i in range(tb_test.shape[1]):
        _, _, _, tb_test_mesh_i = interpolateGridData(ccx_test, ccy_test, tb_test[:, i], mesh_target=uniform_mesh_size,
                                                      interp=interp_method)
        _, _, _, tb_train_mesh_i = interpolateGridData(ccx_train, ccy_train, tb_train[:, i],
                                                       mesh_target=uniform_mesh_size, interp=interp_method)
        # Assign each component to an ensemble
        tb_test_mesh[..., i], tb_train_mesh[..., i] = tb_test_mesh_i, tb_train_mesh_i

    # Swap n_outputs and n_bases axis back so that (grid, n_outputs, n_bases)
    tb_test_mesh = np.swapaxes(tb_test_mesh, -2, -1)
    tb_train_mesh = np.swapaxes(tb_train_mesh, -2, -1)
    t1 = t.time()
    print('\nFinished interpolating mesh data for Tij in {:.4f} s'.format(t1 - t0))

    # Make an individual directory
    tijdir = case.result_paths[time] + '/Tij'
    os.makedirs(tijdir, exist_ok=True)
    comps = ('11', '12', '13', '22', '23', '33')
    tbcomps = ('T_{11}', 'T_{12}', 'T_{13}', 'T_{22}', 'T_{23}', 'T_{33}')
    extent_test = (ccx_test.min(), ccx_test.max(), ccy_test.min(), ccy_test.max())
    extent_train = (ccx_train.min(), ccx_train.max(), ccy_train.min(), ccy_train.max())
    # Go through every component then every basis
    for i in range(tb_test.shape[1]):
        for j in range(tb_test.shape[2]):
            zlabel = '${0}^{1}$ [-]'.format(tbcomps[i], '{(' + str(j + 1) + ')}')
            # Tij is limited individually so that only test and train have the same limit
            tblim = (min(min(tb_test[:, i, j]), min(tb_train[:, i, j])),
                     max(max(tb_test[:, i, j]), max(tb_train[:, i, j])))
            # First, plot test Tij
            figname = 'T' + comps[i] + '_' + str(j + 1) + '_test'
            # Initialize PLot2D_Image object
            tbplot_test = Plot2D_Image(val=tb_test_mesh[..., i, j],
                                       name=figname,
                                       xlabel=xlabel, ylabel=ylabel, val_label=zlabel,
                                       save=save_fig,
                                       show=show,
                                       figdir=tijdir,
                                       rotate_img=True,
                                       extent=extent_test,
                                       figwidth='1/3',
                                       val_lim=tblim,
                                       figheight_multiplier=figheight_multiplier)
            tbplot_test.initializeFigure()
            tbplot_test.plotFigure()
            tbplot_test.axes.ticklabel_format(scilimits=(-4, 5))
            tbplot_test.axes.add_patch(next(patches))
            tbplot_test.finalizeFigure()

            # Then plot train Tij
            figname = 'T' + comps[i] + '_' + str(j + 1) + '_train'
            tbplot_train = Plot2D_Image(val=tb_train_mesh[..., i, j],
                                       name=figname,
                                       xlabel=xlabel, ylabel=ylabel, val_label=zlabel,
                                       save=save_fig,
                                       show=show,
                                       figdir=tijdir,
                                       rotate_img=True,
                                       extent=extent_train,
                                       figwidth='1/3',
                                       val_lim=tblim,
                                        figheight_multiplier=figheight_multiplier)
            tbplot_train.initializeFigure()
            tbplot_train.plotFigure()
            tbplot_train.axes.add_patch(next(patches))
            tbplot_train.finalizeFigure()

"""
Visualize g Using Trained TBDT
"""
if plot_g:
    # Load trained TBDT
    regressor = load(case.result_paths[time] + estimator_name + '.joblib')
    # Predict g
    g_test = regressor.predict(x_test)
    g_train = regressor.predict(x_train)
    # Interpolate to mesh
    _, _, _, g_test_mesh = interpolateGridData(ccx_test, ccy_test, g_test,
                                                                              mesh_target=uniform_mesh_size,
                                                                              interp=interp_method)
    _, _, _, g_train_mesh = interpolateGridData(ccx_train, ccy_train, g_train,
                                                           mesh_target=uniform_mesh_size, interp=interp_method)
    # Make an individual directory
    gdir = case.result_paths[time] + '/g'
    os.makedirs(gdir, exist_ok=True)
    extent_test = (ccx_test.min(), ccx_test.max(), ccy_test.min(), ccy_test.max())
    extent_train = (ccx_train.min(), ccx_train.max(), ccy_train.min(), ccy_train.max())
    # glim = (-10000, 10000)
    glim = None
    # g diff = |g_train - g_test|
    g_diff_mesh = np.abs(g_train_mesh - g_test_mesh)
    # Go through every basis
    for i in range(g_test.shape[1]):
        zlabel = "$|g'^{}|$ [-]".format('{(' + str(i + 1) + ')}')
        cbticks = [0, 100, 10000, 1e6, 1e8, 1e10]
        # g diff = |g_train - g_test|
        figname = 'g' + str(i + 1) + '_diff'
        gplot_diff = Plot2D_Image(val=g_diff_mesh[:, i],
                                  name=figname,
                                  xlabel=xlabel, ylabel=ylabel, val_label=zlabel,
                                  save=save_fig,
                                  show=show,
                                  figdir=gdir,
                                  rotate_img=True,
                                  extent=extent_test,
                                  figwidth='1/3',
                                  val_lim=glim,
                                  figheight_multiplier=figheight_multiplier)
        gplot_diff.initializeFigure()
        gplot_diff.plotFigure(norm='symlog')
        gplot_diff.axes.add_patch(next(patches))
        gplot_diff.finalizeFigure(cbticks=cbticks)

        # # First test g
        # zlabel = '$g^{}$'.format('{(' + str(i + 1) + ')}')
        # # First test g
        # figname = 'g' + str(i + 1) + '_test'
        # gplot_test = Plot2D_Image(val=g_test_mesh[:, i],
        #                           name=figname,
        #                           xLabel=xlabel, yLabel=ylabel, zLabel=zlabel,
        #                           save=save_fig,
        #                           show=show,
        #                           figDir=gdir,
        #                           rotate_img=True,
        #                           extent=extent_test,
        #                           figWidth='1/3',
        #                           zlim=glim)
        # gplot_test.initializeFigure()
        # gplot_test.plotFigure(norm='symlog')
        # gplot_test.axes.add_patch(next(patches))
        # gplot_test.finalizeFigure()
        #
        # # Then train g
        # figname = 'g' + str(i + 1) + '_train'
        # gplot_train = Plot2D_Image(val=g_train_mesh[:, i],
        #                           name=figname,
        #                           xLabel=xlabel, yLabel=ylabel, zLabel=zlabel,
        #                           save=save_fig,
        #                           show=show,
        #                           figDir=gdir,
        #                           rotate_img=True,
        #                           extent=extent_test,
        #                           figWidth='1/3',
        #                           zlim=glim)
        # gplot_train.initializeFigure()
        # gplot_train.plotFigure(norm='symlog')
        # gplot_train.axes.add_patch(next(patches))
        # gplot_train.finalizeFigure()

if plot_x:
    xdir = case.result_paths[time] + '/X'
    os.makedirs(xdir, exist_ok=True)
    # Interpolate to mesh
    _, _, _, x_test_mesh = interpolateGridData(ccx_test, ccy_test, x_test,
                                               mesh_target=uniform_mesh_size,
                                               interp=interp_method)
    _, _, _, x_train_mesh = interpolateGridData(ccx_train, ccy_train, x_train,
                                                mesh_target=uniform_mesh_size, interp=interp_method)

    if fs == 'grad(TKE)_grad(p)':
        labels = ['S^2', 'S^3', '\Omega^2', '\Omega^2S', '\Omega^2S^2', '\Omega^2S\Omega S^2',
                  'A_k^2', 'A_k^2S', 'A_k^2S^2', 'A_k^2SA_kS^2', '\Omega A_k', '\Omega A_kS', '\Omega A_kS^2',
                  '\Omega^2A_kS', 'A_k^2\Omega S', '\Omega^2A_kS^2', 'A_k^2\Omega S^2', '\Omega^2SA_kS^2', 'A_k^2S\Omega S^2',
                  'A_p^2', 'A_p^2S', 'A_p^2S^2', 'A_p^2SA_pS^2', '\Omega A_p', '\Omega A_pS', '\Omega A_pS^2',
                  '\Omega^2A_pS', 'A_p^2\Omega S', '\Omega^2A_pS^2', 'A_p^2\Omega S^2', '\Omega^2SA_pS^2', 'A_p^2S\Omega S^2',
                  'A_kA_p', 'A_kA_pS', 'A_kA_pS^2', 'A_k^2A_pS', 'A_p^2A_kS', 'A_k^2A_pS^2', 'A_p^2A_kS^2',
                  'A_k^2SA_pS^2', 'A_p^2SA_kS^2',
                  '\Omega A_kA_p', '\Omega A_kA_pS', '\Omega A_pA_kS', '\Omega A_kA_pS^2', '\Omega A_pA_kS^2', '\Omega A_kSA_pS^2']
        zlabel = ['1.1', '1.2', '1.3', '1.4', '1.5', '1.6',
                  '2.1', '2.2', '2.3', '2.4', '2.5', '2.6', '2.7', '2.8', '2.9', '2.10', '2.11', '2.12', '2.13',
                  '3.1', '3.2', '3.3', '3.4', '3.5', '3.6', '3.7', '3.8', '3.9', '3.10', '3.11', '3.12', '3.13',
                  '4.1', '4.2', '4.3', '4.4', '4.5', '4.6', '4.7', '4.8', '4.9', '4.10', '4.11', '4.12', '4.13', '4.14', '4.15']
    elif 'grad(TKE)' in fs:
        labels = ['S^2', 'S^3', '\Omega^2', '\Omega^2S', '\Omega^2S^2', '\Omega^2S\Omega S^2',
                  'A_k^2', 'A_k^2S', 'A_k^2S^2', 'A_k^2SA_kS^2', '\Omega A_k', '\Omega A_kS', '\Omega A_kS^2',
                  '\Omega^2A_kS', 'A_k^2\Omega S', '\Omega^2A_kS^2', 'A_k^2\Omega S^2', '\Omega^2SA_kS^2', 'A_k^2S\Omega S^2']
        zlabel = ['1.1', '1.2', '1.3', '1.4', '1.5', '1.6',
                  '2.1', '2.2', '2.3', '2.4', '2.5', '2.6', '2.7', '2.8', '2.9', '2.10', '2.11', '2.12', '2.13']
    else:
        labels = ['S^2', 'S^3', '\Omega^2', '\Omega^2S', '\Omega^2S^2', '\Omega^2S\Omega S^2']
        zlabel = ['1.1', '1.2', '1.3', '1.4', '1.5', '1.6']

    for i in range(len(zlabel)):
        zlabel[i] = 'Feature ' + zlabel[i] + ' [-]'
        labels[i] = '$' + labels[i] + '$ [-]'

    for i in range(x_test.shape[1]):
        # For test features
        figname = 'X' + str(i + 1) + 'test'
        xlim = min(min(x_train[:, i]), min(x_test[:, i])), max(max(x_train[:, i]), max(x_test[:, i]))
        xtest_plot = Plot2D_Image(val=x_test_mesh[..., i],
                                  name=figname,
                                  xlabel=xlabel, ylabel=ylabel, val_label=labels[i],
                                  save=save_fig,
                                  show=show,
                                  figdir=xdir,
                                  rotate_img=True,
                                  extent=extent_test,
                                  figwidth='1/3',
                                  val_lim=xlim,
                                  figheight_multiplier=figheight_multiplier)
        xtest_plot.initializeFigure()
        xtest_plot.plotFigure()
        xtest_plot.axes.add_patch(next(patches))
        xtest_plot.finalizeFigure()

        # For train features
        figname = 'X' + str(i + 1) + 'train'
        xtrain_plot = Plot2D_Image(val=x_train_mesh[..., i],
                                  name=figname,
                                  xlabel=xlabel, ylabel=ylabel, val_label=labels[i],
                                  save=save_fig,
                                  show=show,
                                  figdir=xdir,
                                  rotate_img=True,
                                  extent=extent_train,
                                  figwidth='1/3',
                                  val_lim=xlim,
                                  figheight_multiplier=figheight_multiplier)
        xtrain_plot.initializeFigure()
        xtrain_plot.plotFigure()
        xtrain_plot.axes.add_patch(next(patches))
        xtrain_plot.finalizeFigure()

# # First, take Frobenius norm on Tij to collapse n_output axis
# # so that Tij becomes shape (n_samples, n_bases)
# # 4 times the mesh grid, ignoring b13 and b23
# x_test_contour, y_test_contour = (ccx_test_mesh,)*4, (ccy_test_mesh,)*4
# x_train_contour, y_train_contour = (ccx_train_mesh,)*4, (ccy_train_mesh,)*4
#
# # Go through each basis
# for i in range(1):#tb_test.shape[2]):
#     # Same with interpolated Tij of shape (nx, ny, n_bases, n_outputs), ignore b13, b23
#     vals_test = (tb_test_mesh[..., i, 0], tb_test_mesh[..., i, 1], tb_test_mesh[..., i, 3], tb_test_mesh[..., i, 5])
#     vals_train = (tb_train_mesh[..., i, 0], tb_train_mesh[..., i, 1], tb_train_mesh[..., i, 3], tb_train_mesh[..., i, 5])
#
#     figname = 'Tij_' + str(i + 1)
#     tbplot_test = PlotContourSlices3D(x_test_contour, y_test_contour, vals_test, np.arange(4), contourLvl=contour_lvl, gradientBg=False, equalAxis=True, zDir='x',
#                                       name=figname, xLabel=xlabel, yLabel=ylabel, save=save_fig, show=show,
#                                       figDir=case.result_paths[time],
#                                       viewAngles=(45, 90),
#                                       cbarOrientate='vertical',
#                                       zLabel=None)
#     tbplot_test.initializeFigure(figSize=(2, 1))
#     tbplot_test.plotFigure()
#     for j in range(4):
#         patch_i = next(patches)
#         tbplot_test.axes.add_patch(patch_i)
#         art3d.pathpatch_2d_to_3d(patch_i, z=j, zdir='x')
#
#     # plt.show()
#     tbplot_test.finalizeFigure()
#     # plt.savefig(tbplot_test.figDir + '/' + tbplot_test.name + '.png',
#     #             dpi=1000)
#     print('\nFigure ' + tbplot_test.name + '.png saved in ' + tbplot_test.figDir)
#
#     # tbplot_train = PlotContourSlices3D(x_train_contour, y_train_contour, vals_train, np.arange(4), contourLvl=contour_lvl, gradientBg=False, equalAxis=True,
#     #                              name=figname, xLabel=xlabel, yLabel=ylabel,save=save_fig, show=show, figDir=case.result_paths[time])
#     # tbplot_train.initializeFigure()
#     # tbplot_train.plotFigure()
#     # tbplot_train.axes.add_patch(next(patches))
#     # tbplot_train.finalizeFigure()

max_stars_repo_name: YuyangL/TurbulenceMachineLearning

max_issues_repo_name: YuyangL/TurbulenceMachineLearning

max_forks_repo_name: YuyangL/TurbulenceMachineLearning

reason: import numpy,from scipy