text:
import numpy as np
import pickle
import os

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import style

from PIL import Image, ImageDraw, ImageFont
from sklearn.cluster import DBSCAN

from scipy.spatial import distance

#pickle.dump(data_for_classifier_flattened, open(raw_path, 'wb'))

# OnNotOn #######################################################################################################
# zl path
# radarData_path = 'F:/onNotOn_data/072819_zl_onNotOn/f_data-2019-07-28_22-11-01.258054_zl_onNotOn_rnn/f_data.p'
# videoData_path = 'F:/onNotOn_data/072819_zl_onNotOn/v_data-2019-07-28_22-10-32.249041_zl_onNotOn_rnn/cam1'
# mergedImg_path = 'F:/config_detection/figures/zl_onNotOn_x03y03z03_clustered_esp02ms4'
# raw_path = 'F:/onNotOn_raw/zl_onNoton_raw_flattened.p'

# ag path
# radarData_path = 'F:/onNotOn_data/072819_ag_onNotOn/f_data-2019-07-28_21-44-17.102820_ag_onNotOn_rnn/f_data.p'
# videoData_path = 'F:/onNotOn_data/072819_ag_onNotOn/v_data-2019-07-28_21-44-08.514321_ag_onNotOn_rnn/cam1'
# mergedImg_path = 'F:/config_detection/figures/ag_onNotOn_x03y03z03_clustered_esp02ms4'
# raw_path = 'F:/onNotOn_raw/ag_onNoton_raw_flattened.p'

# zy path
# radarData_path = 'F:/onNotOn_data/072919_zy_onNotOn/f_data.p'
# videoData_path = 'F:/onNotOn_data/072919_zy_onNotOn/v_data-2019-07-29_11-40-34.810544_zy_onNotOn/cam1'
# mergedImg_path = 'F:/config_detection/figures/zy_onNotOn_x03y03z03_clustered_esp02ms4'
# raw_path = 'F:/onNotOn_raw/zy_onNoton_raw_flattened.p'

# Tracking #######################################################################################################

# Palmpad Test ####################################################################################################
# Angled Path
# radarData_path = 'F:/palmpad/f_data_zy_ABC_angled/f_data.p'
# videoData_path = 'F:/palmpad/v_data_zy_ABC_angled/cam2'
# mergedImg_path = 'F:/palmpad/figures/angled'
# Vertical Path
# radarData_path = 'F:/palmpad/f_data_zy_ABC_vertical/f_data.p'
# videoData_path = 'F:/palmpad/v_data_zy_ABC_vertical/cam2'
# mergedImg_path = 'F:/palmpad/figures/vertical'
# Flat Index Path
# radarData_path = 'F:/palmpad/f_data_zy_ABC_flat_index/f_data.p'
# videoData_path = 'F:/palmpad/v_data_zy_ABC_flat_index/cam2'
# mergedImg_path = 'F:/palmpad/figures/flat_index'
# Flat Thumb Path
# radarData_path = 'F:/palmpad/f_data_zy_ABC_flat_thumb/f_data.p'
# videoData_path = 'F:/palmpad/v_data_zy_ABC_flat_thumb/cam2'
# mergedImg_path = 'F:/palmpad/figures/flat_thumb'
# Writing Classification ####################################################################################################
radarData_path = 'F:/palmpad/f_data_ya_A/f_data.p'
videoData_path = 'F:/palmpad/v_data_ya_A/cam2'
mergedImg_path = 'F:/palmpad/figures/ya_A'

# utility directory to save the pyplots
radar_3dscatter_path = 'F:/palmpad/figures/utils/radar_3dscatter'

radar_data = list(pickle.load(open(radarData_path, 'rb')).items())
radar_data.sort(key=lambda x: x[0])  # sort by timestamp
videoData_list = os.listdir(videoData_path)
videoData_timestamps = list(map(lambda x: float(x.strip('.jpg')), videoData_list))

style.use('fivethirtyeight')
white_color = 'rgb(255, 255, 255)'
black_color = 'rgb(0, 0, 0)'
red_color = 'rgb(255, 0, 0)'

DBSCAN_esp = 0.2
DBSCAN_minSamples = 3

# input data for the classifier that has the shape n*4*100, n being the number of samples
num_padding = 100
data_for_classifier = np.zeros((len(radar_data), num_padding, 4))
data_for_classifier_flattened = np.zeros((len(radar_data), 4 * num_padding + 1 + 1 + 1))  # + 1 + 1 for the timestamp as integer ratio

fnt = ImageFont.truetype("arial.ttf", 16)

for i, radarFrame in enumerate(radar_data):

    timestamp, fData = radarFrame
    print('Processing ' + str(i + 1) + ' of ' + str(len(radar_data)))

    closest_video_timestamp = min(videoData_timestamps,
                                  key=lambda x: abs(x - timestamp))
    closest_video_path = os.path.join(videoData_path, str(closest_video_timestamp) + '.jpg')
    closest_video_img = Image.open(closest_video_path)

    # plot the radar scatter
    ax1 = plt.subplot(2, 2, 1, projection='3d')
    ax1.set_xlim((-0.3, 0.3))
    ax1.set_ylim((-0.3, 0.3))
    ax1.set_zlim((-0.3, 0.3))
    ax1.set_xlabel('X', fontsize=10)
    ax1.set_ylabel('Y', fontsize=10)
    ax1.set_zlabel('Z', fontsize=10)
    ax1.set_title('Detected Points', fontsize=10)
    # plot the detected points
    ax1.scatter(fData['x'], fData['y'], fData['z'], c=fData['doppler'], marker='o')

    # Do DBSCAN cluster ###############
    # Do cluster ###############
    # map the points to their doppler value, this is for retrieving the doppler value after clustering
    data = np.asarray([fData['x'], fData['y'], fData['z'], fData['doppler']]).transpose()
    doppler_dict = {}
    for point in data:
        doppler_dict[tuple(point[:3])] = point[3:]
    # get rid of the doppler for clustering TODO should we consider the doppler in clustering?
    data = data[:, :3]

    db = DBSCAN(eps=DBSCAN_esp, min_samples=DBSCAN_minSamples).fit(data)
    core_samples_mask = np.zeros_like(db.labels_, dtype=bool)
    core_samples_mask[db.core_sample_indices_] = True
    labels = db.labels_
    # Number of clusters in labels, ignoring noise if present.
    n_clusters_ = len(set(labels)) - (1 if -1 in labels else 0)
    n_noise_ = list(labels).count(-1)

    ax2 = plt.subplot(2, 2, 2, projection='3d')
    ax2.set_xlim((-0.3, 0.3))
    ax2.set_ylim((-0.3, 0.3))
    ax2.set_zlim((-0.3, 0.3))
    ax2.set_xlabel('X', fontsize=10)
    ax2.set_ylabel('Y', fontsize=10)
    ax2.set_zlabel('Z', fontsize=10)
    ax2.set_title('Clustered Points', fontsize=10)

    unique_labels = set(labels)
    colors = [plt.cm.Spectral(each)
              for each in np.linspace(0, 1, len(unique_labels))]

    clusters = []

    for k, col in zip(unique_labels, colors):
        if k == -1:
            # Black used for noise.
            col = [0, 0, 0, 1]
        class_member_mask = (labels == k)
        xyz = data[class_member_mask & core_samples_mask]
        # ax.scatter(xyz[:, 0], xyz[:, 1], xyz[:, 2], 'o', c=np.array([col]), s=28, marker='o')
        if xyz.any():  # in case there are none objects
            clusters.append(xyz)  # append this cluster data to the cluster list
        # each cluster is a 3 * n matrix
        xyz = data[class_member_mask & ~core_samples_mask]
        ax2.scatter(xyz[:, 0], xyz[:, 1], xyz[:, 2], 'o', c=np.array([col]), s=12, marker='X')  # plot the noise

    # find the center for each cluster
    clusters_centers = list(
        map(lambda xyz: np.array([np.mean(xyz[:, 0]), np.mean(xyz[:, 1]), np.mean(xyz[:, 2])]), clusters))
    clusters.sort(key=lambda xyz: distance.euclidean((0.0, 0.0, 0.0), np.array(
        [np.mean(xyz[:, 0]), np.mean(xyz[:, 1]), np.mean(xyz[:, 2])])))

    # plot the clusters
    for xyz, col in zip(clusters, colors):
        ax2.scatter(xyz[:, 0], xyz[:, 1], xyz[:, 2], 'o', c=np.array([col]), s=28,
                    marker='o')  # plot the cluster points

    #############################
    # center normalize hand cluster
    # clear the hand cluster
    hand_cluster = []

    if len(clusters) > 0:
        hand_cluster = clusters[0]

        xmean = np.mean(hand_cluster[:, 0])
        xmin = np.min(hand_cluster[:, 0])
        xmax = np.max(hand_cluster[:, 0])

        ymean = np.mean(hand_cluster[:, 1])
        ymin = np.min(hand_cluster[:, 1])
        ymax = np.max(hand_cluster[:, 1])

        zmean = np.mean(hand_cluster[:, 2])
        zmin = np.min(hand_cluster[:, 2])
        zmax = np.max(hand_cluster[:, 2])

        # append back the doppler
        # doppler array for this frame
        point_num = hand_cluster.shape[0]

        doppler_array = np.zeros((point_num, 1))
        for j in range(point_num):
            doppler_array[j:, ] = doppler_dict[tuple(hand_cluster[j, :3])]
        hand_cluster = np.append(hand_cluster, doppler_array,
                                 1)  # TODO this part needs validation, are the put-back dopplers correct?

        # Do the Mean Normalization
        # avoid division by zero, check if all the elements in a column are the same
        if np.all(hand_cluster[:, 0][0] == hand_cluster[:, 0]) or xmin == xmax:
            hand_cluster[:, 0] = np.zeros((point_num))
        else:
            hand_cluster[:, 0] = np.asarray(list(map(lambda x: (x - xmean) / (xmax - xmin), hand_cluster[:, 0])))

        if np.all(hand_cluster[:, 1][0] == hand_cluster[:, 1]) or ymin == ymax:
            hand_cluster[:, 1] = np.zeros((point_num))
        else:
            hand_cluster[:, 1] = np.asarray(list(map(lambda y: (y - ymean) / (ymax - ymin), hand_cluster[:, 1])))

        if np.all(hand_cluster[:, 2][0] == hand_cluster[:, 2]) or zmin == zmax:
            hand_cluster[:, 2] = np.zeros((point_num))
        else:
            hand_cluster[:, 2] = np.asarray(list(map(lambda z: (z - zmean) / (zmax - zmin), hand_cluster[:, 2])))
        # pad to 50

        hand_cluster_padded = np.pad(hand_cluster, ((0, num_padding - point_num), (0, 0)), 'constant',
                                 constant_values=0)
    else:
        hand_cluster_padded = np.zeros((num_padding, 4))

    # flatten hand_cluster and add timestamp information
    hand_cluster_padded_flattened = hand_cluster_padded.reshape(( -1))
    hand_cluster_padded_flattened = np.insert(hand_cluster_padded_flattened, 0, timestamp.as_integer_ratio()[1])
    hand_cluster_padded_flattened = np.insert(hand_cluster_padded_flattened, 0, timestamp.as_integer_ratio()[0])
    hand_cluster_padded_flattened = np.insert(hand_cluster_padded_flattened, 0, timestamp)

    data_for_classifier[i] = hand_cluster_padded
    data_for_classifier_flattened[i] = hand_cluster_padded_flattened

    # plot the normalized closest cluster
    ax3 = plt.subplot(2, 2, 3, projection='3d')
    ax3.set_xlim((-1.0, 1.0))
    ax3.set_ylim((-1.0, 1.0))
    ax3.set_zlim((-1.0, 1.0))
    ax3.set_xlabel('X', fontsize=10)
    ax3.set_ylabel('Y', fontsize=10)
    ax3.set_zlabel('Z', fontsize=10)
    ax3.set_title('CLosest Cluster', fontsize=10)

    ax3.scatter(hand_cluster_padded[:, 0], hand_cluster_padded[:, 1], hand_cluster_padded[:, 2], c=hand_cluster_padded[:, 3], marker='o')

    #############################
    # Combine the three images
    plt.savefig(os.path.join(radar_3dscatter_path, str(timestamp) + '.jpg'))
    radar_3dscatter_img = Image.open(os.path.join(radar_3dscatter_path, str(timestamp) + '.jpg'))

    images = [closest_video_img, radar_3dscatter_img]  # add image here to arrange them horizontally
    widths, heights = zip(*(i.size for i in images))
    total_width = sum(widths)
    max_height = max(heights)
    new_im = Image.new('RGB', (total_width, max_height))
    x_offset = 0
    for im in images:
        new_im.paste(im, (x_offset, 0))
        x_offset += im.size[0]

    timestamp_difference = abs(float(timestamp) - float(closest_video_timestamp))
    draw = ImageDraw.Draw(new_im)

    # draw the timestamp difference on the image
    (x, y) = (20, 10)
    message = "Timestamp Difference, abs(rt-vt): " + str(timestamp_difference)
    draw.text((x, y), message, fill=white_color, font=fnt)
    # draw the timestamp
    (x, y) = (20, 30)
    message = "Timestamp: " + str(timestamp)
    draw.text((x, y), message, fill=white_color, font=fnt)

    # draw the number of points
    (x, y) = (20, 60)
    message = "Number of detected points: " + str(xyz.shape[0])
    draw.text((x, y), message, fill=white_color, font=fnt)

    # draw the number of clusters and number of noise point on the clutter plot
    (x, y) = (20, 80)
    message = "Number of clusters: " + str(n_clusters_)
    draw.text((x, y), message, fill=white_color, font=fnt)
    (x, y) = (20, 100)
    message = "Number of outliers: " + str(n_noise_)
    draw.text((x, y), message, fill=white_color, font=fnt)

    # save the combined image
    new_im.save(os.path.join(mergedImg_path, str(timestamp) + '_' + str(timestamp.as_integer_ratio()[0]) + '_' + str(timestamp.as_integer_ratio()[1]) + '.jpg'))
    plt.close('all')

