text:
#!/usr/bin/env python
##############################################################################
#
# diffpy.structure  by DANSE Diffraction group
#                   Simon J. L. Billinge
#                   (c) 2007 trustees of the Michigan State University.
#                   All rights reserved.
#
# File coded by:    Pavol Juhas
#
# See AUTHORS.txt for a list of people who contributed.
# See LICENSE_DANSE.txt for license information.
#
##############################################################################

"""This module defines class Structure.
"""

import copy as copymod
import numpy
import codecs
import six

from diffpy.structure.lattice import Lattice
from diffpy.structure.atom import Atom
from diffpy.structure.utils import _linkAtomAttribute, atomBareSymbol
from diffpy.structure.utils import isiterable

# ----------------------------------------------------------------------------

class Structure(list):
    """Structure --> group of atoms

    Structure class is inherited from Python list.  It contains
    a list of Atom instances.  Structure overloads setitem and setslice
    methods so that the lattice attribute of atoms get set to lattice.

    Data members:
        title   -- structure description
        lattice -- coordinate system (instance of Lattice)
        pdffit  -- None or a dictionary of PDFFit-related metadata
    """

    # default values for instance attributes

    title = ''
    _lattice = None
    pdffit = None

    def __init__(self, atoms=None, lattice=None, title=None,
            filename=None, format=None):
        """define group of atoms in a specified lattice.

        atoms    -- list of Atom instances to be included in this Structure.
                    When atoms argument is an existing Structure instance,
                    the new Structure is its copy.
        lattice  -- instance of Lattice defining coordinate systems, property.
        title    -- string description of the structure
        filename -- optional, name of a file to load the structure from.
                    Overrides atoms argument when specified.
        format   -- optional structure format of the loaded filename.  By default
                    all structure formats are tried one by one.  Ignored when
                    filename has not been specified.

        Structure(stru)     create a copy of Structure instance stru.

        Because Structure is inherited from a list it can use list expansions,
        for example:
            oxygen_atoms = [ for a in stru if a.element == "O" ]
            oxygen_stru = Structure(oxygen_atoms, lattice=stru.lattice)
        """
        # if filename is specified load it and return
        if filename is not None:
            if any((atoms, lattice, title)):
                emsg = "Cannot use filename and atoms arguments together."
                raise ValueError(emsg)
            readkwargs = (format is not None) and {'format' : format} or {}
            self.read(filename, **readkwargs)
            return
        # copy initialization, must be first to allow lattice, title override
        if isinstance(atoms, Structure):
            Structure.__copy__(atoms, self)
        # assign arguments:
        if title is not None:
            self.title = title
        if lattice is not None:
            self.lattice = lattice
        elif self.lattice is None:
            self.lattice = Lattice()
        # insert atoms unless already done by __copy__
        if not len(self) and atoms is not None:
            self.extend(atoms)
        return

    def copy(self):
        '''Return a copy of this Structure object.
        '''
        return copymod.copy(self)

    def __copy__(self, target=None):
        '''Create a deep copy of this instance.

        target   -- optional target instance for copying, useful for
                    copying a derived class.  Defaults to new instance
                    of the same type as self.

        Return a duplicate instance of this object.
        '''
        if target is None:
            target = Structure()
        elif target is self:
            return target
        # copy attributes as appropriate:
        target.title = self.title
        target.lattice = Lattice(self.lattice)
        target.pdffit = copymod.deepcopy(self.pdffit)
        # copy all atoms to the target
        target[:] = self
        return target

    def __str__(self):
        """simple string representation"""
        s_lattice = "lattice=%s" % self.lattice
        s_atoms = '\n'.join([str(a) for a in self])
        return s_lattice + '\n' + s_atoms

    def addNewAtom(self, *args, **kwargs):
        """Add new Atom instance to the end of this Structure.

        All arguments are forwarded to Atom constructor.

        No return value.
        """
        kwargs['lattice'] = self.lattice
        a = Atom(*args, **kwargs)
        self.append(a, copy=False)
        return

    def getLastAtom(self):
        """Return Reference to the last Atom in this structure.
        """
        last_atom = self[-1]
        return last_atom

    def assignUniqueLabels(self):
        """Set a unique label string for each atom in this structure.
        The label strings are formatted as "%(baresymbol)s%(index)i",
        where baresymbol is the element right-stripped of "[0-9][+-]".

        No return value.
        """
        elnum = {}
        # support duplicate atom instances
        islabeled = set()
        for a in self:
            if a in islabeled:  continue
            baresmbl = atomBareSymbol(a.element)
            elnum[baresmbl] = elnum.get(baresmbl, 0) + 1
            a.label = baresmbl + str(elnum[baresmbl])
            islabeled.add(a)
        return

    def distance(self, aid0, aid1):
        """Distance between 2 atoms, no periodic boundary conditions.

        aid0 -- zero based index of the first atom or a string label
                such as "Na1"
        aid1 -- zero based index or string label of the second atom.

        Return float.
        Raise IndexError for invalid arguments.
        """
        # lookup by labels
        a0, a1 = self[aid0, aid1]
        return self.lattice.dist(a0.xyz, a1.xyz)

    def angle(self, aid0, aid1, aid2):
        """The bond angle at the second of three atoms in degrees.

        aid0 -- zero based index of the first atom or a string label
                such as "Na1"
        aid1 -- index or string label for the second atom,
                where the angle is formed
        aid2 -- index or string label for the third atom

        Return float.
        Raise IndexError for invalid arguments.
        """
        a0, a1, a2 = self[aid0, aid1, aid2]
        u10 = a0.xyz - a1.xyz
        u12 = a2.xyz - a1.xyz
        return self.lattice.angle(u10, u12)

    def placeInLattice(self, new_lattice):
        """place structure into new_lattice coordinate system

        sets lattice to new_lattice and recalculate fractional coordinates
        of all atoms so their absolute positions remain the same

        return self
        """
        Tx = numpy.dot(self.lattice.base, new_lattice.recbase)
        Tu = numpy.dot(self.lattice.normbase, new_lattice.recnormbase)
        for a in self:
            a.xyz = numpy.dot(a.xyz, Tx)
            if a.anisotropy:
                a.U = numpy.dot(numpy.transpose(Tu), numpy.dot(a.U, Tu))
        self.lattice = new_lattice
        return self

    def read(self, filename, format='auto'):
        """Load structure from a file, any original data become lost.

        filename -- file to be loaded
        format   -- all structure formats are defined in parsers submodule,
                    when format == 'auto' all parsers are tried one by one

        Return instance of data Parser used to process file.  This
        can be inspected for information related to particular format.
        """
        import diffpy.structure
        import diffpy.structure.parsers
        getParser = diffpy.structure.parsers.getParser
        p = getParser(format)
        new_structure = p.parseFile(filename)
        # reinitialize data after successful parsing
        # avoid calling __init__ from a derived class
        Structure.__init__(self)
        if new_structure is not None:
            self.__dict__.update(new_structure.__dict__)
            self[:] = new_structure
        if not self.title:
            import os.path
            tailname = os.path.basename(filename)
            tailbase = os.path.splitext(tailname)[0]
            self.title = tailbase
        return p

    def readStr(self, s, format='auto'):
        """Load structure from a string, any original data become lost.

        s        -- string with structure definition
        format   -- all structure formats are defined in parsers submodule,
                    when format == 'auto' all parsers are tried one by one

        Return instance of data Parser used to process input string.  This
        can be inspected for information related to particular format.
        """
        from diffpy.structure.parsers import getParser
        p = getParser(format)
        new_structure = p.parse(s)
        # reinitialize data after successful parsing
        # avoid calling __init__ from a derived class
        Structure.__init__(self)
        if new_structure is not None:
            self.__dict__.update(new_structure.__dict__)
            self[:] = new_structure
        return p

    def write(self, filename, format):
        """Save structure to file in the specified format

        No return value.

        Note: available structure formats can be obtained by:
            from parsers import formats
        """
        from diffpy.structure.parsers import getParser
        p = getParser(format)
        p.filename = filename
        s = p.tostring(self)
        with codecs.open(filename, 'w', encoding='UTF-8') as fp:
            fp.write(s)
        return

    def writeStr(self, format):
        """return string representation of the structure in specified format

        Note: available structure formats can be obtained by:
            from parsers import formats
        """
        from diffpy.structure.parsers import getParser
        p = getParser(format)
        s = p.tostring(self)
        return s

    def tolist(self):
        '''Return atoms in this Structure as a standard Python list.
        '''
        rv = [a for a in self]
        return rv

    # Overloaded list Methods and Operators ----------------------------------

    def append(self, a, copy=True):
        """Append atom to a structure and update its lattice attribute.

        a    -- instance of Atom
        copy -- flag for appending a copy of a.
                When False, append a and update a.lattice.

        No return value.
        """
        adup = copy and Atom(a) or a
        adup.lattice = self.lattice
        super(Structure, self).append(adup)
        return

    def insert(self, idx, a, copy=True):
        """Insert atom a before position idx in this Structure.

        idx  -- position in atom list
        a    -- instance of Atom
        copy -- flag for inserting a copy of a.
                When False, append a and update a.lattice.

        No return value.
        """
        adup = copy and copymod.copy(a) or a
        adup.lattice = self.lattice
        super(Structure, self).insert(idx, adup)
        return

    def extend(self, atoms, copy=None):
        """Extend Structure with an iterable of atoms.

        Update the `lattice` attribute of all added atoms.

        Parameters
        ----------
        atoms : iterable
            The `Atom` objects to be appended to this Structure.
        copy : bool, optional
            Flag for adding copies of Atom objects.
            Make copies when `True`, append `atoms` unchanged when ``False``.
            The default behavior is to make copies when `atoms` are of
            `Structure` type or if new atoms introduce repeated objects.
        """
        adups = (copymod.copy(a) for a in atoms)
        if copy is None:
            if isinstance(atoms, Structure):
                newatoms = adups
            else:
                memo = set(id(a) for a in self)
                nextatom = lambda a: (a if id(a) not in memo
                                      else copymod.copy(a))
                mark = lambda a: (memo.add(id(a)), a)[-1]
                newatoms = (mark(nextatom(a)) for a in atoms)
        elif copy:
            newatoms = adups
        else:
            newatoms = atoms
        setlat = lambda a: (setattr(a, 'lattice', self.lattice), a)[-1]
        super(Structure, self).extend(setlat(a) for a in newatoms)
        return

    def __getitem__(self, idx):
        """Get one or more atoms in this structure.

        idx  -- atom identifier, which can be integer, string or iterable.
                When integer use standard list lookup.  For iterables use
                numpy lookup, this supports integer or boolean flag arrays.
                For string or string-containing iterables lookup the atoms
                by string label.

        Return an Atom instance for integer or string index or a substructure
        in all other cases.  Raise IndexError for invalid index or for
        non-unique atom label.

        Examples:

        stru[0]  -->  first atom in the Structure
        stru[stru.element == 'Na']  -->  substructure of all Na atoms
        stru['Na3']  -->  atom with a unique label 'Na3'
        stru['Na3', 2, 'Cl2']  -->  substructure of three atoms, lookup by
            label is more efficient when done for several atoms at once.
        """
        if isinstance(idx, slice):
            rv = self.__emptySharedStructure()
            lst = super(Structure, self).__getitem__(idx)
            rv.extend(lst, copy=False)
            return rv
        try:
            rv = super(Structure, self).__getitem__(idx)
            return rv
        except TypeError:
            pass
        # check if there is any string label that should be resolved
        scalarstringlabel = isinstance(idx, six.string_types)
        hasstringlabel = scalarstringlabel or (
            isiterable(idx) and
            any(isinstance(ii, six.string_types) for ii in idx))
        # if not, use numpy indexing to resolve idx
        if not hasstringlabel:
            idx1 = idx
            if type(idx) is tuple:
                idx1 = numpy.r_[idx]
            indices = numpy.arange(len(self))[idx1]
            rhs = [list.__getitem__(self, i) for i in indices]
            rv = self.__emptySharedStructure()
            rv.extend(rhs, copy=False)
            return rv
        # here we need to resolve at least one string label
        # build a map of labels to indices and mark duplicate labels
        duplicate = object()
        labeltoindex = {}
        for i, a in enumerate(self):
            labeltoindex[a.label] = (
                duplicate if a.label in labeltoindex else i)
        def _resolveindex(aid):
            aid1 = aid
            if type(aid) is str:
                aid1 = labeltoindex.get(aid, None)
                if aid1 is None:
                    raise IndexError("Invalid atom label %r." % aid)
                if aid1 is duplicate:
                    raise IndexError("Atom label %r is not unique." % aid)
            return aid1
        # generate new index object that has no strings
        if scalarstringlabel:
            idx2 = _resolveindex(idx)
        # for iterables preserve the tuple object type
        else:
            idx2 = [_resolveindex(i) for i in idx]
            if type(idx) is tuple:
                idx2 = tuple(idx2)
        # call this function again and hope there is no recursion loop
        rv = self[idx2]
        return rv

    def __setitem__(self, idx, value, copy=True):
        """Assign self[idx] atom to value.

        idx  -- index of atom in this Structure or a slice
        value -- instance of Atom or an iterable.
        copy -- flag for making a copy of the value.  When False, update
                the `lattice` attribute of Atom objects present in value.

        No return value.
        """
        # handle slice assignment
        if isinstance(idx, slice):
            def _fixlat(a):
                a.lattice = self.lattice
                return a
            v1 = value
            if copy:
                keep = set(super(Structure, self).__getitem__(idx))
                v1 = (a if a in keep else Atom(a) for a in value)
            vfinal = filter(_fixlat, v1)
        # handle scalar assingment
        else:
            vfinal = Atom(value) if copy else value
            vfinal.lattice = self.lattice
        super(Structure, self).__setitem__(idx, vfinal)
        return

    def __add__(self, other):
        '''Return new Structure object with appended atoms from other.

        other    -- sequence of Atom instances

        Return new Structure with a copy of Atom instances.
        '''
        rv = copymod.copy(self)
        rv += other
        return rv

    def __iadd__(self, other):
        '''Extend this Structure with atoms from other.

        other    -- sequence of Atom instances

        Return self.
        '''
        self.extend(other, copy=True)
        return self

    def __sub__(self, other):
        '''Return new Structure that has atoms from the other removed.

        other    -- sequence of Atom instances

        Return new Structure with a copy of Atom instances.
        '''
        otherset = set(other)
        keepindices = [i for i, a in enumerate(self) if not a in otherset]
        rv = copymod.copy(self[keepindices])
        return rv

    def __isub__(self, other):
        '''Remove other atoms if present in this structure.

        other    -- sequence of Atom instances

        Return self.
        '''
        otherset = set(other)
        self[:] = [a for a in self if a not in otherset]
        return self

    def __mul__(self, n):
        '''Return new Structure with n-times concatenated atoms from self.
        Atoms and lattice in the new structure are all copies.

        n    -- integer multiple

        Return new Structure.
        '''
        rv = copymod.copy(self[:0])
        rv += n * self.tolist()
        return rv

    # right-side multiplication is the same as left-side
    __rmul__ = __mul__

    def __imul__(self, n):
        '''Concatenate this Structure to n-times more atoms.
        For positive multiple the current Atom objects remain at the
        beginning of this Structure.

        n    -- integer multiple

        Return self.
        '''
        if n <= 0:
            self[:] = []
        else:
            self.extend((n - 1) * self.tolist(), copy=True)
        return self

    # Properties -------------------------------------------------------------

    # lattice

    def _get_lattice(self):
        return self._lattice

    def _set_lattice(self, value):
        for a in self:  a.lattice = value
        self._lattice = value
        return

    lattice = property(_get_lattice, _set_lattice, doc =
        "Coordinate system for this Structure.")

    # composition

    def _get_composition(self):
        rv = {}
        for a in self:
            rv[a.element] = rv.get(a.element, 0.0) + a.occupancy
        return rv

    composition = property(_get_composition,
        doc="Dictionary of chemical symbols and their total occupancies.")

    # linked atom attributes

    element = _linkAtomAttribute('element',
        '''Character array of atom types.  Assignment updates
        the element attribute of the respective atoms.''',
        toarray=numpy.char.array)

    xyz = _linkAtomAttribute('xyz',
        '''Array of fractional coordinates of all atoms.
        Assignment updates xyz attribute of all atoms.''')

    x = _linkAtomAttribute('x',
        '''Array of all fractional coordinates x.
        Assignment updates xyz attribute of all atoms.''')

    y = _linkAtomAttribute('y',
        '''Array of all fractional coordinates y.
        Assignment updates xyz attribute of all atoms.''')

    z = _linkAtomAttribute('z',
        '''Array of all fractional coordinates z.
        Assignment updates xyz attribute of all atoms.''')

    label = _linkAtomAttribute('label',
        '''Character array of atom names.  Assignment updates
        the label attribute of all atoms.''',
        toarray=numpy.char.array)

    occupancy = _linkAtomAttribute('occupancy',
        '''Array of atom occupancies.  Assignment updates the
        occupancy attribute of all atoms.''')

    xyz_cartn = _linkAtomAttribute('xyz_cartn',
        '''Array of absolute Cartesian coordinates of all atoms.
        Assignment updates the xyz attribute of all atoms.''')

    anisotropy = _linkAtomAttribute('anisotropy',
        '''Boolean array for anisotropic thermal displacement flags.
        Assignment updates the anisotropy attribute of all atoms.''')

    U = _linkAtomAttribute('U',
        '''Array of anisotropic thermal displacement tensors.
        Assignment updates the U and anisotropy attributes of all atoms.''')

    Uisoequiv = _linkAtomAttribute('Uisoequiv',
        '''Array of isotropic thermal displacement or equivalent values.
        Assignment updates the U attribute of all atoms.''')

    U11 = _linkAtomAttribute('U11',
        '''Array of U11 elements of the anisotropic displacement tensors.
        Assignment updates the U and anisotropy attributes of all atoms.''')

    U22 = _linkAtomAttribute('U22',
        '''Array of U22 elements of the anisotropic displacement tensors.
        Assignment updates the U and anisotropy attributes of all atoms.''')

    U33 = _linkAtomAttribute('U33',
        '''Array of U33 elements of the anisotropic displacement tensors.
        Assignment updates the U and anisotropy attributes of all atoms.''')

    U12 = _linkAtomAttribute('U12',
        '''Array of U12 elements of the anisotropic displacement tensors.
        Assignment updates the U and anisotropy attributes of all atoms.''')

    U13 = _linkAtomAttribute('U13',
        '''Array of U13 elements of the anisotropic displacement tensors.
        Assignment updates the U and anisotropy attributes of all atoms.''')

    U23 = _linkAtomAttribute('U23',
        '''Array of U23 elements of the anisotropic displacement tensors.
        Assignment updates the U and anisotropy attributes of all atoms.''')

    Bisoequiv = _linkAtomAttribute('Bisoequiv',
        '''Array of Debye-Waller isotropic thermal displacement or equivalent
        values.  Assignment updates the U attribute of all atoms.''')

    B11 = _linkAtomAttribute('B11',
        '''Array of B11 elements of the Debye-Waller displacement tensors.
        Assignment updates the U and anisotropy attributes of all atoms.''')

    B22 = _linkAtomAttribute('B22',
        '''Array of B22 elements of the Debye-Waller displacement tensors.
        Assignment updates the U and anisotropy attributes of all atoms.''')

    B33 = _linkAtomAttribute('B33',
        '''Array of B33 elements of the Debye-Waller displacement tensors.
        Assignment updates the U and anisotropy attributes of all atoms.''')

    B12 = _linkAtomAttribute('B12',
        '''Array of B12 elements of the Debye-Waller displacement tensors.
        Assignment updates the U and anisotropy attributes of all atoms.''')

    B13 = _linkAtomAttribute('B13',
        '''Array of B13 elements of the Debye-Waller displacement tensors.
        Assignment updates the U and anisotropy attributes of all atoms.''')

    B23 = _linkAtomAttribute('B23',
        '''Array of B23 elements of the Debye-Waller displacement tensors.
        Assignment updates the U and anisotropy attributes of all atoms.''')

    # Private Methods --------------------------------------------------------

    def __emptySharedStructure(self):
        '''Return empty Structure with standard attributes same as in self.
        '''
        rv = Structure()
        rv.__dict__.update([(k, getattr(self, k)) for k in rv.__dict__])
        return rv

    # Python 2 Compatibility -------------------------------------------------

    if six.PY2:

        def __getslice__(self, lo, hi):
            return self.__getitem__(slice(lo, hi))

        def __setslice__(self, lo, hi, sequence):
            self.__setitem__(slice(lo, hi), sequence)
            return

    # ------------------------------------------------------------------------

