text:
import numpy as np
from multiprocessing import Pool, cpu_count
import statsmodels.api as sm
from scipy.stats import norm
from statsmodels.gam.api import GLMGam, BSplines
from tqdm import tqdm
from itertools import product
import pandas as pd
import torch
from statsmodels.stats.proportion import proportion_confint
import rpy2.robjects.packages as packages

# Import rpy for testing with R functions
import rpy2.robjects as robj
CondIndTests = packages.importr("CondIndTests")
GeneralisedCovarianceMeasure = packages.importr("GeneralisedCovarianceMeasure")
dHSIC = packages.importr("dHSIC")

# Load files from parent folders
import os
import sys
try: sys.path.append(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
except NameError: print("Cannot load testing module")
from wrapper_resampler import ShiftedTester

# Help functions
alpha = 3
def e(n=1, mu=0): return np.random.normal(size=(n, 1), loc=mu)  # Gaussian noise
def coin(n, p=0.5): return 1*(np.random.uniform(size=(n,1)) < 0.5)
def cb(*args): return np.concatenate(args, axis=1)  # Col bind
def to_r(x): return robj.FloatVector(x)

np.random.seed(1)

# Simulate data from a Gaussian SCM
def scm(n, causal_effect=1, causal_exponent=2):
    X1 = e(n, mu=4*coin(n)-2) # Bimodal Gaussian distribution
    X2 = -X1**2 + 2*e(n)
    X3 = np.sin(X2) + causal_effect*X1**causal_exponent + 2*e(n)
    return cb(X1, X2, X3)

p = norm.pdf

# Weight
def weight(X):
    return p(X[:,1], loc = X[:,1].mean(), scale=np.sqrt(X[:,1].var()))/p(X[:,1], loc=-(X[:,0])**2, scale=2)

# Fitted weight
def weight_fit(X):
    mod1 = GLMGam(X[:,1], smoother=BSplines(X[:,0], df = 10, degree=3)).fit()
    return p(X[:,1], loc = X[:,1].mean(), scale=np.sqrt(X[:,1].var()))/p(X[:,1], loc=mod1.fittedvalues, scale=np.sqrt(mod1.scale))

# Test function: LinReg
def T(X):
    return (sm.OLS(X[:, 2], sm.tools.add_constant(X[:, 0])).fit().pvalues[1] < 0.05)*1

# HSIC test statistic
def THSIC(X):
    pval = dHSIC.dhsic_test(X=to_r(X[:,0]), Y=to_r(X[:,2])).rx2("p.value")[0]
    return 1*(pval < 0.05)

# Kernel Conditional Independence
def KCI(X):
    pval = (CondIndTests.KCI(to_r(X[:,0]), to_r(X[:,2]), to_r(X[:,1])).rx2('pvalue'))[0]
    return 1*(pval < 0.05)

# GeneralisedCovarianceMeasure
def GCM(X):
    pval = (GeneralisedCovarianceMeasure.gcm_test(X=to_r(X[:,0]), Z=to_r(X[:,1]), Y=to_r(X[:,2]), regr_method="gam").rx2('p.value'))[0]
    return 1*(pval < 0.05)

# Loop parameters
causal_effects = np.linspace(0, 1.5, num=7)
causal_exponents = [1, 2]
n = 150
tests = {"HSIC": THSIC, "LinReg": T, "HSICfit": THSIC, "GCM": GCM, "KCI": KCI}

# Compute for goodness-of-fit quantile
cutoff = np.quantile(np.random.uniform(size=(1000, 15)).min(axis=1), 0.005)

# Conditional distribution for goodness-of-fit
def get_cond(loc, scale):
    def cond(X, Y):
        slope = torch.tensor([0])
        linear = loc
        return -(Y.view(-1)-linear)**2/(2.0*scale**2)
    return cond

combinations = list(product(causal_effects, causal_exponents, tests.keys()))

## Wrap as function to apply multiprocessing
def conduct_experiment(i=None):
    out = []
    for c_e, c_exp, test in combinations:
        X = scm(n, causal_effect=c_e, causal_exponent=c_exp)
        if not test in ["GCM", "KCI"]:
            try:
                if test == "HSICfit":
                    psi = ShiftedTester(weight_fit, tests[test], rate=None, replacement="REPL-reject", reject_retries=500)
                else:
                    psi = ShiftedTester(weight, tests[test], rate=None, replacement="REPL-reject", reject_retries=500)
                m = psi.tune_m(X, j_x=[0], j_y=[1], gaussian=False,
                               cond=get_cond(loc=X[:,1].mean(), scale=np.sqrt(X[:,1].var())),
                               m_factor=1.3, p_cutoff=cutoff, replacement=False,
                               m_init=int(np.sqrt(n)), repeats=15)
                out.append(psi.test(X, m=m))
            except:
                # Catch errors from test statistic
                print(f"Error occurred {test}, {n}")
                out.append(np.nan)
        else:
            try:
                out.append(tests[test](X))
            except:
                print(f"Error occurred {test}, {n}")
                out.append(np.nan)
    return(out)

## Conduct multiple experiments with multiprocessing and export to R for plotting:
if __name__ == '__main__':
    repeats = 500

    # Multiprocess
    pool = Pool(cpu_count()-2)
    res = np.array(
        list(tqdm(pool.imap_unordered(conduct_experiment, range(repeats)), total=repeats)))
    pool.close()

    # Count non-nas, to be used for binomial confidence intervals
    counts = (~np.isnan(res)).sum(axis=0)
    res = np.nansum(res, axis=0)

    # Pack as data frame
    df = pd.DataFrame(
        [(x/c, *v, *proportion_confint(x, c, method="binom_test"), c) for x, v, c in zip(res, combinations, counts)],
        columns=["alpha", "CausalEffect", "EffectOrder", "Test", "Lower", "Upper", "Count"])

    # Export to R for ggplotting
    df['alpha'] = df["alpha"].replace(np.NaN, "NA")
    df.to_csv("experiment-cit.csv")

reason: import numpy,from scipy,import statsmodels,from statsmodels


text:
# __protected__ from numba import njit
# __protected__ @njit(cache=True, fastmath=True)

def rk2_step0(state_spect_n12, state_spect, tendencies_n, diss2, dt):
    # transonic block (
    #     complex128[][][] state_spect_n12, state_spect,
    #                      tendencies_n;
    #     float64[][] diss2;
    #     float dt
    # )
    state_spect_n12[:] = (state_spect + ((dt / 2) * tendencies_n)) * diss2

arguments_blocks = {
    "rk2_step0": ["state_spect_n12", "state_spect", "tendencies_n", "diss2", "dt"]
}
__transonic__ = ("0.4.7",)

ext: py

max_stars_repo_name: fluiddyn/transonic

max_issues_repo_name: fluiddyn/transonic

max_forks_repo_name: fluiddyn/transonic

reason: from numba


text:
import os
import pickle
import numpy as np

## Plotting modules
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()

def get_os_friendly_name(path):
    path = path.lstrip("/")
    path = path.replace("/", "_")
    return path

#############################################################################################
# def human_format(num):
#     magnitude = 0
#     while abs(num) >= 1000:
#         magnitude += 1
#         num /= 1000.0
#     # add more suffixes if you need them
#     return '%.2f%s' % (num, ['', 'K', 'M', 'G', 'T', 'P'][magnitude])

# print('the answer is %s' % human_format(7436313))

################################
# kilo  (k) 10^3   thousand    #
# mega  (M) 10^6   million     #
# giga  (G) 10^9   billion     #
# tera  (T) 10^12  trillion    #
# peta  (P) 10^15  quadrillion #
# exa   (E) 10^18  quintillion #
# zetta (Z) 10^21  sextillion  #
# yotta (Y) 10^24  septillion  #
################################
# 1 kilobyte (KiB)	1,024 bytes
# 1 megabyte (MiB)	1,048,576 bytes
# 1 gigabyte (GiB)	1,073,741,824 bytes
# 1 terabyte (TiB)	1,099,511,627,776 bytes
# 1 petabyte (PiB)	1,125,899,906,842,624 bytes
################################
# var cutoff = new SortedList<long, string> {
#        {59, "{3:S}" },
#        {60, "{2:M}" },
#        {60*60-1, "{2:M}, {3:S}"},
#        {60*60, "{1:H}"},
#        {24*60*60-1, "{1:H}, {2:M}"},
#        {24*60*60, "{0:D}"},
#        {Int64.MaxValue , "{0:D}, {1:H}"}
#      };
#############################################################################################

class BasePlotter:
    def __init__(self, loaders, output_dir=None, **options):
        self.loaders = loaders

        # Get output path
        if output_dir:
            self.output_dir = output_dir
        else:
            self.output_dir = self.loaders[0][0].getPlotsPath
        # Create output_dir:
        if not os.path.exists(self.output_dir):
            os.makedirs(self.output_dir)

        self.options = options

    def plot(self, keyx=None, keyy=None):
        pass

    def init_plot(self, nrows=1, ncols=1):
        ############################
        ### Overall plot options ###
        ############################
        context = self.options.get("context", "paper") # notebook | paper
        font_scale = self.options.get("font_scale", 2)
        line_width = self.options.get("line_width", 2.5)
        width  = self.options.get("width", 10)
        height = self.options.get("height", 8)

        sns.set_context(context, font_scale=font_scale, rc={"lines.linewidth": line_width})
        # sns.set_context("paper")
        # sns.set(font_scale=1.5)
        fig, ax = plt.subplots(nrows=nrows, ncols=ncols, figsize=(width, height))  # create figure & 1 axis
        return fig, ax

    def close_plot(self, fig, ax, path, name, save_pickle=True, bbox_extra_artists=[]):
        """This function saves the plot in `.pkl` and other formats.

        To later reload the `.pkl` file, one can use the following (to change format, labels, etc.):

        Example:

            import matplotlib.pyplot as plt
            %matplotlib notebook
            import pickle
            ax = pickle.load( open( "path_to_pkl_figure.pkl", "rb" ) )
            ax.set_title("Alternative title")
            plt.show()

        See:
            https://stackoverflow.com/a/12734723
        """

        ####################################
        ### Save figure to file as it is ###
        ####################################
        ## Make layout tight
        # plt.tight_layout(rect=[0, 0, 0.8, 1])
        plt.tight_layout()

        # Extension can be png/svg/pdf
        dpi = self.options.get("dpi", 300)
        ext = self.options.get("ext", 'svg')
        if isinstance(ext, list):
            filename = []
            for e in ext:
                filename += [os.path.join(path, name+"."+e)]
        else:
            filename = [os.path.join(path, name+"."+ext)]

        if bbox_extra_artists:
            for f in filename:
                fig.savefig(f, dpi=dpi, bbox_extra_artists=bbox_extra_artists, bbox_inches='tight')
        else:
            for f in filename:
                fig.savefig(f, dpi=dpi)
        print("Saved plot:", filename)

        # We do not use bbox_inches to make all figure sizes consistent.
        # fig.savefig(filename, bbox_inches='tight')
        # fig.savefig(filename, bbox_inches='tight', dpi=300)

        if save_pickle:
            # Save pickle files. They will be very easy to modify later.
            pkl_file = os.path.join(path, name+".pkl")
            pickle.dump(ax, open(pkl_file,'wb'))
            print("Saved file:", pkl_file)

        # Close the figure we have created.
        plt.close(fig)
        print()

# get_os_name(key)

reason: import numpy


text:
import numpy as np

from ms_deisotope.scoring import g_test_scaled
from ms_deisotope.averagine import glycan, PROTON, mass_charge_ratio
from ms_peak_picker.peak_set import FittedPeak
from brainpy import isotopic_variants

from glycan_profiling.chromatogram_tree import Unmodified

from .base import ScoringFeatureBase, epsilon

def envelope_to_peak_list(envelope):
    return [FittedPeak(e[0], e[1], 0, 0, 0, 0, 0, 0, 0) for e in envelope]

def scale_theoretical_isotopic_pattern(eid, tid):
    total = sum(p.intensity for p in eid)
    for p in tid:
        p.intensity *= total

def get_nearest_index(query_mz, peak_list):
    best_index = None
    best_error = float('inf')

    for i, peak in enumerate(peak_list):
        error = abs(peak.mz - query_mz)
        if error < best_error:
            best_error = error
            best_index = i
    return best_index

def align_peak_list(experimental, theoretical):
    retain = []
    for peak in experimental:
        retain.append(theoretical[get_nearest_index(peak.mz, theoretical)])
    return retain

def unspool_nodes(node):
    yield node
    for child in (node).children:
        for x in unspool_nodes(child):
            yield x

class IsotopicPatternConsistencyFitter(ScoringFeatureBase):
    feature_type = "isotopic_fit"

    def __init__(self, chromatogram, averagine=glycan, charge_carrier=PROTON, weighted=True):
        self.chromatogram = chromatogram
        self.averagine = averagine
        self.charge_carrier = charge_carrier
        self.scores = []
        self.intensity = []
        self.mean_fit = None
        self.weighted = weighted

        if chromatogram.composition is not None:
            if chromatogram.elemental_composition is not None:
                self.composition = chromatogram.elemental_composition
            else:
                raise Exception(chromatogram.composition)
        else:
            self.composition = None

        self.fit()

    def __repr__(self):
        return "IsotopicPatternConsistencyFitter(%s, %0.4f)" % (self.chromatogram, self.mean_fit)

    def generate_isotopic_pattern(self, charge, node_type=Unmodified):
        if self.composition is not None:
            tid = isotopic_variants(
                self.composition + node_type.composition,
                charge=charge, charge_carrier=self.charge_carrier)
            out = []
            total = 0.
            for p in tid:
                out.append(p)
                total += p.intensity
                if total >= 0.95:
                    break
            return out
        else:
            tid = self.averagine.isotopic_cluster(
                mass_charge_ratio(
                    self.chromatogram.neutral_mass + node_type.mass,
                    charge, charge_carrier=self.charge_carrier),
                charge,
                charge_carrier=self.charge_carrier)
            return tid

    def score_isotopic_pattern(self, deconvoluted_peak, node_type=Unmodified):
        eid, tid = self.prepare_isotopic_patterns(deconvoluted_peak, node_type)
        return g_test_scaled(None, eid, tid)

    def prepare_isotopic_patterns(self, deconvoluted_peak, node_type=Unmodified):
        tid = self.generate_isotopic_pattern(deconvoluted_peak.charge, node_type)
        eid = envelope_to_peak_list(deconvoluted_peak.envelope)
        scale_theoretical_isotopic_pattern(eid, tid)
        tid = align_peak_list(eid, tid)
        return eid, tid

    def fit(self):
        for node in self.chromatogram.nodes.unspool():
            for peak in node.members:
                score = self.score_isotopic_pattern(peak, node.node_type)
                self.scores.append(score)
                self.intensity.append(peak.intensity)
        self.intensity = np.array(self.intensity)
        self.scores = np.array(self.scores)
        if self.weighted:
            self.mean_fit = np.average(self.scores, weights=self.intensity)
        else:
            self.mean_fit = np.mean(self.scores)

    @classmethod
    def score(cls, chromatogram, *args, **kwargs):
        return max(1 - cls(chromatogram).mean_fit, epsilon)

class IsotopicPatternConsistencyModel(ScoringFeatureBase):
    feature_type = "isotopic_fit"

    def __init__(self, averagine=glycan, charge_carrier=PROTON):
        self.averagine = averagine
        self.charge_carrier = charge_carrier

    def fit(self, chromatogram):
        fitter = IsotopicPatternConsistencyFitter(
            chromatogram, averagine=self.averagine,
            charge_carrier=self.charge_carrier)
        return fitter

    def score(self, chromatogram, *args, **kwargs):
        return max(1 - self.fit(chromatogram).mean_fit, epsilon)

max_stars_repo_name: mstim/glycresoft

max_issues_repo_name: mstim/glycresoft

max_forks_repo_name: mstim/glycresoft


text:
import json
import time
import _thread
import numpy as np
from os.path import basename, join
from face_plus_plus import process
import os

"""
process image to get landmarks
"""

def png2json(image, json_path=None, suf_num=4):
    base_name = basename(image)
    json_name = base_name[:-suf_num] + '.json'
    json_file = json_name
    if json_path is not None:
        json_file = join(json_path, json_name)
    return json_file

def get_key():
    key = ['k1', 'k2', 'k3']
    skey = ['s1', 's2', 's3']
    return key, skey

def multi_process(id, images, key, secret, suf_num=4):
    finished = 0
    total = len(images)
    print('-- {} -- starting process, totally {} files !'.format(id, total))
    for image in images:
        json_file = png2json(image, json_path, suf_num=suf_num)
        json_data = process(image, key, secret)
        if json_data is not None:
            with open(json_file, 'w') as f:
                json.dump(json_data, f)
            finished = finished + 1
        time.sleep(1)
        if finished % 100 == 0 and finished > 0:
            print('-- {} -- [{}/{}] finish '.format(id, finished, total))

if __name__ == '__main__':
    keys, skeys = get_key()

    # TODO change the path to your local path
    img_path = 'path/to/image'
    json_path = 'path/to/save json'
    suf_num = 7  # 7 for multipie and 4 for lfw and other dataset, 7 means to remove _xx.png and 4 means to remove .png/.jpg
    while True:
        images = os.listdir(img_path)
        images = [join(img_path, img) for img in images]
        # deduplicate
        images = [img for img in images if not os.path.exists(png2json(img, json_path, suf_num=suf_num))]

        if len(images) == 0:
            break
        print('Starting process, totally {} files !'.format(len(images)))

        key_num = len(keys)
        skip = int(len(images) / key_num)
        index = np.arange(0, key_num + 1) * skip
        index[-1] = len(images)
        print(index)

        try:
            for i in range(key_num - 1):
                _thread.start_new_thread(multi_process, (i, images[index[i]:index[i + 1]], keys[i], skeys[i], suf_num))
        except Exception as e:
            print('Error: ', e)
        multi_process(key_num - 1, images[index[-2]:], keys[0], skeys[0], suf_num)

max_issues_repo_name: Yuxiang-Wei/FFWM

max_forks_repo_name: Yuxiang-Wei/FFWM


text:
"""
Test to compare normal computation methods
Author: Jeff Mahler
"""
import cv2
import numpy as np
import pcl
import os
import sys

import IPython

from autolab_core import RigidTransform
from visualization import Visualizer2D as vis2d
from visualization import Visualizer3D as vis3d

from perception import DepthImage, CameraIntrinsics

KSIZE = 9

if __name__ == '__main__':
    depth_im_filename = sys.argv[1]
    camera_intr_filename = sys.argv[2]

    camera_intr = CameraIntrinsics.load(camera_intr_filename)
    depth_im = DepthImage.open(depth_im_filename, frame=camera_intr.frame)

    depth_im = depth_im.inpaint()

    point_cloud_im = camera_intr.deproject_to_image(depth_im)
    normal_cloud_im = point_cloud_im.normal_cloud_im(ksize=KSIZE)

    vis3d.figure()
    vis3d.points(point_cloud_im.to_point_cloud(), scale=0.0025)

    alpha = 0.025
    subsample = 20
    for i in range(0, point_cloud_im.height, subsample):
        for j in range(0, point_cloud_im.width, subsample):
            p = point_cloud_im[i,j]
            n = normal_cloud_im[i,j]
            n2 = normal_cloud_im_s[i,j]
            if np.linalg.norm(n) > 0:
                points = np.array([p, p + alpha*n])
                vis3d.plot3d(points, tube_radius=0.001, color=(1,0,0))

                points = np.array([p, p + alpha*n2])
                vis3d.plot3d(points, tube_radius=0.001, color=(1,0,1))

    vis3d.show()

max_stars_repo_name: iamlab-cmu/perception

max_issues_repo_name: iamlab-cmu/perception

max_forks_repo_name: iamlab-cmu/perception


text:
import sys, os, re
import logging
import collections
from collections import defaultdict
import cluster_simulation_protos_pb2
import numpy

logging.basicConfig(level=logging.DEBUG)

def usage():
    print "usage: generate-neural-txt-from-dir.py <input_directory_name> <optional: base name for output files. (defaults to inputdir)>"
    sys.exit(1)

logging.debug("len(sys.argv): " + str(len(sys.argv)))

if len(sys.argv) < 2:
    logging.error("Not enough arguments provided.")
    usage()

try:
    input_dir_name = sys.argv[1]
    # Start optional args.
    if len(sys.argv) == 3:
        outfile_name_base = str(sys.argv[2])
    else:
        #make the output files the same as the input but add .txt to end
        outfile_name_base = os.path.join(input_dir_name, "neural.txt")

except:
    usage()

#input_dir_name = "/Users/damianfernandez/IdeaProjects/cluster-scheduler-simulator/experiment_results/prueba-4"
#outfile_name_base = os.path.join(input_dir_name, "neural.txt")

logging.info("Input directory: %s" % input_dir_name)

protobufs_to_iterate = {}

for filename in os.listdir(input_dir_name):
    if filename.endswith(".protobuf"):
        if "mesos" in filename:
            protobufs_to_iterate["mesos"] = os.path.join(input_dir_name, filename)
        elif "omega" in filename:
            protobufs_to_iterate["omega"] = os.path.join(input_dir_name, filename)
        elif "monolithic" in filename:
            protobufs_to_iterate["monolithic"] = os.path.join(input_dir_name, filename)
        continue
    else:
        continue

result_dict = collections.OrderedDict()

for scheduler, protobuf in protobufs_to_iterate.iteritems():

    # Read in the ExperimentResultSet.
    experiment_result_set = cluster_simulation_protos_pb2.ExperimentResultSet()
    infile = open(protobuf, "rb")
    experiment_result_set.ParseFromString(infile.read())
    infile.close()

    logging.debug("Processing %d experiment envs."
                  % len(experiment_result_set.experiment_env))
    for env in experiment_result_set.experiment_env:

        logging.debug("Handling experiment env (%s %s)."
                      % (env.cell_name, env.workload_split_type))
        logging.debug("Processing %d experiment results."
                      % len(env.experiment_result))
        prev_l_val = -1.0
        for exp_result in env.experiment_result:
            logging.debug("Handling experiment result with C = %f and L = %f."
                          % (exp_result.constant_think_time,
                             exp_result.per_task_think_time))
            params = {}
            params["resource_utilization"] = exp_result.cell_state_avg_cpu_utilization
            params["inter_arrival"] = exp_result.avg_job_interarrival_time
            logging.debug("Processing %d workload stats."
                          % len(exp_result.workload_stats))
            for workload_stat in exp_result.workload_stats:
                params[workload_stat.workload_name + "_" + "num_jobs_total"] = workload_stat.num_jobs
                params[workload_stat.workload_name + "_" + "num_jobs_scheduled"] = workload_stat.num_jobs_scheduled
                params[workload_stat.workload_name + "_" + "queue_first"] = workload_stat.avg_job_queue_times_till_first_scheduled
                params[workload_stat.workload_name + "_" + "queue_fully"] = workload_stat.avg_job_queue_times_till_fully_scheduled
                params[workload_stat.workload_name + "_" + "makespan"] = workload_stat.avg_makespan
                params[workload_stat.workload_name + "_" + "makespan_90p"] = workload_stat.makespan_90_percentile
                params[workload_stat.workload_name + "_" + "cpu_util"] = workload_stat.avg_cpu_util
                params[workload_stat.workload_name + "_" + "mem_util"] = workload_stat.avg_mem_util
                params[workload_stat.workload_name + "_" + "duration"] = workload_stat.avg_duration
                params[workload_stat.workload_name + "_" + "inter_arrival"] = workload_stat.avg_interarrival
                params[workload_stat.workload_name + "_" + "tasks"] = workload_stat.avg_tasks

            for param_name, value in params.iteritems():
                if param_name in result_dict:
                    if scheduler in result_dict[param_name]:
                        result_dict[param_name][scheduler].append(value)
                    else:
                        #dict1 = result_dict[param_name]
                        #dict1[scheduler] = [value]
                        result_dict[param_name][scheduler] = [value]
                else:
                    sched_dict = {scheduler: [value]}
                    result_dict[param_name] = sched_dict

#print(str(result_dict))
# Create output files.
# One output file for each unique (cell_name, scheduler_name, metric) tuple.
no_numpty_list = []
val_matrix = None
header_row = ""
body_rows = ""
for param_name, dict in result_dict.iteritems():
    for scheduler_name, values in dict.iteritems():
        header_row += param_name+"_"+scheduler_name+" "
        if val_matrix is None:
            val_matrix = numpy.empty([len(result_dict)*len(dict), len(values)])
        numpy.append(val_matrix, values)
        no_numpty_list.append(values)
header_row += "\n"
results_mat = map(list, zip(*no_numpty_list))
for r in results_mat:
    for c in r:
        body_rows += str(c) + " "
    body_rows += "\n"
# val_matrix = val_matrix.transpose()
#
# for x in xrange(val_matrix.shape[0]):
#     for y in xrange(val_matrix.shape[1]):
#         body_rows += str(val_matrix[x][y])+" "
#     body_rows += "\n"

logging.info("Creating output file: %s" % outfile_name_base)
outfile = open(outfile_name_base, "w")

outfile.write(header_row)
outfile.write(body_rows)
outfile.close()

max_stars_repo_name: DamianUS/sphere

max_issues_repo_name: DamianUS/sphere

max_forks_repo_name: DamianUS/sphere


text:
#! /usr/bin/env python3

# ============================== #
#   This script is employed to   #
#     compute (evolve) the Ks    #
#   parameters that are used in  #
#     the NN controller formula  #
# ============================== #

import os
import time
import random
import argparse

import subprocess
import numpy as np
import multiprocessing

import vie_measure # ViE minimization version with constraints

from deap import base, creator, tools

# ============ EVOLUTION CONFIGURATION ============ #

# length of the command input step
STEP_LEN = 10000
INIT_REC_0 = STEP_LEN*5
INIT_REC_180 = STEP_LEN*11 # skip init + 0 + transition values

# mutation configuration
MUT_PROB = 0.4

MIN_CONSTRAINS = [1.0]*5
MAX_CONSTRAINS = [2.0]*5

TARGET_COMPRESSION = 0.4 # stiffness: 0.1

SIM_SEED = 42
# ================================================= #

def evaluate_sol(solution):
    # prog <robot_conf> <K1> ... <K5> <noiseType> <noiseLevel>
    exec_string = '{} {} {} {} {}'.format('../apps/robotMeasurements',
                                          '../test_inputs/stiff_single/single_0.1.txt',
                                          str(solution),
                                          0,
                                          0,
                                          SIM_SEED)
    cproc = subprocess.run(exec_string.split(' '), capture_output=True)
    raw_data = cproc.stdout.decode("utf-8").strip().split('è,')[1:-1]

    data = np.asarray([[[float(d) for d in e.split(';')] for e in line.split(',')[1:]]
                       for line in raw_data])[:, 1:-1]

    # extract from the data only the initial compression values and the final ones
    initial_coms = np.median(data[list(range(INIT_REC_0-1, (INIT_REC_0+STEP_LEN*3)))], axis=0)    # command 0
    final_coms = np.median(data[list(range(INIT_REC_180-2, (INIT_REC_180+STEP_LEN*3)-2))], axis=0)  # command 180

    init_face_com_1 = np.mean(initial_coms[[1, 3, 4]], axis=0)
    init_face_com_5 = np.mean(initial_coms[[0, 2, 5]], axis=0)
    init_comp = np.linalg.norm(init_face_com_1 - init_face_com_5)

    fin_face_com_1 = np.mean(final_coms[[1, 3, 4]], axis=0)
    fin_face_com_5 = np.mean(final_coms[[0, 2, 5]], axis=0)
    fin_comp = np.linalg.norm(fin_face_com_1 - fin_face_com_5)

    # relative compression wrt to initial position
    current_compression = fin_comp/init_comp

    # consider invalid solutions the ones that compress the module more than the requirements
    return (current_compression, ) if current_compression >= TARGET_COMPRESSION else (1, )

def mutate_sol(solution):
    # gaussian mutation
    for i in range(len(solution.ks)):
        if random.random() < MUT_PROB:
            solution.ks[i] += np.random.normal(0, 0.12)

    return solution,

def cross(sol1, sol2):
    return sol1, sol2

class ModuleEvoConf:
    def __init__(self, min_constraints, max_constraints):
        self.min_constraints = min_constraints
        self.max_constraints = max_constraints
        self.ks = [min_constraints[i] + np.abs(np.random.normal(0, 0.4))
                     for i in range(len(min_constraints))]

    def __str__(self):
        return ' '.join(map(str, self.ks))

    def __repr__(self):
        return self.__str__()

    def __lt__(self, other):
        return self.fitness.values[0] < other.fitness.values[0]

    def __eq__(self, other):
        return (self.ks == other.ks
                and self.fitness.values[0] == other.fitness.values[0])

    def check_viability(self):
        for i in range(len(self.ks)):
            if self.ks[i] < self.min_constraints[i] or self.ks[i] > self.max_constraints[i]:
                return False

        return True

def gen_sol(create, min_constraints, max_constraints):
    """ Generate a random individual for the population. """
    return create(min_constraints, max_constraints)

creator.create('FitnessMin', base.Fitness, weights=(-1.0,))
creator.create('Configuration', ModuleEvoConf, fitness=creator.FitnessMin)

def main(tbox):
    # ======= SIMULATION CONFIGURATION ======= #
    max_num_sims = 12000
    pop_size = 32
    num_mutants = 32
    conv_rate = 0.05
    min_update = 0.001
    seeds = [42, 7240, 213, 96542, 10537996]

    # define output files and directories
    results_file = os.path.join('outputs', 'module_evo_results.csv')
    if not os.path.exists(os.path.dirname(results_file)):
        os.makedirs(os.path.dirname(results_file), exist_ok=True)
    # ======================================== #

    bests_history = []
    start_time = time.time()

    for seed in seeds:
        random.seed(seed)
        np.random.seed(seed)
        pop = tbox.population(pop_size)

        print('Start ViE execution...')
        [pop, bests] = vie_measure.run(pop, tbox, max_num_sims,
                                       conv_rate, num_mutants, min_update)
        print('ViE execution terminated!\n')

        bests_history.append(bests)

        for bh in bests_history:
            print('seed: {}'.format(seed))
            for b in sorted(bh):
                print('\t', b.fitness.values[0], '\t', b)

    with open(results_file, 'w') as out_file:
        out_file.write('fitness,K1,K4,K6,K7,K10\n')
        for bh in bests_history:
            for b in sorted(bh):
                out_file.write(str(b.fitness.values[0]) + ',' + ','.join(str(b).split(' ')) + '\n')

    end_time = int(time.time() - start_time)
    # record how much time it needed for current experiment
    elap_time = '{:02d}:{:02d}:{:02d}'.format(end_time // 3600,
                                              (end_time % 3600 // 60),
                                              end_time % 60)

    print('Computation time: {}'.format(elap_time))

if __name__ == '__main__':
    assert 0.0 <= TARGET_COMPRESSION <= 1.0, 'TARGET_COMPRESSION must be within the interval [0, 1].'

    parser = argparse.ArgumentParser(description='Script for evolving module parameters (Ks)')
    args = parser.parse_args()

    # prepare the structure of the EA algorithm
    toolbox = base.Toolbox()
    toolbox.register('conf', gen_sol, creator.Configuration, MIN_CONSTRAINS, MAX_CONSTRAINS)
    toolbox.register('population', tools.initRepeat, list, toolbox.conf)

    toolbox.register('evaluate', evaluate_sol)
    toolbox.register('mutate', mutate_sol)
    toolbox.register('mate', cross)

    toolbox.register('selectBest', tools.selBest)
    toolbox.register('select', tools.selRoulette)
    toolbox.register('fill', tools.selRandom)

    pool = multiprocessing.Pool()
    toolbox.register('map', pool.map)

    # execute the evolutionary simulation
    main(toolbox)
    pool.close()


text:
from __future__ import print_function, division, absolute_import, unicode_literals
__author__ = 'Alexandre Drouin'

import h5py as h
import json

from numpy import float

from ..spectrum import Spectrum
from ..utils import _is_mz_equal, _is_mz_precision_equal, _is_metadata_type_dict

def hdf5_save(file_name, spectra, compression_type="gzip", compression_opts=4):
    """
    Saves a list of spectra to a HDF5 file.

    Parameters:
    -----------
    file_name: str
        The path to the output file.

    spectra: list of Spectrum
        The list of spectra to save.

    compression_type: str, default="gzip"
        The type of compression to use for compressing the datasets. Refer to H5py's documentation.

    compression_opts: variable, default=4
        The compression options for the specified compression type. Refer to H5py's documentation. In the case of gzip
        compression, this specifies the level of compression used from 0 to 9.

    Note:
    -----
    * The spectra must have identical m/z values.
    * The spectra must have the same m/z precision.
    """
    if not _is_mz_equal(spectra[0].mz_values, spectra):
        raise ValueError("The spectra m/z values must be identicial to save to HDF5.")

    if not _is_mz_precision_equal(spectra[0].mz_precision, spectra):
        raise ValueError("The spectra m/z precisions must be equal to save to HDF5.")

    if not _is_metadata_type_dict(spectra):
        raise ValueError("The spectra metadata must be of type dict or None to save to HDF5.")

    file = h.File(file_name, "w")

    n_spectra = len(spectra)
    n_mz_values = len(spectra[0])

    file.create_dataset("precision", data=spectra[0].mz_precision)

    file.create_dataset("mz", data=spectra[0].mz_values, compression=compression_type,
                        compression_opts=compression_opts)

    spectra_intensity_dataset = file.create_dataset("intensity",
                                                    shape=(n_spectra, n_mz_values),
                                                    dtype=float,
                                                    chunks=(1, n_mz_values),
                                                    compression=compression_type,
                                                    compression_opts=compression_opts)

    for i, spectrum in enumerate(spectra):
        spectra_intensity_dataset[i] = spectrum.intensity_values

    dt = h.special_dtype(vlen=str)
    spectra_metadata_dataset = file.create_dataset("metadata", shape=(n_spectra,), dtype=dt)

    for i, spectrum in enumerate(spectra):
        spectra_metadata_dataset[i] = json.dumps(spectrum.metadata)

    file.close()

def hdf5_load(file_name, metadata=True):
    """
    Loads spectra from a HDF5 file.

    Parameters:
    -----------
    file_name: str
        The path to the file to load.

    metadata: boolean
        Defaults to True. Boolean to check if we load the metadata along with the spectrum data.

    Returns:
    -------
    spectra: list of Spectrum
        The list of spectra extracted form the file.
    """
    file = h.File(file_name, "r")

    mz_precision = file["precision"][...]
    mz_values = file["mz"][...]
    spectra_intensity_dataset = file["intensity"]

    if metadata and "metadata" in file:
        spectra_metadata_dataset = file['metadata']
    else:
        spectra_metadata_dataset = [None] * spectra_intensity_dataset.shape[0]

    spectra = []
    for spectrum_intensity_values, spectrum_metadata in zip(spectra_intensity_dataset, spectra_metadata_dataset):
        try:
            spectrum_metadata = json.loads(spectrum_metadata.decode("utf-8")) if not spectrum_metadata is None else None
        except AttributeError:
            spectrum_metadata = json.loads(spectrum_metadata) if not spectrum_metadata is None else None
        spectra.append(Spectrum(mz_values=mz_values, intensity_values=spectrum_intensity_values,
                                mz_precision=mz_precision, metadata=spectrum_metadata))
    file.close()

    return spectra

max_stars_repo_name: plpla/msvlm

max_issues_repo_name: plpla/msvlm

max_forks_repo_name: plpla/msvlm

reason: from numpy


text:
# IMPORT LIBRARIES
# =============================================================================
import seaborn as sns
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import warnings
# =============================================================================

# VERSION CHECK
# =============================================================================
pd.set_option('display.float_format', lambda x: '%.3f' % x)
sns.set(style='white', context='notebook', palette='deep')
warnings.filterwarnings('ignore')
sns.set_style('white')
# =============================================================================

# SET UP DATA
df = pd.read_csv("./hw04_data_set.csv")

# SPLIT DATA TRAIN AND TEST
# =============================================================================
num_train = 150
num_test = 122
train = df[0:num_train]
test = df[num_train:]
# =============================================================================

# SORTING BY ERUPTIONS
# =============================================================================
train = train.sort_values(by = ["eruptions"])
train = train.reset_index(drop=True)
test = test.sort_values(by = ["eruptions"])
test = test.reset_index(drop=True)
train_np_erup = np.array(train["eruptions"])
train_np_wait = np.array(train["waiting"])
test_np_erup = np.array(test["eruptions"])
test_np_wait = np.array(test["waiting"])
# =============================================================================

# CONFIGS
# =============================================================================
widht_param = 0.37
origin_param = 1.5
num_bins = 10
end_point = (origin_param) + (widht_param*num_bins)
bins = np.linspace(origin_param, end_point, num_bins+1)
# =============================================================================

# FUNCTIONS
# =============================================================================
# FUNCTION FOR EVALUATION
def evaluate(data,data1,data2,data3,data4):
    # ROOT MEAN SQUARE ERROR RMSE CALCULATION
    differences = kernel_smooth(data1,data2,data3,data4) - data["waiting"]
    differences_squared = differences ** 2
    mean_of_differences_squared = differences_squared.mean()
    rmse_val = np.sqrt(mean_of_differences_squared)
    rmse_val = float("{0:.4f}".format(rmse_val))
    # OUTPUT
    print("Kernel Smoother => RMSE is ", rmse_val, " when h is ", widht_param)

# GAUSSIAN KERNEL FUNCTION TO USE AS KERNEL METHOD
def gaussian_kernel(x, x0, h):
    return np.exp(- 0.5 * np.power((x - x0) / h, 2) )

# IMPLEMENTATION OF KERNEL SMOOTHING
def kernel_smooth(x_test, x_train, y_train, bandwidth, kernel_func=gaussian_kernel):
    return np.array([(kernel_func(x_train, x0, bandwidth).dot(y_train) ) /
                     kernel_func(x_train, x0, bandwidth).sum() for x0 in x_test])
# =============================================================================

if __name__ == '__main__':

    # PLOTTING
    # =============================================================================
    a = plt.scatter(train["eruptions"], train["waiting"], edgecolors='b')
    b = plt.scatter(test["eruptions"], test["waiting"], color='r')
    plt.plot(test["eruptions"], kernel_smooth(test_np_erup,train_np_erup,train_np_wait,widht_param), linewidth = 3, color = 'k')
    plt.legend((a,b),
       ('train', 'test'),
       scatterpoints=1,
       loc='upper left',
       ncol=3,
       fontsize=10)
    plt.xlabel('Eruption time (min)')
    plt.ylabel('Waiting time to next eruption (min)')
    plt.title('h = 0.37')
    plt.show()
    # =============================================================================

    # EVALUATION WITH TEST DATA
    # =============================================================================
    evaluate(test,test_np_erup,train_np_erup,train_np_wait,widht_param)
    # =============================================================================


text:
"""
Similarity functions
"""
import math
import numpy as np
from rtree import index

#pylint: disable=invalid-name

def dot(p1, p2, debug = False):
    """Dot product between two points

    Args:
        p1 ([float, float]): x and y coordinates
        p2 ([float, float]): x and y coordinates
    Returns:
        float
    """
    return p1[0] * p2[0] + p1[1] * p2[1]

def normalize(p, debug = False):
    """Normalizes a point/vector

    Args:
        p ([float, float]): x and y coordinates
    Returns:
        float
    """
    l = math.sqrt(p[0]**2 + p[1]**2)
    return [0.0, 0.0] if l == 0 else [p[0]/l, p[1]/l]

def angle(p1, p2, debug = False):
    """Angle between two points

    Args:
        p1 ([float, float]): x and y coordinates
        p2 ([float, float]): x and y coordinates
    Returns:
        float
    """
    return dot(p1, p2, debug)

def angle_similarity(l1, l2, debug = False):
    """Computes the similarity between two lines, based
    on their angles

    Args:
        l1 ([float, float]): x and y coordinates
        l2 ([float, float]): x and y coordinates
    Returns:
        float
    """
    return angle(l1, l2, debug)

def line(p1, p2, debug = False):
    """Creates a line from two points

    From http://stackoverflow.com/a/20679579

    Args:
        p1 ([float, float]): x and y coordinates
        p2 ([float, float]): x and y coordinates
    Returns:
        (float, float, float): x, y and _
    """
    A = (p1[1] - p2[1])
    B = (p2[0] - p1[0])
    C = (p1[0]*p2[1] - p2[0]*p1[1])
    return A, B, -C

def intersection(L1, L2, debug = False):
    """Intersects two line segments

    Args:
        L1 ([float, float]): x and y coordinates
        L2 ([float, float]): x and y coordinates
    Returns:
        bool: if they intersect
        (float, float): x and y of intersection, if they do
    """
    D = L1[0] * L2[1] - L1[1] * L2[0]
    Dx = L1[2] * L2[1] - L1[1] * L2[2]
    Dy = L1[0] * L2[2] - L1[2] * L2[0]
    if D != 0:
        x = Dx / D
        y = Dy / D
        return x, y
    else:
        return False

def distance(a, b, debug = False):
    """Distance between two points

    Args:
        a ([float, float]): x and y coordinates
        b ([float, float]): x and y coordinates
    Returns:
        float
    """
    return math.sqrt((b[0]-a[0])**2 + (b[1]-a[1])**2)

def distance_tt_point(a, b, debug = False):
    """ Euclidean distance between two (tracktotrip) points

    Args:
        a (:obj:`Point`)
        b (:obj:`Point`)
    Returns:
        float
    """
    return math.sqrt((b.lat-a.lat)**2 + (b.lon-a.lon)**2)

def closest_point(a, b, p, debug = False):
    """Finds closest point in a line segment

    Args:
        a ([float, float]): x and y coordinates. Line start
        b ([float, float]): x and y coordinates. Line end
        p ([float, float]): x and y coordinates. Point to find in the segment
    Returns:
        (float, float): x and y coordinates of the closest point
    """
    ap = [p[0]-a[0], p[1]-a[1]]
    ab = [b[0]-a[0], b[1]-a[1]]
    mag = float(ab[0]**2 + ab[1]**2)
    proj = dot(ap, ab, debug)
    if mag ==0 :
        dist = 0
    else:
        dist = proj / mag
    if dist < 0:
        return [a[0], a[1]]
    elif dist > 1:
        return [b[0], b[1]]
    else:
        return [a[0] + ab[0] * dist, a[1] + ab[1] * dist]

def distance_to_line(a, b, p, debug = False):
    """Closest distance between a line segment and a point

    Args:
        a ([float, float]): x and y coordinates. Line start
        b ([float, float]): x and y coordinates. Line end
        p ([float, float]): x and y coordinates. Point to compute the distance
    Returns:
        float
    """
    return distance(closest_point(a, b, p, debug), p, debug)

CLOSE_DISTANCE_THRESHOLD = 1.0
def distance_similarity(a, b, p, T=CLOSE_DISTANCE_THRESHOLD, debug = False):
    """Computes the distance similarity between a line segment
    and a point

    Args:
        a ([float, float]): x and y coordinates. Line start
        b ([float, float]): x and y coordinates. Line end
        p ([float, float]): x and y coordinates. Point to compute the distance
    Returns:
        float: between 0 and 1. Where 1 is very similar and 0 is completely different
    """
    d = distance_to_line(a, b, p, debug)
    r = (-1/float(T)) * abs(d) + 1

    return r if r > 0 else 0

def line_distance_similarity(p1a, p1b, p2a, p2b, T=CLOSE_DISTANCE_THRESHOLD, debug = False):
    """Line distance similarity between two line segments

    Args:
        p1a ([float, float]): x and y coordinates. Line A start
        p1b ([float, float]): x and y coordinates. Line A end
        p2a ([float, float]): x and y coordinates. Line B start
        p2b ([float, float]): x and y coordinates. Line B end
    Returns:
        float: between 0 and 1. Where 1 is very similar and 0 is completely different
    """
    d1 = distance_similarity(p1a, p1b, p2a, T=T, debug=debug)
    d2 = distance_similarity(p1a, p1b, p2b, T=T, debug=debug)
    return abs(d1 + d2) * 0.5

def line_similarity(p1a, p1b, p2a, p2b, T=CLOSE_DISTANCE_THRESHOLD, debug = False):
    """Similarity between two lines

    Args:
        p1a ([float, float]): x and y coordinates. Line A start
        p1b ([float, float]): x and y coordinates. Line A end
        p2a ([float, float]): x and y coordinates. Line B start
        p2b ([float, float]): x and y coordinates. Line B end
    Returns:
        float: between 0 and 1. Where 1 is very similar and 0 is completely different
    """
    d = line_distance_similarity(p1a, p1b, p2a, p2b, T=T, debug=debug)
    a = abs(angle_similarity(normalize(line(p1a, p1b, debug), debug), normalize(line(p2a, p2b, debug), debug), debug))
    return d * a

def bounding_box_from(points, i, i1, thr, debug = False):
    """Creates bounding box for a line segment

    Args:
        points (:obj:`list` of :obj:`Point`)
        i (int): Line segment start, index in points array
        i1 (int): Line segment end, index in points array
    Returns:
        Rect(float, float, float, float): with bounding box min x, min y, max x and max y
    """
    pi = points[i]
    pi1 = points[i1]

    min_lat = min(pi.lat, pi1.lat)
    min_lon = min(pi.lon, pi1.lon)
    max_lat = max(pi.lat, pi1.lat)
    max_lon = max(pi.lon, pi1.lon)

    return min_lat-thr, min_lon-thr, max_lat+thr, max_lon+thr

def segment_similarity(A, B, T=CLOSE_DISTANCE_THRESHOLD, debug = False):
    """Computes the similarity between two segments

    Args:
        A (:obj:`Segment`)
        B (:obj:`Segment`)
    Returns:
        float: between 0 and 1. Where 1 is very similar and 0 is completely different
    """
    l_a = len(A.points)
    l_b = len(B.points)

    idx = index.Index()
    dex = 0
    for i in range(l_a-1):
        idx.insert(dex, bounding_box_from(A.points, i, i+1, T), obj=[A.points[i], A.points[i+1]])
        dex = dex + 1

    prox_acc = []

    for i in range(l_b-1):
        ti = B.points[i].gen2arr()
        ti1 = B.points[i+1].gen2arr()
        bb = bounding_box_from(B.points, i, i+1, T, debug)
        intersects = idx.intersection(bb, objects=True)
        n_prox = []
        i_prox = 0
        a = 0
        for x in intersects:
            a = a + 1
            pi = x.object[0].gen2arr()
            pi1 = x.object[1].gen2arr()
            prox = line_similarity(ti, ti1, pi, pi1, T, debug)
            i_prox = i_prox + prox
            n_prox.append(prox)

        if a != 0:
            #prox_acc.append(i_prox / a)
            prox_acc.append(max(n_prox))
        else:
            prox_acc.append(0)

    return np.mean(prox_acc), prox_acc

def sort_segment_points(Aps, Bps, debug = False):
    """Takes two line segments and sorts all their points,
    so that they form a continuous path

    Args:
        Aps: Array of tracktotrip.Point
        Bps: Array of tracktotrip.Point
    Returns:
        Array with points ordered
    """
    mid = []
    j = 0
    mid.append(Aps[0])
    for i in range(len(Aps)-1):
        dist = distance_tt_point(Aps[i], Aps[i+1], debug)
        for m in range(j, len(Bps)):
            distm = distance_tt_point(Aps[i], Bps[m], debug)
            if dist > distm:
                direction = dot(normalize(line(Aps[i].gen2arr(), Aps[i+1].gen2arr(), debug), debug), normalize(Bps[m].gen2arr(), debug), debug)
                if direction > 0:
                    j = m + 1
                    mid.append(Bps[m])
                    break

        mid.append(Aps[i+1])
    for m in range(j, len(Bps)):
        mid.append(Bps[m])
    return mid


text:
"""This module defines plotting functionality for instantiated PulseTemplates using matplotlib.

Classes:
    - PlottingNotPossibleException.
Functions:
    - plot: Plot a pulse using matplotlib.
"""

from typing import Dict, Tuple, Any, Optional, Set, List, Union
from numbers import Real

import numpy as np
import warnings
import operator
import itertools

from qupulse._program import waveforms
from qupulse.utils.types import ChannelID, MeasurementWindow, has_type_interface
from qupulse.pulses.pulse_template import PulseTemplate
from qupulse.pulses.parameters import Parameter
from qupulse._program.waveforms import Waveform
from qupulse._program._loop import Loop, to_waveform

__all__ = ["render", "plot", "PlottingNotPossibleException"]

def render(program: Union[Loop],
           sample_rate: Real = 10.0,
           render_measurements: bool = False,
           time_slice: Tuple[Real, Real] = None,
           plot_channels: Optional[Set[ChannelID]] = None) -> Tuple[np.ndarray, Dict[ChannelID, np.ndarray],
                                                                    List[MeasurementWindow]]:
    """'Renders' a pulse program.

        Samples all contained waveforms into an array according to the control flow of the program.

        Args:
            program: The pulse (sub)program to render. Can be represented either by a Loop object or the more
                old-fashioned InstructionBlock.
            sample_rate: The sample rate in GHz.
            render_measurements: If True, the third return value is a list of measurement windows.
            time_slice: The time slice to be rendered. If None, the entire pulse will be shown.
            plot_channels: Only channels in this set are rendered. If None, all will.

        Returns:
            A tuple (times, values, measurements). times is a numpy.ndarray of dimensions sample_count where
            containing the time values. voltages is a dictionary of one numpy.ndarray of dimensions sample_count per
            defined channel containing corresponding sampled voltage values for that channel.
            measurements is a sequence of all measurements where each measurement is represented by a tuple
            (name, start_time, duration).
        """
    if has_type_interface(program, Loop):
        waveform, measurements = _render_loop(program, render_measurements=render_measurements)
    else:
        raise ValueError('Cannot render an object of type %r' % type(program), program)

    if waveform is None:
        return np.array([]), dict(), measurements

    if plot_channels is None:
        channels = waveform.defined_channels
    else:
        channels = waveform.defined_channels & plot_channels

    if time_slice is None:
        start_time, end_time = 0, waveform.duration

    elif time_slice[1] < time_slice[0] or time_slice[0] < 0 or time_slice[1] < 0:
        raise ValueError("time_slice is not valid.")

    else:
        start_time, end_time, *_ = time_slice

        # filter measurement windows
        measurements = [(name, begin, length)
                        for name, begin, length in measurements
                        if begin < end_time and begin + length > start_time]

    sample_count = (end_time - start_time) * sample_rate + 1
    if sample_count < 2:
        raise PlottingNotPossibleException(pulse=None,
                                           description='cannot render sequence with less than 2 data points')
    if not round(float(sample_count), 10).is_integer():
        warnings.warn("Sample count {sample_count} is not an integer. Will be rounded (this changes the sample rate).".format(sample_count=sample_count))

    times = np.linspace(float(start_time), float(end_time), num=int(sample_count), dtype=float)
    times[-1] = np.nextafter(times[-1], times[-2])

    voltages = {ch: waveforms._ALLOCATION_FUNCTION(times, **waveforms._ALLOCATION_FUNCTION_KWARGS)
                for ch in channels}
    for ch, ch_voltage in voltages.items():
        waveform.get_sampled(channel=ch, sample_times=times, output_array=ch_voltage)

    return times, voltages, measurements

def _render_loop(loop: Loop,
                 render_measurements: bool,) -> Tuple[Waveform, List[MeasurementWindow]]:
    """Transform program into single waveform and measurement windows.
    The specific implementation of render for Loop arguments."""
    waveform = to_waveform(loop)

    if render_measurements:
        measurement_dict = loop.get_measurement_windows()
        measurement_list = []
        for name, (begins, lengths) in measurement_dict.items():
            measurement_list.extend(zip(itertools.repeat(name), begins, lengths))
        measurements = sorted(measurement_list, key=operator.itemgetter(1))
    else:
        measurements = []

    return waveform, measurements

def plot(pulse: PulseTemplate,
         parameters: Dict[str, Parameter]=None,
         sample_rate: Real=10,
         axes: Any=None,
         show: bool=True,
         plot_channels: Optional[Set[ChannelID]]=None,
         plot_measurements: Optional[Set[str]]=None,
         stepped: bool=True,
         maximum_points: int=10**6,
         time_slice: Tuple[Real, Real]=None,
         **kwargs) -> Any:  # pragma: no cover
    """Plots a pulse using matplotlib.

    The given pulse template will first be turned into a pulse program (represented by a Loop object) with the provided
    parameters. The render() function is then invoked to obtain voltage samples over the entire duration of the pulse which
    are then plotted in a matplotlib figure.

    Args:
        pulse: The pulse to be plotted.
        parameters: An optional mapping of parameter names to Parameter
            objects.
        sample_rate: The rate with which the waveforms are sampled for the plot in
            samples per time unit. (default = 10)
        axes: matplotlib Axes object the pulse will be drawn into if provided
        show: If true, the figure will be shown
        plot_channels: If specified only channels from this set will be plotted. If omitted all channels will be.
        stepped: If true pyplot.step is used for plotting
        plot_measurements: If specified measurements in this set will be plotted. If omitted no measurements will be.
        maximum_points: If the sampled waveform is bigger, it is not plotted
        time_slice: The time slice to be plotted. If None, the entire pulse will be shown.
        kwargs: Forwarded to pyplot. Overwrites other settings.
    Returns:
        matplotlib.pyplot.Figure instance in which the pulse is rendered
    Raises:
        PlottingNotPossibleException if the sequencing is interrupted before it finishes, e.g.,
            because a parameter value could not be evaluated
        all Exceptions possibly raised during sequencing
    """
    from matplotlib import pyplot as plt

    channels = pulse.defined_channels

    if parameters is None:
        parameters = dict()

    program = pulse.create_program(parameters=parameters,
                                   channel_mapping={ch: ch for ch in channels},
                                   measurement_mapping={w: w for w in pulse.measurement_names})

    if program is not None:
        times, voltages, measurements = render(program,
                                               sample_rate,
                                               render_measurements=bool(plot_measurements),
                                               time_slice=time_slice)
    else:
        times, voltages, measurements = np.array([]), dict(), []

    duration = 0
    if times.size == 0:
        warnings.warn("Pulse to be plotted is empty!")
    elif times.size > maximum_points:
        # todo [2018-05-30]: since it results in an empty return value this should arguably be an exception, not just a warning
        warnings.warn("Sampled pulse of size {wf_len} is lager than {max_points}".format(wf_len=times.size,
                                                                                max_points=maximum_points))
        return None
    else:
        duration = times[-1]

    if time_slice is None:
        time_slice = (0, duration)

    legend_handles = []
    if axes is None:
        # plot to figure
        figure = plt.figure()
        axes = figure.add_subplot(111)

    if plot_channels is not None:
        voltages = {ch: voltage
                    for ch, voltage in voltages.items()
                    if ch in plot_channels}

    for ch_name, voltage in voltages.items():
        label = 'channel {}'.format(ch_name)
        if stepped:
            line, = axes.step(times, voltage, **{**dict(where='post', label=label), **kwargs})
        else:
            line, = axes.plot(times, voltage, **{**dict(label=label), **kwargs})
        legend_handles.append(line)

    if plot_measurements:
        measurement_dict = dict()
        for name, begin, length in measurements:
            if name in plot_measurements:
                measurement_dict.setdefault(name, []).append((begin, begin+length))

        color_map = plt.cm.get_cmap('plasma')
        meas_colors = {name: color_map(i/len(measurement_dict))
                       for i, name in enumerate(measurement_dict.keys())}
        for name, begin_end_list in measurement_dict.items():
            for begin, end in begin_end_list:
                poly = axes.axvspan(begin, end, alpha=0.2, label=name, edgecolor='black', facecolor=meas_colors[name])
            legend_handles.append(poly)

    axes.legend(handles=legend_handles)

    max_voltage = max((max(channel, default=0) for channel in voltages.values()), default=0)
    min_voltage = min((min(channel, default=0) for channel in voltages.values()), default=0)

    # add some margins in the presentation
    axes.set_xlim(-0.5+time_slice[0], time_slice[1] + 0.5)
    voltage_difference = max_voltage-min_voltage
    if voltage_difference>0:
        axes.set_ylim(min_voltage - 0.1*voltage_difference, max_voltage + 0.1*voltage_difference)
    axes.set_xlabel('Time (ns)')
    axes.set_ylabel('Voltage (a.u.)')

    if pulse.identifier:
        axes.set_title(pulse.identifier)

    if show:
        axes.get_figure().show()
    return axes.get_figure()

class PlottingNotPossibleException(Exception):
    """Indicates that plotting is not possible because the sequencing process did not translate
    the entire given PulseTemplate structure."""

    def __init__(self, pulse, description = None) -> None:
        super().__init__()
        self.pulse = pulse
        self.description = description
    def __str__(self) -> str:
        if self.description is None:
            return "Plotting is not possible. There are parameters which cannot be computed."
        else:
            return "Plotting is not possible: %s." % self.description

max_stars_repo_name: eendebakpt/qupulse

max_issues_repo_name: eendebakpt/qupulse

max_forks_repo_name: eendebakpt/qupulse


text:
import json

import gym
import numpy as np
import zmq
from typing import Any, Dict, Tuple
from gym import spaces

class WalkingRobotOpenLoopEnv(gym.Env):
    def __init__(self, mech_name: str, objective: str):
        context = zmq.Context()
        self.socket = context.socket(zmq.REQ)
        self.socket.connect("tcp://localhost:5572")
        self.objective = objective

        response = self._send_request({
            "type": "Initial",
            "mechName": mech_name,
        })

        self.id = response["id"]

    def step(self, action: dict):

        response = self._send_request({
            "id": self.id,
            "type": "Command",
            "objective": self.objective,
            "hyperParams": action,
        })

        return response, 0, True, None

    def reset(self):
        return None

    def render(self, mode="human"):
        ...

    def _send_request(self, request: Dict[str, Any]) -> Dict[str, Any]:
        request_message = json.dumps(request)
        self.socket.send(request_message.encode("UTF-8"))
        response = json.loads(self.socket.recv())
        return response

    def close(self):
        request = {
            "type": "Stop",
            "id": self.id,
        }
        try:
            self._send_request(request)
            self.socket.close()
        except zmq.error.ZMQError:
            pass

# Test the environment by doing 20 random steps in the game
if __name__ == "__main__":
    import time
    import gym
    import gym_space_engineers

    while True:
        env = gym.make("SpaceEngineers-WalkingRobot-OpenLoop-v0", mech_name="v2")
        env.reset()

        observation, _, _, _ = env.step({
            "duration": 50, "x": 5.250909035640978, "yRange": 6.6994406917744693, "yMin": -8.947822208617133,
            "zMin": -1.2499644968318842, "zRange": 2.005817017551913, "numberOfCycles": 20
        })

        print(observation)

        env.close()

max_stars_repo_name: GoodAI/SpaceEngineersGym

max_issues_repo_name: GoodAI/SpaceEngineersGym

max_forks_repo_name: GoodAI/SpaceEngineersGym


text:
import math
import statistics
import warnings

import numpy as np
from hmmlearn.hmm import GaussianHMM
from sklearn.model_selection import KFold
from asl_utils import combine_sequences

class ModelSelector(object):
    '''
    base class for model selection (strategy design pattern)
    '''

    def __init__(self, all_word_sequences: dict, all_word_Xlengths: dict, this_word: str,
                 n_constant=3,
                 min_n_components=2, max_n_components=10,
                 random_state=14, verbose=False):
        self.words = all_word_sequences
        self.hwords = all_word_Xlengths
        self.sequences = all_word_sequences[this_word]
        self.X, self.lengths = all_word_Xlengths[this_word]
        self.this_word = this_word
        self.n_constant = n_constant
        self.min_n_components = min_n_components
        self.max_n_components = max_n_components
        self.random_state = random_state
        self.verbose = verbose

    def select(self):
        raise NotImplementedError

    def base_model(self, num_states):
        # with warnings.catch_warnings():
        warnings.filterwarnings("ignore", category=DeprecationWarning)
        # warnings.filterwarnings("ignore", category=RuntimeWarning)
        try:
            hmm_model = GaussianHMM(n_components=num_states, covariance_type="diag", n_iter=1000,
                                    random_state=self.random_state, verbose=False).fit(self.X, self.lengths)
            if self.verbose:
                print("model created for {} with {} states".format(self.this_word, num_states))
            return hmm_model
        except:
            if self.verbose:
                print("failure on {} with {} states".format(self.this_word, num_states))
            return None

class SelectorConstant(ModelSelector):
    """ select the model with value self.n_constant

    """

    def select(self):
        """ select based on n_constant value

        :return: GaussianHMM object
        """
        best_num_components = self.n_constant
        return self.base_model(best_num_components)

class SelectorBIC(ModelSelector):
    """ select the model with the lowest Bayesian Information Criterion(BIC) score

    http://www2.imm.dtu.dk/courses/02433/doc/ch6_slides.pdf
    Bayesian information criteria: BIC = -2 * logL + p * logN
    """

    def select(self):
        """ select the best model for self.this_word based on
        BIC score for n between self.min_n_components and self.max_n_components

        :return: GaussianHMM object
        """
        warnings.filterwarnings("ignore", category=DeprecationWarning)

        tgtScore, tgtModel = float("inf"), None

        for components in range(self.min_n_components, self.max_n_components + 1):
            try:
                chModel = self.base_model(components)
                chModelLen = chModel.score(self.X, self.lengths)

                numCharct = self.X.shape[1]
                numVrs = components * (components - 1) + 2 * numCharct * components

                lgNum = np.log(self.X.shape[0])
                bicResult = -2 * chModelLen + numVrs * lgNum

                if bicResult < tgtScore:
                    tgtScore, tgtModel = bicResult, chModel

            except Exception as _:
                continue

        return tgtModel if tgtModel is not None else self.base_model(self.n_constant)

class SelectorDIC(ModelSelector):
    ''' select best model based on Discriminative Information Criterion

    Biem, Alain. "A model selection criterion for classification: Application to hmm topology optimization."
    Document Analysis and Recognition, 2003. Proceedings. Seventh International Conference on. IEEE, 2003.
    http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.58.6208&rep=rep1&type=pdf
    https://pdfs.semanticscholar.org/ed3d/7c4a5f607201f3848d4c02dd9ba17c791fc2.pdf
    DIC = log(P(X(i)) - 1/(M-1)SUM(log(P(X(all but i))
    '''

    def select(self):
        warnings.filterwarnings("ignore", category=DeprecationWarning)

        palabrsNvs = []
        modelos = []
        puntAcumul = []

        for palabra in self.words:
            if palabra != self.this_word:
                palabrsNvs.append(self.hwords[palabra])

        try:
            for stateNm in range(self.min_n_components, self.max_n_components + 1):
                markvModel = self.base_model(stateNm)

                topWrd = markvModel.score(self.X, self.lengths)
                modelos.append((topWrd, markvModel))

        except Exception as e:
            pass

        for indicePos, model in enumerate(modelos):
            topWrd, markvModel = model

            wordPos = [model[1].score(palabr[0], palabr[1]) for palabr in palabrsNvs]

            score_dic = topWrd - np.mean(wordPos)
            puntAcumul.append(tuple([score_dic, model[1]]))

            compMax = max(puntAcumul, key = lambda x: x[0])
        return compMax[1] if puntAcumul else None

class SelectorCV(ModelSelector):
    ''' select best model based on average log Likelihood of cross-validation folds

    '''

    def select(self):
        warnings.filterwarnings("ignore", category=DeprecationWarning)

        mediasArray = []
        foldSequences = KFold()

        try:
            for componentN in self.n_components:
                compBaseModel = self.base_model(componentN)
                foldValues = []

                for _, seqObs in foldSequences.split(self.sequences):
                    seqObsX, seqObsL = combine_sequences(seqObs, self.sequences)
                    foldValues.append(compBaseModel.score(seqObsX, seqObsL))

                mediasArray.append(np.mean(foldValues))

        except Exception as _:
            pass

        retStates = self.n_components[np.argmax(mediasArray)] if mediasArray else self.n_constant
        return self.base_model(retStates)

max_stars_repo_name: pablomateo/SignLanguageRecognizer

max_issues_repo_name: pablomateo/SignLanguageRecognizer

max_forks_repo_name: pablomateo/SignLanguageRecognizer


text:
'''
    This code is used to resample the firing ratemaps of neurons
'''

from __future__ import print_function, division

import numpy as np
try:
    from itertools import zip_longest as zip_longest
except:
    from itertools import izip_longest as zip_longest

# import external dictionaries
from relevance import *
from misc import *

def resample_ratemap(spike_data):
    small_bins = spike_data['binning_time']
    smaller_spike = spike_data['spikes']
    speed_filter = spike_data['speed_filter']
    masked_x_smaller = spike_data['pos_x']
    masked_y_smaller = spike_data['pos_y']
    hd_smaller = spike_data['theta']
    spatial_rate = spike_data['spatial_rate']
    x_edges = spike_data['x_edges']
    y_edges = spike_data['y_edges']
    HD_rate = spike_data['HD_rate']
    xmax, xmin = spike_data['x_limits']
    ymax, ymin = spike_data['y_limits']
    N_bins = spike_data['N_bins']
    th_N_bins = spike_data['HD_bins']
    n = spike_data['n']
    calc_option = spike_data['calc_option']

    unfiltered_spikes = np.copy(smaller_spike)
    smaller_spike[speed_filter] = 0

    Code = [(i,chr(ord('A')+i)) for i in np.arange(25)] + [(i+25,chr(ord('a')+i)) for i in np.arange(25)]

    # Get spatial trajectory codes
    x_edges = x_edges[:-1] + np.diff(x_edges)/2.
    y_edges = y_edges[:-1] + np.diff(y_edges)/2.
    x_edges, y_edges = np.meshgrid(x_edges, y_edges, indexing='ij')
    x_edges = x_edges.flatten()
    y_edges = y_edges.flatten()

    step = (xmax-xmin)/N_bins;
    x_alloc_codes = np.floor((x_edges-xmin)/step); x_codes = [Code[int(i)][1] for i in x_alloc_codes]
    y_alloc_codes = np.floor((y_edges-xmin)/step); y_codes = [Code[int(i)][1] for i in y_alloc_codes]
    binning_codes = np.array([x_codes[i]+y_codes[i] for i in np.arange(len(y_codes))]).astype("str")
    x_trajectory = (np.floor((x_smaller-xmin)/step)).astype("int"); x_trajectory_codes = [Code[int(i)][1] for i in x_trajectory]
    y_trajectory = (np.floor((y_smaller-xmin)/step)).astype("int"); y_trajectory_codes = [Code[int(i)][1] for i in y_trajectory]
    trajectory_codes = np.array([np.where(binning_codes==x_trajectory_codes[i]+y_trajectory_codes[i])[0] for i in np.arange(len(x_trajectory))])

    # Get head direction trajectory codes
    th_max = 2.*np.pi; th_min = 0; th_step = (th_max-th_min)/th_N_bins
    hd_trajectory = (np.floor((hd_smaller-th_min)/th_step)).astype("int")

    returned_data = []

    if calc_option[0]:
        r_trajectory = np.array([ratemap[trajectory_codes[i]]*small_bins for i in np.arange(len(trajectory_codes))])
        synthetic = np.array([np.random.binomial(1, min(r_trajectory[i],1.0), 100) for i in np.arange(len(r_trajectory))])
        synthetic[speed_filter] = 0; synthetic = synthetic.T

        resampled_spikes = []
        for neuron_number in np.arange(synthetic.shape[0]):
            s = [synthetic[neuron_number][n-i::n] for i in np.arange(n)]
            resampled_spikes.append(np.array([sum(r) for r in zip_longest(*s, fillvalue=0)]).astype("int"))
        resampled_spikes = np.array(resampled_spikes)

        resampled_relevance = np.zeros(resampled_spikes.shape[0])
        for i in np.arange(resampled_spikes.shape[0]):
            resampled_relevance[i] = parallelized_total_relevance((resampled_spikes.shape[1], resampled_spikes[i]))

        returned_data.append(resampled_relevance)

    if calc_option[1]:
        hd_trajectory = np.array([smoothened_firingrates[hd_trajectory[i]]*small_bins for i in np.arange(len(trajectory_codes))])
        synthetic_hd = np.array([np.random.binomial(1, min(hd_trajectory[i],1.0), 100) for i in np.arange(len(hd_trajectory))])
        synthetic_hd[speed_filter] = 0; synthetic_hd = synthetic_hd.T

        resampled_hd_spikes = []

        for neuron_number in np.arange(synthetic_hd.shape[0]):
            s = [synthetic_hd[neuron_number][n-i::n] for i in np.arange(n)]
            resampled_hd_spikes.append(np.array([sum(r) for r in zip_longest(*s, fillvalue=0)]).astype("int"))
        resampled_hd_spikes = np.array(resampled_hd_spikes)

        resampled_hd_relevance = np.zeros(resampled_hd_spikes.shape[0])
        for i in np.arange(resampled_spikes.shape[0]):
            resampled_hd_relevance[i] = parallelized_total_relevance((resampled_hd_spikes.shape[1], resampled_hd_spikes[i]))

        returned_data.append(resampled_hd_relevance)

    if calc_option[2]:
        rhd_trajectory = np.array([ratemap[trajectory_codes[i]]*smoothened_firingrates[hd_trajectory[i]]*small_bins for i in np.arange(len(trajectory_codes))])
        synthetic_rhd = np.array([np.random.binomial(1, min(rhd_trajectory[i],1.0), 100) for i in np.arange(len(rhd_trajectory))])
        synthetic_rhd[speed_filter] = 0; synthetic_rhd = synthetic_rhd.T

        resampled_rhd_spikes = []

        for neuron_number in np.arange(synthetic_rhd.shape[0]):
            s = [synthetic_rhd[neuron_number][n-i::n] for i in np.arange(n)]
            resampled_rhd_spikes.append(np.array([sum(r) for r in zip_longest(*s, fillvalue=0)]).astype("int"))
        resampled_rhd_spikes = np.array(resampled_rhd_spikes)

        resampled_rhd_relevance = np.zeros(resampled_hd_spikes.shape[0])
        for i in np.arange(resampled_spikes.shape[0]):
            resampled_rhd_relevance[i] = parallelized_total_relevance((resampled_rhd_spikes.shape[1], resampled_rhd_spikes[i]))

        returned_data.append(resampled_rhd_relevance)

    if len(returned_data)==1: return returned_data[0]
    else: return returned_data

max_stars_repo_name: rcubero/MSR

max_issues_repo_name: rcubero/MSR

max_forks_repo_name: rcubero/MSR


text:
import numpy as np
from lib.layers.linear import Linear
from lib.layers.sigmoid import Sigmoid
from lib.loss_functions.mse import MSELoss
from lib.models.sequential_model import SequentialModel
import matplotlib.pyplot as plt
from lib.parameter_updates import SGD, Momentum

rate = 0.09
mu = 0.8
data_size = 50000

model = SequentialModel(rate, MSELoss())

model.add_layer(Linear(1, 10, parameter_update=Momentum(mu), bias=True))
model.add_layer(Sigmoid())
model.add_layer(Linear(10, 1, parameter_update=Momentum(mu), bias=True))

data_x = np.array([np.random.rand(data_size)*180])
data_x = np.deg2rad(data_x)
data_y = np.sin(data_x)

errors = model.train(data_x, data_y, batch_size=5, error=False)

# print(errors)

test_x = np.arange(0, np.pi, 0.05).T
test_y = np.sin(test_x).T

pred = [model.forward(np.array([[x]])) for x in test_x]

plt.plot(test_x, test_y)
plt.scatter(test_x, pred)

plt.show()


text:
from enum import Enum
import numpy as np
from PIL import Image

from aoc.challenge_base import ChallengeBase

class MapTile(Enum):
    """
    Values for map tiles
    """
    SAND = 0
    CLAY = 1
    FLOWING_WATER = 2
    STILL_WATER = 3

class Challenge(ChallengeBase):
    """
    Day 16 challenges
    """
    def parse_input(self):
        """
        Parse input lines
        """
        # First, parse the clay ranges
        vert_ranges = []
        horz_ranges = []
        for line in self.lines:
            # Get the basic tokens
            first, second = line.strip().split(", ")
            first_axis, first_dist = first.split("=")
            second_axis, second_range = second.split("=")
            second_start, second_end = second_range.split("..")
            if first_axis == second_axis:
                raise RuntimeError("Expecting two different axes!")

            # Record as (first_axis_pos, range_start, range_end) tuple
            range_tuple = (int(first_dist), int(second_start), int(second_end))
            if first_axis == "x":
                vert_ranges.append(range_tuple)
            elif first_axis == "y":
                horz_ranges.append(range_tuple)
            else:
                raise RuntimeError("Expecting x or y axis")

        # Now use these to build our cross-section map. We want to clip to X and Y range given.
        self.min_x = min([r[0] for r in vert_ranges] + [r[1] for r in horz_ranges])
        self.min_y = min([r[0] for r in horz_ranges] + [r[1] for r in vert_ranges])

        # Build map as 2d array with 0 as sand
        # A real genius would RLE encode the whole thing for super efficiency, but for an idiot like me,
        # that seems like the water filling algorithms would become a bit complicated...
        max_x = max([r[0] for r in vert_ranges] + [r[2] for r in horz_ranges])
        max_y = max([r[0] for r in horz_ranges] + [r[2] for r in vert_ranges])

        # Could leave an extra top row as all sand, as we may need to flow water along it.
        # But we don't need to as the puzzle input kindly leaves some extra little nubbins
        # within and above the top container so it's not an issue in practice
        # However we do need a padding column in X before and after to allow waterfalls off the edge!
        self.map = np.zeros((max_y - self.min_y + 1, max_x - self.min_x + 1 + 2), dtype=np.uint8)

        # Draw the clay
        for x, y_start, y_end in vert_ranges:
            self.map[(y_start - self.min_y):(y_end - self.min_y + 1), x - self.min_x + 1] = MapTile.CLAY.value

        for y, x_start, x_end in horz_ranges:
            self.map[(y - self.min_y), (x_start - self.min_x + 1):(x_end - self.min_x + 2)] = MapTile.CLAY.value

        # Normalize the water source X location too, then set the first visible tile under it to water
        self.flowing_water = [(0, 500 - self.min_x + 1)]
        self.still_water = []
        self.map[self.flowing_water[0]] = MapTile.FLOWING_WATER.value

    def flow_sideways(self, coord, new_flowing_water):
        """
        Flow water sideways
        coord: the coordinate to attempt to fill (to side of existing flowing water)
        new_flowing_water: list to add any new filled coords to
        return: True if more flowing can happen, False if it's time to stop
        """
        if self.map[coord] == MapTile.CLAY.value:
            # Bounded on this side
            return False

        if self.map[coord] != MapTile.FLOWING_WATER.value:
            # Not already flowing
            self.map[coord] = MapTile.FLOWING_WATER.value
            new_flowing_water.append(coord)

        if self.map[coord[0] + 1, coord[1]] not in (MapTile.CLAY.value, MapTile.STILL_WATER.value):
            # Fallen off the end of this level
            return False

        return True

    def check_bounded_water(self, coord):
        """
        Check for a bounded water volume, starting from flowing water
        coord: the starting flowing water coord to look left and right from
        return: bounded, min_x, max_x where bounded is True if this can be converted to still
        """
        min_x = coord[1]
        max_x = coord[1]
        bounded_left = False
        for x in range(coord[1] - 1, -1, -1):
            tile_type = self.map[coord[0], x]
            if tile_type == MapTile.CLAY.value:
                bounded_left = True
                min_x = x + 1
                break
            elif tile_type not in (MapTile.FLOWING_WATER.value, MapTile.STILL_WATER.value):
                break

            # Also, if the tile underneath isn't clay or still water, we are not bounded
            if self.map[coord[0] + 1, x] not in (MapTile.CLAY.value, MapTile.STILL_WATER.value):
                break

        bounded_right = False
        for x in range(coord[1] + 1, self.map.shape[1]):
            tile_type = self.map[coord[0], x]
            if tile_type == MapTile.CLAY.value:
                bounded_right = True
                max_x = x - 1
                break
            elif tile_type not in (MapTile.FLOWING_WATER.value, MapTile.STILL_WATER.value):
                break

            # Also, if the tile underneath isn't clay or still water, we are not bounded
            if self.map[coord[0] + 1, x] not in (MapTile.CLAY.value, MapTile.STILL_WATER.value):
                break

        return (bounded_left and bounded_right), min_x, max_x

    def challenge1(self):
        """
        Day 16 challenge 1
        """
        self.parse_input()

        # Start the water flowing (as 2). When no new tiles are filled, we're done
        did_something = True
        try:
            while did_something:
                did_something = False

                # Iteratively propagate the water! First, continue any downward flow
                # This is simply those tiles with flowing water above that are also sand
                new_flowing_water = []
                for coord in self.flowing_water:
                    coord_below = (coord[0] + 1, coord[1])
                    if coord_below[0] < self.map.shape[0] and self.map[coord_below] == MapTile.SAND.value:
                        self.map[coord_below] = MapTile.FLOWING_WATER.value
                        new_flowing_water.append(coord_below)
                        did_something = True

                # Add any new flowing water
                self.flowing_water += new_flowing_water

                # Now do any sideways filling
                # This is simply those tiles beside flowing water above clay or water which are also above clay or water
                new_flowing_water = []
                fw_indices_to_remove = set()
                for fw_index, coord in enumerate(self.flowing_water):
                    coord_below = (coord[0] + 1, coord[1])
                    if coord_below[0] < self.map.shape[0] and self.map[coord_below] in (MapTile.CLAY.value, MapTile.STILL_WATER.value):
                        # Propagate left
                        current_new_flowing_water = []
                        for x in range(coord[1] - 1, -1, -1):
                            if not self.flow_sideways((coord[0], x), current_new_flowing_water):
                                break

                        # Propagate right
                        for x in range(coord[1] + 1, self.map.shape[1]):
                            if not self.flow_sideways((coord[0], x), current_new_flowing_water):
                                break

                        # Do we need to convert this line to still water?
                        is_bounded, min_x, max_x = self.check_bounded_water(coord)
                        if is_bounded:
                            # Convert the line to still water
                            for x in range(min_x, max_x + 1):
                                self.map[coord[0], x] = MapTile.STILL_WATER.value

                            self.still_water += current_new_flowing_water
                            fw_indices_to_remove.add(fw_index)
                            did_something = True
                        elif current_new_flowing_water:
                            # Add to our new flowing water list
                            new_flowing_water += current_new_flowing_water
                            did_something = True

                cropped_flowing_water = []
                for fw_index, coord in enumerate(self.flowing_water):
                    if fw_index not in fw_indices_to_remove:
                        if coord[0] < self.map.shape[0] - 1:
                            coord_below = (coord[0] + 1, coord[1])
                            coord_left = (coord[0], coord[1] - 1)
                            coord_right = (coord[0], coord[1] + 1)

                            # Remove flowing water from our "to-check" list if it has clay or still water below
                            # and has water either side. This is merely a speed optimization.
                            if (self.map[coord_below] in (MapTile.CLAY.value, MapTile.STILL_WATER.value)
                                and self.map[coord_left] == MapTile.FLOWING_WATER.value
                                and self.map[coord_right] == MapTile.FLOWING_WATER.value):
                                continue

                        cropped_flowing_water.append(coord)

                self.flowing_water = cropped_flowing_water

                # Add any new flowing water
                self.flowing_water += new_flowing_water
                print(f"Last flowing coord: {self.flowing_water[-1]}")

        except RuntimeError as ex:
            print(ex)

        # Debugging output
        im = Image.fromarray(self.map * 255 // 3, mode="L")
        im.save("out.png")

        # Count total water tiles (after first)
        water_tile_count = np.count_nonzero(self.map > 1)
        print(f"Water tile count: {water_tile_count}")

    def challenge2(self):
        """
        Day 16 challenge 2
        """
        # Count total water tiles (after first)
        water_tile_count = np.count_nonzero(self.map > 2)
        print(f"Still water tile count: {water_tile_count}")

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
REPORT_INTERVAL = 60
#CLEANUP_INTERVAL = 60
#CLEANUP_EXPIRED = 60
#DISK_CACHE_SIZE_ENV = "HTTPFS_DISK_CACHE_SIZE"
#DISK_CACHE_DIR_ENV = "HTTPFS_DISK_CACHE_DIR"

from errno import EIO, ENOENT
import time
import numpy
from fuse import FuseOSError, LoggingMixIn, Operations

if True:

    def read(self, path, size, offset, fh):
        # TODO rename size to read_size
        # TODO rename offset to read_offset

        if not path in self.patched_paths:
            # passthrough
            os.lseek(fh, offset, os.SEEK_SET)
            return os.read(fh, size)

        print(f"+ read    {path} + offset {offset} + size {size}")

        t1 = time.time()

        #self.logger.debug("read %s %s %s", path, offset, size)

        if t1 - self.last_report_time > REPORT_INTERVAL:
            """
            self.logger.info(
                "lru hits: {} lru misses: {} disk hits: {} total_requests: {}".format(
                    self.lru_hits,
                    self.lru_misses,
                    self.disk_hits,
                    self.disk_misses,
                    self.total_requests,
                )
            )
            """
            pass
        try:
            self.total_requests += 1

            # FIXME handle github paths. map path to blob url
            # TODO refactor copy-paste code: self.tree_data vs self.blob_data
            if not path in self.blob_sha:
                print(f"  read: path not in self.blob_sha")
                path_parts = path.split("/")
                # find parent in cache: bottom-up
                parent_path = None
                parent_depth = None
                for path_len in range((len(path_parts) - 1), 1, -1):
                    maybe_parent_path = "/".join(path_parts[0:path_len])
                    print(f"  read: path {path} -> parent? {parent_path}")
                    if maybe_parent_path in self.tree_data:
                        parent_path = maybe_parent_path
                        parent_depth = path_len
                        break
                if parent_path == None:
                    #print("  parent NOT found")
                    parent_path = '/'
                    parent_depth = 1
                print(f"  read: path {path} -> parent! {parent_path}")
                # fetch entries: top-down
                #print(f"fetch parents from {parent_depth + 1} to {len(path_parts)}")
                for path_len in range((parent_depth + 1), (len(path_parts) + 1)):
                    parent_path = "/".join(path_parts[0:path_len])

                    if not parent_path in self.tree_sha:
                        raise FuseOSError(ENOENT) # no such file
                        # TODO also in getattr

                    #print(f"  fetch path {parent_path} @ path_len {path_len}")
                    print(f"  read: dl_tree {parent_path} {self.tree_sha[parent_path]}")

                    url = f"https://api.github.com/repos/{self.owner}/{self.repo}/git/trees/{self.tree_sha[parent_path]}"
                    response = self.requests.get(url)
                    self.tree_data[parent_path] = response.json()
                    #print("  set tree_data %s" % parent_path)
                    verify_github_api.verify_tree(self, self.tree_sha[parent_path], self.tree_data[parent_path], parent_path)

                    # add subtree shas
                    # TODO refactor
                    path_prefix = '/' if parent_path == '/' else parent_path + '/'
                    #path_prefix = path + '/'
                    #if path == '/':
                    #    path_prefix = '/'
                    tree = self.tree_data[parent_path]['tree']
                    for tree_item in tree:
                        item_path = path_prefix + tree_item['path']
                        if tree_item['type'] == 'tree':
                            self.tree_sha[item_path] = tree_item['sha']
                        elif tree_item['type'] == 'blob':
                            self.blob_sha[item_path] = tree_item['sha']

            # now path should be in blob_sha
            if not path in self.blob_sha:
                print(f"  read: no such file")
                raise FuseOSError(ENOENT) # no such file

            sha = self.blob_sha[path]

            print(f"  read: path in self.blob_sha -> {sha}")

            #path_parts = path.split("/")
            #parent_path = "/".join(path_parts[:-1])

            attr = self.getattr(path) # TODO verify
            blob_size = attr['st_size']

            #url = "{}:/{}".format(self.schema, path[:-2])
            #url = f"https://api.github.com/repos/{self.owner}/{self.repo}/git/blobs/{sha}"
            #url = f"https://github.com/{self.owner}/{self.repo}/git/blobs/{sha}"
            url = f"https://raw.githubusercontent.com/{self.owner}/{self.repo}/{self.commit}{path}"

            """
            self.logger.debug("read url: {}".format(url))
            self.logger.debug(
                "offset: {} - {} request_size (KB): {:.2f} block: {}".format(
                    offset,
                    offset + size - 1,
                    size / 2 ** 10,
                    offset // self.chunk_size,
                )
            )
            """

            # size != blob_size
            output = numpy.zeros((size,), numpy.uint8)

            t1 = time.time()

            # nothing fetched yet
            last_fetched = -1
            curr_start = offset

            while last_fetched < offset + size:
                chunk_num = curr_start // self.chunk_size
                chunk_start = self.chunk_size * chunk_num

                chunk_id = (url, chunk_num)
                while chunk_id in self.getting:
                    # block following calls to read
                    time.sleep(0.05)

                self.getting.add(chunk_id)
                print(f"  read: call getchunk {url} {chunk_num}")
                chunk_data = self.getchunk(path, url, chunk_num)
                print(f"  read: getchunk -> {len(chunk_data)} bytes")
                self.getting.remove(chunk_id)

                """
                print("  first 128 bytes of chunk_data:")
                hexdump_canonical(chunk_data[:128])
                print("  last 128 bytes of chunk_data:")
                hexdump_canonical(chunk_data[-128:])
                """

                data_start = (
                    curr_start - chunk_num * self.chunk_size
                )

                data_end = min(self.chunk_size, offset + size - chunk_start)
                data = chunk_data[data_start:data_end]
                print(f"  read: data = {len(data)} bytes")

                """
                print("  first 128 bytes of data:")
                hexdump_canonical(data[:128])
                print("  last 128 bytes of data:")
                hexdump_canonical(data[-128:])
                """

                print(f"  chunk {chunk_num} from {chunk_start} to {chunk_start + self.chunk_size}")
                print(f"  data  {chunk_num} from {data_start} to {data_end}")

                d_start = curr_start - offset
                print(f"  read: output[{d_start} : {d_start + len(data)}] = data")
                output[d_start : d_start + len(data)] = data

                last_fetched = curr_start + (data_end - data_start)
                curr_start += data_end - data_start

            _bytes = bytes(output)

            """
            print(f"  read    {path} + offset {offset} + size {size}")
            print("  first 128 bytes of output:")
            hexdump_canonical(_bytes[:128])
            print("  last 128 bytes of output:")
            hexdump_canonical(_bytes[-128:])
            """

            """
            + read    /pkgs/top-level/all-packages.nix + offset 917504 + size 131072
            read: path in self.blob_sha -> daa8d9cd38564425c72e7f53fd72c5a6a5464348
            + getattr /pkgs/top-level/all-packages.nix -> mode 100644 size 1083238 [cache hit]
            + read    /pkgs/top-level/all-packages.nix + offset 1048576 + size 36864
            read: path in self.blob_sha -> daa8d9cd38564425c72e7f53fd72c5a6a5464348
            + getattr /pkgs/top-level/all-packages.nix -> mode 100644 size 1083238 [cache hit]
            read: call getchunk https://raw.githubusercontent.com/TLATER/nixpkgs/c1b3f029a39f39e621eb0f9ab4c18acb2e7f74d0/pkgs/top-level/all-packages.nix 1
            + getchunk /pkgs/top-level/all-packages.nix
            getchunk: cache miss @ lru
            getchunk: cache hit @ disk
            getchunk: cache hit @ disk -> 34662 bytes + type <class 'numpy.ndarray'>
            read: getchunk -> 34662 bytes
            read: data = 34662 bytes
            read: output[0 : 34662] = data
            first 128 bytes of output:
            hexdump: 128 bytes
            00000000  70 20 3d 20 63 61 6c 6c  50 61 63 6b 61 67 65 20  |p = callPackage |
            00000010  2e 2e 2f 74 6f 6f 6c 73  2f 70 61 63 6b 61 67 65  |../tools/package|
            00000020  2d 6d 61 6e 61 67 65 6d  65 6e 74 2f 6e 69 78 2d  |-management/nix-|

            output file: note the ^@ -> error byte
            'p = callPackage ' is start of chunk

            nix-template-rpm = callPackage ../build-support/templaterpm { inherit (pythonPackages) python toposort; };

            nix-t^@p = callPackage ../tools/package-management/nix-top { };

            nix-tree = haskell.lib.justStaticExecutables (haskellPackages.nix-tree);

            -> off by one error
            chunk 0 ends at 'nix-t'
            chunk 1 starts at 'p = callPackage ../tools/package-management/nix-top { };'

            ^@ is NULL byte

            """

            return _bytes

        except Exception as ex:
            self.logger.exception(ex)
            raise

max_stars_repo_name: milahu/python-fuse-githubfs

max_issues_repo_name: milahu/python-fuse-githubfs

max_forks_repo_name: milahu/python-fuse-githubfs


text:
import numpy
from PIL import Image, ImageFont, ImageDraw
from osr2mp4.ImageProcess import imageproc
from osr2mp4 import logger

def round_corner(radius, fill):
	"""Draw a round corner"""
	corner = Image.new('RGBA', (radius, radius), (0, 0, 0, 0))
	draw = ImageDraw.Draw(corner)
	draw.pieslice((0, 0, radius * 2, radius * 2), 180, 270, fill=fill)
	return corner

def round_rectangle(size, radius, fill):
	"""Draw a rounded rectangle"""
	width, height = size
	rectangle = Image.new('RGBA', size, fill)
	corner = round_corner(radius, fill)
	rectangle.paste(corner, (0, 0))
	rectangle.paste(corner.rotate(90), (0, height - radius))  # Rotate the corner and paste it
	rectangle.paste(corner.rotate(180), (width - radius, height - radius))
	rectangle.paste(corner.rotate(270), (width - radius, 0))
	return rectangle

def get_font(name, size=10):
	try:
		font = ImageFont.truetype(name, size=size)
	except Exception as e:
		font = ImageFont.load_default()
		logger.error(f"Error while loading beatmap font at {name}: {repr(e)}")
	return font

def prepare_pp_counter(settings, bg, beatmap):
	PP_CARD_WIDTH = 550
	PP_CARD_HEIGHT = 300
	PP_CARD_LOWER_HEIGHT = 75
	PP_CARD_SIZE = (PP_CARD_WIDTH, PP_CARD_HEIGHT)
	FONT_BIG = get_font(settings.ppsettings["Font"], size = 50)
	FONT_HEADER = get_font(settings.ppsettings["Font"], size = 25)
	FONT_DEFAULT_TITLE = get_font(settings.ppsettings["Font"], size = 20)

	try:
		bg = Image.open(bg).convert("RGBA")
		bg = bg.resize(PP_CARD_SIZE, Image.ANTIALIAS)
		imageproc.changealpha(bg, 0.8)
	except Exception as e:
		logger.error("Error while loading beatmap BG: " + repr(e))
		bg = Image.new("RGBA", PP_CARD_SIZE, (0, 0, 0, 205))

	# create mask
	mask = round_rectangle((PP_CARD_WIDTH*5, PP_CARD_HEIGHT*5), 150, (0,0,0,255))
	mask = mask.resize(PP_CARD_SIZE, Image.ANTIALIAS)

	# create counter in card
	counter = Image.new('RGBA', (PP_CARD_WIDTH,PP_CARD_LOWER_HEIGHT), "#222f3e")
	bg.paste(counter, box = (0,PP_CARD_HEIGHT - PP_CARD_LOWER_HEIGHT))
	text_draw = ImageDraw.Draw(bg)

	# title
	title_font = FONT_DEFAULT_TITLE
	beatmap_title = f'{beatmap.meta["Artist"]} - {beatmap.meta["Title"]}'
	for fs in range(15,50):
	    f = get_font(settings.ppsettings["Font"], size=fs)
	    title_size = text_draw.textsize(beatmap_title, f)
	    if (title_size[0] >= PP_CARD_WIDTH*0.85) or (title_size[1] >= PP_CARD_HEIGHT*0.12):
	        title_font = f
	        logger.info(f"PP Counter Title size: {fs}")
	        break
	text_draw.text((int((PP_CARD_WIDTH-title_size[0])/2),10), beatmap_title, "white", font=title_font)

	# pp
	text_draw.text((int(PP_CARD_WIDTH*0.25), int(PP_CARD_HEIGHT - PP_CARD_LOWER_HEIGHT/1.2)), "pp", "#c8d6e5", font=FONT_BIG)

	# hit counter
	text_draw.text((int(PP_CARD_WIDTH*0.47), int(PP_CARD_HEIGHT - PP_CARD_LOWER_HEIGHT/1.2)), "100", "#1dd1a1", font=FONT_HEADER)
	text_draw.text((int(PP_CARD_WIDTH*0.67), int(PP_CARD_HEIGHT - PP_CARD_LOWER_HEIGHT/1.2)), "50", "#54a0ff", font=FONT_HEADER)
	text_draw.text((int(PP_CARD_WIDTH*0.84), int(PP_CARD_HEIGHT - PP_CARD_LOWER_HEIGHT/1.2)), "miss", "#ff6b6b", font=FONT_HEADER)

	# round corners
	mask.paste(bg, mask = mask)
	mask.save(settings.temp + "pptemplate.png")

	return mask