import pandas as pd

data_for_classifier_flattened = pd.DataFrame(data_for_classifier_flattened)
# data_for_classifier_flattened.to_csv('F:/config_detection/csv/*.csv')

max_stars_repo_name: ApocalyVec/mGesf

max_issues_repo_name: ApocalyVec/mGesf

max_forks_repo_name: ApocalyVec/mGesf

reason: import numpy,from scipy


text:
#!/usr/bin/env python3
# Copyright (c) 2020 Performance Portability authors
# SPDX-License-Identifier: MIT

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import argparse
from statistics import mean, harmonic_mean, median
import math
import sys

def geomean(n):
    # geometric_mean was not added to statistics until Python 3.8
    if (sys.version_info.major > 3) or (
            sys.version_info.major == 3 and sys.version_info.minor >= 8):
        return statistics.geometric_mean(n)
    else:
        return np.prod(n)**(1.0 / float(len(n)))

def harmean(n):
    # Python harmonic_mean returns 0 if there is a 0 in the input
    # Return NaN instead to distinguish from pp
    if 0 in n:
        return np.nan
    return harmonic_mean(n)

def pp(n):
    if 0 in n:
        return 0
    return harmonic_mean(n)

# Set up argument parsing
parser = argparse.ArgumentParser(
    description="Produce table of \"average\" efficiencies")
parser.add_argument(
    'input_file',
    help="CSV file containing performance data")
parser.add_argument(
    'output_file',
    help="Output TeX file")
parser.add_argument(
    '--calc-efficiency',
    action="store_true",
    help="Calculate application efficiency")
parser.add_argument(
    '--input-is-throughput',
    action="store_true",
    help="If calculating application efficiency, then treat the data as throughput (higher is better)")
parser.add_argument(
    '--sort',
    action="store_true",
    help="Sort columns according to performance portability")

args = parser.parse_args()

print('Performance portability metrics')
print()
print('Input file: {}'.format(args.input_file))
print()

# Read in the CSV file as a Pandas DataFrame
data = pd.read_csv(
    args.input_file,
    skipinitialspace=True,
    sep=r',\s+',
    delimiter=',',
    na_values='X')

# In the case of trailing whitespace, the X don't get converted.
# This replaces anything starting with an X to a NaN
data = data.replace(r'^X', np.nan, regex=True)

# Make sure the data is all floating point
data[list(data.columns[1:])] = data[list(data.columns[1:])].apply(pd.to_numeric)

print(data)

# Save a version where NaN is set to 0
data_nona = data.fillna(float(0.0))

if (args.calc_efficiency):
    print("Calculating application efficiency...")

    # Calculate application efficiency
    if (not args.input_is_throughput):
        minimums = data.min(axis=1, skipna=True)
        for col in list(data.columns[1:]):
            data_nona[col] = 100.0 * minimums[:] / data_nona[col]
        data_nona = data_nona.replace([np.inf, -np.inf, np.nan], 0.0)

    else:
        maximums = data.max(axis=1, skipna=True)
        for col in list(data.columns[1:]):
            data_nona[col] = data_nona[col] / maximums[:] * 100.0
        data_nona = data_nona.replace([np.inf, -np.inf, np.nan], 0.0)

else:
    print("Warning: using input data as efficiencies")

# Display data information
print('Number of data items:')
print(data.count())
print()

print(data_nona)

# Compute "averages" for each implementation after discarding non-numeric data
data_nona = data_nona.drop('Device', axis=1)
averages = {"Minimum": min,
            "Arithmetic Mean": mean,
            "Geometric Mean": geomean,
            "Harmonic Mean": harmean,
            "Median": median,
            "Performance Portability": pp}
results = pd.DataFrame()
for (name, f) in averages.items():
    avg = data_nona.apply(f, raw=True).copy()
    avg.name = name
    results = results.append(avg, ignore_index=False)

# Sort columns according to their PP value
# sort_index is not supported by old Pandas
if args.sort:
    order = sorted([col for col in results.columns],
                   key=lambda col: results[col]['Performance Portability'])
    results = results.reindex(order, axis=1)

print()
print(results)

# Write table to LaTeX file
results.to_latex(args.output_file, float_format="%.2f")

print(80 * '-')
print()
print()

ext: py

max_stars_repo_name: UoB-HPC/performance-portability

max_issues_repo_name: UoB-HPC/performance-portability

max_forks_repo_name: UoB-HPC/performance-portability

reason: import numpy


text:
from __future__ import division
import numpy as np
import copy

def SPSA(fnc, initparams, a_par = 1e-6, c_par = .01, args = (), \
         bounds = None, param_tol = None, ftol = 1e-8, maxiter = 10000,
         alpha = .602, gamma = .101,print_iters = 100.):

    def calc(ps): return fnc(ps, *args)

    big_a = maxiter / 10.

    n_iter = 0
    initparams = np.array(initparams)
    n_params = initparams.shape[0]

    p = initparams
    old_val = calc(p)
    saved_p = p*100.

    if bounds is not None:
        if len(bounds) != n_params: raise Exception('Malformed bounds')
        minvals = [x[0] for x in bounds]
        maxvals = [x[1] for x in bounds]
    else:
        minvals = None
        maxvals = None

    while (np.linalg.norm(saved_p - p) / np.linalg.norm(saved_p)) > ftol and n_iter < maxiter:
        saved_p = copy.copy(p)
        ak = a_par / (n_iter + 1 + big_a)**alpha
        ck = c_par / (n_iter + 1)**gamma

        delta = (np.random.randint(0,2,n_params) * 2 - 1)
        p_plus = p + ck*delta
        p_minus = p - ck*delta
        if maxvals:
            p_plus = np.minimum(p_plus,maxvals)
            p_minus = np.minimum(p_minus,maxvals)
        if minvals:
            p_plus = np.maximum(p_plus,minvals)
            p_minus = np.maximum(p_minus,minvals)

        try:
            grad = (calc(p_plus) - calc(p_minus)) / (2.*ck*delta)
        except:
            print p_plus, p_minus, ck, delta
            raise Exception('Stopping')

        # Scale parameter movement
        this_ak = np.ones(n_params)*ak
        new_ps = p
        running = True
        while running:
            ntry = new_ps - this_ak*grad
            oob = np.where(np.logical_or(ntry > maxvals, ntry < minvals))[0]
            new_ps = p - this_ak*grad
            if len(oob) == 0:
                p -= this_ak*grad
                running = False
            else:
                this_ak[oob] = this_ak[oob]/2.

        new_val = calc(p)
        old_val = new_val
        if n_iter % print_iters == 0:
            print "\tIter %05d" % n_iter, new_val, ak, ck, p

        if param_tol is not None:
            pdiff = np.abs(p - saved_p)
            if not np.all(pdiff < param_tol):
                p = saved_p
                continue
        n_iter += 1

    return (p, new_val,n_iter)

max_stars_repo_name: kasmith/btWorld

max_issues_repo_name: kasmith/btWorld

max_forks_repo_name: kasmith/btWorld


text:
#!/usr/bin/python3

# generic imports
import sys
import time

# ROS imports
import rospy

# MoveIt imports
import moveit_commander
import moveit_msgs.msg

# message imports
from dh_gripper_msgs.msg import GripperState

# vision imports
import cv2
import pyrealsense2 as rs

# math related imports
from math import pi
import numpy as np
from rospy.core import rospyinfo
import tf
from tf.transformations import quaternion_matrix
from tf.transformations import translation_matrix
from scipy.spatial.transform import Rotation as R

# speech recognition imports
import speech_recognition as sr
import pyttsx3

class RecognizeColor:
	"""
	A utility class to recognize color name and obtain the coordinates of
	the centroid of the object's contour. It uses speech recognition and
	keyboard input mutually exclusively.
	"""

	def __init__(self, use_speech_recognition = False):
		"""
		Initialize the tools for recognition
		"""
		self.use_speech_recognition = use_speech_recognition
		if use_speech_recognition:
			self.r = sr.Recognizer()
			self.m = sr.Microphone()

		# TODO: add parameterized constructor

		# setting the camera limits
		self.depth_resolution_x = 640
		self.depth_resolution_y = 480
		self.color_resolution_x = 640
		self.color_resolution_y = 480
		self.depth_fps = 30
		self.color_fps = 60

		# minimum and maximum depth to perceive
		self.depth_min = 0.11
		self.depth_max = 1.0

		self.low_H = 0
		self.low_S = 0
		self.low_V = 112
		self.high_H = 180
		self.high_S = 74
		self.high_V = 255

	def get_color_from_speech(self):
		"""
		Returns the color value from speech

		1 -> Red | 2 -> Green | 3 -> Yellow
		"""
		color = 0
		if self.use_speech_recognition:
			print("A moment of silence, please...")
			with self.m as source: self.r.adjust_for_ambient_noise(source)
			print("Set minimum energy threshold to {}".format(self.r.energy_threshold))
			while color==0:
				print("Say something!")
				engine = pyttsx3.init()
				engine.say("Please specify the color")
				engine.runAndWait()
				with self.m as source: audio = self.r.listen(source)
				print("Got it! Now to recognize it...")
				try:
					# recognize speech using Google Speech Recognition
					value = self.r.recognize_google(audio)

					# we need some special handling here to correctly print unicode characters to standard output
					if str is bytes:  # this version of Python uses bytes for strings (Python 2)
						print(u"You said {}".format(value).encode("utf-8"))
						if value == "red":
							color = 1
						elif value == "green":
							color = 2
						elif value == "yellow":
							color = 3
					else:  # this version of Python uses unicode for strings (Python 3+)
						print("You said {}".format(value))
				except sr.UnknownValueError:
					print("Oops! Didn't catch that")
				except sr.RequestError as e:
					print("Uh oh! Couldn't request results from Google Speech Recognition service; {0}".format(e))
			return color
		else:
			print("Cannot use speech recognition")
			return color

	def get_color_from_keyboard(self):
		"""
		Utility method to aid providing color value from keyboard

		1 -> Red | 2 -> Green | 3 -> Yellow
		"""
		if not self.use_speech_recognition:
			return input("Enter the color code: [1] Red | [2] Green | [3] Yellow : ")
		else:
			print("Cannot use keyboard input")
			return 0

	def get_camera_coordinates(self):
		"""
		Returns the 3D camera coordinate of the centroid of the object
		"""
		color_code = 0

		if self.use_speech_recognition:
			color_code = self.get_color_from_speech()
		else:
			color_code = self.get_color_from_keyboard()

		if color_code == 1:
			self.low_H = 0
			self.low_S = 0
			self.low_V = 112
			self.high_H = 180
			self.high_S = 74
			self.high_V = 255
		elif color_code == 2:
			self.low_H = 0
			self.low_S = 81
			self.low_V = 0
			self.high_H = 180
			self.high_S = 152
			self.high_V = 51

		# creating a pipeline to get depth and color images
		pipeline = rs.pipeline()
		config = rs.config()
		config.enable_stream(rs.stream.depth, self.depth_resolution_x, self.depth_resolution_y, rs.format.z16, self.depth_fps)
		config.enable_stream(rs.stream.color, self.color_resolution_x, self.color_resolution_y, rs.format.bgr8, self.color_fps)
		profile = pipeline.start(config)
		depth_sensor = profile.get_device().first_depth_sensor()
		depth_scale = depth_sensor.get_depth_scale()

		camera_coord = []

		while True:
			# wait for a new frame and then get the depth and color frame
			frames = pipeline.wait_for_frames()
			depth_frame = frames.get_depth_frame()
			color_frame = frames.get_color_frame()
			if not depth_frame or not color_frame:
				continue

			# temp var to store camera coordinate; this coordinate might be discarded by the logic in this code that follows
			temp_camera_coord = []

			# create numpy array of depth and color frames
			depth_image = np.asanyarray(depth_frame.get_data())
			color_image = np.asanyarray(color_frame.get_data())

			# cv2.imwrite("colour.jpg", color_image)
			depth_image = depth_image		# TODO: What does this do?

			cv2.imshow('org', color_image)  						# displaying the image

			# check is the requested color is present in the image
			in_range_mask = cv2.inRange(color_image, (self.low_H, self.low_S, self.low_V), (self.high_H, self.high_S, self.high_V))

			# displaying the binary image
			# cv2.imshow('Binary', in_range_mask)

			# obtain a contour of the mask
			contours, _ = cv2.findContours(in_range_mask, cv2.RETR_TREE, cv2.CHAIN_APPROX_TC89_L1)

			# loop through all contours
			for c in contours:
				if cv2.contourArea(c) < 300:  # ignore the contour area when it is less than
					continue
				elif cv2.contourArea(c) > 7000:  # ignore the contour area when it is greater than
					continue
				else:
					M = cv2.moments(c)  	# image moment of the contour
					if M["m00"] != 0:
						px = int(M["m10"] / M["m00"])  # x coordinate centroid
						py = int(M["m01"] / M["m00"])  # y coordinate centroid
						cv2.drawContours(color_image, [c], -1, (0, 255, 0), 2)
						depth_intrin = profile.get_stream(rs.stream.depth).as_video_stream_profile().get_intrinsics()
						color_intrin = profile.get_stream(rs.stream.color).as_video_stream_profile().get_intrinsics()

						depth_to_color_extrin =  profile.get_stream(rs.stream.depth).as_video_stream_profile().get_extrinsics_to( profile.get_stream(rs.stream.color))
						color_to_depth_extrin =  profile.get_stream(rs.stream.color).as_video_stream_profile().get_extrinsics_to( profile.get_stream(rs.stream.depth))

						depth = depth_image[py, px].astype(float)
						# z = depth * depth_scale

						# obtain the depth pixel at the centroid of the contour
						depth_point = rs.rs2_project_color_pixel_to_depth_pixel(
													depth_frame.get_data(), depth_scale, self.depth_min, self.depth_max,
													depth_intrin, color_intrin, depth_to_color_extrin, color_to_depth_extrin, [px,py])

						# centeroid in camera coordinate
						# TODO: What does this logic do?
						temp_camera_coord = rs.rs2_deproject_pixel_to_point(depth_intrin, depth_point, depth)
						temp_camera_coord = np.reshape(temp_camera_coord,(-1,1))
						temp_camera_coord = np.concatenate([temp_camera_coord,np.reshape([1],(-1,1))], axis=0)
					else:
						px, py = 0, 0

					cv2.imshow('Contours', color_image)  # displaying the color image with contours

			# to eliminate the cases when no camera coordinate is obtained
			if len(temp_camera_coord)!=0:
				# to eliminate the cases when depth estimated is more than the allowable range
				if temp_camera_coord[3][0] >= self.depth_max * 1000:	# depth_max parameter is in meters while coordinates are in mm
					cv2.waitKey(1)
					continue
				# eliminating condition that camera is touching the table
				if temp_camera_coord[3][0] <= 0:
					cv2.waitKey(1)
					continue
			else:
				cv2.waitKey(1)
				continue
			print("Press 'q' to perform next step")

			if cv2.waitKey(30)& 0xFF == ord('q'):
				camera_coord = temp_camera_coord
				break

		cv2.destroyAllWindows()		# after pressing any key close all the window

		# the camera coordinates are obtained in mm while the transformation matrix values
		# are obtained in m. Hence, the coordinates are converted from mm to m while
		# preserving the last element in the vector.
		camera_coord = camera_coord/1000
		camera_coord[3][0] = 1

		return camera_coord