# End of class Structure


text:
#!/usr/bin/python

# Copyright (c) 2017 Ivan Iakoupov
#
# Permission is hereby granted, free of charge, to any person obtaining a
# copy of this software and associated documentation files (the "Software"),
# to deal in the Software without restriction, including without limitation
# the rights to use, copy, modify, merge, publish, distribute, sublicense,
# and/or sell copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included
# in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
# DEALINGS IN THE SOFTWARE.

import os.path

import matplotlib.pyplot as p
import numpy as np
from scipy.optimize import minimize
from scipy.integrate import quad
from matplotlib_util import extract_params_from_file_name,\
                            read_column_names
from common_typeset_info import usetex, common_line_width

#Use matplotlib pre-2.0 version style
p.style.use('classic')

prefix = 'plot_data'

data_file_list = [
    'cphase_fidelity_OmegaScattering_and_OmegaStorageRetrieval_aE_a_opt_lambda_sagnac_g1d_0.05_N_10000_kd_0.5.txt',
    'cphase_fidelity_OmegaScattering_and_OmegaStorageRetrieval_aE_a_opt_dualv_sym_sagnac_g1d_0.05_N_10000_kd_0.266.txt',
]

usetex()
p.figure(figsize=(3.3*2,2.5))
data = []
column_dic = []
for file_name in data_file_list:
    param_dict = extract_params_from_file_name(file_name)
    full_path = os.path.join(prefix, file_name)
    if not os.path.exists(full_path):
        print('Path {} doesn\'t exist'.format(full_path))
        continue
    data.append(np.loadtxt(full_path, dtype=np.float64, delimiter=';',
                      unpack=True, skiprows=1))
    column_names = read_column_names(full_path)
    column_dic.append(dict(zip(column_names, range(len(column_names)))))

