text:
#!/usr/bin/env python

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

# Top level for EM supervoxel classifier that learns merge / no merge classification based on some feature set.
# Featured Region Adjacency Graph (dpFRAG.py) creates the graph of neighboring supervoxels and their features.
#
# There is no iteration in this procedure, typically referred to as "flat learning" in the context of merging
#   supervoxels.
# After investigating flat learning, added a simple iterative procedure that agglomerates some percentage of the most
#   confident mergers and continues for set number of iterations.

import numpy as np
import time
import argparse
import os
#import sys
import importlib
from io import BytesIO

from configobj import ConfigObj, flatten_errors
from validate import Validator, ValidateError
import dill
from io import StringIO
from sklearn.preprocessing import scale
from sklearn.utils import shuffle
from scipy import ndimage as nd

from sklearn.neighbors import KNeighborsClassifier
from sklearn.svm import SVC
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier, AdaBoostClassifier
from sklearn.naive_bayes import GaussianNB
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis
from sklearn.linear_model import LogisticRegression

# xxx - how to make imports as optional? make plotting class?
from matplotlib import pylab as pl
from matplotlib import pyplot as plt
#import matplotlib as mpl
#from matplotlib import colors
#from cycler import cycler

from dpLoadh5 import dpLoadh5
# made dpFRAG import dynamic in init
#from dpFRAG import dpFRAG
#from dpFRAGc import dpFRAG
dpFRAG = []
from metrics import pixel_error_fscore
from utils import print_cpu_info_linux