class UR5Control:
	"""
	A utility class containing methods that aid in moving the UR5 arm
	as well as the AG95 gripper. It also contains methods to estimate
	the next position to reach.
	"""

	def __init__(self, gripper_offset=200.0, z_safety=200.0):
		"""
			Parameters
			==========

			- gripper_offset: int : Distance of gripper centre point from tool0 link (in mm)
			- z_safety: int : Minimum distance of tool0 from the table (in mm)
		"""
		moveit_commander.roscpp_initialize(sys.argv)
		rospy.init_node('moveToObj', anonymous=True)

		self.gripper_pub = rospy.Publisher('gripper_pose_1', GripperState, queue_size=1)
		self.gripper_msg = GripperState()

		self.robot = moveit_commander.RobotCommander()
		self.scene = moveit_commander.PlanningSceneInterface()
		group_name = "manipulator"
		self.group = moveit_commander.MoveGroupCommander(group_name)

		# self.scene.add_

		print("Pausing (init)...")
		time.sleep(1)
		print("...done!")

		self.display_trajectory_publisher = rospy.Publisher('/move_group/display_planned_path',
												moveit_msgs.msg.DisplayTrajectory,
												queue_size=20)

		self.z_safety = z_safety
		self.gripper_offset = gripper_offset

	def get_base_link_coordinates(self, camera_coord):
		"""
		Returns the coordinate of the object's centeroid wrt base_link from the camera coordinate.
		"""
		# translation
		trans = 0

		# quaternion representing the rotation
		quat = 0

		# object to listens to the transformations
		listener = tf.TransformListener()

		while not rospy.is_shutdown():
			try:
				# We obtain a vector describing the translation and a quaternion describing the rotation
				# as soon as something is returned, break the loop
				trans, quat = listener.lookupTransform('/base_link', '/camera_link2', rospy.Time(0))
				break
			except:
				continue

		# Obtain homogenous matrices from translation vector and the quaternion
		# Multiply the homogenous transformations to get a combined HT
		combinedHT = np.matmul(translation_matrix(trans), quaternion_matrix(quat))

		base_link_coord = np.matmul(np.array(combinedHT), camera_coord)

		return base_link_coord

	def gripper_open(self):
		"""
		Opens the gripper
		"""
		if not rospy.is_shutdown():
			self.gripper_msg.target_position = 1000
			self.gripper_msg.target_force = 20

			self.gripper_pub.publish(self.gripper_msg)

			# after sending the open gripper command, sleep for 2 seconds
			print("Pausing (gripper_open)...")
			time.sleep(2)
			print("...done!")

	def gripper_close(self):
		"""
		Closes the gripper
		"""
		if not rospy.is_shutdown():
			self.gripper_msg.target_position = 20
			self.gripper_msg.target_force = 20

			self.gripper_pub.publish(self.gripper_msg)

			# after sending the close gripper command, sleep for 2 seconds
			print("Pausing (gripper_close)...")
			time.sleep(2)
			print("...done!")

	def get_goal_pose(self, base_link_coord):
		"""
		Converts the base link coordinate location to goal pose.
		A new parameter maybe added to this method to provide orientation information

		Parameters
		=====
		- base_coord : vector describing the goal bose coordinates

		Returns
		---
		Goal pose using input parameters
		"""
		pose_goal = self.group.get_current_pose().pose

		pose_goal.position.x = base_link_coord[0][0]
		pose_goal.position.y = base_link_coord[1][0]
		pose_goal.position.z = base_link_coord[2][0] + (self.gripper_offset/1000)

		if pose_goal.position.z < (self.z_safety/1000):
			pose_goal.position.z = (self.z_safety/1000)
			print("Z pose set to Z-safety limit!")

		pose_goal.orientation.x = -0.926984315176
		pose_goal.orientation.y = 0.374568844855
		pose_goal.orientation.z = -0.00467163233836
		pose_goal.orientation.w = 0.019401951777

		return pose_goal

	def move_to_location(self, goal_pose):
		"""
		Move the EE to the goal pose

		Parameters
		---
		- goal_pose: pose to which path must be planned
		"""
		current_pose = self.group.get_current_pose().pose
		print("Planning path from pose (current pose):  ")
		print(current_pose)

		print("To pose (goal pose): ")
		print(goal_pose)

		self.group.set_pose_target(goal_pose)

		plan = self.group.plan()

		print("Planning the path. Check RViz for visualization. Path will be executed after 5 seconds...")
		rospy.sleep(5)
		print("Executing path...")

		self.group.go(wait=True)
		self.group.stop()
		self.group.clear_pose_targets()

	def get_place_location_pose(self):
		"""
		Returns
		---
		Returns the pose at which the object must to placed.
		"""
		pose_goal = self.group.get_current_pose().pose

		pose_goal.position.x = 0.697925644962
		pose_goal.position.y = 0.102895313812
		pose_goal.position.z = 0.225430132411

		pose_goal.orientation.x = 0.720109681727
		pose_goal.orientation.y = -0.692563162952
		pose_goal.orientation.z = -0.0219631139694
		pose_goal.orientation.w = 0.0362757939835

		return pose_goal

if __name__ == '__main__':
	ur_control = UR5Control(200.0)

	ur_control.gripper_open()

	recognize_color = RecognizeColor()

	camera_coord  = recognize_color.get_camera_coordinates()

	print("Camera coordinates: ", camera_coord)

	base_link_coord = ur_control.get_base_link_coordinates(camera_coord)

	print("Base link coordinates: ", base_link_coord)

	initial_pose = ur_control.group.get_current_pose().pose

	ur_control.move_to_location(ur_control.get_goal_pose(base_link_coord))

	ur_control.gripper_close()

	ur_control.move_to_location(ur_control.get_place_location_pose())

	ur_control.gripper_open()

	ur_control.move_to_location(initial_pose)

	print("Done!")


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

from __future__ import (print_function, unicode_literals,
                        division, absolute_import)

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

def peval(x, A, k, theta):
    return A * np.sin((2 * np.pi * k * x) + theta)

if __name__ == '__main__':
    x = np.arange(0, 6e-2, 6e-2 / 30)
    p = A, k, theta = 10, 1.0 / 3e-2, np.pi / 6
    y_true = peval(x, *p)
    y_meas = y_true + 2 * np.random.randn(len(x))

    p0 = [8, 1 / 2.3e-2, np.pi / 3]
    print('Initial guess for A, k , theta =', np.array(p0))

    popt, pcov = opt.curve_fit(peval, x, y_meas, p0)
    print('Best fit values of A, k, theta =', popt)
    print('True values of A, k, theta =', np.array([A, k, theta]))

    plt.plot(x, peval(x, *popt), x, y_meas, 'o', x, y_true)
    plt.title('Least-squares fit to noisy data')
    plt.legend(['Fit', 'Noisy', 'True'])
    plt.show()

max_stars_repo_name: lneisenman/pyhints

max_issues_repo_name: lneisenman/pyhints

max_forks_repo_name: lneisenman/pyhints

reason: import numpy,import scipy


text:
import numpy as np
import cv2

cap = cv2.VideoCapture(0)
count = 1
while count != 1001:
    ret, frame = cap.read()
    # gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
    # blur = cv2.GaussianBlur(gray, (5, 5), 2)
    # th3 = cv2.adaptiveThreshold(blur, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY_INV, 11, 2)
    # ret, res = cv2.threshold(th3, minValue, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU)

    value = (33, 33)
    hsv = cv2.cvtColor(frame,cv2.COLOR_BGR2HLS)
    blur = cv2.GaussianBlur(hsv,value,0)
    # blurred = cv2.GaussianBlur(grey, value, 0)
    lower_green = np.array([20, 0.8 * 255, 0.6 * 255])
    upper_green = np.array([255, 0.8 * 255, 0.6 * 255])
    mask = cv2.inRange(hsv, lower_green, upper_green)
    gaussian = cv2.GaussianBlur(mask, (11,11), 0)
    # image, contours, hierarchy = cv2.findContours(thresh1.copy(),cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE)
    # cnt = max(contours, key = lambda x: cv2.contourArea(x))
    # x, y, w, h = cv2.boundingRect(cnt)
    # cv2.rectangle(crop_img, (x, y), (x+w, y+h), (0, 255, 255), 0)
    erosion = cv2.erode(mask, None, iterations = 1)
    dilated = cv2.dilate(erosion,None,iterations = 1)
    median = cv2.medianBlur(dilated, 7)
    cv2.imshow('cropped', frame)
    cv2.imshow('mask', median)
    # #
    # write_img = cv2.resize(median, (50,50))
    # cv2.imwrite('images_data/peace/'+str(count)+'.jpg',write_img)
    # print count
    # count += 1
    k = cv2.waitKey(5) & 0xFF
    if k == 27:
        break

cv2.destroyAllWindows()
cap.release()


text:
from .. import datasets
from . import scorers
from . import models
from . import methods
from .. import __version__
import numpy as np
import sklearn
import os
import pickle
import sys
import time
import subprocess
from multiprocessing import Pool
import itertools
from shap.benchmark.methods import linear_regress, tree_regress, deep_regress

