text:
import random
from copy import deepcopy
from dataclasses import dataclass
from typing import List
from PIL import Image, ImageDraw, ImageFont
import numpy as np

WIDTH = 1100
SIZE = (0, WIDTH)
N = 6
K = 3

@dataclass
class Peak:
    x: int
    y: int

    def __eq__(self, other):
        return self.x == other.x and self.y == other.y

    def __hash__(self):
        return hash((self.x, self.y))

@dataclass
class Edge:
    first: Peak
    second: Peak
    weight: int = 0

    def __eq__(self, other):
        return (
            self.first == other.first
            and self.second == other.second
            or self.first == other.second
            and self.second == other.first
        )

@dataclass
class Graph:
    peaks: List[Peak]
    edges: List[Edge]

def generate_random_graph(peaks_count: int):
    radius = WIDTH / 2
    x_c = [
        SIZE[1] // 2 + radius * np.cos(2 * np.pi * j / peaks_count)
        for j in range(peaks_count)
    ]
    y_c = [
        SIZE[1] // 2 + radius * np.sin(2 * np.pi * j / peaks_count)
        for j in range(peaks_count)
    ]
    peaks = [Peak(int(x_c[i]), int(y_c[i])) for i in range(peaks_count)]
    edges = []
    for i in range(len(peaks)):
        for j in range(i, len(peaks)):
            if peaks[i] == peaks[j]:
                continue
            edges.append(Edge(peaks[i], peaks[j], random.randint(*SIZE)))
    return Graph(peaks, edges)

def find_min_edge_to_min_path(min_path, left_edges):
    left_edges.sort(key=lambda e: e.weight)
    if len(min_path) == 0:
        return left_edges[0]
    insulated_peaks = []
    for e in min_path:
        insulated_peaks.append(e.first)
        insulated_peaks.append(e.second)

    for i in range(len(left_edges)):
        if (
            left_edges[i].first in insulated_peaks
            or left_edges[i].second in insulated_peaks
        ):
            return left_edges[i]

def find_max_edge(edges: List[Edge]):
    return max(edges, key=lambda x: x.weight)

def knp(graph: Graph):
    min_path = []
    left_edges = deepcopy(graph.edges)
    while True:
        insulated_peaks = set()
        for e in min_path:
            insulated_peaks.add(e.first)
            insulated_peaks.add(e.second)
        if len(insulated_peaks) == len(graph.peaks):
            break
        edge = find_min_edge_to_min_path(min_path, left_edges)
        left_edges.remove(edge)
        min_path.append(edge)
    make_clusters(min_path)

def make_clusters(edges):
    edges.sort(key=lambda e: -1 * e.weight)
    im, d = get_peaks_image(edges)
    draw(edges[K - 1 :], im=im, d=d)

def get_peaks_image(edges: List[Edge]):
    im = Image.new("RGBA", (WIDTH, WIDTH), (255, 255, 255, 255))
    d = ImageDraw.Draw(im)
    for edge in edges:
        d.ellipse(
            (
                edge.first.x - 15,
                edge.first.y - 15,
                edge.first.x + 15,
                edge.first.y + 15,
            ),
            fill=(255, 0, 0),
            outline=(0, 0, 0),
        )
        d.ellipse(
            (
                edge.second.x - 15,
                edge.second.y - 15,
                edge.second.x + 15,
                edge.second.y + 15,
            ),
            fill=(255, 0, 0),
            outline=(0, 0, 0),
        )
    return im, d

def draw(edges: List[Edge], im=None, d=None):
    im = im or Image.new("RGBA", (WIDTH, WIDTH), (255, 255, 255, 255))
    d = d or ImageDraw.Draw(im)
    font = ImageFont.truetype("arial.ttf", size=20)
    for edge in edges:
        d.line(
            (edge.first.x, edge.first.y, edge.second.x, edge.second.y),
            fill=128,
            width=5,
        )
        d.text(
            (
                (edge.first.x + edge.second.x) // 2 + random.randint(-30, 30),
                (edge.first.y + edge.second.y) // 2 + random.randint(-30, 30),
            ),
            str(edge.weight),
            fill=(255, 0, 0),
            font=font,
        )
    for edge in edges:
        d.ellipse(
            (
                edge.first.x - 15,
                edge.first.y - 15,
                edge.first.x + 15,
                edge.first.y + 15,
            ),
            fill=(255, 0, 0),
            outline=(0, 0, 0),
        )
        d.ellipse(
            (
                edge.second.x - 15,
                edge.second.y - 15,
                edge.second.x + 15,
                edge.second.y + 15,
            ),
            fill=(255, 0, 0),
            outline=(0, 0, 0),
        )
    im.show()

if __name__ == "__main__":
    graph = generate_random_graph(N)
    print(graph.edges)
    draw(graph.edges)
    knp(graph)


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

import numpy as np
from scipy.integrate import quad
from scipy.optimize import curve_fit
from scipy.special import gamma

class zdist(object):
    """Parameterized redshift probability distribution."""
    def __init__(self, alpha, beta, z0, zmax=None):
        """Continuous model of a redshift probability distribution.

        The parameterized distribution is given by

            n(z) = A * [(z / z0)^alpha] * exp[-(z / z0)^beta],

        where A is the normalization factor.

        Parameters
        ----------
        alpha, beta, z0 : float
            Parameters of the distribution.
        zmax : float, optional
            Maximum redshift of the distribution. If given, the function
            returns zero for redshifts larger than `zmax`, and the
            normalization is computed by numerical integration. If None,
            np.inf is used, and the normalization is given by the gamma
            function. Default is None.

        Notes
        -----
        A good fit to the MICE simulated galaxy catalog are the parameters
        (alpha, beta, z0) = (0.88, 1.4, 0.78).

        Examples
        --------
        >>> nz = zdist(0.88, 1.4, 0.78)
        >>> nz.A
        2.0124103366443378
        >>> nz.pdf([-0.5, 0, 0.5, 1.0, 1.5])
        array([0.        , 0.        , 0.79567909, 0.60772375, 0.29427866])

        >>> nz = zdist(0.88, 1.4, 0.78, zmax=1.4)
        >>> nz.A
        2.4324379840968535
        >>> nz.pdf([-0.5, 0, 0.5, 1.0, 1.5])
        array([0.        , 0.        , 0.96175219, 0.73456706, 0.        ])

        """
        # Check inputs
        assert alpha > 0 and beta > 0 and z0 > 0, "Inputs must be positive."

        self.alpha = alpha
        self.beta = beta
        self.z0 = z0
        self.zmax = zmax
        self.A = self._normalize(self.zmax)

    def _normalize(self, zmax):
        """Compute the normalization factor A."""
        if zmax is None:
            # Use the gamma function
            norm = self.beta / self.z0 / gamma((1. + self.alpha) / self.beta)
            self.zmax = np.inf
        else:
            # Use numerical integration
            integrand = (lambda z: ((z / self.z0)**self.alpha) *
                         np.exp(-((z / self.z0)**self.beta)))
            norm = 1. / quad(integrand, 0, zmax)[0]

        return norm

    def pdf(self, z):
        """Compute the probability density at a given redshift.

        Parameters
        ----------
        z : array_like
            Redshift values.

        Returns
        -------
        float or numpy array
            Probability density at each `z`.

        """
        # Check input
        z = np.atleast_1d(z)
        z[(z < 0) | (z > self.zmax)] = 0

        # Compute density
        x = z / self.z0
        result = self.A * (x**self.alpha) * np.exp(-(x**self.beta))

        # Clean up
        if len(result) == 1:
            result = result[0]

        return result

    def cdf(self, z):
        """Compute the cumulative probability up to a given redshift.

        Parameters
        ----------
        z : array_like
            Redshift values.

        Returns
        -------
        float or numpy array
            Cumulative probability up to each `z`.

        """
        # Check input
        z = np.atleast_1d(z)

        # Compute cumulative probability
        result = np.array([quad(lambda x: self.pdf(x), 0, zz)[0] for zz in z])

        # Clean up
        if len(result) == 1:
            result = result[0]

        return result

    def fit(self, samples, zmax=None):
        """Fit the distribution to a set of samples.

        Parameters
        ----------
        samples : array_like
            Redshift samples.
        zmax : float, optional
            Maximum redshift cutoff of the distribution. Default is None.

        Notes
        -----
        The current parameters (alpha, beta, z0) will be replaced by those
        of the fit.

        """
        # Bin samples
        hist, bin_edges = np.histogram(samples, bins=512, density=True)
        zvals = (bin_edges[:-1] + bin_edges[1:]) / 2

        # Find best fit parameters
        fitfunc = (lambda z, A, alpha, beta, z0:
                   A * ((z / z0)**alpha) * np.exp(-((z / z0)**beta)))
        params, cov = curve_fit(fitfunc, zvals, hist, p0=[1., 0.75, 1., 0.75])

        # Set new parameters
        self.alpha = params[1]
        self.beta = params[2]
        self.z0 = params[3]
        self.zmax = zmax
        self.A = self._normalize(self.zmax)

ext: py

max_stars_repo_name: sfarrens/lenspack