class dpSupervoxelClassifier():

    # Constants
    LIST_ARGS = ['test_chunks', 'label_subgroups', 'label_subgroups_out', 'iterate_merge_perc',
        'thresholds', 'threshold_subgroups']

    def __init__(self, args):

        # save command line arguments from argparse, see definitions in main or run with --help
        for k, v in vars(args).items():
            # do not override any values that are already set as a method of allowing inherited classes to specify
            if hasattr(self,k): continue
            if type(v) is list and k not in self.LIST_ARGS:
                if len(v)==1:
                    setattr(self,k,v[0])  # save single element lists as first element
                elif type(v[0]) is int:   # convert the sizes and offsets to numpy arrays
                    setattr(self,k,np.array(v,dtype=np.int32))
                else:
                    setattr(self,k,v)   # store other list types as usual (floats)
            else:
                setattr(self,k,v)

        # save the command line argument dict as a string
        out = StringIO(); print( vars(args), file=out )
        self.arg_str = out.getvalue(); out.close()

        print('dpSupervoxelClassifier: config file ''%s''' % self.cfgfile)
        # retrieve / save options from ini files, see definitions dpSupervoxelClassifier.ini
        opts = dpSupervoxelClassifier.get_options(self.cfgfile)
        d = vars(self)
        for k, v in opts.items():
            # do not import if have a "True" (non-empty) value from command line
            if k in d and d[k]: continue
            if type(v) is list and k not in self.LIST_ARGS:
                if len(v)==1:
                    setattr(self,k,v[0])  # save single element lists as first element
                elif len(v)>0 and type(v[0]) is int:   # convert the sizes and offsets to numpy arrays
                    setattr(self,k,np.array(v,dtype=np.int32))
                else:
                    setattr(self,k,v)   # store other list types as usual (floats, empties)
            else:
                setattr(self,k,v)

        # save the command line argument dict as a string
        out = StringIO(); print( opts, file=out )
        self.ini_str = out.getvalue(); out.close()

        # Options / Inits

        # dynamic imports
        dpFRAGl = importlib.import_module('dpFRAGc') if self.useFRAGc else importlib.import_module('dpFRAG')
        globals().update({'dpFRAG':dpFRAGl.dpFRAG})

        # these are so standard cubeIter inputs can be used from command line to override from .ini
        if (self.chunk >= 0).all():
            self.chunk_range_beg = self.chunk
            self.chunk_range_end = []
            self.offset_list = self.offset
            self.size_list = self.size

        self.doplots = (self.show_plots or self.export_plots)

        # default is to have sklearn calculate the priors
        self.priors = None
        if self.merge_prior > 0: self.priors = np.array([1-self.merge_prior, self.merge_prior],dtype=np.double)

        # keep the list of merge percentages, use original to store current iteration value
        self.iterate_merge_perc_list = self.iterate_merge_perc

        # initialize for "chunkrange" or "chunklist" mode if these parameters are not empty.
        # this code was modified from that in the em data parser for cuda-convnets2.
        self.chunk_range_beg = self.chunk_range_beg.reshape(-1,3); self.use_chunk_range = False
        self.nchunk_list = self.chunk_range_beg.shape[0]
        self.nchunks = self.nchunk_list

        # offsets are either per chunk or per range, depending on above mode (whether chunk_range_end empty or not)
        self.offset_list = self.offset_list.reshape(-1,3)   # list must have multiple of 3 elements for volumes
        if self.offset_list.shape[0] > 1:
            assert( self.offset_list.shape[0] == self.nchunk_list )
        else:
            self.offset_list = np.ones_like(self.chunk_range_beg)*self.offset_list

        # sizes are either per chunk or per range, depending on above mode (whether chunk_range_end empty or not).
        self.size_list = self.size_list.reshape(-1,3)   # list must have multiple of 3 elements for volumes
        if self.size_list.shape[0] > 1:
            assert( self.size_list.shape[0] == self.nchunk_list )
        else:
            self.size_list = np.ones_like(self.chunk_range_beg)*self.size_list

        if len(self.chunk_range_end) > 0:
            # "chunkrange" mode, chunks are selected based on defined beginning and end of ranges in X,Y,Z
            # range is open ended (python-style, end is not included in range)
            self.chunk_range_end = self.chunk_range_end.reshape(-1,3);
            assert( self.chunk_range_end.shape[0] == self.nchunk_list )
            self.chunk_range_rng = self.chunk_range_end - self.chunk_range_beg
            assert( (self.chunk_range_rng >= 0).all() )     # some bad ranges

            # this allows for iterating in chunk_range mode with sizes greater than chunksize
            # just read the h5info for the rawfile to get the chunksize
            rawh5 = dpLoadh5.readInith5(self.rawfile, self.raw_dataset, [0,0,0], [0,0,0], [128,128,128], 'uint8')
            self.chunksize = rawh5.chunksize
            self.chunkscale = np.ones_like(self.size_list)
            # only allow this if the sizes are a multiple of chunksize
            if (self.size_list % self.chunksize == 0).all():
                scale = self.size_list // self.chunksize
                # only allow this if the ranges are a multiple of the sizes
                if (self.chunk_range_rng % scale == 0).all():
                    self.chunkscale = scale
            self.chunk_range_rng //= self.chunkscale

            self.chunk_range_size = self.chunk_range_rng.prod(axis=1)
            self.chunk_range_cumsize = np.concatenate((np.zeros((1,),dtype=self.chunk_range_size.dtype),
                self.chunk_range_size.cumsum()))
            self.chunk_range_nchunks = self.chunk_range_cumsize[-1]
            self.use_chunk_range = True; self.nchunks = self.chunk_range_nchunks

        # print out all initialized variables in verbose mode
        if self.dpSupervoxelClassifier_verbose:
            # print out info for chunklist / chunkrange modes so that input data is logged
            print(('Using %d ' % self.nchunk_list) + ('ranges' if self.use_chunk_range else 'chunks') + ':')
            fh = BytesIO()
            if self.use_chunk_range:
                np.savetxt(fh, np.concatenate((np.arange(self.nchunk_list).reshape((self.nchunk_list,1)),
                    self.chunk_range_beg, self.chunk_range_end, self.chunk_range_size.reshape((self.nchunk_list,1)),
                    self.size_list, self.offset_list), axis=1),
                    fmt='\t(%d) range %d %d %d to %d %d %d (%d chunks), size %d %d %d, offset %d %d %d',
                    delimiter='', newline='\n', header='', footer='', comments='')
            else:
                np.savetxt(fh, np.concatenate((np.arange(self.nchunk_list).reshape((self.nchunk_list,1)),
                    self.chunk_range_beg, self.size_list, self.offset_list), axis=1),
                    fmt='\t(%d) chunk %d %d %d, size %d %d %d, offset %d %d %d',
                    delimiter='', newline='\n', header='', footer='', comments='')
            cstr = fh.getvalue(); fh.close(); print(cstr.decode('UTF-8'))
            print('Test chunks: %s' % ' '.join([str(x) for x in self.test_chunks]))

            #print('dpSupervoxelClassifier, verbose mode:\n'); print(vars(self))
            print_cpu_info_linux() # for debugging runtime variance on biowulf

        # inits for iterative prior mode
        self.iterative_mode = (self.iterate_count > 0)
        self.iterative_frag = [None] * self.nchunks
        self.iterative_mode_count = 0
        if self.outfile:
            if len(self.iterate_save_ranges) > 0:
                assert(len(self.iterate_save_ranges) % 3 == 0) # must specify as python-style ranges
                n = len(self.iterate_save_ranges)//3
                self.iterate_save_mask = np.zeros((self.iterate_count,), dtype=bool)
                for i in range(n):
                    start, stop, step = self.iterate_save_ranges[3*i:3*i+3]
                    self.iterate_save_mask[range(start,stop,step)] = 1
                self.iterate_save_mask[stop:] = 1
            else:
                self.iterate_save_mask = np.ones((self.iterate_count,), dtype=bool)
        else:
            self.iterate_save_mask = np.zeros((self.iterate_count,), dtype=bool)

        if self.iterative_mode:
            # expand the list out to the number of iterations by repeating the last merge percentage
            tmp = np.zeros((self.iterate_count,), dtype=np.double)
            tmp[:len(self.iterate_merge_perc_list)] = self.iterate_merge_perc_list
            tmp[len(self.iterate_merge_perc_list):] = self.iterate_merge_perc_list[-1]
            self.iterate_merge_perc_list = tmp

        # xxx - clean this up at some point, basically hack for legacy metric comparison scripts
        if len(self.threshold_subgroups) == 0:
            if self.iterative_mode:
                self.threshold_subgroups = np.arange(1,self.iterate_count+1,dtype=np.double)
            else:
                self.threshold_subgroups = np.sort(self.thresholds)[::-1]
        else:
            self.threshold_subgroups = np.sort(self.threshold_subgroups)[::-1]
            if self.iterative_mode:
                assert( self.iterate_count == len(self.threshold_subgroups) )
            else:
                assert( len(self.thresholds) == len(self.threshold_subgroups) )

        # new mode where saving classifiers separately for iterative mode
        if self.classifierout and self.iterative_mode:
            [self.classifierout_path, tmp] = os.path.split(self.classifierout)
            [self.classifierout_name, self.classifierout_ext] = os.path.splitext(tmp)
            # make a subdirectory of the same name without ext and put all the iterative classifier dills into it
            self.classifierout_path = os.path.join(self.classifierout_path,self.classifierout_name)
            os.makedirs(self.classifierout_path, exist_ok=True)

        # for loading new mode where each classifier stored in separate dill file.
        # if classifierin is not a directory then use old mode where all classifiers stored in a single file.
        if self.classifierin and self.iterative_mode and not os.path.isfile(self.classifierin):
            [tmp, self.classifierin_name] = os.path.split(self.classifierin)

        # dpFRAG features are now selectable as sets
        #self.FEATURES, self.FEATURES_NAMES, self.nfeatures = dpFRAG.make_FEATURES()
        d = dpFRAG.make_features(self.feature_set, self.has_ECS)
        for k in ['features','features_names','nfeatures']: setattr(self,k,d[k])

        # xxx - make this more systematic for other libraries?
        #   typically using anaconda, for which this should work since built on mkl
        try:
            import mkl
            print('Setting mkl num_threads to %d' % (self.nthreads,))
            mkl.set_num_threads(self.nthreads)
        except ImportError:
            pass

    def train(self):

        if self.dpSupervoxelClassifier_verbose: print('\nTRAIN')

        if self.trainin and (not self.classifierin or self.doplots):
            if self.dpSupervoxelClassifier_verbose:
                print('Loading training data')
            with open(self.trainin, 'rb') as f: data = dill.load(f)
            target = data['target']; fdata = data['data']
            ntargets = target.size; nfeatures = fdata.shape[1]
            assert( nfeatures == self.nfeatures )

        elif not self.classifierin:
            #dict_keys(['feature_names', 'DESCR', 'target_names', 'target', 'data'])
            nalloc = self.nchunks*self.nalloc_per_chunk
            nfeatures = self.nfeatures
            target = np.zeros((nalloc,), dtype=np.int64)
            fdata = np.zeros((nalloc,nfeatures), dtype=np.double)

            # accumulate training data from all training chunks
            cnt_targets = 0; ntargets = np.zeros((self.nchunks,),dtype=np.int64)
            for chunk in range(self.nchunks):
                cchunk, chunk_list_index, chunk_range_index = self.get_chunk_inds(chunk)
                offset = self.offset_list[chunk_list_index,:]; size = self.size_list[chunk_list_index,:]

                if chunk_list_index in self.test_chunks: continue
                print('Appending training data for chunk %d,%d,%d' % tuple(cchunk.tolist()))

                if self.iterative_mode:
                    if self.iterative_frag[chunk] is None:
                        frag = dpFRAG.makeBothFRAG(self.labelfile, cchunk, size, offset,
                            [self.probfile, self.probaugfile], [self.rawfile, self.rawaugfile],
                            self.raw_dataset, self.gtfile, self.outfile, self.label_subgroups,
                            ['agglomeration_training'], progressBar=self.progress_bar, feature_set=self.feature_set,
                            has_ECS=self.has_ECS, chunk_subgroups=self.chunk_subgroups, no_agglo_ECS=self.no_agglo_ECS,
                            pad_prob_svox_perim = not self.prob_svox_context, neighbor_only=self.neighbor_only,
                            verbose=self.dpSupervoxelClassifier_verbose)
                        frag.isTraining = True; self.iterative_frag[chunk] = frag
                    else:
                        frag = self.iterative_frag[chunk]
                else:
                    frag = dpFRAG.makeTrainingFRAG(self.labelfile, cchunk, size, offset,
                        [self.probfile, self.probaugfile], [self.rawfile, self.rawaugfile],
                        self.raw_dataset, self.gtfile, self.label_subgroups, feature_set=self.feature_set,
                        has_ECS=self.has_ECS, chunk_subgroups=self.chunk_subgroups, neighbor_only=self.neighbor_only,
                        pad_prob_svox_perim = not self.prob_svox_context, progressBar=self.progress_bar,
                        no_agglo_ECS=self.no_agglo_ECS, verbose=self.dpSupervoxelClassifier_verbose)
                frag.createFRAG(update = self.iterative_mode)
                #frag.createFRAG(update = False)
                data = frag.createDataset()
                ntargets[chunk] = data['target'].shape[0]
                target[cnt_targets:cnt_targets+ntargets[chunk]] = data['target']
                fdata[cnt_targets:cnt_targets+ntargets[chunk],:] = data['data']
                cnt_targets += ntargets[chunk]
            target = target[:cnt_targets]; fdata = fdata[:cnt_targets,:]

            if self.trainout:
                #dict_keys(['feature_names', 'DESCR', 'target_names', 'target', 'data'])
                descr = 'Training data from dpFRAG.py with command line:\n' + self.arg_str
                descr = ('With ini file "%s":\n' % (self.cfgfile,)) + self.ini_str
                data['data'] = fdata; data['target'] = target; data['DESCR'] = descr
                with open(self.trainout, 'wb') as f: dill.dump(data, f)

        if not self.classifierin or self.doplots:
            # everyone wants to be norml
            sdata = scale(fdata)   # normalize for the classifiers
            # shuffle data for training as well
            target, sdata = shuffle(target, sdata)

        if self.classifierin:
            if self.dpSupervoxelClassifier_verbose:
                print('\nLoading classifier:'); t = time.time()

            with open(self.classifierin, 'rb') as f: d = dill.load(f)
            self.clf = d['classifier'];
        else:
            if self.dpSupervoxelClassifier_verbose:
                print('\nTraining classifier %s with %d examples and %d features:' % (self.classifier,
                    cnt_targets, nfeatures)); t = time.time()

            # train a classifier
            if self.classifier == 'lda':
                self.clf = LinearDiscriminantAnalysis(solver='svd', store_covariance=False, priors=self.priors)
                #self.clf = LinearDiscriminantAnalysis(solver='eigen', store_covariance=True, priors=self.priors)
            elif self.classifier == 'qda':
                self.clf = QuadraticDiscriminantAnalysis(priors=self.priors)
            elif self.classifier == 'rf':
                # the gala parameters
                #self.clf = RandomForestClassifier(n_estimators=100, criterion='entropy', max_depth=20,
                #    bootstrap=False, random_state=None)
                #self.clf = RandomForestClassifier(n_estimators=5*nfeatures,n_jobs=self.nthreads,max_depth=10)
                self.clf = RandomForestClassifier(n_estimators=256,n_jobs=self.nthreads,max_depth=16)
            elif self.classifier == 'svm':
                self.clf = SVC(kernel='rbf',probability=True,cache_size=2000)
            elif self.classifier == 'nb':
                self.clf = GaussianNB()
            elif self.classifier == 'kn':
                self.clf = KNeighborsClassifier(n_neighbors=10,n_jobs=self.nthreads)
            elif self.classifier == 'dc':
                self.clf = DecisionTreeClassifier(max_depth=10)
            elif self.classifier == 'ada':
                self.clf = AdaBoostClassifier()
            elif self.classifier == 'lr':
                self.clf = LogisticRegression(penalty='l2',dual=False,solver='sag',n_jobs=self.nthreads)
            else:
                assert(False)   # i never try anything, i just do it

            # train to the normalized data and merge or no merge targets
            self.clf.fit(sdata, target)

            if self.classifierout and not self.iterative_mode:
                with open(self.classifierout, 'wb') as f: dill.dump({'classifier':self.clf}, f)

        if self.dpSupervoxelClassifier_verbose:
            print('\tdone in %.4f s' % (time.time() - t))

        # do the agglomeration to use as input to next iteration for iterative mode
        thr = -1
        if self.iterative_mode:
            cnt_targets = 0
            for chunk in range(self.nchunks):
                if self.iterative_frag[chunk] is None or not self.iterative_frag[chunk].isTraining: continue

                # get the feature data for current training chunk only
                cdata = sdata[cnt_targets:cnt_targets+ntargets[chunk],:]
                cnt_targets += ntargets[chunk]

                # merge based on current classifier
                self.iterative_frag[chunk].subgroups_out[-1] = '%.8f' \
                    % self.threshold_subgroups[self.iterative_mode_count]
                clf_predict, thr = self.get_merge_predict_thr(cdata)
                self.iterative_frag[chunk].agglomerate(clf_predict,
                    doWrite=self.iterate_save_mask[self.iterative_mode_count])

                # make the next training iteration load from the current agglomerated supervoxels
                self.iterative_frag[chunk].srcfile = self.iterative_frag[chunk].outfile
                self.iterative_frag[chunk].subgroups = self.iterative_frag[chunk].subgroups_out

        if self.doplots:
            return self.createPlots(target,sdata,self.clf,self.export_plots,
                name=self.classifier + '_train_' + '_'.join([str(x) for x in self.test_chunks]) + \
                '_iter_' + str(self.iterative_mode_count), thr=thr, plot_features=self.plot_features)

    def test(self):
        if self.dpSupervoxelClassifier_verbose: print('\nTEST')

        for chunk in range(self.nchunks):
            cchunk, chunk_list_index, chunk_range_index = self.get_chunk_inds(chunk)

            if chunk_list_index not in self.test_chunks: continue
            data,sdata,thr = self._test(chunk, cchunk, chunk_list_index, chunk_range_index)

        # xxx - this only plots the LAST test chunk (currently this is only used for leave-one-out cross-validations)
        if self.doplots:
            assert(len(self.test_chunks) == 1)  # exporting plots/test data only works for one test chunk
            return self.createPlots(data['target'],sdata,self.clf,self.export_plots,
                name=self.classifier + '_test_' + '_'.join([str(x) for x in self.test_chunks]) + \
                '_iter_' + str(self.iterative_mode_count), thr=thr, plot_features=self.plot_features)

    def _test(self, ichunk, cchunk, chunk_list_index, chunk_range_index):
        print('Exporting testing data for chunk %d,%d,%d' % tuple(cchunk.tolist()))
        offset = self.offset_list[chunk_list_index,:]; size = self.size_list[chunk_list_index,:]

        FRAG = None
        if len(self.test_chunks) == 1 and self.testin:
            if self.dpSupervoxelClassifier_verbose:
                print('Loading testing data')
            with open(self.testin, 'rb') as f: data = dill.load(f)
            FRAG = data['FRAG']; data = data['data']

        frag = self.iterative_frag[ichunk] if self.iterative_mode else None
        subgroups_out = list(self.label_subgroups_out)

        if frag is None:
            if self.doplots:
                frag = dpFRAG.makeBothFRAG(self.labelfile, cchunk, size, offset,
                    [self.probfile, self.probaugfile], [self.rawfile, self.rawaugfile],
                    self.raw_dataset, self.gtfile, self.outfile, self.label_subgroups, subgroups_out,
                    G=FRAG, progressBar=self.progress_bar, feature_set=self.feature_set, has_ECS=self.has_ECS,
                    chunk_subgroups=self.chunk_subgroups, neighbor_only=self.neighbor_only,
                    pad_prob_svox_perim = not self.prob_svox_context, no_agglo_ECS=self.no_agglo_ECS,
                    verbose=self.dpSupervoxelClassifier_verbose)
            else:
                frag = dpFRAG.makeTestingFRAG(self.labelfile, cchunk, size, offset,
                    [self.probfile, self.probaugfile], [self.rawfile, self.rawaugfile],
                    self.raw_dataset, self.outfile, self.label_subgroups, subgroups_out, G=FRAG,
                    progressBar=self.progress_bar, feature_set=self.feature_set, has_ECS=self.has_ECS,
                    chunk_subgroups=self.chunk_subgroups, neighbor_only=self.neighbor_only,
                    pad_prob_svox_perim = not self.prob_svox_context, no_agglo_ECS=self.no_agglo_ECS,
                    verbose=self.dpSupervoxelClassifier_verbose)

        if self.iterative_mode and self.iterative_frag[ichunk] is None:
            frag.isTraining = False; self.iterative_frag[ichunk] = frag

            # write out the starting labels as as subgroup 0
            frag.subgroups_out[-1] = '%.8f' % 0.0
            frag.data_cube = frag.supervoxels_noperim
            verbose = frag.dpWriteh5_verbose; frag.dpWriteh5_verbose = frag.dpFRAG_verbose;
            #self.data_attrs['types_nlabels'] = [ncomps]
            frag.writeCube(); frag.dpWriteh5_verbose = verbose

        if not (len(self.test_chunks) == 1 and self.testin):
            frag.createFRAG(update = self.iterative_mode)
            #frag.createFRAG(update = False)
            data = frag.createDataset(train=self.doplots)

            if self.testout:
                if self.dpSupervoxelClassifier_verbose:
                    print('Dumping testing data')
                descr = 'Testing data from dpFRAG.py with command line:\n' + self.arg_str
                descr = ('With ini file "%s":\n' % (self.cfgfile,)) + self.ini_str
                data['DESCR'] = descr
                with open(self.testout, 'wb') as f: dill.dump({'data':data,'FRAG':frag.FRAG}, f)

        sdata = scale(data['data'])     # normalize for the classifiers

        thr = -1
        if self.iterative_mode:
            # merge based on current classifier
            frag.subgroups_out[-1] = '%.8f' % self.threshold_subgroups[self.iterative_mode_count]
            clf_predict, thr = self.get_merge_predict_thr(sdata)
            frag.agglomerate(clf_predict, doWrite=self.iterate_save_mask[self.iterative_mode_count])

            # make the next training iteration load from the current agglomerated supervoxels
            self.iterative_frag[ichunk].srcfile = self.iterative_frag[ichunk].outfile
            self.iterative_frag[ichunk].subgroups = self.iterative_frag[ichunk].subgroups_out
        else:
            try:
                # predict merge or not on testing cube and write outputs at specified probability thresholds
                frag.threshold_agglomerate(self.clf.predict_proba(sdata), self.thresholds, self.threshold_subgroups)
                # there's an issue here if there are no mergers left, would be better to just copy the current
                #   agglomeration, xxx - deal with this later. handled this explicitly in other locations.
                #except AttributeError:
            except:
                # if the classifier doesn't do probabilities just export single prediction
                frag.subgroups_out[-1] = ('single_' + self.classifier)
                frag.agglomerate(self.clf.predict(sdata))

        return data,sdata,thr

    def iterative_classify(self):
        assert(self.iterative_mode)

        if self.classifierin and os.path.isfile(self.classifierin):
            # this is a legacy mode where all classifiers are stored in same dill
            if self.dpSupervoxelClassifier_verbose:
                print('\nLoading iterative classifiers:'); t = time.time()

            with open(self.classifierin, 'rb') as f: d = dill.load(f)
            self.clfs = d['classifiers'];
            print('\tdone in %.4f s' % (time.time() - t))
        else:
            self.clfs = [None] * self.iterate_count

        # save the ROC-style metrics for later analysis / plotting
        train_metrics = [None] * self.iterate_count; test_metrics = [None] * self.iterate_count

        # each iteration loop trains on the current supervoxels and
        #   then performs a merge (test) based with a normal sklearn predict based on a small merge prior.
        self.trainin = ''; self.trainout = ''; self.testin = ''; self.testout = ''
        if self.test_only:
            print('Test-only iterative mode with merge prior %.4f' % (self.merge_prior,))
            for chunk in range(self.nchunks):
                cchunk, chunk_list_index, chunk_range_index = self.get_chunk_inds(chunk)
                #if chunk_list_index not in self.test_chunks: continue

                for i in range(self.iterate_count):
                    self.iterate_merge_perc = self.iterate_merge_perc_list[i]; self.iterative_mode_count = i
                    print('Iteration %d, merge perc %.4f' % (self.iterative_mode_count+1,self.iterate_merge_perc))

                    if self.clfs[i] is None and self.classifierin and not os.path.isfile(self.classifierin):
                        with open(os.path.join(self.classifierin,self.classifierin_name+'_'+\
                            str(self.iterative_mode_count)+'.dill'), 'rb') as f: d = dill.load(f)
                        self.clf = d['classifier']
                    else:
                        self.clf = self.clfs[i]

                    assert( self.clf is not None) # classifierin must be specified for test only
                    self._test(chunk, cchunk, chunk_list_index, chunk_range_index)

                # let the previous chunk get garbage collected to save memory.
                # that is essentially the purpose of this special mode (iteration loop as inner loop).
                self.iterative_frag[chunk] = None
        else:
            print('Normal iterative mode with merge prior %.4f' % (self.merge_prior,))
            for i in range(self.iterate_count):
                self.iterate_merge_perc = self.iterate_merge_perc_list[i]; self.iterative_mode_count = i
                print('Iteration %d, merge perc %.4f' % (self.iterative_mode_count+1,self.iterate_merge_perc))

                if self.clfs[i] is None and self.classifierin and not os.path.isfile(self.classifierin):
                    with open(os.path.join(self.classifierin,self.classifierin_name+'_'+\
                        str(self.iterative_mode_count)+'.dill'), 'rb') as f: d = dill.load(f)
                    self.clf = d['classifier']
                else:
                    self.clf = self.clfs[i]

                if self.clf is None:
                    train_metrics[i] = self.train()
                test_metrics[i] = self.test()

                if self.classifierout:
                    # new mode that stores all classifiers in separate dill files
                    fn = os.path.join(self.classifierout_path,self.classifierout_name+'_'+\
                        str(self.iterative_mode_count)+self.classifierout_ext)
                    with open(fn, 'wb') as f:
                        dill.dump({'classifier':self.clf,'train_metrics':train_metrics,'test_metrics':test_metrics},f)

    def get_merge_predict_thr(self,data):
        thr = -1; #ntargets = data.shape[0]
        # without merge percentage specified just use the classifier predict().
        clf_predict = self.clf.predict(data)
        # use the merge percentage to merge only this percentage of the most confident mergers.
        # do not use any potential mergers that are against the classifier predict().
        if self.iterate_merge_perc > 0:
            clf_probs = self.clf.predict_proba(data)
            # xxx - figure out how to handle no mergers left more gracefully, just stop at this iteration?
            assert(clf_probs.ndim > 1 and clf_probs.shape[1] > 1) # ran out of mergers
            # use merge percentage as percentage of estimated yes merges remaining
            nmerge = (clf_predict==1).sum(dtype=np.int64); clf_probs = clf_probs[:,1]
            if nmerge > 0:
                thr = np.sort(clf_probs)[::-1][int(nmerge*self.iterate_merge_perc)]
                clf_predict = np.logical_and(clf_predict, clf_probs >= thr)
        return clf_predict, thr

    def get_chunk_inds(self, chunk):
        if self.use_chunk_range:
            chunk_list_index = np.nonzero(chunk >= self.chunk_range_cumsize)[0][-1]
            chunk_range_index = chunk - self.chunk_range_cumsize[chunk_list_index]
            cchunk = np.unravel_index(chunk_range_index, self.chunk_range_rng[chunk_list_index,:]) \
                *self.chunkscale[chunk_list_index,:] + self.chunk_range_beg[chunk_list_index,:]
        else:
            chunk_list_index = chunk; cchunk = self.chunk_range_beg[chunk,:]; chunk_range_index = None
        return cchunk, chunk_list_index, chunk_range_index

    # this method assumes binary classification
    # xxx - plotting code got really messy as per usual, clean this up
    def createPlots(self,target,sdata,clf,export_path,plot_features=False,name='clf',figno=100,thr=-1):
        ntargets = target.size; nfeatures = sdata.shape[1]
        bins = [np.arange(-5,5.1,0.1)] * nfeatures

        pl.figure(figno); plt.clf();
        if plot_features:
            self.plotFeatures(target,sdata,export_path=None,show_plot=False,bins=bins,figno=figno)
            nx, ny = 100, 100
            for x in range(nfeatures):
                for y in range(x+1,nfeatures):
                    pl.subplot(nfeatures-1,nfeatures-1,x*(nfeatures-1)+y)

                    # class 0 and 1 : get boundaries between classes, depending on classifier
                    x_min, x_max = plt.gca().get_xlim()
                    y_min, y_max = plt.gca().get_ylim()
                    xx, yy = np.meshgrid(np.linspace(x_min, x_max, nx),
                                         np.linspace(y_min, y_max, ny))
                    Y = np.zeros((nx*ny,nfeatures)); Y[:,[y,x]] = np.c_[xx.ravel(), yy.ravel()]
                    ZZ = clf.predict(Y)
                    try:
                        # use predict_proba with target merge percentage method to plot boundaries (if available)
                        Z = clf.predict_proba(Y)
                        if thr > 0:
                            Z = np.logical_and(ZZ,(Z[:,1] > thr)).astype(np.double).reshape(xx.shape)
                        else:
                            Z = Z[:, 1].reshape(xx.shape)
                    except:
                        # otherwise just use regular predict to plot boundaries
                        Z = ZZ.astype(np.double).reshape(xx.shape)

                    if nfeatures==2:
                        # if there are only two features, overlay probabilities from above with the feature densities
                        img = 1-np.abs(Z-0.5);
                        pl.imshow(img,interpolation='nearest',extent=(x_min,x_max,y_min,y_max),
                            aspect=(y_max-y_min)/(x_max-x_min), origin='lower', alpha=0.3, cmap='gray',)

                    # overlay the boundaries between the features
                    plt.contour(xx, yy, Z, [0.5], linewidths=2., colors='w')

                    try:
                        # plot the class means if available in the classifier (LDA)
                        plt.scatter(clf.means_[0][y], clf.means_[0][x], s=2, color='r', edgecolors='w')
                        plt.scatter(clf.means_[1][y], clf.means_[1][x], s=2, color='g', edgecolors='w')
                    except AttributeError:
                        pass

        # calculate ROC / PR style metrics using classifier predict (with target merge percentage method) and print
        clf_probs = None; clf_preds = clf.predict(sdata)
        if thr > 0:
            clf_probs = clf.predict_proba(sdata)
            # avoid error if there are no predictions left that predict yes merge
            if clf_probs.ndim > 1 and clf_probs.shape[1] > 1:
                clf_preds = np.logical_and(clf_preds,(clf_probs[:,1] > thr))
        yesmerge = (target==1); notmerge = (target==0);
        nyes = yesmerge.sum(dtype=np.int64); nnot = notmerge.sum(dtype=np.int64)
        fScore, tpr_recall, precision, pixel_error, fpr, tp, tn, fp, fn = pixel_error_fscore( target.astype(bool),
            clf_preds.astype(bool) )
        print('p=%d, n=%d, tp=%d, tn=%d, fp=%d, fn=%d, rec=%.4f, prec=%.4f, fscore=%.4f' % (nyes,nnot,tp,tn,fp,fn,
            tpr_recall,precision,fScore))

        # some plots that only work with certain classifiers, skip for non-relevant classifiers

        pl.figure(figno+1); plt.clf()
        try:
            # plot if the classifier implements predict_proba
            axes = pl.subplot(1,2,1)
            pbins = np.arange(0,1.01,0.01); binw = (pbins[1]-pbins[0])/2; cbins = pbins[:-1]+binw
            if clf_probs is None: clf_probs = clf.predict_proba(sdata)
            tnhist,tmp = np.histogram(1-clf_probs[np.logical_and(notmerge,clf_preds==0),0],bins=pbins)
            tnhist = tnhist.astype(np.double)/nnot
            fnhist,tmp = np.histogram(1-clf_probs[np.logical_and(yesmerge,clf_preds==0),0],bins=pbins)
            fnhist = fnhist.astype(np.double)/nyes
            # avoid error if there are no predictions left that predict yes merge
            if clf_probs.ndim > 1 and clf_probs.shape[1] > 1:
                tphist,tmp = np.histogram(clf_probs[np.logical_and(yesmerge,clf_preds==1),1],bins=pbins)
                tphist = tphist.astype(np.double)/nyes
                fphist,tmp = np.histogram(clf_probs[np.logical_and(notmerge,clf_preds==1),1],bins=pbins)
                fphist = fphist.astype(np.double)/nnot
            else:
                tphist = np.zeros_like(tnhist); fphist = np.zeros_like(fnhist)
            pl.plot(cbins,tphist,'g'); pl.plot(cbins,fphist,'g--')
            pl.plot(cbins,tnhist,'r'); pl.plot(cbins,fnhist,'r--')
            plt.xlim([-0.05,1.05]); plt.ylim([-0.05,1.05])
            plt.xlabel('P(merge)'); plt.ylabel('density'); plt.title('tp/fn Y(%d), tn/fp N(%d)' % (nyes,nnot))

            # plot for LDA classifier only
            axes = pl.subplot(1,2,2)
            plt.xticks(range(nfeatures),rotation=45)
            #FEATURES, FEATURES_NAMES, nfeatures = dpFRAG.make_FEATURES()    # xxx - meh
            axes.set_xticklabels(self.features_names)
            pl.plot(clf.scalings_[:,0],'k')
            d = clf.decision_function(clf.means_)   # distance from means to the decision boundary
            pl.title('dist: not %.5f yes %.5f' % (d[0],d[1]))
        except (AttributeError, TypeError) as e:
            pass

        if export_path:
            figna = [x % (name,) for x in ['merge_features_%s.png', 'merge_probs_eigen_%s.png']]
            nfigna = len(figna)
            for f,i in zip(range(figno, figno+nfigna), range(nfigna)):
                pl.figure(f)
                figure = plt.gcf() # get current figure
                figure.set_size_inches(20, 20)
                plt.savefig(os.path.join(export_path,figna[i]), dpi=72)
        else:
            pl.show()

        # return the ROC style metrics
        metrics = {}
        for i in ('nyes', 'nnot', 'tp', 'tn', 'fp', 'fn', 'tpr_recall', 'precision', 'fScore', 'pixel_error'):
            metrics[i] = locals()[i]
        return metrics

    def plotFeatures(self,target,fdata,export_path,show_plot=True,bins=None,figno=100):
        nfeatures = fdata.shape[1]; #ntargets = target.size
        yesmerge = (target==1); notmerge = (target==0);
        nyes = yesmerge.sum(dtype=np.int64); nnot = notmerge.sum(dtype=np.int64)
        if not bins:
            bins = [np.arange(0,4.6,0.125), np.arange(0,6.15,0.15), np.arange(0.75,3.1,0.075), np.arange(0,257,8),
                np.arange(0,1.025,0.025), np.arange(0,1.025,0.025), np.arange(0,1.025,0.025)]
        nbins = [x.size-1 for x in bins]; print(nbins)
        binw = [(x[1]-x[0])/2 for x in bins]; cbins = [x[:-1]+y for x,y in zip(bins,binw)]
        pl.figure(figno)
        if show_plot: plt.clf()
        for x in range(nfeatures):
            for y in range(x+1,nfeatures):
                pl.subplot(nfeatures-1,nfeatures-1,x*(nfeatures-1)+y)

                img = np.zeros((nbins[x],nbins[y],3),dtype=np.double)
                img[:,:,1] = np.histogram2d(fdata[yesmerge, x], fdata[yesmerge, y], bins=(bins[x],bins[y]))[0]/nyes
                img[:,:,0] = np.histogram2d(fdata[notmerge, x], fdata[notmerge, y], bins=(bins[x],bins[y]))[0]/nnot
                sel = (img > 0); img[sel] = -np.log10(img[sel]); sel = (img > 0); img[sel] = 1-img[sel]/img.max()
                #sel = (img > 0.05);
                bnd = nd.measurements.find_objects(sel);
                #xlim = (bins[x][bnd[0][0].start]+binw[x]/10, bins[x][bnd[0][0].stop]-binw[x]/10)
                #ylim = (bins[y][bnd[0][1].start]+binw[y]/10, bins[y][bnd[0][1].stop]-binw[y]/10)
                if len(bnd) == 0:
                    xlim = [0,1]; ylim = [0,1]
                else:
                    xlim = (cbins[x][bnd[0][0].start], cbins[x][bnd[0][0].stop-1])
                    ylim = (cbins[y][bnd[0][1].start], cbins[y][bnd[0][1].stop-1])

                # imshow uses f-order so x/y are flipped
                pl.imshow(img,interpolation='nearest',extent=(cbins[y][0],cbins[y][-1],cbins[x][0],
                    cbins[x][-1]), aspect=(cbins[y][-1]-cbins[y][0])/(cbins[x][-1]-cbins[x][0]), origin='lower')
                plt.xlim(ylim); plt.ylim(xlim)

                if y==x+1:
                    #FEATURES, FEATURES_NAMES, nfeatures = dpFRAG.make_FEATURES()    # xxx - meh
                    plt.xlabel(self.features_names[y])
                    plt.ylabel(self.features_names[x])
                else:
                    plt.gca().axes.get_xaxis().set_visible(False)
                    plt.gca().axes.get_yaxis().set_visible(False)

        if show_plot:
            if export_path:
                plt.gcf().set_size_inches(20, 20)
                plt.savefig(os.path.join(export_path,'merge_features.png'), dpi=72)
            else:
                pl.show()

    @staticmethod
    def get_options(cfgfile):
        config = ConfigObj(cfgfile,
            configspec=os.path.join(os.path.dirname(os.path.realpath(__file__)),'dpSupervoxelClassifier.ini'))

        # Validator handles missing / type / range checking
        validator = Validator()
        results = config.validate(validator, preserve_errors=True)
        if results != True:
            for (section_list, key, err) in flatten_errors(config, results):
                if key is not None:
                    if not err:
                        print('The "%s" key is missing in the following section(s):%s ' \
                            % (key, ', '.join(section_list)))
                        raise ValidateError
                    else:
                        print('The "%s" key in the section(s) "%s" failed validation' \
                            % (key, ', '.join(section_list)))
                        raise err
                elif section_list:
                    print('The following section(s) was missing:%s ' % ', '.join(section_list))
                    raise ValidateError

        return config

    @staticmethod
    def addArgs(p):
        p.add_argument('--cfgfile', nargs=1, type=str, default='', help='Path/name of ini config file')
        p.add_argument('--trainin', nargs=1, type=str, default='', help='Input file for loading training data (dill)')
        p.add_argument('--classifier', nargs=1, type=str, default='lda', help='Which sklearn classifier to use')
        p.add_argument('--classifierin', nargs=1, type=str, default='',
            help='Input file for loading trained classifier(s) (dill)')
        p.add_argument('--classifierout', nargs=1, type=str, default='',
            help='Output file for saving trained classifier(s) (dill)')
        p.add_argument('--testin', nargs=1, type=str, default='', help='Input file for loading testing data (dill)')
        p.add_argument('--test-chunks', nargs='*', type=int, default=[],
            metavar='CHUNKS', help='Chunks to use for test (override from .ini)')
        p.add_argument('--show-plots', action='store_true', help='Show various plots')
        p.add_argument('--export-plots', nargs=1, type=str, default='',
            help='Export various plots to this path (default no plots)')
        p.add_argument('--plot-features', action='store_true', help='If plotting, whether to include feature plots')
        p.add_argument('--outfile', nargs=1, type=str, default='', help='Override output file for agglomerations')
        p.add_argument('--labelfile', nargs=1, type=str, default='', help='Override input label (supervoxel) file')
        p.add_argument('--probfile', nargs=1, type=str, default='', help='Override input prob file')
        p.add_argument('--probaugfile', nargs=1, type=str, default='', help='Override input prob augment file')
        p.add_argument('--progress-bar', action='store_true', help='Enable progress bar if available')
        p.add_argument('--feature-set', nargs=1, type=str, default='',
            help='Option to control which FRAG features are calculated (override from .ini)')
        # these are so standard cubeIter inputs can be used from command line to override from .ini
        p.add_argument('--chunk', nargs=3, type=int, default=[-1,-1,-1], metavar=('X', 'Y', 'Z'),
            help='Corner chunk to parse out of hdf5')
        p.add_argument('--offset', nargs=3, type=int, default=[0,0,0], metavar=('X', 'Y', 'Z'),
            help='Offset in chunk to read')
        p.add_argument('--size', nargs=3, type=int, default=[256,256,128], metavar=('X', 'Y', 'Z'),
            help='Size in voxels to read')
        p.add_argument('--nthreads', nargs=1, type=int, default=[8],
            help='Number of parallel threads to set for scipy and scikit-learn')
        p.add_argument('--chunk-subgroups', action='store_true',
            help='This mode is for probs and labels that have overlapping context so need to be stored separately.')
        p.add_argument('--neighbor-only', dest='neighbor_only', action='store_true',
            help='Only use boundary voxels labeled with neighboring supervoxels (no background / non-neighbor voxels)')
        p.add_argument('--prob-svox-context', dest='prob_svox_context', action='store_true',
            help='Use context for loading probs and supervoxels along edge faces')
        p.add_argument('--no-agglo-ECS', action='store_true', help='Do not agglomerate ECS supervoxels')
        p.add_argument('--useFRAGc', action='store_true', help='Use the C-optimized version of FRAG')

        p.add_argument('--dpSupervoxelClassifier-verbose', action='store_true',
            help='Debugging output for dpSupervoxelClassifier')

