text:
# -*- coding: utf-8 -*-
"""
Created on Mon May 29 16:04:03 2017

@author: Kristen
"""
import math
import numpy

#finds the distance between two atoms
def distance(position1, position2):
    return math.sqrt(math.pow(position1[0] - position2[0],2) +
                     math.pow(position1[1] - position2[1],2) +
                     math.pow(position1[2] - position2[2],2))

#finds if a triplet could have an Si atom between them
def dists(positions, dist):

    #if there were not enough close to make a triplet, return none
    if len(positions)<3:
        return[""]
    #if there is a triplet and they are close enough to have a Si,
    #return the triplet, else return blank
    if len(positions)==3:
        if distance(positions[1],positions[2])<=dist:
            return positions
        else:
            return[""]
    numbers=[]

    if len(positions)==5:
        print(1)
    #if there are more then 2 close enough to have a Si between them,findthe
    # one that could not given the other two
    for i in range(len(positions)):
        numbers.append(0)
    for i in range(1,len(positions)-1):
        for j in range(1,len(positions)-i):
            #if two positions are not close enough, add a counter to both. If they
            # are close enough, remove a counter from both
            if distance(positions[i],positions[i+j])>dist:
                numbers[i]+=1
                numbers[i+j]+=1
            else:
                numbers[i]-=1
                numbers[i+j]-=1

    #removetheonewiththemostcounters
    del positions[numbers.index(max(numbers))]

    #if these still are not close enough to have a triplet between them, return
    # none. If they are close enough, return the new triplet
    if distance(positions[1],positions[2])<=dist:
        return positions
    else:
        return[""]

#finds four membered rings and returns a list of lists of their locaitons
def find_four(opositions, far):
    rings=[[]]
    remov=[]
    #for each oxygen
    for i in range(len(opositions)):
        rings.append([""])
        rings[i]=[opositions[i]]
        #for each oxygen with an x position higher than the current
        for j in range(1,len(opositions)-i):
            #if th exposition is less than the possible distance between two
            # oxygenatoms(variableinclusionradius)
            if abs(opositions[i][0]-opositions[i+j][0])<=far:
                #if the distance between the two oxygens is less than the
                # characteristic distance(variable inclusion radius)
                if distance(opositions[i],opositions[i+j])<=far:
                    rings[i].append(opositions[i+j])
        rem=0
        if len(rings[i])<4:
            rem=1
        elif len(rings[i])>4:
            while len(rings[i])!=4:
                distances=[]
                for k in range(len(rings[i])):
                    tot_len = 0
                    for l in range(1,len(rings[i])-k):
                        tot_len+=distance(rings[i][k],rings[i][k+l])
                    distances.append(tot_len)
                del rings[i][distances.index(max(distances))]
        if len(rings[i])==4:
            distances=[]
            for n in range(len(rings[i])-1):
                for m in range(1,len(rings[i])-n):
                    distances.append(distance(rings[i][n],rings[i][n+m]))
            for n in range(2):
                del distances[distances.index(max(distances))]
            for n in range(4):
                for m in range(1,len(distances)-n):
                    if abs(distances[n]-distances[n+m])>.03:
                        rem=1
        if rem==1:
            remov.insert(0,i)

    for n in range(len(remov)):
        del rings[remov[n]]

    return rings

#finds the area of triangle
def triarea(p1,p2,p3):
    a=distance(p1,p2)
    b=distance(p2,p3)
    c=distance(p1,p3)
    s=(a+b+c)/2
    return math.sqrt(s*(s-a)*(s-b)*(s-c))

def ringarea(corners):
    n=len(corners)
    area=0.0
    for i in range(n):
      j=(i+1)%n
      area += corners[i][0]*corners[j][1]
      area -= corners[j][0]*corners[i][1]
    area=abs(area)/2.0
    return float(area)

#finds if the silicon atom is within a 4 membered ring
def rem4(rings,si):
    deletes=[]
    for i in range(len(rings)):
        triangles=0
        distances=[]
        locations=[]
        for n in range(len(rings[i])-1):
            for m in range(1,len(rings[i])-n):
                distances.append(distance(rings[i][n],rings[i][n+m]))
                locations.append([n,n+m])
        locations.append(len(rings[i]))
        for n in range(2):
            del locations[distances.index(max(distances))]
            del distances[distances.index(max(distances))]
        for n in range(len(locations)):
            triangles += triarea(rings[i][locations[n][0]],
                                 rings[i][locations[n][1]],si)
        if ringarea(rings[i])==triangles:
            return"n"
    return"y"

#finds the position of a Si given a triplet of oxygen
def si_finder(opositions):

    #characteristic distance
    dist=1.6*math.pow(10,-1)

    #sets up the translation to happen around a basepoint(the first point in
    # the positions)
    trans=[[0,0,0],[opositions[1][0]-opositions[0][0],
            opositions[1][1]-opositions[0][1],
            opositions[1][2]-opositions[0][2]],
    [opositions[2][0]-opositions[0][0],opositions[2][1]-opositions[0][1],
     opositions[2][2]-opositions[0][2]]]

    #finds vector perpendicular to the plane of the three points
    v=numpy.matrix([numpy.linalg.det([[trans[1][1],trans[2][1]],
                                      [trans[1][2],trans[2][2]]]),
    numpy.linalg.det([[trans[1][0],trans[2][0]],[trans[1][2],
                       trans[2][2]]]),
    numpy.linalg.det([[trans[1][0],trans[2][0]],[trans[1][1],trans[2][1]]])])

    #sets up first rotation matrix about the x axis
    theta=math.atan2(v.item(1),v.item(2))
    xmatr=numpy.matrix([[1,0,0],[0,math.cos(theta),-math.sin(theta)],
                        [0,math.sin(theta),math.cos(theta)]])
    trans1=numpy.matrix(trans)
    rot1=numpy.matrix.dot(trans1,xmatr)
    v1=numpy.matrix.dot(v,xmatr)

    #second rotation matrix about the y axis
    rho=math.atan2(v1.item(0),v1.item(2))
    ymatr=numpy.matrix([[math.cos(rho),0,math.sin(rho)],[0,1,0],
                         [-math.sin(rho),0,math.cos(rho)]])
    rot2=numpy.matrix.dot(rot1,ymatr)

    #should be in the xy plane now. Have to rotate such that two points
    # are on the x axis
    alph=math.atan2(rot2.item(4),rot2.item(3))
    bet=math.atan2(rot2.item(7),rot2.item(6))
    r1=math.sqrt(math.pow(rot2.item(3),2)+math.pow(rot2.item(4),2))
    r2=math.sqrt(math.pow(rot2.item(6),2)+math.pow(rot2.item(7),2))
    rot3=numpy.matrix([[rot2.item(0),rot2.item(1),rot2.item(2)],[r1,0,0],
                        [r2*math.cos(bet-alph),r2*math.sin(bet-alph),0]])
    x=r1/2
    y=r2*(1-math.cos(bet-alph))/(2*math.sin(bet-alph))
    z=math.sqrt(abs(math.pow(dist,2)-math.pow(x,2)-math.pow(y,2)))
    si_pos=numpy.matrix([x,y,z])

    #rotate back to originial position
    init=math.atan2(si_pos.item(1),si_pos.item(0))
    r=math.sqrt(math.pow(si_pos.item(0),2)+math.pow(si_pos.item(1),2))
    x=r*math.cos(init+alph)
    y=r*math.sin(init+alph)
    si_pos=numpy.matrix([x,y,z])

    #undo second rotation matrix
    iymatr=numpy.linalg.inv(ymatr)
    si_pos=numpy.matrix.dot(si_pos,iymatr)

    #undo first rotation matrix
    ixmatr=numpy.linalg.inv(xmatr)
    si_pos=numpy.matrix.dot(si_pos,ixmatr)

    #translate back so there is no point at the origin
    si_pos=[si_pos.item(0)+opositions[0][0],si_pos.item(1)+opositions[0][1],
            si_pos.item(2)+opositions[0][2]]

    return si_pos

#locates all possiblee triplets
def o_locator(opositions):

    dist = 1.6*math.pow(10,-1)

    #assumed oxygens are ordered by increasing x values
    #used to collect all the found oxygens close enough to have a single Si
    # between them
    found=[[""]]
    #for each oxygen
    for i in range(len(opositions)):
        found[i]=[opositions[i]]
        #for each oxygen with an x position higher than the current
        for j in range(1,len(opositions)-i):
            #if the x position is less than the possible distance between two
            # oxygenatoms(variableinclusionradius)
            if abs(opositions[i][0]-opositions[i+j][0])<=3.45*math.pow(10,-1):
            #if the distance between the two oxygens is less than the
                # characteristic distance(variable inclusion radius)
                if distance(opositions[i],opositions[i+j])<=3.45*math.pow(10,-1):
                    found[i].append(opositions[i+j])
        found.append([""])

    #removes last appended empty list
    del found[len(found)-1]

    #remove all those too far apart using dist function (variable inclusion
    # radius)
    for n in range(len(found)):
        found[n]=dists(found[n], .345)

    #createanarrayforpositionstoremove
    remov=[]
    #for all atoms with found oxygens
    for n in range(len(found)):
        #add empties to a list for removal
        if found[n]==[""]:
            remov.insert(0,n)

    #remove those in the remove list
    for m in range(len(remov)):
        del found[remov[m]]

    #return the list of those oxygen that have a possible Si between them
    return found

def main():

    #convert data in file into floats and append to a position list
    with open("OfC Positions 120106_008 Python Output_2.txt") as f:
        content=f.readline()

    string = ""

    locations = []

    for i in range(len(content)):
        if content[i] == " ":
            locations.append(float(string))
            string = ""
        else:
            string += content[i]

    locations.append(float(string))

    positions=[[""]]

    for i in range(len(locations)):
        if i%3==0:
            positions[i/3]=[locations[i]]
            positions.append("")
        else:
            positions[i/3].append(locations[i])

    del positions[len(positions)-1]

    #find triplets
    triples=o_locator(positions)

     #find Si positions
    si_locations=[]
    for j in range(len(triples)):
         si_locations.append(si_finder(triples[j]))

     #rings = find four(positions,.35)

    delete=[]

     #for i in range(len(silocations)):
         #if rem4(rings,silocations[i])=="n":
             #delete.append(0,i)

    for i in range(len(delete)):
        del si_locations[delete[i]]

    print(si_locations)

    #write Si positions to an outfile
    out = open("test.txt", "w")
    out.write(str(si_locations))
    out.write("\n")

if __name__=="__main__": main()

max_stars_repo_path: Misc Parts/Silicon from Oxygen Main_KC.py

max_stars_repo_name: lgearty/Silica-Coding-Project

max_issues_repo_name: lgearty/Silica-Coding-Project

max_forks_repo_name: lgearty/Silica-Coding-Project


text:
import os
import sys
import json
import random
import tempfile
import pkg_resources

import numpy as np
import tensorflow as tf

from aimodelshare.aws import get_s3_iam_client, run_function_on_lambda, get_aws_client

def export_reproducibility_env(seed, directory, mode="gpu"):
  # Change the output into json.dumps
  # Argument single seed for all inputs & mode
  data = {
    "global_seed_code": [
      "os.environ['PYTHONHASHSEED'] = '{}'".format(seed),
      "random.seed({})".format(seed),
      "tf.random.set_seed({})".format(seed),
      "np.random.seed({})".format(seed),
    ]
  }

  # Ignore this part for now
  # Local seed codes are tensorflow code that are
  # not affected by the global seed and sometimes require us
  # to define the seed in the function call
  data["local_seed_code"] = [
    "train_ds = tf.keras.preprocessing.image_dataset_from_directory(data_dir, validation_split=0.2, subset='training', seed={}, image_size=(img_height, img_width), batch_size=batch_size)".format(seed),
    "val_ds = tf.keras.preprocessing.image_dataset_from_directory(data_dir, validation_split=0.2, subset='validation', seed={}, image_size=(img_height, img_width), batch_size=batch_size)".format(seed),
  ]

  if mode == "gpu":
    data["gpu_cpu_parallelism_ops"] = [
      "os.environ['TF_DETERMINISTIC_OPS'] = '1'",
      "os.environ['TF_CUDNN_DETERMINISTIC'] = '1'"
    ]
  elif mode == "cpu":
    data["gpu_cpu_parallelism_ops"] = [
      "tf.config.threading.set_inter_op_parallelism_threads(1)"
    ]
  else:
    raise Exception("Error: unknown 'mode' value, expected 'gpu' or 'cpu'")

  installed_packages = pkg_resources.working_set
  installed_packages_list = sorted(["%s==%s" % (i.key, i.version)
    for i in installed_packages])

  data["session_runtime_info"] = {
    "installed_packages": installed_packages_list,
    "python_version": sys.version,
  }

  with open(os.path.join(directory, "reproducibility.json"), "w") as fp:
    json.dump(data, fp)

  return print("Your reproducibility environment is now saved to 'reproducibility.json'")