max_issues_repo_name: sfarrens/lenspack

max_forks_repo_name: sfarrens/lenspack

reason: import numpy,from scipy


text:
"""Stratified sampling for probabilistic program analysis."""
import time
from typing import Optional

from absl import logging
import jax.numpy as np
import numpy as onp
from scipy import stats

from sympais import constraint
from sympais import logger as logger_lib
from sympais import profiles
from sympais import tasks
from sympais.realpaver_solver import RealPaverSolver

def _compute_box_prob_mass(profile, b):
  area = 1
  for name, int_ in b.intervals.items():
    target_dist = profile[name]
    area *= target_dist.cdf(int_.upper) - target_dist.cdf(int_.lower)
  return area

def truncnorm(loc, scale, lower, upper, num_samples, random_state=None):
  a = (lower - loc) / scale
  b = (upper - loc) / scale
  return stats.truncnorm(
      loc=loc, scale=scale, a=a, b=b).rvs(
          num_samples, random_state=random_state)

def run_stratified(task: tasks.Task,
                   seed: int,
                   num_samples: int,
                   batch_size: Optional[int] = None,
                   box_precision: float = 0.1,
                   logger: Optional[logger_lib.Logger] = None):
  """Run the stratified sampler for probabilistic program analysis.

  This resembles a simplified implementation qCoral that uses interval
  boxes produced by RealPaver to perform stratified sampling. This is equivalent
  to the sampling scheme used by qCoral for analyzing the sliced constraints for
  a given path.
  """

  profile = task.profile
  assert isinstance(profile, profiles.IndependentProfile)
  profile = profile.components
  domains = task.domains

  # Compute the mass of the truncated cdf
  assert isinstance(profile, dict)

  pcs = task.constraints
  boxes = RealPaverSolver().solve(pcs, precision=box_precision, domains=task.domains)
  constraint_fn = constraint.build_constraint_fn(pcs, domains)

  outer_boxes = boxes.outer_boxes
  inner_boxes = boxes.inner_boxes

  path_inner_area = sum([_compute_box_prob_mass(profile, b) for b in inner_boxes])
  outer_box_area = [_compute_box_prob_mass(profile, b) for b in outer_boxes]
  path_tot_area = sum([_compute_box_prob_mass(profile, b) for b in boxes])
  path_outer_area = path_tot_area - path_inner_area

  logging.info("inner vol %f, total vol %f", path_inner_area, path_tot_area)

  def _run_batch(rng, num_samples):
    estimate = path_inner_area
    var = 0.0
    # Only sample the outer boxes
    outer_boxes = boxes.outer_boxes
    # Compute the number of samples allocated for each box
    assert len(num_samples) == len(outer_boxes)
    for box_idx, b in enumerate(outer_boxes):
      area = outer_box_area[box_idx]
      num_samples_box = num_samples[box_idx]
      if num_samples_box == 0:
        continue
      proposals = _sample_box(rng, b, num_samples_box)
      box_mean = np.mean(constraint_fn(proposals))
      estimate += area * box_mean
      var += np.square(area) * (1 - box_mean) * (box_mean) / num_samples_box
    return {"mean": estimate, "var": var}

  def _sample_box(rng, box, num_samples):
    proposals = {}
    for k in profile:
      mean = profile[k].loc
      std = profile[k].scale
      if k not in box.intervals:
        # Some intervals may not have a particular variable,
        # which implies that the variable is not constrained
        lower = domains[k][0]
        upper = domains[k][1]
      else:
        interval = box.intervals[k]
        lower = interval.lower
        upper = interval.upper
      proposals[k] = truncnorm(
          loc=mean,
          scale=std,
          lower=lower,
          upper=upper,
          num_samples=num_samples,
          random_state=rng)
    return proposals

  def _allocate_samples(budget):
    outer_boxes = boxes.outer_boxes
    total_area = path_outer_area
    num_samples_boxes = [0] * len(outer_boxes)
    remainders = [0] * len(outer_boxes)
    allocated = 0

    for n in range(len(outer_boxes)):
      area = outer_box_area[n]
      proportion = area / total_area
      quota = budget * proportion
      num_samples_boxes[n] = int(onp.floor(quota))
      remainders[n] = quota - num_samples_boxes[n]
      allocated += num_samples_boxes[n]
    # Allocate remainders
    if allocated < budget:
      to_allocate = budget - allocated
      sort_index = onp.argsort(onp.array(remainders))[::-1]
      sort_index = sort_index[:to_allocate]
      for i in sort_index:
        num_samples_boxes[i] += 1
    return num_samples_boxes

  rng = onp.random.RandomState(seed)
  if batch_size is None:
    batch_size = num_samples
  num_batches = num_samples // batch_size
  assert num_samples % batch_size == 0
  num_samples_per_path = batch_size
  path_samples_sizes = _allocate_samples(num_samples_per_path)

  start_time = time.time()
  mean = 0.0
  var = 0.0
  n = 0
  for batch_idx in range(num_batches):
    tic = time.time()
    batch_result = _run_batch(rng, path_samples_sizes)
    batch_mean = batch_result["mean"]
    batch_var = batch_result["var"]
    mean = (batch_size * batch_idx * mean + batch_mean * batch_size) / (
        (batch_idx + 1) * batch_size)
    var = (
        batch_size / (n + batch_size) * batch_var + n / (n + batch_size) * var +
        (batch_size * n / ((batch_size + n)**2) * np.square(mean - batch_mean)))
    n += batch_size
    toc = time.time()
    logging.info(
        "mean %f, batch mean %f, num_samples %d, elapsed %.4f seconds",
        mean,
        batch_mean,
        (batch_idx + 1) * batch_size,
        toc - tic,
    )
    logs = {
        "mean": mean,
        "var": var,
        "cov": np.sqrt(var) / mean,
        "sample_count": (batch_idx + 1) * batch_size,
        "batch_mean": batch_mean,
        "batch_var": batch_var,
        "time": time.time() - start_time
    }
    if logger is not None:
      logger.write(logs)
  logging.info("Finished sampling, time %f seconds", time.time() - start_time)
  # if logger is not None:
  #   logger.write({"mean": mean})
  return {"mean": mean, "var": var}

max_stars_repo_name: ethanluoyc/sympais

max_issues_repo_name: ethanluoyc/sympais

max_forks_repo_name: ethanluoyc/sympais

reason: import numpy,from scipy,import jax


text:
from flask import g, current_app
from google.auth import credentials, default as default_creds
from google.cloud import bigtable

import numpy as np

# Hack the imports for now
from pychunkedgraph.backend import chunkedgraph

class DoNothingCreds(credentials.Credentials):
    def refresh(self, request):
        pass

def get_client(config):
    project_id = config.get('project_id', 'pychunkedgraph')

    if config.get('emulate', False):
        credentials = DoNothingCreds()
    else:
        credentials, project_id = default_creds()

    client = bigtable.Client(admin=True,
                             project=project_id,
                             credentials=credentials)
    return client

def get_cg():
    if 'cg' not in g:
        table_id = current_app.config['CHUNKGRAPH_TABLE_ID']
        client = get_client(current_app.config)
        g.cg = chunkedgraph.ChunkedGraph(table_id=table_id,
                                         client=client)
    return g.cg

def tobinary(ids):
    """ Transform id(s) to binary format

    :param ids: uint64 or list of uint64s
    :return: binary
    """
    return np.array(ids).tobytes()

max_stars_repo_name: fcollman/PyChunkedGraph

max_issues_repo_name: fcollman/PyChunkedGraph

max_forks_repo_name: fcollman/PyChunkedGraph

reason: import numpy


text:
# -*- coding: UTF-8 -*-
# !/usr/bin/python
# @time     :2019/4/4 10:00
# @author   :Mo
# @function :chatbot based search, encode sentence_vec by char

from conf.path_config import w2v_model_char_path
from conf.path_config import matrix_ques_part_path_char
from utils.text_tools import txtRead, txtWrite, getChinese
from conf.path_config import projectdir, chicken_and_gossip_path
from numpy import float32 as numpy_type
from collections import Counter
import pickle, jieba, os, re
import jieba.posseg as pseg
from gensim import matutils
from math import log
import numpy as np
import gensim
import jieba

def load_word2vec_model(path, bin=False, limit=None):
    word2vec_model = gensim.models.KeyedVectors.load_word2vec_format(path, limit=limit, binary=bin, unicode_errors='ignore')
    return word2vec_model

def encoding_question(w2v_model, char_list):
    '''    生成句子向量
    :param wordlist: 分词list
    :param is_replaced: 是否替换default true
    :param debug_mode: default false
    :return: array句子的向量 len=300
    '''
    try:
        sentence_vec = w2v_model.wv[word2vec_model.index2word[1]] * 0
    except:
        sentence_vec = w2v_model.wv[word2vec_model.index2word[0]] * 0

    for k in range(len(char_list)):
        char_list_one = char_list[k]
        if type(char_list_one) == str:
            try:
                sentence_vec = sentence_vec + w2v_model.wv[char_list_one]
            except Exception as e:
                print(str(e))
                if char_list_one not in [' ', '']:
                    sentence_vec = sentence_vec + 1
    return sentence_vec

