text:
#!/usr/bin/env python3
# Copyright 2020 Graphcore Ltd.
import argparse
import sys
import numpy as np
import tensorflow as tf
from tensorflow.python.ipu import ipu_compiler, scopes, utils
from tensorflow.python.framework import errors

def device_connection_type(value):
    dcts = {"ALWAYS": utils.DeviceConnectionType.ALWAYS,
            "ON_DEMAND": utils.DeviceConnectionType.ON_DEMAND,
            "NEVER": utils.DeviceConnectionType.NEVER}
    return dcts.get(value)

def my_graph(pa, pb, pc):
    o1 = pa + pb
    o2 = pa + pc
    result = o1 + o2
    return result

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("--connection_type",
                        choices=['ALWAYS', 'ON_DEMAND', 'NEVER'],
                        help="Specify connection type")
    parser.set_defaults(connection_type='ALWAYS')
    opts = parser.parse_args()

    with tf.device("cpu"):
        pa = tf.compat.v1.placeholder(np.float32, [2], name="a")
        pb = tf.compat.v1.placeholder(np.float32, [2], name="b")
        pc = tf.compat.v1.placeholder(np.float32, [2], name="c")

    # Create the IPU section of the graph.
    with scopes.ipu_scope("/device:IPU:0"):
        out = ipu_compiler.compile(my_graph, [pa, pb, pc])

    # Define the feed_dict input data.
    fd = {pa: [1., 1.], pb: [0., 1.], pc: [1., 5.]}

    # Connection type from options.
    connection_type = device_connection_type(opts.connection_type)

    cfg = utils.create_ipu_config()
    cfg = utils.auto_select_ipus(cfg, 1)
    cfg = utils.set_ipu_connection_type(cfg,
                                        connection_type,
                                        1)
    utils.configure_ipu_system(cfg)

    # Run the session.
    # If running with DeviceConnectionType.NEVER then anticipate the
    # specific exception with message "configured for compilation only".
    with tf.compat.v1.Session() as sess:
        try:
            result = sess.run(out, fd)
            print(result)
        except tf.errors.InvalidArgumentError as invalid_arg_exception:
            if (connection_type == utils.DeviceConnectionType.NEVER) and \
               ("configured for compilation only" in invalid_arg_exception.message):
                print("Compiled")
                pass
            else:
                print("ERROR: {}".format(invalid_arg_exception.message))
        except:
            general_exception = sys.exc_info()[0]
            print("ERROR: {}".format(general_exception))

if __name__ == '__main__':
    main()


text:
# Uses python3
import sys
from numpy import random

def binary_search(a, x):
    left, right = 0, len(a)-1

    while left <= right:
        mid = (left + right) // 2
        if x == a[mid]:
            return mid
        elif x < a[mid]:
            right = mid - 1
        else:
            left = mid + 1

    return -1

def linear_search(a, x):
    for i in range(len(a)):
        if a[i] == x:
            return i
    return -1

def stress_test():
    flag_correct = True
    i = 0
    while flag_correct:
        N = random.randint(1, 10)
        x = random.randint(1, 10)
        a = random.randint(0, 10, size=N)
        a.sort()

        print(i, a, x, N)

        ls_res = linear_search(a, x)
        bs_res = binary_search(a, x)

        if ls_res != bs_res:
            print(i, a, x, N)
            print('Linear search:', ls_res)
            print('Binary search:', bs_res)
            flag_correct = False

        i += 1

if __name__ == '__main__':
    input = sys.stdin.read()

    # stress_test()

    data = list(map(int, input.split()))

    n = data[0]
    m = data[n + 1]
    a = data[1:n + 1]

    for x in data[n + 2:]:
        # replace with the call to binary_search when implemented
        print(binary_search(a, x), end=' ')

ext: py

max_stars_repo_path: 1 Algorithmic Toolbox/Homeworks/4 Divide-and-Conquer/1 Binary Search/Python/binary_search.py

reason: from numpy


text:
from vectorhub.encoders.audio.pytorch.wav2vec import Wav2Vec
import numpy as np

def test_fairseq_initialize():
    """
    Testing for fairseq initialize
    """

    client = Wav2Vec()
    assert True

def test_fairseq_encoder():
    """
    Testing for fairseq single encode
    """
    client = Wav2Vec()
    audio = client.read(
        'https://vecsearch-bucket.s3.us-east-2.amazonaws.com/voices/common_voice_en_2.wav')
    output = client.encode(audio)
    assert len(output) == 512

reason: import numpy


text:
"""
Import Atyeo data, tensor formation, plotting raw data.
"""
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

"""Data Import"""
CovDataUrl = 'https://raw.githubusercontent.com/meyer-lab/sysSerol-COVID/main/Data%20File/CovData.csv'
Cov = pd.read_csv(CovDataUrl)

def tensorView():
    """Tensor Formation (Patients x Receptors x Antigens)"""
    _, c = Cov.shape
    SMat = Cov.iloc[:,13:c-1:3]
    RMat = Cov.iloc[:,14:c:3]
    NMat = Cov.iloc[:,15:c:3]
    Tensor = np.stack((SMat, RMat, NMat))
    return(SMat, RMat, NMat, Tensor)

def dataPlot():
    "Plotting raw data, color coded by Convalescent/Deceased"
    Plot1 = sns.scatterplot(y="S IgG3", x="S FcRg2A", data=Cov, hue="Outcome")
    Plot1.set(yscale="log")


text:
#!/usr/bin/env python

import argparse
import numpy as np
import pymesh

"""
Highlight self-intersecting faces.
"""

def parse_args():
    parser = argparse.ArgumentParser(description=__doc__);
    parser.add_argument("input_mesh", help="input mesh");
    parser.add_argument("output_mesh", help="output mesh");
    return parser.parse_args();

def main():
    args = parse_args();
    mesh = pymesh.load_mesh(args.input_mesh);

    intersecting_faces = pymesh.detect_self_intersection(mesh);
    intersection_marker = np.zeros(mesh.num_faces);
    for i,face_pair in enumerate(intersecting_faces):
        intersection_marker[face_pair] = i+1;

    mesh.add_attribute("intersecting_faces");
    mesh.set_attribute("intersecting_faces", intersection_marker);
    mesh.add_attribute("is_self_intersecting");
    mesh.set_attribute("is_self_intersecting", intersection_marker > 0);
    pymesh.save_mesh(args.output_mesh, mesh, "intersecting_faces",
            "is_self_intersecting");

if __name__ == "__main__":
    main();


text:
import os
import unittest
import numpy as np

from food_classifier import inference

__author__ = 'Junior Teudjio'

class TestInference(unittest.TestCase):
    @classmethod
    def setUpClass(cls):
        cls.model = inference.Model()

    def test_predict(self):
        sandwich_pred =self.model.predict('../tests/images/sandwich/sandwich.jpg')
        sushi_pred =self.model.predict('../tests/images/sushi/sushi.jpg')

        self.assertEqual(sandwich_pred['label'], 'sandwich')
        self.assertEqual(sushi_pred['label'], 'sushi')

    def test_get_best_model_path(self):
        ckpts_path = '../tests/fake_model_ckpts'
        best_base_model_path = inference._get_best_model_path('base_model',ckpts_path)
        best_mobnet_model_path = inference._get_best_model_path('mobilenet_model', ckpts_path)

        self.assertEqual(best_base_model_path,
                         os.path.join(ckpts_path, 'base_model.ckpt-03-10.64-0.70.hdf5') )
        self.assertEqual(best_mobnet_model_path,
                         os.path.join(ckpts_path, 'mobilenet_model.ckpt-01-10.64-0.90.hdf5') )

    def test_predict_on_batch(self):
        true_labels, pred_labels = self.model.predict_on_batch(folder_path='../tests/images/')
        self.assertTrue(np.allclose(true_labels, pred_labels))
        self.assertTrue(np.allclose(np.array([0., 1.]), true_labels))

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


text:
#!/usr/bin/env python2.7
'''
Copyright 2018 Google LLC

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    https://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
'''

from __future__ import print_function

import numpy as np
import os
import random
from builtins import range
from functools import partial
import grpc

from tensorrtserver.api import api_pb2
from tensorrtserver.api import grpc_service_pb2
from tensorrtserver.api import grpc_service_pb2_grpc
import tensorrtserver.api.model_config_pb2 as model_config

from PIL import Image

def model_dtype_to_np(model_dtype):
  if model_dtype == model_config.TYPE_BOOL:
    return np.bool
  elif model_dtype == model_config.TYPE_INT8:
    return np.int8
  elif model_dtype == model_config.TYPE_INT16:
    return np.int16
  elif model_dtype == model_config.TYPE_INT32:
    return np.int32
  elif model_dtype == model_config.TYPE_INT64:
    return np.int64
  elif model_dtype == model_config.TYPE_UINT8:
    return np.uint8
  elif model_dtype == model_config.TYPE_UINT16:
    return np.uint16
  elif model_dtype == model_config.TYPE_FP16:
    return np.float16
  elif model_dtype == model_config.TYPE_FP32:
    return np.float32
  elif model_dtype == model_config.TYPE_FP64:
    return np.float64
  elif model_dtype == model_config.TYPE_STRING:
    return np.dtype(object)
  return None

def parse_model(status, model_name, batch_size, verbose=False):
  """
  Check the configuration of a model to make sure it meets the
  requirements for an image classification network (as expected by
  this client)
  """
  server_status = status.server_status
  if model_name not in server_status.model_status.keys():
    raise Exception("unable to get status for '" + model_name + "'")

  status = server_status.model_status[model_name]
  config = status.config

  if len(config.input) != 1:
    raise Exception("expecting 1 input, got {}".format(len(config.input)))
  if len(config.output) != 1:
    raise Exception("expecting 1 output, got {}".format(len(config.output)))

  input = config.input[0]
  output = config.output[0]

  if output.data_type != model_config.TYPE_FP32:
    raise Exception("expecting output datatype to be TYPE_FP32, model '" +
                    model_name + "' output type is " +
                    model_config.DataType.Name(output.data_type))

  # Output is expected to be a vector. But allow any number of
  # dimensions as long as all but 1 is size 1 (e.g. { 10 }, { 1, 10
  # }, { 10, 1, 1 } are all ok).
  non_one_cnt = 0
  for dim in output.dims:
    if dim > 1:
      non_one_cnt += 1
      if non_one_cnt > 1:
        raise Exception("expecting model output to be a vector")

  # Model specifying maximum batch size of 0 indicates that batching
  # is not supported and so the input tensors do not expect an "N"
  # dimension (and 'batch_size' should be 1 so that only a single
  # image instance is inferred at a time).
  max_batch_size = config.max_batch_size
  if max_batch_size == 0:
    if batch_size != 1:
      raise Exception("batching not supported for model '" + model_name + "'")
  else:  # max_batch_size > 0
    if batch_size > max_batch_size:
      raise Exception(
        "expecting batch size <= {} for model '{}'".format(max_batch_size, model_name))

  # Model input must have 3 dims, either CHW or HWC
  if len(input.dims) != 3:
    raise Exception(
      "expecting input to have 3 dimensions, model '{}' input has {}".format(
        model_name, len(input.dims)))

  if ((input.format != model_config.ModelInput.FORMAT_NCHW) and
      (input.format != model_config.ModelInput.FORMAT_NHWC)):
    raise Exception("unexpected input format " + model_config.ModelInput.Format.Name(input.format) +
                    ", expecting " +
                    model_config.ModelInput.Format.Name(model_config.ModelInput.FORMAT_NCHW) +
                    " or " +
                    model_config.ModelInput.Format.Name(model_config.ModelInput.FORMAT_NHWC))

  if input.format == model_config.ModelInput.FORMAT_NHWC:
    h = input.dims[0]
    w = input.dims[1]
    c = input.dims[2]
  else:
    c = input.dims[0]
    h = input.dims[1]
    w = input.dims[2]

  return (input.name, output.name, c, h, w, input.format, model_dtype_to_np(input.data_type))