xdata_str = 'OmegaScattering'

p.subplot(1,2,1)
p.plot(data[0][column_dic[0][xdata_str]],
       data[0][column_dic[0]['F_CJ_conditional_tNoInteraction_one']],
       color='#000099', linestyle='-',
       label=r'$F_{\rm CJ,cond}$, $\Lambda$-type',
       linewidth=common_line_width)
p.plot(data[1][column_dic[1][xdata_str]],
       data[1][column_dic[1]['F_CJ_conditional_tNoInteraction_one']],
       color='#990000', linestyle='--',
       label=r'$F_{\rm CJ,cond}$, dual-V',
       linewidth=common_line_width)
p.plot(data[0][column_dic[0][xdata_str]],
       data[0][column_dic[0]['F_CJ_tNoInteraction_one']],
       color='#009900', linestyle=':',
       label=r'$F_{\rm CJ}$, $\Lambda$-type',
       linewidth=common_line_width)
p.plot(data[1][column_dic[1][xdata_str]],
       data[1][column_dic[1]['F_CJ_tNoInteraction_one']],
       'k-.',
       label=r'$F_{\rm CJ}$, dual-V',
       linewidth=common_line_width)
p.xlabel(r"$\Omega_0/\Gamma$")
p.ylim(0.4,1)
p.title(r"(a)")
p.legend(loc='center left')

p.subplot(1,2,2)
p.plot(data[0][column_dic[0][xdata_str]],
       -data[0][column_dic[0]['Deltac']],
       color='#000099', linestyle='-',
       label=r'$-\Delta_{\rm c}/\Gamma$, $\Lambda$-type',
       linewidth=common_line_width)
p.plot(data[1][column_dic[1][xdata_str]],
       -data[1][column_dic[1]['Deltac']],
       color='#990000', linestyle='--',
       label=r'$-\Delta_{\rm c}/\Gamma$, dual-V',
       linewidth=common_line_width)
p.plot(data[0][column_dic[0][xdata_str]],
       data[0][column_dic[0]['delta']],
       color='#009900', linestyle=':',
       label=r'$\delta_{\rm res}/\Gamma$, $\Lambda$-type',
       linewidth=common_line_width)
p.plot(data[1][column_dic[1][xdata_str]],
       data[1][column_dic[1]['delta']],
       'k-.',
       label=r'$\delta_{\rm res}/\Gamma$, dual-V',
       linewidth=common_line_width)
p.xlabel(r"$\Omega_0/\Gamma$")
p.title(r"(b)")
p.legend(loc='upper left')
p.tight_layout(pad=0.1)

p.savefig('figS09.eps')

ext: py

reason: import numpy,from scipy


text:
import matplotlib
import matplotlib.pyplot as plt
import numpy
import unittest
from unittest.mock import patch, Mock

import CrystFEL_Jupyter_utilities.histogram as histogram

class TestHistogram(unittest.TestCase):

    @patch('matplotlib.cbook.silent_list')
    @patch('matplotlib.pyplot')
    @patch('matplotlib.axes.Axes')
    def setUp(self, mock_ax, mock_plt, mock_patch):
        self.fig = mock_plt.figure()
        self.fig.add_subplot(1, 1, 1)
        self.mock_ax = mock_ax()
        self.patch = mock_patch()
        self.mock_ax.get_xlim.return_value = "current_xlim"
        self.mock_ax.hist.return_value = [None, None, self.patch]
        self.data = {'P': [1, 2, 3, 4, 5, 10], 'B': [2, 2, 3, 4, 4, 4],
                     'I': [6, 6, 5, 4, 5, 3], 'A': [1, 1, 1, 1]}
        self.dict_color = {'P': 'gray', 'A': 'cyan', 'B': 'darkblue',
                           'C': 'royalblue', 'I': 'lime',
                           'F': "magenta", 'H': "firebrick", 'R': 'olive'}
        self.title = "test_hist"
        self.xlabel = "test_label"
        self.hist = histogram.Histogram(axs=self.mock_ax, xlabel=self.xlabel,
                                        data_to_histogram=self.data,
                                        colors=self.dict_color,
                                        bins=10, name=self.title)

    def tearDown(self):
        self.fig.reset_mock()
        self.mock_ax.reset_mock()
        self.patch.reset_mock()

    def test_plt(self):
        assert self.fig.add_subplot.called
        self.fig.add_subplot.assert_called_with(1, 1, 1)
        self.fig.reset_mock()

    def test_hist(self):
        assert self.mock_ax.hist.called
        self.mock_ax.hist.assert_called_with(
            alpha=0.9, bins=10, density=1,
            color=['gray', 'cyan', 'darkblue', 'royalblue', 'lime',
                   'magenta', 'firebrick', 'olive', 'lightgray'],
            histtype='stepfilled', range=(1, 10),
            stacked=True, x=[[1, 2, 3, 4, 5, 10], [1, 1, 1, 1],
                             [2, 2, 3, 4, 4, 4], [], [6, 6, 5, 4, 5, 3],
                             [], [], [], []])
        self.assertEqual(self.patch, self.hist.patches)
        assert self.mock_ax.get_xlim.called

    def test_bool_crystal_exluded_green_space(self):
        self.assertEqual(self.hist.bool_crystal_exluded_green_space(1111),
                         False)
        self.hist.range_green_space = [0, 2]
        self.assertEqual(self.hist.bool_crystal_exluded_green_space(3), True)
        self.hist.range_green_space = [None, None]

    def test_set_range_green_space(self):
        self.assertEqual(self.hist.range_green_space, [None, None])
        self.hist.range_green_space = 1, 4
        self.assertEqual(self.hist.range_green_space, (1, 4))

    def test_draw_green_space(self):
        self.hist.range_green_space = 1, 3
        self.hist.draw_green_space()
        assert self.mock_ax.axvspan.called
        self.mock_ax.axvspan.assert_called_with(1, 3, facecolor='#2ca02c',
                                                alpha=0.5)

    def test_set_was_clicked_before(self):
        self.assertEqual(self.hist.was_clicked_before, False)
        self.hist.was_clicked_before = True
        self.assertEqual(self.hist.was_clicked_before, True)

    def test_get_was_clicked_before(self):
        self.hist.was_clicked_before = True
        self.assertEqual(self.hist.was_clicked_before, True)

    def test_set_bins(self):
        self.hist.bins = 15
        self.assertEqual(self.hist.bins, 15)
        self.hist.bins = 22
        self.assertEqual(self.hist.bins, 22)

    def test_set_name(self):
        self.hist.name = "New title"
        self.assertEqual(self.hist.name, "New title")

    def test_update(self):
        self.mock_ax.reset_mock()
        self.data = {'P': [2, 3, 4, 5], 'A': [3, 4, 5, 6, 6]}
        self.data_excluded = [1, 1, 2, 3, 2]
        self.hist.update(self.data, self.data_excluded)
        self.assertListEqual(self.hist.list_data[0], [2, 3, 4, 5])
        self.assertListEqual(self.hist.list_data[1], [3, 4, 5, 6, 6])
        self.assertListEqual(self.hist.list_data[-2], [])
        self.assertListEqual(self.hist.data_excluded, [1, 1, 2, 3, 2])
        assert self.mock_ax.clear.called
        assert self.mock_ax.set_title.called
        self.mock_ax.set_title.assert_called_with("Histogram of " +
                                                  self.hist.name)
        assert self.mock_ax.set_xlabel.called
        self.mock_ax.set_xlabel.assert_called_with(self.xlabel)
        assert self.mock_ax.hist.called

    def test_set_color(self):
        self.dict_color = {'P': 'magenta', 'A': 'magenta', 'B': 'magenta',
                           'C': 'olive', 'I': 'olive', 'F': "olive",
                           'H': "magenta", 'R': 'firebrick'}
        self.hist.list_colors = self.dict_color
        list_test = ['magenta', 'magenta', 'magenta', 'olive', 'olive',
                     'olive', 'magenta', 'firebrick', 'lightgray']
        self.assertListEqual(self.hist.list_colors, list_test)

    def test_current_xlim(self):
        self.assertEqual(self.hist.current_xlim, 'current_xlim')
        self.mock_ax.get_xlim.reset_mock()
        self.hist.update_current_xlim()
        assert self.mock_ax.get_xlim.called

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

reason: import numpy


text:
""" Calculate the shape of a net's output's based on
    its configuration files.
    author: Ioannis Athansiadis (supernlogn)
"""
import numpy as np

__all__ = ["_get_output_shape"]

def _get_output_shape(mc):
  """ Given the image shape from the model configuration,
      this function produces the final feature output shape
      @param mc model configuration dictionary
      @Note channel num is not always right
  """

  # # # shapes are (H, W, C) using the standard tensorflow paradigm
  if(mc.NET == "squeezeDet+"):
    ret_shape = get_squeezeDetPlus_shape(mc)
  elif(mc.NET == "squeezeDet"):
    ret_shape = get_squeezeDet_shape(mc)
  else:
    ret_shape = None

  return (int(ret_shape[0]), int(ret_shape[1]), int(ret_shape[2]))

def get_conv2d_out_shape(in_dims, filter_sizes, strides, padding):
  """ @param in_dims = input (height, width, channels)
      @param filter_sizes = (kernel height, kernel width, number of filters)
      @param strides = (filter height stride, filter width stride)
  """
  if(padding == 'SAME'):
    return (np.ceil(float(in_dims[0]) / float(strides[0])),
            np.ceil(float(in_dims[1]) / float(strides[1])),
            float(filter_sizes[2]))
  elif(padding == 'VALID'):
    return (np.ceil(float(in_dims[0] - filter_sizes[0] + 1) / float(strides[0])),
            np.ceil(float(in_dims[1] - filter_sizes[1] + 1) / float(strides[1])),
            float(filter_sizes[2]))

def get_pool_out_shape(in_dims, filter_sizes, strides, padding):
  if(padding == 'SAME'):
    return (np.ceil(float(in_dims[0]) / float(strides[0])),
            np.ceil(float(in_dims[1]) / float(strides[1])),
            in_dims[2])
  elif(padding == 'VALID'):
    return (np.ceil(float(in_dims[0] - filter_sizes[0] + 1) / float(strides[0])),
            np.ceil(float(in_dims[1] - filter_sizes[1] + 1) / float(strides[1])),
            in_dims[2])

