text:
#!/usr/bin/python3

# Functions to handle Input
#############################################################################################
def read_csv(file):
  # simple function to read data from a file
  data = pd.read_csv(file, sep=';')
  return data

# Functions to handle string/value conversion
#############################################################################################
# function converts format (DD-)HH:MM:SS to seconds
def ave2sec(x):
  if ( '-' in x ):
    vals = x.split('-')
    times = vals[1].split(':')
    sec = 24*3600*int(vals[0])+3600*int(times[0])+60*int(times[1])+int(times[2])
  else:
    times = x.split(':')
    sec = 3600*int(times[0])+60*int(times[1])+int(times[2])
  return (sec)

def scalingref(x):
  # returns reference scaling factor for MPI jobs based on 1.5 factor:
  # doubling cores should make parformance x1.5 (or better)
  if int(x) == 1:
    ref = 1
  else:
    ref = np.power(1/1.5,np.log2(int(x)))

  return ref

def rss2g(x):
  return int(float(x[:-1]))/1024

def reqmem2g(x):
  return int(float(x[:-2]))/1024

# Functions to handle DataFrames
#############################################################################################
def parse_df(data):
  # convert multi-line DataFrame to more compact form for analysis
  import datetime
  from dateutil import parser

  data[['id','subid']] = data.JobID.str.split('_',1, expand=True)
  data.drop(['subid'],axis=1, inplace=True)

  df=pd.DataFrame()
  df=data[~data['JobID'].str.contains("\.")]
  df.rename(columns={'State': 'Parentstate'}, inplace=True)

  data2=data.shift(-1).dropna(subset=['JobID'])
  df2=data2[data2['JobID'].str.contains("\.batch")]

  data2=data.shift(-2).dropna(subset=['JobID'])
  df3=data2[data2['JobID'].str.contains("\.0")]

  df.update(df2.MaxRSS)
  df.update(df3.MaxRSS)
  df.update(df2.AveCPU)
  df.update(df3.AveCPU)
  df=df.join(df2[['State']])
  df.update(df3.State)

  # drop columns that all all nan
  df.dropna(axis=1, inplace=True, how='all')

  # drop rows where any element is nan (errors in the data)
  df.dropna(axis=0, inplace=True, how='any')
  df.reset_index(inplace=True)
  df.loc[:,'State']=df.State.apply(str)

  # reduce data point, 5min accuracy instead of seconds just fine
  df['Submit']=pd.to_datetime(df['Submit']).dt.round('5min')
  df['Start']=pd.to_datetime(df['Start']).dt.round('5min')
  df['End']=pd.to_datetime(df['End']).dt.round('5min')

  # add extra columns to df for analysis purposes
  df.insert(len(df.columns),'runtime',(pd.to_datetime(df['End'])-pd.to_datetime(df['Start']))/np.timedelta64(1,'s'))
  df.insert(len(df.columns),'waittime',(pd.to_datetime(df['Start'])-pd.to_datetime(df['Submit']))/np.timedelta64(1,'s'))

  # NOTE: this value ok only if State=Completed, for State=timeout values close to RawCPUTime
  df.insert(len(df.columns),'cpuavesec',df['AveCPU'].apply(ave2sec) )

  # scaling reference number for MPI jobs
  df.insert(len(df.columns),'cpufactor',df['cpuavesec']/df['runtime'] )

  df.loc[:,'MaxRSS']=df['MaxRSS'].apply(rss2g)
  df.insert(len(df.columns),'MemReqTotal', df['ReqMem'].apply(reqmem2g)*df['NCPUS'] )
  df.insert(len(df.columns),'MemRatio', df['MaxRSS']/df['MemReqTotal'] )
  df.insert(len(df.columns),'TimeRatio', df['runtime']/df['Timelimit'].apply(ave2sec)  )
  df.drop(['AveCPU','CPUTimeRAW','JobID','Elapsed','Account','ReqMem'], axis=1, inplace=True)

  # pick up only gpu and gpushort partitions and release memory and jobs > 5min
  df2=df[(df['Partition'].str.contains("gpu")) & (df['runtime'] > 300) ]
  del df
  df2.to_csv('df2.csv', sep=';', index=True)

  return df2

def times_df(df):
  start = df.Submit.min().floor('1H')
  end = df.End.max().floor('1H')+np.timedelta64(1,'h')
  steps = int((end-start).total_seconds()/3600)
  df2 = pd.DataFrame(index=pd.date_range(start,end,freq='1H'), columns=['Running'])
  df2.Running = 0

  for i in range(0,steps):
    s1=start+i*np.timedelta64(1,'h')
    e1=s1+np.timedelta64(1,'h')
    dfr=df[((df['Start'] >= s1) & (df['Start'] <= e1)) | ((df['End'] >= s1) & (df['End'] <= e1)) ]
    dfq=df[((df['Submit'] >= s1) & (df['Submit'] <= e1)) | ((df['Start'] >= s1) & (df['Start'] <= e1)) ]
    df2.loc[s1,'Running']=int(len(dfr))
    df2.loc[s1,'Queue']=int(len(dfq))

  return df2

# Functions to handle Output summary
#############################################################################################
def print_data(data):
  print('============================')
  print(data)

def writecsv(data):
  data.to_csv('result.csv', sep=';', index=True)

if __name__ == '__main__':
    from argparse import ArgumentParser
    import sys
    import pandas as pd
    import numpy as np
    import time

    parser = ArgumentParser(description='Create histogram stats from sacct data')
    parser.add_argument("file", help="input csv-file")
    args = parser.parse_args()

    times=[]
    comment=[]

    times.append([time.process_time(),'start'])
    data=read_csv(args.file)
    times.append([time.process_time(),'read csv'])
    data=parse_df(data)
    times.append([time.process_time(),'parse_df'])
    dfh=times_df(data)
    times.append([time.process_time(),'times_df'])

    #print_data(dfh)
    writecsv(dfh)

    prev=times[0][0]
    #print(times)
    for i in times:
        print("%s : %.2f seconds" % (i[1], i[0]-prev) )
        prev=i[0]


text:
# ワーシャルフロイド
from scipy.sparse.csgraph import csgraph_from_dense, floyd_warshall

H, W = map(int, input().split())
c = [list(map(int, input().split())) for _ in range(10)]
c = floyd_warshall(csgraph_from_dense(c))

result = 0
for _ in range(H):
    for i in map(int, input().split()):
        if i == -1:
            continue
        result += c[i][1]
print(int(result))

ext: py

max_stars_repo_name: c-yan/atcoder

max_issues_repo_name: c-yan/atcoder

max_forks_repo_name: c-yan/atcoder

reason: from scipy


text:
"""
Copyright 2019 Wietse de Vries and Martijn Bartelds

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.
"""

import numpy as np
import json
from glob import glob
import os
import time
import argparse
import pickle
import codecs
import warnings

import pandas as pd
import scipy as sp
from conll_df import conll_df

from sklearn.base import BaseEstimator, TransformerMixin, ClassifierMixin
from sklearn.pipeline import Pipeline, FeatureUnion
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.preprocessing import MaxAbsScaler
from sklearn.decomposition import TruncatedSVD
from sklearn.calibration import CalibratedClassifierCV
from sklearn.svm import LinearSVC
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import f1_score, precision_recall_fscore_support
from sklearn.metrics.classification import UndefinedMetricWarning

warnings.filterwarnings('ignore', category=UndefinedMetricWarning)

FLOAT_DTYPES = (np.float64, np.float32, np.float16)

def read_problem_info(path, problem_name):
    problem_path = path + os.sep + problem_name
    probleminfo_path = problem_path + os.sep + 'problem-info.json'

    with open(probleminfo_path, 'r') as f:
        probleminfo = json.load(f)

    def list_files(candidate): return glob(problem_path + os.sep + candidate + os.sep + '*.txt')

    candidates = {c['author-name']: list_files(c['author-name']) for c in probleminfo['candidate-authors']}
    return candidates, list_files(probleminfo['unknown-folder'])

def read_collection_info(path):
    with open(path + os.sep + 'collection-info.json', 'r') as f:
        collectioninfo = json.load(f)

    for probleminfo in collectioninfo:
        candidates, unk_folder = read_problem_info(path, probleminfo['problem-name'])
        probleminfo['candidates'] = candidates
        probleminfo['unknown'] = unk_folder

    return collectioninfo

def save_output(outpath, problem_name, test_names, predictions):
    os.makedirs(outpath, exist_ok=True)
    out_data = []

    for unknown, prediction in zip(test_names, predictions):
        unknown_name = unknown.split('/')[-1]
        out_data.append(
            {'unknown-text': unknown_name, 'predicted-author': prediction})

    with open(outpath + os.sep + 'answers-' + problem_name + '.json', 'w') as f:
        json.dump(out_data, f, indent=4)

    print(' > answers saved to file', 'answers-' + problem_name + '.json')

def get_model_inputs(candidates):
    filenames, labels = [], []
    for c, c_files in candidates.items():
        filenames.extend(c_files)
        labels.extend([c] * len(c_files))

    return filenames, labels

def get_test_labels(filenames, collectioninfo, path):
    filenames = [f.replace(path + os.sep, '').split(os.sep) for f in filenames]
    labels = []
    truths = {}
    for problem, _, filename in filenames:
        if problem not in truths:
            with open(path + os.sep + problem + '/ground-truth.json', 'r') as f:
                truth = json.load(f)['ground_truth']
                truths[problem] = {t['unknown-text']: t['true-author'] for t in truth}
        labels.append(truths[problem][filename])

    return labels

class FeatureExtractor(BaseEstimator, TransformerMixin):
    """ Based on https://scikit-learn.org/0.18/auto_examples/hetero_feature_union.html """

    def fit(self, x, y=None):
        return self

    def _read_text(self, filename):
        """ Read raw texts from .txt file """
        with codecs.open(filename, 'r', encoding='utf-8') as f:
            return f.read()

    def _read_words(self, filename):
        """ Read tokens and pos tags from .txt.pos file """
        with codecs.open(filename + '.pos', 'r', encoding='utf-8') as f:
            tokens, pos_tags = [], []
            for line in f.readlines():
                word = line.rstrip().split()
                if len(word) == 2:
                    tokens.append(word[0])
                    pos_tags.append(word[1])
            return tokens, pos_tags

    def _read_syntax(self, filename):
        dep_df = conll_df(filename + '.dep', file_index=False)
        if dep_df.index.levels[1].dtype != np.int64:
            dep_df.reset_index(inplace=True)
            dep_df = dep_df[dep_df['i'].str.isnumeric()]
            dep_df['i'] = dep_df['i'].astype(int)
            dep_df.set_index(['s', 'i'], inplace=True)
        return dep_df.rename(columns={'w': 'word', 'l': 'lemma', 'f': 'dep'})

    def transform(self, filenames):
        features = np.recarray(shape=(len(filenames),),
                               dtype=[
                                   ('text', object),
                                   ('tokens', object),
                                   ('pos_tags', object),
                                   ('syntax', object)
        ])
        for i, filename in enumerate(filenames):
            features['text'][i] = self._read_text(filename)
            features['tokens'][i], features['pos_tags'][i] = self._read_words(filename)
            features['syntax'][i] = self._read_syntax(filename)
        return features

class FeatureSelector(BaseEstimator, TransformerMixin):
    def __init__(self, key):
        self.key = key

    def fit(self, x, y=None):
        return self

    def transform(self, data):
        return data[self.key]

