text:
"""Group by values of a dependent variable."""

from microbepy.common import constants as cn
from microbepy.common.dataframe_sorter import DataframeSorter
from microbepy.data.model_data_provider import ModelDataProvider
from microbepy.common.study_context import StudyContext
from microbepy.common.range_constraint   \
    import RangeConstraintVector, RangeConstraint
from microbepy.common import util

import copy
import matplotlib.pyplot as plt
import numpy as np
import os
import pandas as pd
from sklearn.cluster import KMeans

OUTPUT_DIRECTORY = util.getDataModelPath(None)
DEFAULT_RC = RangeConstraint(lower=-np.inf, upper=np.inf)
DEFAULT_RC_VECTOR = RangeConstraintVector({
    cn.RATE: DEFAULT_RC,
    cn.YIELD: DEFAULT_RC,
    })

class IsolateMutationStatistics(object):
  """
  Information about isolates and mutations.
  """

  def __init__(self, provider):
    self._provider = provider
    #
    self.mutations = None  # list-str
    self.isolate_pairs = None  # list-str
    self.ser_count = None  # pd.Series

  def do(self):
    self.mutations = self._provider.df_X.columns.tolist()
    self.isolate_pairs = self._provider.df_X.index.values
    self.ser_count = self._provider.df_X.sum()

class DepvarGroup(object):
  """
  Create and manage groups formed by being
  in the same range of values of a dependent variable.
  Comparisons are made between isolates found within
  a constrained range and those that are unconstrained.
  """

  def __init__(self,
      depvar,
      mutation_column=cn.GGENE_ID,
      constraints=None,
      is_plot=True,
      output_directory=OUTPUT_DIRECTORY):
    """
    :param str depvar:
    :param str mutation_column:
    :param ModelDataProvider provider:
    :param list-BooleanFunction constraints:
    :param str output_directory:
    """
    self._depvar = depvar
    self._mutation_column = mutation_column
    self._constraints = constraints
    self._output_directory = output_directory
    self._is_plot = is_plot
    #
    self._provider = self.getData()  # Unconstrained by depvar
    self._statistics = IsolateMutationStatistics(self._provider)
    self._statistics.do()

  def getData(self, rc_vector=None):
    context = StudyContext(depvar=self._depvar,
        mutation_column=self._mutation_column)
    provider = ModelDataProvider(context,
        constraints=self._constraints,
        rc_vector=rc_vector)
    provider.do()
    return provider

  def makeStatisticsDF(self, lower, upper, steps):
    """
    Makes a statistics dataframe.
    :param float lower: lower end of range in units of std
    :param float upper: upper end of range in units of std
    :param float steps: number of steps between lower and upper
    :return pd.DataFrame:
        cn.LOWER: lower part of range of dependent variable
        cn.UPPER: upper part of range of dependent variable
        cn.CENTER: center of the range
        cn.ISOLATES: list of isolates in range
        cn.MUTATIONS: list of mutations in the range
        cn.FRACTION: Series with fraction of isolates with each mutation
    """
    step_size = (upper - lower)/(1.0*steps)
    result_dict = {
        cn.LOWER: [],
        cn.UPPER: [],
        cn.CENTER: [],
        cn.ISOLATES: [],
        cn.MUTATIONS: [],
        cn.FRACTION: [],
        }
    for step in range(steps):
      rc_vector = copy.deepcopy(DEFAULT_RC_VECTOR)
      new_lower = (lower + step_size*step)
      new_upper = new_lower + step_size
      rc_vector.update(self._depvar, RangeConstraint(
          lower=new_lower,
          upper=new_upper,
          ))
      provider = self.getData(rc_vector=rc_vector)
      statistics = IsolateMutationStatistics(provider)
      statistics.do()
      #
      result_dict[cn.LOWER].append(new_lower)
      result_dict[cn.UPPER].append(new_upper)
      result_dict[cn.CENTER].append((new_upper+new_lower)/2.0)
      result_dict[cn.ISOLATES].append(statistics.isolate_pairs)
      result_dict[cn.MUTATIONS].append(statistics.mutations)
      ser = statistics.ser_count / (1.0*len(statistics.isolate_pairs))
      result_dict[cn.FRACTION].append(ser)
    df_result = pd.DataFrame(result_dict)
    return df_result

  def makeHeatmap(self, steps, lower=-3, upper=3):
    """
    Creates a heatmap with depvar as x-axis and mutations as y.
    :param int steps:
    :param float lower: lower value for depvar
    :param float upper: upper value for depvar
    """
    max_major_labels = 5
    df_statistics = self.makeStatisticsDF(lower, upper, steps)
    # Construct a dataframe of fractions for each value of depvar
    df_plot = pd.DataFrame([r[cn.FRACTION] for _, r in
        df_statistics.iterrows()])
    xvalues = df_statistics[cn.LOWER].tolist()
    df_plot.index = xvalues
    df_plot = df_plot.transpose()
    sorter = DataframeSorter(df_plot)
    df_plot = sorter.orderRows()
    # Plot
    fig = plt.figure(figsize=(16, 14))
    for idx, species in  \
        enumerate([cn.SPECIES_MIX_DVH, cn.SPECIES_MIX_MMP]):
      df_plot_species = df_plot.copy()
      ax = fig.add_subplot(1, 2, idx+1)
      mutations = [i for i in df_plot_species.index.values
          if i[0] == species]
      df_plot_species = df_plot_species.loc[mutations]
      plot = ax.pcolor(df_plot_species, cmap='jet')
      ax.set_yticks(np.arange(0.5, len(mutations), 1))
      ax.set_yticklabels(mutations, fontdict={'fontsize': 7})
      #
      num_minor_labels = max(1, len(xvalues)/max_major_labels)
      xlabels = [str(x) if i % num_minor_labels == 0 else ""
           for i, x in enumerate(xvalues)]
      ax.set_xticks(np.arange(0.5, len(xvalues), 1))
      ax.set_xticklabels(xlabels)
      ax.set_xlabel(self._depvar)
      ax.set_title(cn.SPECIES_DICT[species])
      if idx == 1:
        fig.colorbar(plot, cmap='jet')
      else:
        ax.set_ylabel("Mutations")
    if self._is_plot:
      plt.show()

max_stars_repo_name: ScienceStacks/MicrobEPy

max_issues_repo_name: ScienceStacks/MicrobEPy

max_forks_repo_name: ScienceStacks/MicrobEPy


text:
import numpy as np
from skimage.filters import threshold_otsu
from skimage.measure import label as label_components

from ..utils.dataUtils import getCoordsWithinSphere, getCoordsOutsideSphere, morphDilation, morphErosion

def morphologicalDilation(img, dilate_size):
  return morphDilation(img, size=dilate_size)

def removeSmallCCs(img, estimated_background_thr=None, min_size_to_preserve=0.1, dilate_size=3):

  if estimated_background_thr is None:
    estimated_background, thr= estimateBackground(img)
    print( estimated_background, thr)
    assert estimated_background < thr, "Error, background should have smaller intensities than signal. Automatic cleaning failed"

  minVal= np.min(img)
  bin_img= morphDilation(img, size=dilate_size)
  bin_img= morphErosion(bin_img, size=dilate_size)

  bin_img[bin_img<thr]=0
  bin_img[bin_img >=thr] = 1
  bin_img= bin_img.astype(bool)
  print("Post-processing step: ", end=".")

  labels, n_labels= label_components(bin_img, return_num=True, background=estimated_background)

  print(".", end="")
  labels_and_num=getLabels2NumVoxels(labels, n_labels)
  print(".", end="")
  allComponents_size= np.sum(labels>0)

  keepIfBiggerThr= min_size_to_preserve*allComponents_size
  final_img= np.ones_like(img)*minVal
  for l, count in labels_and_num:
    if count > keepIfBiggerThr:
      selection= (labels==l)
      final_img[selection]= img[selection]
  print("!")
  return final_img

def getLabels2NumVoxels(labels, n_labels):
  label2Size= np.zeros(n_labels+1, dtype=np.int64)
  flattenLabels= labels.flatten()
  for i in range(flattenLabels.shape[0]):
    label2Size[flattenLabels[i]]+=1

  return [ (l, size) for l, size in enumerate(label2Size) ][1:]

def estimateBackground(img):
  coordsNotInSphere= getCoordsOutsideSphere(img)
  background= np.mean( img[coordsNotInSphere])
  coordsInSmallerSphere= getCoordsWithinSphere(img,0.33)
  notBackground= np.mean( img[coordsInSmallerSphere])
  otsu_thr= threshold_otsu(img)
  thr_select= .5*(notBackground+otsu_thr)
  return background, thr_select

try:
  import numba
  getLabels2NumVoxels= numba.jit(getLabels2NumVoxels, nopython=True, cache=True, parallel=True)
except ImportError:
  pass

def test():
  from ..utils.ioUtils import loadVol
  fname= "devel_code/tests/data/EMD-7023_predLocscale.mrc"
  x= loadVol( fname )
  x_clean= removeSmallCCs(x)

if __name__=="__main__":
  test()

ext: py

max_stars_repo_name: rsanchezgarc/deepEMhancer

max_issues_repo_name: rsanchezgarc/deepEMhancer

max_forks_repo_name: rsanchezgarc/deepEMhancer

reason: import numpy,import numba


text:
# (C) Copyright 1996- ECMWF.
#
# This software is licensed under the terms of the Apache Licence Version 2.0
# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
# In applying this licence, ECMWF does not waive the privileges and immunities
# granted to it by virtue of its status as an intergovernmental organisation
# nor does it submit to any jurisdiction.

import unittest

import numpy as np
from context import data_file

import thermofeel.thermofeel as tfc

# import pytest

# combining the pytest library with the functionality of np testing for use with np.array file type

class TestThermalCalculator(unittest.TestCase):
    def setUp(self):
        # variables for use in thermalindexcalculator
        t = np.genfromtxt(
            data_file("thermofeel_testcases.csv"),
            delimiter=",",
            names=True,
        )
        self.t2m = t["t2m"]
        self.ssr = t["ssr"]
        self.td = t["td"]
        self.va = t["va"]
        self.mrt = t["mrt"]
        self.ssrd = t["ssrd"]
        self.strd = t["strd"]
        self.fdir = t["fdir"]
        self.strr = t["strr"]
        self.cossza = t["cossza"]

        # indices
        tr = np.genfromtxt(
            data_file("thermofeel_test_results.csv"),
            delimiter=",",
            names=True,
        )
        self.rh = tfc.calculate_saturation_vapour_pressure(self.t2m)
        self.rhpercent = tfc.calculate_relative_humidity_percent(self.t2m, self.td)
        self.heatindex = tr["heatindex"]
        self.heatindexadjusted = tr["heatindexadjust"]
        self.utci = tr["utci"]
        self.apparenttemperature = tr["apparenttemperature"]
        self.wbgts = tr["wbgts"]
        self.wbgt = tr["wbgt"]
        self.net = tr["net"]
        self.humidex = tr["humidex"]
        self.windchill = tr["windchill"]
        self.mrtr = tr["mrt"]
        self.mrtw = tr["mrtw"]

    def assert_equal(self, result, calculation):
        self.assertequal = self.assertIsNone(
            np.testing.assert_array_almost_equal(result, calculation, decimal=6)
        )

    def assert_equal_less_precise(self, result, calculation):
        self.assertequal = self.assertIsNone(
            np.testing.assert_array_almost_equal(result, calculation, decimal=1)
        )

    def test_relative_humidity(self):
        self.assert_equal(self.rh, tfc.calculate_saturation_vapour_pressure(self.t2m))

    def test_relative_humidity_percent(self):
        self.assert_equal(
            self.rhpercent, tfc.calculate_relative_humidity_percent(self.t2m, self.td)
        )

    def test_heat_index(self):
        self.assert_equal_less_precise(
            self.heatindex, tfc.calculate_heat_index_simplified(self.t2m)
        )

    def test_mean_radiant_temperature(self):
        self.assert_equal_less_precise(
            self.mrtr,
            tfc.calculate_mean_radiant_temperature(
                self.ssrd / 3600,
                self.ssr / 3600,
                self.fdir / 3600,
                self.strd / 3600,
                self.strr / 3600,
                self.cossza / 3600,
            ),
        )

    def test_utci(self):
        rh_pc = tfc.calculate_relative_humidity_percent(self.t2m, self.td)
        ehPa = tfc.calculate_saturation_vapour_pressure(self.t2m) * rh_pc / 100.0
        self.assert_equal_less_precise(
            self.utci,
            tfc.calculate_utci(
                t2_k=self.t2m, va_ms=self.va, mrt_k=self.mrt, e_hPa=ehPa
            ),
        )

    def test_apparent_temperature(self):
        self.assert_equal_less_precise(
            self.apparenttemperature,
            tfc.calculate_apparent_temperature(self.t2m, self.va),
        )

    def test_wbgts(self):
        self.assert_equal_less_precise(self.wbgts, tfc.calculate_wbgts(self.t2m))

    def test_net(self):
        self.assert_equal_less_precise(
            self.net,
            tfc.calculate_net_effective_temperature(self.t2m, self.va, self.td),
        )

    def test_humidex(self):
        self.assert_equal(self.humidex, tfc.calculate_humidex(self.t2m, self.td))

    def test_wind_chill(self):
        self.assert_equal(self.windchill, tfc.calculate_wind_chill(self.t2m, self.va))

    def test_heat_index_adjusted(self):
        self.assert_equal(
            self.heatindexadjusted, tfc.calculate_heat_index_adjusted(self.t2m, self.td)
        )

    def test_wbgt(self):
        self.assert_equal_less_precise(
            self.wbgt, tfc.calculate_wbgt(self.t2m, self.va, self.mrt)
        )

    def test_mrt_from_wbgt(self):
        wbgt_k = tfc.celcius_to_kelvin(self.wbgt)
        mrt = tfc.calculate_mrt_from_wbgt(self.t2m, wbgt_k, self.va)
        self.assert_equal_less_precise(self.mrtw, mrt)

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