def most_similar_sentence_vec(vec_ques, matrix_org, top_vec=20):
    """
      最相似的句子，句向量与矩阵点乘
    :param vec:
    :param matrix:
    :param keys:
    :param topn:
    :return:
    """
    # 首先对句向量矩阵标号
    matrix_org_index = list(range(len(matrix_org)))
    # Scale a vector to unit length. The only exception is the zerovector, which is returned back unchanged.
    vec_ques_mean = matutils.unitvec(np.array([vec_ques]).mean(axis=0)).astype(numpy_type)
    # matrix_org单位化
    matrix_org_norm = (matrix_org / np.sqrt((matrix_org ** 2).sum(-1))[..., np.newaxis]).astype(numpy_type)
    # 计算两个向量之间的相似度，使用numpy的dot函数，矩阵点乘
    matrix_vec_dot = np.dot(matrix_org_norm, vec_ques_mean)
    # 防止top_vec越界
    top_vec = min(len(matrix_org), top_vec)
    # 相似度排序
    most_similar_sentence_vec_sort = matutils.argsort(matrix_vec_dot, topn=top_vec, reverse=True)

    index_score = []
    for t in most_similar_sentence_vec_sort[:top_vec]:
        index_score.append([matrix_org_index[t], float(matrix_vec_dot[t])])
    return index_score

def create_matrix_org_np(sen_count, word2vec_model, qa_path, matrix_ques_path):
    """
      创建问题句向量
    :param sen_count: int
    :param word2vec_model: gensim model
    :param qa_path: str
    :param matrix_ques_path:str
    :return: None
    """
    if os.path.exists(matrix_ques_path):
        file_matrix_ques = open(matrix_ques_path, 'rb')
        matrix_ques = pickle.load(file_matrix_ques)
        return matrix_ques
    print('create_matrix_org_pkl start!')
    qa_dail = txtRead(qa_path, encodeType='utf-8')
    # questions = []
    matrix_ques = []
    count = 0
    for qa_dail_one in qa_dail:
        ques = getChinese(qa_dail_one.split('\t')[0])
        char_list = [ques_char for ques_char in ques]
        sentence_vec = encoding_question(word2vec_model, char_list)
        matrix_ques.append(sentence_vec)
        if len(matrix_ques)%sen_count == 0 and len(matrix_ques) != 0:
            print("count: " + str(count))
            count += 1
            np.savetxt(projectdir + "/Data/sentence_vec_encode_char/" + str(count)+".txt", matrix_ques)
            matrix_ques = []
            break

    # count += 1
    # np.savetxt(projectdir + "/Data/sentence_vec_encode_char/" + str(count)+".txt", matrix_ques)

    print('create_matrix_org_pkl ok!')
    # return matrix_ques

if __name__ == '__main__':

    # 读取问答语料
    syn_qa_dails = txtRead(chicken_and_gossip_path, encodeType='utf-8')
    # 读取词向量
    word2vec_model = load_word2vec_model(w2v_model_char_path, limit=None)
    # 创建标准问答中问题的句向量，存起来，到matrix_ques_path， 10万条，可自己设置，这里需要耗费点时间
    if not os.path.exists(matrix_ques_part_path_char):
        # matrix_ques = create_matrix_org_np(sen_count=100000, word2vec_model=word2vec_model, qa_path=chicken_and_gossip_path, matrix_ques_path=matrix_ques_part_path_char)
        create_matrix_org_np(sen_count=100000, word2vec_model=word2vec_model, qa_path=chicken_and_gossip_path, matrix_ques_path=matrix_ques_part_path_char)
    # 重载
    matrix_ques = np.loadtxt(matrix_ques_part_path_char)
    print("np.loadtxt(matrix_ques_part_path_char) ok!")
    while True:
        print("你问: ")
        ques_ask = input()
        ques_clean = getChinese(ques_ask)
        char_list = [ques_char for ques_char in ques_clean]
        sentence_vic = encoding_question(word2vec_model, char_list)
        top_20_qid = most_similar_sentence_vec(sentence_vic, matrix_ques, top_vec=20)
        try:
            print("小姜机器人: " + syn_qa_dails[top_20_qid[0][0]].strip().split("\t")[1])
            print([(syn_qa_dails[top_20_qid[i][0]].strip().split("\t")[0], syn_qa_dails[top_20_qid[i][0]].strip().split("\t")[1]) for i in range(len(top_20_qid))])
        except Exception as e:
            # 有的字符可能打不出来
            print(str(e))

reason: import numpy,from numpy


text:
import importlib
import logging
import os

import numpy as np
import torch
from torch.utils.data import Dataset

from args import args

def create_dir(dir_path):
    if not os.path.isdir(dir_path):
        os.makedirs(dir_path)

class ForecastingData:
    def __init__(self):
        raw_X = self.load_data()
        assert raw_X.ndim == 2
        valid_set = int((1-args.valid_ratio-args.test_ratio)*raw_X.shape[0])
        test_set = int((1-args.test_ratio)*raw_X.shape[0])
        self.raw_Xs = np.split(raw_X, [valid_set, test_set])  # trn, val, tst split

        if args.local_norm:
            axis = 0
        else:
            axis = None

        if args.norm_type == 'none':
            self.sc = np.ones((1, 1))
            self.mn = np.zeros((1, 1))
        elif args.norm_type == 'minmax':
            self.mn = self.raw_Xs[0].min(axis=axis, keepdims=True)
            self.sc = self.raw_Xs[0].max(axis=axis, keepdims=True) - self.mn
        elif args.norm_type == 'standard':
            self.mn = self.raw_Xs[0].mean(axis=axis, keepdims=True)
            self.sc = self.raw_Xs[0].std(axis=axis, keepdims=True)

    def load_data(self):
        data_path = os.path.join(args.data_dir, f'{args.dataset}.npy')
        raw_X = np.load(data_path, allow_pickle=True)
        args.n_series = raw_X[0].shape[-1]
        return raw_X

    def get_dataset(self, i):
        return ForecastingDataset(self.raw_Xs[i], self.sc, self.mn)

class ForecastingDataset(Dataset):
    def __init__(self, raw_X, sc, mn):
        self.sc = sc
        self.mn = mn
        self.rse = np.sum((raw_X[args.series_len-1:] - np.mean(raw_X[args.series_len-1:]))**2)
        self.X = self.norm(raw_X).astype(np.float32)
        self.avg = self.X.mean(axis=0).astype(np.float32)

    def norm(self, X):
        return (X - self.mn) / self.sc

    def renorm(self, X):
        return X * self.sc + self.mn

    def __len__(self):
        return self.X.shape[0] - args.series_len

    def __getitem__(self, idx):
        return (self.X[idx:idx+args.series_len], self.X[idx+args.series_len])

max_stars_repo_name: Daikon-Sun/AdjustAutocorrelation

max_issues_repo_name: Daikon-Sun/AdjustAutocorrelation

max_forks_repo_name: Daikon-Sun/AdjustAutocorrelation


text:
import pydicom as dicom
from astropy.io import fits
from astropy.table import Table
import numpy as np

from ..meta.info import dicom_fieldnames

def dicom2fits(in_dcm_path, out_fits_path):

    ds = dicom.dcmread(in_dcm_path)
    arr = np.array(ds.pixel_array)

    hdu = fits.PrimaryHDU(arr)
    hdr = hdu.header

    rows = []
    for field in dicom_fieldnames:
        try:
            if ds.data_element(field) is not None:
                x = str(ds.data_element(field)).replace("'", "")
                y = x.find(":")
                value = x[y+2:]
                # FITS strings have a max size of 68 characters
                hdr['hierarch {}'.format(field)] = value[:40]
                rows.append(dict(key=field, value=value))
        except Exception as e:
            continue

    # write the meta-data also out as a table in a second HDU
    # there shouldn't be an practical limit on this info
    tbl = Table(rows=rows)
    tbl.meta['comments'] = ["Written by dicom2fits"]
    metadata = fits.BinTableHDU(data=tbl)

    hdu_l = fits.HDUList([hdu, metadata])

    hdu_l.writeto(out_fits_path)

max_stars_repo_name: ejeschke/dicomconvert

max_issues_repo_name: ejeschke/dicomconvert

max_forks_repo_name: ejeschke/dicomconvert

reason: import numpy,from astropy


text:
import struct
import time
import math

import matplotlib.pyplot as plt
from matplotlib import gridspec
import numpy as np

MAGIC_VALUE = [ord('C'), ord('M'), ord('Z')]

QUANTIZATION_BINS_DEFAULT = np.round(np.linspace(0.0, 1.0, 256), decimals=4)

PACK_FLOAT = struct.Struct('=f')
PACK_FLOAT_SIZE = 4
PACK_SHAPE = struct.Struct('=i')
PACK_SHAPE_SIZE = 4
PACK_BOOL = struct.Struct('=?')
PACK_BOOL_SIZE = 1
PACK_QUANT = struct.Struct('=c')
PACK_QUANT_SIZE = 1

