text:
import inspect
from math import floor
import re

import numpy as np
import visa

class BaseInstr:
    def __init__(self, rsrc, idn=None, timeout_sec=5, reset=True):
        """
        >>> b = BaseInstr(None)
        debug mode (BaseInstr): skip BaseInstr.__init__.

        >>> rm = visa.ResourceManager()
        ... # rm = visa.ResourceManager('C:\\Windows\\system32\\visa32.dll')
        >>> rm.list_resources()

        >>> rsrc = rm.open_resource('visa://169.254.136.196/GPIB0::4::INSTR')
        >>> rsrc.timeout = 60*1000  # 60s
        >>> idn = 'LSCI,MODEL331S,333342,061404'
        >>> b = BaseInstr(rsrc, idn, reset=True)

        :param rsrc: Instrument VISA resource
        :param timeout_sec: second
        :type timeout_sec: int
        :type reset: bool
        :param idn: can be a regular expression.
        :type idn: str
        """
        self._debug_mode = True if rsrc is None else False
        if self._debug_mode:
            self._dbg_print('skip BaseInstr.__init__.')
        else:
            self._rsrc = rsrc
            self._rsrc.timeout = timeout_sec * 1000  # millisec
            resp_idn = self.q('*IDN?')
            match = re.search(idn, resp_idn)
            if match is None:
                raise RuntimeError('Failed to connect to a instrument. '
                                   '*IDN? response: {}'.format(resp_idn))
            if reset:
                self.reset()

    def _dbg_print(self, *objects, sep=' ', end='\n', file=None, flush=False):
        """
        Parameters: same as built-in print() function.

        >>> b = BaseInstr(None)
        debug mode (BaseInstr): skip BaseInstr.__init__.

        >>> b._dbg_print('Hello ', end='World!\\n')
        debug mode (BaseInstr): Hello World!
        """
        if not self._debug_mode:
            raise RuntimeError('Not debug mode.')
        print('debug mode ({}): '.format(self.__class__.__name__), end='')
        print(*objects, sep=sep, end=end, file=file, flush=flush)

    def _dbg_print_skip_method(self):
        """
        >>> b = BaseInstr(None)
        debug mode (BaseInstr): skip BaseInstr.__init__.
        >>> b._dbg_print_skip_method()
        debug mode (BaseInstr): skip <module>.
        >>> b.func = lambda : b._dbg_print_skip_method()
        >>> b.func()
        debug mode (BaseInstr): skip <lambda>.
        """
        stk = inspect.stack()
        self._dbg_print('skip {}.'.format(stk[1][3]))

    def q(self, query_str, chkerr=False):
        """
        VISA query.

        >>> b = BaseInstr(None)
        debug mode (BaseInstr): skip BaseInstr.__init__.

        >>> b.q('*IDN?')
        debug mode (BaseInstr): q: return empty string.
        ''

        :param query_str:
        :param chkerr:
        :return: debug mode -> ''
        :rtype: str
        """

        if self._debug_mode:
            self._dbg_print('q: return empty string.')
            res = ''
        else:
            res = self._rsrc.query(query_str)

        if chkerr:
            self.check_error()
        return res

    def w(self, write_str, chkerr=False):
        """
        Write.
        Do nothing if debug mode.

        >>> b = BaseInstr(None)
        debug mode (BaseInstr): skip BaseInstr.__init__.
        >>> b.w('*RST', True)
        debug mode (BaseInstr): skip w.
        debug mode (BaseInstr): skip check_error.
        """
        if self._debug_mode:
            self._dbg_print_skip_method()
        else:
            self._rsrc.write(write_str)

        if chkerr:
            self.check_error()

    # Abstract methods ---------------------------------------------------------
    def check_error(self):
        """
        Call super if debug mode.

        >>> b = BaseInstr(None)
        debug mode (BaseInstr): skip BaseInstr.__init__.

        >>> b.check_error()
        debug mode (BaseInstr): skip check_error.
        """
        if self._debug_mode:
            self._dbg_print_skip_method()
            return
        else:
            raise NotImplementedError

    def reset(self):
        """
        Call super if debug mode.
        check_error() after reset.
        """
        if self._debug_mode:
            self._dbg_print_skip_method()
            self.check_error()
        else:
            raise NotImplementedError()

class SourceMeter(BaseInstr):
    @staticmethod
    def _v_points_to_step(v_start, v_end, v_points):
        """
        :return: v_step (Always positive value.)
        :rtype: float

        >>> s = SourceMeter(None)
        debug mode (SourceMeter): skip BaseInstr.__init__.
        >>> s._v_points_to_step(0, 1, 11)
        0.1
        >>> s._v_points_to_step(-1, -2, 11)
        0.1
        """
        return abs(v_end - v_start)/(v_points - 1)

    @staticmethod
    def _v_step_to_points(v_start, v_end, v_step):
        """
        >>> s = SourceMeter(None)
        debug mode (SourceMeter): skip BaseInstr.__init__.
        >>> s._v_step_to_points(0, 1, 0.1)
        11
        >>> s._v_step_to_points(0, 1, 0.101)
        10
        >>> s._v_step_to_points(-1, -2, 0.1)
        11

        :param v_step: Must be positive.
        """
        if v_step < 0:
            raise ValueError('v_step must be positive.')
        return floor(abs(v_end - v_start)/v_step) + 1

    # Abstract methods ---------------------------------------------------------
    def iv_sweep(self, v_start, v_end, v_step=1e-3, v_points=None,
                 i_limit=1e-3):
        """
        Call super().iv_sweep if debug mode.
        debug mode: Simulate 1kOhm resistance with noise ~1nA.

        >>> s = SourceMeter(None)
        debug mode (SourceMeter): skip BaseInstr.__init__.

        >>> np.random.seed(1)
        >>> s.iv_sweep(1e-3, 3e-3, 1e-3)  # doctest: +NORMALIZE_WHITESPACE
        debug mode (SourceMeter): iv_sweep: Return dummy data.
        (array([[  1.00000000e-03,   1.00162435e-06],
               [  2.00000000e-03,   2.00162435e-06],
               [  3.00000000e-03,   3.00162435e-06]]), False)

        >>> np.random.seed(1)
        >>> s.iv_sweep(1, -1, 1)  # doctest: +NORMALIZE_WHITESPACE
        debug mode (SourceMeter): iv_sweep: Return dummy data.
        (array([[  1.00000000e+00,   1.00000162e-03],
               [  0.00000000e+00,   1.62434536e-09],
               [ -1.00000000e+00,  -9.99998376e-04]]), False)

        >>> np.random.seed(1)
        >>> s.iv_sweep(0, 1, 2)
        debug mode (SourceMeter): iv_sweep: Return dummy data.
        (array([[  0.00000000e+00,   1.62434536e-09]]), False)

        :param i_limit: positive. (absolute value)
        :type v_start: float
        :type v_end: float
        :type v_step: float
        :param v_points: overrides v_step
        :type v_points: int
        :return: [(v0, i0), (v1, i1), ...], is_aborted
        :rtype: np.ndarray, bool
        """
        if v_points is None:
            v_points = self._v_step_to_points(v_start, v_end, v_step)
        if self._debug_mode:
            self._dbg_print('iv_sweep: Return dummy data.')
            vs = np.linspace(v_start, v_end, v_points)
            Is = vs / 1000 + np.random.normal(0, 1e-9)
            vis = np.array([vs, Is]).transpose()
            return vis, False
        else:
            raise NotImplementedError()

    def iv_sweep_double(self, v_max, v_step=1e-3, v_points=None,
                        i_limit=1e-3):
        """
        Call super if debug mode.
        Sweep 0V -> v_max -> 0V

        >>> s = SourceMeter(None)
        debug mode (SourceMeter): skip BaseInstr.__init__.

        >>> np.random.seed(1)
        >>> s.iv_sweep_double(0.1, 0.1)  # doctest: +NORMALIZE_WHITESPACE
        debug mode (SourceMeter): iv_sweep: Return dummy data.
        debug mode (SourceMeter): iv_sweep: Return dummy data.
        (array([[  0.00000000e+00,   1.62434536e-09],
               [  1.00000000e-01,   1.00001624e-04],
               [  1.00000000e-01,   9.99993882e-05],
               [  0.00000000e+00,  -6.11756414e-10]]), False)

        >>> np.random.seed(1)
        >>> s.iv_sweep_double(0.1, v_points=2)
        debug mode (SourceMeter): iv_sweep: Return dummy data.
        debug mode (SourceMeter): iv_sweep: Return dummy data.
        (array([[  0.00000000e+00,   1.62434536e-09],
               [  1.00000000e-01,   1.00001624e-04],
               [  1.00000000e-01,   9.99993882e-05],
               [  0.00000000e+00,  -6.11756414e-10]]), False)

        :param i_limit: positive. (even if negative sweep)
        :param v_max: can be negative.
        :param v_points: overrides v_step
        :return: [(v0, i0), (v1, i1), ... (v0, i2n)], is_aborted
        :rtype: (np.ndarray, bool)
        """
        # Implementation example
        vis1, aborted = self.iv_sweep(0, v_max, v_step=v_step,
                                      v_points=v_points, i_limit=i_limit)
        if aborted:
            return vis1, aborted
        vis2, aborted = self.iv_sweep(v_max, 0, v_step=v_step,
                                      v_points=v_points, i_limit=i_limit)
        vis = np.concatenate((vis1, vis2))
        return vis, aborted

class TemperatureController(BaseInstr):
    def read_temp(self):
        raise NotImplementedError

if __name__ == '__main__':
    import doctest
    doctest.testmod()

max_stars_repo_name: wataash/SUSS

max_issues_repo_name: wataash/instr

max_forks_repo_name: wataash/instr


text:
import numpy as np
import xarray
from brainio_collection.packaging import package_data_assembly
import h5py
from brainio_base.assemblies import DataAssembly

DATA_DIR = '/braintree/data2/active/users/tmarques/bs_datasets/FreemanZiemba2013_V1V2data.mat'
ASSEMBLY_NAME = 'movshon.FreemanZiemba2013_V1_properties'
TEXTURE_STIM_NAME = 'movshon.FreemanZiemba2013_properties'

PROPERTY_NAMES = ['texture_modulation_index', 'absolute_texture_modulation_index', 'texture_selectivity',
                  'noise_selectivity', 'texture_sparseness', 'noise_sparseness', 'variance_ratio', 'sample_variance',
                  'family_variance', 'max_texture', 'max_noise', ]

def collect_data(data_dir):
    response_file = h5py.File(data_dir, 'r')

    texture_modulation_index, absolute_texture_modulation_index, texture_selectivity, noise_selectivity, \
    texture_sparseness, noise_sparseness, variance_ratio, sample_variance, family_variance, max_texture, max_noise = \
        calculate_texture_properties(response_file, area='v1')

    # Bins
    max_texture_bins = np.logspace(-1, 3, 13, base=10)
    max_noise_bins = np.logspace(-1, 3, 13, base=10)
    texture_modulation_index_bins = np.linspace(-1, 1, num=25)
    absolute_texture_modulation_index_bins = np.linspace(0, 1, num=13)
    texture_selectivity_bins = np.linspace(0, 1, num=11)
    noise_selectivity_bins = np.linspace(0, 1, num=11)
    texture_sparseness_bins = np.linspace(0, 1, num=11)
    noise_sparseness_bins = np.linspace(0, 1, num=11)
    variance_ratio_bins = np.logspace(-2, 1, num=13)
    sample_variance_bins = np.linspace(0, 1, num=11)
    family_variance_bins = np.linspace(0, 1, num=11)

    # Create DataAssembly with single neuronal properties and bin information
    assembly = np.concatenate((texture_modulation_index, absolute_texture_modulation_index, texture_selectivity,
                               noise_selectivity, texture_sparseness, noise_sparseness, variance_ratio, sample_variance,
                               family_variance, max_texture, max_noise), axis=1)

    assembly = DataAssembly(assembly, coords={'neuroid_id': ('neuroid', range(assembly.shape[0])),
                                              'region': ('neuroid', ['V1'] * assembly.shape[0]),
                                              'neuronal_property': PROPERTY_NAMES},
                            dims=['neuroid', 'neuronal_property'])

    assembly.attrs['number_of_trials'] = 20

    for p in assembly.coords['neuronal_property'].values:
        assembly.attrs[p+'_bins'] = eval(p+'_bins')

    return assembly