class ProbabilityThreshold(BaseEstimator, ClassifierMixin):
    """ Sets uncertain predictions to <UNK> """

    def __init__(self, estimator, pt):
        self.estimator = estimator
        self.pt = pt

    def fit(self, X, y):
        self.classes_ = list(set(y))
        self.estimator.fit(X, y)

    def predict(self, X):
        proba = self.estimator.predict_proba(X)
        y = self.estimator.classes_[np.argmax(proba, axis=1)]

        proba.sort(axis=1)
        diff = proba[:, -1] - proba[:, -2]

        y[diff < self.pt] = '<UNK>'
        return y

def identity(x):
    return x

def _document_group_frequency(X, y=None):
    X = X.toarray()
    X[X > 0] = 1

    if y is None:
        return X.shape[0], X.sum(axis=0)

    group_df = pd.DataFrame(X).groupby(y).sum()

    # in_df = np.vstack([group_df.loc[label] for label in y])
    df = np.vstack([np.asarray(X.sum(axis=0)).flatten() - group_df.loc[label].values for label in y]) + 1

    n_samples = np.vstack([np.sum(y != label) for label in y]) + 1
    return n_samples, df

# class GroupTfidfTransformer(TfidfTransformer):

#     def fit(self, X, y=None):
#         assert y is not None

#         super(GroupTfidfTransformer, self).fit(X, y)

#         dtype = X.dtype if X.dtype in FLOAT_DTYPES else np.float64

#         n_samples, n_features = X.shape
#         n_samples, df = self._document_frequency(X, dtype=dtype)
#         n_group_samples, group_df = self._document_frequency(X, y, dtype=dtype)

#         # perform idf smoothing if required
#         df += int(self.smooth_idf)
#         group_df += int(self.smooth_idf)
#         n_samples += int(self.smooth_idf)
#         n_group_samples += int(self.smooth_idf)

#         # log+1 instead of log makes sure terms with zero idf don't get
#         # suppressed entirely.

#         idf = np.log(n_samples / df) + 1
#         self._idf_diag = sp.sparse.diags(idf, offsets=0,
#                                          shape=(n_features, n_features),
#                                          format='csr',
#                                          dtype=dtype)

#         group_idf = np.log(n_group_samples / group_df) + 1
#         self._group_idf = group_idf

#         return self

#     def transform(self, X, y=None, copy=True):
#         if y is not None:
#             if sp.sparse.issparse(X):
#                 X = X.toarray()
#             return X * self._group_idf

#         return super(GroupTfidfTransformer, self).transform(X, copy=copy)

class GroupTfidfVectorizer(TfidfVectorizer):
    def __init__(self, input='content', encoding='utf-8',
                 decode_error='strict', strip_accents=None, lowercase=True,
                 preprocessor=None, tokenizer=None, analyzer='word',
                 stop_words=None, token_pattern=r"(?u)\b\w\w+\b",
                 ngram_range=(1, 1), max_df=1.0, min_df=1,
                 max_features=None, vocabulary=None, binary=False,
                 dtype=np.float64, norm='l2', use_idf=True, smooth_idf=True,
                 sublinear_tf=False,
                 min_group_df=1):

        super(GroupTfidfVectorizer, self).__init__(
            input=input, encoding=encoding, decode_error=decode_error,
            strip_accents=strip_accents, lowercase=lowercase,
            preprocessor=preprocessor, tokenizer=tokenizer, analyzer=analyzer,
            stop_words=stop_words, token_pattern=token_pattern,
            ngram_range=ngram_range, max_df=max_df, min_df=min_df,
            max_features=max_features, vocabulary=vocabulary, binary=binary,
            dtype=dtype, norm=norm, use_idf=use_idf, smooth_idf=smooth_idf,
            sublinear_tf=sublinear_tf)

        self.min_group_df = min_group_df

    #     if use_gidf:
    #         self._tfidf = GroupTfidfTransformer(norm=norm, use_idf=use_idf,
    #                                             smooth_idf=smooth_idf,
    #                                             sublinear_tf=sublinear_tf)
    #     else:
    #         self._tfidf = TfidfTransformer(norm=norm, use_idf=use_idf,
    #                                        smooth_idf=smooth_idf,
    #                                        sublinear_tf=sublinear_tf)

    def fit_transform(self, raw_documents, y=None):
        self._check_params()
        X = super(TfidfVectorizer, self).fit_transform(raw_documents)
        self._tfidf.fit(X)
        X = self._tfidf.transform(X, copy=False)
        return self._limit_group_features(X, y)

    def _limit_group_features(self, X, y):
        n_samples, y_df = _document_group_frequency(X, y)

        if sp.sparse.issparse(X):
            X = X.tolil()

        if self.min_group_df > 1:
            X[y_df < self.min_group_df] = 0

        return X

class SyntacticGroupTfidfVectorizer(GroupTfidfVectorizer):
    """ Uses ngrams based on dependency syntax hierarchy """

    def _dep_ngrams(self, dep_df):
        """ analyzer can be one of [word, lemma, dep] """

        min_n, max_n = self.ngram_range
        ngrams = []

        unigrams = []

        sent_ids = dep_df['dep'].index.get_level_values(0).unique()
        for sent_id in sent_ids:
            sent_df = dep_df.loc[sent_id]

            tokens = sent_df[self.analyzer]
            heads = sent_df['g']
            for idx, head in heads.items():
                ngram = [tokens[idx]]
                unigrams.append((ngram, idx, head))

                if min_n == 1:
                    ngrams.append(tuple(ngram))

                while len(ngram) < max_n and head > 0 and head in tokens:
                    ngram.append(tokens[head])
                    head = heads[head]

                    if len(ngram) >= min_n:
                        ngrams.append(tuple(ngram))

        if len(ngrams) == 0:
            raise ValueError(str(sent_ids) + str(unigrams) + str(dep_df))

        return ngrams

    def build_analyzer(self):
        return self._dep_ngrams