if __name__ == '__main__':
    parser = argparse.ArgumentParser(description='Flattened Supervoxel Classifier',
        formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    dpSupervoxelClassifier.addArgs(parser)
    args = parser.parse_args()

    svoxClass = dpSupervoxelClassifier(args)
    if svoxClass.iterative_mode:
        svoxClass.iterative_classify()
    else:
        svoxClass.train()
        svoxClass.test()

max_stars_repo_name: erjel/emdrp

max_issues_repo_name: erjel/emdrp

max_forks_repo_name: erjel/emdrp

reason: import numpy,from scipy


text:
"""
In this file we define initializers, the classes
that can be used to initialize various parameters of
networks, such as weights, biases and hidden states.
"""

# Standard libraries
from abc import ABC, abstractmethod

# Third-party libraries
import numpy as np

class Initializer(ABC):
    """
    Abstract base class for initializers.
    """
    @staticmethod
    @abstractmethod
    def initialize(shape):
        """
        This function returns a randomly initialized
        parameter with the given shape.
        :param shape: Shape of the tensor.
        :return: Randomly initialized tensor of shape
        'shape'.
        """

class StandardGaussian(Initializer):
    """
    This class implements standard gaussian
    initialization.
    """
    @staticmethod
    def initialize(shape):
        return np.random.randn(shape[0], shape[1])

class Xavier(Initializer):
    """
    This class implements the Xavier
    initialization method.
    """
    @staticmethod
    def initialize(shape):
        return np.random.randn(shape[0], shape[1]) / np.sqrt((shape[0]))

class XavierReLU(Initializer):
    """
    This class implements the Xavier
    initialization method for ReLUs.
    """
    @staticmethod
    def initialize(shape):
        return np.random.randn(shape[0], shape[1]) / np.sqrt((shape[0] / 2))

class NullInitializer(Initializer):
    """
    This class implements null initialization.
    """
    @staticmethod
    def initialize(shape):
        return np.zeros(shape, dtype="float64")

ext: py

max_stars_repo_name: BendeguzToth/NeuralLanguageModel

max_issues_repo_name: BendeguzToth/NeuralLanguageModel

max_forks_repo_name: BendeguzToth/NeuralLanguageModel

reason: import numpy


text:
from .compare_results import compare
from . import stats
from . import timeclass

def init(base_directory,obs_directory=None,timeperiods=None):
	global BASE_DIRECTORY,OBSERVED_DIRECTORY,TIME_PERIODS
	BASE_DIRECTORY = base_directory
	OBSERVED_DIRECTORY = obs_directory
	TIME_PERIODS = timeperiods

def _globFilenamesOnly(directory,search):
    from glob import glob
    import os
    return [fn.split(os.path.sep)[-1] for fn in glob(os.path.join(directory,search))]

def available():
    return _globFilenamesOnly(BASE_DIRECTORY,'*')

def value_diff(dfL,dfR):
    import numpy as np
    oCols = list(dfL.columns[np.where(dfL.dtypes=='object')])
    numericCols = list(set(dfL.columns)-set(oCols))
    delta = dfL.copy()
    for col in numericCols:
        delta[col] = dfL[col] - dfR[col]
    return delta

def Veneer(*args,**kwargs):
    from veneer import Veneer
    result = Veneer(*args,**kwargs)

    def configureOptions(self,options):
        self.configure_options(options)

    def configure_options(self,options):
        lines = ["# Generated Script","from Dynamic_SedNet.PluginSetup import DSScenarioDetails"]
        lines += ["DSScenarioDetails.%s = %s"%(k,v) for (k,v) in options.items()]
        script = '\n'.join(lines)
        #print(script)
        res = self.model._safe_run(script)

    def set_run_name(self,name):
        self.model.set('scenario.CurrentConfiguration.runName',name,literal=True)

    def set_output_path(self,name):
        self.model.set('scenario.CurrentConfiguration.ModelOutputPath',name.replace('\\','\\\\'),literal=True)

    def run_contributor(self,results,fn):
        import os
        script = """
import Dynamic_SedNet.Results.Contributor.RegionalReportingContributor as RegionalReportingContributor
import Dynamic_SedNet.Tools.ToolsGeneral as ToolsGeneral
contrib = RegionalReportingContributor()
contrib.UseResults('%s')
contrib.Scenario = scenario
contrib.runTimeStep()
if contrib.Success:
  result = None
  ToolsGeneral.ExportTableToCsvString(contrib.contribsTab,True,'%s')
else:
  result = contrib.Error
"""%(results.path,os.path.join(results.path,'%s.csv'%fn))
        self.model.runScript(script)

    result.configureOptions = configureOptions.__get__(result,result.__class__)
    result.configure_options = configure_options.__get__(result,result.__class__)
    result.set_run_name = set_run_name.__get__(result,result.__class__)
    result.set_output_path = set_output_path.__get__(result,result.__class__)
    result.run_contributor = run_contributor.__get__(result,result.__class__)
    return result

class Results(object):
    def __init__(self,run_name):
        from .queries import ResultsQueries
        import os
        self._time_periods = None
        self.run_name = run_name
        self.path = os.path.join(BASE_DIRECTORY,run_name)
        self.runDetails = RunDetails(self.path)
        self.queries = ResultsQueries(self)
        self._loaded = {}

    def _apply_time_series_helpers(self,dataframe):
        from pandas import DataFrame
#        import types
        def of_timeclass(df_self,timeclass):
            time_periods = self.time_periods().ix[df_self.index]
            if isinstance(timeclass,str):
                result = df_self[time_periods['cat']==timeclass]
            else:
                result = df_self[time_periods['cat'].apply(timeclass)]
            self._apply_time_series_helpers(result)
            return result

        def by_timeclass(df_self):
            time_periods = self.time_periods().ix[df_self.index]
            result = df_self.groupby(time_periods['cat'])
            #self._apply_time_series_helpers(result)
            return result

        def by_wateryear(df_self):
            time_periods = self.time_periods().ix[df_self.index]
            result = df_self.groupby(time_periods.Water_year)
            #self._apply_time_series_helpers(result)
            return result

        def of_month(df_self,month):
            result = df_self[df_self.index.month==month]
            self._apply_time_series_helpers(result)
            return result

        def by_month(df_self):
            result = df_self.groupby(df_self.index.month)
#            self._apply_time_series_helpers(result)
            return result

        def of_year(df_self,year):
            result = df_self[df_self.index.year==year]
            self._apply_time_series_helpers(result)
            return result

        def by_year(df_self):
            result = df_self.groupby(df_self.index.year)
#            self._apply_time_series_helpers(result)
            return result

        meths = {'of_timeclass':of_timeclass,
            'by_timeclass':by_timeclass,
            'by_wateryear':by_wateryear,
            'of_month':of_month,
            'by_month':by_month,
            'of_year':of_year,
            'by_year':by_year,
        }

        dataframe.__class__ = type('JoelsDataFrame',(DataFrame,),meths)
#        dataframe.__setattr__('of_timeclass',types.MethodType(of_timeclass,dataframe))
#        dataframe.by_timeclass = types.MethodType(by_timeclass,dataframe)
#        dataframe.by_wateryear = types.MethodType(by_wateryear,dataframe)
#        dataframe.of_month = types.MethodType(of_month,dataframe)
#        dataframe.by_month = types.MethodType(by_month,dataframe)
#        dataframe.of_year = types.MethodType(of_year,dataframe)
#        dataframe.by_year = types.MethodType(by_year,dataframe)

#        dataframe.__orig_dir__ = dataframe.__dir__
#        def new_dir(df_self):
#            return df_self.__orig_dir__() + ['of_timeclass','of_month','of_year',
#                'by_year', 'by_month', 'by_wateryear', 'by_timeclass']

#        dataframe.__setattr__('__dir__',types.MethodType(new_dir,dataframe))

    def available(self):
        return [fn[0:-4] for fn in _globFilenamesOnly(self.path,'*.csv')]

    def get(self,results_type):
        if not results_type in self._loaded:
            import pandas as pd
            import os
            self._loaded[results_type] = pd.DataFrame.from_csv(os.path.join(self.path,results_type+'.csv'))
        return self._loaded[results_type]

    def get_node_ts(self,site,constituents,header=True):
        import pandas as pd
        import os
        if isinstance(constituents,str):
            constituents = [constituents]

        def _fn(con):
            return os.path.join(self.path,'TimeSeries',con,"%s_%s_kilograms.csv"%(con,site))
        def _read(con):
            header_rows = 1 if header else 0
            return pd.read_csv(_fn(con),index_col=0,parse_dates=True,dayfirst=True,names = ['Date', con],skiprows=1)
        result = pd.DataFrame([_read(c)[c] for c in constituents]).transpose()
        self._apply_time_series_helpers(result)
        return result

    def load_observed_ts(self,site,fn=None,location=None):
        """
        Load a time series file of observed loads for a given site.
        """
        import pandas as pd
        if fn is None:
            if location is None:
                location = OBSERVED_DIRECTORY
            import os
            from glob import glob
            matching_files = glob(os.path.join(location,"*%s*.csv"%site))
            fn = matching_files[0]
        result = pd.read_csv(fn,index_col=0,parse_dates=True,dayfirst=True)
        self._apply_time_series_helpers(result)
        return result

    def time_periods(self,fn=None):
        if self._time_periods is None:
            if fn is None:
                fn = TIME_PERIODS
            import pandas as pd
            self._time_periods = pd.read_csv(fn,index_col = 1, parse_dates = [1], dayfirst = True, names=["Water_year","date","cat"], skiprows =1)
        return self._time_periods

class DifferenceResults(object):
    def __init__(self,run1,run2):
        from .queries import ResultsQueries
        self._run1 = Results(run1)
        self._run2 = Results(run2)
        self.run_name = '('+self._run1.run_name + ' - ' + self._run2.run_name+')'
        self.path = None
        self.runDetails = self._run1.runDetails
        self.queries = ResultsQueries(self)

    def available():
        return self._run1.available()

    def get(self,results_type):
        return value_diff(self._run1.get(results_type),self._run2.get(results_type))

    def get_node_ts(self,site,constituents):
        return value_diff(self._run1.get_node_ts(site,constituents),self._run2.get(site,constituents))

    def load_observed_ts(self,site,fn=None,location=None):
        return self._run1.load_observed_ts(site,fn,location)

class RunDetails(object):
    def __init__(self,path):
        import os
        import xml.etree.ElementTree
        self._data = xml.etree.ElementTree.parse(os.path.join(path,'DSScenarioRunInfo.xml')).getroot()
        self.start = self._extractDate('startDate')
        self.startRecording = self._extractDate('startRecDate')
        self.end = self._extractDate('endDate')
        simulationLength = self.end - self.startRecording
        self.yearsOfRecording = simulationLength.days / 365.25

    def _extractDate(self,element):
        from datetime import datetime
        return datetime.strptime(self._data.find(element).text,'%Y-%M-%dT00:00:00')

class MyTest(object):
    def __dir__(self):
        return ['Hi']

max_stars_repo_name: flowmatters/dsed-py

max_issues_repo_name: flowmatters/dsed-py

max_forks_repo_name: flowmatters/dsed-py


text:
import numpy as np
import copy

from frozenlake import FrozenLakeEnv
from plot_utils import plot_values

def policy_improvement(env, V, gamma=1):
    policy = np.zeros([env.nS, env.nA]) / env.nA
    for s in range(env.nS):
        q = q_from_v(env, V, s, gamma)

        # OPTION 1: construct a deterministic policy
        # policy[s][np.argmax(q)] = 1

        # OPTION 2: construct a stochastic policy that puts equal probability on maximizing actions
        best_a = np.argwhere(q == np.max(q)).flatten()
        policy[s] = np.sum([np.eye(env.nA)[i] for i in best_a], axis=0) / len(best_a)

    return policy

def q_from_v(env, V, s, gamma=1):
    q = np.zeros(env.nA)
    for a in range(env.nA):
        for prob, next_state, reward, done in env.P[s][a]:
            q[a] += prob * (reward + gamma * V[next_state])
    return q

def truncated_policy_evaluation(env, policy, V, max_it=1, gamma=1):
    num_it=0
    while num_it < max_it:
        for s in range(env.nS):
            v = 0
            q = q_from_v(env, V, s, gamma)
            for a, action_prob in enumerate(policy[s]):
                v += action_prob * q[a]
            V[s] = v
        num_it += 1
    return V

def truncated_policy_iteration(env, max_it=1, gamma=1, theta=1e-8):
    V = np.zeros(env.nS)
    policy = np.zeros([env.nS, env.nA]) / env.nA
    while True:
        policy = policy_improvement(env, V)
        old_V = copy.copy(V)
        V = truncated_policy_evaluation(env, policy, V, max_it, gamma)
        if max(abs(V-old_V)) < theta:
            break;
    return policy, V

env = FrozenLakeEnv()

policy_tpi, V_tpi = truncated_policy_iteration(env, max_it=2)

# print the optimal policy
print("\nOptimal Policy (LEFT = 0, DOWN = 1, RIGHT = 2, UP = 3):")
print(policy_tpi,"\n")

# plot the optimal state-value function
plot_values(V_tpi)

max_stars_repo_name: RaoulMa/NeuralNets

max_issues_repo_name: RaoulMa/BasicNeuralNets

max_forks_repo_name: RaoulMa/BasicNeuralNets


text:
import numpy as np

import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.projections import PolarAxes
from matplotlib.transforms import Affine2D, Transform
from matplotlib.testing.decorators import image_comparison

from mpl_toolkits.axes_grid1.parasite_axes import ParasiteAxes
from mpl_toolkits.axisartist import SubplotHost
from mpl_toolkits.axes_grid1.parasite_axes import host_subplot_class_factory
from mpl_toolkits.axisartist import angle_helper
from mpl_toolkits.axisartist.axislines import Axes
from mpl_toolkits.axisartist.grid_helper_curvelinear import \
    GridHelperCurveLinear

@image_comparison(['custom_transform.png'], style='default', tol=0.2)
def test_custom_transform():
    class MyTransform(Transform):
        input_dims = output_dims = 2

        def __init__(self, resolution):
            """
            Resolution is the number of steps to interpolate between each input
            line segment to approximate its path in transformed space.
            """
            Transform.__init__(self)
            self._resolution = resolution

        def transform(self, ll):
            x, y = ll.T
            return np.column_stack([x, y - x])

        transform_non_affine = transform

        def transform_path(self, path):
            ipath = path.interpolated(self._resolution)
            return Path(self.transform(ipath.vertices), ipath.codes)

        transform_path_non_affine = transform_path

        def inverted(self):
            return MyTransformInv(self._resolution)

    class MyTransformInv(Transform):
        input_dims = output_dims = 2

        def __init__(self, resolution):
            Transform.__init__(self)
            self._resolution = resolution

        def transform(self, ll):
            x, y = ll.T
            return np.column_stack([x, y + x])

        def inverted(self):
            return MyTransform(self._resolution)

    fig = plt.figure()

    SubplotHost = host_subplot_class_factory(Axes)

    tr = MyTransform(1)
    grid_helper = GridHelperCurveLinear(tr)
    ax1 = SubplotHost(fig, 1, 1, 1, grid_helper=grid_helper)
    fig.add_subplot(ax1)

    ax2 = ParasiteAxes(ax1, tr, viewlim_mode="equal")
    ax1.parasites.append(ax2)
    ax2.plot([3, 6], [5.0, 10.])

    ax1.set_aspect(1.)
    ax1.set_xlim(0, 10)
    ax1.set_ylim(0, 10)

    ax1.grid(True)

@image_comparison(['polar_box.png'], style='default', tol=0.04)
def test_polar_box():
    # Remove this line when this test image is regenerated.
    plt.rcParams['text.kerning_factor'] = 6

    fig = plt.figure(figsize=(5, 5))

    # PolarAxes.PolarTransform takes radian. However, we want our coordinate
    # system in degree
    tr = Affine2D().scale(np.pi / 180., 1.) + PolarAxes.PolarTransform()

    # polar projection, which involves cycle, and also has limits in
    # its coordinates, needs a special method to find the extremes
    # (min, max of the coordinate within the view).
    extreme_finder = angle_helper.ExtremeFinderCycle(20, 20,
                                                     lon_cycle=360,
                                                     lat_cycle=None,
                                                     lon_minmax=None,
                                                     lat_minmax=(0, np.inf))

    grid_locator1 = angle_helper.LocatorDMS(12)
    tick_formatter1 = angle_helper.FormatterDMS()

    grid_helper = GridHelperCurveLinear(tr,
                                        extreme_finder=extreme_finder,
                                        grid_locator1=grid_locator1,
                                        tick_formatter1=tick_formatter1)

    ax1 = SubplotHost(fig, 1, 1, 1, grid_helper=grid_helper)

    ax1.axis["right"].major_ticklabels.set_visible(True)
    ax1.axis["top"].major_ticklabels.set_visible(True)

    # let right axis shows ticklabels for 1st coordinate (angle)
    ax1.axis["right"].get_helper().nth_coord_ticks = 0
    # let bottom axis shows ticklabels for 2nd coordinate (radius)
    ax1.axis["bottom"].get_helper().nth_coord_ticks = 1

    fig.add_subplot(ax1)

    ax1.axis["lat"] = axis = grid_helper.new_floating_axis(0, 45, axes=ax1)
    axis.label.set_text("Test")
    axis.label.set_visible(True)
    axis.get_helper().set_extremes(2, 12)

    ax1.axis["lon"] = axis = grid_helper.new_floating_axis(1, 6, axes=ax1)
    axis.label.set_text("Test 2")
    axis.get_helper().set_extremes(-180, 90)

    # A parasite axes with given transform
    ax2 = ParasiteAxes(ax1, tr, viewlim_mode="equal")
    assert ax2.transData == tr + ax1.transData
    # Anything you draw in ax2 will match the ticks and grids of ax1.
    ax1.parasites.append(ax2)
    ax2.plot(np.linspace(0, 30, 50), np.linspace(10, 10, 50))

    ax1.set_aspect(1.)
    ax1.set_xlim(-5, 12)
    ax1.set_ylim(-5, 10)

    ax1.grid(True)

@image_comparison(['axis_direction.png'], style='default', tol=0.07)
def test_axis_direction():
    # Remove this line when this test image is regenerated.
    plt.rcParams['text.kerning_factor'] = 6

    fig = plt.figure(figsize=(5, 5))

    # PolarAxes.PolarTransform takes radian. However, we want our coordinate
    # system in degree
    tr = Affine2D().scale(np.pi / 180., 1.) + PolarAxes.PolarTransform()

    # polar projection, which involves cycle, and also has limits in
    # its coordinates, needs a special method to find the extremes
    # (min, max of the coordinate within the view).

    # 20, 20 : number of sampling points along x, y direction
    extreme_finder = angle_helper.ExtremeFinderCycle(20, 20,
                                                     lon_cycle=360,
                                                     lat_cycle=None,
                                                     lon_minmax=None,
                                                     lat_minmax=(0, np.inf),
                                                     )

    grid_locator1 = angle_helper.LocatorDMS(12)
    tick_formatter1 = angle_helper.FormatterDMS()

    grid_helper = GridHelperCurveLinear(tr,
                                        extreme_finder=extreme_finder,
                                        grid_locator1=grid_locator1,
                                        tick_formatter1=tick_formatter1)

    ax1 = SubplotHost(fig, 1, 1, 1, grid_helper=grid_helper)

    for axis in ax1.axis.values():
        axis.set_visible(False)

    fig.add_subplot(ax1)

    ax1.axis["lat1"] = axis = grid_helper.new_floating_axis(
        0, 130,
        axes=ax1, axis_direction="left")
    axis.label.set_text("Test")
    axis.label.set_visible(True)
    axis.get_helper().set_extremes(0.001, 10)

    ax1.axis["lat2"] = axis = grid_helper.new_floating_axis(
        0, 50,
        axes=ax1, axis_direction="right")
    axis.label.set_text("Test")
    axis.label.set_visible(True)
    axis.get_helper().set_extremes(0.001, 10)

    ax1.axis["lon"] = axis = grid_helper.new_floating_axis(
        1, 10,
        axes=ax1, axis_direction="bottom")
    axis.label.set_text("Test 2")
    axis.get_helper().set_extremes(50, 130)
    axis.major_ticklabels.set_axis_direction("top")
    axis.label.set_axis_direction("top")

    grid_helper.grid_finder.grid_locator1.set_params(nbins=5)
    grid_helper.grid_finder.grid_locator2.set_params(nbins=5)

    ax1.set_aspect(1.)
    ax1.set_xlim(-8, 8)
    ax1.set_ylim(-4, 12)

    ax1.grid(True)


text:
# TensorFlow Input Pipelines for Large Data Sets
# ischlag.github.io
# TensorFlow 0.11, 07.11.2016

import tensorflow as tf
import numpy as np
import threading

# Generating some simple data
r = np.arange(0.0,100003.0)
raw_data = np.dstack((r,r,r,r))[0]
raw_target = np.array([[1,0,0]] * 100003)

# are used to feed data into our queue
queue_input_data = tf.placeholder(tf.float32, shape=[20, 4])
queue_input_target = tf.placeholder(tf.float32, shape=[20, 3])

queue = tf.FIFOQueue(capacity=50, dtypes=[tf.float32, tf.float32], shapes=[[4], [3]])

enqueue_op = queue.enqueue_many([queue_input_data, queue_input_target])
dequeue_op = queue.dequeue()

# tensorflow recommendation:
# capacity = min_after_dequeue + (num_threads + a small safety margin) * batch_size
data_batch, target_batch = tf.train.batch(dequeue_op, batch_size=15, capacity=40)
# use this to shuffle batches:
# data_batch, target_batch = tf.train.shuffle_batch(dequeue_op, batch_size=15, capacity=40, min_after_dequeue=5)

def enqueue(sess):
  """ Iterates over our data puts small junks into our queue."""
  under = 0
  max = len(raw_data)
  while True:
    print("starting to write into queue")
    upper = under + 20
    print("try to enqueue ", under, " to ", upper)
    if upper <= max:
      curr_data = raw_data[under:upper]
      curr_target = raw_target[under:upper]
      under = upper
    else:
      rest = upper - max
      curr_data = np.concatenate((raw_data[under:max], raw_data[0:rest]))
      curr_target = np.concatenate((raw_target[under:max], raw_target[0:rest]))
      under = rest

    sess.run(enqueue_op, feed_dict={queue_input_data: curr_data,
                                    queue_input_target: curr_target})
    print("added to the queue")
  print("finished enqueueing")

# start the threads for our FIFOQueue and batch
sess = tf.Session()
enqueue_thread = threading.Thread(target=enqueue, args=[sess])
enqueue_thread.isDaemon()
enqueue_thread.start()

coord = tf.train.Coordinator()
threads = tf.train.start_queue_runners(coord=coord, sess=sess)

# Fetch the data from the pipeline and put it where it belongs (into your model)
for i in range(5):
  run_options = tf.RunOptions(timeout_in_ms=4000)
  curr_data_batch, curr_target_batch = sess.run([data_batch, target_batch], options=run_options)
  print(curr_data_batch)

# shutdown everything to avoid zombies
sess.run(queue.close(cancel_pending_enqueues=True))
coord.request_stop()
coord.join(threads)
sess.close()
Blog

max_stars_repo_name: kafkasl/contextualLSTM

max_issues_repo_name: kafkasl/contextualLSTM

max_forks_repo_name: kafkasl/contextualLSTM


text:
# -*- coding: utf-8 -*-
import sys
sys.path.append('/notebooks')

import glob
import numpy as np
import os.path

input_files = glob.glob('/data/input/*encoded.npy')

in_data = []
output = []
for filename in input_files:
    data = np.load(filename)
    seq = len(data)
    for i in range(seq - 45):
        in_data.append(data[i:i+45])

in_data = np.array(in_data)
np.random.shuffle(in_data)

np.save('/data/input/composed_random_input', in_data)


text:
# -*- coding: utf-8 -*-
"""input_image.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1n4mOXiDrCxQL0abSwwIbBP5YXaFn_Axj
"""

import torch
from torchvision import transforms
from google.colab.patches import cv2_imshow
import cv2
import numpy as np
from random import *
from imgaug import augmenters as iaa
import math
# Commented out IPython magic to ensure Python compatibility.
# %cd /content/drive/My Drive/image_cutout/

####Randomly cutout paricular portion of the bounding box
path='이미지 상대경로 및 절대경로'

cutout=cv2.imread(path)
label=open('/content/drive/My Drive/label_cutout/__95bDVirnJgqPmjbHjQAA.txt','r')
list=label.readline()
while(list):
    list=list.split()
    x_center=float(list[1])
    y_center=float(list[2])
    height=float(list[4])*832
    width=float(list[3])*1664
    x_center_new=x_center*1664
    y_center_new=y_center*832
    print(x_center_new,y_center_new)
    #radomly set x1 which is greater than leftmost x coordinates and less than center+1/4 of width
    x1=uniform(x_center_new-width/2,x_center_new+width/4)
    #radomly set y1 which is greater than leftmost y coordinates and less than center+1/4 of height
    y1=uniform(y_center_new-height/2,y_center_new+height/4)
    mark=width/2 #2 is empirical number
    x2=x1+mark
    y2=y1+mark
    x1=int(x1)
    x2=int(x2)
    y1=int(y1)
    y2=int(y2)
    print(x1," ",y1," ",x2," ",y2)
    cutoout = cv2.rectangle(cutout, (x1, y1), (x2,y2), (0, 0, 0), cv2.FILLED)
    list=label.readline()

cv2_imshow(cutout)


text:
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate
import tecplot as tp
import wake_config as conf

import matplotlib.ticker as ticker

width = 3
#from helpers.setup_plot import *
#plot_style = 'THESIS'
#line_style, markers = setup_plot('THESIS', width, width*0.8)

def lumley_map():
    #rechts: axisymmetry
    Jx1=np.linspace(0,2./9.,300)
    Jy1=(3./2.) * (Jx1 *(4./3.))**(2./3.)
    #% links: axisymmetry
    Jx2=np.linspace(0,-1./36.,50)
    Jy2=3./2. * (-Jx2  *4.0/3.0)**(2./3.)
    #% oben: 2-component-turbulence
    Jx3=np.linspace(-1./36.,2./9.,300)
    Jy3=2./9. + Jx3*2
    return Jx1, Jy1, np.flipud(Jx2), np.flipud(Jy2), Jx3, Jy3

def draw_map(x, z, xlin, zlin, invar2, invar3):
#xi, zi = np.linspace(x0, x0, num), np.linspace(z0, z1, num)
    I2 = scipy.interpolate.griddata((x, z), invar2, (xlin, zlin), method='linear')
    I3 = scipy.interpolate.griddata((x, z), invar3, (xlin, zlin), method='linear')

    print('line x is from ' + str(xlin[0]) + ' to ' + str(xlin[-1]))

    Jx1, Jy1, Jx2, Jy2, Jx3, Jy3  = lumley_map()
    width = 3
    line_style, markers = setup_plot('THESIS', width, width*0.8)

    fig, ax = plt.subplots(1,1)
    plt.plot(Jx1, Jy1, Jx2, Jy2, Jx3, Jy3, color='k')
    pcm = plt.scatter(I3, I2, s=2, c = xlin, cmap = plt.cm.viridis)
    #cbar = plt.colorbar(pcm, ticks=[1.05, 1.65])
    #cbar.set_label('$(x-x_{TE})/c_{local}$', labelpad=-5)
    #cbar.ax.set_yticklabels(['0', '3'])
    adjustprops = dict(left=0.12, bottom=0.12, right=0.97, top=0.97, wspace=0.2, hspace=0.2)       # Subplot properties
    adjustprops = dict(left=0.14, bottom=0.2, right=0.97, top=0.97, wspace=0.2, hspace=0.2)       # Subplot properties
    plt.subplots_adjust(**adjustprops)
    plt.xlabel(r'$III_a$', labelpad = 5)
    plt.ylabel(r'$II_a$', labelpad=-9)

    img_name = 'wake_centerline_anisotropy_lumley.pdf'
    plt.savefig(img_name)
    print('written ' + img_name)

    plt.close()

def draw_barycentric(C, xb, yb, x, z, xlin, zlin):
    zlin = np.squeeze(zlin)

    do_annotate = False
    xbi = scipy.interpolate.griddata((x, z), xb, (xlin, zlin), method='linear')
    ybi = scipy.interpolate.griddata((x, z), yb, (xlin, zlin), method='linear')
    print(x.shape)
    print(xlin.shape)
    print(zlin.shape)
    print(xbi.shape)

    width = 2
    fig, ax = plt.subplots(1,1)
#connectpoints()
    x1, y1 = [0, 0.5], [0, np.sqrt(3.0)/2.0]
    x2, y2 = [0, 1], [0, 0]
    x3, y3 = [1, 0.5], [0, np.sqrt(3.0)/2.0]
    plt.plot(x1, y1, x2, y2, x3, y3, color='k')
#plt.xlim([0,1])
#plt.ylim([0,1])
    ax.get_yaxis().set_visible(False)
    ax.get_xaxis().set_visible(False)
    plt.axis('off')
    adjustprops = dict(left=0.05, bottom=0.05, right=0.97, top=0.97, wspace=0.2, hspace=0.2)       # Subplot properties
    plt.subplots_adjust(**adjustprops)

    plt.scatter(xbi, ybi, s=2, c = xlin, cmap=plt.cm.viridis)

    if do_annotate:
        ax.annotate('upstream', xy=(xbi[0], ybi[0]), xytext=(-20,-20),
                            textcoords='offset points', ha='center', va='top',
                            bbox=dict(boxstyle='round,pad=0.2', fc='white', alpha=1),
                            arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=0.5',
                                                                        color='black'), fontsize=11)

        ax.annotate('downstream', xy=(xbi[-1], ybi[-1]), xytext=(-40,40),
                            textcoords='offset points', ha='center', va='top',
                            bbox=dict(boxstyle='round,pad=0.2', fc='white', alpha=1),
                            arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=-0.5',
                                                                        color='black'), fontsize=11)