def preprocess(img, format, dtype, c, h, w):
  """
  Pre-process an image to meet the size, type and format
  requirements specified by the parameters.
  """
  # np.set_printoptions(threshold='nan')

  if c == 1:
    sample_img = img.convert('L')
  else:
    sample_img = img.convert('RGB')

  resized_img = sample_img.resize((w, h), Image.BILINEAR)
  resized = np.array(resized_img)
  if resized.ndim == 2:
    resized = resized[:, :, np.newaxis]

  typed = resized.astype(dtype)

  scaled = (typed / 255) - 0.5

  # Channels are in RGB order. Currently model configuration data
  # doesn't provide any information as to other channel orderings
  # (like BGR) so we just assume RGB.
  return scaled

def postprocess(results, filenames, batch_size):
  """
  Post-process results to show classifications.
  """
  if len(results) != 1:
    raise Exception("expected 1 result, got {}".format(len(results)))

  batched_result = results[0].batch_classes
  if len(batched_result) != batch_size:
    raise Exception("expected {} results, got {}".format(batch_size, len(batched_result)))
  if len(filenames) != batch_size:
    raise Exception("expected {} filenames, got {}".format(batch_size, len(filenames)))

  label, score = [], []
  # batch size is always 1 here, need to modify if were to larger batch_size
  for (index, result) in enumerate(batched_result):
    print("Image '{}':".format(filenames[index]))
    for cls in result.cls:
      label.append(cls.label)
      score += [{"index": cls.label, "val": cls.value}]
      print("    {} ({}) = {}".format(cls.idx, cls.label, cls.value))
  return label[0], score

def requestGenerator(input_name, output_name, c, h, w, format, dtype, model_name, model_version, image_filename,
                     result_filenames):
  # Prepare request for Infer gRPC
  # The meta data part can be reused across requests
  request = grpc_service_pb2.InferRequest()
  request.model_name = model_name
  if model_version is None:
    request.model_version = -1
  else:
    request.model_version = model_version
  # optional pass in a batch size for generate requester over a set of image files, need to refactor
  batch_size = 1
  request.meta_data.batch_size = batch_size
  output_message = api_pb2.InferRequestHeader.Output()
  output_message.name = output_name
  # Number of class results to report. Default is 10 to match with demo.
  output_message.cls.count = 10
  request.meta_data.output.extend([output_message])

  filenames = []
  if os.path.isdir(image_filename):
    filenames = [os.path.join(image_filename, f)
                 for f in os.listdir(image_filename)
                 if os.path.isfile(os.path.join(image_filename, f))]
  else:
    filenames = [image_filename, ]

  filenames.sort()

  # Preprocess the images into input data according to model
  # requirements
  image_data = []
  for filename in filenames:
    img = Image.open(filename)
    image_data.append(preprocess(img, format, dtype, c, h, w))

  request.meta_data.input.add(name=input_name)

  # Send requests of batch_size images. If the number of
  # images isn't an exact multiple of batch_size then just
  # start over with the first images until the batch is filled.
  image_idx = 0
  last_request = False
  while not last_request:
    input_bytes = None
    input_filenames = []
    del request.raw_input[:]
    for idx in range(batch_size):
      input_filenames.append(filenames[image_idx])
      if input_bytes is None:
        input_bytes = image_data[image_idx].tobytes()
      else:
        input_bytes += image_data[image_idx].tobytes()

      image_idx = (image_idx + 1) % len(image_data)
      if image_idx == 0:
        last_request = True

    request.raw_input.extend([input_bytes])
    result_filenames.append(input_filenames)
    yield request

def get_prediction(image_filename, server_host='localhost', server_port=8001,
                   model_name="bolt", model_version=None):
  """
  Retrieve a prediction from a TensorFlow model server

  :param image:       a bolt image
  :param server_host: the address of the TensorRT inference server
  :param server_port: the port used by the server
  :param model_name: the name of the model
  :param timeout:     the amount of time to wait for a prediction to complete
  :return 0:          the integer predicted in the bolt image
  :return 1:          the confidence scores for all classes
  """
  channel = grpc.insecure_channel(server_host + ':' + str(server_port))
  grpc_stub = grpc_service_pb2_grpc.GRPCServiceStub(channel)

  # Prepare request for Status gRPC
  request = grpc_service_pb2.StatusRequest(model_name=model_name)
  # Call and receive response from Status gRPC
  response = grpc_stub.Status(request)
  # Make sure the model matches our requirements, and get some
  # properties of the model that we need for preprocessing
  batch_size = 1
  verbose = False
  input_name, output_name, c, h, w, format, dtype = parse_model(
    response, model_name, batch_size, verbose)

  filledRequestGenerator = partial(requestGenerator, input_name, output_name, c, h, w, format, dtype, model_name,
                                   model_version, image_filename)

  # Send requests of batch_size images. If the number of
  # images isn't an exact multiple of batch_size then just
  # start over with the first images until the batch is filled.
  result_filenames = []
  requests = []
  responses = []

  # Send request
  for request in filledRequestGenerator(result_filenames):
    responses.append(grpc_stub.Infer(request))

  # For async, retrieve results according to the send order
  for request in requests:
    responses.append(request.result())

  idx = 0
  for response in responses:
    print("Request {}, batch size {}".format(idx, batch_size))
    label, score = postprocess(response.meta_data.output, result_filenames[idx], batch_size)
    idx += 1

  return label, score

def random_bolt(img_path='/workspace/web_server/static/small_bolt'):
  """
  Pull a random image out of the small bolt dataset

  :param savePath: the path to save the file to. If None, file is not saved
  :return 0: file selected
  :return 1: label selelcted
  """
  random_dir = random.choice(os.listdir(img_path))
  random_file = random.choice(os.listdir(img_path + '/' + random_dir))

  return img_path + '/' + random_dir + '/' + random_file, random_dir, 'static/small_bolt' + '/' + random_dir + '/' + random_file

max_stars_repo_name: mengdong/code-snippets

max_issues_repo_name: mengdong/code-snippets

max_forks_repo_name: mengdong/code-snippets


text:
from common import small_buffer

import numpy as np

import vaex

def test_unique_arrow():
    ds = vaex.from_arrays(x=vaex.string_column(['a', 'b', 'a', 'a', 'a', 'b', 'b', 'b', 'b', 'a']))
    with small_buffer(ds, 2):
        assert set(ds.unique(ds.x).to_pylist()) == {'a', 'b'}
        values, index = ds.unique(ds.x, return_inverse=True)
        assert np.array(values)[index].tolist() == ds.x.tolist()

def test_unique():
    ds = vaex.from_arrays(colors=['red', 'green', 'blue', 'green'])
    with small_buffer(ds, 2):
        assert set(ds.unique(ds.colors).to_pylist()) == {'red', 'green', 'blue'}
        values, index = ds.unique(ds.colors, return_inverse=True)
        assert np.array(values)[index].tolist() == ds.colors.tolist()

    ds = vaex.from_arrays(x=['a', 'b', 'a', 'a', 'a', 'b', 'b', 'b', 'b', 'a'])
    with small_buffer(ds, 2):
        assert set(ds.unique(ds.x).to_pylist()) == {'a', 'b'}
        values, index = ds.unique(ds.x, return_inverse=True)
        assert np.array(values)[index].tolist() == ds.x.tolist()

def test_unique_f4():
    x = np.array([np.nan, 0, 1, np.nan, 2, np.nan], dtype='f4')
    df = vaex.from_arrays(x=x)
    assert list(sorted(df.x.unique()))[1:] == [np.nan, 0, 1, 2][1:]

def test_unique_nan():
    x = [np.nan, 0, 1, np.nan, 2, np.nan]
    df = vaex.from_arrays(x=x)
    assert list(sorted(df.x.unique()))[1:] == [np.nan, 0, 1, 2][1:]
    with small_buffer(df, 2):
        values, indices = df.unique(df.x, return_inverse=True)
        values = values[indices]
        mask = np.isnan(values)
        assert values[~mask].tolist() == df.x.values[~mask].tolist()
        # assert indices.tolist() == [0, 1, 2, 0, 3, 0]

def test_unique_missing():
    # Create test databn
    x = np.array([None, 'A', 'B', -1, 0, 2, '', '', None, None, None, np.nan, np.nan, np.nan, np.nan])
    df = vaex.from_arrays(x=x)
    uniques = df.x.unique(dropnan=True).tolist()
    assert set(uniques) == set(['', 'A', 'B', -1, 0, 2, None])

def test_unique_string_missing():
    x = ['John', None, 'Sally', None, '0.0']
    df = vaex.from_arrays(x=x)
    result = df.x.unique().to_pylist()

    assert len(result) == 4
    assert'John' in result
    assert None in result
    assert 'Sally'


text:
import json
from re import S
import matplotlib.pyplot as plt
import numpy as np