class StratifiedSplitter():
    """ Split data based on task specifications """

    def __init__(self, n_problems=20, n_candidates=9, n_train_docs=7, n_test_docs=100, unk_ratio=0.5,
                 random_state=None):
        assert n_problems % 4 == 0

        self.n_problems = n_problems
        self.n_candidates = n_candidates
        self.n_train_docs = n_train_docs
        self.n_test_docs = n_test_docs
        self.unk_ratio = unk_ratio
        self.random_state = random_state

    def split(self, X, y=None, groups=None):
        labels, languages = y, groups
        rnd = np.random.RandomState(self.random_state)

        lang_set = sorted(set(languages))
        for lang in lang_set:
            lang_index = np.argwhere(languages == lang)[:, 0]

            for _ in range(self.n_problems // 4):
                candidates = labels[lang_index]
                candidate_set = list(set(candidates) - {'<UNK>', 'problem00020_candidate00000'})
                candidate_subset = rnd.choice(candidate_set, self.n_candidates, replace=False)

                train_index, test_index = [], []
                for candidate in candidate_subset:
                    candidate_index = lang_index[np.argwhere(candidates == candidate)[:, 0]]
                    rnd.shuffle(candidate_index)

                    train_index.extend(candidate_index[:self.n_train_docs])
                    test_index.extend(candidate_index[self.n_train_docs:self.n_train_docs + self.n_test_docs])

                unk_index = lang_index[np.argwhere(np.isin(labels[lang_index], candidate_subset, invert=True))[:, 0]]
                rnd.shuffle(unk_index)
                test_index.extend(unk_index[:int(len(test_index) * self.unk_ratio)])

                rnd.shuffle(train_index)
                rnd.shuffle(test_index)
                yield train_index, test_index

    def get_n_splits(self, X=None, y=None, groups=None):
        return self.n_problems

def merge_probleminfo(collectioninfo, path):
    filenames, labels, languages = [], [], []

    for probleminfo in collectioninfo:
        lang = probleminfo['language']
        prob_name = probleminfo['problem-name']

        prob_filenames, prob_labels = get_model_inputs(probleminfo['candidates'])
        prob_filenames.extend(probleminfo['unknown'])
        prob_labels.extend(get_test_labels(probleminfo['unknown'], collectioninfo, path))

        filenames.extend(prob_filenames)
        labels.extend([prob_name + '_' + f if f != '<UNK>' else '<UNK>' for f in prob_labels])
        languages.extend([lang] * len(prob_filenames))

    return np.array(filenames), np.array(labels), np.array(languages)

def f1_macro(estimator, X, y_true):
    warnings.filterwarnings('ignore', category=UndefinedMetricWarning)
    y_pred = estimator.predict(X)
    return f1_score(y_true, y_pred,
                    labels=estimator.classes_,
                    average='macro')

class VisualTransformer(BaseEstimator, TransformerMixin):  # TODO Fix ngram_range (2,2) and (3,3)
    def __init__(self, ngram_range=(1, 1)):
        self.vocab = None
        self.ngram_range = ngram_range

    def _ngrams(self, x):
        punct_set = ['.', ',', ':', ';', '-', '?', '!', '(', ')', '\'', '"']
        x_puncts = [t for t in x if t in punct_set]

        min_n, max_n = self.ngram_range

        ngrams = []
        for n in range(min_n, max_n + 1):
            ngrams.extend([tuple(x_puncts[i:i+n]) for i in range(len(x_puncts) - n + 1)])

        return ngrams

    def fit(self, X, y=None):
        ngrams = np.unique([ngram for x in X for ngram in self._ngrams(x)])
        self.vocab = list(ngrams)
        return self

    def transform(self, X):
        X_ngrams = [self._ngrams(x) for x in X]
        return [[np.sum([n == ngram for n in x]) for ngram in self.vocab] for x in X_ngrams]

pipeline = Pipeline([
    ('raw_features', FeatureExtractor()),
    ('feature_union', FeatureUnion(
        [
            # Select best character n-gram features  char_ngrams__text__ngram_range
            ('visual', Pipeline([
                ('text', FeatureSelector('tokens')),
                ('features', VisualTransformer(ngram_range=(1, 2))),
                ('scaler', MaxAbsScaler())
            ])),

            # Select best character n-gram features  char_ngrams__text__ngram_range
            ('char_ngrams', Pipeline([
                ('text', FeatureSelector('text')),
                ('tfidf', GroupTfidfVectorizer(analyzer='char', lowercase=False,
                                               ngram_range=(2, 4), use_idf=True, min_group_df=1)),
                ('scaler', MaxAbsScaler()),
                ('best', TruncatedSVD(n_components=150))
            ])),

            # Select best token n-gram features
            ('token_ngrams', Pipeline([
                ('text', FeatureSelector('tokens')),
                ('tfidf', GroupTfidfVectorizer(analyzer='word', ngram_range=(2, 2), min_df=5,
                                               lowercase=False, tokenizer=identity, use_idf=False)),
                ('scaler', MaxAbsScaler()),
                ('best', TruncatedSVD(n_components=150))
            ])),

            # Select best POS tag n-gram features
            ('pos_ngrams', Pipeline([
                ('text', FeatureSelector('pos_tags')),
                ('tfidf', GroupTfidfVectorizer(analyzer='word', tokenizer=identity, lowercase=False,
                                               ngram_range=(1, 4), use_idf=False, min_group_df=1)),
                ('scaler', MaxAbsScaler()),
                ('best', TruncatedSVD(n_components=150))
            ])),

            # Select best dependency tag sn-gram features
            ('dep_sngrams', Pipeline([
                ('text', FeatureSelector('syntax')),
                ('tfidf', SyntacticGroupTfidfVectorizer(analyzer='dep',
                                                        ngram_range=(1, 2), use_idf=False, min_group_df=3)),
                ('scaler', MaxAbsScaler()),
                ('best', TruncatedSVD(n_components=150))
            ])),

            # Select best word sn-gram features
            ('token_sngrams', Pipeline([
                ('text', FeatureSelector('syntax')),
                ('tfidf', SyntacticGroupTfidfVectorizer(analyzer='word',
                                                        ngram_range=(2, 3), use_idf=False, min_group_df=3)),
                ('scaler', MaxAbsScaler()),
                ('best', TruncatedSVD(n_components=150))
            ]))
        ],
    )),
    ('svc', ProbabilityThreshold(
        CalibratedClassifierCV(LinearSVC(), cv=5),
        pt=0.1)
     ),
])

def predict_all_problems(path, outpath, features, evaluate=False):
    start_time = time.time()

    lang_scores = {'en': [], 'fr': [], 'it': [], 'sp': []}

    collectioninfo = read_collection_info(path)
    for probleminfo in collectioninfo:
        print('Solving problem {} [{}]'.format(probleminfo['problem-name'], probleminfo['language']))

        train_filenames, train_labels = get_model_inputs(probleminfo['candidates'])
        test_filenames = probleminfo['unknown']

        pipeline.set_params(
            feature_union__transformer_weights={
                'visual': 1 if 'visual' in features else 0,
                'char_ngrams': 1 if 'char_ngrams' in features else 0,
                'token_ngrams': 1 if 'token_ngrams' in features else 0,
                'pos_ngrams': 1 if 'pos_ngrams' in features else 0,
                'dep_sngrams': 1 if 'dep_sngrams' in features else 0,
                'token_sngrams': 1 if 'token_sngrams' in features else 0
            }
        )

        pipeline.fit(train_filenames, train_labels)
        predictions = list(pipeline.predict(test_filenames))

        if evaluate:
            test_labels = get_test_labels(test_filenames, collectioninfo, path)
            prob_scores = precision_recall_fscore_support(test_labels, predictions,
                                                          labels=list(probleminfo['candidates'].keys()),
                                                          average='macro')
            lang_scores[probleminfo['language']].append(prob_scores)
            print(*prob_scores)

        save_output(outpath, probleminfo['problem-name'], test_filenames, predictions)

    if evaluate:
        p, r, f = [], [], []
        for lang, scores in lang_scores.items():
            lp, lr, lf = [], [], []
            for prob_scores in scores:
                lp.append(prob_scores[0])
                lr.append(prob_scores[1])
                lf.append(prob_scores[2])
            print(lang + ':', np.mean(lp), np.mean(lr), np.mean(lf))
            p.append(lp)
            r.append(lr)
            f.append(lf)
        print('Overall:', np.mean(p), np.mean(r), np.mean(f))

    print('elapsed time:', time.time() - start_time)

def grid_search(path, outpath, features, evaluate=False):
    collectioninfo = read_collection_info(path)
    filenames, labels, languages = merge_probleminfo(collectioninfo, path)

    n_problems = 40

    param_grids = []

    param_grids.append({
        'feature_union__visual__features__ngram_range': [(1, 1), (1, 2), (1, 3), (1, 4)]
    })

    for f in ['char_ngrams', 'pos_ngrams', 'dep_sngrams', 'token_sngrams']:
        ngram_range = [(i, j) for i in [1, 2] for j in [1, 2, 3, 4] if j >= i]
        param_grids.append({
            'feature_union__{}__tfidf__ngram_range'.format(f):  ngram_range,
            'feature_union__{}__tfidf__use_idf'.format(f):      [True, False],
            'feature_union__{}__tfidf__min_group_df'.format(f): range(1, 5),
        })

    param_grids.append({
        'svc__pt': [0.01, 0.05, 0.1, 0.15, 0.2]
    })

    for i, param_grid in enumerate(param_grids):
        print('Starting grid search {}'.format(i))

        cv = GridSearchCV(pipeline, param_grid, f1_macro, n_jobs=-1, iid=False,
                          cv=StratifiedSplitter(n_problems=n_problems, random_state=123),
                          error_score=np.nan, return_train_score=False, refit=False, verbose=3)
        cv.fit(X=filenames, y=labels, groups=languages)

        print('Search finished. Best mean macro f1-score is {} [{}]\n'.format(cv.best_score_, cv.best_index_))

        with open(outpath + os.sep + 'cv_results_{}.pkl'.format(i), 'wb') as f:
            pickle.dump(cv.cv_results_, f)

        with open(outpath + os.sep + 'best_params_{}.json'.format(i), 'w') as f:
            f.write(json.dumps(cv.best_params_))

        pipeline.set_params(**cv.best_params_)

def preprocess_input(raw_path, input_path):
    print(os.system('./scripts/runAll.sh {} {}'.format(raw_path, input_path)))

def main():
    parser = argparse.ArgumentParser(description='SVM Authorship Attribution Method')
    parser.add_argument('-r', type=str, help='Path to the main folder of the preprocessed input data')
    parser.add_argument('-i', type=str, help='Path to the main folder of a collection of attribution problems')
    parser.add_argument('-o', type=str, help='Path to an output folder')
    parser.add_argument('-f', type=str, default='visual,char_ngrams,token_ngrams,pos_ngrams,dep_sngrams,token_sngrams',
                        help='Features [visual,char_ngrams,token_ngrams,pos_ngrams,dep_sngrams,token_sngrams]')
    parser.add_argument('--eval', dest='eval', action='store_const',
                        const=True, default=False,
                        help='Evaluate results')
    parser.add_argument('--grid-search', dest='action', action='store_const',
                        const=grid_search, default=predict_all_problems,
                        help='Perform grid search')
    args = parser.parse_args()

    if not args.i:
        print('ERROR: The input folder is required')
        parser.exit(1)
    if not args.o:
        print('ERROR: The output folder is required')
        parser.exit(1)

    if args.r:
        preprocess_input(args.r, args.i)

    features = args.f.split(',')

    args.action(args.i, args.o, features, args.eval)

if __name__ == '__main__':
    main()

reason: import numpy,import scipy


text:
"""Source code copied and adjusted from gfootball.

GitHub:
https://github.com/google-research/football

License:
Copyright 2019 Google LLC
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.
"""

import numpy as np

def _get_cv2():
    try:
        import cv2  # pylint: disable=import-outside-toplevel
        return cv2
    except ImportError:
        raise ImportError(
            'Could not import cv2. Install GFootball and its dependencies '
            'according to https://github.com/google-research/football.'
        )

class _TextWriter:

    def __init__(
        self, frame, x, y=0, field_coords=False, color=(255, 255, 255)
    ):
        self._frame = frame
        if field_coords:
            x = 400 * (x + 1) - 5
            y = 695 * (y + 0.43)
        self._pos_x = int(x)
        self._pos_y = int(y) + 20
        self._color = color

    def write(self, text, scale_factor=1):
        cv2 = _get_cv2()
        font = cv2.FONT_HERSHEY_SIMPLEX  # pylint: disable=no-member
        text_pos = (self._pos_x, self._pos_y)
        font_scale = 0.5 * scale_factor
        line_type = 1
        cv2.putText(self._frame, text, text_pos, font, font_scale, self._color,  # pylint: disable=no-member
                    line_type)
        self._pos_y += int(20 * scale_factor)

def get_frame(trace):
    """Visualizes a state based on a GFootball trace dict."""
    cv2 = _get_cv2()
    frame = np.uint8(np.zeros((600, 800, 3)))
    corner1 = (0, 0)
    corner2 = (799, 0)
    corner3 = (799, 599)
    corner4 = (0, 599)
    line_color = (0, 255, 255)
    cv2.line(frame, corner1, corner2, line_color)  # pylint: disable=no-member
    cv2.line(frame, corner2, corner3, line_color)  # pylint: disable=no-member
    cv2.line(frame, corner3, corner4, line_color)  # pylint: disable=no-member
    cv2.line(frame, corner4, corner1, line_color)  # pylint: disable=no-member
    cv2.line(frame, (399, 0), (399, 799), line_color)  # pylint: disable=no-member
    writer = _TextWriter(
        frame,
        trace['ball'][0],
        trace['ball'][1],
        field_coords=True,
        color=(255, 0, 0))
    writer.write('B')
    for player_idx, player_coord in enumerate(trace['left_team']):
        writer = _TextWriter(
            frame,
            player_coord[0],
            player_coord[1],
            field_coords=True,
            color=(0, 255, 0))
        letter = 'H'
        if 'left_agent_controlled_player' in trace and player_idx in trace[
            'left_agent_controlled_player']:
            letter = 'X'
        writer.write(letter)
    for player_idx, player_coord in enumerate(trace['right_team']):
        writer = _TextWriter(
            frame,
            player_coord[0],
            player_coord[1],
            field_coords=True,
            color=(255, 255, 0))
        letter = 'A'
        if (
            'opponent_active' in trace and
            player_idx in trace['opponent_active']
        ):
            letter = 'Y'
        elif 'right_agent_controlled_player' in trace and player_idx in trace[
            'right_agent_controlled_player']:
            letter = 'Y'
        writer.write(letter)
    return frame

max_stars_repo_name: shoot-tree-search/sts

max_issues_repo_name: shoot-tree-search/sts

max_forks_repo_name: shoot-tree-search/sts

reason: import numpy


text:
# -*- coding: utf-8 -*-

# %% IMPORTS
# Built-in imports
from os import path

# Package imports
import numpy as np
import pytest

# PUMAS Cube imports
from pumas_cube import (
    calc_flux, export_to_txt, make_flux_plot, make_hist, make_scatter,
    read_cube_HDF5, run_multi_cube)

# %% GLOBALS
# Save the path to this directory
dirpath = path.dirname(__file__)

# Obtain paths to materials directory and input files
input_par = path.join(dirpath, 'data', 'input.par')
materials_dir = path.join(dirpath, 'data', 'materials')
MDF_file = path.join(materials_dir, 'materials.xml')
model_file = path.join(dirpath, 'data', 'test_model.txt')

# Create fixture for providing a temporary directory with input parameters
@pytest.fixture(scope='module')
def working_dir(tmpdir_factory):
    tmpdir = tmpdir_factory.mktemp('test_multi_cube')
    input_par_file = copy_input_par(tmpdir)
    return({'output_dir': tmpdir,
            'input_par': input_par_file})

# %% HELPER FUNCTIONS
# This function creates a test Rubik's cube model with 3 layers
def create_Rubiks_cube_model(filename):     # pragma: no cover
    # Define the number of unique X, Y and Z coordinates
    nx = 15
    ny = 15
    nz = 15

    # Define the boundaries between the layers
    f_layers = np.array([1/3, 2/3])

    # Determine the X, Y and Z coordinate ranges
    x_range = np.arange(0, nx*5, 5)
    y_range = np.arange(0, ny*5, 5)
    z_range = np.arange(0, -nz*5, -5)

    # Determine the actual grid
    Z, Y, X = np.meshgrid(z_range, y_range, x_range, indexing='ij')
    x_arr = X.ravel()
    y_arr = Y.ravel()
    z_arr = Z.ravel()

    # Determine the indices of the layers
    idx = np.array(np.rint(f_layers*nz)*nx*ny, dtype=int)

    # Create a rock_id array
    rock_id_arr = np.empty(nx*ny*nz, dtype=int)
    rock_id_arr[:idx[0]] = 0
    rock_id_arr[idx[0]:idx[1]] = 1
    rock_id_arr[idx[1]:] = 2

    # Create file
    with open(filename, 'w') as file:
        # Loop over all coordinates and rock_ids
        for x, y, z, rock_id in zip(x_arr, y_arr, z_arr, rock_id_arr):
            # Write line to file using a default density of 0.00
            file.write("{:6.2f} {:6.2f} {:6.2f} 0.00 {}\n".format(
                x, y, z, rock_id))

# This function copies the input parameters file to a specific location
def copy_input_par(dest_dir):
    # Read in the input parameters file
    data = np.genfromtxt(input_par, dtype=str, comments="%")
    data_dct = dict(data)

    # Obtain absolute path to provided dest_dir
    dest_dir = path.abspath(dest_dir)

    # Modify paths
    data_dct['output_dir'] = dest_dir
    data_dct['cube_model_0'] = model_file
    data_dct['MDF_filename'] = MDF_file
    data_dct['DEDX_dir'] = materials_dir

    # Create a new data_lst list
    data_lst = [' '.join([key, value]) for key, value in data_dct.items()]

    # Convert to a single string
    data_str = '\n'.join(data_lst)

    # Save to file
    dest_file = path.join(dest_dir, 'input.par')
    with open(dest_file, 'w') as file:
        file.write(data_str)

    # Return path to created file
    return(dest_file)

# %% PYTEST CLASSES AND FUNCTIONS
# Pytest class for run_multi_cube()-function
@pytest.mark.run_on_pass
class Test_run_multi_cube(object):
    # Test default way of calling this function
    def test_default(self, working_dir):
        # Run default model
        run_multi_cube(working_dir['input_par'],
                       N=100,
                       az_rng=(0, 45),
                       el_rng=(88, 90),
                       logE_rng=(2, 2.5))

    # Test using a new range of elevations that overlaps
    def test_el_rng_overlap(self, working_dir):
        # Run default model
        run_multi_cube(working_dir['input_par'],
                       N=100,
                       az_rng=(0, 45),
                       el_rng=(85, 90),
                       logE_rng=(2, 2.5))

    # Test using a new range of elevations that does not overlap
    def test_el_rng_no_overlap(self, working_dir):
        # Run default model
        run_multi_cube(working_dir['input_par'],
                       N=100,
                       az_rng=(0, 45),
                       el_rng=(81, 83),
                       logE_rng=(2, 2.5))

    # Test if providing an invalid input_par value raises an error
    def test_invalid_input_par(self):
        with pytest.raises(OSError):
            run_multi_cube('this_is_not_a_valid_file.txt',
                           N=100,
                           az_rng=(0, 45),
                           el_rng=(85, 90),
                           logE_rng=(2, 2.5))

# Pytest class for read_cube_HDF5()-function
class Test_read_cube_HDF5(object):
    # Test default way of calling this function
    def test_default(self, working_dir):
        # Obtain data
        data = read_cube_HDF5(output_dir=working_dir['output_dir'],
                              az_rng=None,
                              elevation=85,
                              logE_rng=None)

        # Check if certain expected data can be found in this dict
        assert (0, 2.0, 2.1) in data
        assert (0, 2.4, 2.5) in data
        assert (40, 2.0, 2.1) in data
        assert (40, 2.4, 2.5) in data
        assert (60, 2.0, 2.1) not in data
        assert (0, 2.0, 2.5) not in data

    # Test if specific azimuth and logE ranges can be used
    def test_ranges(self, working_dir):
        # Obtain data
        data = read_cube_HDF5(output_dir=working_dir['output_dir'],
                              az_rng=(0, 20),
                              elevation=85,
                              logE_rng=(2, 2.3))

        # Check if certain expected data can be found in this dict
        assert (0, 2.0, 2.1) in data
        assert (0, 2.4, 2.5) not in data
        assert (15, 2.0, 2.1) in data
        assert (15, 2.4, 2.5) not in data
        assert (40, 2.0, 2.1) not in data
        assert (0, 2.0, 2.5) not in data

    # Test if specific azimuth value can be used
    def test_az_value(self, working_dir):
        # Obtain data
        data = read_cube_HDF5(output_dir=working_dir['output_dir'],
                              az_rng=5,
                              elevation=85,
                              logE_rng=(2, 2.5))

        # Check if certain expected data can be found in this dict
        assert (0, 2.0, 2.1) not in data
        assert (0, 2.4, 2.5) not in data
        assert (5, 2.0, 2.1) in data
        assert (5, 2.4, 2.5) in data
        assert (40, 2.0, 2.1) not in data
        assert (0, 2.0, 2.5) not in data

    # Test if specific datasets can be requested
    def test_request_dsets(self, working_dir):
        # Obtain data
        data = read_cube_HDF5('position_xf', 'position_yf', 'position_zf',
                              output_dir=working_dir['output_dir'],
                              az_rng=5,
                              elevation=85,
                              logE_rng=(2, 2.5))

        # Check if certain expected data can be found in this dict
        assert 'position_xf' in data[(5, 2.0, 2.1)]
        assert 'energy_f' not in data[(5, 2.0, 2.1)]

    # Test if no datasets can be requested
    def test_no_dsets(self, working_dir):
        # Obtain data
        data = read_cube_HDF5(None,
                              output_dir=working_dir['output_dir'],
                              az_rng=5,
                              elevation=85,
                              logE_rng=(2, 2.5))

        # Check if certain expected data can be found in this dict
        assert list(data[(5, 2.0, 2.1)].keys()) == ['attrs']

    # Test if giving an invalid directory raises an error
    def test_invalid_dir(self):
        with pytest.raises(OSError):
            read_cube_HDF5(output_dir='invalid_dir',
                           az_rng=(0, 20),
                           elevation=85,
                           logE_rng=(2, 2.3))

    # Test if giving an invalid elevation raises an error
    def test_invalid_elevation(self, working_dir):
        with pytest.raises(ValueError):
            read_cube_HDF5(output_dir=working_dir['output_dir'],
                           az_rng=(0, 20),
                           elevation=80,
                           logE_rng=(2, 2.3))

    # Test if requested an invalid dataset raises an error
    def test_invalid_dset(self, working_dir):
        with pytest.raises(KeyError):
            read_cube_HDF5('invalid_dset_name',
                           output_dir=working_dir['output_dir'],
                           az_rng=(0, 20),
                           elevation=85,
                           logE_rng=(2, 2.3))

# Pytest class for calc_flux()-function
class Test_calc_flux(object):
    # Test default
    def test_default(self, working_dir):
        # Calculate flux
        _ = calc_flux(output_dir=working_dir['output_dir'],
                      az_rng=None,
                      el_rng=85,
                      logE_rng=None)

# Pytest class for make_flux_plot()-function
class Test_make_flux_plot(object):
    # Test default
    def test_default(self, working_dir):
        # Create plot
        make_flux_plot(output_dir=working_dir['output_dir'],
                       az_rng=None,
                       el_rng=85,
                       logE_rng=None)

    # Test if specific azimuth, elevation and logE ranges can be used
    def test_ranges(self, working_dir):
        # Create plot
        make_flux_plot(output_dir=working_dir['output_dir'],
                       az_rng=(0, 20),
                       el_rng=(85, 90),
                       logE_rng=(2, 2.3))

    # Test if specific azimuth value can be used
    def test_az_value(self, working_dir):
        # Create plot
        make_flux_plot(output_dir=working_dir['output_dir'],
                       az_rng=5,
                       el_rng=85,
                       logE_rng=(2, 2.5))

    # Test if the figure can be saved
    def test_savefig(self, working_dir):
        # Create plot
        filename = path.join(working_dir['output_dir'], 'test_flux.png')
        make_flux_plot(output_dir=working_dir['output_dir'],
                       az_rng=(0, 45),
                       el_rng=(85, 90),
                       logE_rng=None,
                       savefig=filename)

# Pytest class for make_hist()-function
class Test_make_hist(object):
    # Test default
    def test_default(self, working_dir):
        # Create plot
        make_hist('energy_f',
                  output_dir=working_dir['output_dir'],
                  az_rng=None,
                  el_rng=85,
                  logE_rng=None)

    # Test if specific azimuth, elevation and logE ranges can be used
    def test_ranges(self, working_dir):
        # Create plot
        make_hist('energy_f',
                  output_dir=working_dir['output_dir'],
                  az_rng=(0, 20),
                  el_rng=(80, 90),
                  logE_rng=(2, 2.3))

    # Test if specific azimuth value can be used
    def test_az_value(self, working_dir):
        # Create plot
        make_hist('energy_f',
                  output_dir=working_dir['output_dir'],
                  az_rng=5,
                  el_rng=85,
                  logE_rng=(2, 2.5))

    # Test if the figure can be saved
    def test_savefig(self, working_dir):
        # Create plot
        filename = path.join(working_dir['output_dir'], 'test_hist.png')
        make_hist('energy_f',
                  output_dir=working_dir['output_dir'],
                  az_rng=None,
                  el_rng=85,
                  logE_rng=None,
                  savefig=filename)

# Pytest class for make_scatter()-function
class Test_make_scatter(object):
    # Test default
    def test_default(self, working_dir):
        # Create plot
        make_scatter(output_dir=working_dir['output_dir'],
                     az_rng=None,
                     el_rng=85,
                     logE_rng=None)

    # Test if specific azimuth, elevation and logE ranges can be used
    def test_ranges(self, working_dir):
        # Create plot
        make_scatter(output_dir=working_dir['output_dir'],
                     az_rng=(0, 20),
                     el_rng=(80, 90),
                     logE_rng=(2, 2.3))

    # Test if specific azimuth value can be used
    def test_az_value(self, working_dir):
        # Create plot
        make_scatter(output_dir=working_dir['output_dir'],
                     az_rng=5,
                     el_rng=85,
                     logE_rng=(2, 2.5))

    # Test if the figure can be saved
    def test_savefig(self, working_dir):
        # Create plot
        filename = path.join(working_dir['output_dir'], 'test_scatter.png')
        make_scatter(output_dir=working_dir['output_dir'],
                     az_rng=None,
                     el_rng=85,
                     logE_rng=None,
                     savefig=filename)

# Pytest class for export_to_txt()-function
class Test_export_to_txt(object):
    # Test default
    def test_default(self, working_dir):
        # Export to file
        filename = path.join(working_dir['output_dir'], 'test.txt')
        export_to_txt(filename,
                      output_dir=working_dir['output_dir'],
                      az_rng=None,
                      el_rng=(80, 90),
                      logE_rng=None)

    # Test if specific azimuth, elevation and logE ranges can be used
    def test_ranges(self, working_dir):
        # Export to file
        filename = path.join(working_dir['output_dir'], 'test.txt')
        export_to_txt(filename,
                      output_dir=working_dir['output_dir'],
                      az_rng=(0, 20),
                      el_rng=(80, 90),
                      logE_rng=(2, 2.3))

    # Test if specific azimuth value can be used
    def test_az_value(self, working_dir):
        # Export to file
        filename = path.join(working_dir['output_dir'], 'test.txt')
        export_to_txt(filename,
                      output_dir=working_dir['output_dir'],
                      az_rng=5,
                      el_rng=85,
                      logE_rng=(2, 2.5))


text:
import tensorflow as tf
import matplotlib.pyplot as plt
from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
import numpy as np

def get_fixed_random(config, num_to_generate=100):
    seed_cont = tf.random.truncated_normal([num_to_generate, config.num_cont_noise])
    seed_cat = tf.math.mod(tf.range(0, num_to_generate), config.num_classes)
    return seed_cont, seed_cat

def generate_images(generator, z_input, c_input, config):
    if not config.conditional:
        c_input = None
    predictions = generator(z_input, c_input, training=False)
    gen_img = _data2plot(predictions, config)

    return gen_img

def _fig2data(fig):
    """
    @brief Convert a Matplotlib figure to a 3d numpy array with RGB channels and return it
    @param fig a matplotlib figure
    @return a numpy 3D array of RGB values
    """
    # draw the renderer
    canvas = FigureCanvas(fig)
    canvas.draw()

    width, height = fig.get_size_inches() * fig.get_dpi()

    image = np.fromstring(canvas.tostring_rgb(), dtype=np.uint8).reshape(height.astype(np.int32), width.astype(np.int32), 3)

    return image

def _data2plot(array, config):
    fig = plt.figure(figsize=(10, 10))
    if config.dataset in ['mnist', 'fashion_mnist']:  # color channel of the dataset
        for i in range(array.shape[0]):
            plt.subplot(10, 10, i + 1)  # 10*10 subplots
            plt.imshow(array[i, :, :, 0] * 127.5 + 127.5, cmap='gray')
            plt.axis('off')
    else:
        for i in range(array.shape[0]):
            plt.subplot(10, 10, i + 1)  # 10*10 subplots
            plt.imshow(array[i, :, :, :])
            plt.axis('off')
    # return fig
    return _fig2data(fig)


text:
from suspect import MRSData

import os
import numpy

spar_types = {
    "floats": ["ap_size", "lr_size", "cc_size", "ap_off_center", "lr_off_center",
               "cc_off_center", "ap_angulation", "lr_angulation", "cc_angulation",
               "image_plane_slice_thickness", "slice_distance", "spec_col_lower_val",
               "spec_col_upper_val", "spec_row_lower_val", "spec_row_upper_val"],
    "integers": ["samples", "rows", "synthesizer_frequency", "offset_frequency",
                 "sample_frequency", "echo_nr", "mix_number", "t0_mul_direction",
                 "echo_time", "repetition_time", "averages", "volumes",
                 "volume_selection_method", "nr_of_slices_for_multislice",
                 "spec_num_col", "spec_num_row", "num_dimensions", "TSI_factor",
                 "spectrum_echo_time", "spectrum_inversion_time", "image_chemical_shift",
                 "t0_mu1_direction"],
    "strings": ["scan_id", "scan_date", "patient_name", "patient_birth_date",
                "patient_position", "patient_orientation", "nucleus",
                "volume_selection_enable", "phase_encoding_enable", "t1_measurement_enable",
                "t2_measurement_enable", "time_series_enable", "Spec.image in plane transf",
                "spec_data_type", "spec_sample_extension", "spec_col_extension",
                "spec_row_extension", "echo_acquisition", "resp_motion_comp_technique",
                "de_coupling", "equipment_sw_verions", "examination_name"],
}

def load_sdat(sdat_filename, spar_filename=None):
    # if the spar filename is not supplied, assume it is in the same folder as
    # the sdat and only differs in the extension
    if spar_filename is None:
        path, ext = os.path.splitext(sdat_filename)
        # match the capitalisation of the sdat extension
        if ext == ".SDAT":
            spar_filename = path + ".SPAR"
        elif ext == ".sdat":
            spar_filename = path + ".spar"

    with open(spar_filename, 'r') as fin:
        parameter_dict = {}
        for line in fin:
            # ignore empty lines and comments starting with !
            if line != "\n" and not line.startswith("!"):
                key, value = map(str.strip, line.split(":", 1))
                if key in spar_types["floats"]:
                    parameter_dict[key] = float(value)
                elif key in spar_types["integers"]:
                    parameter_dict[key] = int(value)
                elif key in spar_types["strings"]:
                    parameter_dict[key] = value
                else:
                    pass
                    #print("{} : {}".format(key, value))

    with open(sdat_filename, 'rb') as fin:
        raw_bytes = fin.read()

    floats = _vax_to_ieee_single_float(raw_bytes)
    data_iter = iter(floats)
    complex_iter = (complex(r, -i) for r, i in zip(data_iter, data_iter))
    raw_data = numpy.fromiter(complex_iter, "complex64")
    raw_data = numpy.reshape(raw_data, (parameter_dict["rows"], parameter_dict["samples"])).squeeze()
    return MRSData(raw_data,
                   5e-4,
                   parameter_dict["synthesizer_frequency"] * 1e-6,
                   te=parameter_dict["echo_time"])

def _vax_to_ieee_single_float(data):
    """Converts a float in Vax format to IEEE format.

    data should be a single string of chars that have been read in from
    a binary file. These will be processed 4 at a time into float values.
    Thus the total number of byte/chars in the string should be divisible
    by 4.

    Based on VAX data organization in a byte file, we need to do a bunch of
    bitwise operations to separate out the numbers that correspond to the
    sign, the exponent and the fraction portions of this floating point
    number

    role :      S        EEEEEEEE      FFFFFFF      FFFFFFFF      FFFFFFFF
    bits :      1        2      9      10                               32
    bytes :     byte2           byte1               byte4         byte3

    """
    f = []
    nfloat = int(len(data) / 4)
    for i in range(nfloat):

        byte2 = data[0 + i*4]
        byte1 = data[1 + i*4]
        byte4 = data[2 + i*4]
        byte3 = data[3 + i*4]

        # hex 0x80 = binary mask 10000000
        # hex 0x7f = binary mask 01111111

        sign = (byte1 & 0x80) >> 7
        expon = ((byte1 & 0x7f) << 1) + ((byte2 & 0x80) >> 7)
        fract = ((byte2 & 0x7f) << 16) + (byte3 << 8) + byte4

        if sign == 0:
            sign_mult = 1.0
        else:
            sign_mult = -1.0

        if 0 < expon:
            # note 16777216.0 == 2^24
            val = sign_mult * (0.5 + (fract/16777216.0)) * pow(2.0, expon - 128.0)
            f.append(val)
        elif expon == 0 and sign == 0:
            f.append(0)
        else:
            f.append(0)
            # may want to raise an exception here ...

    return f

max_stars_repo_name: bennyrowland/suspect

max_issues_repo_name: bennyrowland/suspect

max_forks_repo_name: bennyrowland/suspect


text:
#!/usr/bin/env python3
"""From RAFFT output, build a kinetic model. Starting from only unfolded
structures, it generates a folding kinetic trajectory.

Usage:
rafft_kin rafft.out --plot

"""

from rafft.utils import paired_positions, parse_rafft_output
from numpy import array, zeros, exp, diag
from scipy.linalg import eig, inv
import numpy as np
import matplotlib.pyplot as plt
from copy import deepcopy

def plot_traj(trajectory, struct_list, times, font_size, width, height,
              show_thres, out_file=None):
    """plot trajectory
    """
    trajectory = array(trajectory).real

    plt.rcParams["font.family"] = "serif"
    plt.rcParams["font.size"] = font_size
    plt.rcParams["figure.figsize"] = width, height
    left, width = 0.10, 0.88
    bottom, height = 0.10, 0.88
    rect_scatter = [left, bottom, width, height]
    fig = plt.figure(1)
    kin_f = fig.add_axes(rect_scatter)
    kin_f.grid(True, color="grey", linestyle="--", linewidth=0.2)

    kin_f.set_xlim([times[0], times[-1]])

    for si, struct in enumerate(struct_list):
        if any(trajectory[:, si] > show_thres):
            kin_f.plot(times, trajectory[:, si], alpha=0.8, label=si)

    kin_f.set_xscale("log")
    kin_f.legend(ncol=2, fontsize=int(font_size * 0.8))
    if out_file is not None:
        plt.savefig(out_file, dpi=300, transparent=True)
    else:
        plt.show()

def get_connected_prev(cur_struct, prev_pos):
    "get the connected structures"
    cur_pairs = set(paired_positions(cur_struct.str_struct))
    res = []
    for si, struct in enumerate(prev_pos):
        pairs = set(paired_positions(struct.str_struct))
        if len(pairs - cur_pairs) == 0:
            res += [si]
    return res

def get_connected_next(cur_struct, prev_pos):
    "get the connected structures"
    cur_pairs = set(paired_positions(cur_struct.str_struct))
    for si, struct in enumerate(prev_pos):
        pairs = set(paired_positions(struct.str_struct))
        if len(cur_pairs - pairs) == 0:
            return True

def get_transition_mat(fast_paths, nb_struct, struct_map):
    transition_mat = zeros((nb_struct, nb_struct), dtype=np.longdouble)
    KT = 0.61

    for step_i, fold_step in enumerate(fast_paths):
        for struct in fold_step:
            # get all structures connected from previous step
            lprev_co = get_connected_prev(struct, fast_paths[step_i - 1])
            map_cur, cur_nrj = struct_map[struct.str_struct]

            for si in lprev_co:
                prev_st = fast_paths[step_i-1][si]
                map_prev, prev_nrj = struct_map[prev_st.str_struct]
                delta_nrj = cur_nrj - prev_nrj
                if map_cur != map_prev:
                    # # M_ij = k(i -> j)
                    transition_mat[map_prev, map_cur] = min(1.0, exp(-delta_nrj/KT))
                    transition_mat[map_cur, map_prev] = min(1.0, exp(delta_nrj/KT))

    # normalize input and output flows
    for si in range(nb_struct):
        transition_mat[si, si] = -transition_mat[si,:].sum()

    return transition_mat

def kinetics(fast_paths, max_time, n_steps, initial_pop=None):
    """
    input:
    fast_paths = list of list of (structure, energy)
    initial_pop = list of float giving the initial population of all/few structures

    output:
    trajectory = trajectory of population
    times = time steps
    struct_list = list of unique structures
    str_equi_pop = final population
    """
    seen = set()
    struct_list = []
    for el in fast_paths:
        for struct in el:
            if struct.str_struct not in seen:
                seen.add(struct.str_struct)
                struct_list += [struct]

    # create a non-redundant mapping
    struct_map = {struct.str_struct: (si, struct.energy) for si, struct in enumerate(struct_list)}
    nb_struct = len(struct_list)
    transition_mat = get_transition_mat(fast_paths, nb_struct, struct_map)

    # initialize the kinetic
    if initial_pop is None:
        # starts with the unfolded state
        init_pop = array([1.0] + [0.0 for _ in range(nb_struct-1)], dtype=np.longdouble)
    else:
        init_pop = array([0.0 for _ in range(nb_struct)], dtype=np.longdouble)
        for p, w in initial_pop:
            init_pop[p] = w

    trajectory = [deepcopy(init_pop)]

    # diagonalize the transition matrix to solve the system
    V, W = eig(transition_mat.T, check_finite=True)
    iW = inv(W)

    time_step = max_time / n_steps

    times = [exp(-4)]
    for st in range(n_steps):
        time = exp(time_step * st-4)
        times += [time]
        tmp_pop = W @ diag(exp(V * time)) @ (iW @ init_pop)
        trajectory += [tmp_pop.real/tmp_pop.real.sum()]

    equi_pop = trajectory[-1]

    # g = nx.from_numpy_matrix(transition_mat)
    # nx.draw(g)
    # plt.show()

    str_equi_pop = [(struct.str_struct, struct.energy, ep, struct_map[struct.str_struct][0]) for struct, ep in zip(struct_list, equi_pop.real)]
    return trajectory, times, struct_list, str_equi_pop

max_stars_repo_name: strevol-mpi-mis/RAFFT

max_issues_repo_name: strevol-mpi-mis/RAFFT

max_forks_repo_name: strevol-mpi-mis/RAFFT

reason: import numpy,from numpy,from scipy


text:
"""
A script to collect a batch of human demonstrations that can be used
to generate a learning curriculum (see `demo_learning_curriculum.py`).

The demonstrations can be played back using the `playback_demonstrations_from_pkl.py`
script.
"""

import os
import shutil
import time
import argparse
import datetime
import h5py
from glob import glob
import numpy as np
import json

import robosuite as suite
from robosuite import load_controller_config
from robosuite.wrappers import DataCollectionWrapper, VisualizationWrapper
from robosuite.utils.input_utils import input2action

def collect_human_trajectory(env, device, arm, env_configuration):
    """
    Use the device (keyboard or SpaceNav 3D mouse) to collect a demonstration.
    The rollout trajectory is saved to files in npz format.
    Modify the DataCollectionWrapper wrapper to add new fields or change data formats.

    Args:
        env (MujocoEnv): environment to control
        device (Device): to receive controls from the device
        arms (str): which arm to control (eg bimanual) 'right' or 'left'
        env_configuration (str): specified environment configuration
    """

    env.reset()

    # ID = 2 always corresponds to agentview
    env.render()

    is_first = True

    task_completion_hold_count = -1  # counter to collect 10 timesteps after reaching goal
    device.start_control()

    # Loop until we get a reset from the input or the task completes
    while True:
        # Set active robot
        active_robot = env.robots[0] if env_configuration == "bimanual" else env.robots[arm == "left"]

        # Get the newest action
        action, grasp = input2action(
            device=device,
            robot=active_robot,
            active_arm=arm,
            env_configuration=env_configuration
        )

        # If action is none, then this a reset so we should break
        if action is None:
            break

        # Run environment step
        env.step(action)

        if is_first:
            is_first = False

            # We grab the initial model xml and state and reload from those so that
            # we can support deterministic playback of actions from our demonstrations.
            # This is necessary due to rounding issues with the model xml and with
            # env.sim.forward(). We also have to do this after the first action is
            # applied because the data collector wrapper only starts recording
            # after the first action has been played.
            initial_mjstate = env.sim.get_state().flatten()
            xml_str = env.sim.model.get_xml()
            env.reset_from_xml_string(xml_str)
            env.sim.reset()
            env.sim.set_state_from_flattened(initial_mjstate)
            env.sim.forward()

        env.render()

        # Also break if we complete the task
        if task_completion_hold_count == 0:
            break

        # state machine to check for having a success for 10 consecutive timesteps
        if env._check_success():
            if task_completion_hold_count > 0:
                task_completion_hold_count -= 1  # latched state, decrement count
            else:
                task_completion_hold_count = 10  # reset count on first success timestep
        else:
            task_completion_hold_count = -1  # null the counter if there's no success

    # cleanup for end of data collection episodes
    env.close()

def gather_demonstrations_as_hdf5(directory, out_dir, env_info):
    """
    Gathers the demonstrations saved in @directory into a
    single hdf5 file.

    The strucure of the hdf5 file is as follows.

    data (group)
        date (attribute) - date of collection
        time (attribute) - time of collection
        repository_version (attribute) - repository version used during collection
        env (attribute) - environment name on which demos were collected

        demo1 (group) - every demonstration has a group
            model_file (attribute) - model xml string for demonstration
            states (dataset) - flattened mujoco states
            actions (dataset) - actions applied during demonstration

        demo2 (group)
        ...

    Args:
        directory (str): Path to the directory containing raw demonstrations.
        out_dir (str): Path to where to store the hdf5 file.
        env_info (str): JSON-encoded string containing environment information,
            including controller and robot info
    """

    hdf5_path = os.path.join(out_dir, "demo.hdf5")
    f = h5py.File(hdf5_path, "w")

    # store some metadata in the attributes of one group
    grp = f.create_group("data")

    num_eps = 0
    env_name = None  # will get populated at some point

    for ep_directory in os.listdir(directory):

        state_paths = os.path.join(directory, ep_directory, "state_*.npz")
        states = []
        actions = []

        for state_file in sorted(glob(state_paths)):
            dic = np.load(state_file, allow_pickle=True)
            env_name = str(dic["env"])

            states.extend(dic["states"])
            for ai in dic["action_infos"]:
                actions.append(ai["actions"])

        if len(states) == 0:
            continue

        # Delete the first actions and the last state. This is because when the DataCollector wrapper
        # recorded the states and actions, the states were recorded AFTER playing that action.
        del states[-1]
        del actions[0]

        num_eps += 1
        ep_data_grp = grp.create_group("demo_{}".format(num_eps))

        # store model xml as an attribute
        xml_path = os.path.join(directory, ep_directory, "model.xml")
        with open(xml_path, "r") as f:
            xml_str = f.read()
        ep_data_grp.attrs["model_file"] = xml_str

        # write datasets for states and actions
        ep_data_grp.create_dataset("states", data=np.array(states))
        ep_data_grp.create_dataset("actions", data=np.array(actions))

    # write dataset attributes (metadata)
    now = datetime.datetime.now()
    grp.attrs["date"] = "{}-{}-{}".format(now.month, now.day, now.year)
    grp.attrs["time"] = "{}:{}:{}".format(now.hour, now.minute, now.second)
    grp.attrs["repository_version"] = suite.__version__
    grp.attrs["env"] = env_name
    grp.attrs["env_info"] = env_info

    f.close()

if __name__ == "__main__":
    # Arguments
    parser = argparse.ArgumentParser()
    parser.add_argument(
        "--directory",
        type=str,
        default=os.path.join(suite.models.assets_root, "demonstrations"),
    )
    parser.add_argument("--environment", type=str, default="Lift")
    parser.add_argument("--robots", nargs="+", type=str, default="Panda", help="Which robot(s) to use in the env")
    parser.add_argument("--config", type=str, default="single-arm-opposed",
                        help="Specified environment configuration if necessary")
    parser.add_argument("--arm", type=str, default="right", help="Which arm to control (eg bimanual) 'right' or 'left'")
    parser.add_argument("--camera", type=str, default="agentview", help="Which camera to use for collecting demos")
    parser.add_argument("--controller", type=str, default="OSC_POSE",
                        help="Choice of controller. Can be 'IK_POSE' or 'OSC_POSE'")
    parser.add_argument("--device", type=str, default="keyboard")
    parser.add_argument("--pos-sensitivity", type=float, default=1.5, help="How much to scale position user inputs")
    parser.add_argument("--rot-sensitivity", type=float, default=1.5, help="How much to scale rotation user inputs")
    args = parser.parse_args()

    # Get controller config
    controller_config = load_controller_config(default_controller=args.controller)

    # Create argument configuration
    config = {
        "env_name": args.environment,
        "robots": args.robots,
        "controller_configs": controller_config,
    }

    # Check if we're using a multi-armed environment and use env_configuration argument if so
    if "TwoArm" in args.environment:
        config["env_configuration"] = args.config

    # Create environment
    env = suite.make(
        **config,
        has_renderer=True,
        has_offscreen_renderer=False,
        render_camera=args.camera,
        ignore_done=True,
        use_camera_obs=False,
        reward_shaping=True,
        control_freq=20,
    )

    # Wrap this with visualization wrapper
    env = VisualizationWrapper(env)

    # Grab reference to controller config and convert it to json-encoded string
    env_info = json.dumps(config)

    # wrap the environment with data collection wrapper
    tmp_directory = "/tmp/{}".format(str(time.time()).replace(".", "_"))
    env = DataCollectionWrapper(env, tmp_directory)

    # initialize device
    if args.device == "keyboard":
        from robosuite.devices import Keyboard

        device = Keyboard(pos_sensitivity=args.pos_sensitivity, rot_sensitivity=args.rot_sensitivity)
        env.viewer.add_keypress_callback("any", device.on_press)
        env.viewer.add_keyup_callback("any", device.on_release)
        env.viewer.add_keyrepeat_callback("any", device.on_press)
    elif args.device == "spacemouse":
        from robosuite.devices import SpaceMouse

        device = SpaceMouse(pos_sensitivity=args.pos_sensitivity, rot_sensitivity=args.rot_sensitivity)
    else:
        raise Exception(
            "Invalid device choice: choose either 'keyboard' or 'spacemouse'."
        )

    # make a new timestamped directory
    t1, t2 = str(time.time()).split(".")
    new_dir = os.path.join(args.directory, "{}_{}".format(t1, t2))
    os.makedirs(new_dir)

    # collect demonstrations
    while True:
        collect_human_trajectory(env, device, args.arm, args.config)
        gather_demonstrations_as_hdf5(tmp_directory, new_dir, env_info)

max_stars_repo_name: wangcongrobot/robosuite

max_issues_repo_name: wangcongrobot/robosuite

max_forks_repo_name: wangcongrobot/robosuite


text:
"""
SABR Volatility calibration

"""

import copy
import math
import numpy as np
from optionmodels.analyticalmethods import AnalyticalMethods
from optionmodels.models_params import sabr_params_dict
from optionmodels.utils import Utils
# pylint: disable=invalid-name

class SABRVolatility():
    """
    Stochastic, Alpha, Beta, Rho model

    Extension of Black 76 model to include an easily implementable
    stochastic volatility model

    Beta will typically be chosen a priori according to how traders
    observe market prices:
        e.g. In FX markets, standard to assume lognormal terms, Beta = 1
             In some Fixed Income markets traders prefer to assume
             normal terms, Beta = 0

    Alpha will need to be calibrated to ATM volatility

    """

    def __init__(self, **kwargs):

        # Import dictionary of default parameters
        self.default_dict = copy.deepcopy(sabr_params_dict)

        # Store initial inputs
        inputs = {}
        for key, value in kwargs.items():
            inputs[key] = value

        # Initialise system parameters
        params = Utils.init_sabr_params(inputs)

        self.params = params

    def calibrate(self, **kwargs):
        """
        Run the SABR calibration

        Returns
        -------
        Float
            Black-76 equivalent SABR volatility and price.

        """
        if 'refresh' not in kwargs:
            kwargs['refresh'] = True

        # Update pricing input parameters to default if not supplied
        if 'refresh' in kwargs and kwargs['refresh']:
            params = Utils.init_sabr_params(kwargs)
            F = params['F']
            K = params['K']
            T = params['T']
            r = params['r']
            atmvol = params['atmvol']
            beta = params['beta']
            volvol = params['volvol']
            rho = params['rho']
            output_flag = params['output_flag']
            option = params['option']

        black_vol = self._alpha_sabr(
            F, K, T, beta, volvol, rho, self._find_alpha(
                F=F, T=T, atmvol=atmvol, beta=beta, volvol=volvol, rho=rho,))

        black_price = AnalyticalMethods.black_76(
            F=F, K=K, T=T, r=r, sigma=black_vol, option=option, refresh=True)

        output_dict = {
            'vol':black_vol,
            'price':black_price,
            'all':{
                'Price':black_price,
                'Vol':black_vol
                }
            }

        return output_dict.get(output_flag, "Please enter a valid output flag")

#        if output_flag == 'vol':
#            return black_vol

#        elif output_flag == 'price':
#            return black_price

#        elif output_flag == 'all':
#            return {'Price':black_price,
#                    'Vol':black_vol}

    @staticmethod
    def _alpha_sabr(F, K, T, beta, volvol, rho, alpha):
        """
        The SABR skew vol function

        Parameters
        ----------
        Alpha : Float
            Alpha value.

        Returns
        -------
        result : Float
            Black-76 equivalent SABR volatility.

        """

        dSABR = np.zeros(4)
        dSABR[1] = (
            alpha
            / ((F * K) ** ((1 - beta) / 2)
               * (1
                  + (((1 - beta) ** 2) / 24)
                  * (np.log(F / K) ** 2)
                  + ((1 - beta) ** 4 / 1920)
                  * (np.log(F / K) ** 4))))

        if abs(F - K) > 10 ** -8:
            sabrz = (volvol / alpha) * (F * K) ** (
                (1 - beta) / 2) * np.log(F / K)
            y = (np.sqrt(1 - 2 * rho * sabrz + (
                sabrz ** 2)) + sabrz - rho) / (1 - rho)
            if abs(y - 1) < 10 ** -8:
                dSABR[2] = 1
            elif y > 0:
                dSABR[2] = sabrz / np.log(y)
            else:
                dSABR[2] = 1
        else:
            dSABR[2] = 1

        dSABR[3] = (1 + ((((1 - beta) ** 2 / 24) * alpha ** 2 / (
            (F * K) ** (1 - beta))) + (
                0.25 * rho * beta * volvol * alpha) / (
                    (F * K) ** ((1 - beta) / 2)) + (
                        2 - 3 * rho ** 2) * volvol ** 2 / 24) * T)

        result = dSABR[1] * dSABR[2] * dSABR[3]

        return result

    @classmethod
    def _find_alpha(cls, F, T, atmvol, beta, volvol, rho):
        """
        Find alpha feeding values to _cube_root method.

        Returns
        -------
        result : Float
            Smallest positive root.

        """
        # Alpha is a function of atm vol etc

        alpha = cls._cube_root(
            ((1 - beta) ** 2 * T / (24 * F ** (2 - 2 * beta))),
            (0.25 * rho * volvol * beta * T / (F ** (1 - beta))),
            (1 + (2 - 3 * rho ** 2) / 24 * volvol ** 2 * T),
            (-atmvol * F ** (1 - beta)))

        return alpha

    @classmethod
    def _cube_root(cls, cubic, quadratic, linear, constant):
        """
        Finds the smallest positive root of the input cubic polynomial
        algorithm from Numerical Recipes

        Parameters
        ----------
        cubic : Float
            3rd order term of input polynomial.
        quadratic : Float
            2nd order term of input polynomial.
        linear : Float
            Linear term of input polynomial.
        constant : Float
            Constant term of input polynomial.

        Returns
        -------
        result : Float
            Smallest positive root.

        """
        a = quadratic / cubic
        b = linear / cubic
        C = constant / cubic
        Q = (a ** 2 - 3 * b) / 9
        r = (2 * a ** 3 - 9 * a * b + 27 * C) / 54
        roots = np.zeros(4)

        if r ** 2 - Q ** 3 >= 0:
            cap_A = -np.sign(r) * (abs(r) + np.sqrt(
                r ** 2 - Q ** 3)) ** (1 / 3)
            if cap_A == 0:
                cap_B = 0
            else:
                cap_B = Q / cap_A
            result = cap_A + cap_B - a / 3
        else:
            theta = cls._arccos(r / Q ** 1.5)

            # The three roots
            roots[1] = - 2 * np.sqrt(Q) * math.cos(
                theta / 3) - a / 3
            roots[2] = - 2 * np.sqrt(Q) * math.cos(
                theta / 3 + 2.0943951023932) - a / 3
            roots[3] = - 2 * np.sqrt(Q) * math.cos(
                theta / 3 - 2.0943951023932) - a / 3

            # locate that one which is the smallest positive root
            # assumes there is such a root (true for SABR model)
            # there is always a small positive root

            if roots[1] > 0:
                result = roots[1]
            elif roots[2] > 0:
                result = roots[2]
            elif roots[3] > 0:
                result = roots[3]

            if roots[2] > 0 and roots[2] < result:
                result = roots[2]

            if roots[3] > 0 and roots[3] < result:
                result = roots[3]

        return result

    @staticmethod
    def _arccos(y):
        """
        Inverse Cosine method

        Parameters
        ----------
        y : Float
            Input value.

        Returns
        -------
        result : Float
            Arc Cosine of input value.

        """
        result = np.arctan(-y / np.sqrt(-y * y + 1)) + 2 * np.arctan(1)

        return result

max_stars_repo_name: GBERESEARCH/optionmodels

max_issues_repo_name: GBERESEARCH/optionmodels

max_forks_repo_name: GBERESEARCH/optionmodels


text:
from __future__ import absolute_import
from __future__ import division
from os import sep

from architecture.crowd_count import CrowdCounter
import architecture.network as network
from manage_data.utils import mkdir_if_missing as mkdir_if_missing
from architecture.utils import save_results as save_results
from architecture.cam import cam as apply_cam
import numpy as np
import torch
import pandas as pd
import os.path as osp

def evaluate_model(trained_model, data_loader, model, save_test_results = False, plot_save_dir = "/tmp/", den_scale_factor = 1e3, cam=False):
    net = CrowdCounter(model = model)
    network.load_net(trained_model, net)
    net.cuda()
    net.eval()
    mae = 0.0
    mse = 0.0
    if cam:
        net.net.cam = True
    res = {'img': [], 'et': [], 'gt': [], 'abs_err': [], 'err^2': [], 'err/gt': []}

    for blob in data_loader:
        im_data = blob['data']
        gt_data = blob['gt_density']
        idx_data = blob['idx']
        im_data_norm = im_data / 127.5 - 1.
        density_map = net(im_data_norm)
        if cam in ['grad-cam', 'all']:
            apply_cam(net, blob, density_map, osp.join(plot_save_dir, 'grad-cam'), 'grad-cam')
        if cam in ['grad-cam++', 'all']:
            apply_cam(net, blob, density_map, osp.join(plot_save_dir, 'grad-cam++'), 'grad-cam++')
        #if cam in ['eigen-cam', 'all']:
        #    apply_cam(net, blob, density_map, osp.join(plot_save_dir, 'eigen-cam'), 'eigen-cam')
        if cam in ['layer-cam', 'all']:
            apply_cam(net, blob, density_map, osp.join(plot_save_dir, 'layer-cam'), 'layer-cam')
        density_map = density_map.data.cpu().numpy()
        density_map /= den_scale_factor
        gt_count = np.sum(gt_data)
        et_count = np.sum(density_map)
        mae += abs(gt_count-et_count)
        mse += ((gt_count-et_count)*(gt_count-et_count))

        res['img'].append(idx_data)
        res['et'].append(et_count)
        res['gt'].append(gt_count)
        res['abs_err'].append(abs(gt_count-et_count))
        res['err^2'].append((gt_count-et_count)*(gt_count-et_count))
        res['err/gt'].append(abs(gt_count-et_count)/gt_count)

        if save_test_results:
            print("Plotting results")
            mkdir_if_missing(plot_save_dir)
            save_results(im_data, gt_data, density_map, idx_data, plot_save_dir)
    mae = mae/data_loader.get_num_samples()
    mse = np.sqrt(mse/data_loader.get_num_samples())
    df = pd.DataFrame(res).sort_values('err/gt', axis=0, ascending=False)
    df.to_csv(osp.join(plot_save_dir, 'error_by_img.csv'), sep=',', index=False, float_format='%.5f')

    return mae,mse


text:
"""
	Make a file with Hi-C interactions for every TAD. Each side of the interaction is considered a separate regulatory element that can be gained or lost.

	1. Filter for all interactions that do not take place entirely within a TAD
	2. Write these interactions per TAD to a file.

"""

import sys
import numpy as np

hiCFile = sys.argv[1]
tadFile = sys.argv[2]
outFile = sys.argv[3]

#1. Read the TADs

tads = []
with open(tadFile, 'r') as inF:

	for line in inF:
		line = line.strip()
		splitLine = line.split("\t")

		tads.append([splitLine[0], int(splitLine[1]), int(splitLine[2])])

	tads = np.array(tads, dtype="object")

#2. Read the Hi-C interactions and link to TADs.
tadInteractions = dict()
with open(hiCFile, 'r') as inF:
	lineCount = 0
	for line in inF:

		line = line.strip()
		splitLine = line.split("\t")

		region1 = splitLine[0]
		region2 = splitLine[1]

		regionChr = region1.split("_")[0]
		region1Pos = region1.split("_")[1]
		region2Pos = region2.split("_")[1]

		#1. Map both regions to a TAD.

		tadChrSubset = tads[tads[:,0] == regionChr]
		#Match region 1
		matchingStart = tadChrSubset[:,1] <= int(region1Pos)
		matchingEnd = tadChrSubset[:,2] >= int(region1Pos)

		matchingPos = matchingStart * matchingEnd
		if len(tadChrSubset[matchingPos]) < 1:
			continue
		matchingTad1 = tadChrSubset[matchingPos][0]

		#Match region 1
		matchingStart = tadChrSubset[:,1] <= int(region2Pos)
		matchingEnd = tadChrSubset[:,2] >= int(region2Pos)

		matchingPos = matchingStart * matchingEnd
		if len(tadChrSubset[matchingPos]) < 1:
			continue
		matchingTad2 = tadChrSubset[matchingPos][0]

		#Only consider interactions that are within the same TAD. We assume that interactions always take place within TADs, so otherwise this doesn't fit with our data.

		if matchingTad1[0] == matchingTad2[0] and matchingTad1[1] == matchingTad2[1] and matchingTad1[2] == matchingTad2[2]:
			tadStr = matchingTad1[0] + "_" + str(matchingTad1[1]) + "_" + str(matchingTad1[2])
			if tadStr not in tadInteractions:
				tadInteractions[tadStr] = []
			tadInteractions[tadStr].append(region1Pos)
			tadInteractions[tadStr].append(region2Pos)
		lineCount += 1

#Write the Hi-c interaction indices to a file grouped by TADs
with open(outFile, 'w') as outF:
	for tad in tadInteractions:

		line = tad
		line += "\t"
		for interaction in tadInteractions[tad]:
			line += str(interaction)
			line += ","
		line += "\n"
		outF.write(line)

max_stars_repo_name: UMCUGenetics/svMIL

max_issues_repo_name: UMCUGenetics/svMIL

max_forks_repo_name: UMCUGenetics/svMIL


text:
from typing import Tuple, Dict

import random
import numpy as np

import torch
from torchvision import datasets, transforms

from sklearn.metrics.pairwise import cosine_distances
from matplotlib import pyplot as plt

try:
    from torch.utils.tensorboard import SummaryWriter
except ImportError:
    from tensorboardX import SummaryWriter

CIFAR10_ANNOTATION = {
    0: 'airplane',
    1: 'automobile',
    2: 'bird',
    3: 'cat',
    4: 'deer',
    5: 'dog',
    6: 'frog',
    7: 'horse',
    8: 'ship',
    9: 'truck'
}

def plot_cifar_image(image, label=""):
    plt.title(label)
    plt.imshow(image.permute(1, 2, 0).numpy())
    plt.show()

class AccumulateStats:

    def __enter__(self):
        pass

    def __exit__(self):
        pass

    def __call__(self):
        pass

class AverageMeter(object):
    """
    Computes and stores the average and current value
    """

    def __init__(self):
        self.val = 0
        self.avg = 0
        self.sum = 0
        self.count = 0

    def reset(self):
        self.val = 0
        self.avg = 0
        self.sum = 0
        self.count = 0

    def update(self, val, n=1):
        self.val = val
        self.sum += val * n
        self.count += n
        self.avg = self.sum / self.count

class MeterLogger:

    def __init__(self, meters: Tuple[str], writer: SummaryWriter):
        self.average_meters: Dict[str, AverageMeter] = {k: AverageMeter() for k in meters}
        self._writer = writer

    def update(self, name: str, val, n=1):
        self.average_meters[name].update(val, n)

    def reset(self):
        for meter in self.average_meters.values():
            meter.reset()

    def write(self, step, prefix):
        for name, meter in self.average_meters.items():
            tag = prefix + '/' + name
            self._writer.add_scalar(tag, meter.avg, step)

class ImageLogger:

    def __init__(self, writer: SummaryWriter, mean=None, std=None):
        self._writer = writer
        self.mean = mean
        self.std = std

        if self.mean is not None:
            self.mean = torch.tensor(self.mean).reshape(1, 3, 1, 1)

        if self.std is not None:
            self.std = torch.tensor(self.std).reshape(1, 3, 1, 1)

    def write(self, images, reconstruction, step, prefix):
        images = images.cpu()
        reconstruction = reconstruction.cpu()

        if self.mean is not None and self.std is not None:
            images = images * self.std + self.mean
            reconstruction = reconstruction * self.std + self.mean

        image_tag = prefix + '/' + 'image'
        self._writer.add_images(image_tag, images, step)

        reconstruction_tag = prefix + '/' + 'reconstruction'
        self._writer.add_images(reconstruction_tag, reconstruction, step)

class VQEmbeddingLogger:

    def __init__(self, writer: SummaryWriter):
        self._writer = writer

    def write(self, embeddings, step):
        embeddings = embeddings.detach().cpu().numpy()
        sim = cosine_distances(embeddings)

        self._writer.add_image('cos_sim_vq_embeddings', sim, step, dataformats='HW')

def double_soft_orthogonality(weights: torch.Tensor):
    a = torch.norm(weights @ weights.t() - torch.eye(weights.shape[0]).to(weights.device)) ** 2
    b = torch.norm(weights.t() @ weights - torch.eye(weights.shape[1]).to(weights.device)) ** 2

    return a + b

def set_random_seed(seed: int, cuda: bool = False):
    random.seed(seed)
    np.random.seed(seed)
    torch.manual_seed(seed)

    if cuda:
        torch.cuda.manual_seed(seed)
        torch.cuda.manual_seed_all(seed)
        torch.cuda.random.manual_seed(seed)
        torch.backends.cudnn.deterministic = True


text:
import torch.nn as nn
import numpy as np
from abc import abstractmethod
from tensorboardX import SummaryWriter

class MetaTemplate(nn.Module):
  def __init__(self, model_func, n_way, n_support, flatten=True, leakyrelu=False, tf_path=None, change_way=True):
    super(MetaTemplate, self).__init__()
    self.n_way      = n_way
    self.n_support  = n_support
    self.n_query    = -1 #(change depends on input)
    self.feature    = model_func(flatten=flatten, leakyrelu=leakyrelu)
    self.feat_dim   = self.feature.final_feat_dim
    self.change_way = change_way  #some methods allow different_way classification during training and test
    self.tf_writer = SummaryWriter(log_dir=tf_path) if tf_path is not None else None

  @abstractmethod
  def set_forward(self,x,is_feature):
    pass

  @abstractmethod
  def set_forward_loss(self, x):
    pass

  def forward(self,x):
    out  = self.feature.forward(x)
    return out

  def parse_feature(self,x,is_feature):
    x = x.cuda()
    if is_feature:
      z_all = x
    else:
      x           = x.contiguous().view( self.n_way * (self.n_support + self.n_query), *x.size()[2:])
      z_all       = self.feature.forward(x)
      z_all       = z_all.view( self.n_way, self.n_support + self.n_query, -1)

    z_support   = z_all[:, :self.n_support]
    z_query     = z_all[:, self.n_support:]

    return z_support, z_query

  def correct(self, x):
    scores, loss = self.set_forward_loss(x)
    y_query = np.repeat(range( self.n_way ), self.n_query )

    topk_scores, topk_labels = scores.data.topk(1, 1, True, True)
    topk_ind = topk_labels.cpu().numpy()
    top1_correct = np.sum(topk_ind[:,0] == y_query)
    return float(top1_correct), len(y_query), loss.item()*len(y_query)

  def train_loop(self, epoch, train_loader, optimizer, total_it):
    print_freq = len(train_loader) // 10
    avg_loss=0
    for i, (x,_ ) in enumerate(train_loader):
      self.n_query = x.size(1) - self.n_support
      if self.change_way:
        self.n_way  = x.size(0)
      optimizer.zero_grad()
      _, loss = self.set_forward_loss(x)
      loss.backward()
      optimizer.step()
      avg_loss = avg_loss+loss.item()

      if (i + 1) % print_freq==0:
        print('Epoch {:d} | Batch {:d}/{:d} | Loss {:f}'.format(epoch, i + 1, len(train_loader), avg_loss/float(i+1)))
      if (total_it + 1) % 10 == 0 and self.tf_writer is not None:
        self.tf_writer.add_scalar(self.method + '/query_loss', loss.item(), total_it + 1)
      total_it += 1
    return total_it

  def test_loop(self, test_loader, record = None):
    loss = 0.
    count = 0
    acc_all = []

    iter_num = len(test_loader)
    for i, (x,_) in enumerate(test_loader):
      self.n_query = x.size(1) - self.n_support
      if self.change_way:
        self.n_way  = x.size(0)
      correct_this, count_this, loss_this = self.correct(x)
      acc_all.append(correct_this/ count_this*100  )
      loss += loss_this
      count += count_this

    acc_all  = np.asarray(acc_all)
    acc_mean = np.mean(acc_all)
    acc_std  = np.std(acc_all)
    print('--- %d Loss = %.6f ---' %(iter_num,  loss/count))
    print('--- %d Test Acc = %4.2f%% +- %4.2f%% ---' %(iter_num,  acc_mean, 1.96* acc_std/np.sqrt(iter_num)))

    return acc_mean


text:
#  ===============================================================================================================
#  Copyright (c) 2019, Cornell University. All rights reserved.
#
#  Redistribution and use in source and binary forms, with or without modification, are permitted provided that
#  the following conditions are met:
#
#      * Redistributions of source code must retain the above copyright otice, this list of conditions and
#        the following disclaimer.
#
#      * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and
#        the following disclaimer in the documentation and/or other materials provided with the distribution.
#
#      * Neither the name of Cornell University nor the names of its contributors may be used to endorse or
#        promote products derived from this software without specific prior written permission.
#
#  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED
#  WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
#  A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE
#  FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
#  TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
#  HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
#   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
#  OF SUCH DAMAGE.
#
#  Author: Kai Zhang (kz298@cornell.edu)
#
#  The research is based upon work supported by the Office of the Director of National Intelligence (ODNI),
#  Intelligence Advanced Research Projects Activity (IARPA), via DOI/IBC Contract Number D17PC00287.
#  The U.S. Government is authorized to reproduce and distribute copies of this work for Governmental purposes.
#  ===============================================================================================================

import numpy as np
from osgeo import gdal, gdal_array, osr
import os

def parse_proj_str(proj_str):
    idx1 = proj_str.find('UTM zone')
    idx2 = proj_str.find('",')

    sub_str = proj_str[idx1:idx2]
    hemisphere = sub_str[-1]
    zone_number = int(sub_str[-3:-1])

    return zone_number, hemisphere

def read_dsm_tif(file):
    assert (os.path.exists(file))

    ds = gdal.Open(file)
    geo = ds.GetGeoTransform()
    proj = ds.GetProjection()
    meta = ds.GetMetadata()
    width = ds.RasterXSize
    height = ds.RasterYSize

    assert (ds.RasterCount == 1)     # dsm is only one band
    band = ds.GetRasterBand(1)       # band index is one-based

    data_type = gdal_array.GDALTypeCodeToNumericTypeCode(band.DataType)
    assert (isinstance(data_type(), np.float32))
    image = np.zeros((ds.RasterYSize, ds.RasterXSize), dtype=np.float32)

    image[:, :] = band.ReadAsArray()
    nodata = band.GetNoDataValue()

    # to ease later processing, replace nodata regions with nan
    if nodata is not None:
        mask = np.isclose(image, nodata)
        image[mask] = np.nan

    zone_number, hemisphere = parse_proj_str(proj)
    # return a meta dict
    meta_dict = {
        'geo': geo,
        'proj': proj,
        'meta': meta,
        'img_width': width,
        'img_height': height,
        'nodata': nodata,
        'zone_number': zone_number,
        'hemisphere': hemisphere,
        'ul_easting': geo[0],
        'ul_northing': geo[3],
        'east_resolution': geo[1],
        'north_resolution': abs(geo[5])
    }

    # add other info for convenience
    # lr_easting are left-margin of the last column
    # lr_northing are top-margin of the last row
    meta_dict['lr_easting'] = meta_dict['ul_easting'] + (meta_dict['img_width']-1) * meta_dict['east_resolution']
    meta_dict['lr_northing'] = meta_dict['ul_northing'] - (meta_dict['img_height']-1) * meta_dict['north_resolution']
    meta_dict['area_width'] = meta_dict['lr_easting'] - meta_dict['ul_easting']
    meta_dict['area_height'] = meta_dict['ul_northing'] - meta_dict['lr_northing']
    meta_dict['alt_min'] = float(np.nanmin(image))  # for json serialization
    meta_dict['alt_max'] = float(np.nanmax(image))

    del ds
    return image, meta_dict

def get_driver(file):
    f_ext = os.path.splitext(file)[1]
    for i in range(gdal.GetDriverCount()):
        driver = gdal.GetDriver(i)
        if driver.GetMetadataItem(gdal.DCAP_RASTER):
            d_ext_str = driver.GetMetadataItem(gdal.DMD_EXTENSIONS)
            if d_ext_str is not None:
                for d_ext in d_ext_str.split(' '):
                    if f_ext == '.' + d_ext:
                        return driver
    return None

# out_file: .tif file to write
# geo: (ul_e, ul_n, e_resolution, n_resolution)
# utm_zone: (zone number, N or S)

def write_dsm_tif(image, out_file, geo, utm_zone, nodata_val=None):
    assert (len(image.shape) == 2)      # image should only be 2D

    ul_e, ul_n, e_resolution, n_resolution = geo
    zone_number, hemisphere = utm_zone

    # replace nan with no_data
    if nodata_val is not None:
        image = image.copy()    # avoid modify source data
        image[np.isnan(image)] = nodata_val
    else:
        nodata_val = np.nan

    driver = get_driver(out_file)
    out = driver.Create(out_file, image.shape[1], image.shape[0], 1,
                        gdal_array.NumericTypeCodeToGDALTypeCode(np.float32))
    band = out.GetRasterBand(1)     # one-based index
    band.WriteArray(image.astype(np.float32), 0, 0)
    band.SetNoDataValue(nodata_val)
    band.FlushCache()

    # syntax for geotransform
    # geotransform[0] = top left x
    # geotransform[1] = w-e pixel resolution
    # geotransform[2] = 0
    # geotransform[3] = top left y
    # geotransform[4] = 0
    # geotransform[5] = n-s pixel resolution (negative value)
    out.SetGeoTransform((ul_e, e_resolution, 0, ul_n, 0, -n_resolution))

    srs = osr.SpatialReference();
    srs.SetProjCS('WGS84 / UTM zone {}{}'.format(zone_number, hemisphere));
    srs.SetWellKnownGeogCS('WGS84');
    srs.SetUTM(zone_number, hemisphere=='N');
    out.SetProjection(srs.ExportToWkt())
    out.SetMetadata({'AREA_OR_POINT': 'Area'})

    del out

if __name__ == "__main__":
    gt_tif = '/bigdata/kz298/dataset/core3d_phase1b/Unclassified_Ground_Truth_D4-D9/AOI-D4-Jacksonville/AOI-D4-DSM.tif'
    image, meta_dict = read_dsm_tif(gt_tif)
    print(meta_dict['geo'])
    print(type(meta_dict['geo']))