def set_reproducibility_env(reproducibility_env):
  # Change the input into dict / json
  for global_code in reproducibility_env["global_seed_code"]:
    exec("%s" % (global_code))

  for parallelism_ops in reproducibility_env["gpu_cpu_parallelism_ops"]:
    exec("%s" % (parallelism_ops))

def import_reproducibility_env(reproducibility_env_file):
  with open(reproducibility_env_file) as json_file:
    reproducibility_env = json.load(json_file)
    set_reproducibility_env(reproducibility_env)

  print("Your reproducibility environment is successfully setup")

def import_reproducibility_env_from_model(apiurl):
    if all(["AWS_ACCESS_KEY_ID" in os.environ,
            "AWS_SECRET_ACCESS_KEY" in os.environ,
            "AWS_REGION" in os.environ,
            "username" in os.environ,
            "password" in os.environ]):
        pass
    else:
        return print("'Instantiate Model' unsuccessful. Please provide credentials with set_credentials().")

    aws_client = get_aws_client()
    reproducibility_env_filename = "/runtime_reproducibility.json"

    # Get bucket and model_id for user
    response, error = run_function_on_lambda(
        apiurl, **{"delete": "FALSE", "versionupdateget": "TRUE"}
    )
    if error is not None:
        raise error

    _, bucket, model_id = json.loads(response.content.decode("utf-8"))

    try:
        resp_string = aws_client["client"].get_object(
            Bucket=bucket, Key=model_id + reproducibility_env_filename
        )

        reproducibility_env_string = resp_string['Body'].read()

    except Exception as err:
        print("This model was not deployed with reproducibility support")
        raise err

    # generate tempfile for onnx object
    temp_dir = tempfile.gettempdir()
    temp_path = os.path.join(temp_dir, 'temp_file_name')

    # save onnx to temporary path
    with open(temp_path, "wb") as f:
        f.write(reproducibility_env_string)

    import_reproducibility_env(temp_path)

ext: py

max_stars_repo_name: AIModelShare/AIModelShare

max_issues_repo_name: AIModelShare/AIModelShare

max_forks_repo_name: AIModelShare/aimodelshare

reason: import numpy


text:
import unittest

from quantlib.time.api import (
    Period, Months, Days, Years, Following, Actual365Fixed,
    UnitedStates, today, Unadjusted)
from quantlib.math.matrix import Matrix
import numpy as np
from quantlib.termstructures.volatility.swaption.swaption_vol_matrix \
    import SwaptionVolatilityMatrix
from quantlib.indexes.api import EuriborSwapIsdaFixA
from quantlib.instruments.make_cms import MakeCms
from quantlib.termstructures.yields.api import FlatForward
from quantlib.cashflows.coupon_pricer import set_coupon_pricer
from quantlib.cashflows.cap_floored_coupon import CappedFlooredCmsCoupon
from quantlib.cashflows.conundrum_pricer import (
    YieldCurveModel, NumericHaganPricer, AnalyticHaganPricer)
from quantlib.cashflows.linear_tsr_pricer import LinearTsrPricer
from quantlib.quotes import SimpleQuote
from quantlib.settings import Settings

class CmsFairRateTestCase(unittest.TestCase):
    def setUp(self):
        atm_option_tenors = [Period(1, Months), Period(6, Months)] + \
                            [Period(i, Years) for i in [1, 5, 10, 30]]
        atm_swap_tenors = [Period(1, Years), Period(5, Years), Period(10, Years),
                           Period(30, Years)]

        m = np.array([[.1300, .1560, .1390, .1220],
                      [.1440, .1580, .1460, .1260],
                      [.1600, .1590, .1470, .1290],
                      [.1640, .1470, .1370, .1220],
                      [.1400, .1300, .1250, .1100],
                      [.1130, .1090, .1070, .0930]])

        M = Matrix.from_ndarray(m)

        calendar = UnitedStates()
        self.atm_vol = SwaptionVolatilityMatrix(calendar,
                                                Following,
                                                atm_option_tenors,
                                                atm_swap_tenors,
                                                M,
                                                Actual365Fixed())

        reference_date = calendar.adjust(today())
        Settings().evaluation_date = reference_date
        self.term_structure = FlatForward(reference_date, 0.05, Actual365Fixed())
        self.swap_index = EuriborSwapIsdaFixA(Period(10, Years),
                                              forwarding=self.term_structure)

    def test_fair_rate(self):

        start_date = self.term_structure.reference_date + Period(20, Years)
        payment_date = start_date + Period(1, Years)
        end_date = payment_date
        nominal = 1.0
        gearing = 1.0
        spread = 0.0
        coupon = CappedFlooredCmsCoupon(payment_date, nominal,
                                        start_date, end_date,
                                        self.swap_index.fixing_days,
                                        self.swap_index,
                                        gearing, spread,
                                        ref_period_start=start_date,
                                        ref_period_end=end_date,
                                        day_counter=self.swap_index.day_counter)

        for model in YieldCurveModel:
            pricer = NumericHaganPricer(self.atm_vol, model, SimpleQuote(0.))
            coupon.set_pricer(pricer)
            rate_numeric = coupon.rate
            pricer = AnalyticHaganPricer(self.atm_vol, model, SimpleQuote(0.))
            coupon.set_pricer(pricer)
            rate_analytic = coupon.rate
        pricer = LinearTsrPricer(self.atm_vol, SimpleQuote(0.))
        coupon.set_pricer(pricer)
        rate_lineartsr = coupon.rate
        self.assertAlmostEqual(rate_lineartsr, rate_analytic, 3)

    def test_cms_swap(self):
        swap_lengths = [1, 5, 6, 10]
        spread = 0.
        cms_list = [MakeCms(Period(t, Years),
                            self.swap_index,
                            self.swap_index.ibor_index,
                            spread,
                            Period(10, Days))() for t in swap_lengths]
        for model in YieldCurveModel:
            num_pricer = NumericHaganPricer(self.atm_vol, model, SimpleQuote(0.))
            analytic_pricer = AnalyticHaganPricer(self.atm_vol, model, SimpleQuote(0.))
            if model == YieldCurveModel.NonParallelShifts:
                lineartsr_pricer = LinearTsrPricer(self.atm_vol, SimpleQuote(0.))
            for cms in cms_list:
                set_coupon_pricer(cms[0], num_pricer)
                price_num = cms.npv
                set_coupon_pricer(cms[0], analytic_pricer)
                price_analytic = cms.npv
                self.assertAlmostEqual(price_num, price_analytic, 3)
                if model == YieldCurveModel.NonParallelShifts:
                    set_coupon_pricer(cms[0], lineartsr_pricer)
                    self.assertAlmostEqual(cms.npv, price_analytic, 3)

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

max_stars_repo_name: bpmbank/pyql

max_issues_repo_name: bpmbank/pyql

max_forks_repo_name: bpmbank/pyql


text:
import click
import numpy as np
import pandas as pd

import NaiveDE
import SpatialDE

@click.command()
@click.argument('expression_csv', type=click.Path(exists=True), metavar='<expression csv>')
@click.argument('coordinate_csv', type=click.Path(exists=True), metavar='<cooridnates csv>')
@click.argument('results_csv', type=click.Path(), metavar='<output file>')
@click.option('--model_selection_csv', type=click.Path(), default=None)
def main(expression_csv, coordinate_csv, results_csv, model_selection_csv):
    ''' Perform SpatialDE test on data in input files.

    <expression csv> : A CSV file with expression valies. Columns are genes,
    and Rows are samples

    <coordinates csv> : A CSV file with sample coordinates. Each row is a sample,
    the columns with coordinates must be named 'x' and 'y'. For other formats
    (e.g. 1d or 3d queries), it is recommended to write a custom Python
    script to do the analysis.

    <output file> : P-vaues and other relevant values for each gene
    will be stored in this file, in CSV format.

    '''
    df = pd.read_csv(expression_csv, index_col=0)

    df = df.T[df.sum(0) >= 3].T  # Filter practically unobserved genes

    sample_info = pd.read_csv(coordinate_csv, index_col=0)

    sample_info['total_counts'] = df.sum(1)
    sample_info = sample_info.query('total_counts > 5')  # Remove empty features

    df = df.loc[sample_info.index]
    X = sample_info[['x', 'y']]

    # Convert data to log-scale, and account for depth
    dfm = NaiveDE.stabilize(df.T).T
    res = NaiveDE.regress_out(sample_info, dfm.T, 'np.log(total_counts)').T

    # Perform Spatial DE test with default settings
    results = SpatialDE.run(X, res)

    # Save results and annotation in files for interactive plotting and interpretation
    results.to_csv(results_csv)

    if not model_selection_csv:
        return results

    de_results = results[(results.qval < 0.05)].copy()
    ms_results = SpatialDE.model_search(X, res, de_results)

    ms_results.to_csv(model_selection_csv)

    return results, ms_results

max_stars_repo_name: jiangfuqing/SpatialDE

max_issues_repo_name: jiangfuqing/SpatialDE

max_forks_repo_name: jiangfuqing/SpatialDE


text:
# coding: utf-8
"""
__project_ = 'blog'
__file_name__ = 'tasks'
__author__ = 'admin'
__time__ = '2020-04-13 11:20'
__product_name = PyCharm
"""
import platform

from celery import shared_task
import re
import json
# 词频统计库
import collections
import numpy as np
import jieba
import wordcloud
from PIL import Image

from blog.settings import STATICFILES_DIRS, WINDOWS_FONT_PATH, UBUNTU_FONT_PATH
from storm.models import Article

@shared_task
def word_cloud_handle():
    """
    生成词云
    :param :
    :return:
    """
    print("---开始生成词云---")
    article_content = Article.objects.all().values_list('body', flat=True)
    # 转化为字符串
    article_content_str = json.dumps(list(article_content), ensure_ascii=False)
    # 文本预处理
    # pattern_char = re.compile(u'\t|\n|\.|-|:|;|\)|\(|\?|"')  # 定义正则表达式匹配模式

    # 正则匹配中文
    pattern = re.compile(r'[\u4e00-\u9fa5]+')
    string_data = pattern.findall(article_content_str)
    string_data = json.dumps(string_data, ensure_ascii=False)
    # string_data = re.sub(pattern, '', string_data)  # 将符合模式的字符去除
    # 文本分词
    seg_list_exact = jieba.cut(string_data, cut_all=False)  # 精确模式分词
    object_list = []
    remove_words = [u'的', u'，', u'和', u'是', u'随着', u'对于', u'对', u'等', u'能', u'都', u'。', u' ', u'、', u'中', u'在', u'了',
                    u'通常', u'如果', u'我们', u'需要', u'你', u'"', u'[', u']', u',']  # 自定义去除词库
    # 循环读出每个分词，如果不在去除词库中，分词追加到列表
    for word in seg_list_exact:
        if word not in remove_words:
            object_list.append(word)

    # 词频统计
    word_counts = collections.Counter(object_list)
    # word_counts_top10 = word_counts.most_common(100)  # 获取前10最高频的词
    # print(word_counts_top10)  # 输出检查

    # 判断系统环境，字体文件
    platform_ = platform.system()
    if platform_ == "Windows":
        font_path = WINDOWS_FONT_PATH
    elif platform_ == "Linux":
        font_path = UBUNTU_FONT_PATH
    elif platform_ == "Mac":
        font_path = ''

    # 词频展示
    mask = np.array(Image.open(STATICFILES_DIRS[0] + '/wordcloud.jpg'))  # 定义词频背景
    wc = wordcloud.WordCloud(
        font_path=font_path,  # 设置字体格式
        mask=mask,  # 设置背景图
        max_words=1000,  # 最多显示词数
        max_font_size=40,  # 字体最大值
        height=500,
    )
    # 从字典生成词云
    wc.generate_from_frequencies(word_counts)
    # 从背景图建立颜色方案
    image_colors = wordcloud.ImageColorGenerator(mask)
    # 将词云颜色设置为背景图方案
    wc.recolor(color_func=image_colors)
    # 保存图片
    wc.to_file(STATICFILES_DIRS[0] + "/pywordcloud.png")
    print("---生成词云完成---")

