text:
# Licensed under a 3-clause BSD style license - see LICENSE.rst
from __future__ import division, print_function
from stingray.events import EventList
from stingray.lightcurve import Lightcurve
from stingray.powerspectrum import Powerspectrum, AveragedPowerspectrum
from stingray.powerspectrum import Crossspectrum, AveragedCrossspectrum
import numpy as np
import os
from hendrics.io import load_events, save_events, save_lcurve, load_lcurve
from hendrics.io import save_data, load_data, save_pds, load_pds
from hendrics.io import HEN_FILE_EXTENSION, _split_high_precision_number
from hendrics.io import save_model, load_model, HAS_C256, HAS_NETCDF

import pytest
import glob
from astropy.modeling import models
from astropy.modeling.core import Model

try:
    FileNotFoundError
except NameError:
    FileNotFoundError = IOError

def _dummy_bad(x, z, y=0):
    return

def _dummy(x, y=0):
    return

class TestIO():
    """Real unit tests."""
    @classmethod
    def setup_class(cls):
        cls.dum = 'bubu' + HEN_FILE_EXTENSION

    def test_save_data(self):
        struct = {'a': 0.1, 'b': np.longdouble('123.4567890123456789'),
                  'c': np.longdouble([[-0.5, 3.5]]),
                  'd': 1}
        save_data(struct, self.dum)
        struct2 = load_data(self.dum)
        assert np.allclose(struct['a'], struct2['a'])
        assert np.allclose(struct['b'], struct2['b'])
        assert np.allclose(struct['c'], struct2['c'])
        assert np.allclose(struct['d'], struct2['d'])

    def test_load_and_save_events(self):
        events = EventList([0, 2, 3.], pi=[1, 2, 3], mjdref=54385.3254923845,
                           gti=np.longdouble([[-0.5, 3.5]]))
        events.energy = np.array([3., 4., 5.])
        save_events(events, self.dum)
        events2 = load_events(self.dum)
        assert np.allclose(events.time, events2.time)
        assert np.allclose(events.pi, events2.pi)
        assert np.allclose(events.mjdref, events2.mjdref)
        assert np.allclose(events.gti, events2.gti)
        assert np.allclose(events.energy, events2.energy)

    def test_load_and_save_lcurve(self):
        lcurve = Lightcurve(np.linspace(0, 10, 15), np.random.poisson(30, 15),
                            mjdref=54385.3254923845,
                            gti=np.longdouble([[-0.5, 3.5]]))
        save_lcurve(lcurve, self.dum)
        lcurve2 = load_lcurve(self.dum)
        assert np.allclose(lcurve.time, lcurve2.time)
        assert np.allclose(lcurve.counts, lcurve2.counts)
        assert np.allclose(lcurve.mjdref, lcurve2.mjdref)
        assert np.allclose(lcurve.gti, lcurve2.gti)
        assert lcurve.err_dist == lcurve2.err_dist

    def test_load_and_save_pds(self):
        pds = Powerspectrum()
        pds.freq = np.linspace(0, 10, 15)
        pds.power = np.random.poisson(30, 15)
        pds.mjdref = 54385.3254923845
        pds.gti = np.longdouble([[-0.5, 3.5]])

        save_pds(pds, self.dum)
        pds2 = load_pds(self.dum)
        assert np.allclose(pds.gti, pds2.gti)
        assert np.allclose(pds.mjdref, pds2.mjdref)
        assert pds.m == pds2.m

    def test_load_and_save_xps(self):
        lcurve1 = Lightcurve(np.linspace(0, 10, 150),
                             np.random.poisson(30, 150),
                             mjdref=54385.3254923845,
                             gti=np.longdouble([[-0.5, 3.5]]))
        lcurve2 = Lightcurve(np.linspace(0, 10, 150),
                             np.random.poisson(30, 150),
                             mjdref=54385.3254923845,
                             gti=np.longdouble([[-0.5, 3.5]]))

        xps = AveragedCrossspectrum(lcurve1, lcurve2, 1)

        save_pds(xps, self.dum)
        xps2 = load_pds(self.dum)
        assert np.allclose(xps.gti, xps2.gti)
        assert xps.m == xps2.m
        lag, lag_err = xps.time_lag()
        lag2, lag2_err = xps2.time_lag()
        assert np.allclose(lag, lag2)
        assert hasattr(xps2, 'pds1')

    def test_load_and_save_xps_no_all(self):
        lcurve1 = Lightcurve(np.linspace(0, 10, 150),
                             np.random.poisson(30, 150),
                             mjdref=54385.3254923845,
                             gti=np.longdouble([[-0.5, 3.5]]))
        lcurve2 = Lightcurve(np.linspace(0, 10, 150),
                             np.random.poisson(30, 150),
                             mjdref=54385.3254923845,
                             gti=np.longdouble([[-0.5, 3.5]]))

        xps = AveragedCrossspectrum(lcurve1, lcurve2, 1)

        outfile = 'small_xps' + HEN_FILE_EXTENSION
        save_pds(xps, outfile, save_all=False)
        xps2 = load_pds(outfile)
        assert np.allclose(xps.gti, xps2.gti)
        assert xps.m == xps2.m
        assert not hasattr(xps2, 'pds1')

    def test_load_and_save_xps_quick(self):
        lcurve1 = Lightcurve(np.linspace(0, 10, 150),
                             np.random.poisson(30, 150),
                             mjdref=54385.3254923845,
                             gti=np.longdouble([[-0.5, 3.5]]))
        lcurve2 = Lightcurve(np.linspace(0, 10, 150),
                             np.random.poisson(30, 150),
                             mjdref=54385.3254923845,
                             gti=np.longdouble([[-0.5, 3.5]]))

        xps = AveragedCrossspectrum(lcurve1, lcurve2, 1)

        save_pds(xps, self.dum)
        xps2 = load_pds(self.dum, nosub=True)
        assert np.allclose(xps.gti, xps2.gti)
        assert xps.m == xps2.m

        assert not hasattr(xps2, 'pds1')

    def test_high_precision_split1(self):
        C_I, C_F, C_l, k = \
            _split_high_precision_number("C", np.double(0.01), 8)
        assert C_I == 1
        np.testing.assert_almost_equal(C_F, 0, 6)
        assert C_l == -2
        assert k == "double"

    def test_high_precision_split2(self):
        C_I, C_F, C_l, k = \
            _split_high_precision_number("C", np.double(1.01), 8)
        assert C_I == 1
        np.testing.assert_almost_equal(C_F, np.double(0.01), 6)
        assert C_l == 0
        assert k == "double"

    def test_save_longcomplex(self):
        val = np.longcomplex(1.01 + 2.3j)
        data = {'val': val}
        save_data(data, 'bubu' + HEN_FILE_EXTENSION)
        data_out = load_data('bubu' + HEN_FILE_EXTENSION)

        assert np.allclose(data['val'], data_out['val'])

    @pytest.mark.skipif('not HAS_C256 or not HAS_NETCDF')
    def test_save_longcomplex(self):
        val = np.complex256(1.01 + 2.3j)
        data = {'val': val}
        with pytest.warns(UserWarning) as record:
            save_data(data, 'bubu' + HEN_FILE_EXTENSION)
        assert "complex256 yet unsupported" in record[0].message.args[0]
        data_out = load_data('bubu' + HEN_FILE_EXTENSION)

        assert np.allclose(data['val'], data_out['val'])

    @classmethod
    def teardown_class(cls):
        import shutil
        for dum in glob.glob('bubu*.*'):
            os.unlink(dum)
        shutil.rmtree('bubu')

class TestIOModel():
    """Real unit tests."""

    @classmethod
    def setup_class(cls):
        cls.dum = 'bubu' + HEN_FILE_EXTENSION
        cls.model = models.Gaussian1D() + models.Const1D(amplitude=2)

    def test_save_and_load_Astropy_model(self):
        save_model(self.model, 'bubu_model.p')
        b, kind, _ = load_model('bubu_model.p')
        assert kind == 'Astropy'
        assert isinstance(b, Model)
        assert np.all(self.model.parameters == b.parameters)
        assert np.all(self.model.bounds == b.bounds)
        assert np.all(self.model.fixed == b.fixed)

    def test_save_and_load_callable_model(self):
        constraints0 = {'bounds': ()}
        save_model(_dummy, 'bubu_callable.p', constraints=constraints0)
        b, kind, constraints = load_model('bubu_callable.p')
        assert kind == 'callable'
        assert callable(b)
        assert np.all(_dummy.__code__.co_argcount == b.__code__.co_argcount)
        assert np.all(_dummy.__defaults__ == b.__defaults__)
        assert np.all(constraints == constraints0)

    def test_save_callable_model_wrong(self):
        with pytest.raises(TypeError) as record:
            save_model(_dummy_bad, 'callable_bad.p')
        assert 'Accepted callable models have only' in str(record.value)
        assert not os.path.exists('callable_bad.p')

    def test_save_junk_model(self):
        a = 'g'
        with pytest.raises(TypeError) as record:
            save_model(a, 'bad.p', constraints={'bounds': ()})
        assert 'The model has to be an Astropy model or a callable' \
               in str(record.value)
        assert not os.path.exists('bad.p')

    def test_load_python_model_callable(self):
        modelstring = '''
def model(x, a=2, b=4):
    return x * a + b

constraints = {'fixed': {'a': True}}
'''
        with open('bubu__model__.py', 'w') as fobj:
            print(modelstring, file=fobj)
        b, kind, constraints = load_model('bubu__model__.py')
        assert kind == 'callable'
        assert callable(b)
        assert b.__code__.co_argcount == 3
        assert b.__defaults__[0] == 2
        assert b.__defaults__[1] == 4
        assert np.all(constraints == {'fixed': {'a': True}})

    def test_load_python_model_Astropy(self):
        modelstring = '''
from astropy.modeling import models
model = models.Const1D()
'''
        with open('bubu__model__2__.py', 'w') as fobj:
            print(modelstring, file=fobj)
        b, kind, constraints = load_model('bubu__model__2__.py')
        assert isinstance(b, Model)
        assert kind == 'Astropy'
        assert b.amplitude == 1

    def test_load_model_input_not_string(self):
        '''Input is not a string'''
        with pytest.raises(TypeError) as record:
            b, kind, _ = load_model(1)
        assert 'modelstring has to be an existing file name' \
               in str(record.value)

    def test_load_model_input_file_doesnt_exist(self):
        with pytest.raises(FileNotFoundError) as record:
            b, kind, _ = load_model('dfasjkdaslfj')
        assert 'Model file not found' in str(record.value)

    def test_load_model_input_invalid_file_format(self):
        with open('bubu.txt', 'w') as fobj:
            print(1, file=fobj)
        with pytest.raises(TypeError) as record:
            b, kind, _ = load_model('bubu.txt')
        assert 'Unknown file type' in str(record.value)

    def test_load_data_fails(self):
        with pytest.raises(TypeError) as record:
            load_data('afile.fits')
        assert "The file type is not recognized" in str(record.value)

    @classmethod
    def teardown_class(cls):
        import shutil
        for dum in glob.glob('bubu*.*'):
            os.unlink(dum)