#ax.text(0, 0, r'$x_{2c}$', color='black' , fontsize=12, ha='right', va='top')
#ax.text(x2[1], y2[1], r'$x_{1c}$', color='black' , fontsize=12, ha='left', va='top')
#ax.text(x1[1], y1[1]+0.03, r'$x_{3c}$', color='black' , fontsize=12, ha='center', va='bottom')

    ax.text(0, 0, r'$2c$', color='black' , fontsize=12, ha='right', va='top')
    ax.text(x2[1], y2[1], r'${1c}$', color='black' , fontsize=12, ha='left', va='top')
    ax.text(x1[1], y1[1]+0.03, r'${3c}$', color='black' , fontsize=12, ha='center', va='bottom')

# Rotate angle
    angle = 45
#trans_angle = plt.gca().transData.transform_angles(np.array((45,)), l2.reshape((1, 2)))[0]
# Plot text
    th1 = plt.text(0.2, np.sqrt(3)/4.0+0.05, 'axisymmetric contraction', fontsize=11,  rotation=60, rotation_mode='anchor', ha='center', va='center')
    th2 = plt.text(0.8, np.sqrt(3)/4.0+0.05, 'axisymmetric expansion', fontsize=11,  rotation=-60, rotation_mode='anchor', ha='center', va='center')

    adjustprops = dict(left=0.05, bottom=0.07, right=0.95, top=0.93, wspace=0.2, hspace=0.2)       # Subplot properties
    plt.subplots_adjust(**adjustprops)
    plt.gca().set_aspect('equal')

    img_name = 'wake_centerline_anisotropy_barycentric.pdf'
    plt.savefig(img_name)
    print('written ' + img_name)

    plt.close()