max_stars_repo_name: Iamacode/blog

max_issues_repo_name: Iamacode/blog

max_forks_repo_name: Iamacode/blog


text:
import tensorflow as tf
from tensorflow.python.framework import ops
import numpy as np

def bf16cut_np(xmlp):
  fp32a = np.asarray(xmlp,dtype=np.float32)
  bits = fp32a.view(np.int32)
  # for bf16 8 bit mantis
  b = np.left_shift(np.right_shift(bits,15),15)
  # for bf12 4 bit mantis
  #b = np.left_shift(np.right_shift(bits,19),19)
  # for bf10 2 bit mantis
  #b = np.left_shift(np.right_shift(bits,21),21)
  # for bf9 1 bit mantis
  #b = np.left_shift(np.right_shift(bits,22),22)
  # for bf8 0 bit mantis
  #b = np.left_shift(np.right_shift(bits,23),23)
  return b.view(np.float32)

def bf16cut_tf(xmlp):
  aabs=tf.abs(xmlp)
  rnd2exp=tf.pow(2.0,tf.floor(tf.where(tf.equal(aabs,0.0),aabs, tf.divide(tf.log(aabs),tf.log(2.0))))) # get back the old value rounded to 2 exp
  # for bf16 8 bit mantis
  return  tf.multiply(tf.multiply(tf.divide(tf.floor(tf.multiply(tf.divide(aabs,rnd2exp),256.0)),256.0),rnd2exp),tf.sign(xmlp))
  # for bf12 4 bit mantis
  #return  tf.multiply(tf.multiply(tf.divide(tf.floor(tf.multiply(tf.divide(aabs,rnd2exp),16.0)),16.0),rnd2exp),tf.sign(xmlp))
  # for bf10 2 bit mantis
  #return  tf.multiply(tf.multiply(tf.divide(tf.floor(tf.multiply(tf.divide(aabs,rnd2exp),4.0)),4.0),rnd2exp),tf.sign(xmlp))
  # for bf9 1 bit mantis
  #return  tf.multiply(tf.multiply(tf.divide(tf.floor(tf.multiply(tf.divide(aabs,rnd2exp),2.0)),2.0),rnd2exp),tf.sign(xmlp))
  # for bf8 0 bit mantis
  #return  tf.multiply(tf.multiply(tf.divide(tf.floor(tf.multiply(tf.divide(aabs,rnd2exp),1.0)),1.0),rnd2exp),tf.sign(xmlp))

#def bf16cut_np_merge(xmlp):
#  aabs=np.abs(xmlp)
#  rnd2exp=np.power(2.0,np.floor(np.where(aabs==0.0,0.0,np.log2(aabs)))) # get back the old value rounded to 2 exp
#  return np.multiply(np.multiply(np.divide(np.floor(np.multiply(np.divide(aabs,rnd2exp),256.0)),256.0),rnd2exp),np.sign(xmlp))

def id_ssy(xmlp):
  return xmlp

def id_ssy_grad(op,grad):
    x = op.inputs[0]
    return grad

# Actual gradient:
def id_bf16cut_grad(op, grad):
    x = op.inputs[0]
    # must use tf version or else rise no xla compiler error
    return bf16cut_tf(grad )  # add a "small" error just to see the difference:
    #return grad *  x  # add a "small" error just to see the difference:

# Define custom py_func which takes also a grad op as argument:
def py_func(func, inp, Tout, stateful=True, name=None, grad=None):

    # Need to generate a unique name to avoid duplicates:
    rnd_name = 'PyFuncGrad' + str(np.random.randint(0, 1E+8))

    tf.RegisterGradient(rnd_name)(grad)  # see id_bf16cut_grad for grad example
    g = tf.get_default_graph()
    with g.gradient_override_map({"PyFunc": rnd_name}):
        return tf.py_func(func, inp, Tout, stateful=stateful, name=name)

def id_bf16cut_fp(x, name=None):

    with ops.op_scope([x], name, "id_bf16cut_fp") as name:
        #sqr_x = py_func(np.square,
        sqr_x = py_func(bf16cut_np, # must use np version orelse will rise unsupport  data type, change to id_ssy will be as fast as no custom oprator case, so np instead of H2D/D2H is the reason that slow down custom operator
                        [x],
                        [tf.float32],
                        name=name,
                        grad=id_ssy_grad)  # <-- here's the call to the gradient
        return sqr_x[0]

def id_bf16cut_bp(x, name=None):

    with ops.op_scope([x], name, "id_bf16cut_bp") as name:
        #sqr_x = py_func(np.square,
        sqr_x = py_func(id_ssy, # must use np version orelse will rise unsupport  data type
                        [x],
                        [tf.float32],
                        name=name,
                        grad=id_bf16cut_grad)  # <-- here's the call to the gradient
        return sqr_x[0]

if __name__ == "__main__":
  with tf.Session() as sess:
    x = tf.constant([0,0,1.123456789, 11.123456789,111.123456789,-1.123456789,-11.123456789,-111.123456789],shape=[2,4])
    y = id_bf16cut(x)
    tf.initialize_all_variables().run()

    print(x.eval())
    print(y.eval())
    print(tf.gradients(y, x)[0].eval())

max_stars_repo_name: shengyushen/training

max_issues_repo_name: shengyushen/training

max_forks_repo_name: shengyushen/training


text:
import torch
import numpy as np
from torch.utils.data import Dataset, DataLoader
from pycocotools.coco import COCO
import cv2

class dataset(Dataset):
    def __init__(self, image_folder, annotation_file):
        self.images = image_folder
        self.annotations = annotation_file

        self.coco = COCO(self.annotations)
        self.catId = self.coco.getCatIds(catNms = ['person'])

        self.imgIds = self.coco.getImgIds(catIds = self.catId)

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

    def __getitem__(self, id):
        # print('yahan baar baar nahi aana chahiye')
        print(self.imgIds[id])
        img = self.coco.loadImgs(self.imgIds[id])[0]
        img_name = img['file_name']
        # print(img_name)
        X = np.array(cv2.imread(self.images + str(img_name)))
        height = X.shape[0]
        width = X.shape[1]
        # print(height, width)
        X = cv2.resize(X,(224, 224))
        # print(X.shape)
        X = X.transpose([2,0,1])
        annIds = self.coco.getAnnIds(imgIds=img['id'], catIds=self.catId, iscrowd=None)
        anns = self.coco.loadAnns(annIds)
        y = anns[0]['keypoints']
        y = np.array(y)
        # print(y)
        for i in range(17):
            y[3*i] = 224*1.0/width * y[3*i]
            y[3*i + 1] = 224*1.0/height * y[3*i + 1]
        # print(y)
        return X, y
data_train = dataset('coco/images/', 'coco/annotations/person_keypoints_train2014.json')
train_loader =  DataLoader(data_train, batch_size = 8, shuffle = True)
print(len(train_loader))
data_train.__getitem__(5)


text:
""" Interpolation-related functions and classes.
"""

# p2.6+ compatibility
from __future__ import division, print_function, unicode_literals

import numpy as np
from .utilities import between, indices_from_grid, meshgrid_nd

class CubicSpline(object):
    """ Class to generate a cubic spline through a set of points.

    After initialisation, an instance can be called with an array of
    values xp, and will return the cubic-spline interpolated values yp
    = f(xp).

    The spline can be reset to use a new first and last derivative
    while still using the same initial points by calling the _set_d2()
    method.

    If you want to calculate a new spline using a different set of x
    and y values, you'll have to instantiate a new class.

    The spline generation is based on the NR algorithms (but not their
    code.)
    """
    def __init__(self, x, y, firstderiv=None, lastderiv=None, nochecks=False):
        """
        Parameters
        ----------
        x, y : arrays of floats, shape (N,)
          x and y = f(x).
        firstderiv : float (None)
          Derivative of f(x) at x[0]
        lastderiv : float (None)
          Derivative of f(x) at x[-1]
        nochecks : bool (False)
          If False, check the x array is sorted and unique. Set to True
          for increased speed.
        """
        x = np.asarray(x, dtype=float)
        y = np.asarray(y, dtype=float)
        if not nochecks:
            # check all x values are unique
            if len(x) - len(np.unique(x)) > 0:
                raise Exception('non-unique x values were found!')
            cond = x.argsort()                    # sort arrays
            x = x[cond]
            y = y[cond]

        self.x = x
        self.y = y
        self.npts = len(x)
        self._set_d2(firstderiv, lastderiv)

    def __call__(self,xp):
        """ Given an array of x values, returns cubic-spline
        interpolated values yp = f(xp) using the derivatives
        calculated in _set_d2().
        """
        x = self.x;  y = self.y;  npts = self.npts;  d2 = self.d2

        # make xp into an array
        if not hasattr(xp,'__len__'):  xp = (xp,)
        xp = np.asarray(xp)

        # for each xp value, find the closest x value above and below
        i2 = np.searchsorted(x,xp)

        # account for xp values outside x range
        i2 = np.where(i2 == npts, npts-1, i2)
        i2 = np.where(i2 == 0, 1, i2)
        i1 = i2 - 1

        h = x[i2] - x[i1]
        a = (x[i2] - xp) / h
        b = (xp - x[i1]) / h
        temp = (a**3 - a)*d2[i1] +  (b**3 - b)*d2[i2]
        yp = a * y[i1] + b * y[i2] + temp * h**2 / 6.

        return yp

    def _tridiag(self,temp,d2):
        x, y, npts = self.x, self.y, self.npts
        for i in range(1, npts-1):
            ratio = (x[i]-x[i-1]) / (x[i+1]-x[i-1])
            denom = ratio * d2[i-1] + 2.       # 2 if x vals equally spaced
            d2[i] = (ratio - 1.) / denom       # -0.5 if x vals equally spaced
            temp[i] = (y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1])
            temp[i] = (6.*temp[i]/(x[i+1]-x[i-1]) - ratio * temp[i-1]) / denom
        return temp

    def _set_d2(self, firstderiv=None, lastderiv=None, verbose=False):
        """ Calculates the second derivative of a cubic spline
        function y = f(x) for each value in array x. This is called by
        __init__() when a new class instance is created.

        Parameters
        ----------
        firstderiv : float, (None)
          1st derivative of f(x) at x[0].  If None, then 2nd
          derivative is set to 0 ('natural').
        lastderiv : float (None)
          1st derivative of f(x) at x[-1].  If None, then 2nd
          derivative is set to 0 ('natural').
        """
        if verbose:  print('first deriv,last deriv', firstderiv, lastderiv)
        x, y, npts = self.x, self.y, self.npts
        d2 = np.empty(npts)
        temp = np.empty(npts-1)

        if firstderiv is None:
            if verbose:  print("Lower boundary condition set to 'natural'")
            d2[0] = 0.
            temp[0] = 0.
        else:
            d2[0] = -0.5
            temp[0] = 3./(x[1]-x[0]) * ((y[1]-y[0])/(x[1]-x[0]) - firstderiv)

        temp = self._tridiag(temp,d2)

        if lastderiv is None:
            if verbose:  print("Upper boundary condition set to 'natural'")
            qn = 0.
            un = 0.
        else:
            qn = 0.5
            un = 3./(x[-1]-x[-2]) * (lastderiv - (y[-1]-y[-2])/(x[-1]-x[-2]))

        d2[-1] = (un - qn*temp[-1]) / (qn*d2[-2] + 1.)
        for i in reversed(range(npts-1)):
            d2[i] = d2[i] * d2[i+1] + temp[i]

        self.d2 = d2