def mainstatScript(userInputFileName):
    with open(userInputFileName, 'r') as JsonFile:
        data = json.load(JsonFile)
    amountofComputers = len(data.keys())
    #print("Amount of computers: {}".format(amountofComputers))
    #[dictlist.extend(k,v) for k,v in data.items()]
    Occurances = {}
    portStorage = []
    piLabelKeys = []
    piValues = []
    def func(pct, allvals):
        absolute = int(pct/100.*np.sum(allvals))
        return "{:.1f}%\n{:d} open".format(pct, absolute)
    def countCurrencies(storage, value):
        try:
            storage[value] = storage[value] + 1
        except KeyError as e:
            storage[value] = 1
        return
    for k,v in data.items():
        for port in v['tcp'].keys():

            portStorage.append(port)
    #print(portStorage)
    for data in portStorage:
        countCurrencies(Occurances, data)
    #print(Occurances)
    #for k, v in Occurances.items():
    #    print("Key " + k + " has occurred "  + str(v) + " times")
    portNumber = list(Occurances.keys())
    amountOfTimes = list(Occurances.values())
    def bar():
        c = ['red', 'black', 'gray', 'blue', 'green']
        plt.bar(range(len(Occurances)), amountOfTimes, tick_label=portNumber, color=c)
        plt.xlabel('Port Number', fontweight='bold', color = 'blue', fontsize='12', horizontalalignment='center')
        plt.ylabel('Number Amount Open', fontweight='bold', color = 'blue', fontsize='12', horizontalalignment='center')
        plt.suptitle("Ports Open on Network - Bar\n\n")
        plt.savefig('bar.png')

    def piChart():
        for k,v in Occurances.items():
            piLabelKeys.append(k)
            piValues.append(v)
        fig, ax = plt.subplots(figsize=(6, 5), subplot_kw=dict(aspect="equal"))
        wedges, texts, autotexts = ax.pie(piValues, labels= piLabelKeys, autopct=lambda pct: func(pct, piValues),textprops={'fontsize': 11},shadow=True,radius=1)
        plt.setp(autotexts, size=9, weight="bold")
        ax.set_title("Ports Open on Network - Pie\n")
        plt.axis('equal')
        plt.savefig('pichart.png')
        plt.close()

    piChart()
    bar()
    return amountofComputers


text:
import numpy as np
import pandas as pd
from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt
from functools import partial

np.random.seed(42)

def calc_dist(p1, p2):
    return np.sqrt((p1[0] - p2[0])**2 + (p1[1] - p2[1])**2 + (p1[2] - p2[2])**2)

def make_seg_sphere(segments, center, radius=50):
    #o = [center["Coord X"].iloc[0], center["Coord Y"].iloc[0], center["Coord Z"].iloc[0]]
    o = [center["Coord X"], center["Coord Y"], center["Coord Z"]]
    return segments[calc_dist(o, [segments["Coord X"], segments["Coord Y"], segments["Coord Z"]]) <= radius]

# class SegmentList:
#     def __init__(self, file):
#         segs = pd.read_csv(file)

#     def filter(self, part, func, center):
#         return func(self.segs[self.segs["Type"] == part], center)

class FunctionalGroup:
    def __init__(self, seg_list, center, num_cells, num_clusters, name, start_id, grouping_func, clustering_func):
        self.center = center
        self.n_cells = num_cells
        self.n_clusters = num_clusters
        self.name =  name
        self.start_id = start_id
        #self.grouping_func = grouping_func
        self.clustering_func = clustering_func
        self.group_segs = grouping_func(seg_list, center)

        self.cells = []
        self.clusters = []

        self.make_clusters()
        self.create_cells()

    #Initializes the clusters of the group.
    #Centers are randomly selected (uniquely) from all possible segs.
    #Might want some purposeful spread?
    def make_clusters(self):
        if self.n_clusters > len(self.group_segs):
            raise Exception("The number of segment options should never be bellow n_clusters.")

        center_ids = np.random.choice(range(0, len(self.group_segs)), self.n_clusters, replace=False)
        for id in center_ids:
            self.clusters.append(Cluster(self.group_segs.iloc[id], self.clustering_func, self.group_segs))

    #Creates the group's Cells with randomly assigned number of synapses.
    #Distributes the cell's synapses onto clusters.
    def create_cells(self):
        for i in range(self.n_cells):
            cell = Cell(np.random.randint(low=2, high=8 + 1))
            cluster_ids = np.random.choice(range(0, self.n_clusters), cell.n_syns, replace=False)
            for id in cluster_ids:
                #import pdb; pdb.set_trace()
                cell.add_syn(self.clusters[id].random_seg())

            self.cells.append(cell)

class Cluster:
    def __init__(self, center, clustering_func, group_segs):
        self.center = center
        #import pdb; pdb.set_trace()
        self.n_syns = 0
        self.cluster_segs = clustering_func(group_segs, center)

    #Returns a random segment from cluster segs.
    #Could add count table that spreads synapses.
    def random_seg(self):
        #import pdb; pdb.set_trace()
        self.n_syns += 1
        return self.cluster_segs.iloc[np.random.choice(len(self.cluster_segs))]

#Simple class that contains the number of synapses and the synapse locations for an input cell.
class Cell:
    def __init__(self, n_syns):
        self.n_syns = n_syns
        self.syn_segs = [] #List of Segments where the cell's synapses should be placed.
        self.num_set = 0#Number of synapses already created in bmtk.

    def add_syn(self, seg):
        if len(self.syn_segs) >= self.n_syns:
            raise Exception("Error: too many synapses added to cell.")

        self.syn_segs.append(Segment(seg["BMTK ID"], seg["X"]))

    #Returns the next syn location.
    def get_seg(self):
        if self.num_set >= self.n_syns:
            raise Exception("Too many synapses created for cell.")

        result = self.syn_segs[self.num_set]
        self.num_set += 1
        return result

#Simple class that contains the BMTK morphology id and x of a segment.
class Segment:
    def __init__(self, bmtk_id, x):
        self.bmtk_id = bmtk_id
        self.x = x

def plot_group(segs, group, ax, psoma=[0, 0, 0]):
    #import pdb; pdb.set_trace()
    ax.scatter3D([group.center["Coord X"] + psoma[0]], [group.center["Coord Y"] + psoma[1]], [group.center["Coord Z"] + psoma[2]], color="Black", marker="*", s=50)
    zdata = []
    xdata = []
    ydata = []
    vdata = []

    cell_colors = np.linspace(0, 1, group.n_cells)
    for i, cell in enumerate(group.cells):
        for seg in cell.syn_segs:
            #import pdb; pdb.set_trace()
            s = segs[(segs["X"] == seg.x) & (segs["BMTK ID"] == seg.bmtk_id)].iloc[0]
            xdata.append(s["Coord X"])
            ydata.append(s["Coord Y"])
            zdata.append(s["Coord Z"])
            vdata.append(cell_colors[i])

    ax.scatter3D(np.array(xdata) + psoma[0], np.array(ydata) + psoma[1], np.array(zdata) + psoma[2], c=vdata, cmap="hsv", alpha = 1, s=10)

    zdata = []
    xdata = []
    ydata = []

    for cluster in group.clusters:
        #import pdb; pdb.set_trace()
        #s = df[(df["X"] == seg.x) & (df["BMTK ID"] == seg.bmtk_id)].iloc[0]
        xdata.append(cluster.center["Coord X"])
        ydata.append(cluster.center["Coord Y"])
        zdata.append(cluster.center["Coord Z"])

    ax.scatter3D(np.array(xdata) + psoma[0], np.array(ydata) + psoma[1], np.array(zdata) + psoma[2], color = "Green", s=30)

# df = pd.read_csv("Segments_small.csv")
# dends = df[df["Type"] == "dend"]
# apics = df[df["Type"] == "apic"]
# #import pdb; pdb.set_trace()
# group = FunctionalGroup(dends, dends.sample().iloc[0], 60, 20, "first", partial(make_seg_sphere, radius = 100), partial(make_seg_sphere, radius = 5))

# ax = plt.axes(projection='3d')

# plot_group(dends, group, ax)
# plot_group(apics, FunctionalGroup(apics, apics.sample().iloc[0], 60, 20, "first", partial(make_seg_sphere, radius = 100), partial(make_seg_sphere, radius = 5)), ax)

# # ax.set_xlim([-200, 300])
# # ax.set_ylim([-200, 1200])
# # ax.set_zlim([-200, 50])
# plt.show()


text:
#coding:utf8

import json
label_path =   '/mnt/7/zhihu/ieee_zhihu_cup/data/labels.json'
test_data_path='/mnt/7/zhihu/ieee_zhihu_cup/data/test.npz'
def ensamble(file1,file2,label_path=label_path,     test_data_path=test_data_path,result_csv=None):
    import torch as t
    import numpy as np
    if result_csv is None:
        import time
        result_csv = time.strftime('%y%m%d_%H%M%S.csv')
    a = t.load(file1)
    b = t.load(file2)
    r = 9.0*a+b
    result = r.topk(5,1)[1]

    index2qid = np.load(test_data_path)['index2qid'].item()
    with open(label_path) as f:   label2qid = json.load(f)['id2label']
    rows = range(result.size(0))
    for ii,item in enumerate(result):
        rows[ii] = [index2qid[ii]] + [label2qid[str(_)] for _ in item ]
    import csv
    with open(result_csv,'w') as f:
        writer = csv.writer(f)
        writer.writerows(rows)

if __name__ == '__main__':
    import fire
    fire.Fire()

files = ['CNNText_tmp0.4024_char_test.pth',
 'CNNText_tmp0.4024_char_val.pth',
 'DeepText0.4103_word_test.pth',
 'Inception0.4110_word.pth',
 'LSTMText0.4119_word.pth',
 'LSTMText0.4031_char_test.pth',
 'LSTMText0.4119_word_test.pth',
 'DeepText0.4103_word_val.pth',
 'CNNText_tmp0.4109_word_val.pth',
 'LSTMText0.4119_word_val.pth',
 'RCNN_0.4115_word_test.pth',
 'RCNN_0.4037_char_val.pth',
 'LSTMText0.4031_char_val.pth',
 'RCNN_0.4115_word_val.pth',
 'RCNN_0.4037_char_test.pth',
 'CNNText_tmp0.4109_word_test.pth']

max_stars_repo_name: phasorhand/PyTorchText

max_issues_repo_name: LydonL/PyTorchText

max_forks_repo_name: LydonL/PyTorchText


text:
import trimesh
import numpy as np

def cover(polygon, step):
    """
    Cover a polygon's OBB with a uniform grid.

    Parameters
    -----------
    polygon : shapely.geometry.Polygon
      Polygon to be covered
    step : float
      Distance between grid rows

    Returns
    ------------
    path : (n, 2) float
      Points making up the grid path
    """
    if hasattr(polygon, 'exterior'):
        # polygons
        points = np.array(polygon.exterior.coords)
    elif hasattr(polygon, 'coords'):
        # linestrings etc
        points = np.array(polygon.coords)
    else:
        # numpy arrays
        points = np.array(polygon)

    # find the OBB of the points passed
    rotation, extents = trimesh.bounds.oriented_bounds_2D(points)

    ceil = np.ceil(extents / step) * step
    half = ceil / 2.0

    path = trimesh.transform_points(
        create(bounds=[-half, half], step=step),
        matrix=np.linalg.inv(rotation))

    if False:
        import matplotlib.pyplot as plt
        plt.plot(*path.T, color='r', linestyle='-')
        # plt.scatter(*polygon.T)
        plt.show()

    return path