max_stars_repo_name: astrojuan/HENDRICS

max_issues_repo_name: astrojuan/HENDRICS

max_forks_repo_name: astrojuan/HENDRICS

reason: import numpy,from astropy


text:
import collections
import operator
import pandas as pd
import numpy as np

from .error_check import ErrorChecker
from .requirements import ResourceRequirements, CountableResourceRequirement
from .metric.metric_categories.size import Size

class RequestProcessingInfo:

    """
    Wraps static information on how request is processed in the application,
    e.g. in which order it passes through services, which amount of time spends
    there being processed, and which amount of resources it consumes at each
    simulation step.
    """

    @classmethod
    def build_from_config(cls, request_info : dict):

        request_type = ErrorChecker.key_check_and_load('request_type', request_info, 'request')
        entry_service = ErrorChecker.key_check_and_load('entry_service', request_info, 'request_type', request_type)

        processing_times = {}
        processing_times_raw = ErrorChecker.key_check_and_load('processing_times', request_info, 'request_type', request_type)
        processing_times_unit = ErrorChecker.key_check_and_load('unit', processing_times_raw, 'request_type', request_type)
        processing_times_vals = ErrorChecker.key_check_and_load('values', processing_times_raw, 'request_type', request_type)
        for processing_time_service_entry in processing_times_vals:

            service_name = ErrorChecker.key_check_and_load('service', processing_time_service_entry, 'request_type', request_type)
            upstream_time = ErrorChecker.key_check_and_load('upstream', processing_time_service_entry, 'request_type', request_type)
            downstream_time = ErrorChecker.key_check_and_load('downstream', processing_time_service_entry, 'request_type', request_type)

            processing_times[service_name] = { 'upstream' : cls._update_distribution_parameters(upstream_time, processing_times_unit), \
                                               'downstream' : cls._update_distribution_parameters(downstream_time, processing_times_unit) }

        timeout_raw = ErrorChecker.key_check_and_load('timeout', request_info, 'request_type', request_type)
        timeout_val = ErrorChecker.key_check_and_load('value', timeout_raw, 'request_type', request_type)
        timeout_unit = ErrorChecker.key_check_and_load('unit', timeout_raw, 'request_type', request_type)
        timeout = pd.Timedelta(timeout_val, unit = timeout_unit)
        ErrorChecker.value_check('timeout', timeout, operator.ge, pd.Timedelta(0, unit = timeout_unit), [f'request_type {request_type}'])

        request_size = CountableResourceRequirement.memory_from_dict(ErrorChecker.key_check_and_load('request_size', request_info, 'request_type', request_type))
        response_size = CountableResourceRequirement.memory_from_dict(ErrorChecker.key_check_and_load('response_size', request_info, 'request_type', request_type))

        request_operation_type = ErrorChecker.key_check_and_load('operation_type', request_info, 'request_type', request_type)

        request_processing_requirements = ErrorChecker.key_check_and_load('processing_requirements', request_info, 'request_type', request_type)

        return cls(request_type, entry_service, processing_times, timeout, request_size, response_size, request_operation_type, request_processing_requirements)

    @classmethod
    def _update_distribution_parameters(cls, config : dict, processing_times_unit : str):

        distribution_parameters = {'mean': 0, 'std': 0, 'unit': processing_times_unit}

        if isinstance(config, collections.Iterable):
            distribution_parameters['mean'] = ErrorChecker.key_check_and_load('mean', config, default = None)
            if distribution_parameters['mean'] is None:
                distribution_parameters['mean'] = ErrorChecker.key_check_and_load('value', config)
            distribution_parameters['unit'] = ErrorChecker.key_check_and_load('unit', config, default = processing_times_unit)
            if 'std' in config:
                distribution_parameters['std'] = ErrorChecker.key_check_and_load('std', config)

        else:
            distribution_parameters['mean'] = config

        return distribution_parameters

    def __init__(self,
                 request_type : str,
                 entry_service : str,
                 processing_times : dict,
                 timeout : pd.Timedelta,
                 request_size : CountableResourceRequirement,
                 response_size : CountableResourceRequirement,
                 request_operation_type : str,
                 request_processing_requirements : dict):

        self.request_type = request_type
        self.entry_service = entry_service
        self.processing_times = processing_times
        self.timeout = timeout
        self.request_size = request_size
        self.response_size = response_size
        self.request_operation_type = request_operation_type
        self.resource_requirements = ResourceRequirements.from_dict(request_processing_requirements)

    def get_upstream_processing_time(self, service_name : str) -> pd.Timedelta:

        return self._sample_processing_time(service_name, 'upstream')

    def get_downstream_processing_time(self, service_name : str) -> pd.Timedelta:

        return self._sample_processing_time(service_name, 'downstream')

    def _sample_processing_time(self, service_name : str, direction : str) -> pd.Timedelta:

        processing_time_value = max(np.random.normal(self.processing_times[service_name][direction]['mean'], self.processing_times[service_name][direction]['std']), 0)
        return pd.Timedelta(processing_time_value, unit = self.processing_times[service_name][direction]['unit'])

ext: py

max_stars_repo_name: Remit/autoscaling-simulator

max_issues_repo_name: Remit/autoscaling-simulator

max_forks_repo_name: Remit/autoscaling-simulator

reason: import numpy


text:
# -*- coding: utf-8 -*-
"""
Script for running nulls on simulations IN PARALLEL (n.b., there is NO
parallelization at the level of each null method)
"""

from argparse import ArgumentParser
from pathlib import Path
import time

from joblib import Parallel, delayed
import numpy as np

import threadpoolctl

from brainsmash import mapgen
from brainspace.null_models import moran
from netneurotools import (datasets as nndata,
                           freesurfer as nnsurf,
                           stats as nnstats)
from parspin import burt, simnulls, utils as putils

ROIDIR = Path('./data/raw/rois').resolve()
SPDIR = Path('./data/derivatives/spins').resolve()
DISTDIR = Path('./data/derivatives/geodesic').resolve()
SIMDIR = Path('./data/derivatives/simulated').resolve()

ALPHA = 0.05  # p-value threshold
N_PROC = 24  # number of parallel workers for surrogate generation
N_PERM = 1000  # number of permutations for null models
N_SIM = 1000  # number of simulations to run
SEED = 1234  # reproducibility
SHUFFLE = True  # if we're shuffling sims instead of running paired (r = 0.15)
USE_KNN = False  # whether to use default nearest neigh setting for Burt-2020

def make_surrogates(data, parcellation, scale, spatnull):
    """
    Generates surrogates for `data` using `spatnull` method

    Parameters
    ----------
    data : (N,) pd.DataFrame
    parcellation : {'atl-cammoun2012', 'atl-schaefer2018'}
    scale : str
    spatnull : {'burt2018', 'burt2020', 'moran'}

    Returns
    -------
    surrogates : (N, `N_PERM`) np.ndarray
    """

    if spatnull not in ('burt2018', 'burt2020', 'moran'):
        raise ValueError(f'Cannot make surrogates for null method {spatnull}')

    darr = np.asarray(data)
    dmin = darr[np.logical_not(np.isnan(darr))].min()

    surrogates = np.zeros((len(data), N_PERM))
    for hdata, dist, idx in putils.yield_data_dist(
        DISTDIR, parcellation, scale, data, inverse=(spatnull == 'moran')
    ):

        # handle NaNs before generating surrogates; should only be relevant
        # when using vertex-level data, but good nonetheless
        mask = np.logical_not(np.isnan(hdata))
        surrogates[idx[np.logical_not(mask)]] = np.nan
        hdata, dist, idx = hdata[mask], dist[np.ix_(mask, mask)], idx[mask]

        if spatnull == 'burt2018':
            # Box-Cox transformation requires positive data :man_facepalming:
            hdata += np.abs(dmin) + 0.1
            surrogates[idx] = \
                burt.batch_surrogates(dist, hdata, n_surr=N_PERM, seed=SEED)
        elif spatnull == 'burt2020':
            if parcellation == 'vertex':  # memmap is required for this shit
                index = np.argsort(dist, axis=-1)
                dist = np.sort(dist, axis=-1)
                knn = 1000 if USE_KNN else len(hdata)
                surrogates[idx] = \
                    mapgen.Sampled(hdata, dist, index, knn=knn,
                                   seed=SEED)(N_PERM).T
            else:
                surrogates[idx] = \
                    mapgen.Base(hdata, dist, seed=SEED)(N_PERM, 50).T
        elif spatnull == 'moran':
            mrs = moran.MoranRandomization(joint=True, n_rep=N_PERM,
                                           tol=1e-6, random_state=SEED)
            with threadpoolctl.threadpool_limits(limits=2):
                surrogates[idx] = mrs.fit(dist).randomize(hdata).T

    return surrogates

def _cornblath(x, y, spins, annot):
    """ Calculates p-value using Cornblath method for provided data
    """
    nulls = nnsurf.spin_data(y, version='fsaverage5',
                             lhannot=annot.lh, rhannot=annot.rh,
                             spins=spins, n_rotate=spins.shape[-1])

    return simnulls.calc_pval(x, y, nulls)

def _baum(x, y, spins):
    """ Calculates p-value using Cornblath method for provided data
    """
    nulls = y[spins]
    nulls[spins == -1] = np.nan

    return simnulls.calc_pval(x, y, nulls)

def _genmod(x, y, parcellation, scale, spatnull):
    """ Calculates p-value w/Burt2018/Burt2020/MRS method for provided data
    """
    nulls = make_surrogates(y, parcellation, scale, spatnull)
    return simnulls.calc_pval(x, y, nulls)

def _get_ysim(y, sim):
    """ Gets `sim` column from `y`, accounting for DataFrame vs ndarray
    """
    try:
        return y.iloc[:, sim]
    except AttributeError:
        return y[:, sim]