all_scorers = [
    "runtime",
    "local_accuracy",
    "consistency_guarantees",
    "mask_keep_positive",
    "mask_keep_negative",
    "keep_positive",
    "keep_negative",
    "batch_keep_absolute_r2",
    "mask_remove_positive",
    "mask_remove_negative",
    "remove_positive",
    "remove_negative",
    "batch_remove_absolute_r2"
]

_experiments = []
_experiments += [["corrgroups60", "lasso", m, s] for s in all_scorers for m in linear_regress]
_experiments += [["corrgroups60", "ridge", m, s] for s in all_scorers for m in linear_regress]
_experiments += [["corrgroups60", "decision_tree", m, s] for s in all_scorers for m in tree_regress]
_experiments += [["corrgroups60", "random_forest", m, s] for s in all_scorers for m in tree_regress]
_experiments += [["corrgroups60", "gbm", m, s] for s in all_scorers for m in tree_regress]
_experiments += [["corrgroups60", "ffnn", m, s] for s in all_scorers for m in deep_regress]

def experiments(dataset=None, model=None, method=None, scorer=None):
    for experiment in _experiments:
        if dataset is not None and dataset != experiment[0]:
            continue
        if model is not None and model != experiment[1]:
            continue
        if method is not None and method != experiment[2]:
            continue
        if scorer is not None and scorer != experiment[3]:
            continue
        yield experiment

def run_experiment(experiment, use_cache=True, cache_dir="/tmp"):
    dataset_name, model_name, method_name, scorer_name = experiment

    # see if we have a cached version
    cache_id = "v" + "__".join([__version__, dataset_name, model_name, method_name, scorer_name])
    cache_file = os.path.join(cache_dir, cache_id + ".pickle")
    if use_cache and os.path.isfile(cache_file):
        with open(cache_file, "rb") as f:
            #print(cache_id.replace("__", " ") + " ...loaded from cache.")
            return pickle.load(f)

    # compute the scores
    print(cache_id.replace("__", " ") + " ...")
    sys.stdout.flush()
    start = time.time()
    X,y = getattr(datasets, dataset_name)()
    score = getattr(scorers, scorer_name)(
        X, y,
        getattr(models, dataset_name+"__"+model_name),
        method_name
    )
    print("...took %f seconds.\n" % (time.time() - start))

    # cache the scores
    with open(cache_file, "wb") as f:
        pickle.dump(score, f)

    return score

def run_experiments_helper(args):
    experiment, cache_dir = args
    return run_experiment(experiment, cache_dir=cache_dir)

def run_experiments(dataset=None, model=None, method=None, scorer=None, cache_dir="/tmp", nworkers=1):
    experiments_arr = list(experiments(dataset=dataset, model=model, method=method, scorer=scorer))
    if nworkers == 1:
        out = list(map(run_experiments_helper, zip(experiments_arr, itertools.repeat(cache_dir))))
    else:
        with Pool(nworkers) as pool:
            out = pool.map(run_experiments_helper, zip(experiments_arr, itertools.repeat(cache_dir)))
    return list(zip(experiments_arr, out))

max_stars_repo_name: jmschrei/shap

max_issues_repo_name: jmschrei/shap

max_forks_repo_name: jmschrei/shap


text:
import string
from typing import List, Optional, Callable, Tuple, Dict, Set

import langcodes
import nltk
import numpy as np
import yake
from flair.data import Sentence
from flair.embeddings import WordEmbeddings
from lemmagen import (
    DICTIONARY_SLOVENE,
    DICTIONARY_ENGLISH,
    DICTIONARY_SERBIAN,
    DICTIONARY_ITALIAN,
    DICTIONARY_ROMANIAN,
    DICTIONARY_HUNGARIAN,
    DICTIONARY_FRENCH,
    DICTIONARY_GERMAN,
    DICTIONARY_SPANISH,
    DICTIONARY_CZECH,
    DICTIONARY_BULGARIAN,
    DICTIONARY_ESTONIAN,
)
from lemmagen.lemmatizer import Lemmatizer
from nltk.corpus import stopwords
from sklearn.feature_extraction.text import TfidfVectorizer

from textsemantics.textrank import keywords as textrank_kw
from textsemantics.utils import cos_sim
from textsemantics.utils.udpipe import get_udipipe_lematizer
from textsemantics.utils.word_enrichment import hypergeom_p_values

lemmagen_languages = {
    "slovenian": DICTIONARY_SLOVENE,
    "english": DICTIONARY_ENGLISH,
    "serbian": DICTIONARY_SERBIAN,
    "italian": DICTIONARY_ITALIAN,
    "romanian": DICTIONARY_ROMANIAN,
    "hungarian": DICTIONARY_HUNGARIAN,
    "french": DICTIONARY_FRENCH,
    "german": DICTIONARY_GERMAN,
    "spanish": DICTIONARY_SPANISH,
    "czech": DICTIONARY_CZECH,
    "bulgarian": DICTIONARY_BULGARIAN,
    "estonian": DICTIONARY_ESTONIAN,
}

ft_languages = {
    "slovenian": 'sl',
    "english": 'en',
    "serbian": 'sr',
    "italian": 'it',
    "romanian": 'ro',
    "french": 'fr',
    "german": 'de',
    "spanish": 'es',
    "czech": 'cs',
    "bulgarian": 'bg',
}

def _get_lemmatizer(language: str) -> Callable:
    if language in lemmagen_languages:
        return Lemmatizer(
            dictionary=lemmagen_languages[language.lower()]
        ).lemmatize
    else:
        return get_udipipe_lematizer(language)

def nltk_language(language: str) -> str:
    exceptions = {"slovenian": "slovene"}
    return exceptions.get(language, language)

def _preprocess_corpus(corpus: List[str], language: str) -> List[List[str]]:
    nltk.download("stopwords", quiet=True)
    lemmatizer = _get_lemmatizer(language)
    stop_words = set(stopwords.words(nltk_language(language.lower())))
    tokenizer = nltk.RegexpTokenizer("\w+")

    preprocessed = list()
    for text in corpus:
        text = text.translate(text.maketrans("", "", string.punctuation))
        tokens = tokenizer.tokenize(text.lower())
        tokens = [
            lemmatizer(token)
            for token in tokens
            if token not in stop_words
            and len(token) > 2
            and not token.isnumeric()
        ]
        # lematizer in rare cases produce empty strings - removing them
        tokens = list(filter(lambda a: a != "", tokens))
        preprocessed.append(tokens)
    return preprocessed

def prepare_embeddings(
    tokens_list: List[List[str]], model_name: str = 'sl'
) -> Tuple[
    np.ndarray, Dict[str, np.ndarray], Dict[str, Set[int]], List[Set[str]]
]:
    embedder = WordEmbeddings(model_name)
    word_embs = {}
    doc_embs = list()
    doc2word = list()
    word2doc = dict()

    for i, tokens in enumerate(tokens_list):
        sent = Sentence(" ".join(tokens))
        embedder.embed(sent)
        doc_emb = np.zeros(embedder.embedding_length)
        doc2word.append(set())
        for token in sent.tokens:
            if token.text not in word2doc:
                word2doc[token.text] = set()
            word2doc[token.text].add(i)
            doc2word[i].add(token.text)

            if token.text not in word_embs:
                emb = token.embedding.cpu().detach().numpy()
                word_embs[token.text] = emb
            else:
                emb = word_embs[token.text]

            doc_emb += emb / len(tokens)
        doc_embs.append(doc_emb)

    doc_embs = np.array(doc_embs)

    return doc_embs, word_embs, word2doc, doc2word

def find_corpus_words(
    doc_embs: np.ndarray, word_embs: Dict[str, np.ndarray]
) -> List[List[Tuple[str, float]]]:
    # compute distances
    distances = np.zeros((len(word_embs), doc_embs.shape[0]))
    words = list(word_embs.keys())
    for i, word in enumerate(words):
        for j in range(doc_embs.shape[0]):
            distances[i, j] = 1 - cos_sim(word_embs[word], doc_embs[j, :])

    # compute scores
    doc_desc = []
    for j in range(doc_embs.shape[0]):
        scores = np.zeros(len(word_embs))
        for i in range(len(words)):
            mask = np.full(doc_embs.shape[0], fill_value=True)
            mask[j] = False
            scores[i] = distances[i, j] - np.mean(distances[i, mask])

        idx = np.argsort(scores)
        doc_desc.append([(words[w], scores[w]) for w in idx])
    return doc_desc

def find_document_words(
    doc_embs: np.ndarray,
    word_embs: Dict[str, np.ndarray],
    word2doc: Dict[str, Set[int]],
    doc2word: List[Set[str]],
) -> List[List[Tuple[str, float]]]:
    distances = dict()
    for word, w_emb in word_embs.items():
        for j in word2doc[word]:
            distances[word, j] = 1 - cos_sim(w_emb, doc_embs[j, :])

    doc_desc = list()
    for j in range(doc_embs.shape[0]):
        scores = np.zeros(len(doc2word[j]))
        ind2word = dict(zip(list(range(len(doc2word[j]))), list(doc2word[j])))
        for k, word in enumerate(doc2word[j]):
            # current document and a single word
            sum_of_distances = sum([distances[word, x] for x in word2doc[word]])
            if len(word2doc[word]) > 1:
                mean_distance = (sum_of_distances - distances[word, j]) / (
                    len(word2doc[word]) - 1
                )
            else:
                mean_distance = 0
            scores[k] = distances[word, j] - mean_distance
        idx = np.argsort(scores)
        doc_desc.append([(ind2word[x], scores[x]) for x in idx])
    return doc_desc

def _get_model_name(language: str) -> str:
    lower = language.lower()
    return ft_languages[lower] if lower in ft_languages else lower

def embedding_corpus_keywords(
    texts: Optional[List[str]] = None,
    tokens: Optional[List[List[str]]] = None,
    language: str = "slovenian",
):
    assert bool(texts) != bool(
        tokens
    ), "Parametri naj vsebujejo zgolj besedilo ali zgolj pojavnice"
    if not tokens:
        tokens = _preprocess_corpus(texts, language)
    doc_embs, word_embs, _, _ = prepare_embeddings(tokens, _get_model_name(language))
    return find_corpus_words(doc_embs, word_embs)

def embedding_document_keywords(
    texts: Optional[List[str]] = None,
    tokens: Optional[List[List[str]]] = None,
    language: str = "slovenian",
):
    assert bool(texts) != bool(
        tokens
    ), "Parametri naj vsebujejo zgolj besedilo ali zgolj pojavnice"
    if not tokens:
        tokens = _preprocess_corpus(texts, language)
    doc_embs, word_embs, word2doc, doc2word = prepare_embeddings(tokens, _get_model_name(language))
    return find_document_words(doc_embs, word_embs, word2doc, doc2word)

def enrichment_keywords(
    texts: Optional[List[str]] = None,
    tokens: Optional[List[List[str]]] = None,
    background_texts: Optional[List[str]] = None,
    background_tokens: Optional[List[List[str]]] = None,
    language: str = "slovenian",
):
    assert bool(texts) != bool(
        tokens
    ), "Parametri naj vsebujejo zgolj besedilo ali zgolj pojavnice"
    if not tokens:
        tokens = _preprocess_corpus(texts, language)

    assert (
        bool(background_texts) != bool(background_tokens),
        "Parametri naj vsebujejo zgolj besedilo (background text) ali zgolj "
        "pojavnice (background tokens)",
    )
    if not background_tokens:
        background_tokens = _preprocess_corpus(texts, language)

    temp_tokes = tokens + background_tokens
    joined_texts = [" ".join(tokens) for tokens in temp_tokes]
    vectorizer = TfidfVectorizer()
    X = vectorizer.fit_transform(joined_texts)
    words = vectorizer.get_feature_names()

    keywords = []
    for row in X[: len(tokens)]:
        p_values = hypergeom_p_values(X[len(tokens) :], row)
        order = np.argsort(p_values)
        keywords.append([(words[i], p_values[i]) for i in order])

    return keywords

def tfidf_keywords(
    texts: Optional[List[str]] = None,
    tokens: Optional[List[List[str]]] = None,
    language: str = "slovenian",
):
    assert bool(texts) != bool(
        tokens
    ), "Parametri naj vsebujejo zgolj besedilo ali zgolj pojavnice"
    if not tokens:
        tokens = _preprocess_corpus(texts, language)

    joined_texts = [" ".join(tokens) for tokens in tokens]
    vectorizer = TfidfVectorizer()
    X = vectorizer.fit_transform(joined_texts)
    words = vectorizer.get_feature_names()

    keywords = []
    for row in X:
        features = [(words[i], row[0, i]) for i in row.nonzero()[1]]
        keywords.append(sorted(features, key=lambda tup: tup[1], reverse=True))
    return keywords