class AkimaSpline(object):
    """ Descibes an Akima Spline through a set of points.

    It must be instantiated with a set of `xvals` and `yvals` knot values,
    and then can be called with a new set of x values `x`. This is
    used by `interp_Akima`, see its documentation for more
    information.

    References
    ----------
    "A new method of interpolation and smooth curve fitting based
    on local procedures." Hiroshi Akima, J. ACM, October 1970, 17(4),
    589-602.

    Notes
    -----
    This is adapted from a function written by `Christoph Gohlke
    <http://www.lfd.uci.edu/~gohlke/>`_ under a BSD license:

    Copyright (c) 2007-2012, Christoph Gohlke
    Copyright (c) 2007-2012, The Regents of the University of California
    Produced at the Laboratory for Fluorescence Dynamics
    All rights reserved.
    """
    def __init__(self, xvals, yvals):
        """
        Parameters
        ----------
        xvals, yvals : array_like, shape (N,)
          Reference values. xvals cannot contain duplicates.
        """

        x = np.asarray(xvals, dtype=np.float64)
        y = np.asarray(yvals, dtype=np.float64)
        if x.ndim != 1:
            raise ValueError("x array must be one dimensional")

        n = len(x)
        if n < 3:
            raise ValueError("Array too small")
        if n != len(y):
            raise ValueError("Size of x-array must match data shape")

        dx = np.diff(x)
        if (dx <= 0.0).any():
            isort = np.argsort(x)
            x = x[isort]
            y = y[isort]
            dx = np.diff(x)
            if (dx == 0.).any():
                raise ValueError("x array has duplicate values")

        m = np.diff(y) / dx
        mm = 2. * m[0] - m[1]
        mmm = 2. * mm - m[0]
        mp = 2. * m[n - 2] - m[n - 3]
        mpp = 2. * mp - m[n - 2]

        m1 = np.concatenate(([mmm], [mm], m, [mp], [mpp]))

        dm = np.abs(np.diff(m1))
        f1 = dm[2:n + 2]
        f2 = dm[0:n]
        f12 = f1 + f2

        ids = np.nonzero(f12 > 1e-9 * f12.max())[0]
        b = m1[1:n + 1]

        b[ids] = (f1[ids] * m1[ids + 1] + f2[ids] * m1[ids + 2]) / f12[ids]
        c = (3. * m - 2. * b[0:n - 1] - b[1:n]) / dx
        d = (b[0:n - 1] + b[1:n] - 2. * m) / dx ** 2

        self.xvals, self.yvals, self.b, self.c, self.d = x, y, b, c, d

    def __call__(self, x):
        """
        Parameters
        ----------
        x : array_like, shape (M,)
          Values at which to interpolate.

        Returns
        -------
        vals : ndarray, shape (M,)
           Interpolated values.
        """
        x = np.asarray(x, dtype=np.float64)

        if x.ndim != 1:
            raise ValueError("Array must be one dimensional")

        c0 = x < self.xvals[0]
        c2 = x > self.xvals[-1]
        c1 = ~(c0 | c2)
        if c1.sum() == 0:
            raise ValueError('All points are outside the interpolation range!')
        x1 = x[c1]
        out = np.empty_like(x)
        bins = np.digitize(x1, self.xvals)
        bins = np.minimum(bins, len(self.xvals) - 1) - 1
        b = bins[0:len(x1)]
        wj = x1 - self.xvals[b]
        out[c1] = ((wj * self.d[b] + self.c[b])*wj + self.b[b])*wj + \
                  self.yvals[b]

        # use linear extrapolation for points outside self.xvals
        if c0.any():
            y = out[c1]
            slope = (y[1] - y[0]) / (x1[1] - x1[0])
            intercept = y[0] - slope * x1[0]
            out[c0] = x[c0] *slope + intercept
        if c2.any():
            y = out[c1]
            slope = (y[-2] - y[-1]) / (x1[-2] - x1[-1])
            intercept = y[-1] - slope * x1[-1]
            out[c2] = x[c2] *slope + intercept

        return out

def fit_spline(x, y, bins=4, addknots=None, estimator=np.median):
    """ Find a smooth function that approximates `x`, `y`.

    `bins` is the number of bins into which the sample is split. Returns
    a function f(x) that approximates y from min(x) to max(x).

    Parameters
    ----------
    addknots : sequence of float pairs
       A sequence of (x,y) pairs values that the spline is forced to
       pass through. Default is None.

    Notes
    -----
    The sample is split into bins number of sub-samples with evenly
    spaced x values. The median x and y value within each subsample is
    measured, and a cubic spline is drawn through these subsample
    median points.

    `x` must be sorted lowest -> highest, but need not be evenly spaced.
    """
    x,y = map(np.asarray, (x,y))
    good = ~np.isnan(x)
    good &= ~np.isnan(y)
    x = x[good]
    y = y[good]
    binedges = np.linspace(x.min(), x.max(), bins+1)
    medvals = []
    cbins = []
    for x0,x1 in zip(binedges[:-1], binedges[1:]):
        cond = between(x, x0, x1)
        if not cond.any():
            continue
        cbins.append(estimator(x[cond]))
        medvals.append(estimator(y[cond]))

    if addknots is not None:
        for x,y in addknots:
            cbins.append(x)
            medvals.append(y)
        isort = np.argsort(cbins)
        cbins = np.array(cbins)[isort]
        medvals = np.array(medvals)[isort]

    if len(cbins) < 3:
        raise RuntimeError('Too few bins')
    return AkimaSpline(cbins, medvals)

def interp_Akima(x_new, x, y):
    """Return interpolated data using Akima's method.

    Akima's interpolation method uses a continuously differentiable
    sub-spline built from piecewise cubic polynomials. The resultant
    curve passes through the given data points and will appear smooth
    and natural.

    Parameters
    ----------
    x_new : array_like, shape (M,)
        Values at which to interpolate.
    x, y : array_like, shape (N,)
        Reference values. x cannot contain duplicates.

    Returns
    -------
    vals : ndarray, shape (M,)
       Interpolated values.

    References
    ----------
    "A new method of interpolation and smooth curve fitting based
    on local procedures." Hiroshi Akima, J. ACM, October 1970, 17(4),
    589-602.

    Examples
    --------
    Plot interpolated Gaussian noise:

    >>> x = np.sort(np.random.random(10) * 100)
    >>> y = np.random.normal(0.0, 0.1, size=len(x))
    >>> x2 = np.arange(x[0], x[-1], 0.02)
    >>> y2 = interp_Akima(x2, x, y)
    >>> y3 = interp_spline(x2, x, y)
    >>> from matplotlib import pyplot as plt
    >>> plt.plot(x2, y2, 'b-', label='Akima')
    >>> plt.plot(x2, y3, 'r-', label='Cubic')
    >>> plt.plot(x, y, 'co')
    >>> plt.legend()
    >>> plt.show()
    """
    interpolator = AkimaSpline(x, y)
    return interpolator(x_new)

def interp_spline(x, xvals, yvals, nochecks=False):
    """ Like `numpy.interp`, but using spline instead of linear
    interpolation.

    This is a convenience function wrapped around CubicSpline.
    """
    spl = CubicSpline(xvals, yvals, nochecks=nochecks)
    return spl(x)

def splice(co0, co1, i, j, forced=None):
    """ Join two overlapping curves smoothly using a cubic spline.

    Parameters
    ----------
    co0, co1 : arrays of shape (N,)
      The two curves to be joined. They must have the same length and
      overlap completely.
    i, j : int
      Roughly speaking, co0 values will be retained below i, and co1 values
      will be retained above j.
    forced : int, optional
      The number of pixels and continuum values between i and j that
      continuum will be forced to pass through.

    Returns the new continuum array of shape (N,).
    """
    # go npix to either side of the joining point, and measure slopes
    newco = np.empty_like(co0)

    # derivatives
    d1 = co0[i] - co0[i-1]
    d2 = co1[j] - co1[j-1]

    if forced is not None:
        indices = [i] + list(zip(*forced)[0]) + [j]
        covals = [co0[i]] + list(zip(*forced)[1])+ [co1[j]]
    else:
        indices = [i,j]
        covals = [co0[i],co1[j]]

    indices = np.array(indices, dtype=float)
    spline = CubicSpline(indices, covals,firstderiv=d1, lastderiv=d2)
    newco[:i] = co0[:i].copy()
    newco[i:j] = spline(range(i,j))
    newco[j:] = co1[j:].copy()

    return newco

def _interp3d(ix, iy, iz, a):
    """ Trilinear interpolation.

    Parameters
    ----------
    ix, iy, iz : arrays of floats, shape (M, )
      Coordinates of the points at which to interpolate. Values should
      be floats from 0 - I-1, 0 - J-1 and 0 - K-1.
    a: array of floats, shape (I, J, K)

    Returns
    -------
    output : array of floats, shape (M,)

    Notes
    -----
    This function isn't intended to be used directly. The public
    trilinear interpolation function is `trilinear_interp`.
    """
    assert len(ix) == len(iy) == len(iz)

    output = np.empty(len(ix))

    x0 = ix.astype(int)
    y0 = iy.astype(int)
    z0 = iz.astype(int)
    x1 = x0 + 1
    y1 = y0 + 1
    z1 = z0 + 1

    x1[x1 == a.shape[0]] = a.shape[0] - 1
    y1[y1 == a.shape[1]] = a.shape[1] - 1
    z1[z1 == a.shape[2]] = a.shape[2] - 1

    x = ix - x0
    y = iy - y0
    z = iz - z0
    output = (a[x0,y0,z0] * (1-x) * (1-y) * (1-z) +
              a[x1,y0,z0] * x * (1-y) * (1-z) +
              a[x0,y1,z0] * (1-x) * y * (1-z) +
              a[x0,y0,z1] * (1-x) * (1-y) * z +
              a[x1,y0,z1] * x * (1-y) * z +
              a[x0,y1,z1] * (1-x) * y *z +
              a[x1,y1,z0] * x * y * (1-z) +
              a[x1,y1,z1] * x * y * z)

    return output

def trilinear_interp(x, y, z, xref, yref, zref, vals):
    """ Trilinear interpolation.

    Parameters
    ----------
    x, y, z : arrays of floats, shapes (M,), (N,), (O,)
      Coordinate grid at which to interpolate `vals`.
    xref, yref, zref : array of floats, shapes (I,), (J,), (K,)
      Reference coordinate grid. The grid must be equally spaced along
      each direction, but the spacing can be different between
      directions.
    vals : array of floats, shape (I, J, K)
      Reference values at the reference grid positions.

    Returns
    -------
    output : array of floats, shape (M, N, O)

    See Also
    --------
    barak.plot.arrplot for plotting 2-d slices of the reference and
    interpolated arrays.
    """
    assert (len(xref), len(yref), len(zref)) == vals.shape

    ix = indices_from_grid(x, xref)
    iy = indices_from_grid(y, yref)
    iz = indices_from_grid(z, zref)
    #iX, iY, iZ = (a.ravel() for a in meshgrid_nd(ix, iy, iz))
    iX, iY, iZ = (a.ravel() for a in np.meshgrid(ix, iy, iz))
    out = _interp3d(iX, iY, iZ, vals)

    # Note the index order and transpose!
    return out.reshape(len(z), len(y), len(x)).T