def main():
    assembly = collect_data(DATA_DIR)
    assembly.name = ASSEMBLY_NAME
    print('Packaging assembly')
    package_data_assembly(xarray.DataArray(assembly), assembly_identifier=assembly.name, stimulus_set_identifier=TEXTURE_STIM_NAME,
                          assembly_class='PropertyAssembly', bucket_name='brainio.contrib')

def calculate_texture_properties(response_file, area='v1'):
    textureNumOrder = np.array([327, 336, 393, 402, 13, 18, 23, 30, 38, 48, 52, 56, 60, 71, 99])
    familyNumOrder = np.array([60, 56, 13, 48, 71, 18, 327, 336, 402, 38, 23, 52, 99, 393, 30])
    textureSorted = [np.argwhere(textureNumOrder == familyNumOrder[i])[0][0] for i in range(len(familyNumOrder))]

    responses = response_file[area].value
    # (cellNum) x (timeBin) x (rep) x (sample) x (texType) x (texFamily)
    # (102+)    x (300)     x (20)  x (15)     x (2)       x (15)

    responses = responses[:, :, :, :, :, textureSorted]
    responses = np.moveaxis(responses, 4, 2)  # texType to position 2
    responses = np.moveaxis(responses, 5, 3)  # family to position 3
    responses = np.moveaxis(responses, 5, 4)  # sample to position 4
    # (cellNum) x (timeBin) x (texType) x (texFamily) x (sample) x (rep)

    n_neurons, n_tb, n_type, n_fam, n_smp, n_rep = responses.shape

    responses_stabilized = np.copy(responses)  # copies data for treating as stabilized spike counts

    # Calculates average responses in 10ms timebins (in spike/s)
    responses = responses.reshape(n_neurons, 30, 10, n_type, n_fam, n_smp, n_rep).mean(axis=2) * 1000

    n_tb = responses.shape[1]
    sample_responses = responses.mean(axis=5)  # average across repetitions
    family_responses = sample_responses.mean(axis=4)  # average across samples
    type_responses = family_responses.mean(axis=3)  # average across types
    mean_responses = type_responses.mean(axis=2)  # average across all
    baseline = mean_responses.min(axis=1, keepdims=True)  # baseline response

    latency = calculate_latency(family_responses)  # gets latency using sample averages
    sample_responses = calculate_mean_response(sample_responses, lat=latency)  # mean response during 100ms
    sample_responses = sample_responses - baseline.reshape(-1, 1, 1, 1) # subtracts baseline response and rectifies
    sample_responses[sample_responses < 0] = 0
    family_responses = sample_responses.mean(axis=3)
    max_response = sample_responses.reshape((n_neurons, n_type, -1)).max(axis=2)

    # Sums spikes in 10ms timebins
    responses_stabilized = responses_stabilized.reshape(n_neurons, 30, 10, n_type, n_fam, n_smp, n_rep).sum(axis=2)
    responses_stabilized = calculate_mean_response(responses_stabilized, lat=latency) * 10
    responses_stabilized = np.sqrt(responses_stabilized) + np.sqrt(responses_stabilized + 1)
    responses_stabilized = responses_stabilized.mean(axis=4)

    # Maximum texture response
    max_texture = max_response[:, 1].reshape((-1, 1))

    # Maximum noise response
    max_noise = max_response[:, 0].reshape((-1, 1))

    texture_modulation_index = np.zeros((n_neurons, 1))
    texture_selectivity = np.zeros((n_neurons, 1))
    noise_selectivity = np.zeros((n_neurons, 1))
    texture_sparseness = np.zeros((n_neurons, 1))
    noise_sparseness = np.zeros((n_neurons, 1))
    variance_ratio = np.zeros((n_neurons, 1))
    sample_variance = np.zeros((n_neurons, 1))
    family_variance = np.zeros((n_neurons, 1))

    for n in range(n_neurons):
        # Texture modulation
        texture_modulation_index[n] = calculate_texture_modulation(family_responses[n])[0]
        # Texture selectivity (sparseness over texture families)
        texture_selectivity[n] = calculate_sparseness(family_responses[n, 1])
        # Noise selectivity (sparseness over noise families)
        noise_selectivity[n] = calculate_sparseness(family_responses[n, 0])
        # Texture sparseness (sparseness over texture samples)
        texture_sparseness[n] = calculate_sparseness(sample_responses[n, 1])
        # Noise sparseness (sparseness over noise samples)
        noise_sparseness[n] = calculate_sparseness(sample_responses[n, 0])
        variance_ratio[n], sample_variance[n], family_variance[n] = calculate_variance_ratio(responses_stabilized[n, 1])

    absolute_texture_modulation_index = np.abs(texture_modulation_index)

    return texture_modulation_index, absolute_texture_modulation_index, texture_selectivity, noise_selectivity,\
           texture_sparseness, noise_sparseness, variance_ratio, sample_variance, family_variance, max_texture, \
           max_noise

def calculate_mean_response(response, lat, resp_nbin=10):

    response_shape = response.shape
    n_neur = response_shape[0]
    mean_response = np.zeros(np.concatenate((np.array([n_neur]), np.array(response_shape[2:]))))

    for n in range(n_neur):
        mean_response[n] = response[n, lat[n]:(lat[n] + resp_nbin + 1)].mean(axis=0)

    return mean_response

def calculate_latency(response, t_pk=[7, 22], t_lat_min=2, max_d_lat=10, thrsh=0.15):

    # n_neur, n_tb, n_type, n_fam = meanByFam.shape

    response_shape = response.shape
    n_neur = response_shape[0]

    lat = np.zeros(n_neur).astype(int)
    for n in range(n_neur):
        pref_stim = np.unravel_index(np.argmax(response[n, t_pk[0]:t_pk[1]]), response_shape[1:])
        pk_lat = pref_stim[0]
        if len(pref_stim)==2:
            pref_response = response[n, :, pref_stim[1]]
        elif len(pref_stim)==3:
            pref_response = response[n, :, pref_stim[1], pref_stim[2]]
        elif len(pref_stim) == 4:
            pref_response = response[n, :, pref_stim[1], pref_stim[2], pref_stim[3]]

        pref_response_pk = np.max(pref_response)
        pk_lat = pk_lat + t_pk[0]

        min_lat = np.max([t_lat_min, pk_lat-max_d_lat])
        min_val = np.min(pref_response[min_lat:pk_lat])

        lat_min = np.argmin(pref_response[min_lat:pk_lat]) + min_lat
        range_val = pref_response_pk - min_val

        lat[n] = np.argwhere(pref_response[lat_min:(pk_lat + 1)] <= (min_val + thrsh * range_val))[-1][0] \
                 + lat_min + 1
    return lat

def calculate_texture_modulation(response):
    texMod_fam = (response[1, :] - response[0, :]) / (response[1, :] + response[0, :])
    texMod = np.nanmean(texMod_fam)

    return texMod, texMod_fam

def calculate_sparseness(response):

    response = response.reshape(-1)
    n_stim = response.shape[0]

    spars = (1 - ((response.sum() / n_stim) ** 2) / ((response ** 2).sum() / n_stim)) / (1 - 1 / n_stim)

    return spars

def calculate_variance_ratio(response):

    residual_ms, sample_ms, family_ms = calculate_variance(response)
    response_shape = response.shape

    if len(response_shape) == 3:
        residual_variance = residual_ms
        sample_variance = (sample_ms - residual_ms) / response_shape[2]
        family_variance = (family_ms - sample_ms) / (response_shape[2]*response_shape[1])

    else:
        residual_variance = 0
        sample_variance = sample_ms
        family_variance = (family_ms - sample_ms) / response_shape[1]

    total_variance = residual_variance + sample_variance + family_variance

    variance_ratio = (family_variance / total_variance + 0.02) / (sample_variance / total_variance + 0.02)

    return variance_ratio, sample_variance / total_variance, family_variance / total_variance

def calculate_variance(response):

    response_shape = response.shape

    if len(response_shape) == 3:
        a, b, n = response_shape

        sample_mean = response.mean(axis=2)
        family_mean = sample_mean.mean(axis=1)
        all_mean = family_mean.mean()

        residual_ms = np.sum((response - sample_mean.reshape(a, b, 1)) ** 2) / (a * b * (n - 1))
        sample_ms = n * np.sum((sample_mean - family_mean.reshape(a, 1)) ** 2) / (a * (b - 1))
        family_ms = b*n*np.sum((family_mean - all_mean) ** 2) / (a - 1)
    else:
        a, b = response_shape

        sample_mean = response
        family_mean = sample_mean.mean(axis=1)
        all_mean = family_mean.mean()

        residual_ms = np.nan
        sample_ms = np.sum((sample_mean - family_mean.reshape(a, 1)) ** 2) / (a * (b - 1))
        family_ms = b * np.sum((family_mean - all_mean) ** 2) / (a - 1)

    return residual_ms, sample_ms, family_ms

if __name__ == '__main__':
    main()

ext: py

max_stars_repo_name: dmayo/brain-score

max_issues_repo_name: dmayo/brain-score

max_forks_repo_name: dmayo/brain-score

reason: import numpy


text:
import numpy as np

from neupy import layers
from neupy.utils import asfloat

from base import BaseTestCase