def create(bounds, step):
    """
    Create a grid inside an axis aligned bounding box.

    Parameters
    ----------
    bounds : (2, 2) float
      Axis aligned 2D bounding box
    step : float
      Distance between grid lines
    """
    bounds = np.array(bounds, dtype=np.float64)

    extents = bounds.ptp(axis=0)

    if extents[0] > extents[1]:
        flip = True
        bounds = np.fliplr(bounds)
    else:
        flip = False

    # make sure we've included the whole boundary
    steps = np.ceil(extents / step).astype(int)
    # linearly space between bounds
    sp = np.linspace(*bounds[:, 0], steps[0])

    # do array wangling
    x = (np.tile(sp, (2, 1)).T).ravel()
    y = (np.ones(len(sp) * 2).reshape((len(sp), 2)) * bounds[:, 1])
    # flip every other column
    y[::2] = np.fliplr(y[::2])
    y = y.ravel()

    if flip:
        result = np.column_stack((y, x))
    else:
        result = np.column_stack((x, y))

    return result

max_stars_repo_name: mikedh/pypocketing

max_issues_repo_name: mikedh/pypocketing

max_forks_repo_name: mikedh/pypocketing


text:
#!/usr/bin/env python

import attr
import click
import shutil
from box import Box
from copy import deepcopy
from glob import glob
from pathlib import Path
from argparse import Namespace

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.metrics import roc_curve, auc
from sklearn.model_selection import GroupKFold
from scikitplot.metrics import plot_confusion_matrix

import cv2
import albumentations as A
from albumentations.pytorch.transforms import ToTensorV2

from apex.optimizers import FusedAdam

import torch
from torch import nn
from torch.nn import functional as F
from torch.optim.lr_scheduler import CosineAnnealingWarmRestarts
from torch.utils.data import DataLoader

import pytorch_lightning as pl
from pytorch_lightning.callbacks import ModelCheckpoint
from pytorch_lightning.loggers import NeptuneLogger
from efficientnet_pytorch import EfficientNet
from catalyst.data.sampler import BalanceClassSampler

pl.seed_everything(1337)

ROOT_PATH = Path("/home/konodyuk/alaska2/")
DATASET_PATH = ROOT_PATH / "input/alaska2"
CKPT_PATH = ROOT_PATH / "experiments"
LABELS = ["Cover", "JMiPOD", "JUNIWARD", "UERD"]

class AlaskaDataset():
    def __init__(self, frame, transforms):
        self.frame = frame
        self.transforms = transforms

    def __getitem__(self, idx):
        img = cv2.imread(self.frame.path.values[idx])
        img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
        img = self.transforms(image=img)["image"]
        item = dict()
        item["image"] = img
        item["target"] = self.frame.target.values[idx]
        return item

    def __len__(self):
        return len(self.frame)

    @property
    def labels(self):
        return list(self.frame.target.values)

class AlaskaTTADataset():
    def __init__(self, frame, transforms_list):
        self.frame = frame
        self.transforms_list = transforms_list

    def __getitem__(self, idx):
        frame_idx = idx % len(self.frame)
        tfm_idx = idx // len(self.frame)
        img = cv2.imread(self.frame.path.values[frame_idx])
        img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
        img = self.transforms_list[tfm_idx](image=img)["image"]
        item = dict()
        item["image"] = img
        item["target"] = self.frame.target.values[frame_idx]
        return item

    def __len__(self):
        return len(self.frame) * len(self.transforms_list)

class AlaskaModel(pl.LightningModule):
    def __init__(self, hparams):
        super().__init__()
        self.hparams = hparams
        self.lr = self.hparams.lr
        self.batch_size = self.hparams.batch_size

        self.model = EfficientNet.from_pretrained(
            self.hparams.model_name,
            advprop=self.hparams.model_advprop,
            num_classes=len(LABELS),
        )

        self.criterion = _LabelSmoothing(smoothing=self.hparams.smoothing)
        self.train_aucs = [_MetricAccumulator(
            _alaska_weighted_auc,
            period=self.hparams.metric_period,
            default_values=([0, 1], [0.5, 0.5]),
            window=window,
        ) for window in [10000, 100000]]
        self.valid_auc = _MetricAccumulator(
            _alaska_weighted_auc,
            period=-1,
            default_values=([0, 1], [0.5, 0.5]),
        )
        self.train_confusion = _MetricAccumulator(_confusion_matrix, period=-1)
        self.valid_confusion = _MetricAccumulator(_confusion_matrix, period=-1)

    def forward(self, images):
        return self.model(images)

    def training_step(self, batch, batch_idx):
        images = batch["image"]
        targets = batch["target"]
        logits = self.forward(images)
        loss = self.criterion(logits, targets)

        probas = F.softmax(logits, dim=1).cpu().detach().numpy()
        targets = targets.cpu().detach().numpy()

        y_pred = 1 - probas[:, 0]
        y_true = 1 - (targets.flatten() == 0).astype(int)
        for train_auc in self.train_aucs:
            train_auc.update(y_true, y_pred)

        multiclass_pred = probas.argmax(axis=1)
        multiclass_true = targets.flatten()
        self.train_confusion.update(multiclass_true, multiclass_pred)

        auc_logs = {
            f"train_auc@{train_auc.window}": train_auc.result
            for train_auc in self.train_aucs
        }

        logs = {
            "train_loss": loss,
            "lr": self.cur_lr(),
            **auc_logs,
        }

        return {"loss": loss, "log": logs}

    def cur_lr(self):
        for param_group in self.optimizer.param_groups:
            return param_group['lr']

    def validation_step(self, batch, batch_idx):
        images = batch["image"]
        targets = batch["target"]
        logits = self.forward(images)
        loss = self.criterion(logits, targets)

        probas = F.softmax(logits, dim=1).cpu().detach().numpy()
        targets = targets.cpu().detach().numpy()

        y_pred = 1 - probas[:, 0]
        y_true = 1 - (targets.flatten() == 0).astype(int)
        self.valid_auc.update(y_true, y_pred)

        multiclass_pred = probas.argmax(axis=1)
        multiclass_true = targets.flatten()
        self.valid_confusion.update(multiclass_true, multiclass_pred)

        return {"valid_loss": loss}

    def validation_epoch_end(self, outputs):
        valid_loss_mean = torch.stack([x["valid_loss"] for x in outputs]).mean()
        self.valid_auc.recalc()
        self.valid_confusion.recalc()
        self.train_confusion.recalc()

        try:
            self.logger.experiment.log_image("valid_confusion", self.valid_confusion.result)
            self.logger.experiment.log_image("train_confusion", self.train_confusion.result)
        except:
            pass

        logs = {
            "valid_loss": valid_loss_mean,
            "valid_auc": torch.tensor(self.valid_auc.result),
        }

        for train_auc in self.train_aucs:
            train_auc.reset()
        self.valid_auc.reset()
        self.valid_confusion.reset()
        self.train_confusion.reset()
        return {"log": logs}

    def test_step(self, batch, batch_idx):
        images = batch["image"]
        logits = self.forward(images)

        probas = F.softmax(logits, dim=1).cpu().detach().numpy()
        y_pred = 1 - probas[:, 0]

        return {"y_pred": y_pred}

    def test_epoch_end(self, outputs):
        y_pred = np.concatenate([x["y_pred"] for x in outputs])

        if self.tta:
            y_pred = y_pred.reshape((-1, len(self.tta_transforms())), order='F')
            print("TTA corr:")
            print(pd.DataFrame(y_pred).corr())
            y_pred = y_pred.mean(axis=1)

        res = pd.DataFrame()
        res["Id"] = [f"{i + 1}.jpg".zfill(8) for i in range(len(y_pred))]
        res["Label"] = y_pred
        res.to_csv(self.submission_path, index=False)
        return {}

    def configure_optimizers(self):
        self.optimizer = FusedAdam(self.model.parameters(), lr=self.lr)
        scheduler = {
            "scheduler": CosineAnnealingWarmRestarts(
                self.optimizer,
                T_0=self.hparams.optimizer_period,
            ),
            "interval": "step",
        }
        return [self.optimizer], [scheduler]

    def prepare_data(self):
        input_frame = _create_train_frame()
        input_frame["fold"] = _get_folds(input_frame)
        test_frame = _create_test_frame()
        train_frame = input_frame[input_frame.fold != self.hparams.fold]
        valid_frame = input_frame[input_frame.fold == self.hparams.fold]
        self.train_dataset = AlaskaDataset(
            frame=train_frame,
            transforms=self.train_transforms()
        )
        self.valid_dataset = AlaskaDataset(
            frame=valid_frame,
            transforms=self.valid_transforms()
        )
        self.test_dataset = AlaskaDataset(
            frame=test_frame,
            transforms=self.test_transforms()
        )
        if self.tta:
            self.test_dataset = AlaskaTTADataset(
                frame=test_frame,
                transforms_list=self.tta_transforms()
            )

    def train_dataloader(self):
        return DataLoader(
            dataset=self.train_dataset,
            batch_size=self.batch_size,
            num_workers=8,
            sampler=BalanceClassSampler(self.train_dataset.labels),
        )

    def val_dataloader(self):
        return DataLoader(
            dataset=self.valid_dataset,
            batch_size=self.batch_size,
            num_workers=8,
        )

    def test_dataloader(self):
        return DataLoader(
            dataset=self.test_dataset,
            batch_size=self.batch_size,
            num_workers=8,
        )

    def train_transforms(self):
        return A.Compose([
            self.normalizer(),
            A.PadIfNeeded(
                self.hparams.aug_pad_size,
                self.hparams.aug_pad_size,
                border_mode=self.hparams.aug_border_mode,
                p=1.
            ),
            A.VerticalFlip(p=0.5),
            A.RandomRotate90(p=1.),
            ToTensorV2(),
        ])

    def valid_transforms(self):
        return A.Compose([
            self.normalizer(),
            A.PadIfNeeded(
                self.hparams.aug_pad_size,
                self.hparams.aug_pad_size,
                border_mode=self.hparams.aug_border_mode,
                p=1.
            ),
            ToTensorV2(),
        ])

    def test_transforms(self):
        return A.Compose([
            self.normalizer(),
            A.PadIfNeeded(
                self.hparams.aug_pad_size,
                self.hparams.aug_pad_size,
                border_mode=self.hparams.aug_border_mode,
                p=1.
            ),
            ToTensorV2(),
        ])

    def tta_transforms(self):
        tfms = [
            [],
            [A.VerticalFlip(p=1.)],
            [A.HorizontalFlip(p=1.)],
            [A.VerticalFlip(p=1.), A.HorizontalFlip(p=1.)],
            [A.Transpose(p=1.)],
            [A.Transpose(p=1.), A.VerticalFlip(p=1.)],
            [A.Transpose(p=1.), A.HorizontalFlip(p=1.)],
            [A.Transpose(p=1.), A.VerticalFlip(p=1.), A.HorizontalFlip(p=1.)],
        ]
        return [
            A.Compose([
                self.normalizer(),
                A.PadIfNeeded(
                    self.hparams.aug_pad_size,
                    self.hparams.aug_pad_size,
                    border_mode=self.hparams.aug_border_mode,
                    p=1.
                ),
            ] + tfm + [
                ToTensorV2(),
            ])
            for tfm in tfms
        ]

    def normalizer(self):
        if self.hparams.model_advprop:
            return A.Lambda(image=_advprop_norm)
        else:
            return A.Normalize()