def CloughTocher2d_interpolator(xref, yref, vals):
    """ 2d interpolation (requires Scipy).

    This is a convenience wrapper around
    scipy.interpolation.CloughTocher2Dinterpolator, that saves you
    having to worry about using meshgrid.

    Parameters
    ----------
    xref, yref : array of floats, shapes (J,), (I,)
      Reference coordinate grid. The grid must be equally spaced along
      each direction, but the spacing can be different between
      directions.
    vals : array of floats, shape (I, J)
      Reference values at the reference grid positions.

    Returns
    -------
    interpolator: CloughTocher2DInterpolater instance
      Object that accepts a (y,x) tuple (note reversed order from the
      input to this function!) and returns the interpolated value.

    See Also
    --------
    barak.plot.arrplot for plotting the reference and interpolated arrays.
    """
    try:
        from scipy.interpolate import CloughTocher2DInterpolator
    except ImportError:
        print("Scipy (Version 0.9 or greater) must be installed to "
              "use CloughTocher2d_interp")
        return None

    assert (len(yref), len(xref)) == vals.shape
    XREF,YREF = np.meshgrid(xref, yref)
    interpolator = CloughTocher2DInterpolator((YREF.ravel(), XREF.ravel()),
                                              vals.ravel())

    return interpolator

def CloughTocher2d_interp(x, y, xref, yref, vals):
    """ 2d interpolation (requires Scipy). This is a convenience
    wrapper around CloughTocher2d_interpolator and has a similar
    interface to numpy.interp.

    Parameters
    ----------
    x, y : arrays of floats, shapes (N,), (M,)
      Coordinate grid at which to interpolate `vals`.
    xref, yref : array of floats, shapes (J,), (I,)
      Reference coordinate grid. The grid must be equally spaced along
      each direction, but the spacing can be different between
      directions.
    vals : array of floats, shape (I, J)
      Reference values at the reference grid positions.

    Returns
    -------
    output : array of floats, shape (M, N)

    See Also
    --------
    barak.plot.arrplot for plotting the reference and interpolated arrays.
    """
    interp = CloughTocher2d_interpolator(xref, yref, vals)
    X, Y = np.meshgrid(x,y)
    return interp((Y, X))

class MapCoord_Interpolator:
    """ A wrapper around scipy.indimage.map_coordinates

    It creates an object that can be repeatedly called to perform
    interpolation at different coordinate values.

    Tested in 2, 3 and 4 dimensions.
    """
    def __init__(self, data, vals):
        """
        data : shape (d0, d1, d2, ...)

        vals = [..., x2, x1, x0]

        where

        x0 : shape (d0,)
        x1 : shape (d1,)
        x2 : shape (d2,)

        ...

        """
        assert data.shape == tuple(len(v) for v in vals)
        self.val0 = [v[0] for v in vals]
        self.dval = [v[1] - v[0] for v in vals]
        self.data = data

    def map_xvals_to_coord(self, vals):
        """ Find the coordinates to give to map_coordinates from each
        axis value.
        """
        # Check we have as many coordinates as dimensions
        assert len(vals) == len(self.data.shape)
        vals = map(np.atleast_1d, vals)
        coords = []
        for i,v in enumerate(vals):
            coords.append( (v - self.val0[i]) / self.dval[i])
        return coords

    def __call__(self, vals, order=3, **kwargs):
        """
        note: order of vals is x2,x1,x0

        order : int  (default 3)
          Spline order.

        Possible keywords:
           for extrapolation outside the initial grid

           mode='nearest'
        """
        from scipy.ndimage import map_coordinates
        coords = self.map_xvals_to_coord(vals)
        out = map_coordinates(self.data, coords, order=order, **kwargs)
        return out

max_stars_repo_name: ntejos/Barak

max_issues_repo_name: ntejos/Barak

max_forks_repo_name: ntejos/Barak

reason: import numpy,from scipy


text:
'''
Simple OpenOpt graph stability number example;
requires networkx (http://networkx.lanl.gov)
and FuncDesigner installed.
For maximum clique problem you could use STAB on complementary graph, for networkx it's nx.complement(G)
Unlike networkx maximum_independent_set() we search for *exact* solution.
Limitations on time, cputime, "enough" value, basic GUI features are available.
'''
from openopt import STAB

import networkx as nx
G = nx.path_graph(15) # [(0,1), (1,2), (2,3), ..., (13,14)]
# you can use string and other types, e.g.
# G = nx.Graph(); G.add_node('asdf'); G.add_nodes_from(['a1', 'b2','c3']); G.add_edge('qwerty1','qwerty2')
p = STAB(G, includedNodes = [1, 5, 8], excludedNodes = [0, 4, 10])
# or
# p = STAB(G, includedNodes = ['a1', 'c3'], excludedNodes = ['qwerty1'])

# includedNodes and excludedNodes are optional arguments - nodes that must be present or absent in the solution

# solve by BSD- licensed global nonlinear solver interalg (http://openopt.org/interalg):
r = p.solve('interalg', iprint = 0) # set it to -1 to completely suppress interalg output; may not work for glpk, lpSolve
# or solve by MILP solver:
#r = p.solve('lpSolve') # see http://openopt.org/MILP for other MILP solvers available in OpenOpt

# if your solver is cplex or interalg, you can provide some stop criterion,
# e.g. maxTime, maxCPUTime, fEnough etc, for example
# r = p.solve('cplex', maxTime = 100, fEnough = 10)
# i.e. stop if solution with 10 nodes has been obtained or 100 sec were elapsed
# (gplk and lpSolve has no iterfcn connected and thus cannot handle those criterions)

# also you can use p.manage() to enable basic GUI (http://openopt.org/OOFrameworkDoc#Solving)
# it requires tkinter installed, that is included into PythonXY, EPD;
# for linux use [sudo] easy_install tk or [sodo] apt-get install python-tk
#r = p.manage('cplex')

print(r.ff, r.solution)
'''
Solver:   Time Elapsed = 0.34   CPU Time Elapsed = 0.28
(6, [1, 3, 5, 8, 12, 14])
for the example with text nodes it will be
['qwerty1', 'a1', 'b2', 'c3', 'asdf']
'''

max_stars_repo_name: wangyum/anaconda

max_issues_repo_name: wangyum/anaconda

max_forks_repo_name: wangyum/anaconda

reason: import networkx


text:
#
# This file is part of stcs-mimpc.
#
# Copyright (c) 2020 Adrian Bürger, Angelika Altmann-Dieses, Moritz Diehl.
# Developed at HS Karlsruhe and IMTEK, University of Freiburg.
# All rights reserved.
#
# The BSD 3-Clause License
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# 1. Redistributions of source code must retain the above copyright notice, this
#    list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright notice,
#    this list of conditions and the following disclaimer in the documentation
#    and/or other materials provided with the distribution.
#
# 3. Neither the name of the copyright holder nor the names of its
#    contributors may be used to endorse or promote products derived from
#    this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#

import os
import copy
import datetime as dt
import numpy as np
import casadi as ca

from system import System
from nlpsetup import NLPSetupMPC, NLPSetupMHE

from abc import ABCMeta, abstractmethod

import logging

logger = logging.getLogger(__name__)