max_stars_repo_name: awaldm/flowpost

max_issues_repo_name: awaldm/flowpost

max_forks_repo_name: awaldm/flowpost

reason: import numpy,import scipy


text:
#!/usr/bin/env python
import helper as h
import numpy as np
import pandas as pd
import csv, scipy, json
import warnings
from tqdm import tqdm

warnings.filterwarnings("ignore")

np.set_printoptions(precision=3, suppress=True)

# Parameters
repetitions = 10
datasets = h.datasets()
clfs = h.classifiers().keys()
results = h.results()
p_s = h.p_s()

# Iterating results
for p in p_s:
    print(p)
    for result in tqdm(results, ascii=True):
        # Get params for results
        dbname = result[2]
        repetition = int(result[3][1:])
        cv_method = result[4]
        filename = result[0]
        df = pd.read_csv(filename)

        # Showoff
        # print "Repetition %i with %s on %s (p=%.2f)" % (repetition, cv_method, dbname,1-p)

        # Calculate mean scores and std_s
        mean_scores = df.mean(axis=0).values
        std_scores = df.std(axis=0).values

        # Classifier classification
        adv_w = [0] * len(clfs)
        adv_t = [0] * len(clfs)

        # Going through classifiers
        for i, clf_a in enumerate(clfs):
            clf_acc = mean_scores[i]
            # print "\n---\n%s - ACC = %.3f" % (clf_a, clf_acc)
            leading = np.delete(mean_scores < clf_acc, i)
            # print leading

            acc_vector_leader = df[clf_a].values

            p_ws, p_ts = [], []
            for j, clf_b in enumerate(clfs):
                if j == i:
                    continue
                acc_vector_b = df[clf_b].values

                _, p_w = scipy.stats.wilcoxon(acc_vector_leader, acc_vector_b)
                _, p_t = scipy.stats.ttest_ind(acc_vector_leader, acc_vector_b)
                p_ws.append(p_w)
                p_ts.append(p_t)
            p_ws = np.array(p_ws)
            p_ts = np.array(p_ts)

            mask_p_ws = p_ws < p
            mask_p_ts = p_ts < p
            actual_p_ws = mask_p_ws * leading
            actual_p_ts = mask_p_ts * leading
            a_ws = sum(actual_p_ws)
            a_ts = sum(actual_p_ts)
            """
            print "| Wilcoxon:"
            print p_ws
            print mask_p_ws
            print actual_p_ws
            print a_ws
            print "\n| T-student"
            print p_ts
            print mask_p_ts
            print actual_p_ts
            print a_ts
            """

            adv_w[i] = int(a_ws)
            adv_t[i] = int(a_ts)

        # Collect and write analysis
        data = {
            "parameters": {
                "db": dbname,
                "repetition": repetition,
                "cv_method": cv_method,
            },
            "clfs": list(clfs),
            "mean": list(mean_scores),
            "std": list(std_scores),
            "adv_w": adv_w,
            "adv_t": adv_t,
        }
        jsonname = "jsons/%s_r%i_%s_p%i.json" % (
            dbname,
            repetition,
            cv_method,
            int(p * 100),
        )
        # print jsonname
        with open(jsonname, "w") as outfile:
            json.dump(data, outfile, indent=2, sort_keys=True)