class Patcher:

    def __init__(self, img_x_size, img_y_size, img=None,
                 patch_x_size=8, patch_y_size=8):
        self.img_x_size = img_x_size
        self.img_y_size = img_y_size
        self.img = img
        self.patch_x_size = patch_x_size
        self.patch_y_size = patch_y_size

    def get_patch_coordinates(self):
        x_patches = int(self.img_x_size / self.patch_x_size)
        y_patches = int(self.img_y_size / self.patch_y_size)
        patch_coords = [
            (
                x * self.patch_x_size,
                y * self.patch_y_size,
                (x + 1) * self.patch_x_size,
                (y + 1) * self.patch_y_size
            )
            for y in range(y_patches)
            for x in range(x_patches)
        ]
        return patch_coords

    def build_patches_matrix(self):
        patch_coords = self.get_patch_coordinates()
        self.patches_matrix = np.array([])
        if self.img is not None:
            patch_vectors = []
            for patch_coord in patch_coords:
                # produce crops of original image
                patch_img = self.img.crop(patch_coord)
                # get raster of pixel data from patch image
                patch_vector = np.array(list(patch_img.getdata()))
                # normalize pixels (256 gray levels)
                patch_vector = patch_vector / 255.
                patch_vectors.append(patch_vector)
            self.patches_matrix = np.stack(patch_vectors, axis=0)

def measured_task(message, task_callable, *args, **kwargs):
    print(message, end='', flush=True)
    start = time.perf_counter()
    task_callable(*args, **kwargs)
    end = time.perf_counter()
    print('took', (end-start), 'sec(s)')

def plot_errors(history_dump, plot_output_path):

    # fit plot squarely with subplots (approx N x N) if possible
    # total number of plots
    nplots = len(history_dump)
    # number of columns
    cols = int(math.sqrt(nplots))
    # number of rows
    rows = int(math.ceil(nplots / cols))

    gs = gridspec.GridSpec(rows, cols)
    fig = plt.figure()
    legend = [[], []]
    for i, (name, values) in enumerate(history_dump.items()):
        try:
            options = ERROR_PLOTS[name]
            ax = fig.add_subplot(gs[i])
            line, = ax.plot(
                range(1, len(values)+1), values,
                label=options['label'], color=options['color']
            )
            ax.set_xlabel('epoch')
            ax.set_ylabel(options['label'])
            legend[0].append(line)
            legend[1].append(options['legend_name'])
        except KeyError:
            pass
    lgd = fig.legend(
        legend[0], legend[1], loc='lower center',
        bbox_to_anchor=(0.5, -0.25),
        bbox_transform=fig.transFigure,
    )
    fig.tight_layout()
    # save figure and calculate size with legend outside subplots
    fig.savefig(
        plot_output_path, bbox_extra_artists=(lgd,), bbox_inches='tight'
    )

ERROR_PLOTS = {
    'loss': {
        'color': 'blue',
        'label': 'Loss',
        'legend_name': 'Loss',
    },
    'mean_absolute_error': {
        'color': 'red',
        'label': 'MAE',
        'legend_name': 'Mean Absolute Error',
    },
    'mean_absolute_percentage_error': {
        'color': 'green',
        'label': 'MAPE',
        'legend_name': 'Mean Absolute Percentage Error',
    },
    'mean_squared_error': {
        'color': 'cyan',
        'label': 'MSE',
        'legend_name': 'Mean Squared Error',
    },
    'mean_squared_logarithmic_error': {
        'color': 'magenta',
        'label': 'MSLE',
        'legend_name': 'Mean Squared Logarithmic Error',
    }
}


text:
from keras.models import Sequential, Model
from keras.layers import Dense, Reshape, Input, Bidirectional, multiply
from keras.layers import Conv2D
from keras.layers.core import Activation, Flatten, Lambda, Dropout
from keras.layers.advanced_activations import LeakyReLU, ELU
from keras.layers.normalization import BatchNormalization
from keras.layers.recurrent import GRU as keras_GRU
from keras.layers.recurrent import LSTM as keras_LSTM
from keras.layers.wrappers import TimeDistributed
# from keras.engine import merge
from keras import backend as K

import numpy as np

from network import Network

class GRU_net(Network):
    """Gated recurrent unit network. This inherets from network"""
    def __init__(self, network, filter):
        """
            The bare minimum for initializing this network.

            Keyword arguments:
            network -- the name of the network so we can save the weights.
            filter -- the name of the filter that is being used byt the network
            to do source separation.
        """
        super(GRU_net, self).__init__(network)
        self.network = network
        self.filter = filter

    def model(self):
        """
            The definition of the network. We are using the Fuctional API from
            keras. It inherets most of the things from Network where the
            number of units per layer are defined and the dropout parameters.

            Returns:
            model -- the network architecture.
        """
        mix = Input(shape=(None, self.features), name='input')
        gru1 = keras_GRU(self.units[0], return_sequences=True,
                           dropout=self.DROPOUT,
                           recurrent_dropout=self.RDROPOUT)(mix)
        gru2 = keras_GRU(self.units[1], return_sequences=True,
                           dropout=self.DROPOUT,
                           recurrent_dropout=self.RDROPOUT)(gru1)
        gru3 = keras_GRU(self.units[2], return_sequences=True,
                           dropout=self.DROPOUT,
                           recurrent_dropout=self.RDROPOUT)(gru2)
        gru4 = keras_GRU(self.units[3], return_sequences=True,
                           dropout=self.DROPOUT,
                           recurrent_dropout=self.RDROPOUT)(gru3)

        speaker_1, speaker_2 = self.filter(self.features, mix, gru4)
        model = Model(input=[mix], output=[speaker_1, speaker_2])

        return model

class LSTM_net(Network):
    """LSTM network using standard LSTM units"""
    def __init__(self, network, filter):
        """
            The bare minimum for initializing this network.

            Keyword arguments:
            network -- the name of the network so we can save the weights.
            filter -- the name of the filter that is being used byt the network
            to do source separation.
        """
        super(LSTM_net, self).__init__(network)
        self.network = network
        self.filter = filter

    def model(self):
        """
            The definition of the network. We are using the Fuctional API from
            keras. It inherets most of the things from Network where the
            number of units per layer are defined and the dropout parameters.

            Returns:
            model -- the network architecture.
        """
        mix = Input(shape=(None, self.features), name='input')
        lstm1 = keras_LSTM(self.units[0], return_sequences=True,
                           dropout=self.DROPOUT,
                           recurrent_dropout=self.RDROPOUT)(mix)
        lstm2 = keras_LSTM(self.units[1], return_sequences=True,
                           dropout=self.DROPOUT,
                           recurrent_dropout=self.RDROPOUT)(lstm1)
        lstm3 = keras_LSTM(self.units[2], return_sequences=True,
                           dropout=self.DROPOUT,
                           recurrent_dropout=self.RDROPOUT)(lstm2)
        lstm4 = keras_LSTM(self.units[3], return_sequences=True,
                           dropout=self.DROPOUT,
                           recurrent_dropout=self.RDROPOUT)(lstm3)
        tdd1 = TimeDistributed(Dense(self.features, activation='tanh'), name="dense1")(lstm4)
        tdd2 = TimeDistributed(Dense(self.features, activation='tanh'), name="dense2")(tdd1)

        speaker_1, speaker_2 = self.filter(self.features, mix, tdd2)
        model = Model(input=[mix], output=[speaker_1, speaker_2])

        return model

class LSTM_pobs_net(Network):
    """LSTM using the probabilies of the mixture and tries to learn a clean
        separation based on that. It takes a second input which is the original
        mixture and makes the network look more like a NMF.
    """
    def __init__(self, network, filter):
        """
            The bare minimum for initializing this network.

            Keyword arguments:
            network -- the name of the network so we can save the weights.
            filter -- the name of the filter that is being used byt the network
            to do source separation.
        """
        super(LSTM_pobs, self).__init__(network)
        self.network = network
        self.filter = filter

    def model(self):
        """
            The definition of the network. We are using the Fuctional API from
            keras. It inherets most of the things from Network where the
            number of units per layer are defined and the dropout parameters.

            Returns:
            model -- the network architecture.
        """
        pobs = Input(shape=(None, self.features), name='input')
        mix = Input(shape=(None, self.features), name='true_input')
        lstm1 = keras_LSTM(self.units, return_sequences=True,
                           dropout=self.DROPOUT,
                           recurrent_dropout=self.RDROPOUT)(pobs)
        lstm2 = keras_LSTM(self.units, return_sequences=True,
                           dropout=self.DROPOUT,
                           recurrent_dropout=self.RDROPOUT)(lstm1)
        lstm3 = keras_LSTM(self.units, return_sequences=True,
                           dropout=self.DROPOUT,
                           recurrent_dropout=self.RDROPOUT)(lstm2)
        lstm4 = keras_LSTM(self.units, return_sequences=True,
                           dropout=self.DROPOUT,
                           recurrent_dropout=self.RDROPOUT)(lstm3)
        #  drop = Dropout(0.5)(lstm2)

        speaker_1, speaker_2 = self.filter(self.features, mix, lstm4)
        model = Model(input=[pobs, mix], output=[speaker_1, speaker_2])

        return model