def yake_keywords(
    texts: List[str],
    language: str = "slovenian",
    max_len: int = 1,
    stopwords: List[str] = None,
):
    # yake uses lancodes instead of full language name
    lg = langcodes.find(language).language
    yake_extractor = yake.KeywordExtractor(
        lan=lg, n=max_len, stopwords=stopwords
    )
    return [yake_extractor.extract_keywords(txt) for txt in texts]

def text_rank_keywords(
    texts: Optional[List[str]] = None,
    tokens: Optional[List[List[str]]] = None,
    language: str = "slovenian",
    num_words: int = 20,
    keyphrases: bool = False,
):
    assert bool(texts) != bool(
        tokens
    ), "Parametri naj vsebujejo zgolj besedilo ali zgolj pojavnice"
    if not tokens:
        tokens = _preprocess_corpus(texts, language)

    def text_rank(tokens):
        kw = textrank_kw(
            " ".join(tokens), words=num_words, scores=True, deacc=False
        )
        if not keyphrases:
            kw = [(x, sc) for w, sc in kw for x in w.split()]
        return kw

    return [text_rank(tokens) for tokens in tokens]

max_stars_repo_name: PrimozGodec/text-semantics

max_issues_repo_name: PrimozGodec/text-semantics

max_forks_repo_name: PrimozGodec/text-semantics


text:
#!/usr/bin/env python
"""  A very simple examle script to download, open, and plot an image
     from the Illustris Image host site.

     Author:  P. Torrey (ptorrey@mit.edu) 1/24/15

     Update:  P. Torrey (ptorrey@mit.edu) 6/21/15
		w/ suggestions from Geferson Lucatelli and Connor Bottrell
"""

import numpy as np
import os
# The two lines below make it possible for people to sit in the examples directory and run the scripts.
import sys
sys.path.append('../..')
import sunpy.sunpy__load as sunpy__load		#
import sunpy.sunpy__plot as sunpy__plot

my_api = "enter your own api key here !!!"
dl_base='http://www.illustris-project.org'

try:
    # Depending on the NumPy version, the line below can raise a TypeError even if the file is
    # present.  The try/except block enables us to try loadtxt calls that work for more NumPy
    # versions.
    try:
        catalog = np.loadtxt('directory_catalog_135.txt',
                             dtype={'names'  : ('subdirs', 'galaxy_numbers', 'galaxy_masses'),
                                    'formats': ('S3', 'i10', 'f8')})
    except TypeError:
        catalog = np.loadtxt('directory_catalog_135.txt',
                             dtype={'names'  : ('subdirs', 'galaxy_numbers', 'galaxy_masses'),
                                    'formats': ('S3', 'i8', 'f8')})
except IOError:
    # If we don't already have the file, we have to download it and try again.
    os.system("wget "+dl_base+"/files/directory_catalog_135.txt")
    try:
        catalog = np.loadtxt('directory_catalog_135.txt',
                             dtype={'names'  : ('subdirs', 'galaxy_numbers', 'galaxy_masses'),
                                    'formats': ('S3', 'i10','f8')})
    except TypeError:
        catalog = np.loadtxt('directory_catalog_135.txt',
                             dtype={'names'  : ('subdirs', 'galaxy_numbers', 'galaxy_masses'),
                                    'formats': ('S3', 'i8','f8')})

all_subdirs = catalog['subdirs']
all_galnrs  = catalog['galaxy_numbers']

for index,galnr in enumerate(all_galnrs[:1]):
    cmd = 'wget --content-disposition --header="API-Key: '+my_api+'" "'+dl_base+ \
	'/api/Illustris-1/snapshots/135/subhalos/'+str(galnr)+  \
	'/stellar_mocks/broadband.fits"'

    if( not (os.path.isfile("./broadband_"+str(galnr)+".fits")) ):
	os.system(cmd)
    filename = "./broadband_"+str(galnr)+".fits"

    # retrieve the image.  Could be used for non-parametric fitting, plotting, etc.
    image = sunpy__load.load_broadband_image(filename,band='B_Johnson.res')

    # plot an SDSS gri band idealized image
    sunpy__plot.plot_sdss_gri(filename, savefile='./sdss_gri_'+str(galnr)+'.png')

    os.remove("./"+filename)

max_stars_repo_name: rmandelb/sunpy

max_issues_repo_name: rmandelb/sunpy

max_forks_repo_name: rmandelb/sunpy


text:
#################################################################################
# The Institute for the Design of Advanced Energy Systems Integrated Platform
# Framework (IDAES IP) was produced under the DOE Institute for the
# Design of Advanced Energy Systems (IDAES), and is copyright (c) 2018-2021
# by the software owners: The Regents of the University of California, through
# Lawrence Berkeley National Laboratory,  National Technology & Engineering
# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia University
# Research Corporation, et al.  All rights reserved.
#
# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and
# license information.
#################################################################################

import os
import sys
import shutil
import enum
import copy
import json
import gzip
import numpy as np

import idaes
import pyomo.environ as pyo
from pyomo.common.collections import ComponentSet, ComponentMap
from pyomo.core.expr.visitor import identify_variables
import pyomo.dae as pyodae
from pyomo.common import Executable
from pyomo.dae.flatten import flatten_dae_components
from pyomo.util.subsystems import (
    TemporarySubsystemManager,
    create_subsystem_block,
)
from pyomo.solvers.plugins.solvers.ASL import ASL
from pyomo.opt.solver import SystemCallSolver
from pyomo.common.tempfiles import TempfileManager
from pyomo.common.errors import ApplicationError
from pyomo.util.calc_var_value import calculate_variable_from_constraint
from idaes.core.util.model_statistics import degrees_of_freedom
import idaes.logger as idaeslog
from idaes.core.solvers import get_solver
import idaes.config as icfg

PetscBinaryIOTrajectory = None
PetscBinaryIO = None

def _import_petsc_binary_io():
    global PetscBinaryIOTrajectory
    global PetscBinaryIO
    try:
        import PetscBinaryIOTrajectory
        import PetscBinaryIO
    except ImportError:
        petsc_dir = os.path.join(icfg.bin_directory, "petscpy")
        if not os.path.isdir(petsc_dir):
            return
        sys.path.append(petsc_dir)
        try:
            import PetscBinaryIOTrajectory
            import PetscBinaryIO
        except ImportError:
            pass

_import_petsc_binary_io()

class DaeVarTypes(enum.IntEnum):
    """Enum DAE variable type for suffix"""

    ALGEBRAIC = 0
    DIFFERENTIAL = 1
    DERIVATIVE = 2
    TIME = 3

@pyo.SolverFactory.register("petsc", doc="ASL PETSc interface")
class Petsc(ASL):
    """ASL solver plugin for the PETSc solver.  This adds the option to use an
    alternative executable batch file to run the solver through the WSL."""

    def __init__(self, **kwds):
        self._wsl = kwds.pop("wsl", None)
        super().__init__(**kwds)
        self.options.solver = "petsc"

    def _default_executable(self):
        """In addition to looking for the petsc executable, optionally check for
        a WSL batch file on Windows. Users could potentially also compile a
        cygwin exectable on Windows, so WSL isn't the only option, but it is the
        easiest for Windows."""
        executable = False
        if not self._wsl or self._wsl is None:
            executable = Executable("petsc")
        if sys.platform.startswith("win32") and (not executable):
            # On Windows, if wsl was requested or a normal petsc solver
            # executable was not found, look for a batch file to run it with WSL
            executable = Executable("petsc_wsl.bat")
        if not executable:
            raise RuntimeError("No PETSc executable found.")
        return executable.path()

@pyo.SolverFactory.register("petsc_snes", doc="ASL PETSc SNES interface")
class PetscSNES(Petsc):
    """PETSc solver plugin that sets options for SNES solver.  This turns on
    SNES monitoring, and checks the config for default options"""

    def __init__(self, **kwds):
        if "options" in kwds and kwds["options"] is not None:
            kwds["options"] = copy.deepcopy(kwds["options"])
        else:
            kwds["options"] = {}
        kwds["options"]["--snes_monitor"] = ""
        if "petsc_snes" in idaes.cfg:
            if "options" in idaes.cfg["petsc_snes"]:
                default_options = dict(idaes.cfg["petsc_snes"]["options"])
                default_options.update(kwds["options"])
                kwds["options"] = default_options
        super().__init__(**kwds)
        self.options.solver = "petsc_snes"

@pyo.SolverFactory.register("petsc_ts", doc="ASL PETSc TS interface")
class PetscTS(Petsc):
    """PETSc solver plugin that sets options for TS solver.  This turns on
    TS monitoring, sets the DAE flag, and checks the config for default options.
    """

    def __init__(self, **kwds):
        self._ts_vars_stub = kwds.pop("vars_stub", None)
        if "options" in kwds and kwds["options"] is not None:
            kwds["options"] = copy.deepcopy(kwds["options"])
        else:
            kwds["options"] = {}
        kwds["options"]["--dae_solve"] = ""
        kwds["options"]["--ts_monitor"] = ""
        if "petsc_ts" in idaes.cfg:
            if "options" in idaes.cfg["petsc_ts"]:
                default_options = dict(idaes.cfg["petsc_ts"]["options"])
                default_options.update(kwds["options"])
                kwds["options"] = default_options
        super().__init__(**kwds)
        self.options.solver = "petsc_ts"

    def _postsolve(self):
        stub = os.path.splitext(self._soln_file)[0]
        # There is a type file created by the solver to give the variable types
        # this is needed to read the trajectory data, and we want to include it
        # with other tmp files
        typ_file = stub + ".typ"
        TempfileManager.add_tempfile(typ_file)
        # If the vars_stub option was specified, copy the col and typ files to
        # the working directory. These files are needed to get the names and
        # types of variables and to make sense of trajectory data.
        if self._ts_vars_stub is not None:
            try:
                shutil.copyfile(f"{stub}.col", f"{self._ts_vars_stub}.col")
            except:
                pass
            try:
                shutil.copyfile(f"{stub}.typ", f"{self._ts_vars_stub}.typ")
            except:
                pass
        return ASL._postsolve(self)

@pyo.SolverFactory.register("petsc_tao", doc="ASL PETSc TAO interface")
class PetscTAO(Petsc):
    """This is a place holder for optimization solvers"""

    def __init__(self, **kwds):
        raise NotImplementedError(
            "The PETSc TAO interface has not yet been implemented"
        )

def petsc_available(wsl=None):
    """Check if the IDAES AMPL solver wrapper for PETSc is available.

    Args:
        wsl (bool): If True force WSL version, if False force not WSL version,
            if None, try non-WSL version then try WSL version

    Returns (bool):
        True if PETSc is available
    """
    solver = pyo.SolverFactory("petsc", wsl=wsl)
    if solver is not None:
        try:
            return solver.available()
        except RuntimeError:
            return False
    return False

def _copy_time(time_vars, t_from, t_to):
    """PRIVATE FUNCTION:

    This is used on the flattened (only indexed by time) variable
    representations to copy variable values that are unfixed at the "to" time
    from the value at the "from" time. The PETSc DAE solver uses the initial
    variable values as the initial condition, so this is used to copy the
    previous time in as the initial condition for the next step.

    Args:
        time_vars (list): list of variables or references to variables indexed
            only by time
        t_from (float): time point to copy from
        t_to (float): time point to copy to, only unfixed vars will be
            overwritten

    Returns:
        None
    """
    for v in time_vars:
        if not v[t_to].fixed:
            v[t_to].value = v[t_from].value

def find_discretization_equations(m, time):
    """This returns a list of time discretization equations. Since we aren't
    solving the whole time period simultaneously, we'll want to deactivate
    these constraints.

    Args:
        m (Block): model or block to search for constraints
        time (ContinuousSet):

    Returns:
        list of time discretization constraints
    """
    disc_eqns = []
    for var in m.component_objects(pyo.Var):
        if isinstance(var, pyodae.DerivativeVar):
            if time in ComponentSet(var.get_continuousset_list()):
                parent = var.parent_block()
                name = var.local_name + "_disc_eq"
                disc_eq = getattr(parent, name)
                disc_eqns.append(disc_eq)
    return disc_eqns