# ========== CLI ==========

config_option = click.option(
    "--config",
    "config_path",
    default="./config.yaml",
    show_default=True,
    type=click.Path(
        exists=True,
        file_okay=True,
        dir_okay=False,
        resolve_path=True
    ),
)

ckpt_option = click.option(
    "--ckpt",
    "checkpoint_path",
    show_default=True,
    type=click.Path(
        exists=True,
        file_okay=True,
        dir_okay=False,
        resolve_path=True
    ),
)

@click.group()
def cli():
    pass

@cli.command()
@config_option
def train(config_path):
    config = _get_config(config_path)

    trainer, model = from_config(config, create_loggers=True, copy_sources=True)

    trainer.fit(model)

@cli.command()
@config_option
@ckpt_option
def resume(config_path, checkpoint_path):
    config = _get_config(config_path)

    trainer, model = from_config(config, checkpoint_path=checkpoint_path, create_loggers=True)

    trainer.fit(model)

@cli.command()
@config_option
@ckpt_option
@click.option("--tta", is_flag=True)
def predict(config_path, checkpoint_path, tta):
    config = _get_config(config_path)

    trainer, model = from_config(config, checkpoint_path=checkpoint_path, create_loggers=False, tta=tta)
    model.submission_path = _create_subm_path(checkpoint_path, tta)

    trainer.test(model)

def from_config(config, checkpoint_path=None, create_loggers=True, copy_sources=False, tta=False):
    if copy_sources and not create_loggers:
        print("Can't copy sources without creating loggers")
        return

    hparams = Namespace(**config.hparams.to_dict())
    if checkpoint_path is None:
        model = AlaskaModel(hparams)
    else:
        model = AlaskaModel.load_from_checkpoint(checkpoint_path)

    model.tta = tta

    if create_loggers:
        tracked_sources = [str(Path(__file__).absolute()), config.general.config_path]

        neptune_logger = NeptuneLogger(
            api_key=config.neptune.api_token.strip(),
            project_name=config.neptune.project_name,
            tags=list(config.tags),
            upload_source_files=tracked_sources,
            params=config.hparams.to_dict(),
            upload_stdout=False,
            upload_stderr=False
        )

        neptune_experiment = neptune_logger.experiment
        exp_id = neptune_experiment.id
        config.ckpt_path = CKPT_PATH / exp_id
        joined_tags = _join_tags(config)

        checkpoint_filename = f"{exp_id}_{{epoch:02d}}_{{valid_loss:.2f}}_{{valid_auc:.4f}}_{joined_tags}"
        checkpoint_callback = ModelCheckpoint(
            filepath=f"{config.ckpt_path}/{checkpoint_filename}",
            save_top_k=3,
            verbose=True,
            monitor="valid_auc",
            mode="max",
        )

        if copy_sources:
            for file in tracked_sources:
                shutil.copy(file, config.ckpt_path)

    else:
        neptune_logger = None
        checkpoint_callback = None

    trainer = pl.Trainer(
        resume_from_checkpoint=checkpoint_path,

        accumulate_grad_batches=config.hparams.grad_acc,
        deterministic=True,
        amp_level="O1",

        checkpoint_callback=checkpoint_callback,
        logger=neptune_logger,

        gpus=config.general.gpus,
        default_root_dir=ROOT_PATH,

        gradient_clip_val=0.05,

        overfit_batches=200 if config.general.debug else 0,
    )

    return trainer, model

# ========== UTILS ==========

def _get_config(config_path):
    config = Box.from_yaml(filename=config_path)
    if config.general.debug and "debug" not in config.tags:
        config.tags.append("debug")
    config.general.config_path = config_path
    return config

def _join_tags(config):
    tags = "_".join(config.tags)
    return tags

def _create_subm_path(checkpoint_path, tta=False):
    res = checkpoint_path.replace("ckpt", "csv")
    if tta:
        root, sep, file = res.rpartition("/")
        res = f"{root}{sep}TTA_{file}"
    return res

def _create_train_frame():
    entries = []
    for target, label in enumerate(LABELS):
        for path in (DATASET_PATH / label).glob("*"):
            entries.append(dict(path=str(path), target=target, label=label))
    return pd.DataFrame(entries)

def _create_test_frame():
    entries = []
    for path in (DATASET_PATH / "Test").glob("*"):
        entries.append(dict(path=str(path), target=-1, label="Unknown"))
    return pd.DataFrame(entries).sort_values("path")

def _get_folds(train_frame):
    folds = np.zeros(len(train_frame), dtype=int)
    splitter = GroupKFold(6)
    idx = train_frame['path'].apply(lambda x: int(x.split('/')[-1].split('.')[0])).values
    for fold, (_, idx_val) in enumerate(splitter.split(train_frame.target, train_frame.target, idx)):
        folds[idx_val] = fold
    return folds

def _advprop_norm(image, **kw):
    return image.astype(np.float32) / 255. * 2. - 1.

class _LabelSmoothing(nn.Module):
    def __init__(self, smoothing=0.05):
        super().__init__()
        self.confidence = 1.0 - smoothing
        self.smoothing = smoothing

    def forward(self, x, target):
        if self.training:
            target = self._one_hot(x, target)
            logprobs = F.log_softmax(x, dim=-1)
            nll_loss = -logprobs * target
            nll_loss = nll_loss.sum(-1)
            smooth_loss = -logprobs.mean(dim=-1)
            loss = self.confidence * nll_loss + self.smoothing * smooth_loss
            return loss.mean()
        else:
            return F.cross_entropy(x, target.view(-1))

    def _one_hot(self, x, target):
        result = torch.zeros_like(x, dtype=torch.float)
        result.scatter_(1, target.view(-1, 1), 1)
        return result

@attr.s
class _MetricAccumulator:
    metric = attr.ib()
    period = attr.ib(default=1)
    window = attr.ib(default=-1)
    default_values = attr.ib(default=([], []))
    def __attrs_post_init__(self):
        self.reset()

    def reset(self):
        self.y_true, self.y_pred = deepcopy(self.default_values)
        self.result = -1
        self.step = 0

    def update(self, y_true, y_pred):
        self.y_true += list(y_true)
        self.y_pred += list(y_pred)
        if self.window > 0:
            self.y_true = self.y_true[-self.window:]
            self.y_pred = self.y_pred[-self.window:]
        if self.period > 0 and self.step % self.period == 0:
            self.recalc()
        self.step += 1

    def recalc(self):
        try:
            self.result = self.metric(self.y_true, self.y_pred)
        except:
            self.result = -1

def _alaska_weighted_auc(y_true, y_pred):
    tpr_thresholds = [0.0, 0.4, 1.0]
    weights = [2, 1]

    fpr, tpr, thresholds = roc_curve(y_true, y_pred, pos_label=1)

    areas = np.array(tpr_thresholds[1:]) - np.array(tpr_thresholds[:-1])

    normalization = np.dot(areas, weights)

    competition_metric = 0
    for idx, weight in enumerate(weights):
        y_min = tpr_thresholds[idx]
        y_max = tpr_thresholds[idx + 1]
        mask = (y_min < tpr) & (tpr < y_max)

        x_padding = np.linspace(fpr[mask][-1], 1, 100)

        x = np.concatenate([fpr[mask], x_padding])
        y = np.concatenate([tpr[mask], [y_max] * len(x_padding)])
        y = y - y_min
        score = auc(x, y)
        submetric = score * weight
        best_subscore = (y_max - y_min) * weight
        competition_metric += submetric

    return competition_metric / normalization

def _confusion_matrix(y_true, y_pred):
    fig, ax = plt.subplots(figsize=(16, 12))
    plot_confusion_matrix(y_true, y_pred, ax=ax, cmap="Greens")
    return fig

# ========== MAIN ==========

if __name__ == "__main__":
    cli()


text:
import unittest

import numpy as np

from revgraph.core.values.variable import Variable
from revgraph.core.functions.operations.math.sub import Sub