class LSTM_backwards_net(Network):
    """LSTM network with forward and backward connections"""
    def __init__(self, network, filter):
        """
            The bare minimum for initializing this network.

            Keyword arguments:
            network -- the name of the network so we can save the weights.
            filter -- the name of the filter that is being used byt the network
            to do source separation.
        """
        super(LSTM_backwards, self).__init__(network)
        self.network = network
        self.filter = filter

    def model(self):
        """
            The definition of the network. We are using the Fuctional API from
            keras. It inherets most of the things from Network where the
            number of units per layer are defined and the dropout parameters.

            Returns:
            model -- the network architecture.
        """
        mix = Input(shape=(None, self.features), name='input')
        lstm1 = Bidirectional(keras_LSTM(self.units, return_sequences=True,
                              dropout=self.DROPOUT,
                              recurrent_dropout=self.RDROPOUT))(mix)
        lstm2 = Bidirectional(keras_LSTM(self.units, return_sequences=True,
                              dropout=self.DROPOUT,
                              recurrent_dropout=self.RDROPOUT))(lstm1)
        lstm3 = Bidirectional(keras_LSTM(self.units, return_sequences=True,
                              dropout=self.DROPOUT,
                              recurrent_dropout=self.RDROPOUT))(lstm2)
        lstm4 = Bidirectional(keras_LSTM(self.units, return_sequences=True,
                              dropout=self.DROPOUT,
                              recurrent_dropout=self.RDROPOUT))(lstm3)
        # drop = Dropout(0.5)(lstm2)

        speaker_1, speaker_2 = self.filter(self.features, mix, lstm4)
        model = Model(input=[mix], output=[speaker_1, speaker_2])

        return model

class DNN_net(Network):
    """Deep neural network which uses dense layers for the separation"""
    def __init__(self, network, filter):
        """
            The bare minimum for initializing this network.

            Keyword arguments:
            network -- the name of the network so we can save the weights.
            filter -- the name of the filter that is being used byt the network
            to do source separation.
        """
        super(DNN_net, self).__init__(network)
        self.network = network
        self.filter = filter

    def model(self):
        """
            The definition of the network. We are using the Fuctional API from
            keras. It inherets most of the things from Network where the
            number of units per layer are defined and the dropout parameters.

            Returns:
            model -- the network architecture.
        """
        mix = Input(shape=(None, self.features), name='input')
        dnn1 = TimeDistributed(Dense(500, activation='tanh'), name="dense1")(mix)
        dnn2 = TimeDistributed(Dense(300, activation='tanh'), name="dense2")(dnn1)
        dnn3 = TimeDistributed(Dense(self.features, activation='tanh'), name="dense3")(dnn2)
        dnn4 = TimeDistributed(Dense(self.features, activation='tanh'), name="dense4")(dnn3)

        speaker_1, speaker_2 = self.filter(self.features, mix, dnn4)
        model = Model(input=[mix], output=[speaker_1, speaker_2])

        return model

class CNN_net(Network):
    """Convolutional neural network using convolutions to do the job"""
    def __init__(self, network, filter):
        """
            The bare minimum for initializing this network.

            Keyword arguments:
            network -- the name of the network so we can save the weights.
            filter -- the name of the filter that is being used byt the nestwork
            to do source separation.
        """
        super(CNN_net, self).__init__(network)
        self.network = network
        self.filter= filter

    def model(self):
        """
            The definition of the network. We are using the Fuctional API from
            keras. It inherets most of the things from Network where the
            number of units per layer are defined and the dropout parameters.

            Returns:
            model -- the network architecture.
        """
        kernels = 64
        kernel_width = 3
        kernel_height = 3

        mix = Input(batch_shape=(None,1,self.timesteps, self.features), name='input')

        conv1 = Conv2D(kernels, (kernel_height, kernel_width), padding='same', data_format='channels_first')(mix)
        conv2 = Conv2D(kernels, (kernel_height, kernel_width), padding='same', data_format='channels_first')(conv1)
        conv3 = Conv2D(kernels, (kernel_height, kernel_width), padding='same', data_format='channels_first')(conv2)
        trans = Conv2D(1, (kernel_height, kernel_width), padding='same', data_format='channels_first')(conv3)
        resh1 = Reshape((self.timesteps, self.features))(trans)
        dnn1 = TimeDistributed(Dense(500, activation='sigmoid'))(trans)
        dnn2 = TimeDistributed(Dense(800, activation='sigmoid'))(dnn1)
        dnn2 = TimeDistributed(Dense(200, activation='sigmoid'))(dnn2)

        speaker_1, speaker_2 = self.filter(self.features, mix, dnn2)

        model = Model(input=[mix], output=[speaker_1, speaker_2])

        return model

class RCNN_net(Network):
    """recurrent CNN using both convolutions and lstms"""
    def __init__(self, network, filter):
        """
            The bare minimum for initializing this network.

            Keyword arguments:
            network -- the name of the network so we can save the weights.
            filter -- the name of the filter that is being used byt the network
            to do source separation.
        """
        super(RCNN, self).__init__(network)
        self.network = network
        self.filter = filter

    def model(self):
        """
            The definition of the network. We are using the Fuctional API from
            keras. It inherets most of the things from Network where the
            number of units per layer are defined and the dropout parameters.

            Returns:
            model -- the network architecture.
        """
        kernels = 64
        kernel_width = 2
        kernel_height = 2

        mix = Input(shape=(self.timesteps, self.features, self.sequences), name='input')
        conv1 = Convolution2D(kernels, kernel_height, kernel_width, border_mode='same')(mix)
        conv2 = Convolution2D(kernels, kernel_height, kernel_width, border_mode='same')(conv1)
        dnn1 = TimeDistributed(Dense(self.features, activation='relu'))(conv2)
        resh1 = Reshape((self.timesteps, self.features), input_shape=(1,self.sequences, self.timesteps, self.features))(dnn1)
        lstm1 = Bidirectional(keras_LSTM(self.features, return_sequences=True))(resh1)
        lstm2 = Bidirectional(keras_LSTM(self.features, return_sequences=True))(lstm1)
        speaker_1, speaker_2 = self.filter(self.features, mix, lstm2)

        model = Model(input=[mix], output=[speaker_1, speaker_2])

        return model

max_stars_repo_name: TeunKrikke/SourceSepDL

max_issues_repo_name: TeunKrikke/SourceSepDL

max_forks_repo_name: TeunKrikke/SourceSepDL


text:
import numpy as np
from numpy import array

class ClassifierWrapper(object):

    def __init__(self, cls, min_cases_for_training = 30):
        self.cls = cls

        self.min_cases_for_training = min_cases_for_training
        self.hardcoded_prediction = None

    def fit(self, X, y):
        # if there are too few training instances, use the mean
        if X.shape[0] < self.min_cases_for_training:
            self.hardcoded_prediction = np.mean(y)

        # if all the training instances are of the same class, use this class as prediction
        elif len(set(y)) < 2:
            self.hardcoded_prediction = int(y[0])

        else:
            self.cls.fit(X, y)
            return self

    def predict_proba(self, X, y=None):
        if self.hardcoded_prediction is not None:
            return array([self.hardcoded_prediction] * X.shape[0])

        else:
            #preds_pos_label_idx = np.where(self.cls.classes_ == 0)[0][0]
            #preds = self.cls.predict_proba(X)[:,preds_pos_label_idx]
            preds = self.cls.predict(X)
            return preds

    def fit_predict(self, X, y):

        self.fit(X, y)
        return self.predict_proba(X)

max_stars_repo_name: verenich/time-prediction-benchmark

max_issues_repo_name: verenich/time-prediction-benchmark

max_forks_repo_name: verenich/time-prediction-benchmark


text:
import numpy as np

def initial_velocity(parameters):
  """
  Initializes the velocity as a python dictionary with:
              - keys: "dW1", "db1", ..., "dWL", "dbL"
              - values: numpy arrays of zeros of the same shape as the corresponding gradients/parameters.
  Arguments:
  parameters -- python dictionary containing your parameters.
                  parameters['W' + str(l)] = Wl
                  parameters['b' + str(l)] = bl

  Returns:
  v -- python dictionary containing the current velocity.
                  v['dW' + str(l)] = velocity of dWl
                  v['db' + str(l)] = velocity of dbl
  """

  L = len(parameters) // 2 # number of layers in the neural networks
  v = {}

  # Initialize velocity
  for l in range(L):
      v["dW" + str(l+1)] = np.zeros(parameters["W{}".format(l+1)].shape)
      v["db" + str(l+1)] = np.zeros(parameters["b{}".format(l+1)].shape)

  return v