class NLPSolverMPCBaseClass(NLPSetupMPC, metaclass = ABCMeta):

    _SECONDS_TIMEOUT_BEFORE_NEXT_TIME_GRID_POINT = 5.0
    _LOGFILE_LOCATION = "/tmp"

    @property
    def time_grid(self):

        return self._timing.time_grid

    @property
    def x_data(self):

        try:

            return np.asarray(self._x_data)

        except AttributeError:

            msg = "Optimized states not available yet, call solve() first."

            logging.error(msg)
            raise RuntimeError(msg)

    @property
    def x_hat(self):

        return np.asarray(self._predictor.x_hat)

    @property
    def u_data(self):

        try:

            return np.asarray(self._u_data)

        except AttributeError:

            msg = "Optimized continuous controls not available yet, call solve() first."

            logging.error(msg)
            raise RuntimeError(msg)

    @property
    def b_data(self):

        try:

            return np.asarray(self._b_data)

        except AttributeError:

            msg = "Optimized binary controls not available yet, call solve() first."

            logging.error(msg)
            raise RuntimeError(msg)

    @property
    def c_data(self):

        return np.asarray(self._ambient.c_data)

    @property
    def r_data(self):

        try:

            return np.asarray(self._r_data)

        except AttributeError:

            msg = "Optimized residuals not available yet, call solve() first."

            logging.error(msg)
            raise RuntimeError(msg)

    @property
    def s_ac_lb_data(self):

        try:

            return np.asarray(self._s_ac_lb_data)

        except AttributeError:

            msg = "Optimized slacks for minimum AC operation temperatures not available yet, call solve() first."

            logging.error(msg)
            raise RuntimeError(msg)

    @property
    def s_ac_ub_data(self):

        try:

            return np.asarray(self._s_ac_ub_data)

        except AttributeError:

            msg = "Optimized slacks for maximum AC operation temperatures not available yet, call solve() first."

            logging.error(msg)
            raise RuntimeError(msg)

    @property
    def s_x_data(self):

        try:

            return np.asarray(self._s_x_data)

        except AttributeError:

            msg = "Optimized slacks for soft state constraints not available yet, call solve() first."

            logging.error(msg)
            raise RuntimeError(msg)

    @property
    def s_ppsc_fpsc_data(self):

        try:

            return np.asarray(self._s_ppsc_fpsc_data)

        except AttributeError:

            msg = "Optimized slacks for FPSC safety pump speed not available yet, call solve() first."

            logging.error(msg)
            raise RuntimeError(msg)

    @property
    def s_ppsc_vtsc_data(self):

        try:

            return np.asarray(self._s_ppsc_vtsc_data)

        except AttributeError:

            msg = "Optimized slacks for VTSC safety pump speed not available yet, call solve() first."

            logging.error(msg)
            raise RuntimeError(msg)

    @property
    def solver_name(self):

        return self._solver_name

    def _setup_timing(self, timing):

        self._timing = timing

    def _setup_ambient(self, ambient):

        self._ambient = ambient

    def _set_previous_solver(self, previous_solver):

        self._previous_solver = previous_solver

    def _set_predictor(self, predictor):

        self._predictor = predictor

    def _setup_solver_name(self, solver_name):

        self._solver_name = solver_name

    def _setup_general_solver_options(self):

        self._nlpsolver_options = {}

        self._nlpsolver_options["ipopt.linear_solver"] = "mumps"

        self._nlpsolver_options["ipopt.mumps_mem_percent"] = 10000
        self._nlpsolver_options["ipopt.mumps_pivtol"] = 0.001

        self._nlpsolver_options["ipopt.print_level"] = 5

        self._nlpsolver_options["ipopt.output_file"] = os.path.join( \
            self._LOGFILE_LOCATION, self._solver_name + ".log")
        self._nlpsolver_options["ipopt.file_print_level"] = 5

        self._nlpsolver_options["ipopt.max_cpu_time"] = 720.0

    @abstractmethod
    def _setup_additional_nlpsolver_options(self):

        pass

    def __init__(self, timing, ambient, previous_solver, predictor, solver_name):

        logger.debug("Initializing NLP solver " + solver_name + " ...")

        super().__init__(timing)

        self._setup_timing(timing = timing)
        self._setup_ambient(ambient = ambient)

        self._setup_solver_name(solver_name = solver_name)
        self._set_previous_solver(previous_solver = previous_solver)
        self._set_predictor(predictor = predictor)

        self._setup_general_solver_options()
        self._setup_additional_nlpsolver_options()
        self._setup_collocation_options()

        logger.debug("NLP solver "  + solver_name + " initialized.")

    def set_solver_max_cpu_time(self, time_point_to_finish):

        max_cpu_time = (time_point_to_finish - dt.datetime.now(tz = self._timing.timezone) - \
            dt.timedelta(seconds = self._SECONDS_TIMEOUT_BEFORE_NEXT_TIME_GRID_POINT)).total_seconds()

        self._nlpsolver_options["ipopt.max_cpu_time"] = max_cpu_time

        logger.debug("Maximum CPU time for " + self._solver_name + " set to " + \
            str(max_cpu_time) + " s ...")

    def _setup_nlpsolver(self):

        __dirname__ = os.path.dirname(os.path.abspath(__file__))
        path_to_nlp_object = os.path.join(__dirname__, self._PATH_TO_NLP_OBJECT, \
            self._NLP_OBJECT_FILENAME)

        self._nlpsolver = ca.nlpsol(self._solver_name, "ipopt", path_to_nlp_object,
            self._nlpsolver_options)

    def _set_states_bounds(self):

        '''
        The boundary values for the states will later be defined as soft constraints.
        '''

        self.x_min = self.p_op["T"]["min"] * np.ones( \
            (self._timing.N+1, self.nx - self.nx_aux))
        self.x_max = self.p_op["T"]["max"] * np.ones( \
            (self._timing.N+1, self.nx - self.nx_aux))

        self.x_max[:,self.x_index["T_shx_psc"][-1]] = \
            self.p_op["T_sc"]["T_feed_max"]

    def _set_continuous_control_bounds(self):

        self.u_min = np.hstack([

            self.p_op["v_ppsc"]["min_mpc"] * np.ones((self._timing.N, 1)),
            self.p_op["p_mpsc"]["min_mpc"] * np.ones((self._timing.N, 1)),
            self.p_op["v_plc"]["min_mpc"] * np.ones((self._timing.N, 1)),
            self.p_op["v_pssc"]["min_mpc"] * np.ones((self._timing.N, 1)),

            # The upcoming controls are constrained later in the NLP using inequality constraints

            np.zeros((self._timing.N, 1)),
            np.zeros((self._timing.N, 1))])

        self.u_max = np.hstack([

            self.p_op["v_ppsc"]["max"] * np.ones((self._timing.N, 1)),
            self.p_op["p_mpsc"]["max"] * np.ones((self._timing.N, 1)),
            self.p_op["v_plc"]["max"] * np.ones((self._timing.N, 1)),
            self.p_op["v_pssc"]["max"] * np.ones((self._timing.N, 1)),

            np.inf * np.ones((self._timing.N, 1)),
            np.inf * np.ones((self._timing.N, 1))])

    @abstractmethod
    def _set_binary_control_bounds(self):

        pass

    def _set_nlpsolver_bounds_and_initials(self):

        # Optimization variables bounds and initials

        V_min = []
        V_max = []
        V_init = []

        # Constraints bounds

        lbg = []
        ubg = []

        # Time-varying parameters

        P_data = []

        # Initial states

        if self._timing.grid_position_cursor == 0:

            V_min.append(self._predictor.x_hat)
            V_max.append(self._predictor.x_hat)
            V_init.append(self._predictor.x_hat)

        else:

            V_min.append(self._previous_solver.x_data[0,:])
            V_max.append(self._previous_solver.x_data[0,:])
            V_init.append(self._previous_solver.x_data[0,:])

        for k in range(self._timing.N):

            # Collocation equations

            for j in range(1,self.d+1):

                lbg.append(np.zeros(self.nx))
                ubg.append(np.zeros(self.nx))

            if k < self._timing.grid_position_cursor:

                lbg.append(-np.inf * np.ones(self.nx))
                ubg.append(np.inf * np.ones(self.nx))

            else:

                lbg.append(np.zeros(self.nx))
                ubg.append(np.zeros(self.nx))

            # s_ac_lb

            lbg.append(-1e-1 * np.ones(3)) # vanishing constraints smoothened
            ubg.append(np.inf * np.ones(3))

            # s_ac_ub

            lbg.append(-np.inf * np.ones(3))
            ubg.append(1e-1 * np.ones(3)) # vanishing constraints smoothened

            # Setup objective temperature range condition

            lbg.append(self.p_op["room"]["T_r_a_min"])
            ubg.append(self.p_op["room"]["T_r_a_max"])

            # State limits soft constraints

            lbg.append(self.x_min[k+1,:])
            ubg.append(self.x_max[k+1,:])

            # Assure ppsc is running at high speed when collector temperature is high

            lbg.append(-np.inf * np.ones(4))
            ubg.append(1.0e-1)
            ubg.append(0)
            ubg.append(1.0e-1)
            ubg.append(0)

            # Assure HTS bottom layer mass flows are always smaller or equal to
            # the corresponding total pump flow

            lbg.append(-np.inf * np.ones(2))
            ubg.append(np.zeros(2))

            # SOS1 constraints

            lbg.append(0)
            ubg.append(1)

            # Append new bounds and initials

            for j in range(1,self.d+1):

                V_min.append(-np.inf * np.ones(self.nx))
                V_max.append(np.inf * np.ones(self.nx))
                V_init.append(self._previous_solver.x_data[k,:])

            if  k < self._timing.grid_position_cursor:

                V_min.append(self._previous_solver.b_data[k,:])
                V_max.append(self._previous_solver.b_data[k,:])
                V_init.append(self._previous_solver.b_data[k,:])

            else:

                V_min.append(self.b_min[k,:])
                V_max.append(self.b_max[k,:])
                V_init.append(self._previous_solver.b_data[k,:])

            V_min.append(np.zeros(3))
            V_max.append(np.inf * np.ones(3))
            try:
                V_init.append(self._previous_solver.s_ac_lb_data[k,:])
            except AttributeError:
                V_init.append(np.zeros(3))

            V_min.append(np.zeros(3))
            V_max.append(np.inf * np.ones(3))
            try:
                V_init.append(self._previous_solver.s_ac_ub_data[k,:])
            except AttributeError:
                V_init.append(np.zeros(3))

            V_min.append(-np.inf * np.ones(self.nx-self.nx_aux))
            V_max.append(np.inf * np.ones(self.nx-self.nx_aux))
            try:
                V_init.append(self._previous_solver.s_x_data[k,:])
            except AttributeError:
                V_init.append(np.zeros(self.nx-self.nx_aux))

            V_min.append(np.zeros(2))
            V_max.append(self.p_op["v_ppsc"]["max"] * np.ones(2))
            try:
                V_init.append(self._previous_solver.s_ppsc_fpsc_data[k,:])
            except AttributeError:
                V_init.append(0.0)
            try:
                V_init.append(self._previous_solver.s_ppsc_vtsc_data[k,:])
            except AttributeError:
                V_init.append(0.0)

            if k < self._timing.grid_position_cursor:

                V_min.append(self._previous_solver.u_data[k,:])
                V_max.append(self._previous_solver.u_data[k,:])
                V_init.append(self._previous_solver.u_data[k,:])

            else:

                V_min.append(self.u_min[k,:])
                V_max.append(self.u_max[k,:])
                V_init.append(self._previous_solver.u_data[k,:])

            V_min.append(-np.inf)
            V_max.append(np.inf)

            try:
                V_init.append(self._previous_solver.r_data[k,:])
            except AttributeError:
                V_init.append(0.0)

            if (k+1) == self._timing.grid_position_cursor:

                V_min.append(self._predictor.x_hat)
                V_max.append(self._predictor.x_hat)
                V_init.append(self._predictor.x_hat)

            elif (k+1) < self._timing.grid_position_cursor:

                V_min.append(self._previous_solver.x_data[k+1,:])
                V_max.append(self._previous_solver.x_data[k+1,:])
                V_init.append(self._previous_solver.x_data[k+1,:])

            else:

                V_min.append(-np.inf * np.ones(self.nx))
                V_max.append(np.inf * np.ones(self.nx))
                V_init.append(self._previous_solver.x_data[k+1,:])

            # Append time-varying parameters

            P_data.append(self._ambient.c_data[k,:])
            P_data.append(self._timing.time_steps[k])
            P_data.append(np.zeros(self.nw))

        # Concatenate objects

        self.V_min = ca.veccat(*V_min)
        self.V_max = ca.veccat(*V_max)
        self.V_init = ca.veccat(*V_init)

        self.lbg = np.hstack(lbg)
        self.ubg = np.hstack(ubg)

        self.P_data = ca.veccat(*P_data)

        self._nlpsolver_args = {"p": self.P_data, \
            "x0": self.V_init,
            "lbx": self.V_min, "ubx": self.V_max, \
            "lbg": self.lbg, "ubg": self.ubg}

    def _run_nlpsolver(self):

        logger.info(self._solver_name + ", iter " + \
            str(self._timing.mpc_iteration_count) + ", limit " + \
            str(round(self._nlpsolver_options["ipopt.max_cpu_time"],1)) + " s ...")

        self.nlp_solution = self._nlpsolver(**self._nlpsolver_args)

        if self._nlpsolver.stats()["return_status"] == "Maximum_CpuTime_Exceeded":

            logger.warning(self._solver_name + " returned '" + \
               str(self._nlpsolver.stats()["return_status"]) + "' after " + \
                str(round(self._nlpsolver.stats()["t_wall_total"], 2)) + " s")

        else:

            logger.info(self._solver_name + " returned '" + \
               str(self._nlpsolver.stats()["return_status"]) + "' after " + \
                str(round(self._nlpsolver.stats()["t_wall_total"], 2)) + " s")

    def _collect_nlp_results(self):

        v_opt = np.array(self.nlp_solution["x"])

        x_opt = []
        u_opt = []
        b_opt = []
        r_opt = []
        s_ac_lb_opt = []
        s_ac_ub_opt = []
        s_x_opt = []
        s_ppsc_fpsc_opt = []
        s_ppsc_vtsc_opt = []

        offset = 0

        for k in range(self._timing.N):

            x_opt.append(v_opt[offset:offset+self.nx])

            for j in range(self.d+1):

                offset += self.nx

            b_opt.append(v_opt[offset:offset+self.nb])
            offset += self.nb

            s_ac_lb_opt.append(v_opt[offset:offset+3])
            offset += 3

            s_ac_ub_opt.append(v_opt[offset:offset+3])
            offset += 3

            s_x_opt.append(v_opt[offset:offset+self.nx-self.nx_aux])
            offset += self.nx-self.nx_aux

            s_ppsc_fpsc_opt.append(v_opt[offset:offset+1])
            offset += 1

            s_ppsc_vtsc_opt.append(v_opt[offset:offset+1])
            offset += 1

            u_opt.append(v_opt[offset:offset+self.nu])
            offset += self.nu

            r_opt.append(v_opt[offset:offset+1])
            offset += 1

        x_opt.append(v_opt[offset:offset+self.nx])
        offset += self.nx

        r_opt.append(v_opt[offset:offset+1])
        offset += 1

        self._x_data = ca.horzcat(*x_opt).T
        self._u_data = ca.horzcat(*u_opt).T
        self._b_data = ca.horzcat(*b_opt).T

        self._r_data = ca.horzcat(*r_opt).T

        self._s_ac_lb_data = ca.horzcat(*s_ac_lb_opt).T
        self._s_ac_ub_data = ca.horzcat(*s_ac_ub_opt).T
        self._s_x_data = ca.horzcat(*s_x_opt).T
        self._s_ppsc_fpsc_data = ca.horzcat(*s_ppsc_fpsc_opt).T
        self._s_ppsc_vtsc_data = ca.horzcat(*s_ppsc_vtsc_opt).T

    def solve(self):

        self._setup_nlpsolver()
        self._set_states_bounds()
        self._set_continuous_control_bounds()
        self._set_binary_control_bounds()
        self._set_nlpsolver_bounds_and_initials()
        self._run_nlpsolver()
        self._collect_nlp_results()

    def reduce_object_memory_size(self):

        self._previous_solver = None
        self._predictor = None

    def save_results(self):

        '''
        This function can be used to save the MPC results, possibly including
        solver runtimes, log files etc.
        '''

        pass

class NLPSolverBin(NLPSolverMPCBaseClass):

    def _setup_additional_nlpsolver_options(self):

        self._nlpsolver_options["ipopt.acceptable_tol"] = 0.2
        self._nlpsolver_options["ipopt.acceptable_iter"] = 8
        self._nlpsolver_options["ipopt.acceptable_constr_viol_tol"] = 10.0
        self._nlpsolver_options["ipopt.acceptable_dual_inf_tol"] = 10.0
        self._nlpsolver_options["ipopt.acceptable_compl_inf_tol"] = 10.0
        self._nlpsolver_options["ipopt.acceptable_obj_change_tol"] = 1e-1

        self._nlpsolver_options["ipopt.mu_strategy"] = "adaptive"
        self._nlpsolver_options["ipopt.mu_target"] = 1e-5

    def _set_binary_control_bounds(self):

        self.b_min = self._previous_solver.b_data
        self.b_max = self._previous_solver.b_data