class SubTestCase(unittest.TestCase):
    def setUp(self) -> None:
        self.a = Variable([[1,2,6], [3,4,5]])
        self.b = Variable([[1,1,1], [1,1,1]])
        self.c = Variable([[1,1,1]])
        self.d = Variable([[1], [2]])
        self.e = Variable(2)

    def test_subtracting_same_shape(self):
        # Forward
        op = Sub(self.a, self.b)
        op.register(op)
        op.new_context()
        expected = np.array([[0,1,5], [2,3,4]])
        actual = op.forward()
        self.assertTrue((expected == actual).all())
        # Backward
        gradient = np.ones((2,3))
        op.accumulate(op, gradient)
        expected1 = np.ones((2,3))
        expected2 = np.full((2,3), -1)
        actual1 = self.a.gradient
        actual2 = self.b.gradient
        self.assertTrue((expected1 == actual1).all())
        self.assertTrue((expected2 == actual2).all())

    def test_subtracting_matrix_with_different_number_of_rows(self):
        # Forward
        op = Sub(self.a, self.c)
        op.register(op)
        op.new_context()
        expected = np.array([[0,1,5], [2,3,4]])
        actual = op.forward()
        self.assertTrue((expected == actual).all())
        # Backward
        gradient = np.ones((2,3))
        op.accumulate(op, gradient)
        actual1 = self.a.gradient
        expected2 = np.full((1,3), -2)
        actual2 = self.c.gradient
        self.assertTrue((gradient == actual1).all())
        self.assertTrue((expected2 == actual2).all())

    def test_subtracting_matrix_with_different_number_of_cols(self):
        # Forward
        op = Sub(self.a, self.d)
        op.register(op)
        op.new_context()
        expected = np.array([[0,1,5], [1,2,3]])
        actual = op.forward()
        self.assertTrue((expected == actual).all())
        # Backward
        gradient = np.ones((2,3))
        op.accumulate(op, gradient)
        expected1 = gradient
        actual1 = self.a.gradient
        expected2 = np.full((2,1), -3)
        actual2 = self.d.gradient
        self.assertTrue((expected1 == actual1).all())
        self.assertTrue((expected2 == actual2).all())

    def test_subtracting_matrix_with_scalar(self):
        # Forward
        op = Sub(self.a, self.e)
        op.register(op)
        op.new_context()
        expected = np.array([[-1,0,4], [1,2,3]])
        actual = op.forward()
        self.assertTrue((expected == actual).all())
        # Backward
        gradient = np.ones((2,3))
        op.accumulate(op, gradient)
        expected1 = gradient
        actual1 = self.a.gradient
        expected2 = np.array(-6)
        actual2 = self.e.gradient
        self.assertTrue((expected1 == actual1).all())
        self.assertTrue((expected2 == actual2).all())

    def test_subtracting_nested_operations_with_repeated_nodes(self):
        # Forward
        subexp1 = Sub(self.a, self.b)
        subexp2 = Sub(self.b, self.c)
        op = Sub(subexp1, subexp2)
        op.register(op)
        op.new_context()
        expected = np.array([[0,1,5], [2,3,4]])
        actual = op.forward()
        self.assertTrue((expected == actual).all())
        # Backward
        gradient = np.ones((2,3))
        op.accumulate(op, gradient)
        expected1 = np.ones((2,3))
        expected2 = np.full((2,3), -2)
        expected3 = np.full((1,3), 2)
        self.assertTrue((expected1 == self.a.gradient).all())
        self.assertTrue((expected2 == self.b.gradient).all())
        self.assertTrue((expected3 == self.c.gradient).all())

    def test_subtracting_with_same_node(self):
        # Forward
        op = Sub(Sub(self.a, self.a), Sub(self.a, self.a))
        op.register(op)
        op.new_context()
        expected = np.zeros((2,3))
        actual = op.forward()
        self.assertTrue((expected == actual).all())
        # Backward
        gradient = np.zeros((2,3))
        op.accumulate(op, gradient)
        expected = np.zeros((2,3))
        self.assertTrue((expected == self.a.gradient).all())

max_stars_repo_name: shhoalex/revgraph

max_issues_repo_name: shinghinho/revgraph

max_forks_repo_name: shinghinho/revgraph


text:
from yolact import Yolact
from utils.augmentations import BaseTransform, FastBaseTransform
from utils.functions import SavePath
from layers.output_utils import postprocess
import pycocotools
from pycocotools.coco import COCO
from pycocotools.cocoeval import COCOeval

from data import cfg, set_cfg

import numpy as np
import torch
import torch.backends.cudnn as cudnn
import argparse
import random
import json
import os
import time
import datetime
from pathlib import Path
from PIL import Image
import glob

import matplotlib.pyplot as plt
import cv2

def str2bool(v):
    if v.lower() in ('yes', 'true', 't', 'y', '1'):
        return True
    elif v.lower() in ('no', 'false', 'f', 'n', '0'):
        return False
    else:
        raise argparse.ArgumentTypeError('Boolean value expected.')

def parse_args(argv=None):
    parser = argparse.ArgumentParser(
        description='YOLACT COCO Evaluation')

    parser.add_argument('--mask_det_file', default='results/mask_detections.json', type=str,
                        help='The output file for coco mask results if --coco_results is set.')
    parser.add_argument('--output_web_json', dest='output_web_json', action='store_true',
                        help='If display is not set, instead of processing IoU values, this dumps detections for usage with the detections viewer web thingy.')
    parser.add_argument('--seed', default=None, type=int,
                        help='The seed to pass into random.seed. Note: this is only really for the shuffle and does not (I think) affect cuda stuff.')
    parser.add_argument('--no_crop', default=False, dest='crop', action='store_false',
                        help='Do not crop output masks with the predicted bounding box.')
    parser.add_argument('--score_threshold', default=0, type=float,
                        help='Detections with a score under this threshold will not be considered. This currently only works in display mode.')
    parser.add_argument('--fast_nms', default=True, type=str2bool,
                        help='Whether to use a faster, but not entirely correct version of NMS.')
    parser.add_argument('--cross_class_nms', default=False, type=str2bool,
                        help='Whether compute NMS cross-class or per-class.')
    parser.add_argument('--trained_model',
                        default='weights/ssd300_mAP_77.43_v2.pth', type=str,
                        help='Trained state_dict file path to open. If "interrupt", this will open the interrupt file.')
    parser.add_argument('--config', default=None,
                        help='The config object to use.')
    parser.add_argument('--cuda', default=True, type=str2bool,
                        help='Use cuda to evaulate model')

    parser.set_defaults(no_bar=False, display=False, resume=False, output_coco_json=False, output_web_json=False, shuffle=False,
                        benchmark=False, no_sort=False, no_hash=False, mask_proto_debug=False, crop=True, detect=False, display_fps=False,
                        emulate_playback=False)

    global args
    args = parser.parse_args(argv)

    if args.output_web_json:
        args.output_coco_json = True

    if args.seed is not None:
        random.seed(args.seed)

class Detections:

    def __init__(self):
        self.mask_data = []

    def add_mask(self, image_id:int, category_id:int, segmentation:np.ndarray, score:float):
        """ The segmentation should be the full mask, the size of the image and with size [h, w]. """
        rle = pycocotools.mask.encode(np.asfortranarray(segmentation.astype(np.uint8)))
        rle['counts'] = rle['counts'].decode('ascii') # json.dump doesn't like bytes strings

        self.mask_data.append({
            'image_id': int(image_id),
            'category_id': int(category_id)+1,
            'segmentation': rle,
            'score': float(score)
        })

    def dump(self):
        with open(args.mask_det_file, 'w') as f:
            json.dump(self.mask_data, f)

def evalimage(net:Yolact, path:str):
    frame = torch.from_numpy(cv2.imread(path)).cuda().float()
    batch = FastBaseTransform()(frame.unsqueeze(0))

    pre_time = time.time()
    preds = net(batch)
    aft_time = time.time()
    timee = 0
    timee = aft_time - pre_time

    return preds, timee

def evalimages(net:Yolact):
    times = 0
    cocoGt = COCO("data/test.json")
    detections = Detections()
    for imgid in cocoGt.imgs:
        name = cocoGt.loadImgs(ids=imgid)[0]['file_name']
        h = cocoGt.loadImgs(ids=imgid)[0]['height']
        w = cocoGt.loadImgs(ids=imgid)[0]['width']
        path = "data/test_images/" + name
        predict, timee = evalimage(net, path)
        times += timee
        classes, scores, boxes, masks = postprocess(predict, w, h, crop_masks=args.crop, score_threshold=args.score_threshold)

        classes = list(classes.cpu().numpy().astype(int))
        scores = list(scores.cpu().numpy().astype(float))
        masks = masks.cpu().numpy()
        boxes = boxes.cpu().numpy()

        if len(classes) > 0: # If any objects are detected in this image
            for i in range(masks.shape[0]): # Loop all instances
                # save information of the instance in a dictionary then append on coco_dt list
                if (boxes[i, 3] - boxes[i, 1]) * (boxes[i, 2] - boxes[i, 0]) > 0:
                    detections.add_mask(imgid, classes[i], masks[i,:,:], scores[i])

    detections.dump()
    print(times/100)
    print('Done.')

def evaluate(net:Yolact):
    net.detect.use_fast_nms = args.fast_nms
    net.detect.use_cross_class_nms = args.cross_class_nms
    cfg.mask_proto_debug = args.mask_proto_debug
    evalimages(net)
    return

if __name__ == '__main__':
    parse_args()

    model_path = SavePath.from_str(args.trained_model)
    args.config = model_path.model_name + '_config'
    print('Config not specified. Parsed %s from the file name.\n' % args.config)
    set_cfg(args.config)

    with torch.no_grad():

        if args.cuda:
            cudnn.fastest = True
            torch.set_default_tensor_type('torch.cuda.FloatTensor')
        else:
            torch.set_default_tensor_type('torch.FloatTensor')

        print('Loading model...', end='')
        net = Yolact()
        net.load_weights(args.trained_model)
        net.eval()
        print(' Done.')

        if args.cuda:
            net = net.cuda()

        evaluate(net)


text:
# This code is part of Qiskit.
#
# (C) Copyright IBM 2017, 2020.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"""Piecewise polynomial Chebyshev approximation to a given f(x)."""

from typing import Callable, List, Optional
import numpy as np
from numpy.polynomial.chebyshev import Chebyshev

from qiskit.circuit import QuantumRegister, AncillaRegister
from qiskit.circuit.library.blueprintcircuit import BlueprintCircuit
from qiskit.circuit.exceptions import CircuitError

from .piecewise_polynomial_pauli_rotations import PiecewisePolynomialPauliRotations