def update_parameters(parameters, gradients, v, beta, learning_rate):
  """
  Update parameters using Momentum

  Arguments:
  parameters -- python dictionary containing your parameters:
                  parameters['W' + str(l)] = Wl
                  parameters['b' + str(l)] = bl
  gradients -- python dictionary containing your gradients for each parameters:
                  gradients['dW' + str(l)] = dWl
                  gradients['db' + str(l)] = dbl
  v -- python dictionary containing the current velocity:
                  v['dW' + str(l)] = ...
                  v['db' + str(l)] = ...
  beta -- the momentum hyperparameter, scalar
  learning_rate -- the learning rate, scalar

  Returns:
  parameters -- python dictionary containing your updated parameters
  v -- python dictionary containing your updated velocities
  """

  L = len(parameters) // 2 # number of layers in the neural networks

  for l in range(L):
    # compute velocities
    v["dW" + str(l+1)] = beta*v['dW{}'.format(l+1)] + (1-beta) * gradients['dW{}'.format(l+1)]
    v["db" + str(l+1)] = beta*v['db{}'.format(l+1)] + (1-beta) * gradients['db{}'.format(l+1)]
    # update parameters
    parameters["W" + str(l+1)] = parameters["W" + str(l+1)] - learning_rate * v['dW{}'.format(l+1)]
    parameters["b" + str(l+1)] = parameters["b" + str(l+1)] - learning_rate * v['db{}'.format(l+1)]

  return parameters, v


text:
# -*- coding: utf-8 -*-
"""
Created on Thu Apr 29 17:06:47 2021

@author: arnep

Estimation of the extremal index of data generated by an
an AR(1) model with student-t distributed errors and stable coefficient.

Estimators used are the parametric kappa filter, the stable and rolling horizon
estimator based on Cai (2019), and the sliding block estimators of Berghaus &
Bücher (2018).
"""

import numpy as np
import matplotlib.pyplot as plt
from methods_simulation import ground_truth_ARt, random_ARt, MSE_miss_handler
from methods_simulation import est_ARt_parametric, est_non_parametr_stable, est_non_parametr_rol_hor
from methods_simulation import est_sliding_blocks_Z, est_sliding_blocks_Y, theta_simulator
from methods_extremal_index import MLE, construct_psi_k2, construct_psi_robust_k2
from methods_extremal_index import non_parametric_day_based

## simulation parameters
ph = 0.8
n = 1000
runs = 10
phi= np.ones(n)*ph
horizon_length = 400
dof = 3
func_gt = ground_truth_ARt
func_rand = random_ARt
func_crit = MSE_miss_handler

#object names
name = 'MSE_ARt_stable_phi_' +str(ph)

#estimation parameters
func_non_param = non_parametric_day_based
k = 50
k2 =  int(round(k/n*horizon_length)*1.2)
bn = int(round(n/k)*4)
d_u = 10
func_mle = MLE
func_filter = construct_psi_robust_k2

### est with robust filter
theta_init = [0.8, 0, 0,  3]
func_filter2 = construct_psi_k2
theta_init2 = [ 0.8,  0, 0,  3]

bounds = ((-2,2), #omega
        (-.9999,  .9999),#alpha
        (-.9999,  .9999), #beta
        (2, 10001) #nu
        )

#create simulation object
sim_MSE_ARt = theta_simulator(name, n, runs, func_gt, func_rand, func_crit, dof =dof, phi=phi )

#estimate parametric model
func_est = est_ARt_parametric
print('Estimating parametric model:')
sim_MSE_ARt.simulate(func_est, name = 'Parametric ARt(1) filter, kappa: 10',
                     func_filter = func_filter2, theta_init =theta_init2,
                     bounds = bounds, func_mle = func_mle)

#estimate non-parametric estimator of Cai (2019) (stable and rolling horizon)
func_est = est_non_parametr_stable
print('Estimating stable non-parametric estimator:')
sim_MSE_ARt.simulate(func_est, name = 'Stable non-parametric estimator, k: ' + str(k),
                     d_u = d_u, k = k, func = func_non_param)

func_est = est_non_parametr_rol_hor
print('Estimating rolling horizon non-parametric estimator:')
sim_MSE_ARt.simulate(func_est, name = 'Rolling horizon non-parametric estimator, k: ' + str(k2) +
                     ', horizon length: ' + str(horizon_length),
                     horizon_length= horizon_length, d_u = d_u, k = k2,
                     func = func_non_param)

#estimate sliding blocks estimator of Berghaus & Bucher (2018)
func_est =  est_sliding_blocks_Y
print('Estimating Y-estimator:')
sim_MSE_ARt.simulate(func_est, name = 'Sliding blocks Y-estimator, bn: ' + str(bn),
                     bn = bn)

func_est = est_sliding_blocks_Z
print('Estimating Z-estimator:')
sim_MSE_ARt.simulate(func_est, name = 'Sliding blocks Z-estimator, bn: ' + str(bn),
                     bn = bn)

#get the summary
sim_MSE_ARt.calc_stats()

#plot all the estimators
for key in sim_MSE_ARt.theta.keys():
    sim_MSE_ARt.plot_theta_gt(key)

#plot ground truth and simulation example
sim_MSE_ARt.plot_gt('ARt(1) model with phi = ' + str(ph) )
sim_MSE_ARt.plot_sim_example('ARt(1) model with phi = ' + str(ph))

#plot autoregressive coefficient
plt.figure(figsize=(10, 6), dpi=200)
plt.grid(color = 'black')
plt.plot(phi)
plt.xticks(np.arange(0, n*1.2, n/5))
plt.yticks(np.arange(0, 1.2, .2))
plt.title('$\phi_t$')
plt.show()


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
:mod:`run`
==================

.. module:: run
    :platform: Unix, Windows
    :synopsis:

.. moduleauthor:: hbldh <henrik.blidh@nedomkull.com>

Created on 2015-11-06, 14:24

"""

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

import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import zero_one_loss
from sklearn.tree import DecisionTreeClassifier

from skboost.milboost.classifier import MILBoostClassifier
from skboost.datasets import MUSK1, MUSK2, Hastie_10_2
from skboost.milboost.softmax import *

#ds_train = MUSK2()
#print(ds_train)
#ds_test = MUSK1()
#print(ds_test)
ds_train = Hastie_10_2()
ds_test = Hastie_10_2(random_state=54325)

mil_classifier = MILBoostClassifier(
    base_estimator=DecisionTreeClassifier(max_depth=1),
    #softmax=ISR(),
    #softmax=GeneralizedMean(5.0),
    #softmax=NoisyOR(),
    softmax=LogSumExponential(5.0),
    n_estimators=30,
    learning_rate=1.0,
    verbose=True)
mil_classifier.fit(ds_train.data, ds_train.labels)
print(mil_classifier)

# Calculate step-by-step results per instance on training data.
milboost_train_error = np.zeros((len(mil_classifier.estimators_),))
for i, y_pred in enumerate(mil_classifier.staged_predict(ds_train.data)):
    milboost_train_error[i] = zero_one_loss(y_pred, np.sign(ds_train.labels))

# Calculate step-by-step results per instance on test data.
milboost_test_error = np.zeros((len(mil_classifier.estimators_),))
for i, y_pred in enumerate(mil_classifier.staged_predict(ds_test.data)):
    milboost_test_error[i] = zero_one_loss(y_pred, np.sign(ds_test.labels))

# Calculate step-by-step results per bag on training data.
milboost_bag_train_error = np.zeros((len(mil_classifier.estimators_),))
for i, y_pred in enumerate(mil_classifier.staged_predict(ds_train.data)):
    bag_labels_pred = np.array([np.max(x) for x in np.split(y_pred, ds_train.bag_partitioning)], 'int')
    milboost_bag_train_error[i] = zero_one_loss(bag_labels_pred, ds_train.bag_labels)

milboost_bag_test_error = np.zeros((len(mil_classifier.estimators_),))
for i, y_pred in enumerate(mil_classifier.staged_predict(ds_test.data)):
    bag_labels_pred = np.array([np.max(x) for x in np.split(y_pred, ds_test.bag_partitioning)], 'int')
    milboost_bag_test_error[i] = zero_one_loss(bag_labels_pred, np.sign(ds_test.bag_labels))

plt.plot(np.arange(len(mil_classifier.estimators_)) + 1, milboost_test_error, '{0}'.format('r'),
         label='{0} Test Error'.format('MILBoost, {0}'.format(mil_classifier.softmax_fcn)))
plt.plot(np.arange(len(mil_classifier.estimators_)) + 1, milboost_train_error, '{0}--'.format('g'),
         label='{0} Train Error'.format('MILBoost, {0}'.format(mil_classifier.softmax_fcn)))
plt.plot(np.arange(len(mil_classifier.estimators_)) + 1, milboost_bag_test_error, '{0}'.format('b'),
         label='{0} Bag Test Error'.format('MILBoost, {0}'.format(mil_classifier.softmax_fcn)))
plt.plot(np.arange(len(mil_classifier.estimators_)) + 1, milboost_bag_train_error, '{0}--'.format('c'),
         label='{0} Bag Train Error'.format('MILBoost, {0}'.format(mil_classifier.softmax_fcn)))

plt.gca().set_ylim((0.0, 0.5))
plt.xlabel('n_estimators')
plt.ylabel('error rate')
plt.legend(loc='upper right', fancybox=True)

plt.show()

max_stars_repo_name: TMRert/skboost

max_issues_repo_name: TMRert/skboost

max_forks_repo_name: TMRert/skboost


text:
import numpy as np

def add_tuples(a: tuple[int, ...], b: tuple[int, ...]) -> tuple[int, ...]:
    return tuple(p + q for p, q in zip(a, b))

def float_dict_comparison(d0: dict[any, float], d1: dict[any, float]) -> bool:
    return np.all([np.isclose(d0[x], d1[x]) for x in d0]) and d0.keys() == d1.keys()

max_stars_repo_name: hmeretti/grid-world

max_issues_repo_name: hmeretti/grid-world

max_forks_repo_name: hmeretti/grid-world


text:
"""States are containers for occupancy values.  States hold a
number of numpy arrays which are all of the same length.