def _set_dae_suffixes_from_variables(m, variables, deriv_diff_map):
    """Write suffixes used by the solver to identify different variable types
    and associated derivative and differential variables.

    Args:
        m: model to search for variables and write suffixes to
        variables (list): List of time indexed variables at a specific time
            point
        deriv_diff_map (ComponentMap): Maps DerivativeVar data objects to
            differential variable data objects

    Returns:
        None
    """
    # The dae_suffix provides the solver information about variables types
    # algebraic, differential, derivative, or time, see DaeVarTypes
    m.dae_suffix = pyo.Suffix(
        direction=pyo.Suffix.EXPORT,
        datatype=pyo.Suffix.INT,
    )
    # The dae_link suffix provides the solver a link between the differential
    # and derivative variable, by assigning the pair a unique integer index.
    m.dae_link = pyo.Suffix(
        direction=pyo.Suffix.EXPORT,
        datatype=pyo.Suffix.INT,
    )
    dae_var_link_index = 1
    differential_vars = []
    i = 0
    for var in variables:
        if var in deriv_diff_map:
            deriv = var
            diffvar = deriv_diff_map[deriv]
            m.dae_suffix[diffvar] = int(DaeVarTypes.DIFFERENTIAL)
            m.dae_suffix[deriv] = int(DaeVarTypes.DERIVATIVE)
            m.dae_link[diffvar] = dae_var_link_index
            m.dae_link[deriv] = dae_var_link_index
            i += 1
            dae_var_link_index += 1
            if not diffvar.fixed:
                differential_vars.append(diffvar)
            else:
                raise RuntimeError(
                    f"Problem cannot contain a fixed differential variable and "
                    f"unfixed derivative. Consider either fixing the "
                    f"corresponding derivative or adding a constraint for the "
                    f"differential variable {diffvar} possibly using an "
                    f"explicit time variable."
                )
    return differential_vars

def _get_derivative_differential_data_map(m, time):
    """Get a map from data objects of derivative variables to the
    corresponding data objects of differential variables.

    Args:
        m: Model to search for DerivativeVars
        time: Set with respect to which DerivativeVars must be differentiated

    Returns:
        (ComponentMap): Map from derivative data objects to differential
            data objects
    """
    # Get corresponding derivative and differential data objects,
    # with no attention paid to fixed or active status.
    deriv_diff_list = []
    for var in m.component_objects(pyo.Var):
        if isinstance(var, pyodae.DerivativeVar) and time in ComponentSet(
            var.get_continuousset_list()
        ):
            deriv = var
            diffvar = deriv.get_state_var()
            for idx in var:
                if deriv[idx].fixed and pyo.value(abs(deriv[idx])) > 1e-10:
                    raise RuntimeError(
                        f"{deriv[idx]} is fixed to a nonzero value "
                        f"{pyo.value(deriv[idx])}. This is "
                        f"most likely a modeling error. Instead of fixing the "
                        f"derivative consider adding a constraint like "
                        f"dxdt = constant"
                    )
                deriv_diff_list.append((deriv[idx], diffvar[idx]))

    # Get unfixed variables in active constraints
    active_con_vars = ComponentSet()
    for con in m.component_data_objects(pyo.Constraint, active=True):
        for var in identify_variables(con.expr, include_fixed=False):
            active_con_vars.add(var)

    # Filter out derivatives that are fixed or not in an active constraint
    filtered_deriv_diff_list = []
    for deriv, diff in deriv_diff_list:
        if deriv in active_con_vars:
            filtered_deriv_diff_list.append((deriv, diff))
    return ComponentMap(filtered_deriv_diff_list)

def _sub_problem_scaling_suffix(m, t_block):
    """Copy scaling factors from the full model to the submodel.  This assumes
    the scaling suffixes could be in two places.  First check the parent block
    of the component (typical place for idaes models) then check the top-level
    model.  The top level model will take precedence.
    """
    if not hasattr(t_block, "scaling_factor"):
        # if the subsystem block doesn't already have a scaling suffix, make one
        t_block.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT)
    # first check the parent block for scaling factors
    for c in t_block.component_data_objects((pyo.Var, pyo.Constraint)):
        if hasattr(c.parent_block(), "scaling_factor"):
            if c in c.parent_block().scaling_factor:
                t_block.scaling_factor[c] = c.parent_block().scaling_factor[c]
    # now check the top level model
    if hasattr(m, "scaling_factor"):
        for c in t_block.component_data_objects((pyo.Var, pyo.Constraint)):
            if c in m.scaling_factor:
                t_block.scaling_factor[c] = m.scaling_factor[c]

def petsc_dae_by_time_element(
    m,
    time,
    timevar=None,
    initial_constraints=None,
    initial_variables=None,
    detect_initial=True,
    skip_initial=False,
    snes_options=None,
    ts_options=None,
    wsl=None,
    keepfiles=False,
    symbolic_solver_labels=True,
    vars_stub=None,
    trajectory_save_prefix=None,
):
    """Solve a DAE problem step by step using the PETSc DAE solver.  This
    integrates from one time point to the next.

    Args:
        m (Block): Pyomo model to solve
        time (ContinuousSet): Time set
        timevar (Var): Optional specification of a time variable, which can be
            used to write constraints that are an explicit function of time.
        initial_constraints (list): Constraints to solve in the initial
            condition solve step.  Since the time-indexed constraints are picked
            up automatically, this generally includes non-time-indexed
            constraints.
        initial_variables (list): This is a list of variables to fix after the
            initial condition solve step.  If these variables were originally
            unfixed, they will be unfixed at the end of the solve. This usually
            includes non-time-indexed variables that are calculated along with
            the initial conditions.
        detect_initial (bool): If True, add non-time-indexed variables and
            constraints to initial_variables and initial_constraints.
        skip_initial (bool): Don't do the initial condition calculation step,
            and assume that the initial condition values have already been
            calculated. This can be useful, for example, if you read initial
            conditions from a separately solved steady state problem, or
            otherwise know the initial conditions.
        snes_options (dict): PETSc nonlinear equation solver options
        ts_options (dict): PETSc time-stepping solver options
        wsl (bool): if True use WSL to run PETSc, if False don't use WSL to run
            PETSc, if None automatic. The WSL is only for Windows.
        keepfiles (bool): pass to keepfiles arg for solvers
        symbolic_solver_labels (bool): pass to symbolic_solver_labels argument
            for solvers. If you want to read trajectory data from the
            time-stepping solver, this should be True.
        vars_stub (str or None): Copy the `*.col` and `*.typ` files to the
            working directory using this stub if not None.  These are needed to
            interpret the trajectory data.
        trajectory_save_prefix (str or None): If a string is provided the
            trajectory data will be saved as gzipped json

    Returns:
        List of solver results objects from each solve. If there are initial
        condition constraints and they are not skipped, the first object will
        be from the initial condition solve.  Then there should be one for each
        time element for each TS solve.
    """
    solve_log = idaeslog.getSolveLogger("petsc-dae")
    regular_vars, time_vars = flatten_dae_components(m, time, pyo.Var)
    regular_cons, time_cons = flatten_dae_components(m, time, pyo.Constraint)
    tdisc = find_discretization_equations(m, time)

    solver_snes = pyo.SolverFactory("petsc_snes", options=snes_options, wsl=wsl)
    solver_dae = pyo.SolverFactory(
        "petsc_ts", options=ts_options, wsl=wsl, vars_stub=vars_stub
    )

    if initial_variables is None:
        initial_variables = []
    if initial_constraints is None:
        initial_constraints = []

    if detect_initial:
        rvset = ComponentSet(regular_vars)
        rcset = ComponentSet(regular_cons)
        icset = ComponentSet(initial_constraints)
        ivset = ComponentSet(initial_variables)
        initial_variables = list(ivset | rvset)
        initial_constraints = list(icset | rcset)

    # First calculate the inital conditions and non-time-indexed constraints
    res_list = []
    t0 = time.first()
    if not skip_initial:
        with TemporarySubsystemManager(to_deactivate=tdisc):
            constraints = [
                con[t0] for con in time_cons if t0 in con
            ] + initial_constraints
            variables = [var[t0] for var in time_vars] + initial_variables
            if len(constraints) > 0:
                # if the initial condition is specified and there are no
                # initial constraints, don't try to solve.
                t_block = create_subsystem_block(
                    constraints,
                    variables,
                )
                # set up the scaling factor suffix
                _sub_problem_scaling_suffix(m, t_block)
                with idaeslog.solver_log(solve_log, idaeslog.INFO) as slc:
                    res = solver_snes.solve(t_block, tee=slc.tee)
                res_list.append(res)

    tprev = t0
    count = 1
    fix_derivs = []
    with TemporarySubsystemManager(
        to_deactivate=tdisc,
        to_fix=initial_variables + fix_derivs,
    ):
        # Solver time steps
        deriv_diff_map = _get_derivative_differential_data_map(m, time)
        for t in time:
            if t == time.first():
                # t == time.first() was handled above
                continue
            constraints = [con[t] for con in time_cons if t in con]
            variables = [var[t] for var in time_vars]
            # Create a temporary block with references to original constraints
            # and variables so we can integrate this "subsystem" without
            # altering the rest of the model.
            t_block = create_subsystem_block(constraints, variables)
            differential_vars = _set_dae_suffixes_from_variables(
                t_block,
                variables,
                deriv_diff_map,
            )
            # We need to check if there are derivatives in the problem before
            # sending this to the solver.  We'll assume that if you are using
            # this and don't have any differential equations, you are making a
            # mistake.
            if len(differential_vars) < 1:
                raise RuntimeError(
                    "No differential equations found at t = %s, "
                    "you do not need a DAE solver." % t
                )
            if timevar is not None:
                t_block.dae_suffix[timevar[t]] = int(DaeVarTypes.TIME)
            # Set up the scaling factor suffix
            _sub_problem_scaling_suffix(m, t_block)
            # Take initial conditions for this step from the result of previous
            _copy_time(time_vars, tprev, t)
            with idaeslog.solver_log(solve_log, idaeslog.INFO) as slc:
                res = solver_dae.solve(
                    t_block,
                    tee=slc.tee,
                    keepfiles=keepfiles,
                    symbolic_solver_labels=symbolic_solver_labels,
                    export_nonlinear_variables=differential_vars,
                    options={"--ts_init_time": tprev, "--ts_max_time": t},
                )
            if trajectory_save_prefix is not None:
                tj = PetscTrajectory(
                    stub=vars_stub, delete_on_read=True, unscale=t_block
                )
                tj.to_json(f"{trajectory_save_prefix}_{count}.json.gz")
            tprev = t
            count += 1
            res_list.append(res)
    return res_list

def calculate_time_derivatives(m, time):
    """Calculate the derivative values from the discretization equations.

    Args:
        m (Block): Pyomo model block
        time (ContinuousSet): Time set

    Returns:
        None
    """
    for var in m.component_objects(pyo.Var):
        if isinstance(var, pyodae.DerivativeVar):
            if time in ComponentSet(var.get_continuousset_list()):
                parent = var.parent_block()
                name = var.local_name + "_disc_eq"
                disc_eq = getattr(parent, name)
                for i, v in var.items():
                    try:
                        if disc_eq[i].active:
                            v.value = 0 # Make sure there is a value
                            calculate_variable_from_constraint(v, disc_eq[i])
                    except KeyError:
                        pass  # discretization equation may not exist at first time