class NLPSolverRel(NLPSolverMPCBaseClass):

    def _setup_additional_nlpsolver_options(self):

        self._nlpsolver_options["ipopt.acceptable_tol"] = 0.2
        self._nlpsolver_options["ipopt.acceptable_iter"] = 8
        self._nlpsolver_options["ipopt.acceptable_constr_viol_tol"] = 10.0
        self._nlpsolver_options["ipopt.acceptable_dual_inf_tol"] = 10.0
        self._nlpsolver_options["ipopt.acceptable_compl_inf_tol"] = 10.0
        self._nlpsolver_options["ipopt.acceptable_obj_change_tol"] = 1e-1

        self._nlpsolver_options["ipopt.mu_strategy"] = "adaptive"
        self._nlpsolver_options["ipopt.mu_target"] = 1e-5

    def _set_binary_control_bounds(self):

        self.b_min = np.zeros((self._timing.N, self.nb))
        self.b_max = np.ones((self._timing.N,self.nb))
        self.b_max[:,-1] = 0.0

        if self._timing._remaining_min_up_time > 0:

            locked_time_grid_points = np.where(self._timing.time_grid < self._timing._remaining_min_up_time)[0]

            self.b_min[locked_time_grid_points, :] = np.repeat(self._timing._b_bin_locked, len(locked_time_grid_points), 0)
            self.b_max[locked_time_grid_points, :] = np.repeat(self._timing._b_bin_locked, len(locked_time_grid_points), 0)

class NLPSolverMHE(NLPSetupMHE):

    @property
    def x_data(self):

        try:

            return np.asarray(self._x_data)

        except AttributeError:

            msg = "Optimized states not available yet, call solve() first."

            logging.error(msg)
            raise RuntimeError(msg)

    @property
    def x_hat(self):

        try:

            return np.squeeze(self._x_data[-1,:])

        except AttributeError:

            msg = "Current states estimate not available yet, call solve() first."

            logging.error(msg)
            raise RuntimeError(msg)

    @property
    def w_data(self):

        try:

            return np.asarray(self._w_data)

        except AttributeError:

            msg = "Optimized process noise not available yet, call solve() first."

            logging.error(msg)
            raise RuntimeError(msg)

    @property
    def u_data(self):

        return self._measurement.u_data

    @property
    def b_data(self):

        return self._measurement.b_data

    @property
    def c_data(self):

        return self._ambient.c_data

    @property
    def solver_name(self):

        return self._solver_name

    def _setup_ambient(self, ambient):

        self._ambient = ambient

    def _setup_measurement(self, measurement):

        self._measurement = measurement

    def _setup_solver_name(self, solver_name):

        self._solver_name = solver_name

    def _setup_solver_options(self):

        self._nlpsolver_options = {}
        self._nlpsolver_options["ipopt.linear_solver"] = "mumps"
        self._nlpsolver_options["ipopt.print_level"] = 5
        self._nlpsolver_options["ipopt.max_cpu_time"] = 1.0e6

    def _initialize_arrival_cost_covariance(self):

        self.P_x_arr = copy.deepcopy(self.R_w)

    def _initialize_arrival_cost(self):

        # Initialize arrival cost with initial measurement

        self.x_arr = np.zeros(self.nx)

        self.x_arr[self.x_index["T_hts"][0]] = \
            self._measurement.y_data[0, self.y_index["T_hts"][0]]

        self.x_arr[self.x_index["T_hts"][1]] = \
            self._measurement.y_data[0, self.y_index["T_hts"][1]]

        self.x_arr[self.x_index["T_hts"][2]] = \
            self._measurement.y_data[0, self.y_index["T_hts"][2]]

        self.x_arr[self.x_index["T_hts"][3]] = \
            self._measurement.y_data[0, self.y_index["T_hts"][3]]

        self.x_arr[self.x_index["T_lts"][0]] = \
            self._measurement.y_data[0, self.y_index["T_lts"][0]]

        self.x_arr[self.x_index["T_lts"][1]] = \
            self._measurement.y_data[0, self.y_index["T_lts"][1]]

        self.x_arr[self.x_index["T_fpsc"]] = \
            self._measurement.y_data[0, self.y_index["T_fpsc_s"]]
        self.x_arr[self.x_index["T_fpsc_s"]] = \
            self._measurement.y_data[0, self.y_index["T_fpsc_s"]]

        self.x_arr[self.x_index["T_vtsc"]] = \
            self._measurement.y_data[0, self.y_index["T_vtsc_s"]]
        self.x_arr[self.x_index["T_vtsc_s"]] = \
            self._measurement.y_data[0, self.y_index["T_vtsc_s"]]

        self.x_arr[self.x_index["T_pscf"]] = \
            self._measurement.y_data[0, self.y_index["T_shx_psc"][1]]
        self.x_arr[self.x_index["T_pscr"]] = \
            self._measurement.y_data[0, self.y_index["T_shx_psc"][0]]

        self.x_arr[self.x_index["T_shx_psc"][:2]] = \
            self._measurement.y_data[0, self.y_index["T_shx_psc"][0]]
        self.x_arr[self.x_index["T_shx_psc"][2:]] = \
            self._measurement.y_data[0, self.y_index["T_shx_psc"][1]]
        self.x_arr[self.x_index["T_shx_ssc"][:2]] = \
            self._measurement.y_data[0, self.y_index["T_shx_ssc"][0]]
        self.x_arr[self.x_index["T_shx_ssc"][2:]] = \
            self._measurement.y_data[0, self.y_index["T_shx_ssc"][1]]

        self.x_arr[self.x_index["T_fcu_w"]] = \
            self._measurement.y_data[0, self.y_index["T_fcu_w"]]
        self.x_arr[self.x_index["T_fcu_a"]] = \
            self._measurement.y_data[0, self.y_index["T_r_a"][0]]

        self.x_arr[self.x_index["T_r_c"]] = \
            self._measurement.y_data[0, self.y_index["T_r_c"]]
        self.x_arr[self.x_index["T_r_a"]] = \
            self._measurement.y_data[0, self.y_index["T_r_a"]]

    def _initialize_ekf_for_arrival_cost_update(self):

        self.hfcn = ca.Function("h", [self.x, self.c], [self.h])
        self.H = self.hfcn.jac()

        ode = {"x": self.x, "p": ca.veccat(self.c, self.u, self.b, self.w), \
            "ode": self.f}

        self.phi = ca.integrator("integrator", "cvodes", ode, \
            {"t0": 0.0, "tf": self._timing.dt_day})
        self.Phi = self.phi.jac()

    def __init__(self, timing, ambient, measurement, solver_name):

        logger.info("Initializing NLP solver " + solver_name + " ...")

        super().__init__(timing)

        self._setup_ambient(ambient = ambient)
        self._setup_measurement(measurement = measurement)

        self._setup_solver_name(solver_name = solver_name)

        self._setup_solver_options()
        self._setup_collocation_options()

        self._setup_model()
        self._initialize_arrival_cost_covariance()
        self._initialize_arrival_cost()
        self._initialize_ekf_for_arrival_cost_update()

        logger.info("NLP solver "  + solver_name + " initialized.")

    def _setup_nlpsolver(self):

        __dirname__ = os.path.dirname(os.path.abspath(__file__))
        path_to_nlp_object = os.path.join(__dirname__, self._PATH_TO_NLP_OBJECT, \
            self._NLP_OBJECT_FILENAME)

        self._nlpsolver = ca.nlpsol(self._solver_name, "ipopt", path_to_nlp_object,
            self._nlpsolver_options)

    def _set_nlpsolver_bounds_and_initials(self):

        # Optimization variables bounds and initials

        V_min = []
        V_max = []
        V_init = []

        # Constraints bounds

        lbg = []
        ubg = []

        # Time-varying parameters

        P_data = []

        # Initial states

        V_min.append(-np.inf * np.ones(self.nx-2))
        V_min.append(np.zeros(2))
        V_max.append(np.inf * np.ones(self.nx-2))
        V_max.append(2.5 * np.ones(2))
        try:
            V_init.append(self._x_data[0,:])
        except AttributeError:
            V_init.append(self.x_arr)

        P_data.append(self.x_arr)
        P_data.append(np.linalg.inv(self.P_x_arr))

        for k in range(self._timing.N):

            # Collocation equations

            for j in range(1,self.d+1):

                lbg.append(np.zeros(self.nx))
                ubg.append(np.zeros(self.nx))

            lbg.append(np.zeros(self.nx))
            ubg.append(np.zeros(self.nx))

            # Append new boundaries and initials

            for j in range(1,self.d+1):

                V_min.append(-np.inf * np.ones(self.nx))
                V_max.append(np.inf * np.ones(self.nx))
                try:
                    V_init.append(self._x_data[k,:])
                except AttributeError:
                    V_init.append(self.x_arr)

            P_data.append(self._measurement.b_data[k,:])
            P_data.append(self._measurement.u_data[k,:])
            P_data.append(self._measurement.y_data[k,:])
            P_data.append(self._ambient.c_data[k,:])
            P_data.append(self._timing.time_steps[k])

            V_min.append(-np.inf * np.ones(self.nw))
            V_max.append(np.inf * np.ones(self.nw))
            V_init.append(np.zeros(self.nw))

            V_min.append(-np.inf * np.ones(self.nx-2))
            V_min.append(np.zeros(2))
            V_max.append(np.inf * np.ones(self.nx-2))
            V_max.append(1.0 * np.ones(2))
            try:
                V_init.append(self._x_data[k+1,:])
            except AttributeError:
                V_init.append(self.x_arr)

        # Append time-varying parameters

        P_data.append(self._measurement.y_data[k+1,:])

        # Concatenate objects

        self.V_min = ca.veccat(*V_min)
        self.V_max = ca.veccat(*V_max)
        self.V_init = ca.veccat(*V_init)

        self.lbg = np.hstack(lbg)
        self.ubg = np.hstack(ubg)

        self.P_data = ca.veccat(*P_data)

        self._nlpsolver_args = {"p": self.P_data, \
            "x0": self.V_init,
            "lbx": self.V_min, "ubx": self.V_max, \
            "lbg": self.lbg, "ubg": self.ubg}

    def _run_nlpsolver(self):

        logger.info("Running " + self._solver_name + " at MPC iteration " + \
            str(self._timing.mhe_iteration_count) + " ...")

        self.nlp_solution = self._nlpsolver(**self._nlpsolver_args)

        logger.info("Solver " + self._solver_name + " finished with '" + \
           str(self._nlpsolver.stats()["return_status"]) + "' after " + \
            str(round(self._nlpsolver.stats()["t_wall_total"], 2)) + " s ...")

    def _collect_nlp_results(self):

        v_opt = np.array(self.nlp_solution["x"])

        x_opt = []
        w_opt = []

        offset = 0

        for k in range(self._timing.N):

            x_opt.append(v_opt[offset:offset+self.nx])

            for j in range(self.d+1):

                offset += self.nx

            w_opt.append(v_opt[offset:offset+self.nw])
            offset += self.nw

        x_opt.append(v_opt[offset:offset+self.nx])
        offset += self.nx

        self._x_data = ca.horzcat(*x_opt).T
        self._w_data = ca.horzcat(*w_opt).T

    def _update_arrival_cost(self):

        # cf. CasADi example mhe_spring_damper.py

        H0 = self.H(i0=self.x_data[0,:], i1=self._ambient.c_data[0,:])["Do0Di0"]
        K = ca.mtimes([self.P_x_arr, H0.T, \
            np.linalg.inv(ca.mtimes([H0, self.P_x_arr, H0.T]) + self.W_r)])
        self.P_x_arr = ca.mtimes((ca.DM.eye(self.nx) - ca.mtimes(K, H0)), self.P_x_arr)

        h0 = self.hfcn(self.x_data[0,:], self._ambient.c_data[0,:])
        self.x_arr = np.squeeze(self.x_arr) + ca.mtimes(K, \
            self._measurement.y_data[0,:] - h0 - ca.mtimes(H0, np.squeeze(self.x_arr) - self.x_data[0,:]))
        self.x_arr = np.asarray(self.phi(x0 = self.x_arr, \
            p = ca.veccat(self._ambient.c_data[0,:], self._measurement.u_data[0,:], \
                self._measurement.b_data[0,:], self.w_data[0,:]))["xf"])

        F = self.Phi(x0 = self.x_data[0,:], \
            p = ca.veccat(self._ambient.c_data[0,:], self._measurement.u_data[0,:], \
                self._measurement.b_data[0,:], self.w_data[0,:]))["DxfDx0"]
        self.P_x_arr = np.asarray(ca.mtimes([F, self.P_x_arr, F.T])) + self.R_w

    def solve(self):

        self._setup_nlpsolver()
        self._set_nlpsolver_bounds_and_initials()
        self._run_nlpsolver()
        self._collect_nlp_results()
        self._update_arrival_cost()

    def save_results(self):

        '''
        This function can be used to save the MHE results, possibly including
        solver runtimes, log files etc.
        '''

        pass

    def publish(self):

        '''
        This function can be used to publish the MHE results, e. g., via a
        message broker, for communication with other processes
        '''

        pass