state.q_val holds a numpy array of quantity values

These are the most important values in the system, they
show the quantity of material (in moles) of a given species
in a given compartment.

The rest of the arrays are useful for making selections:

state.species is a numpy (string) array of species IDs
state.compartment is a numpy (string) array of compartment IDs
state.x_pos is a numpy array of x_positions of the compartment centers
state.y_pos is a numpy array of y_positions of the compartment centers
state.z_pos is a numpy array of z_positions of the compartment centers

Note that species values for a given compartment will only be created
if there is either:
1) a Reaction in that compartment which involves that Species
or 2) a connection from that compartment which involves that Species

When initializing the state, an index is constructed to easily
determine an index given a compartment and a species:

index = state.index[compID][specID]
"""

from openrxn import unit
from openrxn.model import FlatModel
import numpy as np
import pandas as pd

class State(object):
    def __init__(self, model=None, dataframe=None, units=[unit.nanometer]*3):
        """State objects can be initialized using either a
        FlatModel or a dataframe object.  At minimum, the
        dataframe needs to have "species" and "compartment"
        columns."""

        self.index = {}
        self.units = units

        if model is not None:
            assert isinstance(model,FlatModel), "Error! A state object needs a FlatModel to initialize."
            self._init_from_model(model)
        elif dataframe is not None:
            if 'species' not in dataframe.columns or 'compartment' not in dataframe.columns:
                raise ValueError("Error! dataframe must contain columns for 'species' and 'compartment'")
            self._init_from_df(dataframe)

        self.size = len(self.compartment)
        self.q_val = np.zeros((self.size))

    def _init_from_model(self, model):

        big_species_list = []
        big_comp_list = []
        big_x_list = []
        big_y_list = []
        big_z_list = []

        running_index = 0
        for c_tag, c in model.compartments.items():

            self.index[c_tag] = {}

            # figure out which species are associated with this compartment
            spec = []
            for other_c, conn in c.connections.items():
                spec += list(conn[1].species_rates.keys())
            for rxn in c.reactions:
                spec += rxn.reactant_IDs
                spec += rxn.product_IDs
            spec_set = set(spec)

            for i,s in enumerate(spec_set):
                self.index[c_tag][s] = running_index + i

            running_index += len(spec_set)

            big_species_list += list(spec_set)
            big_comp_list += [c.ID]*len(spec_set)

            # for x, y and z, average the boundary values
            x = [None,None,None]
            for i in range(len(c.pos)):
                x[i] = 0.5*(c.pos[i][0]+c.pos[i][1]).to(self.units[i]).magnitude

            big_x_list += [x[0]]*len(spec_set)
            big_y_list += [x[1]]*len(spec_set)
            big_z_list += [x[2]]*len(spec_set)

        self.species = np.array(big_species_list)
        self.compartment = np.array(big_comp_list)
        self.x_pos = np.array(big_x_list)
        self.y_pos = np.array(big_y_list)
        self.z_pos = np.array(big_z_list)

    def _init_from_df(self, df):

        # assign columns to self arrays
        self.species = np.array(df['species'])
        self.compartment = np.array(df['compartment'])
        if 'x_pos' in df.columns:
            self.x_pos = np.array(df['x_pos'])
        if 'y_pos' in df.columns:
            self.y_pos = np.array(df['y_pos'])
        if 'z_pos' in df.columns:
            self.z_pos = np.array(df['z_pos'])

        # building self.index dictionary
        for i in range(len(df['species'])):
            if df['compartment'][i] not in self.index:
                self.index[df['compartment'][i]] = {}
            self.index[df['compartment'][i]][df['species'][i]] = i

    def to_dataframe(self):
        df = pd.DataFrame()

        df['species'] = self.species
        df['compartment'] = self.compartment
        df['q_val'] = self.q_val

        if hasattr(self,'x_pos'):
            df['x_pos'] = self.x_pos
        if hasattr(self,'y_pos'):
            df['y_pos'] = self.y_pos
        if hasattr(self,'z_pos'):
            df['z_pos'] = self.z_pos

        return df

    def to_csv(self, filename):
        df = self.to_dataframe()
        df.to_csv(filename)

    def to_csv_no_q(self, filename):
        df = self.to_dataframe()

        cols = [c for c in df.columns if c != 'q_val']
        df.to_csv(filename,columns=cols)

max_stars_repo_name: alexrd/OpenRXN

max_issues_repo_name: ADicksonLab/OpenRXN

max_forks_repo_name: ADicksonLab/OpenRXN


text:
# BSD 3-Clause License; see https://github.com/scikit-hep/awkward-1.0/blob/main/LICENSE

import pytest  # noqa: F401
import numpy as np  # noqa: F401
import awkward as ak  # noqa: F401

def test_na_union():
    one = ak.from_iter([1, None, 3], highlevel=False)
    two = ak.from_iter([[], [1], None, [3, 3, 3]], highlevel=False)
    tags = ak.layout.Index8(np.array([0, 1, 1, 0, 0, 1, 1], dtype=np.int8))
    index = ak.layout.Index64(np.array([0, 0, 1, 1, 2, 2, 3], dtype=np.int64))
    array = ak.Array(
        ak.layout.UnionArray8_64(tags, index, [one, two]), check_valid=True
    )
    assert ak.to_list(array) == [1, [], [1], None, 3, None, [3, 3, 3]]

    assert ak.to_list(ak.is_none(array)) == [
        False,
        False,
        False,
        True,
        False,
        True,
        False,
    ]

class DummyRecord(ak.Record):
    def __repr__(self):
        return f"<{self.x}>"

class DummyArray(ak.Array):
    def __repr__(self):
        return "<DummyArray {}>".format(" ".join(repr(x) for x in self))

class DeepDummyArray(ak.Array):
    def __repr__(self):
        return "<DeepDummyArray {}>".format(" ".join(repr(x) for x in self))

def test_behaviors():
    behavior = {}
    behavior["Dummy"] = DummyRecord
    behavior[".", "Dummy"] = DummyArray
    behavior["*", "Dummy"] = DeepDummyArray

    content = ak.layout.NumpyArray(np.array([1.1, 2.2, 3.3, 4.4, 5.5]))
    recordarray = ak.layout.RecordArray({"x": content})
    recordarray.setparameter("__record__", "Dummy")

    array = ak.Array(recordarray, behavior=behavior, check_valid=True)
    assert repr(array) == "<DummyArray <1.1> <2.2> <3.3> <4.4> <5.5>>"
    assert repr(array[0]) == "<1.1>"

    offsets = ak.layout.Index64(np.array([0, 3, 3, 5], dtype=np.int64))
    listoffsetarray = ak.layout.ListOffsetArray64(offsets, recordarray)

    array2 = ak.Array(listoffsetarray, behavior=behavior, check_valid=True)

    assert array2.layout.parameter("__record__") is None
    assert array2.layout.purelist_parameter("__record__") == "Dummy"

    assert (
        repr(array2)
        == "<DeepDummyArray <DummyArray <1.1> <2.2> <3.3>> <DummyArray > <DummyArray <4.4> <5.5>>>"
    )
    assert repr(array2[0]) == "<DummyArray <1.1> <2.2> <3.3>>"
    assert repr(array2[0, 0]) == "<1.1>"

    recordarray2 = ak.layout.RecordArray({"outer": listoffsetarray})

    array3 = ak.Array(recordarray2, behavior=behavior, check_valid=True)
    assert type(array3) is ak.Array
    assert type(array3["outer"]) is DeepDummyArray
    assert (
        repr(array3["outer"])
        == "<DeepDummyArray <DummyArray <1.1> <2.2> <3.3>> <DummyArray > <DummyArray <4.4> <5.5>>>"
    )

def test_flatten():
    assert ak.to_list(
        ak.flatten(
            ak.Array([[1.1, 2.2, 3.3], [], [4.4, 5.5]], check_valid=True), axis=1
        )
    ) == [1.1, 2.2, 3.3, 4.4, 5.5]

def test_string_equal():
    trials = [
        (
            ["one", "two", "", "three", "four", "", "five", "six", ""],
            ["one", "two", "", "three", "four", "", "five", "six", ""],
        ),
        (
            ["one", "two", "", "three", "four", "", "five", "six"],
            ["one", "two", "", "three", "four", "", "five", "six"],
        ),
        (
            ["one", "two", "", "three", "four", "", "five", "six", ""],
            ["one", "Two", "", "threE", "four", "", "five", "siX", ""],
        ),
        (
            ["one", "two", "", "three", "four", "", "five", "six"],
            ["one", "Two", "", "threE", "four", "", "five", "siX"],
        ),
        (
            ["one", "two", "", "thre", "four", "", "five", "six", ""],
            ["one", "two", "", "three", "four", "", "five", "six", ""],
        ),
        (
            ["one", "two", "", "thre", "four", "", "five", "six"],
            ["one", "two", "", "three", "four", "", "five", "six"],
        ),
        (
            ["one", "two", "", "three", "four", "", "five", "six", ""],
            ["one", "two", ":)", "three", "four", "", "five", "six", ""],
        ),
        (
            ["one", "two", "", "three", "four", "", "five", "six"],
            ["one", "two", ":)", "three", "four", "", "five", "six"],
        ),
        (
            ["one", "two", "", "three", "four", "", "five", "six", ""],
            ["", "two", "", "three", "four", "", "five", "six", ""],
        ),
        (
            ["one", "two", "", "three", "four", "", "five", "six"],
            ["", "two", "", "three", "four", "", "five", "six"],
        ),
    ]

    for left, right in trials:
        assert ak.to_list(
            ak.Array(left, check_valid=True) == ak.Array(right, check_valid=True)
        ) == [x == y for x, y in zip(left, right)]

def test_string_equal2():
    assert ak.to_list(
        ak.Array(["one", "two", "three", "two", "two", "one"], check_valid=True)
        == "two"
    ) == [False, True, False, True, True, False]


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
A benchmark of SMPyBandits Policies, for Single-Players multi-armed Bandits algorithms.

- Homepage for the benchmarks: https://github.com/Naereen/SMPyBandits-benchmarks/
- Homepage: https://SMPyBandits.GitHub.io/
- Author: Lilian Besson and contributors
- License: MIT
- Date: October 2019
"""
from __future__ import division, print_function  # Python 2 compatibility