class PiecewiseChebyshev(BlueprintCircuit):
    r"""Piecewise Chebyshev approximation to an input function.

    For a given function :math:`f(x)` and degree :math:`d`, this class implements a piecewise
    polynomial Chebyshev approximation on :math:`n` qubits to :math:`f(x)` on the given intervals.
    All the polynomials in the approximation are of degree :math:`d`.

    The values of the parameters are calculated according to [1].

    Examples:

        .. jupyter-execute::

            import numpy as np
            from qiskit import QuantumCircuit
            from qiskit.circuit.library.arithmetic.piecewise_chebyshev import PiecewiseChebyshev
            f_x, degree, breakpoints, num_state_qubits = lambda x: np.arcsin(1 / x), 2, [2, 4], 2
            pw_approximation = PiecewiseChebyshev(f_x, degree, breakpoints, num_state_qubits)
            pw_approximation._build()
            qc = QuantumCircuit(pw_approximation.num_qubits)
            qc.h(list(range(num_state_qubits)))
            qc.append(pw_approximation.to_instruction(), qc.qubits)
            print(qc.draw(output='mpl'))

    References:

        [1]: Haener, T., Roetteler, M., & Svore, K. M. (2018).
             Optimizing Quantum Circuits for Arithmetic.
             `arXiv:1805.12445 <http://arxiv.org/abs/1805.12445>`_
    """

    def __init__(self,
                 f_x: Callable[[int], float],
                 degree: Optional[int] = None,
                 breakpoints: Optional[List[int]] = None,
                 num_state_qubits: Optional[int] = None,
                 name: str = 'pw_cheb') -> None:
        r"""
        Args:
            f_x: the function to be approximated.
            degree: the degree of the polynomials.
                Defaults to ``1``.
            breakpoints: the breakpoints to define the piecewise-linear function.
                Defaults to the full interval.
            num_state_qubits: number of qubits representing the state.
            name: The name of the circuit object.
        """
        super().__init__(name=name)

        # define internal parameters
        self._num_state_qubits = None

        # Store parameters
        self._f_x = f_x
        self._degree = degree if degree is not None else 1
        self._breakpoints = breakpoints if breakpoints is not None else [0]

        self._polynomials = None

        self.num_state_qubits = num_state_qubits

    def _check_configuration(self, raise_on_failure: bool = True) -> bool:
        valid = True

        if self._f_x is None:
            valid = False
            if raise_on_failure:
                raise AttributeError('The function to be approximated has not been set.')

        if self._degree is None:
            valid = False
            if raise_on_failure:
                raise AttributeError('The degree of the polynomials has not been set.')

        if self._breakpoints is None:
            valid = False
            if raise_on_failure:
                raise AttributeError('The breakpoints have not been set.')

        if self.num_state_qubits is None:
            valid = False
            if raise_on_failure:
                raise AttributeError('The number of qubits has not been set.')

        if self.num_qubits < self.num_state_qubits + 1:
            valid = False
            if raise_on_failure:
                raise CircuitError('Not enough qubits in the circuit, need at least '
                                   '{}.'.format(self.num_state_qubits + 1))

        return valid

    @property
    def f_x(self) -> Callable[[int], float]:
        """The function to be approximated.

                Returns:
                    The function to be approximated.
        """
        return self._f_x

    @f_x.setter
    def f_x(self, f_x: Optional[Callable[[int], float]]) -> None:
        """Set the function to be approximated.

        Note that this may change the underlying quantum register, if the number of state qubits
        changes.

        Args:
            f_x: The new function to be approximated.
        """
        if self._f_x is None or f_x != self._f_x:
            self._invalidate()
            self._f_x = f_x

            self._reset_registers(self.num_state_qubits)

    @property
    def degree(self) -> int:
        """The degree of the polynomials.

                Returns:
                    The degree of the polynomials.
        """
        return self._degree

    @degree.setter
    def degree(self, degree: Optional[int]) -> None:
        """Set the error tolerance.

        Note that this may change the underlying quantum register, if the number of state qubits
        changes.

        Args:
            degree: The new degree.
        """
        if self._degree is None or degree != self._degree:
            self._invalidate()
            self._degree = degree

            self._reset_registers(self.num_state_qubits)

    @property
    def breakpoints(self) -> List[int]:
        """The breakpoints for the piecewise approximation.

                Returns:
                    The breakpoints for the piecewise approximation.
        """
        breakpoints = self._breakpoints

        # it the state qubits are set ensure that the breakpoints match beginning and end
        if self.num_state_qubits is not None:
            num_states = 2 ** self.num_state_qubits

            # If the last breakpoint is < num_states, add the identity polynomial
            if breakpoints[-1] < num_states:
                breakpoints = breakpoints + [num_states]

            # If the first breakpoint is > 0, add the identity polynomial
            if breakpoints[0] > 0:
                breakpoints = [0] + breakpoints

        return breakpoints

    @breakpoints.setter
    def breakpoints(self, breakpoints: Optional[List[int]]) -> None:
        """Set the breakpoints for the piecewise approximation.

        Note that this may change the underlying quantum register, if the number of state qubits
        changes.

        Args:
            breakpoints: The new breakpoints for the piecewise approximation.
        """
        if self._breakpoints is None or breakpoints != self._breakpoints:
            self._invalidate()
            self._breakpoints = breakpoints

            self._reset_registers(self.num_state_qubits)

    @property
    def polynomials(self) -> List[List[float]]:
        """The polynomials for the piecewise approximation.

                Returns:
                    The polynomials for the piecewise approximation.
        """
        if self.num_state_qubits is None:
            return [[]]

        # note this must be the private attribute since we handle missing breakpoints at
        # 0 and 2 ^ num_qubits here (e.g. if the function we approximate is not defined at 0
        # and the user takes that into account we just add an identity)
        num_intervals = len(self._breakpoints)

        # Calculate the polynomials
        polynomials = []
        for i in range(0, num_intervals - 1):
            # Calculate the polynomial approximating the function on the current interval
            poly = Chebyshev.interpolate(self._f_x, self._degree,
                                         domain=[self._breakpoints[i],
                                                 self._breakpoints[i + 1]])
            # Convert polynomial to the standard basis and rescale it for the rotation gates
            poly = 2 * poly.convert(kind=np.polynomial.Polynomial).coef
            # Convert to list and append
            polynomials.append(poly.tolist())

        # If the last breakpoint is < 2 ** num_qubits, add the identity polynomial
        if self._breakpoints[-1] < 2 ** self.num_state_qubits:
            polynomials = polynomials + [[2 * np.arcsin(1)]]

        # If the first breakpoint is > 0, add the identity polynomial
        if self._breakpoints[0] > 0:
            polynomials = [[2 * np.arcsin(1)]] + polynomials

        return polynomials

    @polynomials.setter
    def polynomials(self, polynomials: Optional[List[List[float]]]) -> None:
        """Set the polynomials for the piecewise approximation.

        Note that this may change the underlying quantum register, if the number of state qubits
        changes.

        Args:
            polynomials: The new breakpoints for the piecewise approximation.
        """
        if self._polynomials is None or polynomials != self._polynomials:
            self._invalidate()
            self._polynomials = polynomials

            self._reset_registers(self.num_state_qubits)

    @property
    def num_state_qubits(self) -> int:
        r"""The number of state qubits representing the state :math:`|x\rangle`.

        Returns:
            The number of state qubits.
        """
        return self._num_state_qubits

    @num_state_qubits.setter
    def num_state_qubits(self, num_state_qubits: Optional[int]) -> None:
        """Set the number of state qubits.

        Note that this may change the underlying quantum register, if the number of state qubits
        changes.

        Args:
            num_state_qubits: The new number of qubits.
        """
        if self._num_state_qubits is None or num_state_qubits != self._num_state_qubits:
            self._invalidate()
            self._num_state_qubits = num_state_qubits

            # Set breakpoints if they haven't been set
            if num_state_qubits is not None and self._breakpoints is None:
                self.breakpoints = [0, 2 ** num_state_qubits]

            self._reset_registers(num_state_qubits)

    def _reset_registers(self, num_state_qubits: Optional[int]) -> None:
        if num_state_qubits is not None:
            qr_state = QuantumRegister(num_state_qubits, 'state')
            qr_target = QuantumRegister(1, 'target')
            self.qregs = [qr_state, qr_target]
            self._ancillas = []
            self._qubits = qr_state[:] + qr_target[:]
            self._qubit_set = set(self._qubits)

            num_ancillas = num_state_qubits
            if num_ancillas > 0:
                qr_ancilla = AncillaRegister(num_ancillas)
                self.add_register(qr_ancilla)

        else:
            self.qregs = []
            self._qubits = []
            self._qubit_set = set()
            self._ancillas = []

    def _build(self):
        """Build the circuit. The operation is considered successful when q_objective is :math:`|1>`
        """
        # do not build the circuit if _data is already populated
        if self._data is not None:
            return

        self._data = []

        # check whether the configuration is valid
        self._check_configuration()

        poly_r = PiecewisePolynomialPauliRotations(self.num_state_qubits,
                                                   self.breakpoints, self.polynomials)

        qr_state = self.qubits[:self.num_state_qubits]
        qr_target = [self.qubits[self.num_state_qubits]]
        qr_ancillas = self.qubits[self.num_state_qubits + 1:]

        # Apply polynomial approximation
        self.append(poly_r.to_instruction(), qr_state[:] + qr_target + qr_ancillas[:])

max_stars_repo_name: ewinston/qiskit-sdk-py

max_issues_repo_name: ewinston/qiskit-sdk-py

max_forks_repo_name: ewinston/qiskit-sdk-py

reason: import numpy,from numpy


text:
"""
This file is part of PUQ
Copyright (c) 2013-2106 PUQ Authors
See LICENSE file for terms.
"""
from __future__ import absolute_import, division, print_function

import sys
import puq
import numpy as np
import copy
import pymc
from scipy.stats.kde import gaussian_kde