text:
# Copyright 2017 Neural Networks and Deep Learning lab, MIPT
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

from typing import List, Iterable, Union

import numpy as np

from deeppavlov.core.models.nn_model import NNModel

class SiameseModel(NNModel):
    """The class implementing base functionality for siamese neural networks.

    Args:
        batch_size: A size of a batch.
        num_context_turns: A number of ``context`` turns in data samples.
        *args: Other parameters.
        **kwargs: Other parameters.
    """

    def __init__(self,
                 batch_size: int,
                 num_context_turns: int = 1,
                 *args, **kwargs) -> None:
        super().__init__(*args, **kwargs)

        self.batch_size = batch_size
        self.num_context_turns = num_context_turns

    def load(self, *args, **kwargs) -> None:
        pass

    def save(self, *args, **kwargs) -> None:
        pass

    def train_on_batch(self, batch: List[List[np.ndarray]], y: List[int]) -> float:
        b = self._make_batch(list(batch))
        loss = self._train_on_batch(b, y)
        return loss

    def __call__(self, batch: Iterable[List[np.ndarray]]) -> Union[np.ndarray, List[str]]:
        y_pred = []
        buf = []
        for j, el in enumerate(batch, start=1):
            context = el[:self.num_context_turns]
            responses = el[self.num_context_turns:]
            buf += [context + [el] for el in responses]
            if len(buf) >= self.batch_size:
                for i in range(len(buf) // self.batch_size):
                    b = self._make_batch(buf[i*self.batch_size:(i+1)*self.batch_size])
                    yp = self._predict_on_batch(b)
                    y_pred += list(yp)
                lenb = len(buf) % self.batch_size
                if lenb != 0:
                    buf = buf[-lenb:]
                else:
                    buf = []
        if len(buf) != 0:
            b = self._make_batch(buf)
            yp = self._predict_on_batch(b)
            y_pred += list(yp)
        y_pred = np.asarray(y_pred)
        if len(responses) > 1:
            y_pred = np.reshape(y_pred, (j, len(responses)))
        return y_pred

    def reset(self) -> None:
        pass

    def _train_on_batch(self, batch: List[np.ndarray], y: List[int]) -> float:
        pass

    def _predict_on_batch(self, batch: List[np.ndarray]) -> np.ndarray:
        pass

    def _predict_context_on_batch(self, batch: List[np.ndarray]) -> np.ndarray:
        pass

    def _predict_response_on_batch(self, batch: List[np.ndarray]) -> np.ndarray:
        pass

    def _make_batch(self, x: List[List[np.ndarray]]) -> List[np.ndarray]:
        b = []
        for i in range(len(x[0])):
            z = [el[i] for el in x]
            b.append(np.asarray(z))
        return b

max_stars_repo_name: ineersa/DeepPavlov

max_issues_repo_name: ineersa/DeepPavlov

max_forks_repo_name: ineersa/DeepPavlov


text:
#!/usr/bin/python3

import unittest
from unittest.mock import patch

import numpy as np
from numpy.linalg.linalg import LinAlgError

from tmc import points

from tmc.utils import load, get_out, patch_helper

module_name = 'src.meeting_planes'
meeting_planes = load(module_name, 'meeting_planes')
ph = patch_helper(module_name)

@points('p03-07.1')
class MeetingPlanes(unittest.TestCase):

    def test_first(self):
        a1=1
        b1=4
        c1=5

        a2=3
        b2=2
        c2=1

        a3=2
        b3=4
        c3=1
        p=(a1,b1,c1,
           a2,b2,c2,
           a3,b3,c3)
        system="(a1=%i, b1=%i, c1=%i, a2=%i, b2=%i, c2=%i, a3=%i, b3=%i, c3=%i)" % p

        x, y, z = meeting_planes(a1, b1, c1,
                                 a2, b2, c2,
                                 a3, b3, c3)

        self.assertAlmostEqual(z, a1*y + b1*x+c1, msg="Solution %f does not satisfy the first equation of system %s!" % (z, system))
        self.assertAlmostEqual(z, a2*y + b2*x+c2, msg="Solution %f does not satisfy the second equation of system %s!" % (z, system))
        self.assertAlmostEqual(z, a3*y + b3*x+c3, msg="Solution %f does not satisfy the third equation of system %s!" % (z, system))

    def test_calls(self):
        a1=1
        b1=4
        c1=5

        a2=3
        b2=2
        c2=1

        a3=2
        b3=4
        c3=1

        with patch(ph("np.linalg.solve"), wraps=np.linalg.solve) as psolve:
            meeting_planes(a1, b1, c1, a2, b2, c2, a3, b3, c3)
            psolve.assert_called_once()

    def test_underdetermined(self):
        a1=1
        b1=4
        c1=2
        p=(a1,b1,c1,
           a1,b1,c1,
           a1,b1,c1)
        system="(a1=%i, b1=%i, c1=%i, a2=%i, b2=%i, c2=%i, a3=%i, b3=%i, c3=%i)" % p
        with self.assertRaises(np.linalg.linalg.LinAlgError,
                               msg="Under determined system %s should raise an exception!" % system):
            meeting_planes(*p)

    def test_inconsistent(self):
        a1=1
        b1=4
        c1=2
        p=(a1,b1,c1,
           a1,b1,c1,
           a1,b1,c1+1)
        system="(a1=%i, b1=%i, c1=%i, a2=%i, b2=%i, c2=%i, a3=%i, b3=%i, c3=%i)" % p
        with self.assertRaises(np.linalg.linalg.LinAlgError,
                               msg="Inconsistent system %s should raise an exception!" % system):
            meeting_planes(*p)

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

reason: import numpy,from numpy


text:
# DEFINICION DE FUNCIONES
# Por: Sebastian Sarasti
# Lista de funciones para el tratamiento de datos de un detector de NaI
# TASA DE CUENTAS
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def tasa(spectranuclei):
    """This function helps to transform a spectra in gross spectra. Time of exposition is the
    second value of data recolected"""
    tim = spectranuclei[1]
    for i in spectranuclei:
        spectra_g = spectranuclei[2:len(spectranuclei)] / tim
    return spectra_g

# CONTEO NETO
def netcount(spectranuclei, spectrabk):
    """
    This function help to estimate net count rate.
    The time is not necessary because the second value of vector
    is time of sampling in seconds. Because of zero-index, Python uses
    the first value
    """
    timnuclei = spectranuclei[1]
    timbk = spectrabk[1]
    spectranuclei = spectranuclei[2:len(spectranuclei)]
    spectrabk = spectrabk[2:len(spectrabk)]
    spectran_g = spectranuclei / timnuclei
    spectrabk_g = spectrabk / timbk
    spectran_n = spectran_g - spectrabk_g
    for i in range(len(spectran_n)):
        if spectran_n[i] < 0:
            spectran_n[i] = 0
        else:
            spectran_n[i] = spectran_n[i]
    return spectran_n

# SUAVIZADO MOVIL 3 VALORES
def smooth_three(spectra):
    spectra_smoothed = np.zeros(len(spectra))
    for i in range(3, len(spectra) - 1):
        spectra_smoothed[i + 1] = (spectra[i] + spectra[i - 1] + spectra[i - 2]) / 3
    return spectra_smoothed

# SUAVIZADO MOVIL 5 VALORES
def smooth_five(spectra):
    spectra_smoothed = np.zeros(len(spectra))
    for i in range(4, len(spectra) - 1):
        spectra_smoothed[i + 1] = (spectra[i] + spectra[i - 1] + spectra[i - 2] - spectra[i - 3] + spectra[i - 4]) / 5
    return spectra_smoothed

# NORMALIZAR HISTOGRAMA SIMULADO
def normalizar_sim(histograma):
    h = np.histogram(histograma, bins=150)
    y = h[0]
    x = h[1]
    x = x[0:len(x) - 1]
    s = max(y)
    y_n = np.zeros(len(y))
    for i in range(len(y_n)):
        y_n[i] = y[i] / s
    return (x, y_n)

# NORMALIZAR ESPECTROS
def normalizar_exp(spectra, nuclei):
    if nuclei == "Cs":
        m = max(spectra[3000:len(spectra)])
        spectra_nor = np.zeros(len(spectra))
        for i in range(len(spectra)):
            spectra_nor[i] = spectra[i] / m
    else:
        m = max(spectra)
        spectra_nor = np.zeros(len(spectra))
        for i in range(len(spectra)):
            spectra_nor[i] = spectra[i] / m
    return (spectra_nor)

# DISTRIBUCION GAUSSIANO
def spectrum(E, osc, sigma):
    """
    Añade la resolucion del detector al espectro simulado
    """
    x = np.linspace(0, 2, 1500, endpoint=True)
    gE = []
    for Ei in x:
        tot = 0
        for Ej, os in zip(E, osc):
            tot += os * np.exp(-((((Ej - Ei) / sigma) ** 2)))
        gE.append(tot)
    x = x[0:1500]
    gE = gE[0:1500]
    return x, gE

# CALCULO DEL SIGMA EN FUNCION DE LA RESOLUCION
def gauss_fit(spec_e_neto,spec_s_bruto,energia, nuclei):
    if nuclei=="Cs":
        spec_e_neto=list(np.float_(spec_e_neto))
        spec_s_bruto=np.histogram(spec_s_bruto,bins=250)
        x=np.linspace(0,spec_s_bruto[1][250], num=1000, endpoint=True)
        sigma=0.001
        res_dif=0
        m=max(spec_e_neto[3000:len(spec_e_neto)])
        indice=spec_e_neto.index(m)
        # se determina la mitad de la altura y se buscan los valores superiores a este
        n=m/2
        p=np.where(spec_e_neto[0:indice]<=n)
        j1=p[0]
        m1=j1[-1]
        q=np.where(spec_e_neto[indice:len(spec_e_neto)]<=n)
        j2=q[0]
        j2=j2+indice
        m2=j2[0]
        # una vez identificados los valores, se busca los indices de este en el arreglo de energia
        x1=energia[m1]
        x2=energia[m2]
        xm=energia[indice]
        # calculo de la resolucion
        res_e=(x2-x1)/xm
        while res_dif<0.01:
            # se suaviza el espectro simulado
            gE=spectrum(spec_s_bruto[1][0:250],spec_s_bruto[0],sigma,x)
            # se determina el maximo y su indice
            m=max(gE)
            indice=gE.index(max(gE))
            # se determina la mitad de la altura y se buscan los valores superiores a este
            n=m/2
            p=np.where(gE[0:indice]<=n)
            j1=p[0]
            m1=j1[-1]
            q=np.where(gE[indice:-1]<=n)
            j2=q[0]
            j2=j2+indice
            m2=j2[0]
            # una vez identificados los valores, se busca los indices de este en el arreglo de energia
            x1=x[m1]
            x2=x[m2]
            xm=x[indice]
            # calculo de la resolucion
            res_s=(x2-x1)/xm
            sigma=sigma+0.001
            res_dif=(res_s-res_e)/res_e
    else:
        spec_e_neto=list(np.float_(spec_e_neto))
        spec_s_bruto=np.histogram(spec_s_bruto,250)
        x=np.linspace(0,spec_s_bruto[1][250], num=1000, endpoint=True)
        sigma=0.001
        res_dif=0
        m=max(spec_e_neto)
        indice=spec_e_neto.index(max(spec_e_neto))
        # se determina la mitad de la altura y se buscan los valores superiores a este
        n=m/2
        p=np.where(spec_e_neto[0:indice]<=n)
        j1=p[0]
        m1=j1[-1]
        q=np.where(spec_e_neto[indice:len(spec_e_neto)]<=n)
        j2=q[0]
        j2=j2+indice
        m2=j2[0]
        # una vez identificados los valores, se busca los indices de este en el arreglo de energia
        x1=energia[m1]
        x2=energia[m2]
        xm=energia[indice]
        # calculo de la resolucion
        res_e=(x2-x1)/xm
        while res_dif<0.01:
            # se suaviza el espectro simulado
            gE=spectrum(spec_s_bruto[1][0:250],spec_s_bruto[0],sigma,x)
            # se determina el maximo y su indice
            m=max(gE)
            indice=gE.index(max(gE))
            # se determina la mitad de la altura y se buscan los valores superiores a este
            n=m/2
            p=np.where(gE[0:indice]<=n)
            j1=p[0]
            m1=j1[-1]
            q=np.where(gE[indice:-1]<=n)
            j2=q[0]
            j2=j2+indice
            m2=j2[0]
            # una vez identificados los valores, se busca los indices de este en el arreglo de energia
            x1=x[m1]
            x2=x[m2]
            xm=x[indice]
            # calculo de la resolucion
            res_s=(x2-x1)/xm
            sigma=sigma+0.001
            res_dif=(res_s-res_e)/res_e
    return (gE,x,sigma,res_s,res_e,res_dif)

# AREA
def area(spectra, energy, a, b):
    """
    This function helps to determinate area under each photopeak through,
    trapeziodal method.
    Where:
        spectra: net count rate vector
        energy: energy vector
        a: lower limit
        b: upper limit
    """
    spec = []
    spec = spectra
    spec = spec[a:b]
    energy = energy[a:b]
    a = 0
    a = np.trapz(spec, energy)
    return a

# INCERTIBUMBRE
def incertidumbre(spectra, bk, energy, a, b, tiempo_s, tiempo_bk):
    a_s = 0
    a_bk = 0
    a_s = area(spectra, energy, a, b)
    a_b = area(bk, energy, a, b)
    incer = (a_s / tiempo_s + a_b / tiempo_bk) ** 0.5
    return incer

def max_peak(spectra, a, b):
    """
    This function helps to determinate area under each photopeak through,
    trapeziodal method.
    Where:
        spectra: net count rate vector
        energy: energy vector
        a: lower limit
        b: upper limit
    """
    spec = []
    spec = spectra
    spec = spec[a:b]
    return max(spec)

def spec_hist(his, ran):
    a = np.histogram(his, bins=ran)
    x = a[1][0:8190]
    y = a[0]
    return (x, y)

def spec_s_act(spec, pa_s, ac):
    a = pa_s / ac
    spec_c = spec[1] / a
    return (spec[0], spec_c)

def spec_s_f(spec, f):
    x = spec[0]
    y = spec[1] / f
    return (x, y)