def run_null(parcellation, scale, spatnull, alpha):
    """
    Runs spatial null models for given combination of inputs

    Parameters
    ----------
    parcellation : str
        Name of parcellation to be used
    scale : str
        Scale of `parcellation` to be used
    spatnull : str
        Name of spin method to be used
    alpha : float
        Spatial autocorrelation parameter to be used
    """

    print(f'{time.ctime()}: {parcellation} {scale} {spatnull} {alpha} ',
          flush=True)

    # filenames (for I/O)
    spins_fn = SPDIR / parcellation / spatnull / f'{scale}_spins.csv'
    pvals_fn = (SIMDIR / alpha / parcellation / 'nulls' / spatnull
                / f'{scale}_nulls.csv')
    perms_fn = pvals_fn.parent / f'{scale}_perms.csv'

    if SHUFFLE:
        pvals_fn = pvals_fn.parent / f'{scale}_nulls_shuffle.csv'
        perms_fn = perms_fn.parent / f'{scale}_perms_shuffle.csv'

    if pvals_fn.exists() and perms_fn.exists():
        return

    # load simulated data
    alphadir = SIMDIR / alpha
    if parcellation == 'vertex':
        x, y = simnulls.load_vertex_data(alphadir, n_sim=N_SIM)
    else:
        x, y = simnulls.load_parc_data(alphadir, parcellation, scale,
                                       n_sim=N_SIM)

    # if we're computing info on SHUFFLED data, get the appropriate random `y`
    if SHUFFLE:
        y = _get_ysim(y, np.random.default_rng(1).permutation(N_SIM))

    # calculate the null p-values
    if spatnull == 'naive-para':
        pvals = nnstats.efficient_pearsonr(x, y, nan_policy='omit')[1]
        perms = np.array([np.nan])
    elif spatnull == 'cornblath':
        fn = SPDIR / 'vertex' / 'vazquez-rodriguez' / 'fsaverage5_spins.csv'
        x, y = np.asarray(x), np.asarray(y)
        spins = simnulls.load_spins(fn, n_perm=N_PERM)
        fetcher = getattr(nndata, f"fetch_{parcellation.replace('atl-', '')}")
        annot = fetcher('fsaverage5', data_dir=ROIDIR)[scale]
        out = Parallel(n_jobs=N_PROC, max_nbytes=None)(
            delayed(_cornblath)(x[:, sim], y[:, sim], spins, annot)
            for sim in putils.trange(x.shape[-1], desc='Running simulations')
        )
        pvals, perms = zip(*out)
    elif spatnull == 'baum':
        x, y = np.asarray(x), np.asarray(y)
        spins = simnulls.load_spins(spins_fn, n_perm=N_PERM)
        out = Parallel(n_jobs=N_PROC, max_nbytes=None)(
            delayed(_baum)(x[:, sim], y[:, sim], spins)
            for sim in putils.trange(x.shape[-1], desc='Running simulations')
        )
        pvals, perms = zip(*out)
    elif spatnull in ('burt2018', 'burt2020', 'moran'):
        xarr = np.asarray(x)
        out = Parallel(n_jobs=N_PROC, max_nbytes=None)(
            delayed(_genmod)(xarr[:, sim], _get_ysim(y, sim),
                             parcellation, scale, spatnull)
            for sim in putils.trange(x.shape[-1], desc='Running simulations')
        )
        pvals, perms = zip(*out)
    else:  # vazquez-rodriguez, vasa, hungarian, naive-nonparametric
        x, y = np.asarray(x), np.asarray(y)
        spins = simnulls.load_spins(spins_fn, n_perm=N_PERM)
        out = Parallel(n_jobs=N_PROC, max_nbytes=None)(
            delayed(simnulls.calc_pval)(x[:, sim], y[:, sim], y[spins, sim])
            for sim in putils.trange(x.shape[-1], desc='Running simulations')
        )
        pvals, perms = zip(*out)

    # save to disk
    putils.save_dir(perms_fn, np.atleast_1d(perms), overwrite=False)
    putils.save_dir(pvals_fn, np.atleast_1d(pvals), overwrite=False)

def main():
    parcellations = putils.get_cammoun_schaefer(data_dir=ROIDIR)

    # get inputs
    args = get_parser()

    # reset some stuff
    for param in ('n_perm', 'n_proc', 'seed', 'n_sim', 'shuffle', 'use_knn'):
        globals()[param.upper()] = args[param]

    print(f'N_PERM: {N_PERM}',
          f'N_PROC: {N_PROC}',
          f'SEED: {SEED}',
          f'SHUFFLE: {SHUFFLE}',
          f'N_SIM: {N_SIM}',
          f'SPATNULLS: {args["spatnull"]}',
          f'ALPHAS: {args["alpha"]}\n', sep='\n')

    if args['show_params']:
        return

    # everyone loves a four-level-deep nested for-loop :man_facepalming:
    for spatnull in args['spatnull']:
        for alpha in args['alpha']:
            if spatnull in simnulls.VERTEXWISE:
                run_null('vertex', 'fsaverage5', spatnull, alpha)
            for parcellation, annotations in parcellations.items():
                for scale in annotations:
                    run_null(parcellation, scale, spatnull, alpha)

def get_parser():
    parser = ArgumentParser()
    parser.add_argument('--show_params', default=False, action='store_true')
    parser.add_argument('--n_perm', default=N_PERM, type=int)
    parser.add_argument('--n_proc', default=N_PROC, type=int)
    parser.add_argument('--seed', default=SEED, type=int)
    parser.add_argument('--shuffle', default=False, action='store_true')
    parser.add_argument('--use_knn', default=False, action='store_true')
    parser.add_argument('--spatnull', choices=simnulls.SPATNULLS,
                        default=simnulls.SPATNULLS, nargs='+')
    parser.add_argument('--alpha', choices=simnulls.ALPHAS,
                        default=simnulls.ALPHAS, nargs='+')
    parser.add_argument('n_sim', type=int, default=N_SIM, nargs='?')
    return vars(parser.parse_args())

if __name__ == "__main__":
    main()


text:
"""Plot scalar and vector fields."""

from itertools import chain

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.patches as mpat
from mpl_toolkits.axes_grid1 import make_axes_locatable

from . import conf, phyvars, _helpers
from .error import NotAvailableError
from .stagyydata import StagyyData

# The location is off for vertical velocities: they have an extra
# point in (x,y) instead of z in the output

def _threed_extract(step, var, walls=False):
    """Return suitable slices and coords for 3D fields."""
    is_vector = not valid_field_var(var)
    hwalls = is_vector or walls
    i_x = conf.field.ix
    i_y = conf.field.iy
    i_z = conf.field.iz
    if i_x is not None or i_y is not None:
        i_z = None
    if i_x is not None or i_z is not None:
        i_y = None
    if i_x is None and i_y is None and i_z is None:
        i_x = 0
    if i_x is not None:
        xcoord = step.geom.y_walls if hwalls else step.geom.y_centers
        ycoord = step.geom.z_walls if walls else step.geom.z_centers
        i_y = i_z = slice(None)
        varx, vary = var + '2', var + '3'
    elif i_y is not None:
        xcoord = step.geom.x_walls if hwalls else step.geom.x_centers
        ycoord = step.geom.z_walls if walls else step.geom.z_centers
        i_x = i_z = slice(None)
        varx, vary = var + '1', var + '3'
    else:
        xcoord = step.geom.x_walls if hwalls else step.geom.x_centers
        ycoord = step.geom.y_walls if hwalls else step.geom.y_centers
        i_x = i_y = slice(None)
        varx, vary = var + '1', var + '2'
    if is_vector:
        data = (step.fields[varx].values[i_x, i_y, i_z, 0],
                step.fields[vary].values[i_x, i_y, i_z, 0])
    else:
        data = step.fields[var].values[i_x, i_y, i_z, 0]
    return (xcoord, ycoord), data

def valid_field_var(var):
    """Whether a field variable is defined.

    This function checks if a definition of the variable exists in
    :data:`~stagpy.phyvars.FIELD` or :data:`~stagpy.phyvars.FIELD_EXTRA`.

    Args:
        var (str): the variable name to be checked.
    Returns:
        bool: True is the var is defined, False otherwise.
    """
    return var in phyvars.FIELD or var in phyvars.FIELD_EXTRA

def get_meshes_fld(step, var, walls=False):
    """Return scalar field along with coordinates meshes.

    Only works properly in 2D geometry and 3D cartesian.

    Args:
        step (:class:`~stagpy._step.Step`): a step of a StagyyData instance.
        var (str): scalar field name.
        walls (bool): consider the walls as the relevant mesh.
    Returns:
        tuple of :class:`numpy.array`: xmesh, ymesh, fld, meta
            2D arrays containing respectively the x position, y position, the
            values and the metadata of the requested field.
    """
    fld, meta = step.fields[var]
    hwalls = (walls or fld.shape[0] != step.geom.nxtot or
              fld.shape[1] != step.geom.nytot)
    if step.geom.threed and step.geom.cartesian:
        (xcoord, ycoord), fld = _threed_extract(step, var, walls)
    elif step.geom.twod_xz:
        xcoord = step.geom.x_walls if hwalls else step.geom.x_centers
        ycoord = step.geom.z_walls if walls else step.geom.z_centers
        fld = fld[:, 0, :, 0]
    else:  # twod_yz
        xcoord = step.geom.y_walls if hwalls else step.geom.y_centers
        ycoord = step.geom.z_walls if walls else step.geom.z_centers
        if step.geom.curvilinear:
            pmesh, rmesh = np.meshgrid(xcoord, ycoord, indexing='ij')
            xmesh, ymesh = rmesh * np.cos(pmesh), rmesh * np.sin(pmesh)
        fld = fld[0, :, :, 0]
    if step.geom.cartesian:
        xmesh, ymesh = np.meshgrid(xcoord, ycoord, indexing='ij')
    return xmesh, ymesh, fld, meta

def get_meshes_vec(step, var):
    """Return vector field components along with coordinates meshes.

    Only works properly in 2D geometry and 3D cartesian.

    Args:
        step (:class:`~stagpy._step.Step`): a step of a StagyyData instance.
        var (str): vector field name.
    Returns:
        tuple of :class:`numpy.array`: xmesh, ymesh, fldx, fldy
            2D arrays containing respectively the x position, y position, x
            component and y component of the requested vector field.
    """
    if step.geom.threed and step.geom.cartesian:
        (xcoord, ycoord), (vec1, vec2) = _threed_extract(step, var)
    elif step.geom.twod_xz:
        xcoord, ycoord = step.geom.x_walls, step.geom.z_centers
        vec1 = step.fields[var + '1'].values[:, 0, :, 0]
        vec2 = step.fields[var + '3'].values[:, 0, :, 0]
    elif step.geom.cartesian and step.geom.twod_yz:
        xcoord, ycoord = step.geom.y_walls, step.geom.z_centers
        vec1 = step.fields[var + '2'].values[0, :, :, 0]
        vec2 = step.fields[var + '3'].values[0, :, :, 0]
    else:  # spherical yz
        pcoord = step.geom.p_walls
        pmesh = np.outer(pcoord, np.ones(step.geom.nrtot))
        vec_phi = step.fields[var + '2'].values[0, :, :, 0]
        vec_r = step.fields[var + '3'].values[0, :, :, 0]
        vec1 = vec_r * np.cos(pmesh) - vec_phi * np.sin(pmesh)
        vec2 = vec_phi * np.cos(pmesh) + vec_r * np.sin(pmesh)
        pcoord, rcoord = step.geom.p_walls, step.geom.r_centers
        pmesh, rmesh = np.meshgrid(pcoord, rcoord, indexing='ij')
        xmesh, ymesh = rmesh * np.cos(pmesh), rmesh * np.sin(pmesh)
    if step.geom.cartesian:
        xmesh, ymesh = np.meshgrid(xcoord, ycoord, indexing='ij')
    return xmesh, ymesh, vec1, vec2