def get_fire_out_shape(in_dims, s1x1, e1x1, e3x3):
  sq1x1_shape = get_conv2d_out_shape(in_dims, (1, 1, s1x1), [1, 1],'SAME')
  ex1x1_shape = get_conv2d_out_shape(sq1x1_shape, (1, 1, e1x1), [1, 1], 'SAME')
  ex3x3_shape = get_conv2d_out_shape(sq1x1_shape, (3, 3, e3x3), [1, 1], 'SAME')
  return (ex3x3_shape[0], ex3x3_shape[1], ex3x3_shape[2] + ex1x1_shape[2])

def get_squeezeDetPlus_shape(mc):
  input_shape = (mc.IMAGE_HEIGHT,
                mc.IMAGE_WIDTH,
                3)
  conv1_shape = get_conv2d_out_shape(input_shape,
                                    (7, 7, 96),
                                    [2, 2], 'VALID')
  pool1_shape = get_pool_out_shape(conv1_shape,
                                  (3, 3),
                                  [2, 2], 'VALID')
  fire2_shape = get_fire_out_shape(pool1_shape,
                                  s1x1=96, e1x1=64,
                                  e3x3=64)
  fire3_shape = get_fire_out_shape(fire2_shape,
                                  s1x1=96, e1x1=64,
                                  e3x3=64)
  fire4_shape = get_fire_out_shape(fire3_shape,
                                  s1x1=192, e1x1=128,
                                  e3x3=128)
  pool4_shape = get_pool_out_shape(fire4_shape,
                                  (3, 3),
                                  [2, 2], 'VALID')
  fire5_shape = get_fire_out_shape(pool4_shape,
                                  s1x1=192, e1x1=128,
                                  e3x3=128)
  fire6_shape = get_fire_out_shape(fire5_shape,
                                  s1x1=288, e1x1=192,
                                  e3x3=192)
  fire7_shape = get_fire_out_shape(fire6_shape,
                                  s1x1=288, e1x1=192,
                                  e3x3=192)
  fire8_shape = get_fire_out_shape(fire7_shape,
                                  s1x1=384, e1x1=256,
                                  e3x3=256)
  pool8_shape = get_pool_out_shape(fire8_shape,
                                  (3, 3),
                                  [2, 2], 'VALID')
  fire9_shape = get_fire_out_shape(pool8_shape,
                                  s1x1=384, e1x1=256,
                                  e3x3=256)
  fire10_shape = get_fire_out_shape(fire9_shape,
                                    s1x1=384, e1x1=256,
                                    e3x3=256)
  fire11_shape = get_fire_out_shape(fire10_shape,
                                    s1x1=384, e1x1=256,
                                    e3x3=256)
  num_output = mc.ANCHOR_PER_GRID * (mc.CLASSES + 1 + 4)
  conv12_shape = get_conv2d_out_shape(fire11_shape,
                                      (3, 3, num_output),
                                      [1, 1], 'SAME')
  return conv12_shape

def get_squeezeDet_shape(mc):
  input_shape = (mc.IMAGE_HEIGHT,
                mc.IMAGE_WIDTH,
                3)
  conv1_shape = get_conv2d_out_shape(input_shape,
                                    (3, 3, 64),
                                    [2, 2], 'SAME')
  pool1_shape = get_pool_out_shape(conv1_shape,
                                  (3, 3),
                                  [2, 2], 'SAME')
  fire2_shape = get_fire_out_shape(pool1_shape,
                                  s1x1=16, e1x1=64,
                                  e3x3=64)
  fire3_shape = get_fire_out_shape(fire2_shape,
                                  s1x1=16, e1x1=64,
                                  e3x3=64)
  pool3_shape = get_pool_out_shape(fire3_shape,
                                   (3, 3),
                                   [2, 2], 'SAME')
  fire4_shape = get_fire_out_shape(pool3_shape,
                                  s1x1=32, e1x1=128,
                                  e3x3=128)
  fire5_shape = get_fire_out_shape(fire4_shape,
                                  s1x1=32, e1x1=128,
                                  e3x3=128)
  pool5_shape = get_pool_out_shape(fire5_shape,
                                   (3, 3),
                                   [2, 2], 'SAME')
  fire6_shape = get_fire_out_shape(pool5_shape,
                                  s1x1=48, e1x1=192,
                                  e3x3=192)
  fire7_shape = get_fire_out_shape(fire6_shape,
                                  s1x1=48, e1x1=192,
                                  e3x3=192)
  fire8_shape = get_fire_out_shape(fire7_shape,
                                  s1x1=64, e1x1=256,
                                  e3x3=256)
  fire9_shape = get_fire_out_shape(fire8_shape,
                                  s1x1=64, e1x1=256,
                                  e3x3=256)
  fire10_shape = get_fire_out_shape(fire9_shape,
                                    s1x1=96, e1x1=384,
                                    e3x3=384)
  fire11_shape = get_fire_out_shape(fire10_shape,
                                    s1x1=96, e1x1=384,
                                    e3x3=384)
  num_output = mc.ANCHOR_PER_GRID * (mc.CLASSES + 1 + 4)
  conv12_shape = get_conv2d_out_shape(fire11_shape,
                                      (3, 3, num_output),
                                      [1, 1], 'SAME')
  return conv12_shape

max_stars_repo_name: supernlogn/squeezeDetTL

max_issues_repo_name: supernlogn/squeezeDetTL

max_forks_repo_name: supernlogn/squeezeDetTL


text:
from unittest import TestCase

from jtes import jaccard_timespan_event_score
import numpy as np

class JTESTest(TestCase):

    def test_both_empty(self):
        self.assertEqual(1, jaccard_timespan_event_score(np.array([]), np.array([])))

    def test_empty_prediction(self):
        y_true = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T01:00:00')),
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T03:00:10'))
        ])
        self.assertEqual(0, jaccard_timespan_event_score(y_true, np.array([])))

    def test_empty_true(self):
        y_pred = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T01:00:00')),
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T03:00:10'))
        ])
        self.assertEqual(0, jaccard_timespan_event_score(np.array([]), y_pred))

    def test_full_overlap(self):
        y = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T01:00:00')),
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T03:00:10'))
        ])

        self.assertEqual(1, jaccard_timespan_event_score(y, y))

    def test_single_half_overlap(self):
        y_true = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T01:00:00')),
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T04:00:00'))
        ])

        y_pred = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T01:00:00')),
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T05:00:00')),
        ])
        self.assertEqual(0.75, jaccard_timespan_event_score(y_true, y_pred))

    def test_double_half_overlap(self):
        y_true = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T01:00:00')),
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T04:00:00'))
        ])

        y_pred = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T02:00:00')),
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T05:00:00')),
        ])
        self.assertEqual(0.5, jaccard_timespan_event_score(y_true, y_pred))

    def test_false_positive(self):
        y_true = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T01:00:00')),
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T03:00:05')),
        ])

        y_pred = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T01:00:00')),
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T03:00:05')),
            (np.datetime64('1901-01-01T03:00:00'), np.datetime64('1901-01-01T03:00:05')),
        ])
        self.assertEqual(2 / 3, jaccard_timespan_event_score(y_true, y_pred))

    def test_false_negative(self):
        y_true = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T01:00:00')),
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T03:00:05')),
            (np.datetime64('1901-01-01T03:00:00'), np.datetime64('1901-01-01T03:00:05')),
        ])

        y_pred = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T01:00:00')),
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T03:00:05')),
        ])
        self.assertEqual(2 / 3, jaccard_timespan_event_score(y_true, y_pred))

    def test_duplicates(self):
        y_true = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T01:00:00')),
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T01:00:00')),  # dub
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T03:00:05')),
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T03:00:05')),  # dub
        ])

        y_pred = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T01:00:00')),
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T03:00:05')),
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T03:00:05')),  # dub
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T03:00:05')),  # dub
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T03:00:05')),  # dub
        ])
        with self.assertRaises(ValueError):
            jaccard_timespan_event_score(y_true, y_pred)

    def test_simple_pred_split(self):
        y_true = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T01:00:00')),
        ])

        y_pred = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T00:20:00')),
            (np.datetime64('1900-01-01T00:20:00'), np.datetime64('1900-01-01T01:00:00')),
        ])
        self.assertEqual(0.5, jaccard_timespan_event_score(y_true, y_pred))

    def test_pred_split(self):
        y_true = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T01:00:00')),
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T04:00:00'))
        ])

        y_pred = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T00:20:00')),
            (np.datetime64('1900-01-01T00:20:00'), np.datetime64('1900-01-01T01:00:00')),
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T03:10:00')),
            (np.datetime64('1900-01-01T03:10:00'), np.datetime64('1900-01-01T03:45:00')),
            (np.datetime64('1900-01-01T03:45:00'), np.datetime64('1900-01-01T04:00:00'))
        ])
        self.assertEqual(((20/60 + 40/60)/2 + (10/60 + 35/60 + 15/60)/3)/2,
                         jaccard_timespan_event_score(y_true, y_pred))

    def test_pred_overlap(self):
        y_true = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T01:00:00')),
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T04:00:00'))
        ])

        y_pred = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T01:00:00')),
            (np.datetime64('1900-01-01T00:30:00'), np.datetime64('1900-01-01T01:00:00')),
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T03:18:00')),
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T04:00:00')),
            (np.datetime64('1900-01-01T03:10:00'), np.datetime64('1900-01-01T03:44:00'))
        ])
        self.assertEqual(((0.5+1)/2 + (18/60+1+34/60)/3)/2, jaccard_timespan_event_score(y_true, y_pred))

        y_pred = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T04:00:00')),
        ])
        self.assertEqual((1/4 + 1/4)/2, jaccard_timespan_event_score(y_true, y_pred))

    def test_t0_lt_t1(self):
        j_wrong = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T01:00:00')),
            (np.datetime64('1900-01-01T03:00:10'), np.datetime64('1900-01-01T03:00:00'))  # t0 > t1
        ])
        y_correct = np.array([
            (np.datetime64('1900-01-01T00:00:00'), np.datetime64('1900-01-01T01:00:00')),
            (np.datetime64('1900-01-01T03:00:00'), np.datetime64('1900-01-01T03:00:10'))
        ])

        with self.assertRaises(ValueError):
            jaccard_timespan_event_score(j_wrong, np.array([]))
        with self.assertRaises(ValueError):
            jaccard_timespan_event_score(np.array([]), j_wrong)
        with self.assertRaises(ValueError):
            jaccard_timespan_event_score(j_wrong, j_wrong)
        with self.assertRaises(ValueError):
            jaccard_timespan_event_score(j_wrong, y_correct)
        with self.assertRaises(ValueError):
            jaccard_timespan_event_score(y_correct, j_wrong)

max_stars_repo_name: DEDDIAG/python-jtes

max_issues_repo_name: DEDDIAG/python-jtes