class PetscTrajectory(object):
    def __init__(
        self,
        stub=None,
        vecs=None,
        json=None,
        pth=None,
        vis_dir="Visualization-data",
        delete_on_read=False,
        unscale=None,
    ):
        """Class to read PETSc TS solver trajectory data.  This can either read
        PETSc output by providing the ``stub`` argument, a trajectory dict by
        providing ``vecs`` or a json file by providing ``json``.

        Args:
            stub (str): file name stub for variable info
            pth (str): path where variable info and trajectory data are stored
                if None, use current working directory
            vis_dir (str): subdirectory where visualization data is stored
            delete_on_read (bool): if true delete trajectory data after reading
            unscale (Block): if not None this is a block to read scale factors
                to unscale the trajectory
        """
        if PetscBinaryIOTrajectory is None:
            raise RuntimeError("PetscBinaryIOTrajectory could not be imported")
        if pth is not None:
            stub = os.path.join(pth, stub)
            vis_dir = os.path.join(pth, vis_dir)
        if stub is not None:
            self.stub = stub
            self.vis_dir = vis_dir
            self.path = pth
            self.unscale = unscale
            self._read()
            if delete_on_read:
                self.delete_files()
            if unscale is not None:
                self._unscale(unscale)
        elif vecs is not None:
            self.vecs = vecs
            self.time = vecs["_time"]
            self.vars = list(vecs.keys())
        elif json is not None:
            self.from_json(json)
        else:
            raise RuntimeError("To read trajectory, either provide stub, vecs, or json")

    def _read(self):
        with open(f"{self.stub}.col") as f:
            names = list(map(str.strip, f.readlines()))
        with open(f"{self.stub}.typ") as f:
            typ = list(map(int, f.readlines()))
        self.vars = [name for i, name in enumerate(names) if typ[i] in [0, 1]]
        (t, v, names) = PetscBinaryIOTrajectory.ReadTrajectory("Visualization-data")
        self.time = t
        self.vecs_by_time = v
        self.vecs = dict.fromkeys(self.vars, None)
        for k in self.vecs.keys():
            self.vecs[k] = [0] * len(self.time)
        self.vecs["_time"] = list(self.time)
        for i, v in enumerate(self.vars):
            for j, vt in enumerate(self.vecs_by_time):
                self.vecs[v][j] = vt[i]

    def _var_name_list(self, vars):
        if isinstance(vars, str):
            vars = [vars]
        if hasattr(vars, "ctype"):
            if issubclass(vars.ctype, pyo.Var):
                vars = [vars]
        vars = list(map(str, vars))
        for name in vars:
            if name not in self.vars:
                raise RuntimeError(f"Variable {name} not found.")
        return vars

    def get_vec(self, var):
        """Return the vector of variable values at each time point for var.

        Args:
            var (str or Var): Variable to get vector for.

        Retruns (list):
            vector of variable values at each time point

        """
        var = str(var)
        return self.vecs[var]

    def get_dt(self):
        """Get a list of time steps

        Args:
            None

        Returns:
            (list)
        """
        dt = [None] * (len(self.time) - 1)
        for i in range(len(self.time) - 1):
            dt[i] = self.time[i + 1] - self.time[i]
        return dt

    def interpolate_vecs(self, times):
        """Create a new vector dictionary interpolated at times. This method
        will also extraplote values outside the original time range, so care
        should be taken not to specify times too far outside the range.

        Args:
            times (list): list of times to interpolate. These must be in
                increasing order.

        Returns (dict):
            Dictionary of vectors for values at interpolated points
        """
        vecs = dict.fromkeys(self.vars, None)
        vecs["_time"] = copy.copy(times)
        for var in vecs:
            vecs[var] = np.interp(vecs["_time"], self.vecs["_time"], self.vecs[var])
        return vecs

    def _unscale(self, m):
        """If variable scale factors are used, the solver will see scaled
        variables, and the scaled trajectory will be written. This function
        uses variable scaling facors from the given model to unscale the
        trajectory.

        Args:
            m (Block): model or block to read scale factors from.

        Returns:
            None
        """
        # Variables might show up more than once because of References
        already_scaled = set()
        for var in m.component_data_objects():
            vname = str(var)
            if vname in self.vecs and vname not in already_scaled:
                s = None
                if hasattr(var.parent_block(), "scaling_factor"):
                    s = var.parent_block().scaling_factor.get(var, s)
                if hasattr(m, "scaling_factor"):
                    s = m.scaling_factor.get(var, s)
                if s is not None:
                    for i, x in enumerate(self.vecs[vname]):
                        self.vecs[vname][i] = x / s
                already_scaled.add(vname)

    def delete_files(self):
        """Delete the trajectory data and variable information files.

        Args:
            None

        Returns:
            None
        """
        shutil.rmtree(self.vis_dir)
        os.remove(f"{self.stub}.col")
        os.remove(f"{self.stub}.typ")

    def to_json(self, pth):
        """Dump the trajectory data to a json file in the form of a dictionary
        with variable name keys and '_time' with vectors of values at each time.

        Args:
            pth (str): path for json file to write

        Returns:
            None
        """
        if pth.endswith(".gz"):
            with gzip.open(pth, "w") as fp:
                fp.write(json.dumps(self.vecs).encode("utf-8"))
        else:
            with open(pth, "w") as fp:
                json.dump(self.vecs, fp)

    def from_json(self, pth):
        """Read the trajectory data from a json file in the form of a dictionary.

        Args:
            pth (str): path for json file to write

        Returns:
            None
        """
        if pth.endswith(".gz"):
            with gzip.open(pth, "r") as fp:
                self.vecs = json.loads(fp.read())
        else:
            with open(pth, "r") as fp:
                self.vecs = json.load(fp)
        self.time = self.vecs["_time"]
        self.vars = list(self.vecs.keys())

max_stars_repo_name: CATER-UCF/idaes-pse

max_issues_repo_name: CATER-UCF/idaes-pse

max_forks_repo_name: CATER-UCF/idaes-pse

reason: import numpy,import pyomo,from pyomo


text:
import fileinput

import numpy as np
import matplotlib.pyplot as plt

selected_dim = 1
p_crit = 0.644701

beta = 0.276486
v_para = 1.733847
v_perp = 1.096854
beta_p = beta
z = v_para/v_perp
alpha = beta_p/v_para

t_max = 1024

sizes = {}
for line in fileinput.input():
    pc, dims, seed, size, mass, lt, volume = \
            [float(i.strip()) for i in line.split()]

    if dims == selected_dim:
        nodes = int(size**dims)
        if nodes not in sizes:
            sizes[nodes] = {}
        if pc not in sizes[nodes]:
            sizes[nodes][pc] = {"mass": [], "lt": [], "volume": []}

        sizes[nodes][pc]["mass"].append(int(mass))
        sizes[nodes][pc]["volume"].append(int(volume))
        sizes[nodes][pc]["lt"].append(int(lt))

survival_fig, survival_ax = plt.subplots()
scaled_survival_fig, scaled_survival_ax = plt.subplots()

for s in sizes:
    ts = np.arange(1, t_max+1)
    lts = np.array(sizes[s][p_crit]["lt"])
    p_surv = np.array([np.count_nonzero(lts >= t) for t in ts])/lts.shape[0]

    survival_ax.plot(ts, p_surv, label=f"N={s} p={p_crit}")
    scaled_survival_ax.plot(
            ts/s**(z/(selected_dim)),
            (ts**alpha)*p_surv,
            label=f"N={s} p={p_crit}")

survival_ax.plot(ts, ts**-alpha, label="$t^{-\\alpha}$", ls='--')
survival_ax.set_xlabel("$t$")
survival_ax.set_ylabel("$P(t)$")
survival_ax.set_xscale("log")
survival_ax.set_yscale("log")
survival_ax.legend()
survival_fig.savefig("figures/survival.pdf")

scaled_survival_ax.set_xlabel("$t/N^{z/d}$")
scaled_survival_ax.set_ylabel("$t^{\\alpha} P(t)$")
scaled_survival_ax.set_xscale("log")
scaled_survival_ax.set_yscale("log")
scaled_survival_ax.set_ylim(1e-4, 4)
scaled_survival_ax.legend()
scaled_survival_fig.savefig("figures/scaled-survival.pdf")

max_stars_repo_name: arashbm/lattice

max_issues_repo_name: arashbm/lattice

max_forks_repo_name: arashbm/lattice


text:
"""
MIT License

Copyright (c) 2019 ming

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
===============================
@File  : speedProfileOptimizer4.py
@Author: ming.ustb@outlook.com
@Date  : 19-6-12
@GitHub: https://github.com/yangmingustb/PTPSim
"""

import pprint
from cvxopt import matrix, solvers
import math
import numpy
import matplotlib.pyplot as plt

s0 = 200.0
t0 = 24.2
start = [s0, t0]
init_v = 55/3.6  # init_speed
init_acc = 0  # init_acceleration
init_jerk = 0  # init_jerk
init_s = s0
init_t = t0
lane_width = 3.75
refLineRho = lane_width*0.5
laneChaneRefLine = lane_width*1.5
static_obs = [ [s0 + 65, laneChaneRefLine - 0.3]]
start_SRho = [s0, laneChaneRefLine, 0.0 * math.pi / 180.0]
max_t = 8  # 8s
max_s = static_obs[0][0]-s0  # 100m
s_end = s0+max_s
print('==================')
print(s_end)
reso_t = 0.1
reso_s = 2
max_v = 20
max_acc = 2.5  # liu chang liu. 2017 IV. speed profile planning
ref_v = 40.0/3.6  # 前方没有障碍物时的参考速度
#   ref_v_static = 0  #   前方有静态障碍物时的速度
max_jerk = 5
n_t = int(max_t / reso_t)

r = 1  # 碰撞检测圆半径
d = 2  # 圆心间距
# obs = [[60, 0]]

obs_threshold = 2  # 静态障碍物膨胀阈值
safeDistance = r + obs_threshold
safeSpeed = 10/3.6      # 安全速度
overtakingSpeed = 60/3.6
stoppingSpeed= 0
flag_obs = 1  # 障碍物标志位，0前方无障碍物，1前方静态障碍物，2前方动态障碍物
rough_st = []
rough_v = []
# guidance_v = ref_v * numpy.ones((1, n_t))
# for i in range(n_t):
# 	# if i == int(55 / max_s * n_t):
# 	# 	for j in numpy.arange(i, 77, 1):
# 	# 		v=ref_v - j
# 	# 		if v <= stoppingSpeed:
# 	# 			v = stoppingSpeed
# 	# 		guidance_v[0][j] = v
# 	if i == int(75/max_s * n_t):
# 		print('i:',i)
# 		for j in numpy.arange(i, n_t, 1):
# 			guidance_v[0][j] = stoppingSpeed
#
# 	else:continue
# print("----------------")
# print(guidance_v)
w_vel = 1
# w_vel = 0.1
w_acc = 1
w_jerk = 10
showAnimation = 1

def cal_H_vel():
	H_vel = numpy.eye(n_t + 1)

	for i in range(n_t + 1):
		if i > 0 and i < n_t:
			H_vel[i][i] = 2
		if i > 0 and i < n_t + 1:
			H_vel[i - 1][i] = -1
			H_vel[i][i - 1] = -1

	#   print("H_vel size:", numpy.shape(H_vel), H_vel)

	return H_vel

def cal_q_vel(guidance_v):
	"""
	暂时使用默认参考速度进行优化，等调试成熟再使用粗糙速度来优化
	:return:
	"""
	q_vel = numpy.zeros((1, n_t + 1))

	if flag_obs == 0:
		q_vel[0][0] = -ref_v
		q_vel[0][n_t] = ref_v

	if flag_obs == 1:
		for i in range(n_t + 1):
			if i < 1:
				q_vel[0][i] = -guidance_v[0][i]
			elif i >= n_t:

				q_vel[0][i] = guidance_v[0][i - 1]

			else:
				q_vel[0][i] = guidance_v[0][i - 1] - guidance_v[0][i]
	#   print('q_vel:', numpy.shape(q_vel), q_vel)

	return q_vel

def cal_H_acc():
	H_acc = numpy.eye(n_t + 2)

	for i in range(n_t + 2):
		if i == 1:
			H_acc[i][i] = 5
			H_acc[i - 1][i] = - 2
			H_acc[i][i - 1] = - 2
		if i > 1 and i < n_t:
			H_acc[i][i] = 6
			H_acc[i - 1][i] = - 4
			H_acc[i - 2][i] = 1
			H_acc[i][i - 1] = - 4
			H_acc[i][i - 2] = 1

		if i == n_t:
			H_acc[i][i] = 5
			H_acc[i - 1][i] = - 4
			H_acc[i - 2][i] = 1
			H_acc[i][i - 1] = - 4
			H_acc[i][i - 2] = 1

		if i == n_t + 1:
			H_acc[i - 1][i] = - 2
			H_acc[i - 2][i] = 1
			H_acc[i][i - 1] = - 2
			H_acc[i][i - 2] = 1

	#   print("H_acc size:", numpy.shape(H_acc), "\n", H_acc)

	return H_acc

def cal_H_jerk():
	H_jerk = numpy.eye(n_t + 3)

	for i in range(n_t + 3):
		if i == 1:
			H_jerk[i][i] = 10
			H_jerk[i - 1][i] = - 3
			H_jerk[i][i - 1] = - 3
		if i == 2:
			H_jerk[i][i] = 19
			H_jerk[i - 1][i] = - 12
			H_jerk[i - 2][i] = 3
			H_jerk[i][i - 1] = - 12
			H_jerk[i][i - 2] = 3

		if i > 2 and i < n_t:
			H_jerk[i][i] = 20
			H_jerk[i - 1][i] = - 15
			H_jerk[i - 2][i] = 6
			H_jerk[i - 3][i] = - 1

			H_jerk[i][i - 1] = - 15
			H_jerk[i][i - 2] = 6
			H_jerk[i][i - 3] = - 1

		if i == n_t:
			H_jerk[i][i] = 19
			H_jerk[i - 1][i] = - 15
			H_jerk[i - 2][i] = 6
			H_jerk[i - 3][i] = - 1

			H_jerk[i][i - 1] = - 15
			H_jerk[i][i - 2] = 6
			H_jerk[i][i - 3] = - 1

		if i == n_t + 1:
			H_jerk[i][i] = 10
			H_jerk[i - 1][i] = - 12
			H_jerk[i - 2][i] = 6
			H_jerk[i - 3][i] = - 1
			H_jerk[i][i - 1] = - 12
			H_jerk[i][i - 2] = 6
			H_jerk[i][i - 3] = - 1

		if i == n_t + 2:
			H_jerk[i - 1][i] = - 3
			H_jerk[i - 2][i] = 3
			H_jerk[i - 3][i] = - 1
			H_jerk[i][i - 1] = - 3
			H_jerk[i][i - 2] = 3
			H_jerk[i][i - 3] = - 1

	#   print("H_jerk size:", numpy.shape(H_jerk), "\n", H_jerk)
	#   print(H_jerk[82])
	return H_jerk