class SliceLayerConnectionsTestCase(BaseTestCase):
    def test_change_output_layer(self):
        network = layers.join(
            layers.Input(10, name='input-1'),
            layers.Relu(5, name='relu-1'),
            layers.Relu(1, name='relu-2'),
        )

        self.assertEqual(network.input_shape, (10,))
        self.assertEqual(network.output_shape, (1,))
        self.assertEqual(len(network), 3)

        relu_1_network = network.end('relu-1')
        self.assertEqual(relu_1_network.input_shape, (10,))
        self.assertEqual(relu_1_network.output_shape, (5,))
        self.assertEqual(len(relu_1_network.layers), 2)

        predict = relu_1_network.compile()
        x_test = asfloat(np.ones((7, 10)))
        y_predicted = predict(x_test)
        self.assertEqual(y_predicted.shape, (7, 5))

    def test_select_network_branch(self):
        network = layers.join(
            layers.Input(10, name='input-1'),
            [[
                layers.Relu(1, name='relu-1'),
            ], [
                layers.Relu(2, name='relu-2'),
            ]]
        )

        self.assertEqual(network.input_shape, (10,))
        self.assertEqual(network.output_shape, [(1,), (2,)])
        self.assertEqual(len(network), 3)

        relu_1_network = network.end('relu-1')
        self.assertEqual(relu_1_network.input_shape, (10,))
        self.assertEqual(relu_1_network.output_shape, (1,))
        self.assertEqual(len(relu_1_network), 2)

        predict = relu_1_network.compile()
        x_test = asfloat(np.ones((7, 10)))
        y_predicted = predict(x_test)
        self.assertEqual(y_predicted.shape, (7, 1))

        relu_2_network = network.end('relu-2')
        self.assertEqual(relu_2_network.input_shape, (10,))
        self.assertEqual(relu_2_network.output_shape, (2,))
        self.assertEqual(len(relu_2_network), 2)

    def test_cut_output_layers_in_sequence(self):
        network = layers.join(
            layers.Input(10, name='input-1'),
            layers.Relu(5, name='relu-1'),
            layers.Relu(1, name='relu-2'),
        )

        self.assertEqual(network.input_shape, (10,))
        self.assertEqual(network.output_shape, (1,))
        self.assertEqual(len(network), 3)

        cutted_network = network.end('relu-1').end('input-1')
        self.assertEqual(cutted_network.input_shape, (10,))
        self.assertEqual(cutted_network.output_shape, (10,))
        self.assertEqual(len(cutted_network), 1)

        predict = cutted_network.compile()
        x_test = asfloat(np.ones((7, 10)))
        y_predicted = predict(x_test)
        self.assertEqual(y_predicted.shape, (7, 10))

    def test_cut_using_layer_object(self):
        relu = layers.Relu(2)
        network = layers.Input(10) > relu > layers.Sigmoid(1)

        self.assertEqual(network.input_shape, (10,))
        self.assertEqual(network.output_shape, (1,))
        self.assertEqual(len(network), 3)

        cutted_network = network.end(relu)
        self.assertEqual(cutted_network.input_shape, (10,))
        self.assertEqual(cutted_network.output_shape, (2,))
        self.assertEqual(len(cutted_network), 2)

    def test_unknown_layer_name_exception(self):
        network = layers.join(
            layers.Input(10, name='input-1'),
            layers.Relu(5, name='relu-1'),
            layers.Relu(1, name='relu-2'),
        )
        with self.assertRaises(NameError):
            network.end('abc')

    def test_change_input_layer(self):
        network = layers.join(
            layers.Input(10, name='input-1'),
            layers.Relu(5, name='relu-1'),
            layers.Relu(1, name='relu-2'),
        )

        self.assertEqual(network.input_shape, (10,))
        self.assertEqual(network.output_shape, (1,))
        self.assertEqual(len(network), 3)

        relu_1_network = network.start('relu-1')
        self.assertEqual(relu_1_network.input_shape, (10,))
        self.assertEqual(relu_1_network.output_shape, (1,))
        self.assertEqual(len(relu_1_network), 2)

        predict = relu_1_network.compile()
        x_test = asfloat(np.ones((7, 10)))
        y_predicted = predict(x_test)
        self.assertEqual(y_predicted.shape, (7, 1))

    def test_cut_input_and_output_layers(self):
        network = layers.join(
            layers.Input(10, name='input-1'),
            layers.Relu(8, name='relu-0'),
            layers.Relu(5, name='relu-1'),
            layers.Relu(2, name='relu-2'),
            layers.Relu(1, name='relu-3'),
        )

        self.assertEqual(network.input_shape, (10,))
        self.assertEqual(network.output_shape, (1,))
        self.assertEqual(len(network), 5)

        cutted_network = network.start('relu-1').end('relu-2')

        self.assertEqual(cutted_network.input_shape, (8,))
        self.assertEqual(cutted_network.output_shape, (2,))
        self.assertEqual(len(cutted_network), 2)

        predict = cutted_network.compile()
        x_test = asfloat(np.ones((7, 8)))
        y_predicted = predict(x_test)
        self.assertEqual(y_predicted.shape, (7, 2))

    def test_cut_input_layers_in_sequence(self):
        network = layers.join(
            layers.Input(10, name='input-1'),
            layers.Relu(5, name='relu-1'),
            layers.Relu(1, name='relu-2'),
        )

        self.assertEqual(network.input_shape, (10,))
        self.assertEqual(network.output_shape, (1,))
        self.assertEqual(len(network), 3)

        cutted_network = network.start('relu-1').start('relu-2')
        self.assertEqual(cutted_network.input_shape, (5,))
        self.assertEqual(cutted_network.output_shape, (1,))
        self.assertEqual(len(cutted_network), 1)

        predict = cutted_network.compile()
        x_test = asfloat(np.ones((7, 5)))
        y_predicted = predict(x_test)
        self.assertEqual(y_predicted.shape, (7, 1))

    def test_connect_cutted_layers_to_other_layers(self):
        network = layers.join(
            layers.Input(10, name='input-1'),
            layers.Relu(8, name='relu-0'),
            layers.Relu(5, name='relu-1'),
            layers.Relu(2, name='relu-2'),
            layers.Relu(1, name='relu-3'),
        )

        self.assertEqual(network.input_shape, (10,))
        self.assertEqual(network.output_shape, (1,))
        self.assertEqual(len(network), 5)

        cutted_network = network.start('relu-1').end('relu-2')
        self.assertEqual(cutted_network.input_shape, (8,))
        self.assertEqual(cutted_network.output_shape, (2,))
        self.assertEqual(len(cutted_network), 2)

        new_network = layers.join(
            layers.Input(8),
            cutted_network,
            layers.Sigmoid(11),
        )
        self.assertEqual(new_network.input_shape, (8,))
        self.assertEqual(new_network.output_shape, (11,))
        self.assertEqual(len(new_network), 4)

        predict = network.compile()
        x_test = asfloat(np.ones((7, 10)))
        y_predicted = predict(x_test)
        self.assertEqual(y_predicted.shape, (7, 1))

        predict = new_network.compile()
        x_test = asfloat(np.ones((7, 8)))
        y_predicted = predict(x_test)
        self.assertEqual(y_predicted.shape, (7, 11))

    def test_get_layer_by_name_from_connection(self):
        network = layers.join(
            layers.Input(10, name='input-1'),
            layers.Relu(8, name='relu-0'),
            layers.Relu(5, name='relu-1'),
        )

        reul0 = network.layer('relu-0')
        self.assertEqual(reul0.output_shape, (8,))

        reul1 = network.layer('relu-1')
        self.assertEqual(reul1.output_shape, (5,))

        with self.assertRaises(NameError):
            network.layer('some-layer-name')

max_stars_repo_name: vishalbelsare/neupy

max_issues_repo_name: vishalbelsare/neupy

max_forks_repo_name: vishalbelsare/neupy


text:
import os
os.environ['CUDA_VISIBLE_DEVICES']='1'
import os
import tqdm
import argparse

import torch
from torchvision.utils import save_image
import torch.nn as nn
# from model import WaveEncoder, WaveDecoder

from utils.core import feature_wct
from utils.core import feature_adin
from utils.core import feature_adin_without_segment
from utils.core import feature_wct_without_segment
from utils.io import Timer, open_image, load_segment, compute_label_info
from xiaokemodel import XiaoKeEncoder, XiaoKeDecoder
import numpy as np
import torchvision.transforms as transforms

from scipy.io import loadmat
from PIL import Image
from scipy.misc import imread, imresize
import cv2
import datetime

IMG_EXTENSIONS = [
    '.jpg', '.JPG', '.jpeg', '.JPEG',
    '.png', '.PNG',
]

def is_image_file(filename):
    return any(filename.endswith(extension) for extension in IMG_EXTENSIONS)

class WCT2:
    def __init__(self, model_path='./model_checkpoints', transfer_at=['encoder', 'skip', 'decoder'], option_unpool='cat5', device='cuda:0', verbose=False):

        self.transfer_at = set(transfer_at)
        assert not(self.transfer_at - set(['encoder', 'decoder', 'skip'])), 'invalid transfer_at: {}'.format(transfer_at)
        assert self.transfer_at, 'empty transfer_at'
        model_path = './xiaoke_video_checkpoints/'
        encoder_path = 'xiaoke_encoder.pth'
        decoder_path = 'xiaoke_decoder_0.0001_4.pth'

        model_path = './xiaoke_checkpoints/'
        encoder_path = 'xiaoke_encoder.pth'
        decoder_path = 'xiaoke_decoder_87.pth'
        self.device = torch.device(device)
        self.verbose = verbose

        self.encoder = XiaoKeEncoder(option_unpool).to(self.device)
        self.decoder = XiaoKeDecoder(option_unpool).to(self.device)
        self.encoder.load_state_dict(torch.load(os.path.join(model_path,encoder_path),map_location=lambda storage, loc: storage))
        self.decoder.load_state_dict(torch.load(os.path.join(model_path,decoder_path),map_location=lambda storage, loc: storage))

    def print_(self, msg):
        if self.verbose:
            print(msg)

    def encode(self, x, skips, level):
        return self.encoder.encode(x, skips, level)

    def decode(self, x, skips, level):
        return self.decoder.decode(x, skips, level)

    def get_all_feature(self, x):
        skips = {}
        feats = {'encoder': {}, 'decoder': {}}
        for level in [1, 2, 3, 4]:
            x = self.encode(x, skips, level)
            if 'encoder' in self.transfer_at:
                feats['encoder'][level] = x

        if 'encoder' not in self.transfer_at:
            feats['decoder'][4] = x
        for level in [4, 3, 2]:
            x = self.decode(x, skips, level)
            if 'decoder' in self.transfer_at:
                feats['decoder'][level - 1] = x
        return feats, skips

    def transfer(self, content):
        content_feat, content_skips = content, {}

        wct2_enc_level = [1, 2, 3, 4]
        wct2_dec_level = [1, 2, 3, 4]
        wct2_skip_level = ['pool1', 'pool2', 'pool3']
        label_set,label_indicator = None, None
        for level in [1, 2, 3, 4]:
            content_feat = self.encode(content_feat, content_skips, level)
        for level in [4, 3, 2, 1]:
            content_feat = self.decode(content_feat, content_skips, level)
        return content_feat

def get_all_transfer():
    ret = []
    for e in ['encoder']:
        for d in ['decoder']:
            for s in ['skip']:
                _ret = set([e, d, s]) & set(['encoder', 'decoder', 'skip'])
                if _ret:
                    ret.append(_ret)
    return ret