def plot_scalar(step, var, field=None, axis=None, **extra):
    """Plot scalar field.

    Args:
        step (:class:`~stagpy._step.Step`): a step of a StagyyData instance.
        var (str): the scalar field name.
        field (:class:`numpy.array`): if not None, it is plotted instead of
            step.fields[var].  This is useful to plot a masked or rescaled
            array.  Note that if conf.scaling.dimensional is True, this
            field will be scaled accordingly.
        axis (:class:`matplotlib.axes.Axes`): the axis objet where the field
            should be plotted.  If set to None, a new figure with one subplot
            is created.
        extra (dict): options that will be passed on to
            :func:`matplotlib.axes.Axes.pcolormesh`.
    Returns:
        fig, axis, surf, cbar
            handles to various :mod:`matplotlib` objects, respectively the
            figure, the axis, the surface returned by
            :func:`~matplotlib.axes.Axes.pcolormesh`, and the colorbar returned
            by :func:`matplotlib.pyplot.colorbar`.
    """
    if step.geom.threed and step.geom.spherical:
        raise NotAvailableError(
            'plot_scalar not implemented for 3D spherical geometry')

    xmesh, ymesh, fld, meta = get_meshes_fld(step, var,
                                             walls=not conf.field.interpolate)
    if conf.field.interpolate and \
       step.geom.spherical and step.geom.twod_yz and \
       fld.shape[0] == step.geom.nytot:
        # add one point to close spherical annulus
        xmesh = np.concatenate((xmesh, xmesh[:1]), axis=0)
        ymesh = np.concatenate((ymesh, ymesh[:1]), axis=0)
        newline = (fld[:1] + fld[-1:]) / 2
        fld = np.concatenate((fld, newline), axis=0)
    xmin, xmax = xmesh.min(), xmesh.max()
    ymin, ymax = ymesh.min(), ymesh.max()

    if field is not None:
        fld = field
    if conf.field.perturbation:
        fld = fld - np.mean(fld, axis=0)
    if conf.field.shift:
        fld = np.roll(fld, conf.field.shift, axis=0)

    fld, unit = step.sdat.scale(fld, meta.dim)

    if axis is None:
        fig, axis = plt.subplots(ncols=1)
    else:
        fig = axis.get_figure()

    if step.sdat.par['magma_oceans_in']['magma_oceans_mode']:
        rcmb = step.sdat.par['geometry']['r_cmb']
        xmax = rcmb + 1
        ymax = xmax
        xmin = -xmax
        ymin = -ymax
        rsurf = xmax if step.timeinfo['thick_tmo'] > 0 \
            else step.geom.r_walls[0, 0, -3]
        cmb = mpat.Circle((0, 0), rcmb, color='dimgray', zorder=0)
        psurf = mpat.Circle((0, 0), rsurf, color='indianred', zorder=0)
        axis.add_patch(psurf)
        axis.add_patch(cmb)

    extra_opts = dict(
        cmap=conf.field.cmap.get(var),
        vmin=conf.plot.vmin,
        vmax=conf.plot.vmax,
        norm=mpl.colors.LogNorm() if var == 'eta' else None,
        rasterized=conf.plot.raster,
        shading='gouraud' if conf.field.interpolate else 'flat',
    )
    extra_opts.update(extra)
    surf = axis.pcolormesh(xmesh, ymesh, fld, **extra_opts)

    cbar = None
    if conf.field.colorbar:
        cax = make_axes_locatable(axis).append_axes(
            'right', size="3%", pad=0.15)
        cbar = plt.colorbar(surf, cax=cax)
        cbar.set_label(meta.description +
                       (' pert.' if conf.field.perturbation else '') +
                       (f' ({unit})' if unit else ''))
    if step.geom.spherical or conf.plot.ratio is None:
        axis.set_aspect('equal')
        axis.set_axis_off()
    else:
        axis.set_aspect(conf.plot.ratio / axis.get_data_ratio())

    axis.set_adjustable('box')
    axis.set_xlim(xmin, xmax)
    axis.set_ylim(ymin, ymax)
    return fig, axis, surf, cbar

def plot_iso(axis, step, var, **extra):
    """Plot isocontours of scalar field.

    Args:
        axis (:class:`matplotlib.axes.Axes`): the axis handler of an
            existing matplotlib figure where the isocontours should
            be plotted.
        step (:class:`~stagpy._step.Step`): a step of a StagyyData instance.
        var (str): the scalar field name.
        extra (dict): options that will be passed on to
            :func:`matplotlib.axes.Axes.contour`.
    """
    xmesh, ymesh, fld, _ = get_meshes_fld(step, var)
    if conf.field.shift:
        fld = np.roll(fld, conf.field.shift, axis=0)
    extra_opts = dict(linewidths=1)
    if 'cmap' not in extra and conf.field.isocolors:
        extra_opts['colors'] = conf.field.isocolors.split(',')
    elif 'colors' not in extra:
        extra_opts['cmap'] = conf.field.cmap.get(var)
    if conf.plot.isolines:
        extra_opts['levels'] = sorted(conf.plot.isolines)
    extra_opts.update(extra)
    axis.contour(xmesh, ymesh, fld, **extra_opts)