max_stars_repo_name: milanavuckovic/thermofeel

max_issues_repo_name: milanavuckovic/thermofeel

max_forks_repo_name: milanavuckovic/thermofeel

reason: import numpy


text:
# Copyright (c) 2020 Mujammil Ali A S. All Rights Reserved.
#
# This module is Licensed under the Apache License, Version 2.0 with Commons
# Clause (the "License") as defined below, subject to the following condition.
#
#     https://commonsclause.com/
#
# “Commons Clause” License Condition v1.0
#
# Without limiting other conditions in the License, the grant of rights
# under the License will not include, and the License does not grant to you,
# the right to Sell the Software.

"""Streamlit based Evento Deployment Audit Module

This module acts as a web-based user interface and does the following stuff
    - Connects with Redis server and it is Evento protocol aware.
    - Capable of configuring ROI settings and checking the alert stream.
    - Capabel of handling five ROI regions.
    - Capable of handling three different event Type
    - Reads the ROI settings file with last saved settings and save with the
    latest change in ROI settings.

This module contains the following functions:
    - runStreamlitServer(): Runs the streamlit server and does either ROI
        settings configuration or alert stream based on selected application
        mode and event type.
    - connectRedis(): Connect with redis server instance once in the entire
        streamlit session by caching
    - configureROI(): Reads the ROI settings file and save it again with new
        information if altered. If the selected application mode is "Event
        Alert Stream", then this function is invoked.
    - streamAlert(): Once ROI settings configuration is in place, it streams
        the processed frame from Evento Video Analytics Pipeline Module based
        on the selected event type If the selected application mode is "ROI
        Configuration", then this function is invoked.
"""

import streamlit as st
import time
import numpy as np
import redis
import base64
import cv2
import json
import pandas as pd
import os
import argparse

contourROI = {}
ROIColor = {"ROI 1": (0,255,0), "ROI 2": (255,0,255), "ROI 3": (255,255,0), "ROI 4": (255,0,0), "ROI 5": (50,30,70)}

my_json = st.empty()

contourROI = {}

settingsFileName = {""}

@st.cache(allow_output_mutation=True)
def connectRedis():
    """This function connects with the Redis Server and cache the connection for the
    entire streamlit session.

    Args:
        Nothing
    Returns:
        prd (redis): Returns the redis server instance.
    """

    prd = redis.Redis(host='localhost', port=6379,charset="utf-8", decode_responses=True)

    if not prd.ping():
        raise Exception('Redis unavailable')

    return prd

def configureROI(prd, contourROI, eventType, FLAGS):
    """This function configures the ROI settings based on eventType. It parses the keys that obeys the Redis
    Evento protocol to get the latest or relevant frame/scene for which the ROI settings has to be done. It
    helps in creating new point to form contour for the selected ROI, delete the last point (x,y), or deletes
    the entire config using the slider and buttons. It finally stores the changes as json file for each event.

    Args:
        prd (redis): Returns the redis server instance.
        contourROI (dict): contains the last saved/currently modified json settings value as dict.
        eventType (int): Mapped integer value for type of the event.
        FLAGS: contains the command line arguments that holds the path of config file to be stored.

    Returns:
        Nothing. However, it stores the ROI settings as json file for each event.
    """

    image_placeholder = st.empty()

    #resp = prd.xrevrange(camera,count=1)
    #image_string = resp[0][1]['image']

    image_string = prd.get(eventType)

    img_data = base64.b64decode(image_string)
    nparr = np.fromstring(img_data, np.uint8)
    img_np = cv2.imdecode(nparr, cv2.IMREAD_COLOR)
    height, width, channel = img_np.shape
    overlay_np = img_np.copy()
    final_np = img_np.copy()

    whichROI = st.sidebar.selectbox('Select RoI',('ROI 1', 'ROI 2', 'ROI 3', 'ROI 4', 'ROI 5'))

    x = st.sidebar.slider('X Direction', 0, width-1)
    y = st.sidebar.slider('Y Direction', 0, height-1)

    alpha = st.sidebar.slider('Transperancy', 0.0, 1.0, 0.15, 0.005)

    if st.sidebar.button('<<<<<<Create Points>>>>>>'):
        if [x, y] not in list(contourROI[whichROI]):
            contourROI[whichROI].append([x, y])

    if st.sidebar.button('<<<<<Delete last Points>>>>'):
        if not list(contourROI[whichROI]) == []:
            del contourROI[whichROI][-1]

    if st.sidebar.button('<<<<<<<<Clear All>>>>>>>>'):
        contourROI = {'ROI 1':[],'ROI 2':[],'ROI 3':[], 'ROI 4':[], 'ROI 5':[]}

    for index, whichROI in enumerate(contourROI.keys()):
        if not list(contourROI[whichROI]) == []:

            ctr = np.array(contourROI[whichROI]).reshape((-1,1,2)).astype(np.int32)

            cv2.drawContours(overlay_np,[ctr],0,ROIColor[whichROI],-1)
            cv2.drawContours(overlay_np,[ctr],0,(255,255,255),4)
            cv2.addWeighted(overlay_np, alpha, final_np, 1 - alpha, 0, final_np)
            moment = cv2.moments(ctr)

            if (moment["m10"] != moment["m00"])  and  (moment["m01"] != moment["m00"]):

                centerX = int(moment["m10"] / moment["m00"])
                centerY = int(moment["m01"] / moment["m00"])
                cv2.putText(final_np, str(index + 1), (centerX, centerY), cv2.FONT_HERSHEY_SIMPLEX, 0.7, (0,0,0), 2)

    with open(os.path.join(FLAGS.app_config_dir, eventType + '.json'), 'w') as json_file:
        json.dump(contourROI, json_file)

    img_np = cv2.circle(final_np, (x, y),10,(255,0,0),-1)

    image_placeholder.image(final_np, channels="BGR")

    #df = pd.DataFrame.from_dict(contourROI)
    df = pd.DataFrame(dict([(k,pd.Series(v)) for k,v in contourROI.items() ]))

    st.subheader('ROI Points')

    st.write(df)

def streamAlert(app_mode, prd, contourROI):
    """This function stream the processed video as images to the image container. It gets the video feed from Redis
    Stream of Evento video Analytics Pipeline.

    Args:
        app_mode (int): Integer value mapped to application mode. 1 - ROI Configuration and 2 - Event Alert Stream
        prd (redis): Returns the redis server instance.
        contourROI (dict): contains the last saved/currently modified json settings value as dict.

    Returns:
        Nothing. It streams forever until the application mode is changed to "ROI Configuration". Yes, it abuses the
         while loop usage which breaks the Steamlit execution model at the time of developint this.
    """

    image_placeholder_ias = st.empty()
    st.write(" ")
    st.write(" ")

    videoResolution = st.sidebar.selectbox('Select Resolution',('Original', '360*240'))

    st.sidebar.title(" ")
    st.sidebar.title(" ")
    st.sidebar.title(" ")
    st.sidebar.title(" ")
    st.sidebar.title(" ")
    st.sidebar.title(" ")

    while True:

        if app_mode != "Event Alert Stream":
            break

        resp = prd.xrevrange("camera:0",count=1)
        image_string = resp[0][1]['image']

        appMetadata = json.loads(resp[0][1]['appMetadata'])

        img_data = base64.b64decode(image_string)
        nparr = np.fromstring(img_data, np.uint8)
        img_np = cv2.imdecode(nparr, cv2.IMREAD_COLOR)

        if videoResolution == "360*240":
            img_np = cv2.resize(img_np,(360,240))

        image_placeholder_ias.image(img_np, channels="BGR")

        time.sleep(0.025)