__author__ = "Lilian Besson"
__version__ = "0.1"

# Write the benchmarking functions here.
# See "Writing benchmarks" in the asv docs for more information.

import numpy as np

import SMPyBandits
from SMPyBandits import Arms
from SMPyBandits.Environment import MAB
from SMPyBandits import Policies

# Tries to know number of CPU
try:
    from multiprocessing import cpu_count
    CPU_COUNT = cpu_count()  #: Number of CPU on the local machine
except ImportError:
    CPU_COUNT = 1

# ------------------------------------------------
# ----------------- For Policies -----------------

min_arm, max_arm = 0.1, 0.9
make_MAB = lambda nbArms: SMPyBandits.Environment.MAB({'arm_type': Arms.Bernoulli, 'params': tuple(np.linspace(min_arm, max_arm, nbArms))})

algorithm_map = {
    "Uniform": Policies.Uniform,
    "Exp3": Policies.EpsilonDecreasing,
    "UCB": Policies.UCB,
    "Thompson": Policies.Thompson,
    "kl-UCB": Policies.klUCB,
    "MOSS-Anytime": Policies.MOSSAnytime,
    "AdBandits": Policies.AdBandits,
    "ApproximatedFHGittins": Policies.ApproximatedFHGittins,
}
if CPU_COUNT >= 8:
    algorithm_map.update({
        "Softmax": Policies.SoftmaxDecreasing,
        "Exp3++": Policies.Exp3PlusPlus,
        "Discounted-UCB": Policies.DiscountedUCB,
        "Discounted-Thompson": Policies.DiscountedThompson,
        "kl-UCB+": Policies.klUCBPlus,
        "kl-UCB++": Policies.klUCBPlusPlus,
        "kl-UCB-Switch": Policies.klUCBswitchAnytime,
        "BESA": Policies.BESA,
        "RCB": Policies.RCB,
        "PHE": Policies.PHE,
        "OC-UCB": Policies.OCUCB,
        "OSSB": Policies.OSSB,
        "Bayes-UCB": Policies.BayesUCB,
        "Tsallis-Inf": Policies.TsallisInf,
        "UCBoostEpsilon": Policies.UCBoostEpsilon,
    })

values_algorithm = list(algorithm_map.keys())

values_nbArms = [
    2,
]
if CPU_COUNT >= 8:
    values_nbArms += [
        3, 4, 5, 6, 7, 8, 9,
        # 12, 16, 24, 32, 48, 64,  # TODO
    ]
    # max_nbArms = 32  # XXX
    # values_nbArms = list(range(2, max_nbArms + 1))  # XXX

values_horizon = [100, 250, 500, 750]
# values_horizon += [
#     1000, 1250, 1500, 1750,
# ]
if CPU_COUNT >= 8:
    values_horizon += [
        2000, 2500,  # XXX
        3000, 3500,  # XXX
        4000, 4500,  # XXX
        5000, 5500,  # XXX
        6000, 6500,  # XXX
        7000, 7500,  # XXX
        8000, 8500,  # XXX
        9000, 9500,  # XXX
        # 10000, 15000,  # XXX
        # 20000, 25000,  # XXX
        # 30000,  # XXX
    ]

print("values_algorithm =", values_algorithm)  # DEBUG
print("values_nbArms =", values_nbArms)  # DEBUG
print("values_horizon =", values_horizon)  # DEBUG

class SP:
    """
    A benchmark of SMPyBandits Policies. In progress.

    - https://asv.readthedocs.io/en/stable/benchmarks.html#timing-benchmarks
    """
    processes = CPU_COUNT
    # repeat = (10, 200 if CPU_COUNT >= 8 else 20, 4800)
    # # number = 100
    # timeout = 4800

    params = [
        values_algorithm,
        values_nbArms,
        values_horizon,
    ]
    param_names = [
        "algorithm",
        "nbArms",
        "horizon",
    ]

    def setup(self, algname, nbArms, horizon):
        self.MAB = make_MAB(nbArms)
        self.algorithm = algorithm_map[algname]
        self.nbArms = nbArms
        self.horizon = horizon

    # ------- Simulation function -------

    def full_simulation(self, algname, nbArms, horizon):
        MAB = make_MAB(nbArms)
        alg = algorithm_map[algname](nbArms)
        alg.startGame()
        sumReward = 0
        choices = [-1] * horizon
        for t in range(horizon):
            arm = alg.choice()
            reward = MAB.draw(arm)
            sumReward += reward
            alg.getReward(arm, reward)
            choices[t] = arm
        bestArmChoice = len([c for c in choices if c == (nbArms - 1)])
        return alg, sumReward, bestArmChoice

    # ------- Memory benchmarks -------
    # https://asv.readthedocs.io/en/stable/writing_benchmarks.html#memory

    mem_createAlgorithm = full_simulation

    # ------- Peak memory benchmarks -------
    # https://asv.readthedocs.io/en/stable/writing_benchmarks.html#peak-memory

    peakmem_createAlgorithm = full_simulation

    # ------- Timing benchmarks -------
    # https://asv.readthedocs.io/en/stable/writing_benchmarks.html#timing

    time_choice_and_getReward = full_simulation

    # ------- Tracking benchmarks -------
    # https://asv.readthedocs.io/en/stable/writing_benchmarks.html#tracking

    def track_sumReward(self, algname, nbArms, horizon):
        _, sumReward, _ = self.full_simulation(algname, nbArms, horizon)
        return sumReward
    track_sumReward.unit = "Sum reward"

    def track_meanReward(self, algname, nbArms, horizon):
        _, sumReward, _ = self.full_simulation(algname, nbArms, horizon)
        return sumReward / horizon
    track_meanReward.unit = "Mean reward"

    def track_regret(self, algname, nbArms, horizon):
        MAB = make_MAB(nbArms)
        sumReward = self.track_sumReward(algname, nbArms, horizon)
        sumBestReward = MAB.maxArm * horizon
        return max(0, sumBestReward - sumReward)
    track_regret.unit = "Sum regret"

    def track_meanRegret(self, algname, nbArms, horizon):
        return self.track_regret(algname, nbArms, horizon) / horizon
    track_meanRegret.unit = "Mean regret"

    def track_normalizedRegret(self, algname, nbArms, horizon):
        return self.track_regret(algname, nbArms, horizon) / np.log(horizon)
    track_normalizedRegret.unit = "Normalized regret"

    def track_bestArmChoice(self, algname, nbArms, horizon):
        _, _, bestArmChoice = self.full_simulation(algname, nbArms, horizon)
        return bestArmChoice
    track_bestArmChoice.unit = "Best arm selections"

    def track_bestArmChoiceRate(self, algname, nbArms, horizon):
        _, _, bestArmChoice = self.full_simulation(algname, nbArms, horizon)
        return bestArmChoice / horizon
    track_bestArmChoiceRate.unit = "Best arm sample rate"

max_stars_repo_name: Naereen/SMPyBandits-benchmarks

max_issues_repo_name: Naereen/SMPyBandits-benchmarks

max_forks_repo_name: Naereen/SMPyBandits-benchmarks