max_stars_repo_name: adbuerger/stcs-mimpc

max_issues_repo_name: adbuerger/stcs-mimpc

max_forks_repo_name: adbuerger/stcs-mimpc


text:
import numpy as np
import matplotlib.pyplot as plt
import scipy.io as sio
# outs = np.load('Numpy-batch-7-outs.npy')
# masks = np.load('Numpy-batch-7-masks.npy')

final_outs = []
final_masks = []
final_images = []

loss = np.load('./OutMasks-UNet_Loss_bs=4_ep=1_lr=0.001.npy')

for i in range(38):
    outs = np.load(
        './npy-files/out-files/UNetSmall-NonZero-unetsmall-batch-{}-outs.npy'.format(i))
    masks = np.load(
        './npy-files/out-files/UNetSmall-NonZero-unetsmall--batch-{}-masks.npy'.format(i))
    images = np.load(
        './npy-files/out-files/UNetSmall-NonZero-unetsmall--batch-{}-images.npy'.format(i))

    final_outs.append(outs)
    final_masks.append(masks)
    final_images.append(images)
final_outs = np.asarray(final_outs)
final_masks = np.asarray(final_masks)
final_images = np.asarray(final_images)

print(final_images[0].shape)
for i in range(38):
    print(final_images[i].shape)
    plt.imshow(np.squeeze(final_images[i][49, :, :]), cmap='gray')
    plt.show()

# print(final_outs.shape)
# print(final_masks.shape)

# sio.savemat('./mat-files/final_outputs.mat', {'data': final_outs})
# sio.savemat('./mat-files/final_masks.mat', {'data': final_masks})
# sio.savemat('./mat-files/final_images.mat', {'data': final_images})
# for i in range(len(outs)):
#     plt1 = 255 * np.squeeze(outs[i, :, :, :]).astype('uint8')
#     plt2 = 255 * np.squeeze(masks[i, :, :, :]).astype('uint8')
#     print(plt1, plt2)
#     plt.subplot(1, 2, 1)
#     plt.imshow(plt1, cmap='gray')
#     plt.title("UNet Out")
#     plt.subplot(1, 2, 2)
#     plt.imshow(plt2, cmap='gray')
#     plt.title("Mask")

max_stars_repo_name: Rehan-Ahmar/UNet-Zoo

reason: import numpy,import scipy


text:
from __future__ import annotations

import logging
from typing import TYPE_CHECKING, Union

import matplotlib.pyplot as plt
import numpy as np

from .._doc import doc
from ._image import Image
from ._utils import show_image

logger = logging.getLogger(__name__)

if TYPE_CHECKING:
    from .mesh import TriangleMesh

@doc(Image, prefix='Data class for working with 2D image data', shape='(i,j) ')
class Plane(Image, ndim=2):

    def show(self,
             *,
             ax: plt.Axes = None,
             title: str = None,
             **kwargs) -> 'plt.Axes':
        """Plot the image using :mod:`matplotlib`.

        Parameters
        ----------
        ax : matplotlib.axes.Axes, optional
            Axes to use for plotting.
        title : str, optional
            Title for the plot.
        **kwargs
            These parameters are passed to
            :func:`matplotlib.pyplot.imshow`.

        Returns
        -------
        ax : matplotlib.axes.Axes
            Instance of :class:`matplotlib.axes.Axes`
        """
        return show_image(self.image, ax=ax, title=title, **kwargs)

    @doc(show)
    def plot(self, *args, **kwargs):
        return self.show(*args, **kwargs)

    def generate_mesh(self, **kwargs) -> TriangleMesh:
        """Generate mesh from binary (segmented) image.

        Parameters
        ----------
        **kwargs:
            Keyword arguments are passed to
            :func:`nanomesh.plane2mesh`

        Returns
        -------
        mesh : TriangleMesh
            Description of the mesh.
        """
        from nanomesh.image2mesh import plane2mesh
        return plane2mesh(image=self.image, **kwargs)

    def select_roi(self, from_points: np.ndarray = None):
        """Select region of interest in interactive matplotlib figure.

        Parameters
        ----------
        from_points : (n, 2) numpy.ndarray, optional
            List of points that are used as anchors for the roi
            selection.

        Returns
        -------
        roi : `nanomesh.image._roi2d.ROISelector`
            Region of interest object. Bounding box is stored in
            :attr:`roi.bbox`.
        """
        from ._roi2d import ROISelector
        ax = self.show(title='Select region of interest')
        if from_points is not None:
            # reverse columns to match image
            from_points = from_points[:, ::-1]
            ax.scatter(*from_points.T)
        roi = ROISelector(ax, snap_to=from_points)
        return roi

    def crop(self, left: int, top: int, right: int, bottom: int) -> Plane:
        """Crop image to pixel indices.

        Parameters
        ----------
        left, top, right, bottom : int
            Index of pixel delimiting cropping box.

        Returns
        -------
        Plane
            New instance of :class:`Plane`.
        """
        return Plane(self.image[top:bottom, left:right])

    def crop_to_roi(self, bbox: np.ndarray) -> Plane:
        """Crop plane to rectangle defined by bounding box.

        Parameters
        ----------
        bbox : (4,2) numpy.ndarray
            List of points describing region of interest. The bounding box
            may be rotated.

        Returns
        -------
        Plane
            Cropped region as :class:`Plane` object.
        """
        from ._roi2d import extract_rectangle
        cropped = extract_rectangle(self.image, bbox=bbox)
        return Plane(cropped)

    def compare_with_mesh(self, mesh: TriangleMesh) -> 'plt.Axes':
        """Make a plot comparing the image with the given mesh.

        Parameters
        ----------
        mesh : TriangleMesh
            Mesh to compare the image with.

        Returns
        -------
        plt.Axes
        """
        from ..utils import compare_mesh_with_image
        return compare_mesh_with_image(image=self.image, mesh=mesh)

    def compare_with_digitized(self,
                               digitized: Union[np.ndarray, 'Plane'],
                               cmap: str = None,
                               **kwargs) -> 'plt.Axes':
        """Compare image with digitized (segmented) image. Returns a plot with
        the overlay of the digitized image.

        Parameters
        ----------
        digitized : numpy.ndarray, Plane
            Digitized image of the same dimensions to overlay
        cmap : str
            Matplotlib color map for :func:`matplotlib.pyplot.imshow`
        **kwargs
            These parameters are passed to :func:`skimage.color.label2rgb`.

        Returns
        -------
        ax : matplotlib.axes.Axes
        """
        from skimage.color import label2rgb

        if isinstance(digitized, Plane):
            digitized = digitized.image

        # bg_label=0 is default for scikit-image from 0.19 onwards
        kwargs.setdefault('bg_label', 0)
        image_overlay = label2rgb(digitized, image=self.image, **kwargs)

        fig, ax = plt.subplots()

        ax.imshow(image_overlay, interpolation='none', cmap=cmap)
        ax.axis('image')
        ax.set_xticks([])
        ax.set_yticks([])
        ax.set_title('Image comparison')

        return ax

    def compare_with_other(self,
                           other: Union[np.ndarray, 'Plane'],
                           cmap: str = None,
                           **kwargs) -> 'plt.Axes':
        """Compare image with other image.

        Parameters
        ----------
        other : numpy.ndarray, Plane
            Other image of the same dimensions to overlay
        cmap : str
            Matplotlib color map for :func:`matplotlib.pyplot.imshow`
        **kwargs
            These parameters are passed to :func:`skimage.util.compare_images`.

        Returns
        -------
        ax : matplotlib.axes.Axes
        """
        from skimage.util import compare_images

        if isinstance(other, Plane):
            other = other.image

        kwargs.setdefault('method', 'checkerboard')
        kwargs.setdefault('n_tiles', (4, 4))
        comp = compare_images(self.image, other, **kwargs)

        fig, ax = plt.subplots()

        ax.imshow(comp, interpolation='none', cmap=cmap)
        ax.axis('image')
        ax.set_xticks([])
        ax.set_yticks([])
        ax.set_title(f'Image comparison ({kwargs["method"]})')

        return ax

    def clear_border(self, *, object_label: int, fill_val: int,
                     **kwargs) -> Plane:
        """Remove objects at the border of the image.

        Parameters
        ----------
        object_label : int
            Label of the objects to remove.
        fill_val : int
            Cleared objects are set to this value.
        **kwargs
            These parameters are passed to
            :func:`skimage.segmentation.clear_border`.

        Returns
        -------
        Plane
            New instance of :class:`Plane`.
        """
        from skimage import segmentation

        objects = (self.image == object_label).astype(int)
        border_cleared = segmentation.clear_border(objects, **kwargs)
        mask = (border_cleared != objects)

        out = self.image.copy()
        out[mask] = fill_val
        return self.__class__(out)

    def try_all_threshold(self, **kwargs):
        """Produce a plot trying all available thresholds using
        :func:`skimage.filters.try_all_threshold`.

        Parameters
        ----------
        **kwargs
            These parameters are passed to
            :func:`skimage.filters.try_all_threshold`.
        """
        from skimage import filters
        kwargs.setdefault('verbose', False)
        self.apply(filters.try_all_threshold, **kwargs)

max_stars_repo_name: hpgem/nanomesher

max_issues_repo_name: hpgem/nanomesher

max_forks_repo_name: hpgem/nanomesher


text:
import numpy as np

def __print_success_message():
    print('Tests passed')

def test_normalize(normalize):
    test_shape = (np.random.choice(range(1000)), 32, 32, 3)
    test_numbers = np.random.choice(range(256), test_shape)
    normalize_out = normalize(test_numbers)

    assert type(normalize_out).__module__ == np.__name__, \
        'Not Numpy Object'

    assert normalize_out.shape == test_shape, \
        'Incorrect Shape. {} shape found'.format(normalize_out.shape)

    assert normalize_out.max() <= 1 and normalize_out.min() >= 0, \
        'Incorrect Range. {} to {} found'.format(normalize_out.min(), normalize_out.max())

    __print_success_message()

def test_grayscale(grayscale):
    test_shape = (np.random.choice(range(1000)), 32, 32, 3)
    test_numbers = np.random.choice(range(256), test_shape)
    grayscale_out = grayscale(test_numbers)

    assert type(grayscale_out).__module__ == np.__name__, \
        'Not Numpy Object'

    assert len(grayscale_out.shape) < 2,  \
        'Incorrect Shape. {} shape found'.format(grayscale_out.shape)

    __print_success_message()

max_stars_repo_name: sbatururimi/Traffic-Sign-Classifier

max_issues_repo_name: sbatururimi/Traffic-Sign-Classifier

max_forks_repo_name: sbatururimi/Traffic-Sign-Classifier