def runStreamlitServer(FLAGS):
    """This function starts the streamlit server with all UI components and it allows user to change the application
    mode and event Type. It behaves differently based on the changes in corresponding UI components. It reads
    the ROI settings Json file and store it in dict. It also handles the Evento Protocol control signal.

    Args:
        FLAGS: contains the command line arguments that holds the path of config file to be stored.

    Returns:
        Nothing. It streams forever until the application mode is changed to "ROI Configuration". Yes, it abuses the
         while loop usage which breaks the Steamlit execution model at the time of developint this.
    """
    dummy = 0

    st.markdown("
Event Alert Manager
", unsafe_allow_html=True)

    st.subheader('Video Stream')

    st.sidebar.title('Control Panel')

    app_mode = st.sidebar.selectbox("Application Mode",
        ("ROI Configuration", "Event Alert Stream"))

    eventType = st.sidebar.selectbox("Event",
    ("EventStat", "MonitorVehicle", "AlertStranger"))

    configFile = os.path.join(FLAGS.app_config_dir,eventType + '.json')

    if not os.path.exists(configFile):
        contourROI = {'ROI 1':[],'ROI 2':[],'ROI 3':[], 'ROI 4':[], 'ROI 5':[]}
        with open(configFile, 'w') as json_file:
            json.dump(contourROI, json_file)
    else:
        contourROI = json.load(open(configFile))

    prd = connectRedis()

    if app_mode == "ROI Configuration":

        prd.set("app_mode", 1)

        configureROI(prd, contourROI, eventType, FLAGS)

    elif app_mode == "Event Alert Stream":

        prd.set("app_mode", 2)

        if eventType == "EventStat":
            prd.set("event_type", 0)

        elif eventType == "AlertStranger":
            prd.set("event_type", 1)

        elif eventType == "MonitorVehicle":
            prd.set("event_type", 2)

        streamAlert(app_mode, prd, contourROI)

    else:

        prd.set("app_mode", 0)

if __name__ == '__main__':
#   """
#   This Function is the entry point for this module. It collects the
#   user arguments based on argparse library. The list of options are
#
#       --app_config_dir - Directory where the json files that has the ROI
#            settings information. Every event has each file of json.
#
#   How to run in CLI: streamlit run evento/ui_server.py -- --app_config_dir "config/"
#
#   Finally it triggers the streamlit server with web interface.
#   """
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument(
        "--app_config_dir",
        type=str,
        default='config',
        help=("Directory include:'AlertStranger.json', 'EventStat.json', "
              "'MonitorVehicle.json', created by ui_server.py."),
        required=False)

    FLAGS = parser.parse_args()

    runStreamlitServer(FLAGS)

max_stars_repo_name: makemypoc/evento

max_issues_repo_name: makemypoc/evento

max_forks_repo_name: makemypoc/evento


text:
"""Module with a class for Jenoptik ceilometer."""
from collections import namedtuple
import netCDF4
import numpy as np
import numpy.ma as ma
from cloudnetpy.instruments.ceilometer import Ceilometer
from cloudnetpy import utils

instrument_info = namedtuple('instrument_info',
                             ['calibration_factor',
                              'overlap_function_params',
                              'is_range_corrected'])

CEILOMETER_INFO = {
    'punta-arenas': instrument_info(
        calibration_factor=1e-12,
        overlap_function_params=(0, 1),
        is_range_corrected=True),
    'mace-head': instrument_info(
        calibration_factor=4.5e-11,
        overlap_function_params=(500, 200),
        is_range_corrected=False)
}

class JenoptikCeilo(Ceilometer):
    """Class for Jenoptik chm15k ceilometer."""
    def __init__(self, file_name, site_name):
        super().__init__(file_name)
        self.model = 'chm15k'
        self.dataset = netCDF4.Dataset(self.file_name)
        self.variables = self.dataset.variables
        self.noise_params = (70, 2e-14, 0.3e-6, (1e-9, 4e-9))
        self.calibration_info = _read_calibration_info(site_name)

    def read_ceilometer_file(self):
        """Reads data and metadata from Jenoptik netCDF file."""
        self.range = self._calc_range()
        self.time = self._convert_time()
        self.date = self._read_date()
        self.backscatter = self._convert_backscatter()
        self.metadata = self._read_metadata()

    def _calc_range(self):
        """Assumes 'range' means the upper limit of range gate."""
        ceilo_range = self._getvar('range')
        return ceilo_range - utils.mdiff(ceilo_range)/2

    def _convert_time(self):
        time = self.variables['time']
        try:
            assert all(np.diff(time) > 0)
        except AssertionError:
            raise RuntimeError('Inconsistent ceilometer time stamps.')
        if max(time) > 24:
            time = utils.seconds2hours(time)
        return time

    def _read_date(self):
        return [self.dataset.year, self.dataset.month, self.dataset.day]

    def _convert_backscatter(self):
        """Steps to convert Jenoptik SNR to raw beta."""
        beta_raw = self._getvar('beta_raw')
        data_std = self._getvar('stddev')
        normalised_apd = self._get_nn()
        beta_raw *= utils.transpose(data_std / normalised_apd)
        if not self.calibration_info.is_range_corrected:
            beta_raw *= self.range ** 2
        overlap_function = _get_overlap(self.range, self.calibration_info)
        beta_raw /= overlap_function
        beta_raw *= self.calibration_info.calibration_factor
        return beta_raw

    def _get_nn(self):
        """TODO: Should be checked."""
        nn1 = self._getvar('nn1', 'NN1')
        median_nn1 = ma.median(nn1)
        if 120 < median_nn1 < 160:
            step_factor, reference, scale = 1.24, 140, 5
        elif 3200 < median_nn1 < 4000:
            step_factor, reference, scale = 1.035, 3685, 1
        else:
            return 1
        return step_factor ** (-(nn1-reference) / scale)

    def _getvar(self, *args):
        """Reads data of variable (array or scalar) from netcdf-file."""
        for arg in args:
            if arg in self.variables:
                var = self.variables[arg]
                return var[0] if utils.isscalar(var) else var[:]

    def _read_metadata(self):
        meta = {'tilt_angle': self._getvar('zenith')}
        return meta

def _get_overlap(range_ceilo, calibration_info):
    """Approximative overlap function."""
    params = calibration_info.overlap_function_params
    return utils.array_to_probability(range_ceilo, *params)

def _read_calibration_info(site_name):
    if 'punta' in site_name.lower():
        return CEILOMETER_INFO['punta-arenas']
    elif 'mace' in site_name.lower():
        return CEILOMETER_INFO['mace-head']

max_stars_repo_name: tukiains/actris-cloudnet

max_issues_repo_name: tukiains/actris-cloudnet

max_forks_repo_name: tukiains/actris-cloudnet


text:
import timeit
import math
import numpy as np
import gc
import torch
import cProfile

log_output = ""

def log(s):
    print(s)
    global log_output
    log_output = log_output + s + "\n"

# short hands

# torch quat/vector utils

def to_torch(x, dtype=torch.float, device='cuda:0', requires_grad=False):
    return torch.tensor(x, dtype=dtype, device=device, requires_grad=requires_grad)

@torch.jit.script
def quat_mul(a, b):
    assert a.shape == b.shape
    shape = a.shape
    a = a.reshape(-1, 4)
    b = b.reshape(-1, 4)

    x1, y1, z1, w1 = a[:, 0], a[:, 1], a[:, 2], a[:, 3]
    x2, y2, z2, w2 = b[:, 0], b[:, 1], b[:, 2], b[:, 3]
    ww = (z1 + x1) * (x2 + y2)
    yy = (w1 - y1) * (w2 + z2)
    zz = (w1 + y1) * (w2 - z2)
    xx = ww + yy + zz
    qq = 0.5 * (xx + (z1 - x1) * (x2 - y2))
    w = qq - ww + (z1 - y1) * (y2 - z2)
    x = qq - xx + (x1 + w1) * (x2 + w2)
    y = qq - yy + (w1 - x1) * (y2 + z2)
    z = qq - zz + (z1 + y1) * (w2 - x2)

    quat = torch.stack([x, y, z, w], dim=-1).view(shape)

    return quat

@torch.jit.script
def normalize(x, eps: float = 1e-9):
    return x / x.norm(p=2, dim=-1).clamp(min=eps, max=None).unsqueeze(-1)

@torch.jit.script
def quat_apply(a, b):
    shape = b.shape
    a = a.reshape(-1, 4)
    b = b.reshape(-1, 3)
    xyz = a[:, :3]
    t = xyz.cross(b, dim=-1) * 2
    return (b + a[:, 3:] * t + xyz.cross(t, dim=-1)).view(shape)

@torch.jit.script
def quat_rotate(q, v):
    shape = q.shape
    q_w = q[:, -1]
    q_vec = q[:, :3]
    a = v * (2.0 * q_w ** 2 - 1.0).unsqueeze(-1)
    b = torch.cross(q_vec, v, dim=-1) * q_w.unsqueeze(-1) * 2.0
    c = q_vec * \
        torch.bmm(q_vec.view(shape[0], 1, 3), v.view(
            shape[0], 3, 1)).squeeze(-1) * 2.0
    return a + b + c

@torch.jit.script
def quat_rotate_inverse(q, v):
    shape = q.shape
    q_w = q[:, -1]
    q_vec = q[:, :3]
    a = v * (2.0 * q_w ** 2 - 1.0).unsqueeze(-1)
    b = torch.cross(q_vec, v, dim=-1) * q_w.unsqueeze(-1) * 2.0
    c = q_vec * \
        torch.bmm(q_vec.view(shape[0], 1, 3), v.view(
            shape[0], 3, 1)).squeeze(-1) * 2.0
    return a - b + c

@torch.jit.script
def quat_axis(q, axis=0):
    # type: (Tensor, int) -> Tensor
    basis_vec = torch.zeros(q.shape[0], 3, device=q.device)
    basis_vec[:, axis] = 1
    return quat_rotate(q, basis_vec)

@torch.jit.script
def quat_conjugate(a):
    shape = a.shape
    a = a.reshape(-1, 4)
    return torch.cat((-a[:, :3], a[:, -1:]), dim=-1).view(shape)

@torch.jit.script
def quat_unit(a):
    return normalize(a)

@torch.jit.script
def quat_from_angle_axis(angle, axis):
    theta = (angle / 2).unsqueeze(-1)
    xyz = normalize(axis) * theta.sin()
    w = theta.cos()
    return quat_unit(torch.cat([xyz, w], dim=-1))

@torch.jit.script
def normalize_angle(x):
    return torch.atan2(torch.sin(x), torch.cos(x))

@torch.jit.script
def tf_inverse(q, t):
    q_inv = quat_conjugate(q)
    return q_inv, -quat_apply(q_inv, t)

@torch.jit.script
def tf_apply(q, t, v):
    return quat_apply(q, v) + t

@torch.jit.script
def tf_vector(q, v):
    return quat_apply(q, v)

@torch.jit.script
def tf_combine(q1, t1, q2, t2):
    return quat_mul(q1, q2), quat_apply(q1, t2) + t1

@torch.jit.script
def get_basis_vector(q, v):
    return quat_rotate(q, v)

def mem_report():
    '''Report the memory usage of the tensor.storage in pytorch
    Both on CPUs and GPUs are reported'''

    def _mem_report(tensors, mem_type):
        '''Print the selected tensors of type
        There are two major storage types in our major concern:
            - GPU: tensors transferred to CUDA devices
            - CPU: tensors remaining on the system memory (usually unimportant)
        Args:
            - tensors: the tensors of specified type
            - mem_type: 'CPU' or 'GPU' in current implementation '''
        total_numel = 0
        total_mem = 0
        visited_data = []
        for tensor in tensors:
            if tensor.is_sparse:
                continue
            # a data_ptr indicates a memory block allocated
            data_ptr = tensor.storage().data_ptr()
            if data_ptr in visited_data:
                continue
            visited_data.append(data_ptr)

            numel = tensor.storage().size()
            total_numel += numel
            element_size = tensor.storage().element_size()
            mem = numel*element_size /1024/1024 # 32bit=4Byte, MByte
            total_mem += mem
            element_type = type(tensor).__name__
            size = tuple(tensor.size())

            # print('%s\t\t%s\t\t%.2f' % (
            #     element_type,
            #     size,
            #     mem) )
        print('Type: %s Total Tensors: %d \tUsed Memory Space: %.2f MBytes' % (mem_type, total_numel, total_mem) )

    gc.collect()

    LEN = 65
    objects = gc.get_objects()
    #print('%s\t%s\t\t\t%s' %('Element type', 'Size', 'Used MEM(MBytes)') )
    tensors = [obj for obj in objects if torch.is_tensor(obj)]
    cuda_tensors = [t for t in tensors if t.is_cuda]
    host_tensors = [t for t in tensors if not t.is_cuda]
    _mem_report(cuda_tensors, 'GPU')
    _mem_report(host_tensors, 'CPU')
    print('='*LEN)

def grad_norm(params):
    grad_norm = 0.
    for p in params:
        if p.grad is not None:
            grad_norm += torch.sum(p.grad ** 2)
    return torch.sqrt(grad_norm)

def print_leaf_nodes(grad_fn, id_set):
    if grad_fn is None:
        return
    if hasattr(grad_fn, 'variable'):
        mem_id = id(grad_fn.variable)
        if not(mem_id in id_set):
            print('is leaf:', grad_fn.variable.is_leaf)
            print(grad_fn.variable)
            id_set.add(mem_id)

    # print(grad_fn)
    for i in range(len(grad_fn.next_functions)):
        print_leaf_nodes(grad_fn.next_functions[i][0], id_set)

def policy_kl(p0_mu, p0_sigma, p1_mu, p1_sigma):
    c1 = torch.log(p1_sigma/p0_sigma + 1e-5)
    c2 = (p0_sigma**2 + (p1_mu - p0_mu)**2)/(2.0 * (p1_sigma**2 + 1e-5))
    c3 = -1.0 / 2.0
    kl = c1 + c2 + c3
    kl = kl.sum(dim=-1) # returning mean between all steps of sum between all actions
    return kl.mean()

max_stars_repo_name: eanswer/mbrl-lib

max_issues_repo_name: eanswer/mbrl-lib

max_forks_repo_name: eanswer/mbrl-lib


text:
import mmcv
import numpy as np
import pycocotools.mask as mask_util
import torch
import torch.nn as nn
import torch.nn.functional as F
from torch.nn.modules.utils import _pair

from mmdet.core import auto_fp16, force_fp32, mask_target
from mmdet.ops.carafe import CARAFEPack
from ..builder import build_loss
from ..registry import HEADS
from ..utils import ConvModule, build_upsample_layer
import pickle
import time
import pdb

@HEADS.register_module
class MGANHead(nn.Module):

    def __init__(self,
                 num_convs=2,
                 roi_feat_size=7,
                 in_channels=512,
                 conv_out_channels=512,
                 conv_cfg=None,
                 norm_cfg=None,
                 class_agnostic=False,
                 loss_mask=dict(
                     type='CrossEntropyLoss', # todo: binary?
                     use_mask=True,
                     loss_weight=0.3),
                 loss_occ=dict(
                     type='CrossEntropyLoss',
                     use_mask=True,
                     loss_weight=1.0
                 )):
        super(MGANHead, self).__init__()
        self.num_convs = num_convs
        self.roi_feat_size = roi_feat_size
        self.in_channels = in_channels
        self.conv_out_channels = conv_out_channels

        self.conv_cfg = conv_cfg
        self.norm_cfg = norm_cfg
        self.fp16_enabled = False

        self.convs = nn.ModuleList()
        for i in range(self.num_convs):
            in_channels = (
                self.in_channels if i == 0 else self.conv_out_channels)
            self.convs.append(
                ConvModule(
                    in_channels,
                    self.conv_out_channels,
                    3,
                    padding=1,
                    conv_cfg=conv_cfg,
                    norm_cfg=norm_cfg))
        logits_in_channel = self.conv_out_channels
        self.conv_logits = nn.Conv2d(logits_in_channel, 1, 1)
        self.relu = nn.ReLU(inplace=True)
        self.debug_imgs = None
        self.class_agnostic = class_agnostic
        self.loss_mask = build_loss(loss_mask)
        self.loss_occ = build_loss(loss_occ)

    @auto_fp16()
    def forward(self, x):
        for conv in self.convs:
            x = conv(x)
        # x = self.conv_logits(x).sigmoid() * x
        x_fpm = self.conv_logits(x).sigmoid()
        x = x_fpm * x

        # with open('/mmdetection/result/attent_vis/attention_fpm/' + str(time.time()), 'wb') as pkl_file:
        #     pickle.dump(x_fpm, pkl_file)
        return x

    def get_target(self, sampling_results, gt_masks, rcnn_train_cfg):
        pos_proposals = [res.pos_bboxes for res in sampling_results]
        pos_assigned_gt_inds = [
            res.pos_assigned_gt_inds for res in sampling_results
        ]
        mask_targets = mask_target(pos_proposals, pos_assigned_gt_inds,
                                   gt_masks, rcnn_train_cfg)
        return mask_targets

    def get_bbox_target(self, sampling_results, gt_bboxes, gt_labels,
            rcnn_train_cfg):
        pos_proposals = [res.pos_bboxes for res in sampling_results]
        # neg_proposals = [res.neg_bboxes for res in sampling_results]
        pos_gt_bboxes = [res.pos_gt_bboxes for res in sampling_results]
        pos_gt_labels = [res.pos_gt_labels for res in sampling_results]
        reg_classes = 1
        cls_reg_targets = bbox_target(
                pos_proposals,
                # neg_proposals,
                pos_gt_bboxes,
                pos_gt_labels,
                rcnn_train_cfg,
                reg_classes,
                target_means=self.target_means,
                target_stds=self.target_stds)
        return cls_reg_targets

    @force_fp32(apply_to=('mask_pred', ))
    def loss(self, mask_pred, mask_targets, labels):
        loss = dict()
        if self.class_agnostic:
            loss_mask = self.loss_mask(mask_pred, mask_targets,
                                       torch.zeros_like(labels))
        else:
            loss_mask = self.loss_mask(mask_pred, mask_targets, labels)
        loss['loss_mask'] = loss_mask
        return loss

    def loss_occ(self,
             cls_score,
             labels,
             occ_weights,
             ):
        loss = dict()
        loss_sum = []
        batch_size = cls_score.shape[0]
        for i in range(batch_size):
            loss_i = F.cross_entropy(cls_score[i:i+1], labels[i:i+1])
            occ = occ_weights[i]
            occ_level = 1 - sum(sum(occ)) / (occ.shape[0] * occ.shape[1])
            loss_sum.append(occ_level * loss_i)

        loss["occ"] = sum(loss_sum) / batch_size

        return loss

    def get_seg_masks(self, mask_pred, det_bboxes, det_labels, rcnn_test_cfg,
                      ori_shape, scale_factor, rescale):
        """Get segmentation masks from mask_pred and bboxes.

        Args:
            mask_pred (Tensor or ndarray): shape (n, #class+1, h, w).
                For single-scale testing, mask_pred is the direct output of
                model, whose type is Tensor, while for multi-scale testing,
                it will be converted to numpy array outside of this method.
            det_bboxes (Tensor): shape (n, 4/5)
            det_labels (Tensor): shape (n, )
            img_shape (Tensor): shape (3, )
            rcnn_test_cfg (dict): rcnn testing config
            ori_shape: original image size

        Returns:
            list[list]: encoded masks
        """
        if isinstance(mask_pred, torch.Tensor):
            mask_pred = mask_pred.sigmoid().cpu().numpy()
        assert isinstance(mask_pred, np.ndarray)
        # when enabling mixed precision training, mask_pred may be float16
        # numpy array
        mask_pred = mask_pred.astype(np.float32)

        cls_segms = [[] for _ in range(self.num_classes - 1)]
        bboxes = det_bboxes.cpu().numpy()[:, :4]
        labels = det_labels.cpu().numpy() + 1

        if rescale:
            img_h, img_w = ori_shape[:2]
        else:
            img_h = np.round(ori_shape[0] * scale_factor).astype(np.int32)
            img_w = np.round(ori_shape[1] * scale_factor).astype(np.int32)
            scale_factor = 1.0

        for i in range(bboxes.shape[0]):
            if not isinstance(scale_factor, (float, np.ndarray)):
                scale_factor = scale_factor.cpu().numpy()
            bbox = (bboxes[i, :] / scale_factor).astype(np.int32)
            label = labels[i]
            w = max(bbox[2] - bbox[0] + 1, 1)
            h = max(bbox[3] - bbox[1] + 1, 1)

            if not self.class_agnostic:
                mask_pred_ = mask_pred[i, label, :, :]
            else:
                mask_pred_ = mask_pred[i, 0, :, :]

            bbox_mask = mmcv.imresize(mask_pred_, (w, h))
            bbox_mask = (bbox_mask > rcnn_test_cfg.mask_thr_binary).astype(
                np.uint8)

            if rcnn_test_cfg.get('crop_mask', False):
                im_mask = bbox_mask
            else:
                im_mask = np.zeros((img_h, img_w), dtype=np.uint8)
                im_mask[bbox[1]:bbox[1] + h, bbox[0]:bbox[0] + w] = bbox_mask

            if rcnn_test_cfg.get('rle_mask_encode', True):
                rle = mask_util.encode(
                    np.array(im_mask[:, :, np.newaxis], order='F'))[0]
                cls_segms[label - 1].append(rle)
            else:
                cls_segms[label - 1].append(im_mask)

        return cls_segms


text:
from PIL import Image
from PIL.ImageFilter import GaussianBlur
import numpy as np
from scipy.signal import convolve2d

def spatial_kernel(s=1, channels=1, reduce=True):
    l = int(s * 2) | 1  # determine kernel size
    x = np.tile((np.arange(l) - l // 2).reshape(1, l), (l, 1))
    y = np.transpose(x)
    d = x**2 + y**2
    g = np.exp(-d / (2 * (s**2))) / (2 * np.pi * (s**2))
    if reduce:
        g /= np.sum(g.ravel())
    g = np.dstack([g] * channels)
    return g

def range_kernel(I, s=1, reduce=True):
    assert I.shape[0] == I.shape[1]
    l, c = I.shape[0], I.shape[-1]
    I = I.copy().astype(np.float64)
    I = np.abs((I - I[l // 2, l // 2]))  # minus center pixel value
    I = np.max(I, axis=-1)
    k = np.exp(-I / (2 * (s**2))) / (2 * np.pi * (s**2))
    if reduce:
        k /= np.sum(k.ravel())
    k = np.dstack([k] * c)
    return k

def bilateral_filter(img, s1, s2):
    img = np.array(img.convert('RGB'))

    channels = img.shape[-1]
    g = spatial_kernel(s = s1, channels = channels, reduce = False)
    L = g.shape[0]

    h, w = img.shape[0], img.shape[1]
    pad_width = ((L//2, L//2), (L//2, L//2), (0, 0))
    padded_img = np.pad(img, pad_width = pad_width, mode = "edge").astype(np.float64)
    filtered = np.zeros(img.shape, dtype=np.float64)
    for i in range(L // 2, h + L // 2):
        for j in range(L // 2, w + L // 2):
            I = padded_img[i - L // 2: i + L // 2 + 1, j - L // 2: j + L // 2 + 1, :]  # extract patch
            k = range_kernel(I, s = s2, reduce = False)
            W = np.multiply(g, k).reshape(-1, channels)
            I = I.reshape(-1, channels)
            p = np.sum(np.multiply(I, W), axis=0) / np.sum(W, axis=0)
            filtered[i - L // 2, j - L // 2, :] = p
    filtered = filtered.astype(np.uint8)
    return Image.fromarray(filtered)

max_stars_repo_name: mashaashkolnik/image-editor

reason: import numpy,from scipy


text:
import json, os, argparse, sys, time, torch
from unittest import result
import torch2trt
import trt_pose.coco
import trt_pose.models
from torch2trt import TRTModule
import trt_pose.models
import cv2
import torchvision.transforms as transforms
import PIL.Image
from trt_pose.draw_objects import DrawObjects
from trt_pose.parse_objects import ParseObjects

import common
import tensorrt as trt
import numpy as np

def check_path(path):
    if not os.path.exists(path):
        raise Exception("{} is not exists".format(path))

# -----------------------------------------------------------------------------
parser = argparse.ArgumentParser()
parser.add_argument("-e", "--engine", default="resnet18_baseline_att_224x224_A_epoch_249.engine", help="the path of tensorrt engine", type=str)
parser.add_argument("-j", "--json", default="./human_pose.json", type=str)
parser.add_argument("-s", "--size", default="3,224,224",help="the input size of the torch model.", type=str)
parser.add_argument("-b", "--batch", default=1, help="the value of batch size.", type=int)
args = parser.parse_args()

# -----------------------------------------------------------------------------
# define variable

jfile_path = args.json
engine_path = args.engine

# check the path is exists
[ check_path(path) for path in [ jfile_path, engine_path ] ]

# input_size = [ int(i) for i in args.size.split(",")] # chaneel, height, width
if "224" in engine_path:
    input_size = [3, 224, 224]
elif "256" in engine_path:
    input_size = [3, 256, 256]

batch_size = args.batch

# -----------------------------------------------------------------------------
print('Get human pose parser ... ', end="", flush=True)
t0 = time.time()
with open(jfile_path, 'r') as f:
    human_pose = json.load(f)
topology = trt_pose.coco.coco_category_to_topology(human_pose)
num_parts, num_links = len(human_pose['keypoints']), len(human_pose['skeleton'])
print("{}s".format(round(time.time()-t0,3)))

# -----------------------------------------------------------------------------
print('Load model ... ', end="", flush=True)
t0 = time.time()

trt_logger = trt.Logger(trt.Logger.WARNING)     # TensorRT logger singleton
trt.init_libnvinfer_plugins(trt_logger, '')     # We first load all custom plugins shipped with TensorRT
runtime = trt.Runtime(trt_logger)               # Initialize runtime needed for loading TensorRT engine from file

with open(engine_path, 'rb') as serial_engine:                      # read engine from serial data
    engine = runtime.deserialize_cuda_engine(serial_engine.read())

inputs, outputs, bindings, stream = common.allocate_buffers(engine, 1)  # allocate buffer
context = engine.create_execution_context()                             # create excute object

print("{}s".format(round(time.time()-t0,3)))

# -----------------------------------------------------------------------------
print('Do inference with testing data ... ', end="", flush=True)
t0 = time.time()
data = torch.zeros((batch_size, input_size[0], input_size[1], input_size[2]))
np.copyto(inputs[0].host, data.ravel())
results = common.do_inference(context, bindings=bindings, inputs=inputs, outputs=outputs, stream=stream)
print("{}s (FPS:{})".format(round(time.time() - t0,3), int(50.0/(time.time() - t0))))
# -----------------------------------------------------------------------------


text:
"""Tests."""

from math import ceil  # type: ignore
import datetime  # type: ignore
import pytest  # type: ignore
import pandas as pd  # type: ignore
import numpy as np  # type: ignore
import altair as alt  # type: ignore

from src.penn_chime.charts import new_admissions_chart, admitted_patients_chart, chart_descriptions
from src.penn_chime.models import SimSirModel, sir, sim_sir_df, build_admits_df
from src.penn_chime.parameters import Parameters
from src.penn_chime.presentation import display_header
from src.penn_chime.settings import DEFAULTS
from src.penn_chime.defaults import RateLos

PARAM = Parameters(
    current_hospitalized=100,
    doubling_time=6.0,
    known_infected=5000,
    market_share=0.05,
    relative_contact_rate=0.15,
    susceptible=500000,
    hospitalized=RateLos(0.05, 7),
    icu=RateLos(0.02, 9),
    ventilated=RateLos(0.01, 10),
    n_days=60,
)

MODEL = SimSirModel(PARAM)

# set up

# we just want to verify that st _attempted_ to render the right stuff
# so we store the input, and make sure that it matches what we expect
class MockStreamlit:
    def __init__(self):
        self.render_store = []
        self.markdown = self.just_store_instead_of_rendering
        self.latex = self.just_store_instead_of_rendering
        self.subheader = self.just_store_instead_of_rendering

    def just_store_instead_of_rendering(self, inp, *args, **kwargs):
        self.render_store.append(inp)
        return None

    def cleanup(self):
        """
        Call this after every test, unless you intentionally want to accumulate stuff-to-render
        """
        self.render_store = []

st = MockStreamlit()

# test presentation

def test_penn_logo_in_header():
    penn_css = '<link rel="stylesheet" href="https://www1.pennmedicine.org/styles/shared/penn-medicine-header.css">'
    display_header(st, MODEL, PARAM)
    assert len(
        list(filter(lambda s: penn_css in s, st.render_store))
    ), "The Penn Medicine header should be printed"

def test_the_rest_of_header_shows_up():
    random_part_of_header = "implying an effective $R_t$ of"
    assert len(
        list(filter(lambda s: random_part_of_header in s, st.render_store))
    ), "The whole header should render"

def test_mitigation_statement():
    st.cleanup()
    expected_doubling = "outbreak **reduces the doubling time to 7.8** days"
    display_header(st, MODEL, PARAM)
    assert [s for s in st.render_store if expected_doubling in s]
    # assert len((list(filter(lambda s: expected_doubling in s, st.render_store))))
    st.cleanup()
    expected_halving = "outbreak **halves the infections every 51.9** days"
    halving_params = Parameters(
        current_hospitalized=100,
        doubling_time=6.0,
        known_infected=5000,
        market_share=0.05,
        relative_contact_rate=0.7,
        susceptible=500000,
        hospitalized=RateLos(0.05, 7),
        icu=RateLos(0.02, 9),
        ventilated=RateLos(0.01, 10),
        n_days=60,
    )
    halving_model = SimSirModel(halving_params)
    display_header(st, halving_model, halving_params)
    assert [s for s in st.render_store if expected_halving in s]
    #assert len((list(filter(lambda s: expected_halving in s, st.render_store))))
    st.cleanup()

st.cleanup()

@pytest.mark.xfail()
def test_header_fail():
    """
    Just proving to myself that these tests work
    """
    some_garbage = "ajskhlaeHFPIQONOI8QH34TRNAOP8ESYAW4"
    display_header(st, PARAM)
    assert len(
        list(filter(lambda s: some_garbage in s, st.render_store))
    ), "This should fail"
    st.cleanup()

def test_defaults_repr():
    """
    Test DEFAULTS.repr
    """
    repr(DEFAULTS)

# Test the math

def test_sir():
    """
    Someone who is good at testing, help
    """
    sir_test = sir(100, 1, 0, 0.2, 0.5, 1)
    assert sir_test == (
        0.7920792079207921,
        0.20297029702970298,
        0.0049504950495049506,
    ), "This contrived example should work"

    assert isinstance(sir_test, tuple)
    for v in sir_test:
        assert isinstance(v, float)
        assert v >= 0

    # Certain things should *not* work
    with pytest.raises(TypeError) as error:
        sir("S", 1, 0, 0.2, 0.5, 1)
    assert str(error.value) == "can't multiply sequence by non-int of type 'float'"

    with pytest.raises(TypeError) as error:
        sir(100, "I", 0, 0.2, 0.5, 1)
    assert str(error.value) == "can't multiply sequence by non-int of type 'float'"

    with pytest.raises(TypeError) as error:
        sir(100, 1, "R", 0.2, 0.5, 1)
    assert str(error.value) == "unsupported operand type(s) for +: 'float' and 'str'"

    with pytest.raises(TypeError) as error:
        sir(100, 1, 0, "beta", 0.5, 1)
    assert str(error.value) == "bad operand type for unary -: 'str'"

    with pytest.raises(TypeError) as error:
        sir(100, 1, 0, 0.2, "gamma", 1)
    assert str(error.value) == "unsupported operand type(s) for -: 'float' and 'str'"

    with pytest.raises(TypeError) as error:
        sir(100, 1, 0, 0.2, 0.5, "N")
    assert str(error.value) == "unsupported operand type(s) for /: 'str' and 'float'"

    # Zeros across the board should fail
    with pytest.raises(ZeroDivisionError):
        sir(0, 0, 0, 0, 0, 0)

def test_sim_sir():
    """
    Rounding to move fast past decimal place issues
    """
    raw_df = sim_sir_df(5, 6, 7, 0.1, 0.1, 40)

    first = raw_df.iloc[0, :]
    last = raw_df.iloc[-1, :]

    assert round(first.susceptible, 0) == 5
    assert round(first.infected, 2) == 6
    assert round(first.recovered, 0) == 7
    assert round(last.susceptible, 2) == 0
    assert round(last.infected, 2) == 0.18
    assert round(last.recovered, 2) == 17.82

    assert isinstance(raw_df, pd.DataFrame)

def test_new_admissions_chart():
    projection_admits = pd.read_csv("tests/projection_admits.csv")
    chart = new_admissions_chart(alt, projection_admits, PARAM)
    assert isinstance(chart, alt.Chart)
    assert chart.data.iloc[1].hospitalized < 1
    assert round(chart.data.iloc[40].icu, 0) == 25

    # test fx call with no params
    with pytest.raises(TypeError):
        new_admissions_chart()

    empty_chart = new_admissions_chart(alt, pd.DataFrame(), PARAM)
    assert empty_chart.data.empty

def test_admitted_patients_chart():
    census_df = pd.read_csv("tests/census_df.csv")
    chart = admitted_patients_chart(alt, census_df, PARAM)
    assert isinstance(chart, alt.Chart)
    assert chart.data.iloc[1].hospitalized == 1
    assert chart.data.iloc[49].ventilated == 203

    # test fx call with no params
    with pytest.raises(TypeError):
        admitted_patients_chart()

    empty_chart = admitted_patients_chart(alt, pd.DataFrame(), PARAM)
    assert empty_chart.data.empty

def test_model(model=MODEL, param=PARAM):
    # test the Model

    assert model.infected == 40000.0
    assert isinstance(model.infected, float)  # based off note in models.py

    # test the class-calculated attributes
    assert model.detection_probability == 0.125
    assert model.intrinsic_growth_rate == 0.12246204830937302
    assert model.beta == 3.2961405355450555e-07
    assert model.r_t == 2.307298374881539
    assert model.r_naught == 2.7144686763312222
    assert model.doubling_time_t == 7.764405988534983

    # test the things n_days creates, which in turn tests sim_sir, sir, and get_dispositions
    assert len(model.raw_df) == param.n_days + 1 == 61

    raw_df = model.raw_df
    first = raw_df.iloc[0, :]
    second = raw_df.iloc[1, :]
    last = raw_df.iloc[-1, :]

    assert first.susceptible == 500000.0
    assert round(second.infected, 0) == 43735

    assert round(last.susceptible, 0) == 67202
    assert round(raw_df.recovered[30], 0) == 224048

    assert [d[0] for d in model.dispositions.values()] == [100.0, 40.0, 20.0]
    assert [round(d[60], 0) for d in model.dispositions.values()] == [1182.0, 473.0, 236.0]

    # test that admissions are being properly calculated (thanks @PhilMiller)
    admissions = build_admits_df(param.n_days, model.dispositions)
    cumulative_admissions = admissions.cumsum()
    diff = cumulative_admissions["hospitalized"][1:-1] - (
        0.05 * 0.05 * (raw_df.infected[1:-1] + raw_df.recovered[1:-1]) - 100
    )
    assert (diff.abs() < 0.1).all()

def test_chart_descriptions(p=PARAM):
    # new admissions chart
    projection_admits = pd.read_csv('tests/projection_admits.csv')
    chart = new_admissions_chart(alt, projection_admits, p)
    description = chart_descriptions(chart, p.labels)

    hosp, icu, vent, asterisk = description.split("\n\n")  # break out the description into lines

    max_hosp = chart.data['hospitalized'].max()
    assert str(ceil(max_hosp)) in hosp

    max_icu_ix = chart.data['icu'].idxmax()
    assert max_icu_ix + 1 == len(chart.data)
    assert "*" in icu

    # test asterisk
    param = PARAM
    param.n_days = 600

    projection_admits = pd.read_csv('tests/projection_admits.csv')
    # projection_admits = projection_admits.rename(columns={'hospitalized': 'Hospitalized', 'icu': 'ICU', 'ventilated': 'Ventilated'})
    chart = new_admissions_chart(alt, projection_admits, p)
    description = chart_descriptions(chart, p.labels)
    assert "*" not in description

    # census chart
    census_df = pd.read_csv('tests/census_df.csv')
    # census_df = census_df.rename(columns={'hospitalized': 'Hospitalized', 'icu': 'ICU', 'ventilated': 'Ventilated'})
    PARAM.as_date = True
    chart = admitted_patients_chart(alt, census_df, p)
    description = chart_descriptions(chart, p.labels)

    assert str(ceil(chart.data['ventilated'].max())) in description
    assert str(chart.data['icu'].idxmax()) not in description
    assert datetime.datetime.strftime(chart.data.iloc[chart.data['icu'].idxmax()].date, '%b %d') in description

max_stars_repo_name: ellisbjohns/Chime

max_issues_repo_name: ellisbjohns/Chime

max_forks_repo_name: ellisbjohns/Chime


text:
from __future__ import division, print_function
from collections import OrderedDict
from itertools import chain
from functools import partial
import multiprocessing as mp
import os
import sys

try:
    import modshogun as sg
except ImportError:  # new versions just call it shogun
    import shogun as sg

import numpy as np
import pandas as pd
import progressbar as pb
from six.moves import xrange
from scipy.stats.mstats import mquantiles

import generate

if 'OMP_NUM_THREADS' in os.environ:
    num_threads = int(os.environ['OMP_NUM_THREADS'])
else:
    num_threads = mp.cpu_count()
sg.get_global_parallel().set_num_threads(num_threads)

def get_estimates(gen, sigmas=None, n_reps=100, n_null_samps=1000,
                  cache_size=64, rep_states=False, name=None,
                  save_samps=False, thresh_levels=(.2, .1, .05, .01)):
    if sigmas is None:
        sigmas = np.logspace(-1.7, 1.7, num=30)
    sigmas = np.asarray(sigmas)

    mmd = sg.QuadraticTimeMMD()
    mmd.set_num_null_samples(n_null_samps)
    mmd_mk = mmd.multikernel()
    for s in sigmas:
        mmd_mk.add_kernel(sg.GaussianKernel(cache_size, 2 * s**2))

    info = OrderedDict()
    for k in 'sigma rep mmd_est var_est p'.split():
        info[k] = []
    thresh_names = []
    for l in thresh_levels:
        s = 'thresh_{}'.format(l)
        thresh_names.append(s)
        info[s] = []
    if save_samps:
        info['samps'] = []

    thresh_prob = 1 - np.asarray(thresh_levels)

    bar = pb.ProgressBar()
    if name is not None:
        bar.start()
        bar.widgets.insert(0, '{} '.format(name))
    for rep in bar(xrange(n_reps)):
        if rep_states:
            rep = np.random.randint(0, 2**32)
            X, Y = gen(rs=rep)
        else:
            X, Y = gen()
        n = X.shape[0]
        assert Y.shape[0] == n
        mmd.set_p(sg.RealFeatures(X.T))
        mmd.set_q(sg.RealFeatures(Y.T))

        info['sigma'].extend(sigmas)
        info['rep'].extend([rep] * len(sigmas))

        stat = mmd_mk.compute_statistic()
        info['mmd_est'].extend(stat / (n / 2))

        samps = mmd_mk.sample_null()
        info['p'].extend(np.mean(samps >= stat, axis=0))
        if save_samps:
            info['samps'].extend(samps.T)

        info['var_est'].extend(mmd_mk.compute_variance_h1())

        threshes = np.asarray(mquantiles(samps, prob=thresh_prob, axis=0))
        for s, t in zip(thresh_names, threshes):
            info[s].extend(t)

    info = pd.DataFrame(info)
    info.set_index(['sigma', 'rep'], inplace=True)
    return info

def try_int(x):
    try:
        i = int(x)
    except (TypeError, ValueError):
        return x
    else:
        return i if i == x else x

def main():
    import argparse
    parser = argparse.ArgumentParser(description='''
        Evaluate different kernel selection methods for picking the bandwidth
        of an RBF kernel.

        Makes an hdf5 file containing results in OUT_FILE; see fixed_eval.py
        to make a jupyter notebook evaluating the results.
        '''.strip())

    data = parser.add_argument_group('Data')
    generate.add_problem_args(data)
    data.add_argument('-n', type=int, default=500)

    parser.add_argument('--n-reps', type=int, default=100)
    parser.add_argument('--n-null-samps', type=int, default=1000)

    parser.add_argument('--sigma-base', type=float, default=1)
    parser.add_argument('--sigma-max-mult', type=float, default=10**1.7)
    parser.add_argument('--n-sigmas', type=int, default=30)

    io = parser.add_argument_group('I/O')
    io.add_argument('out_file', nargs='?',
                    help="Defaults to res_fixed/PROBLEM/(d|rat)ARG/nN.h5.")
    io.add_argument('out_key', nargs='?', default='df')

    g = io.add_mutually_exclusive_group()
    g.add_argument('--save-samps', action='store_true', default=False)
    g.add_argument('--no-save-samps', action='store_false', dest='save_samps')

    io.add_argument('--levels', type=float, nargs='*',
                        default=[.2, .1, .05, .01])

    args = parser.parse_args()
    if args.out_file is None:
        name = next(chain(
            ('{}/d{}/n{}'.format(base, try_int(d), args.n) for base, d
              in [('sg', args.sg), ('gmd', args.gmd), ('gvd', args.gvd)]
              if d is not None),
            ('{}/rat{}/n{}'.format(base, try_int(r), args.n) for base, r
             in [('blobs', args.blobs)]
             if r is not None),
        ))
        args.out_file = os.path.join(
            os.path.dirname(__file__), 'res_fixed/{}.h5'.format(name))
    else:
        name = os.path.relpath(args.out_file)
        if args.out_key != 'df':
            name += ':{}'.format(args.out_key)

    d = os.path.dirname(args.out_file)
    if d and not os.path.isdir(d):
        os.makedirs(d)

    v = np.log10(args.sigma_max_mult)
    sigmas = np.logspace(-v, v, num=args.n_sigmas) * args.sigma_base

    gen = partial(generate.generate_data, args, args.n, dtype=np.float64)
    info = get_estimates(
        gen, n_reps=args.n_reps, n_null_samps=args.n_null_samps,
        rep_states=True, name=name, sigmas=sigmas,
        thresh_levels=args.levels, save_samps=args.save_samps)

    if args.save_samps:
        # for serialization, need to convert samps to a string :|
        info['samps'] = info.samps.map(lambda x: x.tostring())

    try:
        info.to_hdf(args.out_file, args.out_key,
                    format='table', complevel=9, complib='blosc', append=True)
    except Exception as e:
        if sys.stdin.isatty() and sys.stderr.isatty():
            print(e, file=sys.stderr)
            import IPython; IPython.embed()
        else:
            raise

if __name__ == '__main__':
    main()

max_stars_repo_name: Kensuke-Mitsuzawa/opt-mmd

max_issues_repo_name: Kensuke-Mitsuzawa/opt-mmd

max_forks_repo_name: Kensuke-Mitsuzawa/opt-mmd


text:
import numpy as np
import pandas as pd
import pickle as pkl
from collections import defaultdict
import os
import torch
import argparse

from tqdm.auto import tqdm, trange
from collections import defaultdict
from scipy.stats import pearsonr

import bert_score

wmt18_sys_to_lang_pairs = ['cs-en', 'de-en', 'et-en', 'fi-en', 'ru-en', 'tr-en', 'zh-en']
wmt18_sys_from_lang_pairs = ['en-cs', 'en-de', 'en-et', 'en-fi', 'en-ru', 'en-tr', 'en-zh']
wmt18_sys_lang_pairs = wmt18_sys_to_lang_pairs + wmt18_sys_from_lang_pairs

import argparse

def get_wmt18_seg_data(lang_pair):
    src, tgt = lang_pair.split('-')

    RRdata = pd.read_csv(
        "wmt18/wmt18-metrics-task-package/manual-evaluation/RR-seglevel.csv", sep=' ')
    # RRdata_lang = RRdata[RRdata['LP'] == lang_pair] # there is a typo in this data. One column name is missing in the header
    RRdata_lang = RRdata[RRdata.index == lang_pair]

    systems = set(RRdata_lang['BETTER'])
    systems.update(list(set(RRdata_lang['WORSE'])))
    systems = list(systems)
    sentences = {}
    for system in systems:
        with open("wmt18/wmt18-metrics-task-package/input/wmt18-metrics-task-nohybrids/system-outputs/newstest2018/{}/newstest2018.{}.{}".format(lang_pair, system, lang_pair)) as f:
            sentences[system] = f.read().split("\n")

    with open("wmt18/wmt18-metrics-task-package/input/wmt18-metrics-task-nohybrids/"
              "references/{}".format('newstest2018-{}{}-ref.{}'.format(src, tgt, tgt))) as f:
        references = f.read().split("\n")

    ref = []
    cand_better = []
    cand_worse = []
    for index, row in RRdata_lang.iterrows():
        cand_better += [sentences[row['BETTER']][row['SID']-1]]
        cand_worse += [sentences[row['WORSE']][row['SID']-1]]
        ref += [references[row['SID']-1]]

    return ref, cand_better, cand_worse

def kendell_score(scores_better, scores_worse):
    total = len(scores_better)
    correct = torch.sum(scores_better > scores_worse).item()
    incorrect = total - correct
    return (correct - incorrect)/total

def get_wmt18_seg_bert_score(lang_pair, scorer, cache=False, from_en=True, batch_size=64):
    filename = ''
    if from_en:
        if scorer.idf:
            filename = "cache_score/from_en/18/{}/wmt18_seg_from_{}_{}_idf.pkl".format(scorer.model_type, *lang_pair.split('-'))
        else:
            filename = "cache_score/from_en/18/{}/wmt18_seg_from_{}_{}.pkl".format(scorer.model_type, *lang_pair.split('-'))
    else:
        if scorer.idf:
            filename = "cache_score/to_en/18/{}/wmt18_seg_to_{}_{}_idf.pkl".format(scorer.model_type, *lang_pair.split('-'))
        else:
            filename = "cache_score/to_en/18/{}/wmt18_seg_to_{}_{}.pkl".format(scorer.model_type, *lang_pair.split('-'))
    if os.path.exists(filename):
        with open(filename, "rb") as f:
            return pkl.load(f)
    else:
        refs, cand_better, cand_worse = get_wmt18_seg_data(lang_pair)
        if scorer.idf:
            scorer.compute_idf(refs)
        scores_better = list(scorer.score(cand_better, refs, batch_size=batch_size))
        scores_worse = list(scorer.score(cand_worse, refs, batch_size=batch_size))
        if cache:
            os.makedirs(os.path.dirname(filename), exist_ok=True)
            with open(filename, "wb") as f:
                pkl.dump((scores_better, scores_worse), f)
        return scores_better, scores_worse

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("-d", "--data", default="wmt18", help="path to wmt16 data")
    parser.add_argument("-m", "--model", nargs="+", help="models to tune")
    parser.add_argument("-l", "--log_file", default="wmt18_log.csv", help="log file path")
    parser.add_argument("--idf", action="store_true")
    parser.add_argument("-b", "--batch_size", type=int, default=64)
    parser.add_argument(
        "--lang_pairs",
        nargs="+",
        default=wmt18_sys_to_lang_pairs,
        help="language pairs used for tuning",
    )
    args = parser.parse_args()

    torch.set_grad_enabled(False)

    header = 'model_type'
    for lang_pair in args.lang_pairs + ['avg']:
        header += f',{lang_pair}'
    print(header)
    if not os.path.exists(args.log_file):
        with open(args.log_file, 'w') as f:
            print(header, file=f)

    print(args.model)
    for model_type in args.model:
        scorer = bert_score.scorer.BERTScorer(model_type=model_type, idf=args.idf)
        results = defaultdict(dict)
        for lang_pair in tqdm(args.lang_pairs):
            scores_better, scores_worse = get_wmt18_seg_bert_score(lang_pair, scorer, batch_size=args.batch_size, cache=True, from_en=False)
            for sb, sw, name in zip(scores_better, scores_worse, ["P", "R", "F"]):
                results[lang_pair][f"{model_type} {name}"] = kendell_score(sb, sw)

        for name in ["P", "R", "F"]:
            temp = []
            for lang_pair in args.lang_pairs:
                temp.append(results[lang_pair][f"{model_type} {name}"])
            results["avg"][f"{model_type} {name}"] = np.mean(temp)

            msg = f"{model_type} {name} (idf)" if args.idf else f"{model_type} {name}"
            for lang_pair in args.lang_pairs + ['avg']:
                msg += f",{results[lang_pair][f'{model_type} {name}']}"
            print(msg)
            with open(args.log_file, "a") as f:
                print(msg, file=f)

        del scorer

if __name__ == "__main__":
    main()


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

# ----------------------------------------------------------------------------
#
# TITLE   : illustrate_actionangle_script
# PROJECT : Pal 5 update MW pot constraints
#
# ----------------------------------------------------------------------------

# Docstring
"""Illustrate Actions & Angles.

Make a little movie to illustrate how to compute actions.

Routing Listings
----------------
create_frames
create_movie

Parameters
----------
basefilename : str
    first used to create frames, which are individual plots
    next used to make the movie, from the frames

outputfilename : str
    the location of the movie

Examples
--------
In general
>>> python illustrate_actionangle.py basefilename outputfilename

A specific example
>>> python illustrate_actionangle.py ../movies/aaI ../actionIllustration.mpg

"""

__author__ = "Jo Bovy"
__copyright__ = "Copyright 2016, 2020, "
__maintainer__ = "Nathaniel Starkman"

__all__ = ["create_frames", "create_movie"]

###############################################################################
# IMPORTS

# GENERAL
import sys
import subprocess
import numpy
from tqdm import tqdm

import astropy.units as u

from matplotlib import pyplot, gridspec

# CUSTOM
from galpy.potential import IsochronePotential
from galpy.orbit import Orbit
from galpy.actionAngle import actionAngleIsochrone
from galpy.util import bovy_plot

###############################################################################
# PARAMETERS

_RO, _VO = 8.0, 220.0

###############################################################################
# CODE
###############################################################################

def create_frames(basefilename):
    """Create Frames.

    Parameters
    ----------
    basefilename : str

    Returns
    -------
    None

    """
    pot = IsochronePotential(normalize=1.0, b=1.2)
    aAT = actionAngleIsochrone(ip=pot)
    bmock = 0.8
    aAF = actionAngleIsochrone(b=bmock)

    o = Orbit([1.0, 0.3, 0.9, 0.2, 1.0, 0.1])
    tfac, skip = 10, 1
    ndesired = 30 * 25
    times = numpy.linspace(0.0, 7.0, tfac * ndesired) * u.Gyr

    # Subsample
    otimesIndices = (numpy.arange(len(times)) / float(len(times) - 1)) ** 10 * (
        len(times) - 2
    )
    otimesIndices = numpy.unique(numpy.floor(otimesIndices).astype("int"))
    if len(otimesIndices) > ndesired:
        sfac = int(numpy.floor(len(otimesIndices) / float(ndesired)))
        otimesIndices = otimesIndices[::sfac]
    otimes = times[otimesIndices]
    o.integrate(times, pot)

    # Compute all actions in the wrong potential
    acfs = aAF.actionsFreqsAngles(
        o.R(times), o.vR(times), o.vT(times), o.z(times), o.vz(times), o.phi(times),
    )
    js = (acfs[0], acfs[1], acfs[2])
    danglerI = ((numpy.roll(acfs[6], -1) - acfs[6]) % (2.0 * numpy.pi))[:-1]
    jrC = numpy.cumsum(acfs[0][:-1] * danglerI) / numpy.cumsum(danglerI)

    # And also the actions in the true potential
    jsT = aAT(o.R(times), o.vR(times), o.vT(times), o.z(times), o.vz(times))
    jrT = numpy.median(jsT[0]) * _RO * _VO

    # ---------------------------------------------------------------
    # Plotting

    # Setup gridspec
    gs = gridspec.GridSpec(1, 3, wspace=0.325, bottom=0.2, left=0.075)

    # For each time step, plot: orbit, Jr, <Jr>
    for ii in tqdm(range(len(otimes))):
        bovy_plot.bovy_print(
            fig_width=11.2,
            fig_height=4.0,
            axes_labelsize=17.0,
            text_fontsize=12.0,
            xtick_labelsize=13.0,
            ytick_labelsize=13.0,
        )
        pyplot.figure()
        pyplot.subplot(gs[0])
        minIndx = otimesIndices[ii] - 100

        if minIndx < 0:
            minIndx = 0
        bovy_plot.bovy_plot(
            [o.x(otimes[ii : ii + 1]) * _RO],
            [o.z(otimes[ii : ii + 1]) * _RO],
            "o",
            ms=15.0,
            gcf=True,
            xrange=[-19.0, 19.0],
            yrange=[-19.0, 19.0],
            xlabel=r"$X\,(\mathrm{kpc})$",
            ylabel=r"$z\,(\mathrm{kpc})$",
        )

        if ii > 0:
            bovy_plot.bovy_plot(
                o.x(times[minIndx : otimesIndices[ii] : skip]) * _RO,
                o.z(times[minIndx : otimesIndices[ii] : skip]) * _RO,
                "-",
                overplot=True,
            )

        pyplot.subplot(gs[1])
        bovy_plot.bovy_plot(
            [times[otimesIndices[ii]].value],
            [js[0][otimesIndices[ii]] * _RO * _VO],
            "o",
            ms=15.0,
            gcf=True,
            xrange=[0.0, 1.0 + times[otimesIndices[ii]].value],
            yrange=[0.0, 349.0],
            xlabel=r"$\mathrm{time}\,(\mathrm{Gyr})$",
            ylabel=r"$J_R\,(\mathrm{kpc\,km\,s}^{-1})$",
        )

        if ii > 0:
            bovy_plot.bovy_plot(
                times[: otimesIndices[ii] : skip].value,
                js[0][: otimesIndices[ii] : skip] * _RO * _VO,
                "-",
                overplot=True,
            )

        pyplot.axhline(jrT, ls="--", color="k")
        pyplot.subplot(gs[2])
        bovy_plot.bovy_plot(
            [otimes[ii].value],
            [jrC[otimesIndices[ii]] * _RO * _VO],
            "o",
            ms=15.0,
            gcf=True,
            xrange=[0.0, 1.0 + times[otimesIndices[ii]].value],
            yrange=[
                (otimes[ii] / times[-1]) ** 0.3 * (jrT - 3),
                349.0 + (otimes[ii] / times[-1]) ** 0.3 * (jrT + 3 - 349.0),
            ],
            xlabel=r"$\mathrm{time}\,(\mathrm{Gyr})$",
            ylabel=r"$J_R\,(\mathrm{kpc\,km\,s}^{-1})$",
        )

        if ii > 0:
            bovy_plot.bovy_plot(
                times[: otimesIndices[ii] : skip].value,
                jrC[: otimesIndices[ii] : skip] * _RO * _VO,
                "-",
                overplot=True,
            )

        pyplot.axhline(jrT, ls="--", color="k")
        bovy_plot.bovy_end_print(basefilename + "_%05d.png" % ii)

    return None

# /def

# --------------------------------------------------------------------------

def create_movie(basefilename, outputfilename):
    """Create Movie.

    Parameters
    ----------
    basefilename : str
    outputfilename : str

    Returns
    -------
    None

    """
    framerate = 25
    bitrate = 1000000

    try:
        subprocess.check_call(
            [
                "ffmpeg",
                "-i",
                basefilename + "_%05d.png",
                "-y",
                "-framerate",
                str(framerate),
                "-r",
                str(framerate),
                "-b",
                str(bitrate),
                outputfilename,
            ]
        )

    except subprocess.CalledProcessError:
        print("'ffmpeg' failed")

    return None

# /def

###############################################################################
# Command Line
###############################################################################

if __name__ == "__main__":

    create_frames(sys.argv[1])

    create_movie(sys.argv[1], sys.argv[2])

# /if

###############################################################################
# END

reason: import numpy,import astropy


text:
from transformers import BertForTokenClassification
import torch
from transformers import BertTokenizer
import numpy as np
import nltk.data
nltk.download('punkt')

import torch
from torch.utils.data import TensorDataset, DataLoader, RandomSampler, SequentialSampler
from transformers import BertTokenizer, BertConfig, AutoModelForTokenClassification, AutoConfig

from keras.preprocessing.sequence import pad_sequences
from sklearn.model_selection import train_test_split

from transformers import BertForTokenClassification, AdamW

from transformers import get_linear_schedule_with_warmup

from seqeval.metrics import accuracy_score
from sklearn.metrics import f1_score, classification_report, precision_score, recall_score

import torch.nn as nn

from tqdm import trange
import numpy as np
import matplotlib.pyplot as plt

from nltk.tokenize import sent_tokenize

import os

##########################################################
# import wandb
# from transformers import TrainingArguments, Trainer

# wandb.init(project="project", entity="3rd_year_project")
##########################################################

class NER_BERT(object):
    device = torch.device("cuda")
    #device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

    tag2idx = {'O':0, 'ID':1, 'PHI':2, 'NAME':3, 'CONTACT':4, 'DATE':5, 'AGE':6, 'PROFESSION':7, 'LOCATION':8, 'PAD': 9}
    tag_values = ["O","ID", "PHI", "NAME", "CONTACT", "DATE", "AGE", "PROFESSION", "LOCATION", "PAD"]

    tokenizer = BertTokenizer.from_pretrained('bert-base-cased', num_labels=len(tag2idx), do_lower_case=False)

    MAX_LEN = 75 # max length of sequence, needs for padding
    bs = 32 # batch size

    """Abstract class that other NER plugins should implement"""
    def __init__(self):

        #config = AutoConfig.from_pretrained('ArishkaBelovishka/bert-i2b2')
        #self.model = AutoModelForTokenClassification.from_pretrained('ArishkaBelovishka/bert-i2b2', config = config)

        # Uncomment the following if you want to load your fine-tuned model from Models folder.
        # If you just want to run NER use hugging-face repository where fine-tuned on half of i2b2 data model lives.

        if os.path.exists("Models/BERT_epoch-10.pt"):
            print("Loading model")
            state_dict = torch.load("Models/BERT_epoch-10.pt", map_location=torch.device('cuda'))
            print("Loaded model")

            self.model = BertForTokenClassification.from_pretrained(
                "bert-base-cased",
                state_dict = state_dict,
                num_labels=len(NER_BERT.tag2idx),
                output_attentions = True,
                output_hidden_states = True
            )
        else:
            self.model = BertForTokenClassification.from_pretrained(
                "bert-base-cased",
                num_labels=len(NER_BERT.tag2idx),
                output_attentions = True,
                output_hidden_states = True
            )

    def perform_NER(self,text):
        """Implementation of the method that should perform named entity recognition"""
        # tokenizer to divide data into sentences (thanks, nltk)

        list_of_sents = sent_tokenize(text)

        list_of_tuples_by_sent = []

        for sent in list_of_sents:
            # , truncation=True
            tokenized_sentence = self.tokenizer.encode(sent, truncation=True) # BERT tokenizer is clever, it will internally divide the sentence by words, so all we need to provide there is sentence and it will return an array where each token is either special token/word/subword, refer to BERT WordPiece tokenizer approach
            # truncation=True to comply with 512 length of the sentence
            input_ids = torch.tensor([tokenized_sentence])

            with torch.no_grad():
                # Run inference/classification
                output = self.model(input_ids)
            label_indices = np.argmax(output[0].to("cpu").numpy(), axis=2)
            tokens = self.tokenizer.convert_ids_to_tokens(input_ids.to('cpu').numpy()[0])
            new_tokens, new_labels = [], []
            for token, label_idx in zip(tokens, label_indices[0]):
                if token.startswith("##"):
                    new_tokens[-1] = new_tokens[-1] + token[2:]
                else:
                    new_labels.append(self.tag_values[label_idx])
                    new_tokens.append(token)

            list_of_tuples = []
            for token, label in zip(new_tokens, new_labels):
                list_of_tuples.append((token, label))
                #print("{}\t{}".format(label, token))

            list_of_tuples_by_sent.append(list_of_tuples)

        # remove [CLS] and [SEP] tokens to comply wth xml structure
        for i in range(len(list_of_tuples_by_sent)):
            for tag in self.tag_values:
                if ('[CLS]', tag) in list_of_tuples_by_sent[i]:
                    list_of_tuples_by_sent[i].remove(('[CLS]', tag))

                if ('[SEP]', tag) in list_of_tuples_by_sent[i]:
                    list_of_tuples_by_sent[i].remove(('[SEP]', tag))

        return list_of_tuples_by_sent

    # Needed for transform_sequences
    def tokenize_and_preserve_labels(self, sentence, text_labels):
        tokenized_sentence = []
        labels = []
        for word, label in zip(sentence, text_labels):
            # Tokenize the word and count # of subwords the word is broken into
            tokenized_word = NER_BERT.tokenizer.tokenize(word)
            n_subwords = len(tokenized_word)
            # Add the tokenized word to the final tokenized word list
            tokenized_sentence.extend(tokenized_word)
            # Add the same label to the new list of labels `n_subwords` times
            labels.extend([label] * n_subwords)
        return tokenized_sentence, labels

    def transform_sequences(self,tokens_labels):
        """method that transforms sequences of (token,label) into feature sequences. Returns two sequence lists for X and Y"""
        print("I am in transform seq")
        # result - one document, result[i] is sentence in document, result [i][i] is word in sentence
        tokenized_sentences = []
        labels = []
        for index, sentence in enumerate(tokens_labels):
            text_labels = []
            sentence_to_feed = []
            for word_label in sentence:
                text_labels.append(word_label[1])
                sentence_to_feed.append(word_label[0])
            a, b = self.tokenize_and_preserve_labels(sentence_to_feed, text_labels)
            tokenized_sentences.append(a)
            labels.append(b)

        # Now need to split long tokenized sequences into subsequences of length less than 512 tokens
        # not to loose valuable information in NER, basically not to cut sentences
        # i2b2 docs are very ugly and sentences in them are usually way too long as doctors forget to put full stops...
        # tokenized_sentences AND labels are the same strucutre of 2d arrays

        # I need to take care of the issue if I am going to split beginning of the word and its end, like
        # Arina is tokenized as "Ari" and "##na", thus I cannot separate the two, otherwise it will not make sense

        distributed_tokenized_sentences, distributed_labels = [], []
        for sent, label in zip(tokenized_sentences, labels):
            if len(sent) > NER_BERT.MAX_LEN:
                while len(sent) > NER_BERT.MAX_LEN:
                    #print("I am in while loop to truncate sequence")
                    index = NER_BERT.MAX_LEN - 2
                    for i in range(NER_BERT.MAX_LEN - 2, 0, -1):
                        if sent[i][:2] == "##":
                            index = index - 1
                        else:
                            break

                    new_sent = sent[:index] # 511 because we want to append [SEP] token in the end
                    new_label = label[:index]

                    sent = sent[index:]  # update given sent
                    label = label[index:]

                    distributed_tokenized_sentences.append(new_sent)
                    distributed_labels.append(new_label)

                distributed_tokenized_sentences.append(sent)
                distributed_labels.append(label)
                #print(sent)

            else:
                distributed_tokenized_sentences.append(sent)
                distributed_labels.append(label)

        input_ids = pad_sequences([NER_BERT.tokenizer.convert_tokens_to_ids(txt) for txt in distributed_tokenized_sentences],
                                maxlen=NER_BERT.MAX_LEN, dtype="long", value=0.0,
                                truncating="post", padding="post")

        tags = pad_sequences([[NER_BERT.tag2idx.get(l) for l in lab] for lab in distributed_labels],
                            maxlen=NER_BERT.MAX_LEN, value=NER_BERT.tag2idx["PAD"], padding="post",
                            dtype="long", truncating="post")

        # Result is pair X (array of sentences, where each sentence is an array of words) and Y (array of labels)
        return input_ids, tags

    def learn(self, X_train,Y_train, epochs=1):
        """Function that actually train the algorithm"""
        # if torch.cuda.is_available():
        #     self.model.cuda()

        tr_masks = [[float(i != 0.0) for i in ii] for ii in X_train]

        print("READY TO CREATE SOME TENZORS!!!!!!!!!!!!!!!!!!!!!!!!!!")
        tr_inputs = torch.tensor(X_train).type(torch.long)
        tr_tags = torch.tensor(Y_train).type(torch.long)
        tr_masks = torch.tensor(tr_masks).type(torch.long)

        train_data = TensorDataset(tr_inputs, tr_masks, tr_tags)
        train_sampler = RandomSampler(train_data)
        train_dataloader = DataLoader(train_data, sampler=train_sampler, batch_size=NER_BERT.bs)

        print("READY TO PREPARE OPTIMIZER!!!!!!!!!!!!!!!!!!!!!!!!!!")

        # Weight decay in Adam optimiser (adaptive gradient algorithm) is a regularisation technique which is extensively disucssed in this paper:
        # https://arxiv.org/abs/1711.05101
        # (Like L2 for SGD but different)
        # resularisation of the model objective function in order to prevent overfitting of the model.
        FULL_FINETUNING = True
        if FULL_FINETUNING:
            param_optimizer = list(self.model.named_parameters())
            no_decay = ['bias', 'gamma', 'beta']
            optimizer_grouped_parameters = [
                {'params': [p for n, p in param_optimizer if not any(nd in n for nd in no_decay)],
                'weight_decay_rate': 0.01}, # in AdamW implementation (default: 1e-2)
                {'params': [p for n, p in param_optimizer if any(nd in n for nd in no_decay)],
                'weight_decay_rate': 0.0}
            ]
        else:
            param_optimizer = list(self.model.classifier.named_parameters())
            optimizer_grouped_parameters = [{"params": [p for n, p in param_optimizer]}]

        # TODO: change to new implementation of AdamW: torch.optim.AdamW(...)
        optimizer = AdamW(
            optimizer_grouped_parameters,
            lr=3e-5,
            eps=1e-8
        )

        max_grad_norm = 1.0

        # Total number of training steps is number of batches * number of epochs.
        total_steps = len(train_dataloader) * epochs

        # Create the learning rate scheduler.
        # We need it to adjust learning rate if the accuracy does not change between epochs much,
        # basically pushing the model to learn.
        # https://sajjjadayobi.github.io/blog/markdown/2021/05/23/adamw-warmup.html
        scheduler = get_linear_schedule_with_warmup(
            optimizer,
            num_warmup_steps=0,
            num_training_steps=total_steps
        )

        print("START TRAINING!!!!!!!!!!!!!!!!!!!!!!!!!!!!")
        ## Store the average loss after each epoch so we can plot them.
        loss_values, validation_loss_values = [], []

        # just for intermediate model save naming
        epoch_num = 3

        for _ in trange(epochs, desc="Epoch"):
            # ========================================
            #               Training
            # ========================================
            # Perform one full pass over the training set.
            # clean the cache not to fail with video memory
            # if torch.cuda.is_available():
            #     torch.cuda.empty_cache()

            # just for intermediate model save naming
            epoch_num += 1

            # Put the model into training mode.
            self.model.train()
            # Reset the total loss for this epoch.
            total_loss = 0

            print("Start backprop and optimisation!!! Epoch has passed!!!!!!!!!!!!!!!!!!!!!!!")

            # Training loop
            for step, batch in enumerate(train_dataloader):
                print("We are in the batch!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!")

                # add batch to gpu
                batch = tuple(b.to(NER_BERT.device) for b in batch)
                b_input_ids, b_input_mask, b_labels = batch
                # Always clear any previously calculated gradients before performing a backward pass.
                self.model.zero_grad()
                # forward pass
                # This will return the loss (rather than the model output)
                # because we have provided the `labels`.
                outputs = self.model(b_input_ids, token_type_ids=None,
                                attention_mask=b_input_mask, labels=b_labels)
                # get the loss
                loss = outputs[0]

                # Perform a backward pass to calculate the gradients.
                loss.backward()
                # track train loss
                total_loss += loss.item()
                # Clip the norm of the gradient
                # This is to help prevent the "exploding gradients" problem.
                torch.nn.utils.clip_grad_norm_(parameters=self.model.parameters(), max_norm=max_grad_norm)
                # update parameters
                optimizer.step()
                # Update the learning rate.
                scheduler.step()
                print("We processed one batch!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!")

            # Calculate the average loss over the training data.
            avg_train_loss = total_loss / len(train_dataloader)
            print("Average train loss: {}".format(avg_train_loss))

            # Store the loss value for plotting the learning curve.
            loss_values.append(avg_train_loss)

            # Save intermediate weights of the model, i.e. if computer goes crazy and drops the training or you
            # want to test the performance from different epochs
            torch.save(self.model.state_dict(), os.path.join("Models_intermediate/", 'BERT_epoch-{}.pt'.format(epoch_num)))

                #Plot the learning curve.
        plt.figure()
        plt.plot(loss_values, 'b-o', label="training loss")
        # Label the plot.
        plt.title("Learning curve")
        plt.xlabel("Epoch")
        plt.ylabel("Loss")
        plt.legend()

        plt.show()

    def evaluate(self, X_test,Y_test):
        if torch.cuda.is_available():
            self.model.cuda()
        """Function to evaluate algorithm"""
        val_masks = [[float(i != 0.0) for i in ii] for ii in X_test]
        val_inputs = torch.tensor(X_test).type(torch.long)
        val_tags = torch.tensor(Y_test).type(torch.long)
        val_masks = torch.tensor(val_masks).type(torch.long)

        valid_data = TensorDataset(val_inputs, val_masks, val_tags)
        valid_sampler = SequentialSampler(valid_data)
        valid_dataloader = DataLoader(valid_data, sampler=valid_sampler, batch_size=NER_BERT.bs)
        # seed
        # for _ in range(2):
        #valid_dataloader = DataLoader(valid_data, shuffle=True, batch_size=NER_BERT.bs)

        # for one random seed of valid_dataloader:
        # ...

        # ========================================
        #               Validation
        # ========================================
        # After the completion of each training epoch, measure our performance on
        # our validation set.

        # Put the model into evaluation mode to set dropout and batch normalization layers to evaluation mode to have consistent results
        self.model.eval()
        # Reset the validation loss for this epoch.
        eval_loss, eval_accuracy = 0, 0
        nb_eval_steps, nb_eval_examples = 0, 0
        predictions , true_labels = [], []
        for batch in valid_dataloader:
            batch = tuple(t.to(self.device) for t in batch)
            b_input_ids, b_input_mask, b_labels = batch

            # Telling the model not to compute or store gradients,
            # saving memory and speeding up validation
            with torch.no_grad():
                # Forward pass, calculate logit predictions.
                # This will return the logits rather than the loss because we have not provided labels.
                outputs = self.model(b_input_ids, token_type_ids=None,
                                attention_mask=b_input_mask, labels=b_labels)
            # Move logits and labels to CPU
            logits = outputs[1].detach().cpu().numpy()
            label_ids = b_labels.to('cpu').numpy()

            # Calculate the accuracy for this batch of test sentences.
            eval_loss += outputs[0].mean().item()
            predictions.extend([list(p) for p in np.argmax(logits, axis=2)])
            true_labels.extend(label_ids)

        eval_loss = eval_loss / len(valid_dataloader)
        print("Validation loss: {}".format(eval_loss))
        pred_tags = [NER_BERT.tag_values[p_i] for p, l in zip(predictions, true_labels)
                                    for p_i, l_i in zip(p, l) if NER_BERT.tag_values[l_i] != "PAD"]

        ###############################################################################
        # reconstruct given text for purposes of algorithms' performance comparison
        # our X_test is again a list of sentences, i.e. 2d array
        tokens = [self.tokenizer.convert_ids_to_tokens(sent) for sent in X_test]
        # Unpack tokens into 1d array to be able to go through it with labels
        # [PAD] and not just PAD because that is what BERT actually puts
        tokens_flat = [item for sublist in tokens for item in sublist if item != "[PAD]"]

        #for sentence in tokens:
        new_tokens, new_labels = [], []
        for token, pred in zip(tokens_flat, pred_tags):
            #print("{}\t{}".format(token, pred))
            if token.startswith("##"):
                new_tokens[-1] = new_tokens[-1] + token[2:]
            else:
                new_labels.append(pred)
                new_tokens.append(token)
        ###############################################################################

        valid_tags = [NER_BERT.tag_values[l_i] for l in true_labels
                                    for l_i in l if NER_BERT.tag_values[l_i] != "PAD"]

        print("Validation Accuracy: {}".format(accuracy_score(valid_tags, pred_tags)))  # was other way around, why?
        print("Validation F1-Score: {}".format(f1_score(valid_tags, pred_tags, average='weighted'))) # correct
        print("Validation precision: {}".format(precision_score(valid_tags, pred_tags, average='weighted')))
        print("Validation recall: {}".format(recall_score(valid_tags, pred_tags, average='weighted')))
        labels = ["ID", "PHI", "NAME", "CONTACT", "DATE", "AGE",
                "PROFESSION", "LOCATION"]
        print(classification_report(valid_tags, pred_tags, digits=4, labels=labels))
        print()

        ###############################################################################

        # to evaluate union/intersection of algorithms
        # for t, l in zip(new_tokens, new_labels):
        #     print("{}\t{}".format(t, l))
        return new_labels

        # # Use plot styling from seaborn.
        # sns.set(style='darkgrid')

        # # Increase the plot size and font size.
        # sns.set(font_scale=1.5)
        # plt.rcParams["figure.figsize"] = (12,6)

        # # Plot the learning curve.
        # plt.plot(loss_values, 'b-o', label="training loss")
        # plt.plot(validation_loss_values, 'r-o', label="validation loss")

        # # Label the plot.
        # plt.title("Learning curve")
        # plt.xlabel("Epoch")
        # plt.ylabel("Loss")
        # plt.legend()

        # plt.show()

    def save(self, model_path):
        """
        Function to save model. Models are saved as h5 files in Models directory. Name is passed as argument
        :param model_path: Name of the model file
        :return: Doesn't return anything
        """
        torch.save(self.model.state_dict(), "Models/"+model_path+".pt")
        print("Saved model to disk")