def cal_H_total():
	H_vel = cal_H_vel()
	tmp_y = numpy.zeros((n_t + 1, 2))
	H_vel = numpy.c_[H_vel, tmp_y]
	tmp_x = numpy.zeros((2, n_t + 3))
	H_vel = numpy.r_[H_vel, tmp_x]

	H_acc = cal_H_acc()
	tmp_y = numpy.zeros((n_t + 2, 1))
	H_acc = numpy.c_[H_acc, tmp_y]
	tmp_x = numpy.zeros((1, n_t + 3))
	H_acc = numpy.r_[H_acc, tmp_x]

	H_jerk = cal_H_jerk()
	#   print(H_vel)
	#   print(H_acc)
	#   print(H_jerk)

	H = 2 * w_vel / reso_t * H_vel + 2 * w_acc / (reso_t ** 3) * H_acc + 2 * w_jerk / (reso_t ** 5) * H_jerk
	return H

def cal_q(guidance_v):
	q_vel = cal_q_vel(guidance_v)
	tmp_y = numpy.zeros((1, 2))
	tmp_q = numpy.c_[q_vel, tmp_y]
	q = -2 * w_vel * tmp_q

	#   print(q)
	return q

def constraints_acc():
	h_acc = reso_t ** 2 * max_acc * numpy.ones((n_t, 1))
	g_acc = numpy.zeros((n_t, n_t + 3))
	for i in range(n_t):
		g_acc[i][i] = 1
		g_acc[i][i + 1] = -2
		g_acc[i][i + 2] = 1
	return g_acc, h_acc

def constraints_jerk():
	h_jerk = reso_t ** 3 * max_jerk * numpy.ones((n_t, 1))
	g_jerk = numpy.zeros((n_t, n_t + 3))
	for i in range(n_t):
		g_jerk[i][i] = -1
		g_jerk[i][i + 1] = 3
		g_jerk[i][i + 2] = -3
		g_jerk[i][i + 3] = 1

	return g_jerk, h_jerk

def constraints_speed():
	h_vel = reso_t * max_v * numpy.ones((n_t, 1))
	g_vel = numpy.zeros((n_t, n_t + 3))
	for i in range(n_t):
		g_vel[i][i] = -1
		g_vel[i][i + 1] = 1

	return g_vel, h_vel

def constraints_monotonicity():
	"""
	单调性约束
	:return: 维度, (n_t+2)*(n_t+3)
	"""
	h_inc = numpy.zeros((n_t + 2, 1))
	g_inc = numpy.zeros((n_t + 2, n_t + 3))
	for i in range(n_t + 2):
		g_inc[i][i] = 1
		g_inc[i][i + 1] = -1

	return g_inc, h_inc

def constraints_range():
	h_max = s_end * numpy.ones((n_t, 1))
	h_init = init_s * numpy.ones((n_t, 1))
	g_range = numpy.zeros((n_t, n_t + 3))
	for i in range(n_t):
		g_range[i][i] = 1

	return g_range, h_max, h_init

def static_obs_constraints():
	h_obs1 = (static_obs[0][0] - safeDistance - 2 * d) * numpy.ones((n_t + 3, 1))
	h_obs2 = (-static_obs[0][0] - safeDistance) * numpy.ones((n_t + 3, 1))
	g_obs1 = numpy.eye(n_t + 3)
	g_obs2 = -numpy.eye(n_t + 3)

	return g_obs1, g_obs2, h_obs1, h_obs2

def constraints_equation():
	b = numpy.zeros((6, 1))
	b[0][0] = reso_t * init_v + init_s
	b[1][0] = init_acc * reso_t ** 2 + init_s + 2 * reso_t * init_v
	b[2][0] = init_jerk * reso_t ** 3 + 3 * reso_t * init_v + 3 * init_acc * reso_t ** 2 + init_s

	a = numpy.zeros((6, n_t + 3))
	a[0][0] = 1
	a[1][1] = 1
	a[2][2] = 1

	a[3][n_t - 2] = 1
	a[3][n_t - 1] = -2
	a[3][n_t] = 1

	a[4][n_t - 2] = 2
	a[4][n_t - 1] = -3
	a[4][n_t + 1] = 1

	a[5][n_t - 2] = 3
	a[5][n_t - 1] = -4
	a[5][n_t + 2] = 1

	return a, b

def all_constraints():
	g_acc, h_acc = constraints_acc()
	g_jerk, h_jerk = constraints_jerk()
	g_vel, h_vel = constraints_speed()
	g_inc, h_inc = constraints_monotonicity()
	g_max, h_max, h_init = constraints_range()

	g = numpy.r_[g_acc, -g_acc]
	g = numpy.r_[g, g_jerk]
	g = numpy.r_[g, -g_jerk]
	g = numpy.r_[g, g_vel]
	g = numpy.r_[g, g_inc]
	g = numpy.r_[g, g_max]
	g = numpy.r_[g, -g_max]

	h = numpy.r_[h_acc, h_acc]
	h = numpy.r_[h, h_jerk]
	h = numpy.r_[h, h_jerk]
	h = numpy.r_[h, h_vel]
	h = numpy.r_[h, h_inc]
	h = numpy.r_[h, h_max]
	h = numpy.r_[h, h_init]

	if len(static_obs):
		g_obs1, g_obs2, h_obs1, h_obs2 = static_obs_constraints()
		g = numpy.r_[g, g_obs1]
		# g = numpy.r_[g, g_obs2]
		h = numpy.r_[h, h_obs1]
		# h = numpy.r_[h, h_obs2]

	a, b = constraints_equation()
	return g, h, a, b

def QP_solver(guidance_v):
	H = cal_H_total()
	q = cal_q(guidance_v)
	q = q.T
	g, h, a, b = all_constraints()

	P = matrix(H)
	q = matrix(q)
	G = matrix(g)
	h = matrix(h)
	A = matrix(a)
	b = matrix(b)
	result = solvers.qp(P, q, G, h, A, b)
	s_set = result['x']

	return s_set

def plotTrajectory():
	s_set = QP_solver()
	s_set = numpy.array(s_set).T

	#   print(type(s_set))
	shape = numpy.shape(s_set)
	t_tmp = [init_t]
	s_tmp = [init_s]

	for i in range(shape[1]):
		t_tmp.append(init_t + reso_t * (i + 1))
		s_tmp.append(s_set[0][i])
	return t_tmp, s_tmp

def plot_figure():
	# guidance_v= EMspeedOptimizer()
	guidance_v = ref_v * numpy.ones((1, n_t))
	s_set = QP_solver(guidance_v)
	s_set = numpy.array(s_set).T

	#   print(type(s_set))
	shape = numpy.shape(s_set)
	t_tmp = [init_t]
	s_tmp = [init_s]

	for i in range(shape[1]):
		t_tmp.append(init_t + reso_t * (i + 1))
		s_tmp.append(s_set[0][i])

	print('s_tmp:', len(s_tmp))

	font1 = {'family': 'Times New Roman', 'weight': 'normal', 'size': 10}
	plt.rcParams['font.sans-serif'] = ['Times New Roman']  # 如果要显示中文字体，则在此处设为：SimHei

	# plt.style.use('ggplot')
	plt.figure(figsize=(3.5, 4.2))  # 单位英寸， 3.5

	p1 = [0.18, 0.85, 0.72, 0.1]
	p2 = [0.18, 0.60, 0.72, 0.08]
	p3 = [0.18, 0.35, 0.72, 0.05]
	p4 = [0.18, 0.1, 0.72, 0.05]
	ax1 = plt.axes(p1)  # 2行2列,第一个图
	ax2 = plt.axes(p2)  # 2行2列,第一个图
	ax3 = plt.axes(p3)  # 2行2列,第一个图
	ax4 = plt.axes(p4)  # 2行2列,第一个图

	ax1.plot(t_tmp, s_tmp, lw=1, label='s-t Graph', c='b', alpha=1)
	ax1.grid(linestyle="--", linewidth=0.5, alpha=1)
	ax1.set_title('s-t Graph', font1)
	ax1.set_xlabel('Time (s)', font1)
	ax1.set_ylabel('s (m)', font1)
	# plt.xticks(fontproperties='Times New Roman', fontsize=10)
	# plt.yticks(fontproperties='Times New Roman', fontsize=10)
	# ax1.set_xlim(-1, 9)
	# ax1.set_ylim(190, 280)
	#   plt.grid(True)
	#   plt.title('A simple plot')

	v_tmp = []
	for i in range(len(s_tmp) - 1):
		v = (s_tmp[i + 1] - s_tmp[i]) / reso_t
		v_tmp.append(v)

	ax2.plot(t_tmp[0: (len(t_tmp) - 1)], v_tmp, lw=1, label='Speed Profile', c='b', alpha=0.8)
	ax2.set_title('Speed Profile', font1)
	ax2.grid(linestyle="--", linewidth=0.5, alpha=1)
	# plt.legend(loc=0)  # 图例位置自动
	# ax2.axis('tight')
	ax2.set_xlabel('Time (s)', font1)
	ax2.set_ylabel('v (m/s)', font1)
	# ax2.set_xlim(-1, 9)
	# ax2.set_ylim(-1, 20)

	acc_tmp = []
	for i in range(len(s_tmp) - 2):
		acc = (s_tmp[i + 2] - 2 * s_tmp[i + 1] + s_tmp[i]) / reso_t ** 2
		acc_tmp.append(acc)

	ax3.plot(t_tmp[0: (len(t_tmp) - 2)], acc_tmp, lw=1, label='Acceleration Profile', alpha=0.8, c='b')
	ax3.set_title('Acceleration Profile', font1)
	ax3.grid(linestyle="--", linewidth=0.5, alpha=1)

	# ax3.legend(loc=0)  # 图例位置自动
	# ax3.axis('tight')
	ax3.set_xlabel('Time (s)', font1)
	ax3.set_ylabel(r'a (m/s$\mathrm{^2}$)', font1)
	# ax3.set_xlim(-1, 9)
	# ax3.set_ylim(-4, 4)

	jerk_tmp = []
	for i in range(len(s_tmp) - 3):
		jerk = (s_tmp[i + 3] - 3 * s_tmp[i + 2] + 3 * s_tmp[i + 1] - s_tmp[i]) / reso_t ** 3
		jerk_tmp.append(jerk)

	ax4.plot(t_tmp[0: (len(t_tmp) - 3)], jerk_tmp, lw=1, label='Jerk Profile', c='b', alpha=0.8)
	ax4.set_title('Jerk Profile', font1)
	ax4.grid(linestyle="--", linewidth=0.5, alpha=1)

	# plt.legend(loc=0)  # 图例位置自动
	# ax4.axis('tight')
	ax4.set_xlabel('Time (s)', font1)
	ax4.set_ylabel(r'jerk (m/s$\mathrm{^3}$)', font1)
	# ax4.set_xlim(-1, 9)
	# ax4.set_ylim(-10, 10)
	plt.savefig('../SimGraph/speedOptimizer4_061101.svg')
	plt.savefig('../SimGraph/speedOptimizer4_061101.tiff', dpi=600)
	#   plt.grid(True)
	plt.show()

def EMspeedOptimizer():
	"""
	EM speed optimization
	:return:
	"""
	guidance_v = ref_v * numpy.ones((1, n_t))
	for i in range(n_t):
		# if i == int(55 / max_s * n_t):
		# 	for j in numpy.arange(i, 77, 1):
		# 		v=ref_v - j
		# 		if v <= stoppingSpeed:
		# 			v = stoppingSpeed
		# 		guidance_v[0][j] = v
		if i == int(45/max_s * n_t):
			print('i:',i)
			for j in numpy.arange(i, n_t, 1):
				guidance_v[0][j] = stoppingSpeed

		else:continue
	# print("----------------")
	# print(guidance_v)

	s_set = QP_solver(guidance_v)
	s_set = numpy.array(s_set).T

	#   print(type(s_set))
	shape = numpy.shape(s_set)
	t_tmp = [init_t]
	s_tmp = [init_s]

	for i in range(shape[1]):
		t_tmp.append(init_t + reso_t * (i + 1))
		s_tmp.append(s_set[0][i])

	v_tmp = []
	for i in range(len(s_tmp) - 1):
		v = (s_tmp[i + 1] - s_tmp[i]) / reso_t
		v_tmp.append(v)

	speed = numpy.array([v_tmp[0:80]])
	print("speed:", speed)
	return speed

if __name__ == '__main__':
	plot_figure()

max_stars_repo_name: yangmingustb/PTPSim

max_issues_repo_name: yangmingustb/PTPSim

max_forks_repo_name: yangmingustb/PTPSim