def plot_vec(axis, step, var):
    """Plot vector field.

    Args:
        axis (:class:`matplotlib.axes.Axes`): the axis handler of an
            existing matplotlib figure where the vector field should
            be plotted.
        step (:class:`~stagpy._step.Step`): a step of a StagyyData instance.
        var (str): the vector field name.
    """
    xmesh, ymesh, vec1, vec2 = get_meshes_vec(step, var)
    dipz = step.geom.nztot // 10
    if conf.field.shift:
        vec1 = np.roll(vec1, conf.field.shift, axis=0)
        vec2 = np.roll(vec2, conf.field.shift, axis=0)
    if step.geom.spherical or conf.plot.ratio is None:
        dipx = dipz
    else:
        dipx = step.geom.nytot if step.geom.twod_yz else step.geom.nxtot
        dipx = int(dipx // 10 * conf.plot.ratio) + 1
    axis.quiver(xmesh[::dipx, ::dipz], ymesh[::dipx, ::dipz],
                vec1[::dipx, ::dipz], vec2[::dipx, ::dipz],
                linewidths=1)

def _findminmax(sdat, sovs):
    """Find min and max values of several fields."""
    minmax = {}
    for step in sdat.walk.filter(snap=True):
        for var in sovs:
            if var in step.fields:
                field, meta = step.fields[var]
                field, _ = sdat.scale(field, meta.dim)
                if var in minmax:
                    minmax[var] = (min(minmax[var][0], np.nanmin(field)),
                                   max(minmax[var][1], np.nanmax(field)))
                else:
                    minmax[var] = np.nanmin(field), np.nanmax(field)
    return minmax

def cmd():
    """Implementation of field subcommand.

    Other Parameters:
        conf.field
        conf.core
    """
    sdat = StagyyData()
    lovs = _helpers.list_of_vars(conf.field.plot)
    # no more than two fields in a subplot
    lovs = [[slov[:2] for slov in plov] for plov in lovs]
    minmax = {}
    if conf.plot.cminmax:
        conf.plot.vmin = None
        conf.plot.vmax = None
        sovs = set(slov[0] for plov in lovs for slov in plov)
        minmax = _findminmax(sdat, sovs)
    for step in sdat.walk.filter(snap=True):
        for vfig in lovs:
            fig, axes = plt.subplots(ncols=len(vfig), squeeze=False,
                                     figsize=(6 * len(vfig), 6))
            for axis, var in zip(axes[0], vfig):
                if var[0] not in step.fields:
                    print(f"{var[0]!r} field on snap {step.isnap} not found")
                    continue
                opts = {}
                if var[0] in minmax:
                    opts = dict(vmin=minmax[var[0]][0], vmax=minmax[var[0]][1])
                plot_scalar(step, var[0], axis=axis, **opts)
                if len(var) == 2:
                    if valid_field_var(var[1]):
                        plot_iso(axis, step, var[1])
                    elif valid_field_var(var[1] + '1'):
                        plot_vec(axis, step, var[1])
            if conf.field.timelabel:
                time, unit = sdat.scale(step.timeinfo['t'], 's')
                time = _helpers.scilabel(time)
                axes[0, 0].text(0.02, 1.02, f'$t={time}$ {unit}',
                                transform=axes[0, 0].transAxes)
            oname = '_'.join(chain.from_iterable(vfig))
            plt.tight_layout(w_pad=3)
            _helpers.saveplot(fig, oname, step.isnap)

max_stars_repo_name: StagPython/StagPy

max_issues_repo_name: StagPython/StagPy

max_forks_repo_name: StagPython/StagPy


text:
# Copyright (c) 2018 Eric Kerfoot, see LICENSE file
'''
This script will read in the original cardiac-dig.mat file and produce a .npz file containing multiple clones of the
original data with the spatial deformation applied. Since the augment is way too slow to use in online during training
this is needed to statically generate variations of the original data.
'''

import multiprocessing as mp
from multiprocessing import sharedctypes
import numpy as np

from scipy.interpolate import interp2d
from scipy.ndimage import geometric_transform
from scipy.io import loadmat

def mapping(coords,interx,intery):
    '''Map the coordinates from the original image to the deformed one by adding values from the given interp2d objects.'''
    y,x=coords[:2]
    dx=interx(x,y)
    dy=intery(x,y)

    return (y+dy,x+dx)+tuple(coords[2:])

def deformBothAugment(img,mask,defrange=5):
    '''
    Deform the central parts of the image/mask pair using interpolated randomized deformation. This is a particularly
    slow implementation of this concept, however it produces good results and other versions are still too slow to use
    online during training anyway.
    '''
    grid=np.zeros((4,4,2),int)
    y=np.linspace(0,img.shape[0],grid.shape[0])
    x=np.linspace(0,img.shape[1],grid.shape[1])

    grid[1:3,1:3,:]=2*defrange*np.random.random_sample((2,2,2))-defrange

    interx=interp2d(x,y,grid[...,0],'linear')
    intery=interp2d(x,y,grid[...,1],'linear')
    xargs=(interx,intery)

    return geometric_transform(img,mapping,extra_arguments=xargs),geometric_transform(mask,mapping,extra_arguments=xargs)

def toShared(arr):
    '''Convert the given Numpy array to a shared ctypes object.'''
    carr=np.ctypeslib.as_ctypes(arr)
    return sharedctypes.RawArray(carr._type_, carr)

def fromShared(arr):
    '''Map the given ctypes object to a Numpy array, this is expected to be a shared object from the parent.'''
    return np.ctypeslib.as_array(arr)

data=loadmat('./cardiac-dig.mat')

images_LV=data['images_LV'].astype(np.float32)
endo_LV=data['endo_LV']
epi_LV=data['epi_LV']

numClones=20

output_LV=np.zeros(tuple(images_LV.shape)+(numClones+1,),dtype=images_LV.dtype)
endoout_LV=np.zeros(tuple(endo_LV.shape)+(numClones+1,),dtype=endo_LV.dtype)
epiout_LV=np.zeros(tuple(epi_LV.shape)+(numClones+1,),dtype=epi_LV.dtype)

output_LV[...,0]=images_LV
endoout_LV[...,0]=endo_LV
epiout_LV[...,0]=epi_LV

output_LV=toShared(output_LV)
endoout_LV=toShared(endoout_LV)
epiout_LV=toShared(epiout_LV)

def generate(start,end):
    '''
    This generator is executed in parallel on subprocesses, filling data into the shared arrays created above. This
    relies on fork() semantics so will not work in Windows.
    '''
    loutput_LV=fromShared(output_LV)
    lendoout_LV=fromShared(endoout_LV)
    lepiout_LV=fromShared(epiout_LV)

    for i in range(start,end):
        print(i,start,end)

        im0=images_LV[...,i]
        endo=endo_LV[...,i]
        epi=epi_LV[...,i]
        mask0=np.dstack([endo,epi])

        for c in range(numClones):
            im,mask=deformBothAugment(im0,mask0,5)
            loutput_LV[...,i,c+1]=im
            lendoout_LV[...,i,c+1]=mask[...,0]
            lepiout_LV[...,i,c+1]=mask[...,1]

if __name__=='__main__':
    ranges=[(a[0],a[-1]+1) for a in np.array_split(np.arange(images_LV.shape[2]),mp.cpu_count())] # list of start,end indices

    with mp.Pool(processes=mp.cpu_count()) as pool:
        pool.starmap(generate,ranges)

    np.savez_compressed('./cardiac-dig_cloned_%i.npz'%numClones,numClones=numClones,
                        images_LV=fromShared(output_LV),
                        endo_LV=fromShared(endoout_LV),
                        epi_LV=fromShared(epiout_LV))

reason: import numpy,from scipy


text:
from sklearn.feature_selection import SelectorMixin
from sklearn.base import BaseEstimator
from sklearn.utils import check_array
from sklearn.utils.validation import check_is_fitted

from scipy.stats import pearsonr
import pandas as pd
import numpy as np

class FilterSelector(SelectorMixin, BaseEstimator):
    def __init__(self, r = 0.9, alpha = 0.05):
        '''
        r: 相関係数を高いとする閾値．これより高いと相関が高いとして，削除対象とする．
        alpha: 有意水準．これよりp値が低い時，相関係数の値が信頼に足ると判断．特別な理由がない限りはこれを変えないのがベター．
        '''
        self.r = r
        self.alpha = alpha

    def fit(self, X, y = None):
        '''
        X: pd.DataFrame or 2次元のnp.ndarray or 2次元のlist．すでに作成された特徴量空間．
        y: 不要．
        '''
        # checkかつnp.ndarrayに変換される．
        X = check_array(X)

        # よくつかうので列数を変数で保存
        n_features = X.shape[1]

        # 相関係数とそのp値を入れるリストの鋳型を作成
        # [i, j, 0]が相関係数， [i, j, 1]がp値
        self.corr_ = np.empty((n_features, n_features, 2), dtype = float)

        # i番目の要素は，i番目と相関係数がcorr以上かつp値が有意水準alpha未満のものの集合．これの鋳型を作成
        pair = [set() for _ in range(n_features)]

        for i in range(n_features):
            for j in range(i, n_features):
                if i == j:  # 一緒の時はそもそも相関係数1, p値0は自明なので例外処理
                    self.corr_[i, i] = [1.0, 0.0]
                else:
                    self.corr_[i, j] = list(pearsonr(X[:, i], X[:, j]))
                    self.corr_[j, i] = self.corr_[i, j]
                    if self.corr_[i, j, 0] > self.r and self.corr_[i, j, 1] < self.alpha:    # 相関係数が高く， かつp値が小さい (信頼に足る) とき．
                        # i, j番目の要素の集合にそれぞれj, iを追加
                        pair[i].add(j)
                        pair[j].add(i)

        # 何個相関係数が高い係数が存在するのかの値を求める．
        def _delete_recursive(pair, boolean = np.ones(n_features, dtype = bool)):
            # 相関係数が高いものがいくつあるのかのnp.ndarrayを作成
            order_pair = np.array([len(s) for s in pair])

            # これがすべて0のとき，もう削除は終わってるので終了処理
            if np.sum(order_pair) == 0:
                return boolean

            # 一番相関係数が高いペア数を多く持っているものの要素番号 (index) を取得
            i_max = np.argmax(order_pair)

            # 相関係数高いペアからその番号を削除
            pair = [p.difference({i_max}) for p in pair]

            # その集合は削除かつ，その特徴量はbooleanをFalseに．
            pair[i_max] = set()
            boolean[i_max] = False
            return _delete_recursive(pair, boolean)

        self.support_ = _delete_recursive(pair)

        return self

    def _get_support_mask(self):
        check_is_fitted(self, 'support_')
        return self.support_

if __name__ == '__main__':
    from pdb import set_trace
    from sklearn.datasets import load_boston

    data = load_boston()
    X = pd.DataFrame(data.data, columns = data.feature_names)
    # X = pd.concat([X, pd.Series(np.ones(X.shape[0]))], axis = 1)

    selector = FilterSelector()
    selector.fit(X)

    set_trace()


text:
"""
nodes for scan operations
"""

from __future__ import division, absolute_import
from __future__ import print_function, unicode_literals

import theano

from .. import core
from .. import utils

@core.register_node("scan_input")
class ScanInputNode(core.NodeImpl):

    """
    Node that transfroms a sequence-wise input to an element-wise input,
    for replacement in a ScanNode.
    """

    def compute_output(self, network, in_vw):
        scan_axis = network.find_hyperparameter(["scan_axis"], 1)

        def remove_scan_axis(in_val):
            out_val = list(in_val)
            out_val.pop(scan_axis)
            # assuming in_val is a list or tuple
            if isinstance(in_val, tuple):
                return tuple(out_val)
            else:
                assert isinstance(in_val, list)
                return out_val

        # construct output
        network.create_vw(
            name="default",
            is_shared=False,
            shape=remove_scan_axis(in_vw.shape),
            broadcastable=remove_scan_axis(in_vw.broadcastable),
            tags={"input"},
        )

@core.register_node("scan_state")
class ScanStateNode(core.NodeImpl):

    """
    Container for hidden state, where one specifies an initial value and
    a next value via other nodes in the tree.

    NOTE: initial state MUST have the correct shape (ie. same shape as
    the output of the next_state node)

    initial_scan_state_reference:
    optional reference to a node to take initial state from
    if not given, the default input of the node is used
    """

    hyperparameter_names = ("initial_scan_state_reference",
                            "initial_state_reference",
                            "initial_state",
                            "next_scan_state_reference",
                            "next_state_reference",
                            "next_state")
    input_keys = ("initial_state",)

    def init_state(self, network):
        node_name = network.find_hyperparameter(
            ["initial_scan_state_reference",
             "initial_state_reference",
             "initial_state"],
            None)
        if node_name is None:
            # otherwise set it to the default input of the node
            node_name = network.get_all_input_edges()["default"]
        # add dependency in dag
        network.take_output_from(
            node_name,
            to_key="initial_state")
        # set next state node
        next_state = network.find_hyperparameter(["next_scan_state_reference",
                                                  "next_state_reference",
                                                  "next_state"])
        network.set_data("next_state", next_state)

    def compute_output(self, network, initial_state):
        # copy initial state, to be later referenced
        network.copy_vw(
            name="initial_state",
            previous_vw=initial_state,
            tags={"input"},
        )
        # create a new variable representing the output of this node,
        # so that the scan node can replace it with the node's input at a
        # previous time step
        network.create_vw(
            name="default",
            is_shared=False,
            shape=initial_state.shape,
            broadcastable=initial_state.broadcastable,
            tags={"output"},
        )

@core.register_node("scan")
class ScanNode(core.Wrapper1NodeImpl):

    """
    root node for a scan operation. transforms an element-wise child subtree
    into a sequence-wise subtree, taking into account ScanStateNode's
    """

    hyperparameter_names = ("scan_axis", )
    input_keys = ("default", "final_child_output",)

    def __init__(self, name, *args, **kwargs):
        super(ScanNode, self).__init__(name, *args, **kwargs)
        self._scan_input_node = ScanInputNode(name=self.name + "_input")

    def architecture_children(self):
        return ([self._scan_input_node]
                + super(ScanNode, self).architecture_children())

    def init_state(self, network):
        # send input to scan input node, and take output from child
        super(ScanNode, self).init_state(network)
        # link scan input node to child
        child1, child2 = self.architecture_children()
        network.add_dependency(child1.name, child2.name)

    def compute_output(self, network, sequence_input, element_output):
        scan_axis = network.find_hyperparameter(["scan_axis"], 1)
        # FIXME delete hard coded axis when actually shuffling dimensions
        scan_axis = 0

        # ################################ sequences ##########################

        # for now, only the single ScanInputNode
        input_network = network[self._scan_input_node.name]
        element_input = input_network.get_vw("default").variable
        element_input_vars = [element_input]
        # FIXME shuffle dimensions to use scan_axis
        # TODO assert ndim is right
        input_sequences = [sequence_input.variable]

        # ############################### outputs_info #######################
        # ---
        # sources:
        # 1. all outputs of the subtree
        # 2. next state for scan state node

        # outputs
        # original_outputs are ones whose variable will be replaced with a
        # version with an additional scan axis
        # NOTE: having  variables from nested scans is desired here, because
        # we want to transform the output into a version with a sequence
        # dimension
        original_outputs = network.find_vws_in_subtree(tags={"output"})
        # set all_outputs to a copy of original_outputs
        all_outputs = list(original_outputs)
        # in case the element-wise output of the scan is not tagged as an
        # "output" variable, we don't want to replace it with a sequence
        if element_output not in all_outputs:
            all_outputs.append(element_output)
        # find corresponding variable wrappers
        element_output_vars = [variable_wrapper.variable
                               for variable_wrapper in all_outputs]
        # set initial outputs_info that we do not care about the values
        outputs_info = [None] * len(all_outputs)

        # scan state nodes
        # FIXME do not get state in a nested scan
        scan_state_nodes = network.find_nodes_in_subtree(ScanStateNode)
        scan_state_networks = [network[node.name] for node in scan_state_nodes]
        scan_state_outputs = [net.get_vw("default")
                              for net in scan_state_networks]
        scan_state_vars = [variable_wrapper.variable
                           for variable_wrapper in scan_state_outputs]
        # finding idxs of scan state nodes in
        scan_state_idxs = [all_outputs.index(var)
                           for var in scan_state_outputs]
        # find the order that the states will appear in scan's inputs
        scan_state_order = sorted(range(len(scan_state_idxs)),
                                  key=lambda x: scan_state_idxs[x])
        # get the node name of the next state from each scan state node
        scan_state_next_names = [net.get_data("next_state")
                                 for net in scan_state_networks]
        # each scan state should probably have a unique next state
        # delete the assertion if this assumption does not hold
        # (the code should work just fine, this is just a sanity check)
        assert len(scan_state_next_names) == len(set(scan_state_next_names))
        scan_state_next_networks = [network[name]
                                    for name in scan_state_next_names]
        scan_state_next_vws = [net.get_vw("default")
                               for net in scan_state_next_networks]
        scan_state_next_idxs = [all_outputs.index(var)
                                for var in scan_state_next_vws]
        # finding initial states
        scan_state_initial_vws = [net.get_vw("initial_state")
                                  for net in scan_state_networks]
        # updates outputs_info to contain initial state
        for idx, node, init_vw, next_vw in zip(scan_state_idxs,
                                               scan_state_nodes,
                                               scan_state_initial_vws,
                                               scan_state_next_vws):
            # make sure initial and final shape are the same
            # ---
            # NOTE: node is only passed in for debugging purposes
            assert init_vw.shape == next_vw.shape, dict(
                msg=("Initial and final state from ScanStateNode must be "
                     "the same."),
                node=node,
                init_shape=init_vw.shape,
                next_shape=next_vw.shape,
                init_vw=init_vw,
                next_vw=next_vw,
            )
            outputs_info[idx] = init_vw.variable

        # ############################## non_sequences ########################
        # ---
        # for now, don't specify any non-sequences and hope that theano
        # can optimize the graph
        non_sequence_vars = []
        non_sequences = []

        # ################################### scan ############################

        def step(*scan_vars):
            # calculate number for each type of scan var
            num_inputs = len(input_sequences)
            num_outputs = len([x for x in outputs_info if x is not None])
            num_non_sequences = len(non_sequences)
            assert len(scan_vars) == (num_inputs
                                      + num_outputs
                                      + num_non_sequences)
            assert num_outputs == len(scan_state_nodes)

            # break down scan vars into appropriate categories
            scan_input_vars = scan_vars[:num_inputs]
            scan_output_vars = scan_vars[num_inputs:num_inputs + num_outputs]
            scan_non_sequences = scan_vars[num_inputs + num_outputs:]

            # setup variables for replacement
            to_replace = []
            for_replace = []
            # input vars
            to_replace += element_input_vars
            for_replace += scan_input_vars
            # non sequences
            to_replace += non_sequence_vars
            for_replace += scan_non_sequences
            # scan state nodes
            for state_idx, scan_output_var in zip(scan_state_order,
                                                  scan_output_vars):
                to_replace.append(scan_state_vars[state_idx])
                for_replace.append(scan_output_var)

            assert len(to_replace) == len(for_replace)

            # perform scan
            new_outputs = utils.deep_clone(
                element_output_vars,
                replace=dict(zip(to_replace, for_replace)),
            )

            final_outputs = list(new_outputs)

            # set next state for recurrent state nodes
            for state_idx, next_idx in zip(scan_state_idxs,
                                           scan_state_next_idxs):
                final_outputs[state_idx] = final_outputs[next_idx]

            return final_outputs

        # edit all outputs of subtree
        results, updates = theano.scan(
            fn=step,
            outputs_info=outputs_info,
            sequences=input_sequences,
            non_sequences=non_sequences,
        )

        # ############################# post-processing #######################

        # scan automatically unwraps lists, so rewrap if needed
        if not isinstance(results, list):
            results = [results]
        assert len(results) == len(element_output_vars)
        result_map = dict(zip(element_output_vars, results))

        # TODO store updates in network (to later be used in new_update_deltas)
        # NOTE: before doing this, look into the effects of manipulating the
        # update deltas of random variables - it might not work as expected

        # FIXME unshuffle dimensions for scan axes

        def transform_shape(old_shape):
            tmp = list(old_shape)
            tmp.insert(scan_axis, None)
            # FIXME what if length is < scan_axis
            return tuple(tmp)

        def transform_output(output_variable):
            # FIXME
            return output_variable

        # FIXME mutate current_variables in network to have new updates

        # create final output with the appropriate result
        network.create_vw(
            name="default",
            variable=transform_output(result_map[element_output.variable]),
            shape=transform_shape(element_output.shape),
        )

reason: import theano


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

__author__ = "Lucas Miguel S Ponce"
__email__ = "lucasmsp@gmail.com"

from ddf_library.utils import generate_info

import functools
import pandas as pd
import numpy as np

def aggregation_stage_1(data, settings):
    nfrag = len(data)
    info1 = settings['schema'][0]
    info1['function'] = [_aggregate, settings]

    from .hash_partitioner import hash_partition

    hash_params = {'columns': settings['groupby'],
                   'nfrag': nfrag, 'schema': [info1]}
    output = hash_partition(data, hash_params)['data']
    settings['intermediate_result'] = True
    return output, settings

def _aggregate(data, params):
    """Perform a partial aggregation."""
    columns = params['groupby']
    operation_list = params['operation']

    if columns is '*':
        columns = list(data.columns)

    operations = _generate_agg_operations(operation_list, columns)
    data = data.groupby(columns).agg(**operations)\
        .reset_index()\
        .reset_index(drop=True)

    return data, params

def aggregation_stage_2(data1, params):
    """Combining the aggregation with other fragment."""
    columns = params['groupby']
    operation_list = params['operation']
    frag = params['id_frag']

    operations = _generate_agg_operations(operation_list, columns, True)
    data1 = data1.groupby(columns).agg(**operations).reset_index()\
        .reset_index(drop=True)

    info = generate_info(data1, frag)
    return data1, info

def _generate_agg_operations(operations_list, groupby, replace=False):
    """
    Used to create a pd.NamedAgg list to be used by Pandas to aggregate.
    Replace option is used in to merge the partial result of each partition.
    """
    operations = dict()
    for col, func, target in operations_list:
        if '*' in col:
            target = groupby

        if replace:
            if func in 'count':
                func = 'sum'
            elif func == 'set':
                func = _merge_set
            elif func == 'list':
                func = _merge_list
            col = target

        else:
            if func == 'list':
                func = list
            elif func == 'set':
                func = _collect_set

        operations[target] = pd.NamedAgg(column=col, aggfunc=func)
    return operations

def _merge_list(series):
    return functools\
        .reduce(lambda x, y: list(x) + list(y), series.tolist())

def _collect_set(x):
    """Generate a set of a group."""
    return list(set(x))

def _merge_set(series):
    """Merge set list."""
    return functools\
        .reduce(lambda x, y: list(set(list(x) + list(y))), series.tolist())

max_stars_repo_name: eubr-bigsea/Compss-Python

max_issues_repo_name: eubr-bigsea/Compss-Python

max_forks_repo_name: eubr-bigsea/Compss-Python


text:
import numpy as np
import theano
import theano.tensor as T
import inspect
import collections
from collections import OrderedDict

from autodiff.context import Context
import autodiff.utils as utils
from autodiff.functions import escape, escaped_call

class Symbolic(object):
    """
    A class that converts a Python function into a symbolic function of Theano
    objects. The symbolic function can be used to compile Theano functions of
    the original (Python) function, its gradient, or a Hessian-vector product.
    """

    def __init__(self,
                 pyfn,
                 context=None,
                 force_floatX=False,
                 borrowable=None,
                 ignore=None,
                 infer_updates=False,
                 escape_on_error=False):
        """
        Arguments
        ---------

        borrow : tuple of objects
            If an object in this tuple is encountered while tracing the
            function, then its symbolic representation will alias that object's
            memory location. This means that *inplace* operations on the Python
            (likely NumPy) object will affect the symbolic function.

        """

        if context is None:
            context = Context(borrowable=utils.as_seq(borrowable, tuple),
                              ignore=utils.as_seq(ignore, tuple),
                              force_floatX=force_floatX,
                              infer_updates=infer_updates,
                              escape_on_error=escape_on_error)
        assert isinstance(context, Context)
        self.context = context

        if isinstance(pyfn, Symbolic):
            pyfn = pyfn.pyfn
        self._pyfn = pyfn

        self._symfn = self.context.recompile(self.pyfn)

    def __get__(self, instance, owner=None):
        """
        Necessary descriptor for decorator compatibility.

        At decoration time, methods have not been bound. However, when bound
        methods are accessed, the __get__ method is called, so we can monitor
        that call and bind the method as necessary.
        """
        if instance is not None:
            method = self.pyfn.__get__(instance, owner)
            self._pyfn = method
            self._symfn = self.context.recompile(self.pyfn)
        return self

    def __call__(self, *args, **kwargs):
        return self.trace(*args, **kwargs)[1]

    @property
    def pyfn(self):
        return self._pyfn

    @property
    def symfn(self):
        return self._symfn

    @property
    def cache(self):
        return self._cache

    @property
    def sym_vars(self):
        return self.context.sym_vars

    @property
    def tags(self):
        return self.context.tags

    def get_symbolic(self, x):
        return self.context.get_symbolic(x)

    def trace(self, *args, **kwargs):
        """
        Call the symbolic function on args and kwargs, returning the symbolic
        result and storing all shadowed variables.
        """
        # clean args and kwargs
        c_args, c_kwargs = utils.clean_int_args(*args, **kwargs)

        # call the symfn
        results = self.symfn(*c_args, **c_kwargs)

        # get a tuple of the symbolic inputs
        # but avoid 'self' and 'cls' bound arguments
        all_args = utils.expandedcallargs(self.symfn, *c_args, **c_kwargs)
        if (inspect.ismethod(self.pyfn) or
           (len(all_args) > 0 and type(all_args[0]) is type)):
            all_args = all_args[1:]

        # store the inputs and outputs so they can be accessed later
        # self.s_inputs = tuple(self.get_symbolic(a) for a in all_args)
        # self.s_outputs = utils.as_seq(results, tuple)
        inputs = tuple(self.get_symbolic(a) for a in all_args)
        # outputs = utils.as_seq(results, tuple)

        return inputs, results

    def get_function_compile_args(self, inputs, outputs):
        """
        Helper function: given the symbolic inputs and outputs,
        return the appropriate arguments for theano.function to compile
        a function.
        """
        return dict(inputs=inputs, outputs=outputs)

    def get_gradient_compile_args(self,
                                  inputs,
                                  outputs,
                                  wrt=None,
                                  reduction=None):
        """
        Helper function: given the symbolic inputs and outputs, as well as
        a theano graph and wrt/reduction info, return the appropriate arguments
        for theano.function to compile a gradient.
        """
        wrt = utils.as_seq(wrt)

        if reduction in ['sum', 'max', 'mean', 'min', 'prod', 'std', 'var']:
            reduction = getattr(theano.tensor, reduction)

        if isinstance(reduction, collections.Callable):
            if 'numpy' in reduction.__module__:
                reduction = getattr(theano.tensor, reduction.__name__)
            outputs = [reduction(o) if o.ndim > 0 else o for o in outputs]

        if any([o.ndim != 0 for o in outputs]):
            raise TypeError('Gradient requires either scalar outputs or a '
                            'reduction that returns a scalar.')

        # get wrt variables. If none were specified, use inputs.
        if len(wrt) == 0:
            wrt = [i for i in inputs]
        else:
            wrt = [self.get_symbolic(w) for w in wrt]

        grads = utils.flatten([T.grad(o, wrt=wrt) for o in outputs])

        return dict(inputs=inputs, outputs=utils.as_seq(grads, tuple))

    def get_hessian_vector_compile_args(self,
                                        inputs,
                                        outputs,
                                        wrt=None,
                                        reduction=None):
        """
        Helper function: given the symbolic inputs and outputs, as well as
        a theano graph and wrt/reduction/vectors info, return the appropriate
        argumentsfor theano.function to compile a Hessian-vector product.
        """
        wrt = utils.as_seq(wrt)

        if reduction in ['sum', 'max', 'mean', 'min', 'prod', 'std', 'var']:
            reduction = getattr(theano.tensor, reduction)

        if isinstance(reduction, collections.Callable):
            if 'numpy' in reduction.__module__:
                reduction = getattr(theano.tensor, reduction.__name__)
            outputs = [reduction(o) if o.ndim > 0 else o for o in outputs]

        if any([o.ndim != 0 for o in outputs]):
            raise TypeError('Gradient requires either scalar outputs or a '
                            'reduction that returns a scalar.')

        # get wrt variables. If none were specified, use inputs.
        if len(wrt) == 0:
            wrt = [i for i in inputs]
        else:
            wrt = [self.get_symbolic(w) for w in wrt]

        grads = utils.flatten([T.grad(o, wrt=wrt) for o in outputs])

        sym_vectors = tuple(T.TensorType(
            dtype=w.dtype, broadcastable=[False] * w.ndim)()
            for w in wrt)
        hessian_vectors = utils.as_seq(T.Rop(grads, wrt, sym_vectors), tuple)

        return dict(inputs=inputs + sym_vectors, outputs=hessian_vectors)

    def compile(self,
                function=False,
                gradient=False,
                hessian_vector=False,
                inputs=None,
                outputs=None,
                wrt=None,
                reduction=None,
                allow_input_downcast=True):

        assert isinstance(function, bool)
        assert isinstance(gradient, bool)
        assert isinstance(hessian_vector, bool)

        if not (function or gradient or hessian_vector):
            raise ValueError(
                'At least one of `function`, `gradient`, or `hessian_vector` '
                'must be True when calling `compile()`.')

        sym_inputs = tuple(
          self.get_symbolic(i) for i in utils.as_seq(inputs))
        sym_outputs = tuple(
          self.get_symbolic(o) for o in utils.as_seq(outputs))

        fn_inputs = sym_inputs
        fn_outputs = ()

        if function:
            fn_args = self.get_function_compile_args(inputs=sym_inputs,
                                                     outputs=sym_outputs)
            fn_outputs += fn_args['outputs']

        if gradient:
            g_args = self.get_gradient_compile_args(inputs=sym_inputs,
                                                    outputs=sym_outputs,
                                                    wrt=wrt,
                                                    reduction=reduction)
            fn_outputs += g_args['outputs']

        if hessian_vector:
            hv_args = self.get_hessian_vector_compile_args(inputs=sym_inputs,
                                                           outputs=sym_outputs,
                                                           wrt=wrt,
                                                           reduction=reduction)
            fn_inputs = hv_args['inputs']
            fn_outputs += hv_args['outputs']

        if len(fn_outputs) == 1:
            fn_outputs = fn_outputs[0]

        new_inputs = tuple(i.type() for i in fn_inputs)
        givens = dict(zip(fn_inputs, new_inputs))

        if self.context.infer_updates:
            updates = self.context.updates
        else:
            updates = collections.OrderedDict()

        fn = theano.function(inputs=new_inputs,
                             outputs=fn_outputs,
                             givens=givens,
                             updates=updates,
                             on_unused_input='ignore',
                             allow_input_downcast=allow_input_downcast)

        return fn

    def compile_function(self,
                         inputs=None,
                         outputs=None,
                         allow_input_downcast=True):
        """
        Based on traced variables, compile a Theano function of the inputs that
        returns the outputs.
        """
        fn = self.compile(
            function=True,
            inputs=inputs,
            outputs=outputs,
            allow_input_downcast=allow_input_downcast)
        return fn

    def compile_gradient(self,
                         inputs=None,
                         outputs=None,
                         wrt=None,
                         reduction=None,
                         allow_input_downcast=True):
        """
        Based on traced variables, compile a Theano function of the
        inputs that returns the gradient of the outputs with respect to wrt.
        If wrt is None, it is assumed to be all of the inputs. A reduction may
        be specified (since gradients are defined with respect to scalars); if
        None is supplied, it is assumed to be 'sum'.
        """
        fn = self.compile(
            gradient=True,
            inputs=inputs,
            outputs=outputs,
            wrt=wrt,
            reduction=reduction,
            allow_input_downcast=allow_input_downcast)
        return fn

    def compile_function_gradient(self,
                                  inputs=None,
                                  outputs=None,
                                  wrt=None,
                                  reduction=None,
                                  allow_input_downcast=True):
        """
        Based on traced variables, compile a Theano function of the
        inputs that returns both the outputs and the gradient of the outputs
        with respect to wrt. If wrt is None, it is assumed to be all of the
        inputs. A reduction may be specified (since gradients are defined with
        respect to scalars); if None is supplied, it is assumed to be 'sum'.
        """

        fn = self.compile(
            function=True,
            gradient=True,
            inputs=inputs,
            outputs=outputs,
            wrt=wrt,
            reduction=reduction,
            allow_input_downcast=allow_input_downcast)
        return fn

class Tracer(Symbolic):
    """
    A Symbolic class for tracing variables through multiple functions.
    """

    def __init__(self,
                 context=None,
                 force_floatX=False,
                 infer_updates=False,
                 borrowable=None,
                 ignore=None,
                 escape_on_error=False):
        super(Tracer, self).__init__(pyfn=lambda: None,
                                     context=context,
                                     force_floatX=force_floatX,
                                     infer_updates=infer_updates,
                                     borrowable=borrowable,
                                     ignore=ignore,
                                     escape_on_error=escape_on_error)

    def trace(self, pyfn, *args, **kwargs):
        symbolic = Symbolic(pyfn=pyfn, context=self.context)
        return symbolic.trace(*args, **kwargs)[1]

class Function(Symbolic):
    """
    A Symbolic tracer which is specialized for a specific function, passed at
    initialization.
    """

    def __init__(self,
                 pyfn,
                 context=None,
                 force_floatX=False,
                 borrowable=None,
                 ignore=None,
                 infer_updates=False,
                 escape_on_error=False,
                 use_cache=True):
        super(Function, self).__init__(pyfn=pyfn,
                                       context=context,
                                       force_floatX=force_floatX,
                                       borrowable=borrowable,
                                       ignore=ignore,
                                       infer_updates=infer_updates,
                                       escape_on_error=escape_on_error)

        self._cache = dict()
        self.use_cache = use_cache

    def __call__(self, *args, **kwargs):
        all_args = utils.expandedcallargs(self.symfn, *args, **kwargs)
        if (inspect.ismethod(self.pyfn) or
           (len(all_args) > 0 and type(all_args[0]) is type)):
            all_args = all_args[1:]

        key = tuple(
            (np.asarray(a).ndim, np.asarray(a).dtype) for a in all_args)
        if key not in self.cache or not self.use_cache:
            self.context.reset()
            inputs, outputs = self.trace(*args, **kwargs)
            self.cache[key] = self.get_theano_function(inputs, outputs)
        fn = self.cache[key]
        return fn(*all_args)

    def get_theano_function(self, inputs, outputs):
        fn = self.compile_function(inputs=inputs, outputs=outputs)
        return fn

class Gradient(Function):
    def __init__(self,
                 pyfn,
                 wrt=None,
                 reduction=None,
                 force_floatX=False,
                 infer_updates=False,
                 borrowable=None,
                 ignore=None,
                 escape_on_error=False,
                 context=None,
                 use_cache=True):
        super(Gradient, self).__init__(pyfn=pyfn,
                                       force_floatX=force_floatX,
                                       borrowable=borrowable,
                                       ignore=ignore,
                                       infer_updates=infer_updates,
                                       context=context,
                                       escape_on_error=escape_on_error,
                                       use_cache=use_cache)
        self.wrt = utils.as_seq(wrt, tuple)
        self.reduction = reduction

    def get_theano_function(self, inputs, outputs):
        fn = self.compile_gradient(inputs=inputs,
                                   outputs=outputs,
                                   wrt=self.wrt,
                                   reduction=self.reduction)
        return fn

class HessianVector(Gradient):

    def __call__(self, *args, **kwargs):
        if 'vectors' in kwargs:
            vectors = kwargs.pop('vectors')
        else:
            raise ValueError(
                'HessianVector must be called with the keyword \'vectors\'.')
        vectors = utils.as_seq(vectors, tuple)

        all_args = utils.expandedcallargs(self.symfn, *args, **kwargs)

        key = tuple(np.asarray(a).ndim for a in all_args)
        if key not in self.cache or not self.use_cache:
            self.context.reset()
            inputs, outputs = self.trace(*args, **kwargs)
            self.cache[key] = self.get_theano_function(inputs, outputs)
        fn = self.cache[key]

        if len(self.wrt) > 0 and len(vectors) != len(self.wrt):
            raise ValueError('Expected {0} items in `vectors`; received '
                             '{1}.'.format(len(self.wrt), len(vectors)))
        elif len(self.wrt) == 0 and len(vectors) != len(inputs):
            raise ValueError('Expected {0} items in `vectors`; received '
                             '{1}.'.format(len(inputs), len(vectors)))

        return fn(*(all_args + vectors))

    def get_theano_function(self, inputs, outputs):
        fn = self.compile(hessian_vector=True,
                          inputs=inputs,
                          outputs=outputs,
                          wrt=self.wrt,
                          reduction=self.reduction)
        return fn

class VectorArg(object):

    def __init__(self,
                 pyfn,
                 init_args=None,
                 init_kwargs=None,
                 context=None,
                 force_floatX=False,
                 borrowable=None,
                 ignore=None,
                 infer_updates=False,
                 escape_on_error=False,
                 function=False,
                 gradient=False,
                 hessian_vector=False):

        if isinstance(pyfn, Symbolic):
            pyfn = pyfn.pyfn
        self.pyfn = pyfn

        init_args = utils.as_seq(init_args, tuple)
        init_kwargs = utils.as_seq(init_kwargs, dict)

        self.init_args = utils.expandedcallargs(
            pyfn, *init_args, **init_kwargs)

        def wrapped_function(vector):
            return pyfn(*escaped_call(self.args_from_vector, vector))

        def wrapper(*args, **kwargs):
            vector = self.vector_from_args(args, kwargs)
            v_args = self.args_from_vector(vector)
            return vector, pyfn(*v_args)

        symbolic = Symbolic(pyfn=wrapper,
                            context=context,
                            force_floatX=force_floatX,
                            infer_updates=infer_updates,
                            borrowable=borrowable,
                            ignore=ignore,
                            escape_on_error=escape_on_error)

        _, (sym_vector, result) = symbolic.trace(*init_args, **init_kwargs)

        fn = symbolic.compile(function=function,
                              gradient=gradient,
                              hessian_vector=hessian_vector,
                              inputs=sym_vector,
                              outputs=result)
        self.fn = fn

    def __call__(self, *args, **kwargs):
        return self.fn(*args, **kwargs)

    def vector_from_args(self, args, kwargs):
        if len(args) + len(kwargs) > 1:
            all_args = utils.expandedcallargs(self.pyfn, *args, **kwargs)
            return np.concatenate([np.asarray(a).flatten() for a in all_args])
        elif len(args) > 0:
            return np.asarray(args[0]).flatten()
        elif len(kwargs) > 0:
            return np.asarray(list(kwargs.values())[0]).flatten()
        else:
            return None

    def args_from_vector(self, vector):
        new_args = []
        idx = 0
        for arg in escape(self.init_args):
            # new_args.append(vector[idx: idx + arg.size].reshape(*arg.shape))
            # avoid using arg.size or prod(arg.shape) because Theano's prod
            # doesn't support R-op
            new_args.append(
                vector[idx: idx + safesize(arg)].reshape(arg.shape))
            idx += safesize(arg)
        return new_args

def safesize(arg):
    if arg.ndim == 0:
        size = 1
    else:
        size = arg.shape[0]
        for i in list(range(1, arg.ndim)):
            size *= arg.shape[i]
    return size

max_stars_repo_name: gwtaylor/pyautodiff

max_issues_repo_name: gwtaylor/pyautodiff

max_forks_repo_name: gwtaylor/pyautodiff

reason: import numpy,import theano


text:
import cortex
import numpy as np
import tempfile

from cortex import db, dataset

subj, xfmname, nverts, volshape = "S1", "fullhead", 304380, (31,100,100)

def test_braindata():
    vol = np.random.randn(*volshape)
    tf = tempfile.TemporaryFile(suffix='.png')
    mask = db.get_mask(subj, xfmname, "thick")

    data = dataset.Volume(vol, subj, xfmname, cmap='RdBu_r', vmin=0, vmax=1)
    # quickflat.make_png(tf, data)
    mdata = data.masked['thick']
    assert len(mdata.data) == mask.sum()
    assert np.allclose(mdata.volume[:, mask], mdata.data)

def test_dataset():
    vol = np.random.randn(*volshape)
    stack = (np.ones(volshape[::-1])*np.linspace(0, 1, volshape[0])).T
    mask = db.get_mask(subj, xfmname, "thick")

    ds = dataset.Dataset(randvol=(vol, subj, xfmname), stack=(stack, subj, xfmname))
    ds.append(thickstack=ds.stack.masked['thick'])
    tf = tempfile.NamedTemporaryFile(suffix=".hdf")
    ds.save(tf.name)

    ds = dataset.Dataset.from_file(tf.name)
    assert len(ds['thickstack'].data) == mask.sum()
    assert np.allclose(ds['stack'].data[mask], ds['thickstack'].data)
    return ds

def test_findmask():
    vol = np.random.rand(10, *volshape)
    mask = db.get_mask(subj, xfmname, "thin")
    ds = dataset.Volume(vol[:, mask], subj, xfmname)
    assert np.allclose(ds.volume[:, mask], vol[:, mask])
    return ds

def test_rgb():
    red, green, blue, alpha = [np.random.randn(*volshape) for _ in range(4)]

    rgb = dataset.VolumeRGB(red, green, blue, subj, xfmname)
    assert rgb.volume.shape == tuple([1] + list(volshape) + [4])
    assert rgb.volume.dtype == np.uint8
    assert rgb.volume[..., 3].max() > 0

    rgba = dataset.VolumeRGB(red, green, blue, subj, xfmname, alpha=alpha)
    assert rgba.volume.shape == tuple([1] + list(volshape) + [4])

    data = dataset.Volume.random(subj, xfmname)
    assert data.raw.volume.shape == tuple([1] + list(volshape) + [4])
    data.raw.to_json()

    red, green, blue, alpha = [np.random.randn(nverts) for _ in range(4)]

    rgb = dataset.VertexRGB(red, green, blue, subj)
    assert rgb.vertices.shape == (1, nverts, 4)
    assert rgb.vertices.dtype == np.uint8
    assert rgb.vertices[..., 3].max() > 0

    rgba = dataset.VertexRGB(red, green, blue, subj, alpha=alpha)
    assert rgba.vertices.shape == (1, nverts, 4)

    data = dataset.Vertex.random(subj)
    assert data.raw.vertices.shape == (1, nverts, 4)
    data.raw.to_json()

def test_2D():
    d1 = cortex.Volume.random(subj, xfmname)
    d2 = cortex.Volume.random(subj, xfmname).masked['thick']
    twod = cortex.Volume2D(d1, d2)
    cortex.Volume2D(d1.data, d2.data, subject=subj, xfmname=xfmname, vmin=0, vmax=2, vmin2=1)
    twod.to_json()
    return twod

def test_braindata_hash():
    d = cortex.Volume.random(subj, xfmname)
    hash(d)

def test_dataset_save():
    tf = tempfile.NamedTemporaryFile(suffix=".hdf")
    mrand = np.random.randn(2, *volshape)
    rand = np.random.randn(*volshape)
    ds = cortex.Dataset(test=(mrand, subj, xfmname))
    ds.append(twod=cortex.Volume2D(rand, rand, subj, xfmname))
    ds.append(rgb =cortex.VolumeRGB(rand, rand, rand, subj, xfmname))
    ds.append(vert=cortex.Vertex.random(subj))
    ds.save(tf.name)

    ds = cortex.load(tf.name)
    assert isinstance(ds.test, cortex.Volume)
    assert ds.test.data.shape == mrand.shape
    assert isinstance(ds.twod, cortex.Volume2D)
    assert ds.twod.dim1.data.shape == rand.shape
    assert ds.twod.dim2.data.shape == rand.shape
    assert ds.rgb.volume.shape == tuple([1] + list(volshape) + [4])
    assert isinstance(ds.vert, cortex.Vertex)

def test_mask_save():
    tf = tempfile.NamedTemporaryFile(suffix=".hdf")
    ds = cortex.Dataset(test=(np.random.randn(*volshape), subj, xfmname))
    ds.append(masked=ds.test.masked['thin'])
    data = ds.masked.data
    ds.save(tf.name)

    ds = cortex.load(tf.name)
    assert ds.masked.shape == volshape
    assert np.allclose(ds.masked.data, data)

def test_overwrite():
    tf = tempfile.NamedTemporaryFile(suffix=".hdf")
    ds = cortex.Dataset(test=(np.random.randn(*volshape), subj, xfmname))
    ds.save(tf.name)

    ds.save()
    assert ds.test.data.shape == volshape

def test_pack():
    tf = tempfile.NamedTemporaryFile(suffix=".hdf")
    ds = cortex.Dataset(test=(np.random.randn(*volshape), subj, xfmname))
    ds.save(tf.name, pack=True)

    ds = cortex.load(tf.name)
    pts, polys = cortex.db.get_surf(subj, "fiducial", "lh")
    dpts, dpolys = ds.get_surf(subj, "fiducial", "lh")
    assert np.allclose(pts, dpts)

    overlay_db = cortex.db.get_overlay(subj, None, modify_svg_file=False)
    rois_db = overlay_db.rois.labels.elements.keys()
    # keep the temporary file object in memory to avoid the file being deleted
    temp_file = ds.get_overlay(subj, "rois")
    overlay_ds = cortex.db.get_overlay(subj, temp_file.name,
                                       modify_svg_file=False)
    rois_ds = overlay_ds.rois.labels.elements.keys()
    assert rois_db == rois_ds

    xfm = cortex.db.get_xfm(subj, xfmname)
    assert np.allclose(xfm.xfm, ds.get_xfm(subj, xfmname).xfm)

def test_map():
    dv = cortex.Volume.random(subj, xfmname)
    dv.map("nearest")

def test_convertraw():
    ds = cortex.Dataset(test=(np.random.randn(*volshape), subj, xfmname))
    ds.test.raw

def test_vertexdata_copy():
    vd = cortex.Vertex(np.random.randn(nverts), subj)
    vdcopy = vd.copy(vd.data)
    assert np.allclose(vd.data, vdcopy.data)

def test_vertexdata_set():
    vd = cortex.Vertex(np.random.randn(nverts), subj)
    newdata = np.random.randn(nverts)
    vd.data = newdata
    assert np.allclose(newdata, vd.data)

def test_vertexdata_index():
    vd = cortex.Vertex(np.random.randn(10, nverts), subj)
    assert np.allclose(vd[0].data, vd.data[0])

def test_vertex_rgb_movie():
    r = g = b = np.random.randn(nverts)
    rgb = cortex.VertexRGB(r, g, b, subj)

def test_volumedata_copy():
    v = cortex.Volume(np.random.randn(*volshape), subj, xfmname)
    vc = v.copy(v.data)
    assert np.allclose(v.data, vc.data)

def test_volumedata_copy_with_custom_mask():
    mask = cortex.get_cortical_mask(subj, xfmname, "thick")
    mask[16] = True
    nmask = mask.sum()
    data = np.random.randn(nmask)
    v = cortex.Volume(data, subj, xfmname, mask=mask)
    vc = v.copy(v.data)
    assert np.allclose(v.data, vc.data)

def test_int64_in_dataviewrgb():
    data = np.arange(np.product(volshape)).reshape(volshape, order='C')
    view = cortex.VolumeRGB(data, data + 1, data + 2, subject=subj,
                            xfmname=xfmname)
    cortex.quickshow(view)

    data = np.arange(nverts)
    view = cortex.VertexRGB(data, data + 1, data + 2, subject=subj)
    cortex.quickshow(view)

def test_vmin_none_in_dataview2d():
    data = np.arange(np.product(volshape)).reshape(volshape, order='C')
    view = cortex.Volume2D(data, data + 1, subject=subj, xfmname=xfmname)
    cortex.quickshow(view)

    data = np.arange(nverts)
    view = cortex.Vertex2D(data, data + 1, subject=subj)
    cortex.quickshow(view)

def test_dataset_operators():
    vol = cortex.Volume.random(subj, xfmname)
    array = np.random.randn(*volshape)

    assert np.allclose(vol.data + array, (vol + array).data)
    assert np.allclose(vol.data - array, (vol - array).data)
    assert np.allclose(vol.data * array, (vol * array).data)
    assert np.allclose(vol.data // array, (vol // array).data) # floordiv
    assert np.allclose(vol.data / array, (vol / array).data) # truediv
    # numpy doesn't like fractional powers of negative numbers
    assert np.allclose(vol.data ** array, (vol ** array).data, equal_nan=True)
    assert np.allclose(-vol.data, (-vol).data)
    assert np.allclose(abs(vol.data), abs(vol).data)

max_stars_repo_name: gallantlab/pycortex

max_issues_repo_name: gallantlab/pycortex

max_forks_repo_name: gallantlab/pycortex