max_forks_repo_name: DEDDIAG/python-jtes


text:
# Author: Andrey Boytsov <andrey.boytsov@uni.lu> <andrey.m.boytsov@gmail.com>
# License: BSD 3 clause (C) 2017

# Visualizing 8x8 digits dataset (from sklearn). Visualizing first half, then transforming the rest.
# Feel free to play with starting_y

import sys
import os
from sklearn.datasets import load_digits
import matplotlib.pyplot as plt
import numpy as np

# Importing from parent directory
sys.path.append(os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))))
import lion_tsne

if __name__ == "__main__":
    data = load_digits()
    X_all = data.images.reshape((-1, 64)) # 8x8 image to 64-length vector
    X = X_all[:len(X_all)//2, :]
    X2 = X_all[len(X_all)//2:, :]
    labels_all = data.target
    labels = labels_all[:len(labels_all)//2]
    labels2 = labels_all[len(labels_all)//2:]

    dTSNE = lion_tsne.LionTSNE(perplexity=20)
    y = dTSNE.fit(X, verbose=2, optimizer_kwargs={'momentum': 0.8}, random_seed=0)
    # Option 1. Start with random values
    # starting_y = None
    # Option 2. Start with Ys that corresponded to closest original Xs
    starting_y = 'closest'
    # Option 3. Start with a center of class. Only for testing, it won't be available in training.
    #starting_y = np.zeros((len(labels2), 2))
    #for i in range(len(starting_y)):
    #    starting_y[i,:] = np.mean(y[labels == labels2[i], :], axis=0)
    y2 = dTSNE.transform(X2, y=starting_y, verbose=2,
                         optimizer_kwargs={'momentum': 0.8, 'n_iter': 3000}, random_seed=1)

    color_list = ['blue', 'red', 'green', 'yellow', 'cyan', 'black', 'magenta', 'pink', 'brown', 'orange']

    plt.gcf().set_size_inches(10, 10)
    legend_list = list()
    for l in set(sorted(labels)):
        plt.scatter(y[labels == l, 0], y[labels == l, 1], c=color_list[l])
        legend_list.append(str(data.target_names[l]))
    for l in set(sorted(labels)):
        plt.scatter(y2[labels2 == l, 0], y2[labels2 == l, 1], c=color_list[l], marker='v')
        legend_list.append(str(data.target_names[l])+" embedded")
    # It did not work. Clusters were elsewhere.
    #plt.xlim([-70,90])
    #plt.ylim([-90,70])
    plt.legend(legend_list)
    plt.show()

max_stars_repo_name: andreyboytsov/DynamicTSNE

max_issues_repo_name: andreyboytsov/DynamicTSNE

max_forks_repo_name: andreyboytsov/DynamicTSNE


text:
# -*- coding: utf-8 -*-
"""
Created on Wed Aug  5 00:43:27 2020

@author: jisuk
"""

# %%  import importance module

from __future__ import absolute_import, division, print_function

import tensorflow as tf
from tensorflow.keras import Model, layers
import numpy as np

# %% MNIST dataset parameters
num_classes = 10  # total classes (0-9 digits)

# training parameters
learning_rate = 0.001
training_steps = 500
batch_size = 128
display_step = 10

# network parameters
conv1_filters = 32  # number of filters for 1st conv layer
conv2_filters = 64  # number of filters for 2nd conv layer
fc1_units = 1024  # number of neurons for 1st fully-connected layer

# %% prepare NIMST data
#from tensorflow.keras.datasets import mnist
from tensorflow.keras.datasets import mnist
(x_train, y_train), (x_test, y_test) = mnist.load_data()
# convert to float32
x_train, x_test = np.array(x_train, np.float32), np.array(x_test, np.float32)
# normalize images value from [0 255] to [0 1]
x_train, x_test = x_train / 255.0, x_test/255.0

# %% use tf.data API to shuffle and batch data.
train_data = tf.data.Dataset.from_tensor_slices((x_train, y_train))
train_data = train_data.repeat().shuffle(5000).batch(batch_size).prefetch(1)

# %% Create TF Model

class ConvNet(Model):
    # set layers
    def __init__(self):
        super(ConvNet, self).__init__()
        # convolutional layer with 32 filters and a kernel size o f 5
        self.conv1 = layers.Conv2D(32, kernel_size=5, activation=tf.nn.relu)
        # max pooling (down-sampling) with kernel size of 2 and strides of 2.
        self.maxpool1 = layers.MaxPool2D(2, strides=2)

        # convolution layer with 64 filters and a kernel size of 3
        self.conv2 = layers.Conv2D(64, kernel_size=3, activation=tf.nn.relu)
        # max pooling (down-sampling) with kernel size of 2 and strides of 2.
        self.maxpool2 = layers.MaxPool2D(2, strides=2)

        # Flatten thedata to a 1-D vector for the fully connected layer.
        self.flatten = layers.Flatten()

        # fully connected layer.
        self.fc1 = layers.Dense(1024)
        # apply dropout (tf is_training is False, dropout is not applied)
        self.dropout = layers.Dropout(rate=0.5)

        # output layer class prediction
        self.out = layers.Dense(num_classes)

    # set forward pass
    def call(self, x, is_training=False):
        x = tf.reshape(x, [-1, 28, 28, 1])
        x = self.conv1(x)
        x = self.maxpool1(x)
        x = self.conv2(x)
        x = self.maxpool2(x)
        x = self.flatten(x)
        x = self.fc1(x)
        x = self.dropout(x, training=is_training)
        x = self.out(x)
        if not is_training:
            # tf cross entropy expect logits without softmax so only
            # apply softmax when not training
            x = tf.nn.softmax(x)
        return x

# build neural network model.
conv_net = ConvNet()

# %% cross entrolpy loss
# note that this will apply softmax to the logits

def cross_entropy_loss(x, y):
    # convert labels to int 64 for tf cross-entropy function.
    y = tf.cast(y, tf.int64)
    # apply softmax to logits and compute cross-entropy
    loss = tf.nn.sparse_softmax_cross_entropy_with_logits(labels=y, logits=x)
    # average loss across the batch.
    return tf.reduce_mean(loss)

# accuracy metric

def accuracy(y_pred, y_true):
    # predicted class is the index fo highest score in prediction vector
    correct_prediction = tf. equal(
        tf.argmax(y_pred, 1), tf.cast(y_true, tf.int64))
    return tf.reduce_mean(tf.cast(correct_prediction, tf.float32), axis=-1)

# stochasitc gradient descent optimizer
optimizer = tf.optimizers.Adam(learning_rate)

# %% optimization result

def run_optimization(x, y):
    # wrap computation inside a gradientTape for automatic differentiation
    with tf.GradientTape() as g:
        # forward pass
        pred = conv_net(x, is_training=True)
        # compute loss
        loss = cross_entropy_loss(pred, y)

    # variables to update, i.e. trainable variables
    trainable_variables = conv_net.trainable_variables

    # compute gradients
    gradients = g.gradient(loss, trainable_variables)

    # update W and b following gradients
    optimizer.apply_gradients(zip(gradients, trainable_variables))

# %% run training for the given number of steps
for step, (batch_x, batch_y) in enumerate(train_data.take(training_steps), 1):
    # run the optimization to update W and b values.
    run_optimization(batch_x, batch_y)

    if step % display_step == 0:
        pred = conv_net(batch_x)
        loss = cross_entropy_loss(pred, batch_y)
        acc = accuracy(pred, batch_y)
        print("step: %i, loss: %f, accuracy: %f" % (step, loss, acc))

# %% test model on validation set
pred = conv_net(x_test)
print("test Accuracy: %f" % accuracy(pred, y_test))

# %% visualize predictions
#import matplotlib.pyplot as plt
import matplotlib.pyplot as plt
# %% predict 5 images from validation set
n_images = 5
test_images = x_test[:n_images]
predictions = conv_net(test_images)

# display image and model prediction
for i in range(n_images):
    plt.imshow(np.reshape(test_images[i], [28, 28]), cmap='gray')
    plt.show()
    print("Model prediction: %i" % np.argmax(predictions.numpy()[i]))

max_stars_repo_path: tensorflow examples/3 neural networks/convolutional_network.py


text:
import torch
from config import config as cfg
import backbones
import logging
import random
import losses
import torch
import torch.distributed as dist
import torch.nn.functional as F
import torch.nn as nn
import torch.utils.data.distributed
from utils.utils_logging import AverageMeter, init_logging
import copy
from utils.utils_callbacks import CallBackVerification, CallBackLogging, CallBackModelCheckpoint
from eval_local import CallBack_LocalVerifi
from tqdm import tqdm
import gc
import copy
import os
import time
import pickle
import numpy as np
from functools import reduce

def FedPavg(models,weights):
    aggr = copy.deepcopy(models[0])
    weights = [ w/sum(weights) for w  in weights]
    with torch.no_grad():
        for name in aggr:
            tmp  = 0
            for i in range(len(models)):
                tmp += weights[i] * models[i][name]
            aggr[name] = tmp
    return aggr

def FedAvg_on_FC(pretrain_fc,models,weights,p):
    weights = [ w/sum(weights) for w  in weights]
    aggr = copy.deepcopy(models[0]) * weights[0]
    with torch.no_grad():
        for i in range(1,len(models)):
            aggr += models[i] * weights[i]
    if p == 1:
        pretrain_fc = aggr
    else:
        pretrain_fc = (1-p)*pretrain_fc + p * aggr
    return pretrain_fc

class SpreadOut_Module(nn.Module):
    def __init__(self,all_FC,margin=0.7,local=False,mode='sum'):
        super(SpreadOut_Module,self).__init__()
        self.FC = nn.Parameter(all_FC)
        self.margin = margin
        self.mode = mode

    def forward(self):
        FC = F.normalize(self.FC)
        similarity = torch.matmul(FC,FC.t())
        loss = F.relu(similarity.masked_select(~torch.eye(len(FC), dtype=bool).cuda())-self.margin)
        if self.mode == 'sum':
            loss = torch.sum(loss**2)
        elif self.mode == 'mean':
            loss = torch.mean(loss**2)
        return loss

class Server(object):
    def __init__(self,clients,data,args):
        self.data = data
        self.output_dir = args.output_dir
        self.args = args
        self.local_epoch = args.local_epoch
        self.clients = clients
        self.num_client = args.num_client
        self.csr = args.client_sampled_ratio
        self.logger = logging.getLogger('FL_face.server')
        if args.add_pretrained_data:
            self.public_train_loader = data.public_train_loader
            self.public_test_loader = data.public_test_loader
        ### Create backbone
        self.dropout = 0.4 if cfg.dataset is "webface" else 0
        self.federated_model = eval("backbones.{}".format(args.network))(False, dropout=self.dropout, fp16=cfg.fp16) # on cpu
        self.margin_softmax = eval("losses.{}".format(args.loss))(s=30.0,m=0.4)

        ## Load pretrained backbone
        if self.args.pretrained_root is not '':
            model = torch.load(self.args.pretrained_bb_path,map_location='cpu')
            self.federated_model.load_state_dict(model)
            self.logger.info('Succesfully load model from %s'%(self.args.pretrained_bb_path))
        else:
            self.logger.info('Train from scratch!')
        self.federated_model.eval()
        ###
        self.train_loss = []
        self.global_epoch = 0
        self.global_round = 0

        self.rank = 0 #dist.get_rank()
        self.local_rank = 0 #args.local_rank

        ## 1:1 verify on small dataset (cfp, lfw...)
        self.callback_verification = CallBackVerification(1, self.rank, cfg.val_targets, cfg.val_rec, self.num_client)

        ## 1:1 verify on our local dataset, sample 10
        self.callback_local_veri = CallBack_LocalVerifi(frequent=1, rank=self.rank,data_dir=cfg.local_rec,output_dir=self.output_dir)
        self.local_candidates = sorted(list(np.random.permutation(self.num_client)[:10]))
        print('Local Veri Candidates',self.local_candidates)

        self.callback_checkpoint = CallBackModelCheckpoint(self.rank, self.output_dir)
        self.current_client_list = None

        ## local FC initialization
        self.norm_before_avg = False
        if args.init_fc:
            self.Initialize_local_FC(False)

        ### pretrained FC initialization
        if args.add_pretrained_data:
            if self.args.init_fc:
                self.pretrained_fc, self.pretrained_label = self.Initialize_pretrain_FC(only_labels=False)
            else:
                _,self.pretrained_label = self.Initialize_pretrain_FC(only_labels=True)
                self.pretrained_fc = torch.load(self.args.pretrained_fc_path,map_location='cpu')
                self.logger.info('Use pretrained perfect FC')

        ### BCE branch
        if self.args.BCE_local:
            # self.logger.info('No init BCE')
            self.logger.info('Init BCE')
            ## No matter init_fc or not, BCE is initialized with the same weight
            for i in range(len(self.clients)):
                self.clients[i].bce_module.initialize(self.clients[i].fc_module.fc.data)

    def test(self):
        with torch.no_grad():
            self.federated_model.to(self.local_rank)
            self.federated_model = nn.DataParallel(self.federated_model)
            self.federated_model.eval()
            self.callback_verification(self.global_round, self.federated_model, None, th=0)
            # only rank 0 has highest_acc_list
            if self.callback_verification.highest_acc_list[-1][0] == self.global_round:
                self.callback_checkpoint(self.global_round,self.federated_model,None)
                self.logger.info('Save server model, epoch %d model...'%(self.global_round))
            self.federated_model = self.federated_model.module.cpu()

            if self.global_round >= 0 and self.global_round % 1 == 0:
                torch.save(self.federated_model.state_dict(), os.path.join(self.output_dir, "backbone_%d.pth"%self.global_round))

            # if self.local_rank == 0 and (self.global_epoch+1)%5 == 0:
                ## IJBC
                # torch.save(self.federated_model.state_dict(), os.path.join(self.output_dir, "ijbc_tmp.pth"))
                # os.system('CUDA_VISIBLE_DEVICES=%s sh run_IJBC.sh %s %s %d %s'\
                        # %(self.args.gpus,os.path.join(self.output_dir,'ijbc_tmp.pth'),\
                        #   self.output_dir,self.global_epoch,os.path.join(cfg.val_rec,'IJBC')))

    def Initialize_local_FC(self,save_to_disk=False):
        self.args.pre_init_path = os.path.join(self.args.pretrained_root,'preCos_init_AN.pth')
        ### Load pre-forward features from pretrained model
        if os.path.exists(self.args.pre_init_path):
            self.logger.info('Preload Clients FC init.')
            init_matrix = torch.load(self.args.pre_init_path,map_location='cpu')
            self.logger.info('Clients init shape %r'%list(init_matrix.shape))
            start = 0
            for i in range(len(self.clients)):
                num_classes = self.clients[i].num_classes
                self.clients[i].fc_module.update_from_tensor(init_matrix[start:start+num_classes,:])
                start += num_classes
        ### Infer the pretrained model to generate features
        else:
            collected_fc = []
            for i in range(len(self.clients)):
                self.logger.info('Client %d start initialize!'%(i))
                self.clients[i].data_update_fc(self.federated_model.state_dict(),self.norm_before_avg, \
                                                fc_name='cen_feats_init',save_to_disk=save_to_disk)
                collected_fc.append(self.clients[i].fc_module.fc.data)
            collected_fc = torch.cat(collected_fc,dim=0)
            # torch.save(collected_fc,self.args.pre_init_path)
            gc.collect()
            torch.cuda.empty_cache()

    def Initialize_pretrain_FC(self,load_pth=True,save_pth=False,only_labels=False):
        fc_name = os.path.join(self.args.pretrained_root,'preCos_pretrain_init_AN.pth')
        labels_name = os.path.join(self.args.pretrained_root,'preCos_pretrain_labels.pth')

        if os.path.exists(fc_name) and load_pth == True :
            raw_labels = torch.load(labels_name,map_location='cpu')
            self.logger.info('Preload pretrain labels, shape: %r'%list(raw_labels.shape))
            if only_labels:
                return None,raw_labels

            init_matrix = torch.load(fc_name, map_location='cpu')
            self.logger.info('Preload pretrain fc, shape: %r'%list(init_matrix.shape))
            return init_matrix, raw_labels

        ## forward backbone to create features
        ### Creating model
        backbone = eval("backbones.{}".format(self.args.network))(False, dropout=self.dropout, fp16=cfg.fp16)
        backbone.load_state_dict(self.federated_model.state_dict())
        backbone = nn.DataParallel(backbone)
        backbone.to(self.local_rank)
        backbone.eval()
        # forward public data
        raw_labels = []
        ID_feature = dict()
        self.logger.info('Initialize pretrained fc on server')
        with torch.no_grad():
            pbar = tqdm(total=len(self.public_test_loader),ncols=120)
            for step, (img, label) in enumerate(self.public_test_loader):
                raw_labels.append(label.cpu())
                if not only_labels:
                    features = backbone(img.to(self.local_rank))
                    if self.norm_before_avg:
                        features = F.normalize(features)
                    for i,ID in enumerate(label):
                        ID = ID.item()
                        if ID not in ID_feature:
                            ID_feature[ID] = [torch.zeros_like(features[0]),0]
                        ID_feature[ID][0] += features[i]
                        ID_feature[ID][1] += 1
                pbar.update(1)
            pbar.close()

        raw_labels = torch.cat(raw_labels)
        self.logger.info('Generate pretrain labels, %r'%list(raw_labels.shape))
        if only_labels:
            del backbone
            return None, raw_labels

        init_matrix = torch.zeros(len(ID_feature),512).float()
        for i in range(len(init_matrix)):
            ## average
            init_matrix[i] = (ID_feature[i][0] / ID_feature[i][1]).cpu()
        # torch.save(init_matrix,fc_name)
        # torch.save(raw_labels, raw_name)
        self.logger.info('Generate pretrain fc, %r'%list(init_matrix.shape))
        del ID_feature, backbone
        gc.collect()
        torch.cuda.empty_cache()
        return init_matrix, raw_labels

    def Generate_pretrain_feats(self):
        self.logger.info('Generating pretrained features for HN sampling')
        backbone = eval("backbones.{}".format(self.args.network))(False, dropout=self.dropout, fp16=cfg.fp16)
        backbone.load_state_dict(self.federated_model.state_dict())
        backbone = nn.DataParallel(backbone)
        backbone.to(self.local_rank)
        backbone.eval()
        # forward public data
        raw_feats = []
        with torch.no_grad():
            pbar = tqdm(total=len(self.public_test_loader),ncols=120)
            for step, (img, label) in enumerate(self.public_test_loader):
                features = backbone(img.to(self.local_rank))
                features = F.normalize(features)
                raw_feats.append(features.cpu())
                pbar.update(1)
            pbar.close()
        raw_feats = torch.cat(raw_feats,dim=0)
        # torch.save(raw_feats,'tmp.pth')
        backbone = backbone.cpu()
        del backbone
        return raw_feats

    def train(self):
        models = []
        models_fc = []
        losses = []
        data_sizes =[]

        ## for feature-based HN
        if self.args.add_pretrained_data:
            # self.logger.info('load tmp.pth')
            # self.pretrained_feats = torch.load('tmp.pth')
            self.pretrained_feats = self.Generate_pretrain_feats()

        ### Adjust Local epoch and decay
        if self.args.adaptive_local_epoch and self.global_round != 0:
            self.local_epoch = max(4,self.local_epoch-2)
            cfg.train_decay = max(1,int(3/4*self.local_epoch))

        #### start train
        for idx,i in enumerate(self.current_client_list):
            self.logger.info('Round %d : [%d/%d] Client %d start training!'%(self.global_round,idx+1,len(self.current_client_list),i))
            self.logger.info('Server send backbone to clients')
            self.clients[i].backbone_state_dict = self.federated_model.state_dict()

            ## Adjust local epoch
            self.clients[i].local_epoch = self.local_epoch

            ## Train with local verify test
            if self.clients[i].cid in self.local_candidates:
                if self.args.add_pretrained_data:
                    self.clients[i].train_with_public_data(self.global_epoch,callback_verification=self.callback_local_veri,\
                        public_train_loader=self.public_train_loader,pretrained_fc=self.pretrained_fc,choose_hard_negative=True,\
                        pretrained_label=self.pretrained_label,pretrained_feats=self.pretrained_feats)
                else:
                    self.clients[i].train(self.global_epoch,callback_verification=self.callback_local_veri)
            ## train w/o local test
            else:
                if self.args.add_pretrained_data:
                    self.clients[i].train_with_public_data(self.global_epoch,\
                        public_train_loader=self.public_train_loader,pretrained_fc=self.pretrained_fc,choose_hard_negative=True,\
                        pretrained_label=self.pretrained_label,pretrained_feats=self.pretrained_feats)
                else:
                    self.clients[i].train(self.global_epoch)
            #######################################################
            losses.append(self.clients[i].get_train_loss())

            if self.args.return_all:
                models.append(self.clients[i].get_model())
                models_fc.append(self.clients[i].get_global_fc())
                self.clients[i].fc_module.remove_pretrain()
            else:
                models.append(self.clients[i].get_model())
            data_sizes.append(self.clients[i].get_data_size())

        avg_loss = sum(losses)/len(losses)
        self.logger.info("================")
        self.logger.info("Train Round {}. Avg Train loss among all clients {:.6f}".format(self.global_round,avg_loss))

        if self.args.return_all:
            p = 1.0
            self.logger.info('==========Do Fed FC==========')
            self.pretrain_fc = FedAvg_on_FC(self.pretrained_fc, models_fc, data_sizes, p=p)
            self.logger.info('==========Do FedPavg==========')
            if self.args.aggr_alg in ['FedAvg','FedProx']:
                aggr_state_dict = FedPavg(models,data_sizes)
                if p != 1.0:
                    global_state_dict = self.federated_model.state_dict()
                    for name in global_state_dict:
                        aggr_state_dict[name] = (1-p)*global_state_dict[name] + p * aggr_state_dict[name]
            self.federated_model.load_state_dict(aggr_state_dict)
        else:
            self.logger.info('==========Do FedPavg==========')
            if self.args.aggr_alg in ['FedAvg','FedProx']:
                aggr_state_dict = FedPavg(models,data_sizes)
            self.federated_model.load_state_dict(aggr_state_dict)

    def SpreadOut(self,sp_iter = 5,mode='sum'):
        assert(self.current_client_list is not None)
        FC_over_clients = []
        for idx,i in enumerate(self.current_client_list):
            FC_over_clients.append(self.clients[i].fc_module.fc.data)
        FC_over_clients = torch.cat(FC_over_clients,dim=0)
        self.logger.info('=====Collect FC and cat to a big matrix=====')
        ##### Initialize the SpreadOut module
        SP = SpreadOut_Module(FC_over_clients.detach().clone(),margin=0.4,mode=mode).to(self.local_rank)
        SP_opt = torch.optim.SGD(SP.parameters(), lr=cfg.lr*10,momentum=0.9,weight_decay=cfg.weight_decay)
        self.logger.info('=====SpreadOut Module Create=====')
        for i in range(sp_iter):
            SP_opt.zero_grad()
            loss = SP()
            self.logger.info('- SP iter %d Loss :  %.5e , Start backward'%(i, loss.item()))
            loss.backward()
            SP_opt.step()
        #####
        with torch.no_grad():
            FC = SP.FC.data.cpu()
        for idx,i in enumerate(self.current_client_list):
            self.clients[i].fc_module.fc.data = FC[idx*self.clients[i].num_classes:(idx+1)*self.clients[i].num_classes]
        #####
        # if (self.global_epoch) % self.args.save_fc_iter == 0:
            # torch.save(FC_over_clients,os.path.join(self.output_dir,'Ep%d_allFC_beforeSP.pth'%self.global_epoch))
            # torch.save(FC,os.path.join(self.output_dir,'Ep%d_allFC_afterSP.pth'%self.global_epoch))
        self.logger.info('=====Update FC in partial FC module=====')
        ## delete garbage
        del SP,FC,FC_over_clients
        torch.cuda.empty_cache()
        gc.collect()
        return


text:
from urllib.request import Request, urlopen
import json
import os
import numpy as np
import pandas as pd
from datetime import datetime, timedelta
from matplotlib import pyplot as plt
import windrose
import random
from climvis import core, cfg

def test_get_ts():

    df_cities = pd.read_csv(cfg.world_cities)
    dfi = df_cities.loc[df_cities.Name.str.contains('innsbruck', case=False,
                                                    na=False)].iloc[0]

    df = core.get_cru_timeseries(dfi.Lon, dfi.Lat)
    assert df.grid_point_elevation > 500  # we are in the alps after all
    assert df.distance_to_grid_point < 50000  # we shouldn't be too far

    # It's different data but I wonder how we compare to the
    # Innsbruck climate station we studied a couple of weeks ago?
    url = ('https://raw.githubusercontent.com/fmaussion/'
           'scientific_programming/master/data/innsbruck_temp.json')
    req = urlopen(Request(url)).read()

    # Read the data
    data = json.loads(req.decode('utf-8'))
    for k, v in data.items():
        data[k] = np.array(data[k])

    # select
    t = data['TEMP'][np.nonzero((data['YEAR'] <= 2016))]
    yrs = data['YEAR'][np.nonzero((data['YEAR'] <= 2016))]
    dfs = df.loc[(df.index.year >= yrs.min()) &
                 (df.index.year <= yrs.max())].copy()
    assert len(dfs) == len(yrs)
    dfs['ref'] = t
    dfs = dfs[['tmp', 'ref']]

    # Check that we have good correlations at monthly and annual scales
    assert dfs.corr().values[0, 1] > 0.95
    assert dfs.groupby(dfs.index.year).mean().corr().values[0, 1] > 0.9

    # Check that altitude correction is helping a little
    z_diff = df.grid_point_elevation - dfi.Elevation
    dfs['tmp_cor'] = dfs['tmp'] + z_diff * 0.0065
    dfm = dfs.mean()
    assert np.abs(dfm.ref - dfm.tmp_cor) < np.abs(dfm.ref - dfm.tmp)

def test_get_url():

    df_cities = pd.read_csv(cfg.world_cities)
    dfi = df_cities.loc[df_cities.Name.str.contains('innsbruck', case=False,
                                                    na=False)].iloc[0]

    url = core.get_googlemap_url(dfi.Lon, dfi.Lat)
    assert 'maps.google' in url

def test_write_html(tmpdir):

    df_cities = pd.read_csv(cfg.world_cities)
    dfi = df_cities.loc[df_cities.Name.str.contains('innsbruck', case=False,
                                                    na=False)].iloc[0]

    dir = str(tmpdir.join('html_dir'))
    core.write_html(dfi.Lon, dfi.Lat, directory=dir)
    assert os.path.isdir(dir)

def test_perc_dir_from_data():

    url = core.url_from_input('ellboegen', '3', core.base_url)
    req = urlopen(Request(url)).read()
    data = json.loads(req.decode('utf-8'))
    data['time'] = [datetime(1970, 1, 1) + timedelta(milliseconds=ds) for ds
                    in data['datumsec']]

    wind_dir_percent = core.perc_dir_from_data(data)
    # check if there is no na value in the computation of percentages
    assert np.all(~np.isnan(list(wind_dir_percent.values())))

def test_max_wind_speed():

    url = core.url_from_input('ellboegen', '3', core.base_url)
    req = urlopen(Request(url)).read()
    data = json.loads(req.decode('utf-8'))
    data['time'] = [datetime(1970, 1, 1) + timedelta(milliseconds=ds) for ds
                    in data['datumsec']]

    # I want to be sure that the minutes present are not different than
    # [0, 10, 20, 30, 40, 50]
    data['time'] = test_name_to_data()

    # check if nan values are put in the right positions
    removed_index = random.sample(range(1, len(data['time'])-1), 10)
    time = np.delete(data['time'], removed_index)

    check = []
    for i in range(len(time)):
        check.append(time[i].minute)

    possible_minutes = [0, 10, 20, 30, 40, 50]
    position = possible_minutes.index(check[0])

    # get the sequence that has to repeat from the initial datetime minute
    sequence = []

    for i in range(position, position + 6):
        min = possible_minutes[i % 6]
        sequence.append(min)

    position = 0
    # convert to float to insert nan values
    check = np.array(check, dtype=float)
    sequence = np.array(sequence, dtype=float)

    ff = np.array(data['ff'])
    time_array = np.array(time)

    # check if every minute corresponds to the minute in the sequence
    # if not (so there is a jump of ten minutes in check) a nan is inserted
    # the loop ends when it arrives to the last element of array
    i = 0
    while True:

        if not check[i] == sequence[position]:
            check = np.insert(check, i, np.nan)
            time_array = np.insert(time_array, i, np.nan)
            ff = np.insert(ff, i, np.nan)
        if position < 5:
            position += 1
        else:
            position = 0
        i += 1
        try:
            check[i]
        except IndexError:
            break

    nan_indexes = np.argwhere(np.isnan(ff)).tolist()
    # np.argwhere output is a list of lists
    # itertools allows to make a list out of a list of lists
    import itertools
    nan_indexes = list(itertools.chain.from_iterable(nan_indexes))
    assert sorted(removed_index) == nan_indexes

def test_name_to_data():
    url = core.url_from_input('ellboegen', '3', core.base_url)
    req = urlopen(Request(url)).read()
    data = json.loads(req.decode('utf-8'))
    data['time'] = [datetime(1970, 1, 1) + timedelta(milliseconds=ds) for ds
                    in data['datumsec']]

    check = []
    possible_minutes = [0, 10, 20, 30, 40, 50]

    for i in range(0, len(data['time'])):
        check.append(data['time'][i].minute)

    check = np.array(check)
    boolean = np.isin(check, possible_minutes)

    # check if all the wrong minutes are corrected

    if False in boolean:
        indexes = list(np.where(boolean == False))[0].tolist()
        for i in indexes:
            right_min = data['time'][i].minute * 10
            if right_min < 60:
                data['time'][i] = data['time'][i].replace(minute=right_min)
            else:
                del data['time'][i]
            assert data['time'][i].minute in possible_minutes
    return data['time']

def test_windrose_data():
    url = core.url_from_input('ellboegen', '3', core.base_url)
    req = urlopen(Request(url)).read()
    data = json.loads(req.decode('utf-8'))
    data['time'] = [datetime(1970, 1, 1) + timedelta(milliseconds=ds) for ds
                    in data['datumsec']]

    fig = plt.figure(figsize=(6, 4))
    ax = core.windrose_data(data['dd'], data['ff'], fig)

    assert type(ax) == windrose.windrose.WindroseAxes

def test_write_html_windrose(tmpdir):

    station, days = 'ellboegen', '3'

    dir = str(tmpdir.join('html_windrose_dir'))
    core.write_html_wind_rose(station, days, directory=dir)
    # check if it exists the directory in which index_windrose.html is created
    assert os.path.isdir(dir)


text:
### FIRE DOMAIN ONLY

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import cPickle
from collections import defaultdict
import sys, os, math, re
os.environ['KERAS_BACKEND']='tensorflow'
from keras.preprocessing.sequence import pad_sequences
from keras.layers import LSTM, Dense
from keras.models import Sequential
from sklearn.metrics import mean_squared_error
import matplotlib

from keras.models import model_from_json
from pandas import HDFStore

max_sen_len = 100

dir_list = [ './Data_20_2hyp/', './Data_20_4hyp/', './Data_50_2hyp/', './Data_50_4hyp/', './Data_100_2hyp/', './Data_100_4hyp/' ]
# dir_list = [ 'Data_50_2hyp/', 'Data_50_4hyp/', 'Data_100_2hyp/', 'Data_100_4hyp/' ]

input_flist = [ 'Inputs/fire2.csv', 'Inputs/fire5.csv', 'Inputs/fire10.csv', 'Inputs/fire20.csv', 'Inputs/fire50.csv']
output_flist = [ 'Outputs/res_fire2.npy', 'Outputs/res_fire5.npy', 'Outputs/res_fire10.npy', 'Outputs/res_fire20.npy', 'Outputs/res_fire50.npy']
preds_flist = [ 'Preds/pred_fire2.npy', 'Preds/pred_fire5.npy', 'Preds/pred_fire10.npy', 'Preds/pred_fire20.npy', 'Preds/pred_fire50.npy']

# Iterate through dir_list using dir_name
# Iterate through input_flist using input_fname
# Iterate through output_flist using output_fname
# Iterate through preds_flist using preds_fname

for dir_name in dir_list:

	for index_ in xrange(5):

		input_fname = input_flist[index_]
		output_fname = output_flist[index_]
		preds_fname = preds_flist[index_]

		print (str (dir_name + input_fname))

		dataset_inter = []
		ms = open(str(dir_name + input_fname),'r')
		for line in ms.readlines():
		    val = line.split()
		    val = [float(i) for i in val]
		    dataset_inter.append(val)
		dataset_inter = np.array(dataset_inter)

		dataset = []
		dataset = np.array(dataset_inter)
		dataset = dataset.astype('float32')

		# x=[]
		# for i in xrange(0,100):
		# 	x.append(i)
		# plt.figure('data')
		# for i in dataset:
		# 	plt.plot(x,i)
		# plt.show()

		np.random.seed(7)
		# np.random.shuffle(dataset)
		dataset = np.reshape(dataset, (dataset.shape[0], dataset.shape[1], 1))

		#################################
		# split into train and test sets
		val_split = 0.3
		train_size = int(len(dataset) * (1.0 - val_split))
		test_size = len(dataset) - train_size

		train, test = dataset[0:train_size,:], dataset[train_size:len(dataset),:]

		x_train = train[:,:-1]
		y_train = train[:,1:]

		x_test = test[:,:-1]
		y_test = test[:,1:]
		#################################

		#x, y = dataset[:,:-1], dataset[:,1:]

		############### change model from sequence prediction to last value prediction ###################

		y_train_last = [[y_train[i,-1] for j in xrange(0,y_train.shape[1])] for i in xrange(0,y_train.shape[0])]

		y_test_last = [[y_test[i,-1] for j in xrange(0,y_test.shape[1])] for i in xrange(0,y_test.shape[0])]

		y_train_last = np.array(y_train_last)
		y_test_last = np.array(y_test_last)

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

		model = Sequential()
		model.add (LSTM (1, input_shape = (99,1), return_sequences=True, stateful=False, dtype='float32', dropout=0.3, recurrent_dropout=0.3))
		model.compile(loss="mse", optimizer="adam", metrics=['accuracy'])
		# model.fit (x_train, y_train, batch_size=1, epochs=20, validation_data=(x_test, y_test), shuffle=True)
		model.fit (x_train, y_train_last, batch_size=1, epochs=10, validation_data=(x_test, y_test_last), shuffle=True)
		predictions = model.predict(x_test)

		with open(str(dir_name + preds_fname),'wb') as pr:
			np.save(pr,predictions)

		with open(str(dir_name + output_fname),'wb') as res:
			np.save(res,y_test_last)

		# a = np.load(fp)

		# serialize model to JSON

		'''
		model_json = model.to_json()
		with open("model_fire10.json", "w") as json_file:
		    json_file.write(model_json)
		# serialize weights to HDF5

		'''

		# ######## model.save_weights("model_boul10.h5")
		# print("Saved model to disk")

		# later...x

		# # load json and create model
		# json_file = open('model_boul10.json', 'r')
		# loaded_model_json = json_file.read()
		# json_file.close()
		# loaded_model = model_from_json(loaded_model_json)
		# # load weights into new model
		# loaded_model.load_weights("model_boul10.h5")
		# print("Loaded model from disk")

		# # evaluate loaded model on test data
		# loaded_model.compile(loss='binary_crossentropy', optimizer='rmsprop', metrics=['accuracy'])
		# score = loaded_model.evaluate(X, Y, verbose=0)
		# print("%s: %.2f%%" % (loaded_model.metrics_names[1], score[1]*100))

max_stars_repo_name: PranavEranki/Time-Series-Experiments

max_issues_repo_name: PranavEranki/Time-Series-Experiments

max_forks_repo_name: PranavEranki/Time-Series-Experiments


text:
# -*- coding: utf-8 -*-
"""
Created on Tue Jul 31 16:13:20 2018

@author:    Jonas Hartmann @ Gilmour group @ EMBL Heidelberg

@descript:  Sklearn-like class for removal of outliers.
"""

#------------------------------------------------------------------------------

# IMPORTS

# External general
from __future__ import division
import numpy as np

# External specific
from sklearn.ensemble import IsolationForest

#------------------------------------------------------------------------------

### SKLEARN-LIKE TRANSFORM: OUTLIER REMOVAL

class RemoveOutliers(object):
    """Identifies and removes outliers from feature spaces based on one or
    more approaches/criteria selected from the following:

    - Value of any feature above or below a given threshold value
    - Value of any feature above or below a given percentile
    - Sum of absolute standardized feature values above a given percentile
    - Classified as outlier by sklearn.ensemble.IsolationForest

    Implemented as an sklearn-like transformer with these methods:
        .fit(X) -- determines the threshold on X.
        .transform(X[,y]) -- removes outliers from X [and y].
        .fit_transform(X[,y]) -- combines fit and transform.

    WARNING: Sklearn does not fully support transform operations on y.
             This class therefore cannot be used in sklearn pipelines.
    """

    def __init__(self, method, bounds='upper', threshold=None, percentile=95,
                 isoforest_params='default'):
        """Initialize RemoveOutliers instance.

        Parameters
        ----------
        method : string
            One of the following strings denoting the method to use:
                'absolute_thresh'
                    Remove any sample that has a feature value above/below a
                    given absolute threshold. The threshold can be a single
                    value for all features or a distinct value per feature.
                'percentile_thresh'
                    Remove any sample that has a feature value above/below a
                    given percentile threshold (determined based on training
                    data). The percentile threshold is computed for each
                    feature individually.
                'merged_percentile_thresh'
                    Remove any sample that has overall comparably extreme
                    values. Specifically, remove samples for which the sum
                    of the absolute standardized feature values is above a
                    given percentile threshold (determined based on training
                    data).
                'isolation_forest'
                    Use sklearn.ensemble.IsolationForest to determine and
                    remove outliers.
        bounds : string, optional, default 'upper'
            Relevant if method is 'absolute_thresh' or 'percentile_thresh';
            ignored otherwise.
            Bounds can be either 'upper' or 'lower', designating whether values
            above or below the threshold should be considered outliers. By
            default ('upper'), any samples with values above the threshold are
            removed as outliers.
        threshold : numeric or array of shape (N_features) or None, optional, default None
            Relevant if method is 'absolute_thresh'; ignored otherwise.
            Specifies either a single threshold for all features or an array of
            individual thresholds for each feature.
        percentile : int in range [1,99], optional, default 95
            Relevant if method is 'percentile_thresh' or
            'merged_percentile_thresh'; ignored otherwise.
            Percentile to compute the threshold.
        isoforest_params : dict or 'default', optional, default 'default'
            Parameters for instantiation of sklearn.ensemble.IsolationForest
            class. If 'default', sklearn's defaults are used. Otherwise must be
            a dictionary of keyword arguments.
        """

        # Bind parameters
        self.method           = method
        self.bounds           = bounds
        self.threshold        = threshold
        self.percentile       = percentile
        self.isoforest_params = isoforest_params

        # Check that method exists
        if method not in ['absolute_thresh', 'percentile_thresh',
                          'merged_percentile_thresh', 'isolation_forest']:
            raise NotImplementedError("Method not implemented: " + method)

        # Check that bounds argument makes sense (where relevant)
        if method in ['absolute_thresh', 'percentile_thresh']:
            if bounds not in ['upper', 'lower']:
                raise ValueError("bounds must be 'upper' or 'lower' for " +
                                 "method: " + method)

        # Check that threshold is defined (where relevant)
        if method=='absolute_thresh' and self.threshold is None:
            raise ValueError("threshold must not be None if method is " +
                             "absolute_thresh!")

    def fit(self, X):
        """Determine threshold for outlier detection.

        Parameters
        ----------
        X : array of shape (n_samples, n_features)
            Feature space used for threshold prediction.

        Returns
        -------
        self
        """

        # Absolute threshold
        if self.method == 'absolute_thresh':
            self.threshold_ = self.threshold

        # Percentile threshold
        elif self.method == 'percentile_thresh':
            percentile = np.percentile(X, self.percentile, axis=0)
            self.threshold_ = percentile

        # Merged percentile threshold
        elif self.method == 'merged_percentile_thresh':
            self.X_mean_ = X.mean(axis=0)
            self.X_std_  = X.std(axis=0)
            X_merged = (X - self.X_mean_) / self.X_std_
            X_merged = np.sum(np.abs(X_merged), axis=1)
            percentile = np.percentile(X_merged, self.percentile)
            self.threshold_ = percentile

        # Isolation forest
        elif self.method == 'isolation_forest':
            if self.isoforest_params == 'default':
                self.isoforest_ = IsolationForest()
            else:
                self.isoforest_ = IsolationForest(**self.isoforest_params)
            self.isoforest_.fit(X)

        # Done
        return self

    def transform(self, X, y=None):
        """Apply thresholds determined in fit to remove outliers from given
        feature space.

        Parameters
        ----------
        X : array of shape (n_samples, n_features)
            Outliers will be removed from this array.
        y : array of shape (n_samples, m_features) or list of such arrays, optional, default None
            Outliers determined based on X will also be removed from this array
            or form each array in an iterable of arrays.

        Returns
        -------
        Xt : array of shape (n_samples-n_outliers, n_features)
            Same as X but without outliers.
        yt : array of shape (n_samples-n_outliers, m_features) or list of such arrays
            Same as y but without outliers. Only returned if y is not None.
        """

        # Create inlier mask for absolute or percentile threshold
        if self.method in ['absolute_thresh', 'percentile_thresh']:
            if self.bounds == 'upper':
                outlier_mask = X > self.threshold_
            if self.bounds == 'lower':
                outlier_mask = X < self.threshold_
            self.is_inlier_ = ~ outlier_mask.max(axis=1)

        # Create inlier mask for merged percentile threshold
        elif self.method == 'merged_percentile_thresh':
            X_merged = X_merged = (X - self.X_mean_) / self.X_std_
            X_merged = np.sum(np.abs(X_merged), axis=1)
            outlier_mask = X_merged > self.threshold_
            self.is_inlier_ = ~ outlier_mask

        # Create inlier mask for isolation forest
        elif self.method == 'isolation_forest':
            is_inlier = self.isoforest_.predict(X)
            self.is_inlier_ = is_inlier == 1

        # Remove outliers from X
        Xt = X[self.is_inlier_, :]
        self.X_removed_ = X.shape[0] - Xt.shape[0]

        # Remove outliers from y
        if y is not None:

            if type(y) == list:
                yt = []
                for y_i in y:
                    yt.append(y_i[self.is_inlier_, :])

            else:
                yt = y[self.is_inlier_, :]

        # Return results
        if y is None:
            return Xt
        else:
            return Xt, yt

    def fit_transform(self, X, y=None):
        """Calls fit(X, y) followed by transform(X, y).
        See the respective methods for more information.
        """

        return self.fit(X).transform(X, y)

#------------------------------------------------------------------------------