def run_bulk(in_path,dir_name):
    accurate_segment = True
    device = 'cpu' if config.cpu or not torch.cuda.is_available() else 'cuda:0'
    device = torch.device(device)

    transfer_at = set()
    if config.transfer_at_encoder:
        transfer_at.add('encoder')
    if config.transfer_at_decoder:
        transfer_at.add('decoder')
    if config.transfer_at_skip:
        transfer_at.add('skip')
    # cw, ch =  640,360
    cw, ch =  640,400

    # The filenames of the content and style pair should match
    c_transforms = transforms.Compose([transforms.Resize((ch,cw), interpolation=Image.NEAREST),transforms.CenterCrop((ch // 16 * 16, cw // 16 * 16)),transforms.ToTensor()])

    fnames = os.listdir(in_path)
    fnames.sort()
    sample_fnames = fnames[4:8]
    for fname in tqdm.tqdm(sample_fnames):
        if not is_image_file(fname):
            print('invalid file (is not image), ', fname)
            continue
        # content
        _content = os.path.join(in_path, fname)
        content = Image.open(_content).convert('RGB') # 别忘了这边的to(device)

        content = c_transforms(content).unsqueeze(0).to(device)
        print('current frame {} and shape is {}'.format(fname,content.shape))

        _output = os.path.join(config.output, dir_name+fname)
        print('current in_path_frame ', _content)
        print('current out_path_frame ',_output)
        if not config.transfer_all:
            with Timer('Elapsed time in whole WCT: {}', config.verbose):
                postfix = '_'.join(sorted(list(transfer_at)))
                fname_output = _output.replace('.png', '_{}_{}.png'.format(config.option_unpool, postfix))
                wct2 = WCT2(transfer_at=transfer_at, option_unpool=config.option_unpool, device=device, verbose=config.verbose)
                with torch.no_grad():
                    img = wct2.transfer(content)

                save_image(img.clamp_(0, 1), fname_output, padding=0)
        else:
            for _transfer_at in get_all_transfer():
                with Timer('Elapsed time in whole WCT: {}', config.verbose):
                    postfix = '_'.join(sorted(list(_transfer_at)))
                    fname_output = _output.replace('.png', '_{}_{}.png'.format(config.option_unpool, postfix))
                    wct2 = WCT2(transfer_at=_transfer_at, option_unpool=config.option_unpool, device=device, verbose=config.verbose)
                    with torch.no_grad():
                        starttime = datetime.datetime.now()
                        img = wct2.transfer(content)
                        endtime = datetime.datetime.now()
                        print('xiaoke with adin 运行时间为----',(endtime - starttime))
                    save_image(img.clamp_(0, 1), fname_output, padding=0)
                # break

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('--content', type=str, default='./examples/content')
    parser.add_argument('--content_segment', type=str, default='./examples/content_segment')
    parser.add_argument('--style', type=str, default='./examples/style')
    parser.add_argument('--style_segment', type=str, default='./examples/style_segment')
    parser.add_argument('--output', type=str, default='./outputs')
    parser.add_argument('--image_size', type=int, default=512)
    parser.add_argument('--alpha', type=float, default=1)
    parser.add_argument('--option_unpool', type=str, default='cat5', choices=['sum', 'cat5'])
    parser.add_argument('-e', '--transfer_at_encoder', action='store_true')
    parser.add_argument('-d', '--transfer_at_decoder', action='store_true')
    parser.add_argument('-s', '--transfer_at_skip', action='store_true')
    parser.add_argument('-a', '--transfer_all', action='store_true')
    parser.add_argument('--cpu', action='store_true')
    parser.add_argument('--verbose', action='store_true')
    parser.add_argument('--is_wct',action='store_true')
    parser.add_argument('--label_mapping', type=str, default='ade20k_semantic_rel.npy')
    parser.add_argument('--model_path', help='folder to model path', default='baseline-resnet50_dilated8-ppm_bilinear_deepsup')
    parser.add_argument('--arch_encoder', default='resnet50_dilated8', help="architecture of net_encoder")
    parser.add_argument('--arch_decoder', default='ppm_bilinear_deepsup', help="architecture of net_decoder")
    parser.add_argument('--suffix', default='_epoch_20.pth', help="which snapshot to load")
    parser.add_argument('--fc_dim', default=2048, type=int, help='number of features between encoder and decoder')
    parser.add_argument('--num_class', default=150, type=int, help='number of classes')
    parser.add_argument('--padding_constant', default=8, type=int, help='maxmimum downsampling rate of the network')
    parser.add_argument('--gpu_id', default=0, type=int, help='gpu_id for evaluation')
    parser.add_argument('--imgSize', default=[300, 400, 500, 600], nargs='+', type=int, help='list of input image sizes.' 'for multiscale testing, e.g. 300 400 500')

    # 不能定义两次同样的参数
    config = parser.parse_args()

    transform = transforms.Compose([transforms.Normalize(mean=[102.9801, 115.9465, 122.7717], std=[1., 1., 1.])])
    print(config)
    if not os.path.exists(os.path.join(config.output)):
        os.makedirs(os.path.join(config.output))
    data_path = config.content
    for in_path in os.listdir(data_path):
        run_bulk(os.path.join(data_path,in_path),in_path)

# 树林的图片
# 171124_D1_HD_01
# 170216A_122_ForestTrail_1080
# 170216A_070_LookingUpThroughForest_1080
# 180705_01_0
# 190416_10_Drone1_0
# Forest_15_1_Videv
# Forest_15_4_Videv
# on
# WalkingThroughTreesatSunsetVidev

# 树叶
# Autumn_leaves_in_motion_0
# autumn_leaves
# autumn-leaves-blowing-in-the-wind-H264
# 180705_01_0

# 海浪
# 46234354
# walking_on_the_beac

# 雪山
# 180607_A_00

# 开车
# 180607_A_10

# 飞机
# Airbus_A380_Landing_2__Videv
# Evening_landin
# PlaneLand

# 海边瑜伽
# Ao_Nang_Beach_Yoga_MP4_HDV_1080p25__TanuriX_Stock_Footage_N
# MVI_126

# 水稻
# Barley_3_Videv
# HandStrokin
# wild_gras
# windygrassnoaudi-

# 船
# beach1
# sailing_boa

# 天空
# Becco_di_Filadonna_su_Vall
# Blue_Sky_and_Clouds_Timelapse_0892__Videv

# 老鼠
# CotswoldSequence

# 奶牛
# cow
# Cow_Mother_and_cal
# Cows_
# Limousin_Cows_1__VIdev
# Limousin_Cows_2__Videv

# 日落
# Lonely_tree_at_sunset_CCBY_NatureCli
# MilkyWaywithTreeVidev
# SilhouetteJogge
# Sun_to_Sea_Model__Pan_Down_MP4_HDV_1080p25__TanuriX_Stock_Footage_W
# Sunris
# TimelapseSunse
# Wakeboarding_on_the_Lak

# 马
# Dirty_Hors

# 黑白鸟
# Pigeon-Stock-Vide
# Red_fod
# seagul-H264
# seagulls_on_the_beac
# Weave

# 建筑
# Run_5_wo_metadata_h264420_720p_UH

# 鸭子
# SeaBirdsSwimming_
# Swans__1287_

# 羊
# Shee

'''
CUDA_VISIBLE_DEVICES=6 python transfer.py --content ./examples/content --style ./examples/style --content_segment ./examples/content_segment --style_segment ./examples/style_segment/ --output ./outputs/ --verbose --image_size 512 -a
'''

'''
python xiaoketransfer.py --content ./examples/demo_content/ --style ./examples/demo_style/ -a --output ./examples/demo_stylization --is_wct --image_size 400
CUDA_VISIBLE_DEVICES=1 python xiaoketransfer.py --content ./examples/dataset/alley_2/ --style ./examples/dataset/fangao.png -a --output ./examples/stylization

CUDA_VISIBLE_DEVICES=1 python xiaoketransfer2.py --content ./examples/data/MPI-Sintel-complete/training/clean/temple_2 --style ./examples/data/fangao.png -a --output ./examples/stylization

CUDA_VISIBLE_DEVICES=1 python xiaoketransfer2.py --content ./examples/data/MPI-Sintel-complete/training/clean/mountain_1 --style ./examples/data/fangao.png -a --output ./examples/stylization

CUDA_VISIBLE_DEVICES=1 python xiaoketransfer2.py --content ./examples/data/MPI-Sintel-complete/training/clean/temple_2 --style ./examples/data/fangao.png -a --output ./examples/stylization --is_wct

'''

'''
'../data/video-picture/160825_26_WindTurbines4_1080'
python xiaoketransfer2.py --content ../data/video-picture/160825_26_WindTurbines4_1080 --style ./examples/data/fangao.png -a --output ./examples/160825_26_WindTurbines4_1080_adain

'''

'''
'../data/video-picture/xxx'
python compare_model.py --content ../data/video-picture/ --style ../data/reference/tar0056_orange_forest.png -a --output ./examples/compare_model
'''

max_forks_repo_name: xuzhongyou/Colorization

reason: import numpy,from scipy


text:
#!/usr/bin/env python
import PyDSTool
from time import clock
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pickle
from matplotlib.ticker import FormatStrFormatter,ScalarFormatter

'''This script runs the ODE system of the third in silico experiment
A comparison with the averaged traces of the MCell simulations is done afterwards.
'''

params = {'axes.labelsize': 7,
           'axes.titlesize': 6,
          'legend.fontsize': 6,
           'xtick.labelsize': 6,
           'ytick.labelsize': 6,
            'figure.figsize': (3.2,2.9)}
            #'figure.figsize': (6.4, 6)}
mpl.rcParams.update(params)

Na = 6.02214e23
#vmito = 0.033221e-15 #lt
#vcyto = 0.232549e-15 #
vmito = 0.020706e-15 #lt OM
vcyto = 0.0181e-15 # IM
vcube = 0.0923e-15

no1_ant = 20e3
k7 = 800.0 #0.74
k8 = 0.7 #0.05
k9 = 0.58 #0.37
k10 = 0.48#0.47
no_atp = 3800 #number of atphase
cm = 3e3 #nr.of particles = 12mM #adp+atp = cm
a12 = 24
a21 = 40.0
a23 = 4.0
a32 = 5e3
ac = 9806 #clamp*Na*vcyto
kp  = 1.0
n_porin = 7100

DSargs = PyDSTool.args(name='ex')

DSargs.pars = { 'k1_on':'4.0',#          #KDo
                'k1_off':'100', #         #KDo
                'k2_on': '6.4',#         #KTi
		'k2_off':'40000.0',#     #KTi
		'k5_on':'0.4',           #KTo
		'k5_off':'200.0',           #KTo
		'k6_on':'4.0',           #KDi
		'k6_off':'40000.0',      #KDi
		'k7':k7,           #kp
		'k8':k8,              #kcp
		'k9':k9,             #kt
		'k10':k10,            #kd
		'no_ant': no1_ant,
		'a12':'24', #s-1
		'a21':'40.0',#s-1
		'a65':'924',           #s-1
		'a56':'1e3',             #s-1,
		'a16':'452457',           #s-1
		'a61':'11006',           #s-1
		'a25':'1.17e-12',     #s-1
		'a52':'2.0',          #s-1
		'a54':'1e2',          #s-1
		'a45':'1e2',          #s-1
		'a43':'0.8',         #uM-1s-1
		'a34':'1e2',
		'a32':'5e3',          #s-1
		'a23':'4.0',         #uM-1s-1
		'no1_atp': no_atp, #number of atphase
		'vcyto': vcyto,
		'vcube':vcube,#syn -vmito
		'vmito':vmito,
		'kp':kp,
		'n_porin': n_porin,
		'n_cha':'7e5',
		'k_cha':'1.0',
		#'k_cha_off':'0',
        'k_cha_off':'2e-2',
		'Na':Na,
		'fa':'0.5',
		'ac':ac,
		't_on1':'0.001',
		't_off1':'0.005',
		't_on2':'0.010',
		't_off2':'0.014'}

#DSargs.fnspecs  = {'alb': (['l','al','lb','bl','bla','la','ala','blb','alap','blbp'], 'no_ant-l-al-la-lb-bl-bla-ala-blb-alap-blbp') }
DSargs.fnspecs = {'k_cha_t': (['t'], "(heav(t>=t_on1)*heav(t<t_off1)+heav(t>t_on2)*heav(t<t_off2))*(k_cha-k_cha_off)+k_cha_off")}
#DSargs.fnspecs = {'k_cha_t': (['t'], "(heav(t>=t_on1))*k_cha_off")}

DSargs.varspecs = { 'am':'-k6_on*(1e6*am/(Na*vmito))*bl + k6_off*bla - k6_on*(1e6*am/(Na*vmito))*l + k6_off*la - 2.0*k6_on*fa*(1e6*am/(Na*vmito))*al + k6_off*ala + k6_off*alap + a34*h3es - a43*(1e6*am/(Na*vmito))*h3eo',
                    'bm':'-k2_on*(1e6*bm/(Na*vmito))*al + k2_off*(no_ant-l-al-la-lb-bl-bla-ala-blb-alap-blbp) - k2_on*(1e6*bm/(Na*vmito))*l + k2_off*lb - 2.0*k2_on*fa*(1e6*bm/(Na*vmito))*bl + k2_off*blb + k2_off*blbp -a23*h3e*(1e6*bm/(Na*vmito)) + a32*h3es',
                    'bo':'-kp*(1e6*bo/(Na*vcube))*n_porin + kp*(1e6*bc/(Na*vcyto))*n_porin - n_cha*k_cha_t(t)*(1e6*bo/(Na*vcube))',
                    'bc':'-k5_on*(1e6*bc/(Na*vcyto))*l + k5_off*bl - k5_on*(1e6*bc/(Na*vcyto))*la + k5_off*bla- 2.0*k5_on*fa*(1e6*bc/(Na*vcyto))*lb + k5_off*blb + k5_off*blbp - kp*(1e6*bc/(Na*vcyto))*n_porin + kp*(1e6*bo/(Na*vcube))*n_porin',
                    'l':'-k1_on*(1e6*ac/(Na*vcyto))*l + k1_off*al- k2_on*(1e6*bm/(Na*vmito))*l + k2_off*lb - k5_on*(1e6*bc/(Na*vcyto))*l + k5_off*bl - k6_on*(1e6*am/(Na*vmito))*l + k6_off*la',
                    'al':'k1_on*(1e6*ac/(Na*vcyto))*l - k1_off*al - k2_on*(1e6*bm/(Na*vmito))*al + k2_off*(no_ant-l-al-la-lb-bl-bla-ala-blb-alap-blbp) - k6_on*fa*(1e6*am/(Na*vmito))*al + k6_off*ala - k6_on*fa*(1e6*am/(Na*vmito))*al + k6_off*alap',
                    'la':'k6_on*(1e6*am/(Na*vmito))*l - k6_off*la  - k5_on*(1e6*bc/(Na*vcyto))*la + k5_off*bla - k1_on*fa*(1e6*ac/(Na*vcyto))*la + k1_off*ala - k1_on*fa*(1e6*ac/(Na*vcyto))*la + k1_off*alap',
                    'lb':'k2_on*(1e6*bm/(Na*vmito))*l - k2_off*lb - k1_on*(1e6*ac/(Na*vcyto))*lb + k1_off*(no_ant-l-al-la-lb-bl-bla-ala-blb-alap-blbp) - k5_on*fa*(1e6*bc/(Na*vcyto))*lb + k5_off*blb - k5_on*fa*(1e6*bc/(Na*vcyto))*lb + k5_off*blbp',
                    'bl':'k5_on*(1e6*bc/(Na*vcyto))*l - k5_off*bl- k6_on*(1e6*am/(Na*vmito))*bl + k6_off*bla - k2_on*fa*(1e6*bm/(Na*vmito))*bl + k2_off*blb- k2_on*fa*(1e6*bm/(Na*vmito))*bl + k2_off*blbp',
                    'bla':'k6_on*(1e6*am/(Na*vmito))*bl - k6_off*bla - k8*bla + k7*(no_ant-l-al-la-lb-bl-bla-ala-blb-alap-blbp) + k5_on*(1e6*bc/(Na*vcyto))*la - k5_off*bla',
                     'ala':'k1_on*fa*(1e6*ac/(Na*vcyto))*la - k1_off*ala + k6_on*fa*(1e6*am/(Na*vmito))*al - k6_off*ala -k10*ala + k10*alap',
                     'alap':'k1_on*fa*(1e6*ac/(Na*vcyto))*la - k1_off*alap + k6_on*fa*(1e6*am/(Na*vmito))*al - k6_off*alap - k10*alap + k10*ala',
                     'blb':'k2_on*fa*(1e6*bm/(Na*vmito))*bl - k2_off*blb + k5_on*fa*(1e6*bc/(Na*vcyto))*lb - k5_off*blb + k9*blbp - k9*blb',
                     'blbp':'k2_on*fa*(1e6*bm/(Na*vmito))*bl - k2_off*blbp + k5_on*fa*(1e6*bc/(Na*vcyto))*lb - k5_off*blbp -k9*blbp +k9*blb',
                     'eo':'-a65*eo + a56*(no1_atp-eo-ei-h3eo-h3es-h3e)+ a16*ei - a61*eo',
                     'ei':'-a16*ei + a61*eo - a12*ei + a21*h3e',
                     'h3eo':'-a45*h3eo + a54*(no1_atp-eo-ei-h3eo-h3es-h3e) + a34*h3es - a43*h3eo*(1e6*am/(Na*vmito))',
                   'h3es':'a43*(1e6*am/(Na*vmito))*h3eo - a34*h3es + a23*h3e*(1e6*(bm)/(Na*vmito)) - a32*h3es',
                   'h3e':'-a23*h3e*(1e6*(bm)/(Na*vmito)) + a32*h3es - a25*h3e + a52*(no1_atp-eo-ei-h3eo-h3es-h3e) + a12*ei - a21*h3e'}

# initial conditions
DSargs.ics  = {'am':37739,
	           'bm':500,
	           'bc':200,
	           'bo':0,
                'l':500,
                'al':14.5e3,
                'lb':0,
                'bl':300,
                'bla':100,
                'la':140,
                'ala':2.2e3,
                'blb':0,
                'alap':2.2e3,
                'blbp':0,
                'eo':1000,
                'ei':20,
                'h3eo':33,
                'h3es':80,
                'h3e':1940}
DSargs.tdomain = [0,0.014]

#ode  = PyDSTool.Generator.Vode_ODEsystem(DSargs)    # an instance of the 'Generator' class.
ode  = PyDSTool.Generator.Radau_ODEsystem(DSargs)
start = clock()
traj = ode.compute('polarization')
print ('  ... finished in %.3f seconds.\n' % (clock()-start))
pd   = traj.sample()

#------------average traces are loaded -----------------------

with open('./av_10c_1e9', 'rb') as f: av_ar = pickle.load(f)

with open('./av_10r_1e9', 'rb') as f: av_ra = pickle.load(f)

with open('./av_10rc_1e9', 'rb') as f: av_cu = pickle.load(f)
time = np.loadtxt('./time')

#---------

plt.figure(1)

ax = plt.subplot(4,2,2)
plt.plot(1e3*time-2.5,av_ar[0,:], 'b')
plt.plot(1e3*time-2.5,av_ra[0,:], 'red')
plt.plot(1e3*time-2.5,av_cu[0,:], 'm')
plt.plot(1e3*pd['t']-2.5, ac*np.ones(len(pd['t'])),c='k')
ax.locator_params(axis='y',nbins=3)
plt.ylabel(r'#ADP$_{\rm outside}$')
ax.locator_params(axis='x',nbins=8)
ax.yaxis.set_major_formatter(ScalarFormatter(useMathText=True)) #con 10 como base
plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0)) #notacion cientifica
ax.set_yticks([0.94e4,0.97e4,1e4])
#plt.text(1,29650,'A', ha='center', va='center', fontsize=11)
#plt.ylim(28000,30000)
plt.xlim(0,10)
#plt.arrow(1, 2.8e4, 0, 150,lw= 5, fc='k', ec='k')
#ax.xaxis.set_major_formatter(plt.NullFormatter())

ax2 = plt.subplot(4,2,4)
plt.plot(1e3*time-2.5,av_ar[1,:],'b')
plt.plot(1e3*time-2.5,av_ra[1,:],'red')
plt.plot(1e3*time-2.5,av_cu[1,:],'m')
plt.plot(1e3*pd['t']-2.5, pd['am'],c='k')
ax2.locator_params(axis='y',nbins=3)
ax2.locator_params(axis='x',nbins=8)
ax2.yaxis.set_major_formatter(ScalarFormatter(useMathText=True))
plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
#plt.text(1,68300,'C', ha='center', va='center', fontsize=11)
#plt.arrow(1, 67500, 0, 72,lw= 5, fc='k', ec='k')
plt.xlim(0,10)
ax2.set_yticks([3.75e4,3.80e4,3.85e4])
plt.ylabel(r'#ADP$_{\rm matrix}$')
#ax2.xaxis.set_major_formatter(plt.NullFormatter())

ax3 = plt.subplot(4,2,5)
plt.plot(1e3*time-2.5,av_ar[2,:], 'b')
plt.plot(1e3*time-2.5,av_ra[2,:], 'red')
plt.plot(1e3*time-2.5,av_cu[2,:], 'm')
plt.plot(1e3*pd['t']-2.5,pd['bm'],c='k')
#plt.text(1,1235,'D', ha='center', va='center', fontsize=11)
ax3.locator_params(axis='y',nbins=2)
ax3.yaxis.set_major_formatter(ScalarFormatter(useMathText=True))
plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
plt.xlim(0,10)
ax3.set_yticks([4.5e2,5e2,5.5e2])
plt.ylabel(r'#ATP$_{\rm matrix}$')
#plt.xlabel('Time (msec)')
#ax3.xaxis.set_major_formatter(plt.NullFormatter())
ax3.set_xticks([0,2,4,6,8,10])

ax1 = plt.subplot(4,2,6)

plt.plot(1e3*time-2.5,av_ar[3,:], 'b')
plt.plot(1e3*time-2.5,av_ra[3,:], 'red')
plt.plot(1e3*time-2.5,av_cu[3,:], 'm')
plt.plot(1e3*pd['t']-2.5, pd['bc'],c='k')
ax1.locator_params(axis='y',nbins=3)
plt.ylabel(r'#ATP$_{\rm outside}$')
ax1.yaxis.set_major_formatter(ScalarFormatter(useMathText=True))
plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
#ax1.set_yticks([0e3,3e3,6e3])
#plt.text(1,650,'B', ha='center', va='center', fontsize=11)
plt.xlim(0, 10)
#ax1.xaxis.set_major_formatter(plt.NullFormatter())
ax1.set_yticks([0e2,2e2,4e2])
ax1.set_xticks([0,2,4,6,8,10])
#plt.xlabel('Time (msec)')

ax4 = plt.subplot(4,1,4)

plt.plot(1e3*time-2.5,av_ar[29,:],'b')
plt.plot(1e3*time-2.5,av_ra[29,:],'red')
plt.plot(1e3*time-2.5,av_cu[29,:],'m')
plt.plot(1e3*pd['t']-2.5,pd['bo'],'k')
#plt.text(0.5,120,'E', ha='center', va='center', fontsize=11)#
#plt.xlabel('Time (msec)')
plt.ylabel(r'#ATP$_{\rm syn}$')
plt.subplots_adjust(wspace =0.5)
ax4.locator_params(axis='y',nbins=2)
#ax4.xaxis.set_major_formatter(plt.NullFormatter())
ax4.set_yticks([0,1.1e2,2.2e2])
ax4.set_xticks([0,2,4,6,8,10])
plt.ylim(0,2.2e2)
plt.xlim(0,10)

plt.tight_layout(pad=0.05, w_pad=0.001)
#plt.savefig('syn2_av_10_3x2n.pdf', transparent=True, format = 'pdf',dpi=600,bbox_inches='tight')

plt.show()


text:
import json
import click
import numpy as np
import jenks_natural_breaks

@click.command()
@click.argument('filename', type=click.Path(exists=True))
@click.argument('nr-classes', type=int)
@click.option('--attribute')
def cli(filename, nr_classes, attribute):
    result = {}
    with open(filename, 'r') as myfile:
        geojson = json.loads(myfile.read())
        layername = geojson["name"]
        aantal_list = [ft["properties"]["aantal"] for ft in geojson["features"]]
        np_array = np.array(aantal_list)
        breaks = jenks_natural_breaks.classify(np_array, 5)
        breaks[0] = 0
        result[layername] = {}
        result[layername][attribute] = breaks
    print(json.dumps(result))

if __name__=="__main__":
    cli()

max_stars_repo_name: arbakker/corona-kaart

max_issues_repo_name: arbakker/corona-kaart

max_forks_repo_name: arbakker/corona-map-nl


text:
# -*- coding: utf-8 -*-
"""
Created on Fri Jun  1 11:28:20 2018

@author: Shih-Yang Lin

This script calculates the counterfactual market share when absentee voting
is allowed in all markets.
"""

import numpy as np
import pandas as pd
from BLP_market_share import calculate_n_consumer_market_share as cal_mkt_shr
import matplotlib.pyplot as plt
import sys

# Read data
# =============================================================================

if sys.platform == 'darwin':
    file_path = "/Users/errard/Dropbox/SSQ_Election/RawData20180515.csv"
    market_population_file_name = "/Users/errard/Dropbox/SSQ_Election/adj_census_pop_T.csv"
    xi_file_name = "/Users/errard/Dropbox/SSQ_Election/xi.csv"
else:
    file_path = "C:/SSQ_Election/RawData20180515.csv"
    market_population_file_name = "C:/SSQ_Election/adj_census_pop_T.csv"
    xi_file_name = "C:/SSQ_Election/xi.csv"

def read_csv(file_path):
    '''
    This function read csv files to pandas dataframe.

    input:
        file_path: str.
    output:
        dta: pandas dataframe
    '''
    f = open(file_path)
    dta = pd.read_csv(f)
    f.close()
    return(dta)

dta = read_csv(file_path)
dta.loc[:, 'const'] = np.ones((len(dta), 1))

market_population = read_csv(market_population_file_name)
market_population = market_population.drop(columns = ['depcity', 'arrcity'])

xi = np.loadtxt(xi_file_name)
dta = dta.assign(xi = xi)

dta = pd.merge(dta, market_population,
               left_on = ['depcity', 'arrcity'],
               right_on = ['depcity_code', 'arrcity_code'],
               how = 'left')

# Set up hyper parameters
# =============================================================================
n_consumers = int(1e4)

estimates = [-5.246357, -.6285018, -.8258803, -1.604927, 1.448219,
              0.7957358,  1.093028, -0.2235259, -0.1370856, -6.126835, 1]
sigma = [2.816481]

data = dta[['mkt', 'priced', 'timed', 'pop_dep', 'moninc', 'edu',
            'unemploymentrate_arr', 'unemploymentrate_dep',
            'poll', 'mayor_arr', 'const', 'xi', 'census_pop', 'brand',
            'marketshare1']].copy()

end_var = [0]

seed = 654781324

# =============================================================================
# Create new product -- absentee voting

abs_vot = dta[['mkt', 'priced', 'timed', 'pop_dep', 'moninc', 'edu',
            'unemploymentrate_arr', 'unemploymentrate_dep',
            'poll', 'mayor_arr', 'const', 'xi', 'census_pop',
            'brand', 'marketshare1']].groupby(['mkt']).agg({
               'priced': 'mean', 'timed': 'mean', 'pop_dep': 'mean',
               'moninc': 'mean', 'edu': 'mean', 'unemploymentrate_arr': 'mean',
               'unemploymentrate_dep': 'mean', 'poll': 'mean', 'mayor_arr': 'mean',
               'const': 'mean', 'xi': 'median', 'census_pop': 'mean', 'brand': 'mean',
               'marketshare1': 'mean'})

abs_vot.loc[:, 'priced'] = 0
abs_vot.loc[:, 'timed'] = 0
abs_vot.loc[:, 'brand'] = 0
abs_vot.loc[:, 'marketshare1'] = 0
abs_vot = abs_vot.reset_index()

# Add convenient voting options to each markets
data = data.append(abs_vot, ignore_index = True)

mkt = data['mkt'].copy()

#%% Calculating counterfactual market share.
data.loc[:, 'c_mkt_shr'] = cal_mkt_shr(estimates,
            sigma, data.drop(columns = ['mkt', 'census_pop', 'brand', 'marketshare1']),
            end_var, mkt, n_consumers, seed = seed)

#%%
total_pop = sum(market_population['census_pop'])
observed_voters = sum(dta['marketshare1'] * dta['census_pop'])
counterfactual_n_voters = sum(data['c_mkt_shr']* data['census_pop'])

#%% Calculate traffic volume change for voting
n_voters_difference = round(counterfactual_n_voters - observed_voters)
turnout_difference = round(n_voters_difference*100/total_pop, 2)
observed_turnout = round(observed_voters*100/total_pop, 2)
counterfactual_turnout = observed_turnout + turnout_difference
print('')
print('=========================================================')
print('If convenience voting is allowed:')
print('The number of voters will increase by:')
print(str(n_voters_difference))

print('The turnout rate will increase by:')
print(str(turnout_difference) + '%')

print('The turnout rate before convenience voting is allowed:')
print(str(observed_turnout) + '%')
print('The turnout rate after convenience voting is allowed:')
print(str(counterfactual_turnout) + '%')
print('=========================================================')

#%%
data['b_brand'] = 0
data.loc[(data['brand'] >= 201) & (data['brand'] <= 203), 'b_brand'] = 2
data.loc[(data['brand'] > 203) & (data['brand'] <= 303), 'b_brand'] = 3
data.loc[data['brand'] == 401, 'b_brand'] = 4
data.loc[data['brand'] == 502, 'b_brand'] = 5

#%%
data = data.assign(predicted_voters = round(data['c_mkt_shr'] * data['census_pop']),
                   observed_voters = round(data['marketshare1'] * data['census_pop']))
#%% Produce Table xxx
summary_1 = data.loc[:, ['b_brand', 'predicted_voters', 'observed_voters']].groupby(by = 'b_brand').sum()
summary_1['difference'] = summary_1['predicted_voters'] - summary_1['observed_voters']
summary_1 = summary_1.assign(mkt_shr_before = round(summary_1['observed_voters']*100/total_pop, 2),
                             mkt_shr_after = round(summary_1['predicted_voters']*100/total_pop, 2))

print('The number of non-voters before convenience voting is allowed is:')
print(str(total_pop - sum(summary_1['observed_voters'])))

print('The number of non-voters after convenience voting is allowed is:')
print(str(total_pop - sum(summary_1['predicted_voters'])))

#%% Produce Table xxx

city_location_dict = {1: 'N', 2: 'N', 3: 'N', 8: 'N',
                      4: 'C', 11: 'C', 12: 'C', 13: 'C', 14: 'C',
                      5: 'S', 6: 'S', 9: 'S', 16: 'S'}
dep_region = data['mkt'].copy()
arr_region = data['mkt'].copy()

for i in range(len(dep_region)):
    if len(str(dep_region[i])) == 3:
        dep_region[i] = str(dep_region[i])[0]
        arr_region[i] = str(arr_region[i])[1:3]
    else:
        dep_region[i] = str(dep_region[i])[0:2]
        arr_region[i] = str(arr_region[i])[2:4]

dep_region = dep_region.replace(city_location_dict)
arr_region = arr_region.replace(city_location_dict)

data = data.assign(dep_region = dep_region, arr_region = arr_region)

summary_2 = data.loc[:, ['dep_region', 'arr_region',
                         'predicted_voters', 'observed_voters']].groupby(by = ['dep_region', 'arr_region']).agg({
                         'predicted_voters': 'sum', 'observed_voters': 'sum'})

census_pop_by_region = data.loc[:, ['mkt', 'dep_region', 'arr_region', 'census_pop']].copy()
census_pop_by_region = census_pop_by_region.groupby(by = ['mkt', 'dep_region', 'arr_region']).mean().reset_index()
census_pop_by_region = census_pop_by_region.loc[:, ['dep_region', 'arr_region', 'census_pop']].groupby(by = ['dep_region', 'arr_region']).sum()

summary_2['census_pop'] = census_pop_by_region['census_pop']

summary_2['difference'] = summary_2['predicted_voters'] - summary_2['observed_voters']
summary_2['difference%'] = round(summary_2['difference']*100/summary_2['observed_voters'], 2)

summary_2['turnout_after'] = round(summary_2['predicted_voters']*100/summary_2['census_pop'], 2)
summary_2['turnout_before'] = round(summary_2['observed_voters']*100/summary_2['census_pop'], 2)
summary_2['turnout_diff'] = summary_2['turnout_after'] - summary_2['turnout_before']


text:
import numpy as np
from scipy import optimize

def findgyrostuff(array):
    """
    Find the Larmor radius and gyrocenter for a particle executing many revolutions.
    Assumes the circle will be well filled.
    Parameters
    ----------
    array: Array particle positions in (x,xp,y,yp,z,zp)

    Returns
    -------
    (grho, gcen): Tuple with Larmor radius and gyrocenter

    """
    xmin = np.argmin(array[:, 0])
    xmax = np.argmax(array[:, 0])

    ymin = np.argmin(array[:, 2])
    ymax = np.argmax(array[:, 2])

    grho = (abs(array[ymax, 2] - array[ymin, 2]) + abs(array[xmax, 0] - array[xmin, 0])) / 4.
    gcen = ((array[ymax, 0] + array[ymin, 0]) / 2., (array[xmax, 2] + array[xmin, 2]) / 2.)

    return grho, gcen

def _findR(x, y, xc, yc):
    return np.sqrt((x-xc)**2 + (y-yc)**2)

def _f(c, x, y):
    R = _findR(x, y, *c)
    return R - np.average(R)

def opt_circle(x, y):
    """
    Preferred method for calculating Larmor radius and gyrocenter.
     Least squares fit to coordinates of a circle.
     Based on example from Scipy cookbook
    Parameters
    ----------
    x: Particle x coordinates
    y: Particle y coordinates

    Returns
    -------
    xc: Coordinate in x of circle center
    yc: Coordinate in y of circle center
    R: Circle radius
    residual: Residual of the fit (sum(Ri - R)**2)

    """
    x_m = np.average(x)
    y_m = np.average(y)
    ce = x_m, y_m
    center, ier = optimize.leastsq(_f, ce, args=(x, y))
    xc, yc = center
    Ri = _findR(x, y, *center)
    R = np.average(Ri)
    residual = np.sum((Ri - R)**2)
    return xc, yc, R, residual

def nrGR(B, v):
    """
    Calculate the nonrelativistic Larmor radius.
    Parameters
    ----------
    B: Magnetic field in Tesla
    v: Transverse velocity in m/s

    Returns
    -------
    r: Larmor radius in m

    """
    me = 9.10938356e-31

    p = me * v
    e = 1.6e-19
    r = p / (e * B)
    return r

def relGR(B, v, vl):
    """
    Calculate the relativistically corrected Larmor radius
    Parameters
    ----------
    B: Magnetic field in Tesla
    v: Transverse velocity in m/s
    vl: Total velocity in m/s

    Returns
    -------

    """
    c = 299792458
    me = 9.10938356e-31

    vtot = np.sqrt(v**2 + vl**2)
    gamma = 1/np.sqrt(1-vtot**2/c**2)
    p = gamma * me * v
    r = p / (1.6e-19 * B)
    return r

max_stars_repo_name: radiasoft/mcool

max_issues_repo_name: radiasoft/mcool

max_forks_repo_name: radiasoft/mcool


text:
from keras.backend.tensorflow_backend import set_session
from keras.preprocessing.image import load_img
from keras_vggface.utils import preprocess_input
from keras_vggface.vggface import VGGFace
from scipy.spatial.distance import cosine
import tensorflow as tf
from PIL import Image
import numpy as np
import os

# Set environment variable to make sure that TF does not eat up all GPU memory
os.environ['TF_FORCE_GPU_ALLOW_GROWTH'] = 'true'

class VGGFaceRecognizer:
    """
    Easy Face Recognizer Class, built on top of RC Mali's VGG Face Implementation
    URL (Keras VGG Face): https://github.com/rcmalli/keras-vggface
    License: MIT

    Author: Animikh Aich
    GitHub: https://github.com/animikhaich
    LinkedIn: https://www.linkedin.com/in/animikh-aich/
    """

    def __init__(self, model='senet50'):
        """
        __init__ Model Initializer

        - Sets the flags for GPU memory growth (if available)
        - Selects and Initializes the Model

        Args:
            model (str, optional): Choose Base Feature Extractor: VGG16, ResNet50, and SeNet50. Defaults to 'senet50'.
        """

        # Allow Memory Growth for Tensorflow GPU
        if tf.test.is_gpu_available():
            # For TF 1x
            if int(tf.__version__.split('.')[0]) == 1:
                config = tf.ConfigProto()
                config.gpu_options.allow_growth = True
                session = tf.Session(config=config)
                set_session(session)

            # For TF 2x
            else:
                gpu_devices = tf.config.experimental.list_physical_devices(
                    'GPU'
                )
                tf.config.experimental.set_memory_growth(gpu_devices[0], True)

        # Some Instance Variables
        self.model = model
        self.face_dict = dict()

        # initialize VGG Face using VGGFace Model
        self.recognizer = VGGFace(include_top=False, model=model)

    def calculate_similarity(self, vector_1, vector_2):
        """
        calculate_similarity

        Calculate the Cosine Similarity of Similarity Feature Vectors
        Flattens the Vectors for VGG16 base due to higher high dimension vector

        Args:
            vector_1 (numpy array): Feature Vector of Face 1
            vector_2 (numpy array): Feature Vector of Face 2

        Returns:
            float: Similarity Score, Lower means more similar
        """
        # VGG creates higher dimensional vector, and hence needs to be flattened
        if self.model == 'vgg16':
            vector_1 = vector_1.flatten()
            vector_2 = vector_2.flatten()

        # Remove Extra Dimension
        vector_1 = np.squeeze(vector_1)
        vector_2 = np.squeeze(vector_2)

        # Calculate the cosine similarity
        return cosine(vector_1, vector_2)

    def register_faces(self, face_list, name_list):
        """
        register_faces

        Given a list of faces and names, register them and save them to the database
        These are the faces that a query face will be compared against.
        Every time this is called, it clear the existing database

        Args:
            face_list (list): List of PIL Face Images
            name_list (list): List of Strings containing the name of each face - Unique

        Returns:
            dict: Dictionary containing the name of the face and it's corresponding feature vector
        """
        self.face_dict = dict()
        for face, name in zip(face_list, name_list):
            self.face_dict[name] = self.feature_extractor(face)

        return self.face_dict

    def add_faces(self, face_list, name_list):
        """
        add_faces

        Given a list of faces and names, register them and save them to the database
        These are the faces that a query face will be compared against.
        Every time this is called, it appends to existing faces in the DB

        Args:
            face_list (list): List of PIL Face Images
            name_list (list): List of Strings containing the name of each face - Unique

        Returns:
            dict: Dictionary containing the name of the face and it's corresponding feature vector
        """
        for face, name in zip(face_list, name_list):
            self.face_dict[name] = self.feature_extractor(face)

        return self.face_dict

    def feature_extractor(self, face):
        """
        feature_extractor

        Primary CNN based feature extractor for faces. Uses pre-trained model by RC Mali for the same.
        You can choose the architecture while creating the instance of the class

        Args:
            face (PIL Image): PIL image of a cropped out face

        Returns:
            list: Returns the feature vector that is extracted from the given face
        """
        # Resize the face, face is a PIL image
        face = face.resize((224, 224), Image.ANTIALIAS)
        face = np.asarray(face).astype(np.float)
        face = np.expand_dims(face, axis=0)

        if self.model == 'vgg16':
            face = preprocess_input(face, version=1)
        else:
            face = preprocess_input(face, version=2)

        return self.recognizer.predict(face)

    def recognize(self, face, thresh=0.25):
        """
        recognize

        Given a Query Face, this function compares it against the existing faces and returns a match, if any

        Args:
            face (PIL Image): Cropped face in the form of a PIL image
            thresh (float, optional): Cutoff Threshold used to discard unmatched faces. Defaults to 0.25.

        Returns:
            string: Name of the matched face, if recognized, else None
        """
        # Extract features of the Query Face
        query_features = self.feature_extractor(face)
        temp_sim_dict = dict()

        # Compare the Query Face features with the existing database faces
        for key in self.face_dict.keys():
            db_face_features = self.face_dict[key]
            score = self.calculate_similarity(
                db_face_features, query_features
            )
            temp_sim_dict[key] = score

        try:
            # Return None if no face matches
            if min(temp_sim_dict.values()) > thresh:
                return None
        except:
            return None

        # Find most similar face
        most_similar_face = min(temp_sim_dict, key=temp_sim_dict.get)

        return most_similar_face


text:
import numpy as np

def contains_nan(array):
    """
    Efficiently checks if a numpy array contains a NaN value.

    cf. http://stackoverflow.com/questions/6736590/fast-check-for-nan-in-numpy

    A NaN values may indicate a calibration error. The source finder doesn't
    know how to interpret a NaN value so we reject such an image.

    args:
        array: a Numpy array

    Returns:
        Error string if array contains NaN, False otherwise
    """
    # summing an array containing NaN will result in NaN
    if np.isnan(np.sum(array)):
        return "Image data contains NaN value"
    else:
        return False

max_stars_repo_name: ajstewart/tkp

max_issues_repo_name: ajstewart/tkp

max_forks_repo_name: ajstewart/tkp


text:
import pandas as pd
import numpy as np
import lightgbm as lgb
from sklearn.model_selection import StratifiedKFold
from datetime import datetime
import os
import logging

def train_lgb(X, y, params, folder, rewrite=False, fold=5):
    os.makedirs(folder, exist_ok=rewrite)
    spl = StratifiedKFold(n_splits=fold, random_state=2019, shuffle=True)
    feature_importance_df = pd.DataFrame()
    counter = 0
    for trn_ids, val_ids in spl.split(X, y):
        counter += 1
        print('fold {}'.format(counter))
        X_trn, y_trn = X.iloc[trn_ids], y[trn_ids]
        X_val, y_val = X.iloc[val_ids], y[val_ids]
        lgb_trn = lgb.Dataset(X_trn, y_trn)
        lgb_val = lgb.Dataset(X_val, y_val)
        model = lgb.train(params=params, train_set=lgb_trn,verbose_eval=500,
                          valid_sets=[lgb_trn,lgb_val],valid_names=['train','valid'])
        model.save_model(folder + '/{}_{}'.format(counter, model.best_score['valid']))
        fold_importance_df = pd.DataFrame()
        fold_importance_df["feature"] = X_trn.columns
        fold_importance_df["importance"] = model.feature_importance()
        fold_importance_df["fold"] = counter
        feature_importance_df = pd.concat([feature_importance_df, fold_importance_df], axis=0)
    feature_importance_df = feature_importance_df.groupby('feature', as_index=False)['importance'].mean().sort_values(by='importance', ascending=False)
    feature_importance_df.to_csv(os.path.join(folder, 'feature_importance.csv'), index=False)
    print('done')

def lgb_kfold(train_df,test_df,features,target,cat_features,folds,params):
    oof_preds = np.zeros(train_df.shape[0])
    sub_preds = np.zeros(test_df.shape[0])

    cv_list = []
    for n_fold, (train_idx, valid_idx) in enumerate(folds.split(train_df[features], train_df['salary_round'])):
        print ('FOLD:' + str(n_fold))

        train_x, train_y = train_df[features].iloc[train_idx], train_df[target].iloc[train_idx]
        valid_x, valid_y = train_df[features].iloc[valid_idx], train_df[target].iloc[valid_idx]

        print ('train_x shape:',train_x.shape)
        print ('valid_x shape:',valid_x.shape)

        dtrain = lgb.Dataset(train_x, label=train_y,categorical_feature=cat_features)
        dval = lgb.Dataset(valid_x, label=valid_y, reference=dtrain,categorical_feature=cat_features)
        bst = lgb.train(params, dtrain, num_boost_round=50000,
            valid_sets=[dval,dtrain], verbose_eval=500,early_stopping_rounds=500, )
        new_list = sorted(zip(features, bst.feature_importance('gain')),key=lambda x: x[1], reverse=True)[:]
        for item in new_list:
            print (item)

        oof_preds[valid_idx] = bst.predict(valid_x, num_iteration=bst.best_iteration)
        oof_cv = mean_absolute_error(valid_y,  oof_preds[valid_idx])
        cv_list.append(oof_cv)
        print (cv_list)
        sub_preds += bst.predict(test_df[features], num_iteration=bst.best_iteration) / folds.n_splits

    cv = mean_absolute_error(train_df[target],  oof_preds)
    print('Full OOF MAE %.6f' % cv)

    train_df['lgb_y'] = oof_preds
    test_df['lgb_y'] = sub_preds

    return train_df,test_df


text:
#standard imports
from typing import Tuple,Dict
#scientific imports
import numpy as np
#project imports
from astropy.stats import sigma_clip

from fitter.fit_functions import scipyFit,gaussian
from plotter.plot_handler import plot_sigma_clipping,plot_interpolation,plot_noise_residual
from res.conf_file_str import internal_noise_value,analysis_obs_time_value,internal_mag_value,analysis_upper_mag_limit,internal_multiple_mag,internal_mag_noise
from data_handler.signal_features import compute_periodogram,noise,get_mag

def get_magnitude(data : np.ndarray,kwargs : Dict,mag_check : bool = True ) -> Dict:
    """
    Changes magnitude according to Panda et al. 2018 if noise is added
    :param kwargs: run config
    :return: run config with magnitude changed
    """
    f_data = compute_periodogram(data,kwargs)
    noise_val = noise(f_data)
    mag = get_mag(noise_val)

    upper_limit_mag = 14.5 if analysis_upper_mag_limit not in kwargs else kwargs[analysis_upper_mag_limit]
    np.savetxt("computed_magnitude.txt",np.array([mag]))

    if mag > upper_limit_mag and mag_check:
        raise ValueError(f"Magnitude for noise {kwargs[internal_noise_value]} too high! Computed mag: {'%.2f' % mag}, upper limit: {'%.2f' % upper_limit_mag}")

    kwargs[internal_mag_value] = mag

    return kwargs

def refine_data(data : np.ndarray, kwargs : Dict) -> Tuple[np.ndarray,Dict]:
    """
    Refines the dataset and returns it
    :param data: Dataset that needs refining
    :param kwargs: Run configuration
    :return: Good dataset
    """
    data = normalize_data(data)
    data = reduce_data_to_observation_time(data,kwargs)
    data,kwargs = add_noise(data,kwargs)
    data = remove_stray(data,kwargs)
    data = interpolate(data,kwargs)
    return data,kwargs

def get_total_gap(data : np.ndarray) -> float:
    """
    Computes the total gap in time units for a given dataset
    :param data: datset of the lightcurve
    :return: total gap in units of time
    """
    values,counts,most_common = get_diff_values_counts_most_common(data)

    values[values - most_common < 10**-5] = 0
    return float(np.sum(values*counts))

def compute_duty_cycle(data : np.ndarray) -> float:
    """
    Computes the duty cycle for a given lightcurve.
    :param data: dataset of the lightcurve
    :return: duty cycle in percent
    """
    x = data[0]
    total_length = x[-1] - x[0]

    total_gap = get_total_gap(data)

    return (total_length-total_gap)*100/total_length

def compute_observation_time(data : np.ndarray) -> float:
    """
    Computes the total observation time for a given dataset
    :param data: dataset of the lightcurve
    :return: total observation time
    """
    x = data[0]
    total_length = x[-1] - x[0]
    total_gap = get_total_gap(data)
    return total_length - total_gap

def reduce_data_to_observation_time(data : np.ndarray, kwargs : Dict) -> np.ndarray:
    """
    Reduces the observation to a given time by cutting in the edges
    :param data: dataset of the lightcurve
    :param time: target observation time
    :return: array with given length
    """
    if analysis_obs_time_value in kwargs.keys():
        time = kwargs[analysis_obs_time_value]
        null_obs_time = compute_observation_time(data)
        if time > null_obs_time or null_obs_time < 0:
            raise ValueError(f"Target observation time cannot be reached. Needs to be bigger zero and smaller than "
                             f"zero observation time. 0 time is {null_obs_time}, target obs time is {time}")

        x = data[0]
        y = data[1]
        obs_time = compute_observation_time(np.array((x, y)))
        fringe = 0
        _,_,most_common = get_diff_values_counts_most_common(data)
        while (obs_time > time):
            mask = np.logical_and(data[0] > data[0][int(len(data[0])/2)] - time/2 - fringe/2
                                  ,data[0]< data[0][int(len(data[0])/2)] + time/2 + fringe/2)
            y = data[1][mask]
            x = data[0][mask]
            obs_time = compute_observation_time(np.array((x, y)))
            fringe = obs_time - time

        return normalize_data(np.array((x,y)))
    else:
        return data

def reduce_data_to_duty_cycle(data : np.ndarray, target_duty_cycle : float) -> np.ndarray:
    """
    Adds gaps to the dataset to reach a given duty cycle
    :param data: dataset of the lightcurve
    :param target_duty_cycle: duty cycle that is needed to reach. Needs to be less then the starting duty cycle
    :return: data with given dutycycle
    """

    null_duty_cycle = compute_duty_cycle(data)
    if target_duty_cycle > null_duty_cycle or target_duty_cycle > 100 or target_duty_cycle < 0:
        raise ValueError(f"Target duty cycle can't be reached. It needs to be between 0 and a 100 and smaller than"
                         "the starting duty cycle. Starting duty cycle is {'.%2f'%null_duty_cycle}% "
                         "and target duty cycle is {'.%2f'%target_duty_cycle}%")

    x = data[0]
    y = data[1]
    duty_cycle =compute_duty_cycle(np.array((x,y)))
    fringe = 0
    while(duty_cycle > target_duty_cycle):
        target_indices = len(data[0])*((null_duty_cycle - target_duty_cycle)/100 + fringe)
        mask = np.logical_or(data[0] < data[0][int(len(data[0])/2 - target_indices/2)]
                             ,data[0] >  data[0][int(len(data[0])/2 + target_indices/2)])
        y = data[1][mask]
        x = data[0][mask]
        duty_cycle = compute_duty_cycle(np.array((x, y)))
        fringe = duty_cycle - target_duty_cycle

    return np.array((x,y))

def normalize_data(data:np.ndarray) -> np.ndarray:
    """
    Subtracts the zero point of the time array and removes nans and infs
    :param data: Dataset
    :return: zeroed in dataset
    """
    x = data[0]
    y = data[1]

    for i in [np.inf,-np.inf,np.nan]:
        if i in y:
            n = len(y[y==i])

        x = x[y != i]
        y = y[y != i]

    return np.array(((x - x[0]), y))

def get_diff_values_counts_most_common(data : np.ndarray) -> Tuple[np.ndarray,np.ndarray,float]:
    """
    Returns the difference values and counts for a given dataset in the temporal axis as well as the most common value
    :param data: Dataset of the lightcurve
    :return: Tuple of 2 arrays, consisting of values and counts as well as a float with the most common value
    """
    x = data[0]
    diff = x[1:] - x[:-1]
    (values, counts) = np.unique(diff, return_counts=True)
    most_common = values[np.argmax(counts)]
    return values,counts,most_common

def get_gaps(data:np.ndarray) -> Tuple[int,float]:
    """
    Finds gaps in a given dataset
    :param data: Dataset
    :return: Gap ids and most common difference between two values
    """
    x = data[0]
    #approximation of difference
    diff = np.round(x[1:] - x[0:-1],decimals=2)
    values, counts, most_common = get_diff_values_counts_most_common(data)

    gap_ids = np.where(diff > 3)
    if len(gap_ids[0]) == 0:
        gap_ids = None
    else:
        gap_ids = gap_ids[0][0]

    return gap_ids,most_common

def remove_stray(data:np.ndarray,kwargs : Dict) -> np.ndarray:
    """
    Removes stray values from the dataset (values bigger than 5 sigma of the distribution of datapoints)
    :param data: dataset
    :param kwargs: Run configuration
    :return: dataset with removed strays
    """

    clipped = sigma_clip(data[1],sigma=5,maxiters=10)
    plot_sigma_clipping(data,clipped.mask,kwargs)
    data = np.array([data[0][~clipped.mask],data[1][~clipped.mask]])
    return data

def get_line(x1,x2,y1,y2,most_common):
    """
    Creates a line between two points
    :param x1: x1 point
    :param x2: x2 point
    :param y1: y1 pint
    :param y2: y2 point
    :param most_common: most common value in dataset
    :return:
    """
    b = (y1 - y2)/(x1 -x2)
    a = y1 - b*x1

    delta_t = x2 -x1
    n = int(delta_t//most_common)
    x = np.linspace(x1+most_common,x2,n,endpoint=False)
    y = a+b*x
    return np.array((x,y))

def interpolate(data : np.ndarray,kwargs : Dict)->np.ndarray:
    """
    Interpolates the dataset within the gaps
    :param data: dataset
    :return: interpolated dataset
    """
    gap_arr_ids = []
    while True:
        gap,most_common = get_gaps(data)

        if gap is None or not gap.size:
            break

        line_data = get_line(data[0][gap],data[0][gap+1],data[1][gap],data[1][gap+1],most_common)
        x = np.hstack((data[0][:gap],line_data[0],data[0][gap+1:]))
        y = np.hstack((data[1][:gap], line_data[1], data[1][gap + 1:]))
        ids = np.arange(0,len(line_data[0])) + gap

        data = np.array((x,y))
        gap_arr_ids = gap_arr_ids + ids.tolist()

    plot_interpolation(data,gap_arr_ids,kwargs)

    return data

def add_noise(data : np.ndarray, kwargs : Dict) -> Tuple[np.ndarray,Dict]:
    """
    Adds noise to the signal configured within the conf dict
    :param data: Dataset
    :param kwargs: Run conf
    :return: Noisier signal
    """
    if internal_noise_value in kwargs.keys() or (internal_multiple_mag in kwargs.keys() and kwargs[internal_multiple_mag]):
        range_data = np.amax(data[1]) - np.amin(data[1])
        binsize = int(range_data * np.exp(-range_data / 9500000))  # empirically brought down to a reasonable value

        bins = np.linspace(np.amin(data[1]), np.amax(data[1]), binsize)
        hist = np.histogram(data[1], bins=bins, density=True)[0]

        mids = (bins[1:] + bins[:-1]) / 2
        cen = np.average(mids, weights=hist)
        wid = np.sqrt(np.average((mids - cen) ** 2, weights=hist))

        p0 = [0, cen, wid]
        bins = bins[:-1]

        popt, __ = scipyFit(bins, hist, gaussian, p0)

        (cen, wid) = (popt[1], popt[2])

        x = data[0]
        y = data[1]

        if internal_noise_value in kwargs.keys():
            noise = float(kwargs[internal_noise_value]) * (np.random.normal(cen,100,len(y)))
            plot_noise_residual(data,np.array((x,(y + noise))),kwargs)
            kwargs = get_magnitude(np.array((x,(y + noise))), kwargs)
        else:
            old_mag = kwargs[internal_mag_value]
            kwargs = get_magnitude(np.array((x, y)), kwargs, False)
            noise_val = 0
            noise = 0
            while old_mag - kwargs[internal_mag_value] > 0:
                noise_val +=0.1
                noise = noise_val * (np.random.normal(cen, 100, len(y)))
                kwargs = get_magnitude(np.array((x, (y+noise))), kwargs,False)

            kwargs[internal_mag_value] = old_mag
            kwargs[internal_mag_noise] = noise_val
            plot_noise_residual(data, np.array((x, (y + noise))), kwargs)

        return np.array((x,(y + noise))),kwargs
    else:
        return data,kwargs

max_stars_repo_name: MarcoMuellner/APOLLO

max_issues_repo_name: MarcoMuellner/APOLLO

max_forks_repo_name: MarcoMuellner/APOLLO

reason: import numpy,import scipy,from astropy


text:
"""Copyright © 2020-present, Swisscom (Schweiz) AG.
All rights reserved."""

from pickle import dump, load
import tensorflow as tf
import numpy as np

from tensorflow import keras
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.utils import class_weight

from codi.mlp_codi import MLPCoDi
from codi.codi_utils import generate_plots, \
                            plot_precision_threshold_tradeoff, \
                            filter_nlp_datapoints

from codi.codi_trainer import CoDiTrainer

class NLPTrainer(MLPCoDi, CoDiTrainer):
    """
    Creates a CoDiTrainer instance.

    Parameters:
        name: string
        yaml_model_path: path of yaml containing model params.
        yaml_train_path: path to yaml containing training params.

    Attributes:
        datasets: train_dataset, test_dataset, validation_dataset, unlabelled_dataset
        scaler: the scaler instance ued to scale train, val and test sets.
        train_params: parameters extracted from yaml
        class_weights: class weights computed for imbalanced classes.
    """
    def __init__(self, name='name', yaml_model_path='mlp_codi.yaml',
                 yaml_train_path='train_codi.yaml'):
        super(NLPTrainer, self).__init__(name=name, yaml_path=yaml_model_path)
        self.load_model()
        self.train_params = self.load_yaml(yaml_train_path)
        self.scaler = None
        self.class_weights = None
        self.train_dataset = None
        self.test_dataset = None
        self.validation_dataset = None
        self.unlabelled_dataset = None
        self.x_train = None
        self.y_train = None
        self.x_val = None
        self.y_val = None
        self.x_test = None
        self.y_test = None

    def train(self):
        """
        Trains the model.
        """
        es = keras.callbacks.EarlyStopping(monitor='val_loss',
                                           patience=self.train_params['patience'],
                                           mode='auto',
                                           baseline=None, restore_best_weights=True)
        self.compile(loss=self.train_params['loss'], optimizer=self.train_params['optimizer'],
                     metrics=['accuracy'])

        history = self.fit(self.train_dataset, validation_data=self.validation_dataset,
                           epochs=self.train_params['epochs'],
                           class_weight=self.class_weights, callbacks=[es])

        test_score = self.evaluate(self.test_dataset)

        self.save_weights(self.train_params['path_to_ckpt'])

        test_pred = self.predict(self.test_dataset)

        generate_plots([history.history['accuracy'], history.history['loss'],
                        history.history['val_loss'], history.history['val_accuracy']])
        plot_precision_threshold_tradeoff(test_pred, self.y_test)

        return history.history['accuracy'][-1], test_score[-1]

    def create_labelled_dataset(self, iterator):
        """
        Creates labelled dataset as tf.data.dataset object

        Args:
            iterator: A customed torch iterator.
        """
        features = self.train_params['features']

        for idx, example_batch in enumerate(iterator):
            y_batch = [example.is_correct for example in example_batch[0]]
            X_batch = np.empty((len(example_batch[0]), 0))
            for feature in features:
                feat_arr = [getattr(example, feature) for example in example_batch[0]]
                if len(np.shape(feat_arr)) == 3:
                    feat_arr = np.squeeze(feat_arr, axis=2)
                if len(np.shape(feat_arr)) == 1:
                    feat_arr = np.expand_dims(feat_arr, axis=1)
                X_batch = np.concatenate([X_batch, feat_arr], axis=1)
            if idx == 0:
                X = X_batch
                y = y_batch
            else:
                X = np.append(X, X_batch, axis=0)
                y = np.append(y, y_batch, axis=0)

        y = y.astype(np.float32)

        self.scaler = StandardScaler()

        self.x_train, self.x_test, self.y_train, self.y_test = train_test_split(X, y, test_size=0.2, random_state=42)

        self.x_train, self.x_val, self.y_train, self.y_val = train_test_split(self.x_train, self.y_train, test_size=0.1,
                                                                              random_state=42)
        self.class_weights = dict(enumerate(class_weight.compute_class_weight('balanced',
                                                                              np.unique(self.y_train),
                                                                              self.y_train)))
        self.x_train = self.scaler.fit_transform(self.x_train)
        self.x_val = self.scaler.transform(self.x_val)
        self.x_test = self.scaler.transform(self.x_test)

        self.train_dataset = tf.data.Dataset.from_tensor_slices((self.x_train, self.y_train))
        self.train_dataset = self.train_dataset.shuffle(buffer_size=512).batch(self.train_params['batch_size'])

        self.validation_dataset = tf.data.Dataset.from_tensor_slices((self.x_val, self.y_val))
        self.validation_dataset = self.validation_dataset.batch(self.train_params['test_batch_size'])

        self.test_dataset = tf.data.Dataset.from_tensor_slices((self.x_test, self.y_test))
        self.test_dataset = self.test_dataset.batch(self.train_params['test_batch_size'])

        dump(self.scaler, open(self.train_params['path_to_scaler'], 'wb'))

    def create_unlabelled_dataset(self, iterator):
        """
        Creates unlabelled dataset as tf.data.dataset object.

        Args:
            iterator: A customed torch iterator.
        """
        features = self.train_params['features']

        for example_batch, _ in iterator:
            X_batch = np.empty((len(example_batch), 0))
            for feature in features:
                feat_arr = [getattr(example, feature) for example in example_batch]
                if len(np.shape(feat_arr)) == 3:
                    feat_arr = np.squeeze(feat_arr, axis=2)
                if len(np.shape(feat_arr)) == 1:
                    feat_arr = np.expand_dims(feat_arr, axis=1)
                X_batch = np.concatenate([X_batch, feat_arr], axis=1)

            self.scaler = load(open(self.train_params['path_to_scaler'], 'rb'))

            X_batch = self.scaler.transform(X_batch)
            self.unlabelled_dataset = tf.data.Dataset.from_tensor_slices(X_batch)

            self.unlabelled_dataset = self.unlabelled_dataset.batch(len(X_batch))

            scores = self.test()
            for index, example in enumerate(example_batch):
                example.score = scores[index]

        return filter_nlp_datapoints(iterator, self.train_params['filtering'])

    def test(self):
        """
        Classifies (evaluates) new unlabelled data using the previously trained model.
        """
        self.load_weights(self.train_params['path_to_ckpt']).expect_partial()
        return self.predict(self.unlabelled_dataset)