class Calibrate(object):
    """
    Bayesian calibration of variables.

    Args:
      model(string or function): An expression or function.
      The standard python math functions defined at
      https://docs.python.org/2/library/math.html are implemented.
      'pi' and 'e' are not defined.

      cvars(dict): Dictionary containing calibration variable
        prior PDFs and calibration types.  Example:
        {
            'a': {'prior': 'Uniform(5,20)', 'type': 'D'},
            'b': {'prior': 'Normal(1,100)', 'type': 'S'},
        }

      nvars(dict): A dictionary of variable names
        and data arrays representing observed values.
        Example: {'x': xdata, 'y': ydata} where
        the data arrays are nx2 (two- column) with the first
        column values and the second column deviations (measurement errors).

      outvar(dict): Dictionary with one entry containing the
        observed output data and uncertainty as a two-column array.

    Returns: A Calibrate object. You should use it to call the run method
      to do the actual MCMC calibration.

    -------------------

    Prior PDFs:
      Normal(mean,dev) : For example, "Normal(100,1)"
      Uniform(min,max) : For example, "Uniform(100,120)"

    Calibration Types:
        'D': Deterministic. Calibrate for a fixed value.
        'S': Stochastic. Calibrate the variable mean and deviation.
    """

    # model can be response surface, equation, python function

    def __init__(self, model, cvars, nvars, outvar, MAP='fmin_powell'):

        out_var_name = list(outvar.keys())[0]

        if callable(model):
            # model is a python function
            self.model = model
        else:
            try:
                if model.find('=') != -1:
                    m = model.split('=')
                    if m[0].strip() == out_var_name:
                        model = m[1].strip()
                    elif m[1].strip() == out_var_name:
                        model = m[0].strip()
                    else:
                        raise
                model = sympy.sympify(model, _clash)
            except:
                err = "Model expression must be of form 'varname=expression'\
                and varname must have measured data in the data table."
                raise ValueError(err)

            # find our variable names
            val_names = list(map(str, model.free_symbols))

            # all variables in the model must be defined
            if set(val_names) != set(nvars.keys()).union(set(cvars.keys())):
                missing = list(set(val_names).difference(set(nvars.keys()).union(set(cvars.keys()))))
                if missing == []:
                    errstr = 'Error: Model did not use all parameters'
                else:
                    errstr = "Error: Not all parameters in the model were defined."
                    errstr += "\'%s\' not defined." % missing[0]
                raise ValueError(errstr)

            # turn our symbolic expression into a fast, safe function call
            self.model = lambdify(model.free_symbols, model, dummify=False, modules=['numpy', 'mpmath', 'sympy'])

        out_var = outvar[out_var_name]
        var = {}
        means = {}
        devs = {}
        dlen = out_var.shape[0]
        self.num_samples = 100000
        self.num_burn = 20000
        self.num_thin = 8

        # Calibration variables
        for v in cvars.keys():
            v = str(v)  # pymc doesn't like unicode
            # convert to lowercase and parse
            d = cvars[v]['prior'].lower().replace('(', ',').replace(')', '').split(',')
            d[1] = float(d[1])
            d[2] = float(d[2])
            if cvars[v]['type'] == 'S':
                if d[0] == 'normal':
                    # normal prior with mean d[1] and deviation d[2]
                    means[v] = pymc.Normal(v + '_mean', mu=d[1], tau=1 / d[2] ** 2, value=d[1])
                    values = np.linspace(d[1] - 3 * d[2], d[1] + 3 * d[2], out_var.shape[0])
                    dval = d[2]
                elif d[0] == 'uniform':
                    # uniform prior from d[1] to d[2]
                    means[v] = pymc.Uniform(v + '_mean', d[1], d[2], value=(d[1] + d[2]) / 2.0)
                    values = np.linspace(d[1], d[2], out_var.shape[0])
                    dval = 1
                else:
                    start_val = (d[1] + d[2]) / 2.0
                    values = np.linspace(d[1], d[2], out_var.shape[0])
                    dval = 1
                    means[v] = pymc.Stochastic(name=v + '_mean',
                                               logp=lambda value: -np.log(value),
                                               doc='',
                                               parents={},
                                               value=start_val)

                # Jeffrey prior for deviation
                devs[v] = pymc.Stochastic(name=v + '_dev',
                                          logp=lambda value: -np.log(value),
                                          doc='',
                                          parents={},
                                          value=dval)

                # to get a reliable deviation, we need more samples
                if self.num_samples < 1000000:
                    self.num_samples *= 10
                    self.num_burn *= 10
                    self.num_thin *= 10
                # create a stochastic node for pymc with the mean and
                # dev from above and some initial values to try
                var[v] = pymc.Normal(v, mu=means[v], tau=1.0 / devs[v] ** 2, value=values)
            else:
                if d[0] == 'normal':
                    var[v] = pymc.Normal(v, mu=d[1], tau=1 / d[2] ** 2)
                elif d[0] == 'uniform':
                    var[v] = pymc.Uniform(v, lower=d[1], upper=d[2])
                elif d[0] == 'jeffreys':
                    start_val = (d[1] + d[2]) / 2.0
                    var[v] = pymc.Stochastic(name=v, doc='',
                                             logp=lambda value: -np.log(value),
                                             parents={},
                                             value=start_val)
                else:
                    print('Unknown probability distribution: %s' % d[0])
                    return None
        for v in nvars.keys():
            var[v] = nvars[v][:, 0]

        results = pymc.Deterministic(
            eval=self.model,
            name='results',
            parents=var,
            doc='',
            trace=True,
            verbose=0,
            dtype=float,
            plot=False,
            cache_depth=2)

        mdata = out_var[:, 0]
        mdata_err = out_var[:, 1]

        mcmc_model_out = pymc.Normal('model_out', mu=results, tau=1.0 / mdata_err ** 2, value=mdata, observed=True)
        self.mcmc_model = pymc.Model(list(var.values()) + list(means.values()) + list(devs.values()) + [mcmc_model_out])

        if MAP is not None:
            # compute MAP and use that as start for MCMC
            map_ = pymc.MAP(self.mcmc_model)
            map_.fit(method=MAP)

            print('\nmaximum a posteriori (MAP) using', MAP)
            for v in cvars.keys():
                print('%s=%s' % (v, var[v].value))
            print()

        # NOT calibration variables
        for v in nvars.keys():
            data = nvars[v][:, 0]
            err = nvars[v][:, 1]
            if np.all(err <= 0.0):
                var[v] = data
            else:
                err[err == 0] = 1e-100
                # norm_err = pymc.Normal(v + '_err', mu=0, tau=1.0 / err ** 2)
                # var[v] = data + norm_err
                var[v] = pymc.Normal(v + '_err', mu=data, tau=1.0 / err ** 2)

        results = pymc.Deterministic(
            eval=self.model,
            name='results',
            parents=var,
            doc='',
            trace=True,
            verbose=0,
            dtype=float,
            plot=False,
            cache_depth=2)

        self.mcmc_model = pymc.Model(list(var.values()) + list(means.values()) + list(devs.values()) + [mcmc_model_out])
        self.cvars = cvars
        self.var = var
        self.dlen = dlen
        self.means = means
        self.devs = devs

    def run(self, samples=None, progress=True):
        """
        Perform MCMC calibration. Returns pdfs and MCMC traces.

        Args:
          samples(integer or tuple): A tuple containing the
            number of samples, the number to burn, and the number to thin.
            If samples is an integer, burn will be 20% of the samples and
            thin will be 8.  Default will use between 10000 and 1000000
            samples, depending on the number of stochastic variables
            being calibrated.

          progress(boolean): If True, will display a progress bar.

        Returns(tuple):
          Returns a tuple containing cvars and a pdf.
          cvars is modified to include key 'trace'
          which will be an array.  It will also have a key 'pdf' which
          will be a PDF function.  For GAUSSIAN type, it will also
          include traces 'mtrace' and 'dtrace' and 'jpdf' corresponding
          to the mean and deviation traces and the joint PDF.

        """
        if samples is None:
            num_samples = self.num_samples
            num_burn = self.num_burn
            num_thin = self.num_thin
        else:
            if type(samples) == tuple:
                if len(samples) != 3:
                    raise ValueError("Error: samples should be a number or tuple of length 3.")
                num_samples, num_burn, num_thin = samples
            else:
                num_samples = samples
                num_burn = int(samples * 0.20)
                num_thin = 8

        Calibrate.mcmc = pymc.MCMC(self.mcmc_model)
        Calibrate.mcmc.sample(
            iter=num_samples,
            burn=num_burn,
            thin=num_thin,
            tune_interval=10000,
            tune_throughout=True,
            progress_bar=progress)

        if Calibrate.mcmc is None:
            return None

        for v in self.cvars.keys():
            t = self.var[v].trace[:]
            if len(t.shape) == 2:
                self.cvars[v]['ntraces'] = t.shape[1]
            else:
                self.cvars[v]['ntraces'] = 1
            self.cvars[v]['trace'] = t.ravel()

        for v in self.means.keys():
            self.cvars[v]['mtrace'] = self.means[v].trace[:]
            self.cvars[v]['dtrace'] = self.devs[v].trace[:]

        # collect all the independent variables and compute KDE
        col_count = max([self.cvars[v]['ntraces'] for v in self.cvars])
        for cv in self.cvars.keys():
            if self.cvars[cv]['type'] == 'S':
                data = np.column_stack((self.cvars[cv]['dtrace'], self.cvars[cv]['mtrace']))
                try:
                    self.cvars[cv]['jpdf'] = gaussian_kde(data.T)
                except:
                    self.cvars[cv]['jpdf'] = None
            # multidimensional traces get flattened and others
            # get repeated to match size.
            if self.cvars[cv]['ntraces'] == col_count:
                n = 1
            else:
                n = col_count
            try:
                self.cvars[cv]['pdf'] = gaussian_kde(self.cvars[cv]['trace'].ravel())
            except:
                self.cvars[cv]['pdf'] = None
            self.cvars[cv]['trace'] = self.cvars[cv]['trace'].ravel().repeat(n)

        data = np.column_stack([self.cvars[v]['trace'] for v in sorted(self.cvars.keys())])
        try:
            k = gaussian_kde(data.T)
        except:
            k = None

        return (self.cvars, k)

def calibrate(params, caldata, err, func, num_samples=None):
    """
    Performs Bayesian calibration for a variable or variables.  This
    function provides a convenient interface for use with PUQ.

    Args:
      params (list): PUQ input parameters. Parameters which have caldata
      included are not calibrated. Parameters without caldata must
      have a caltype.

      caldata (list or array): Experimental output values.

      err (float or array of floats): Standard deviation of the measurement
          error in **caldata**.  If this is a scalar, then use the same error
          for every data point.

      func: Response function.

      num_samples:  A tuple containing the number of samples, the number
        to burn, and the number to thin.  If samples is an integer, burn
        will be 20% of the samples and thin will be 8.  Default will use
        between 10000 and 1000000 samples, depending on the number of
        stochastic variables being calibrated.

    Returns:
      A copy of **params** modified with the calibrated variables.
    """
    print("Performing Bayesian Calibration...")

    cvars = {}
    nvars = {}
    ovar = {}

    err = np.broadcast_to(err, caldata.shape)
    ovar['out'] = np.column_stack((caldata, err))

    # For each PUQ parameter
    for p in params:
        p.name = str(p.name)
        if hasattr(p, 'caldata') and p.caldata is not None:
            # noncalibration parameter with measurements and errors
            if len(p.caldata.shape) != 2:
                print('\nWarning: caldata for %s should have two columns.' % p.name)
                print('Column 1 is the data and column 2 is the error.')
                print("Assuming error is zero and continuing.\n")
                err = np.broadcast_to(0, p.caldata.shape)
                p.caldata = np.column_stack((p.caldata, err))
            nvars[p.name] = p.caldata
        else:
            if isinstance(p, puq.NormalParameter):
                prior = 'Normal(%s,%s)' % (p.pdf.mean, p.pdf.dev)
            elif isinstance(p, puq.UniformParameter):
                prior = 'Uniform(%s,%s)' % (p.pdf.range[0], p.pdf.range[1])
            else:
                print('\nWARNING: Only Normal and Uniform priors are currently supported.')
                print('\tContinuing using Uniform\n')
                prior = 'Uniform(%s,%s)' % (p.pdf.range[0], p.pdf.range[1])
            cvars[p.name] = {'prior': prior, 'type': p.caltype}

    c = Calibrate(func, cvars, nvars, ovar)
    cvars, kpdf = c.run(samples=num_samples)

    newparams = copy.copy(params)
    for i, p in enumerate(newparams):
        if hasattr(p, 'caldata') and p.caldata is not None:
            continue

        vals = cvars[p.name]['trace']
        print("Calibrated %s to a PDF with mean=%s and dev=%s" % (p.name, np.mean(vals), np.std(vals)))
        pdf = puq.ExperimentalPDF(vals, fit=True)
        newparams[i] = puq.CustomParameter(newparams[i].name,
                                           newparams[i].description,
                                           pdf=pdf,
                                           use_samples=True)
        try:
            newparams[i].values = copy.copy(p.values)
        except:
            pass
        newparams[i].trace = vals
        newparams[i].original_parameter = p

    return newparams, kpdf

max_stars_repo_name: c-PRIMED/puq

max_issues_repo_name: c-PRIMED/puq

max_forks_repo_name: c-PRIMED/puq

reason: import numpy,from scipy