text:
import json, sys
from os import close
import pandas as pd
import time
from datetime import datetime

import numpy as np
import math
from pandas.plotting import register_matplotlib_converters
import warnings

sys.path.append("../")
# sys.path.insert(1, 'users/grahamlim/Documents/Documents-iMac/my_algos/config_py_files')
from config_py_files import johansen as joh
# import johansen.coint_johansen

# df_object = json.loads(sys.argv[1])
df_object = pd.read_json(sys.argv[1])

# print('printing df_object...')
# print(df_object.head())
# print(df_object.columns)

# print(df_object.iloc[0])

# the now timestmap will also be taken from df_object later
actual_now = datetime.now()
now_readable = actual_now.strftime(("%d_%m_%Y_%H_%M_%S"))
now = str(actual_now)

front_leg_instrument = df_object.iloc[0]["leg"]
middle_leg_instrument = df_object.iloc[1]["leg"]
back_leg_instrument = df_object.iloc[2]["leg"]
tf = int(df_object.iloc[3]["tf"])

# print(f'tf is {tf}')
# print(f'front_leg_instrument = {front_leg_instrument}')

# tick_index = first_df["ticks"]
# close_price = first_df["close"]
# leg_name = first_df["leg"]

# temp_dict = {'ticks': tick_index, leg_name: close_price}

# converted_first_df = pd.DataFrame(temp_dict)
# converted_first_df.set_index('ticks', inplace= True)
# print(converted_first_df.head())
# print(f'length: {len(converted_first_df.index)}')

def json_df_converter(df_object, object_index):
    raw_df = df_object.iloc[object_index]
    tick_index = raw_df["ticks"]
    close_price = raw_df["close"]
    leg_name = raw_df["leg"]

    temp_dict = {'ticks': tick_index, leg_name: close_price}
    converted_df = pd.DataFrame(temp_dict)
    converted_df.set_index('ticks', inplace= True)
    return converted_df

def triple_df_builder(df_object):
    #get leg 1
    try:
        df_front_leg = json_df_converter(df_object, 0)

        df_middle_leg = json_df_converter(df_object, 1)

        df_back_leg = json_df_converter(df_object, 2)

        df = pd.concat([df_front_leg, df_middle_leg, df_back_leg], axis=1)

        #some contracts commence earlier than the others so we want to have the 3 instruments overlap
        #for the butterfly spread to work
        df.dropna(inplace = True)

        # print(f'length of dataframe is {len(df)}')
        return df
    except:
        print('one of the legs are not processing')
        pass

def johansen_test(df):
    # print('running Johansen test...')
    register_matplotlib_converters()

    # For checking cointegration
    sys.path.append("..")
    warnings.filterwarnings('ignore')

    # Run Johansen test on the whole dataset

    # Store the results of Johansen test after applying on the dataframe
    result = joh.coint_johansen(df, 0, 1)

    # Print trace statistics and eigen statistics
    # print ('--------------------------------------------------')
    # print ('--> Trace Statistics')
    # print ('variable statistic Crit-90% Crit-95%  Crit-99%')
    # for i in range(len(result.lr1)):
    #     print ("r <= " + str(i), round(result.lr1[i], 4), round(result.cvt[i, 0],4), round(result.cvt[i, 1],4), round(result.cvt[i, 2],4))
    # print ('--------------------------------------------------')
    # print ('--> Eigen Statistics')
    # print ('variable statistic Crit-90% Crit-95%  Crit-99%')
    # for i in range(len(result.lr2)):
    #     print ("r <= " + str(i), round(result.lr2[i], 4), round(result.cvm[i, 0],4), round(result.cvm[i, 1],4), round(result.cvm[i, 2],4))

    #Print half-life, or time taken for data to revert to mean
    theta = result.eig[0]
    half_life = math.log(2)/theta
    # print ('--------------------------------------------------')
    # print ('--> Half Life')
    # print (f"The expected holding period is {half_life} periods.")

    # Store the results of Johansen test
    half_df_length = len(df)/2
    lookback_johansen = int(round(half_df_length)) #set to half of whatever dataframe, always.
    result = joh.coint_johansen(df[:lookback_johansen], 0, 1)

    # Store the value of eigenvector. Using this eigenvector, you can create the spread
    ev = result.evec

    # Take the transpose and use the first row of eigenvectors as it forms strongest cointegrating spread
    ev = result.evec.T[0]

    # Normalise the eigenvectors by dividing the row with the first value of eigenvector
    ev = ev/ev[0]

    front_vector = round(ev[0], 2)
    middle_vector = round(ev[1],2)
    back_vector = round(ev[2],2)

    # Print the mean reverting spread
    # print(f"\nCurrent Cointegrated Spread for {tf}= {front_vector}*{front_leg_instrument} + ({middle_vector})*{middle_leg_instrument} + ({back_vector})*{back_leg_instrument}")
    # print(f"Johansen test performed on {df.index[-1]}")
    return [front_vector, middle_vector, back_vector]

def butterfly_sl_backtester(df, front_vector, middle_vector, back_vector, lookback, std_dev, sl_flag, sl_std_dev, one_way_fee):
    # print('running new butterfly sl backtester...')
    df['spread'] = df[front_leg_instrument]*front_vector + df[middle_leg_instrument]*middle_vector + df[back_leg_instrument]*back_vector

    # requires x and y plotting to be done in numpy arrays
    # print("Spread for Timeframe {0}, Lookback {1}, Std Deviation {2}".format(tf, lookback, std_dev))
    # spread_timestamps = list(df.index.values)
    # fig.plot(spread_timestamps, df["spread"], label="spread", width=50, height=15)
    # fig.show()

    # Moving Average and Moving Standard Deviation
    df['moving_average'] = df.spread.rolling(lookback).mean()
    df['moving_std_dev'] = df.spread.rolling(lookback).std()

    # Upper band and lower band
    df['upper_band'] = df.moving_average + std_dev*df.moving_std_dev
    df['lower_band'] = df.moving_average - std_dev*df.moving_std_dev

    # Long entries
    df['long_entry'] = df.spread < df.lower_band
    df['long_exit'] = df.spread >= df.moving_average

    # Short entries
    df['short_entry'] = df.spread > df.upper_band
    df['short_exit'] = df.spread <= df.moving_average

    #binary positions tracking
    df['positions_long'] = np.nan
    df.loc[df.long_entry, 'positions_long'] = 1
    df.loc[df.long_exit, 'positions_long'] = 0
    df['positions_short'] = np.nan
    df.loc[df.short_entry, 'positions_short'] = -1
    df.loc[df.short_exit, 'positions_short'] = 0

    # ------------------------- Stop Loss Logic -------------------------------
    if sl_flag:
        # Add a stop loss x standard deviation away from lower band.
        # For example, refer to the above markdown cell
        df['long_sl'] = df.spread < (df.lower_band - sl_std_dev*df.moving_std_dev)

        # Whenever the long_sl is True, set the positions_long_sl to 0
        df.loc[df.long_sl, 'positions_long_sl'] = 0

        # Whenever the price reverts to the mean, set the positions_long_sl to 1
        df.loc[df.long_exit, 'positions_long_sl'] = 1

        # Add a stop loss x standard deviation away from upper band.
        df['short_sl'] = df.spread > (df.upper_band + sl_std_dev*df.moving_std_dev)

        # Similar to above but for short side
        df.loc[df.short_sl, 'positions_short_sl'] = 0
        df.loc[df.short_exit, 'positions_short_sl'] = 1
        df = df.fillna(method='ffill')

        # Multiply with the positions column to include the impact of the stop loss
        df.positions_long = df.positions_long * df.positions_long_sl
        df.positions_short = df.positions_short * df.positions_short_sl

    # Fill NaN values
    df = df.fillna(method='ffill')

    # Consolidate the positions
    df['positions'] = df.positions_long + df.positions_short

    # Calculate spread difference
    df['spread_diff'] = df.spread-df.spread.shift(1)

    # Calculate portfolio value
    df['portfolio_value'] = abs(front_vector)*df[front_leg_instrument] + abs(middle_vector)*df[middle_leg_instrument] + abs(back_vector)*df[back_leg_instrument]

    # Calculate daily returns
    df['daily_returns'] = df.spread_diff / df['portfolio_value'].shift(1)

    # Calculate strategy returns
    df['strategy_returns'] = df.daily_returns  * df.positions.shift(1)

    # Calculate fresh entries
    # fresh longs & shorts, where previously positions = 0.
    # from adapting shift() - https://stackoverflow.com/questions/41399538/comparing-previous-row-values-in-pandas-dataframe

    df['fresh_entry'] = np.where(((df['positions']==-1) | (df['positions']==1)) & (df['positions'].shift()==0), True, False)

    # Calculate fresh exits
    df['fresh_exit'] = np.where(((df['positions'].shift()==-1) | (df['positions'].shift()==1)) & (df['positions']==0), True, False)

    # Calculate one way butterfly fees
    df['one_way_fee'] = np.where(((df['fresh_entry'] == True) | (df['fresh_exit']== True)), one_way_fee, 0)

    # Calculate net strategy returns, after fees
    df['net_strategy_returns'] = df["strategy_returns"] + df["one_way_fee"]

    df = df.dropna()

    # Calculate cumulative strategy returns
    df['cumret'] = ((df.net_strategy_returns+1)).cumprod()

    # print("{0}, {1}, {2} Butterfly Returns with Timeframe {3}, Lookback {4}, Std Deviation {5}".format(front_leg_instrument, middle_leg_instrument, back_leg_instrument, tf,lookback, std_dev ))

    #PLOTTING NOT COMPATIBLE WITH .PY SCRIPT IN TERMINAL
    # df.cumret.plot(figsize=(10,5))
    # plt.title("{0}, {1}, {2} Butterfly Returns with Timeframe {3}, Lookback {4}, Std Deviation {5}".format(front_leg_instrument, middle_leg_instrument, back_leg_instrument, tf,lookback, std_dev ))
    # plt.xlabel('Date')
    # plt.ylabel('Cumulative Returns')
    # plt.legend()
    # plt.show()

    length = len(df.index)
    last_index = len(df.index) - 1
    # print(last_index)
    ROI = df['cumret'].iloc[last_index]
    ROI = ROI.item()

    # print(f'ROI is {df["cumret"].iloc[last_index]} for {tf}')
    # print(f'length of {tf} dataframe is {len(df)}')

    #PLOTTING NOT COMPATIBLE WITH .PY SCRIPT IN TERMINAL
    # df[['spread','moving_average',"lower_band",'upper_band']].plot(figsize=(10,4))
    # plt.grid(True)
    # plt.title("{0}, {1}, {2} Butterfly Bollinger Bands with Timeframe {3}".format(front_leg_instrument, middle_leg_instrument, back_leg_instrument, tf))
    # plt.axis('tight')
    # plt.ylabel('Price')

    output_string = {
        "front_leg": front_leg_instrument,
        "middle_leg": middle_leg_instrument,
        "back_leg" : back_leg_instrument,
        "tf": tf,
        "ROI": ROI,
        "length": length,
        "lookback": lookback,
        "std_dev": std_dev,
        "front_vector": front_vector,
        "middle_vector": middle_vector,
        "back_vector": back_vector,
        "backtest_created_timestamp": now_readable,
        "backtest_alt_timestamp": now,
    }
    # print(df.head())

    return [output_string, df]
    # return df_json

df = triple_df_builder(df_object)
# get vectors
vectors = johansen_test(df)
# print('printing all vectors...')
# print(vectors)
front_vector = vectors[0]
middle_vector = vectors[1]
back_vector = vectors[2]

single_result = butterfly_sl_backtester(df, front_vector, middle_vector, back_vector, 20, 2, True, 1, -0.0015)

# get length of dataframe

select_df = single_result[1]
exportable_df = select_df["cumret"]
df_length = len(exportable_df.index)
# print(df_length)

#to reduce size, only select every nth row to achieve 100 rows only
# based on https://stackoverflow.com/questions/25055712/pandas-every-nth-row
nth_factor = int(df_length/100)-1
# print(nth_factor)

df = exportable_df.iloc[::nth_factor]

json_df = df.to_json(f'./data/{now_readable}_{tf}_cumret.json', orient = 'columns')
# df.to_csv(f'./data/{now}_{tf}_running.csv')

# trying ot read exported json file for sys.stdout.write sending
# https://www.geeksforgeeks.org/read-json-file-using-python/
# f = open(f"./data/{now}_{tf}_running.json")
# print(json_df)
# json_df = json.load(f)

# get output string
# print(single_result[0])
sys.stdout.write(json.dumps(single_result[0]))

# dataframe to json
# https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.to_json.html
# json_df = single_result[1].to_json(orient = "columns")
# print(str(json_df))

# sys.stdout.write(json.dumps(json_df))


text:
import matplotlib.image as mpimg
import matplotlib.pyplot as plt
import numpy as np
import cv2
import glob
import time
import pickle
from sklearn.svm import LinearSVC
from sklearn.preprocessing import StandardScaler
from utilities import extract_features, get_hog_features, find_cars, find_cars_2, add_heat, apply_threshold, draw_labeled_bboxes, draw_res
import matplotlib.image as mpimg
import matplotlib.pyplot as plt
import numpy as np
import cv2
from skimage import img_as_float
from scipy.ndimage.measurements import label
from moviepy.editor import VideoFileClip

class LabelsQueue():
    def __init__(self):
        # Number labels to store
        self.queue_len = 10
        self.queue = []

    # Put new frame
    def put_labels(self, labels):
        if (len(self.queue) > self.queue_len):
            tmp = self.queue.pop(0)
        self.queue.append(labels)

    # Get last N frames hot boxes
    def get_labels(self):
        b = []
        for label in self.queue:
            b.extend(label)
        return b

print('Loading Classifier parameters...')
# load a pe-trained svc model from a serialized (pickle) file
dist_pickle = pickle.load( open("svc_pickle.p", "rb" ) )

# get attributes of our svc object
svc = dist_pickle["svc"]
X_scaler = dist_pickle["scaler"]
orient = dist_pickle["orient"]
pix_per_cell = dist_pickle["pix_per_cell"]
cell_per_block = dist_pickle["cell_per_block"]
spatial_size = dist_pickle["spatial_size"]
hist_bins = dist_pickle["hist_bins"]

last_hot_labels = LabelsQueue()
img = mpimg.imread('test_images/test5.jpg')

def cars_detection(img):
    out_img, win_pos = find_cars_2(img, svc, X_scaler, orient,'YCrCb', pix_per_cell,
                                 cell_per_block, spatial_size, hist_bins)
    # Read in image similar to one shown above

    plt.imshow(out_img)
    plt.savefig('output_car_found_test5.jpg')

    heat = np.zeros_like(img[:, :, 0]).astype(np.float)

    # AVG boxes
    last_hot_labels.put_labels(win_pos)
    win_pos = last_hot_labels.get_labels()

    # Add heat to each box in box list
    heat = add_heat(heat, win_pos)

    plt.imshow(heat, cmap='hot')
    plt.savefig('heatmap_before_test5.jpg')
    # Apply threshold to help remove false positives
    heat = apply_threshold(heat, 2)
    plt.imshow(heat, cmap='hot')
    plt.savefig('heatmap_after_test5.jpg')

    # Visualize the heatmap when displaying
    heatmap = np.clip(heat, 0, 255)
    # Find final boxes from heatmap using label function
    labels = label(heatmap)
    #draw_img = draw_labeled_bboxes(np.copy(img), labels)

    # plot_row2(img, out_img, 'Source', 'All Detections')
    # plot_row2(draw_img, heatmap, 'Car Positions', 'Heat Map')
    return labels, heatmap

def image_pipeline(image):
    # Vehicle detection pipeline
    labels, heatmap = cars_detection(image)

    # Create an output image to draw on and  visualize the result
    out_img = image #np.dstack((image, image, image))*255
    # Draw output
    processed_img = draw_res(image, labels)
    return processed_img

# out_img, xpto = find_cars_2(img, svc, X_scaler, orient, pix_per_cell, cell_per_block, spatial_size,
#                    hist_bins)
# out_img, xpto = find_cars(img,ystart,ystop,scale, svc, X_scaler, orient, pix_per_cell, cell_per_block, spatial_size,
#                     hist_bins)

# last_hot_labels = LabelsQueue()
# output_video = 'video_result.mp4'
# clip1 = VideoFileClip("project_video.mp4")#.subclip(0,3)
# video_clip = clip1.fl_image(image_pipeline)
# video_clip.write_videofile(output_video, audio=False)

out_img = image_pipeline(img)

plt.imshow(out_img)

#print(out_img)
plt.savefig('car_found_test5_processed.jpg')

ext: py

reason: import numpy,from scipy


text:
import os
import sys
from functools import partial
import os
import pickle
import sys
import torch
from copy import deepcopy
import numpy as np
import matplotlib.pyplot as plt
import csv
from collections import defaultdict
import math
import glob
import re

from safe_control_gym.utils.configuration import ConfigFactory
from safe_control_gym.utils.plotting import plot_from_logs, window_func, load_from_log_file, COLORS, LINE_STYLES
from safe_control_gym.utils.registration import make
from safe_control_gym.utils.utils import mkdirs, set_dir_from_config, set_device_from_config, set_seed_from_config, save_video

##############################################################################################
#############################   Testing Utils   ##############################################
##############################################################################################

def make_fixed_init_env(config, env_seed=None, n_episodes=10, init_state=None):
    """"""
    assert init_state is not None, "Must provide an init_state dict for the env."
    config = deepcopy(config_origin)
    # modify env config with fixed init
    config.task_config.randomized_init = False
    config.task_config.init_state_randomization_info = None
    # use deep_set since init_state is None, cannot assign
    deep_set(config, "task_config.init_state", {})
    for state_k, state_v in init_state.items():
        deep_set(config, "task_config.init_state.{}".format(state_k), state_v)
    # make env
    env_func = partial(make, config.task, seed=env_seed, output_dir=config.output_dir, **config.task_config)
    env = env_func()
    env = RecordEpisodeStatistics(env, n_episodes)
    env.add_tracker("constraint_violation", 0, mode="queue")
    env.add_tracker("constraint_values", 0, mode="queue")
    env.add_tracker("mse", 0, mode="queue")
    return env

def run_with_env(config, ctrl, env, render=False, n_episodes=10, verbose=False):
    """Test policy with a given environment and save test trajectories."""
    ep_obs_list, ep_act_list = [], []
    ep_returns, ep_lengths = [], []
    frames = []

    obs, info = env.reset()
    ep_obs, ep_act = [deepcopy(obs[:env.state_dim])], []

    obs = ctrl.obs_normalizer(obs)
    c = info.get("constraint_values", 0)

    while len(ep_returns) < n_episodes:
        with torch.no_grad():
            obs = torch.FloatTensor(obs).to(ctrl.device)
            c = torch.FloatTensor(c).to(ctrl.device)

            if config.algo == "sac" or config.algo == "ddpg":
                action = ctrl.agent.ac.act(obs, deterministic=True)
            elif config.algo == "ppo":
                action = ctrl.agent.ac.act(obs)
            elif config.algo == "safe_explorer_ppo":
                action = ctrl.agent.ac.act(obs, c=c)
            else:
                raise NotImplementedError("ac.act not implemented.")

        obs, reward, done, info = env.step(action)
        ep_obs.append(deepcopy(obs[:env.state_dim]))
        ep_act.append(deepcopy(action))

        if render:
            env.render()
            frames.append(env.render("rgb_array"))
        if verbose:
            print("obs {} | act {}".format(obs, action))

        if done:
            assert "episode" in info
            ep_returns.append(info["episode"]["r"])
            ep_lengths.append(info["episode"]["l"])
            obs, _ = env.reset()

            ep_obs_list.append(np.stack(ep_obs))
            ep_act_list.append(np.stack(ep_act))
            ep_obs, ep_act = [deepcopy(obs[:env.state_dim])], []

        obs = ctrl.obs_normalizer(obs)
        c = info.get("constraint_values", 0)

    # collect evaluation results
    ep_lengths = np.asarray(ep_lengths)
    ep_returns = np.asarray(ep_returns)
    eval_results = {
        "ep_returns": ep_returns,
        "ep_lengths": ep_lengths,
        "ep_obs_list": ep_obs_list,
        "ep_act_list": ep_act_list,
    }
    if len(frames) > 0:
        eval_results["frames"] = frames
    # Other episodic stats from evaluation env.
    if len(env.queued_stats) > 0:
        queued_stats = {k: np.asarray(v) for k, v in env.queued_stats.items()}
        eval_results.update(queued_stats)
    return eval_results

def save_traj_to_csv(config, env, results, save_dir):
    """ """
    ep_obs_list = results["ep_obs_list"]
    ep_act_list = results["ep_act_list"]
    header = ["t", "x", "x_dot", "z", "z_dot", "theta", "theta_dot", "action1", "action2"]
    dim = 9
    state_dim = 6
    action_dim = 2
    n_episodes = len(ep_obs_list)
    os.makedirs(save_dir, exist_ok=True)
    for i in range(n_episodes):
        # (T+1, state_dim)
        states = ep_obs_list[i]
        # (T, action_dim)
        actions = ep_act_list[i]
        Tp1 = len(states)
        x = np.arange(Tp1)
        data = np.zeros((Tp1, dim))
        data[:, 0] = x
        data[:, 1:1+state_dim] = states
        thrust = np.clip(actions, env.action_space.low, env.action_space.high)
        if config.task_config.normalized_rl_action_space:
            thrust = (1 + env.norm_act_scale * thrust) * env.hover_thrust
        data[:-1, -action_dim:] = thrust
        rows = data.tolist()
        csv_path = os.path.join(save_dir, "traj_{}.csv".format(i))
        with open(csv_path, "w") as f:
            writer = csv.writer(f)
            writer.writerow(header)
            writer.writerows(rows)
    print("Test trajectories (time) saved to csv files.")

def plot_time_traj(env, results, output_path="temp.png", max_episode_plot=1):
    """Plot each dimension of state and action in time."""
    ep_obs_list = results["ep_obs_list"]
    ep_act_list = results["ep_act_list"]
    n_episodes = len(ep_obs_list)
    # process trajectories
    # (T, state_dim)
    traj_goal = env.X_GOAL
    if env.QUAD_TYPE == 1:
        state_names = ["z", "z_dot"]
        state_dim = 2
        act_dim = 1
    elif env.QUAD_TYPE == 2:
        state_names = ["x", "x_dot", "z", "z_dot", "theta", "theta_dot"]
        state_dim = 6
        act_dim = 2
    ncols = state_dim + 1
    n_rows = min(n_episodes, max_episode_plot)
    # fig = plt.figure()
    h_size = ncols * 4
    v_size = n_rows * 3
    fig = plt.figure(figsize=(h_size, v_size))
    axes = fig.subplots(nrows=n_rows, ncols=ncols)

    for i in range(n_rows):
        if n_rows == 1:
            axe = axes
        else:
            axe = axes[i]
        # (T+1, state_dim), (T, act_dim)
        traj, acts = ep_obs_list[i], ep_act_list[i]
        x = list(range(len(traj)))

        for j in range(state_dim):
            traj_goal_dim = traj_goal[:, j]
            traj_dim = traj[:, j]
            axe[j].plot(x[1:], traj_goal_dim, color="green")
            axe[j].plot(x, traj_dim, color="black")
            if i == 0:
                axe[j].set_title(state_names[j])

        for k in range(act_dim):
            action_dim = acts[:, k]
            axe[-1].plot(x[:-1], action_dim)
            if i == 0:
                axe[-1].set_title("action")
    # save fig
    plt.tight_layout()
    plt.savefig(output_path)
    plt.show()
    print("Test trajectories (time) plotted.")

def plot_phase_traj(env, results, output_path="temp.png", max_episode_plot=1):
    """Plot the trajectory in phase (state/velocity) spaces."""
    ep_obs_list = results["ep_obs_list"]
    ep_act_list = results["ep_act_list"]
    n_episodes = len(ep_obs_list)
    marker_size = 100
    # process trajectories
    # (T, state_dim)
    traj_goal = env.X_GOAL
    if env.QUAD_TYPE == 1:
        state_names = ["z", "z_dot"]
        state_dim = 2
        ncols = 2
    elif env.QUAD_TYPE == 2:
        state_names = ["x", "x_dot", "z", "z_dot", "theta", "theta_dot"]
        state_dim = 6
        ncols = 3
    n_rows = min(n_episodes, max_episode_plot)
    # fig = plt.figure()
    h_size = ncols * 4
    v_size = n_rows * 3
    fig = plt.figure(figsize=(h_size, v_size))
    axes = fig.subplots(nrows=n_rows, ncols=ncols)

    for i in range(n_rows):
        if n_rows == 1:
            axe = axes
        else:
            axe = axes[i]
        # (T+1, state_dim), (T, act_dim)
        traj = ep_obs_list[i]

        # (T+1 or T, 2)
        if env.QUAD_TYPE == 1:
            z_goal, z_dot_doal = traj_goal[:, 0], traj_goal[:, 1]
            z, z_dot = traj[:, 0], traj[:, 1]
            x = np.zeros_like(z)

            axe[0].plot(x[1:], z_goal, color="green")
            # axe[0].scatter(x[-1], z_goal[-1], color="green", s=marker_size)
            axe[0].plot(x, z, color="black")
            axe[0].scatter(x[-1], z[-1], color="red", s=marker_size//2)
            axe[0].set_xlabel("x")
            axe[0].set_ylabel("z")
            axe[0].set_title("z traj")

            axe[1].plot(x[1:], z_dot_doal, color="green")
            # axe[1].scatter(x[-1], z_dot_doal[-1], color="green", s=marker_size)
            axe[1].plot(x, z_dot, color="black")
            axe[1].scatter(x[-1], z_dot[-1], color="red", s=marker_size//2)
            axe[1].set_xlabel("x_dot")
            axe[1].set_ylabel("z_dot")
            axe[1].set_title("z_dot traj")

        elif env.QUAD_TYPE == 2:
            x_goal, x_dot_goal, z_goal, z_dot_doal = traj_goal[:, 0], traj_goal[:, 1], traj_goal[:, 2], traj_goal[:, 3]
            x, x_dot, z, z_dot = traj[:, 0], traj[:, 1], traj[:, 2], traj[:, 3]
            theta, theta_dot = traj[:, 4], traj[:, 5]

            axe[0].plot(x_goal, z_goal, color="green")
            # axe[0].scatter(x_goal[-1], z_goal[-1], color="green", s=marker_size)
            axe[0].plot(x, z, color="black")
            axe[0].scatter(x[-1], z[-1], color="red", s=marker_size//2)
            axe[0].set_xlabel("x")
            axe[0].set_ylabel("z")
            axe[0].set_title("x-z traj")

            axe[1].plot(x_dot_goal, z_dot_doal, color="green")
            # axe[1].scatter(x_dot_goal[-1], z_dot_doal[-1], color="green", s=marker_size)
            axe[1].plot(x_dot, z_dot, color="black")
            axe[1].scatter(x_dot[-1], z_dot[-1], color="red", s=marker_size//2)
            axe[1].set_xlabel("x_dot")
            axe[1].set_ylabel("z_dot")
            axe[1].set_title("x_dot-z_dot traj")

            axe[2].plot(theta, theta_dot, color="black")
            axe[2].scatter(theta[-1], theta_dot[-1], color="red", s=marker_size//2)
            axe[2].set_xlabel("theta")
            axe[2].set_ylabel("theta_dot")
            axe[2].set_title("theta-theta_dot traj")
    # save fig
    plt.tight_layout()
    plt.savefig(output_path)
    plt.show()
    print("Test trajectories (phase) plotted.")

##############################################################################################
#############################   Plotting Utils   #############################################
##############################################################################################

def reward_to_cost(y):
    """Converts RL reward to control cost (for plotting), 250 is total reward/return upper bound."""
    return 250 - y

def load_stats(legend_dir_specs,
               scalar_names=[],
               window=None,
               x_num_max=None,
               x_rescale_factor=None):
    """Gets all processed statistics for multiple scalars."""
    scalar_stats = {}
    for scalar_name in scalar_names:
        # Get all stats.
        stats = defaultdict(list)
        for l, dirs in legend_dir_specs.items():
            for d in dirs:
                # Pick from either log source (tensorboard or log text files).
                path = os.path.join(d, "logs", scalar_name + ".log")
                _, x, _, y = load_from_log_file(path)
                x, y = np.asarray(x), np.asarray(y)
                # quick hack to convert reward to cost (consistent with other control methods)
                if "return" in scalar_name.lower() or "reward" in scalar_name.lower():
                    y = reward_to_cost(y)
                # rescale the x-axis, e.g. convert to training step to time
                if x_rescale_factor is not None:
                    x *= x_rescale_factor
                # Smoothing.
                if window:
                    x, y = window_func(x, y, window, np.mean)
                stats[l].append([x, y])
        # Post-processing.
        x_max = float("inf")
        for _, runs in stats.items():
            for x, y in runs:
                # Align length of x data (get min across all runs & all algos).
                x_max = min(x_max, len(x))
        if x_num_max:
            x_max = min(x_max, x_num_max)
        processed_stats = {}
        for name, runs in stats.items():
            # Use same x for all runs to an algo.
            x = np.array([x[:x_max] for x, _ in runs])[0]
            # Different y for different runs.
            y = np.stack([y[:x_max] for _, y in runs])
            # Record stats.
            processed_stats[name] = {
                "x": x,
                "mean": np.mean(y, axis=0),
                "std": np.std(y, axis=0),
                "median": np.median(y, axis=0),
                "quantiles": np.quantile(y, [0.25, 0.75], axis=0),
                "min": np.min(y, axis=0),
                "max": np.max(y, axis=0),
            }
        # Copy over stats.
        scalar_stats[scalar_name] = deepcopy(processed_stats)
    return scalar_stats

def plot_from_exps(legend_dir_specs,
                    out_path="temp.jpg",
                    scalar_names=[],
                    title="Traing Curves",
                    sub_titles=["Loss Curve"],
                    xlabel="Epochs",
                    ylabels=["Loss"],
                    window=None,
                    x_num_max=None,
                    x_rescale_factor=None,
                    num_std=1,
                    use_median_quantile=False,
                    cols_per_row=3):
    """Plots 1 statistic figure at a time."""
    # Get all stats.
    scalar_stats = load_stats(legend_dir_specs,
                              scalar_names=scalar_names,
                              window=window,
                              x_num_max=x_num_max,
                              x_rescale_factor=x_rescale_factor)
    # Make plots.
    num_plots = len(scalar_stats)
    num_rows = math.ceil(num_plots / cols_per_row)
    num_cols = min(num_plots, cols_per_row)
    fig = plt.figure(figsize=(12, 6))
    axes = fig.subplots(nrows=num_rows, ncols=num_cols)
    for idx, scalar_name in enumerate(scalar_names):
        row_idx = idx // num_rows
        col_idx = idx % num_cols
        if num_rows > 1:
            ax = axes[row_idx, col_idx]
        elif num_cols > 1:
            ax = axes[col_idx]
        else:
            ax = axes
        processed_stats = scalar_stats[scalar_name]
        for i, name in enumerate(sorted(processed_stats.keys())):
            color_i = (i + 0) % len(COLORS)
            color = COLORS[color_i]
            line_i = (i + 0) % len(LINE_STYLES)
            linestyle = LINE_STYLES[line_i][-1]
            x = processed_stats[name]["x"]
            if use_median_quantile:
                y_median = processed_stats[name]["median"]
                y_quantiles = processed_stats[name]["quantiles"]
                y_quant_1st = y_quantiles[0]
                y_quant_3rd = y_quantiles[1]
                ax.plot(x, y_median, label=name, color=color, alpha=0.7, linestyle=linestyle)
                ax.fill_between(x, y_quant_3rd, y_quant_1st, alpha=0.1, color=color)
            else:
                y_mean = processed_stats[name]["mean"]
                y_std = processed_stats[name]["std"]
                ax.plot(x, y_mean, label=name, color=color, alpha=0.7, linestyle=linestyle)
                ax.fill_between(x, y_mean + num_std * y_std, y_mean - num_std * y_std, alpha=0.1, color=color)
        ax.set_xlabel(xlabel)
        ax.set_ylabel(ylabels[idx])
        # ax.set_ylim((-10, None))
    # Postprocess plot.
    fig.suptitle(title)
    fig.subplots_adjust(bottom=0.15)
    lines = []
    labels = []
    for ax in fig.axes:
        axLine, axLabel = ax.get_legend_handles_labels()
        lines.extend(axLine)
        labels.extend(axLabel)
        break
    fig.legend(lines,
               labels,
               loc='lower center',
               fancybox=True,
               shadow=True,
               borderaxespad=0.1,
               ncol=7)
    plt.savefig(out_path)
    plt.show()
    return scalar_stats

def save_stats_to_csv(scalar_stats,
                      algo_name_map=None,
                      scalar_name_map=None,
                      csv_path="stats.csv"):
    """Saves the queried experiment statistics to a csv file."""
    curves = ["mean", "std", "min", "max", "median", "top_quartile", "bottom_quartile"]
    header = []
    stat_rows = []
    x_already = False
    for scalar_name in scalar_stats:
        stats = scalar_stats[scalar_name]
        true_scalar_name = scalar_name_map[scalar_name] if scalar_name_map else scalar_name
        # Collect stats.
        for algo_name in sorted(stats.keys()):
            true_algo_name = algo_name_map[algo_name] if algo_name_map else algo_name
            stat = stats[algo_name]
            # X.
            if not x_already:
                header.append("x-Step")
                stat_rows.append(stat["x"])
                x_already = True
            # Y.
            header.extend(["y-{}-{}-{}".format(true_scalar_name, true_algo_name, c) for c in curves])
            stat_rows.append(stat["mean"])
            stat_rows.append(stat["std"])
            stat_rows.append(stat["min"])
            stat_rows.append(stat["max"])
            stat_rows.append(stat["median"])
            stat_rows.append(stat["quantiles"][1])
            stat_rows.append(stat["quantiles"][0])
    # Make rows.
    stat_mtx = np.array(stat_rows).transpose()
    rows = stat_mtx.tolist()
    # Write to csv.
    os.makedirs(os.path.dirname(csv_path), exist_ok=True)
    with open(csv_path, "w") as f:
        writer = csv.writer(f)
        writer.writerow(header)
        writer.writerows(rows)

max_stars_repo_name: catgloss/safe-control-gym

max_issues_repo_name: catgloss/safe-control-gym

max_forks_repo_name: catgloss/safe-control-gym

reason: import numpy


text:
# Copyright 2021 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#      http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""Helper functions for building pseudospectral methods."""

from typing import Callable, Tuple

import jax.numpy as jnp
from jax_cfd.base import grids
from jax_cfd.spectral import types as spectral_types

def truncated_rfft(u: spectral_types.Array) -> spectral_types.Array:
  """Applies the 2/3 rule by truncating higher Fourier modes.

  Args:
    u: the real-space representation of the input signal

  Returns:
    Downsampled version of `u` in rfft-space.
  """
  uhat = jnp.fft.rfft(u)
  k, = uhat.shape
  final_size = int(2 / 3 * k) + 1
  return 2 / 3 * uhat[:final_size]

def padded_irfft(uhat: spectral_types.Array) -> spectral_types.Array:
  """Applies the 3/2 rule by padding with zeros.

  Args:
    uhat: the rfft representation of a signal

  Returns:
    An upsampled signal in real space which 3/2 times larger than the input
    signal `uhat`.
  """
  n, = uhat.shape
  final_shape = int(3 / 2 * n)
  smoothed = jnp.pad(uhat, (0, final_shape - n))
  assert smoothed.shape == (final_shape,), "incorrect padded shape"
  return 1.5 * jnp.fft.irfft(smoothed)

def circular_filter_2d(grid: grids.Grid) -> spectral_types.Array:
  """Circular filter which roughly matches the 2/3 rule but is smoother.

  Follows the technique described in Equation 1 of [1]. We use a different value
  for alpha as used by pyqg [2].

  Args:
    grid: the grid to filter over

  Returns:
    Filter mask

  Reference:
    [1] Arbic, Brian K., and Glenn R. Flierl. "Coherent vortices and kinetic
    energy ribbons in asymptotic, quasi two-dimensional f-plane turbulence."
    Physics of Fluids 15, no. 8 (2003): 2177-2189.
    https://doi.org/10.1063/1.1582183

    [2] Ryan Abernathey, rochanotes, Malte Jansen, Francis J. Poulin, Navid C.
    Constantinou, Dhruv Balwada, Anirban Sinha, Mike Bueti, James Penn,
    Christopher L. Pitt Wolfe, & Bia Villas Boas. (2019). pyqg/pyqg: v0.3.0
    (v0.3.0). Zenodo. https://doi.org/10.5281/zenodo.3551326.
    See:
    https://github.com/pyqg/pyqg/blob/02e8e713660d6b2043410f2fef6a186a7cb225a6/pyqg/model.py#L136
  """
  kx, ky = grid.rfft_mesh()
  max_k = ky[-1, -1]

  circle = jnp.sqrt(kx**2 + ky**2)
  cphi = 0.65 * max_k
  filterfac = 23.6
  filter_ = jnp.exp(-filterfac * (circle - cphi)**4.)
  filter_ = jnp.where(circle <= cphi, jnp.ones_like(filter_), filter_)
  return filter_

def vorticity_to_velocity(
    grid: grids.Grid
) -> Callable[[spectral_types.Array], Tuple[spectral_types.Array,
                                            spectral_types.Array]]:
  """Constructs a function for converting vorticity to velocity, both in Fourier domain.

  Solves for the stream function and then uses the stream function to compute
  the velocity. This is the standard approach. A quick sketch can be found in
  [1].

  Args:
    grid: the grid underlying the vorticity field.

  Returns:
    A function that takes a vorticity (rfftn) and returns a velocity vector
    field.

  Reference:
    [1] Z. Yin, H.J.H. Clercx, D.C. Montgomery, An easily implemented task-based
    parallel scheme for the Fourier pseudospectral solver applied to 2D
    Navier–Stokes turbulence, Computers & Fluids, Volume 33, Issue 4, 2004,
    Pages 509-520, ISSN 0045-7930,
    https://doi.org/10.1016/j.compfluid.2003.06.003.
  """
  kx, ky = grid.rfft_mesh()
  two_pi_i = 2 * jnp.pi * 1j
  laplace = two_pi_i ** 2 * (abs(kx)**2 + abs(ky)**2)
  laplace = laplace.at[0, 0].set(1)

  def ret(vorticity_hat):
    psi_hat = -1 / laplace * vorticity_hat
    vxhat = two_pi_i * ky * psi_hat
    vyhat = -two_pi_i * kx * psi_hat
    return vxhat, vyhat

  return ret

def filter_step(step_fn: spectral_types.StepFn, filter_: spectral_types.Array):
  """Returns a filtered version of the step_fn."""
  def new_step_fn(state):
    return filter_ * step_fn(state)
  return new_step_fn

max_stars_repo_name: ngam/jax-cfd

max_issues_repo_name: ngam/jax-cfd

max_forks_repo_name: ngam/jax-cfd

reason: import jax,from jax


text:
from core.calib.kruppa import KruppaSolver, KruppaSolverMC, KruppaSolverRANSAC
from core.calib.focal import FocalSolverStrum

from opt.dist import DistSolver
from core.match import Matcher
from utils import cv_wrap as W
from utils import vmath as M
import viz as V
import cv2
import numpy as np

import autograd.numpy as anp
from autograd import jacobian

def mcheck(x):
    if x is None:
        return False
    if np.any(np.isnan(x)):
        return False
    return True

class FundCov(object):
    def __init__(self):
        self.pt_a_ = None
        self.pt_b_ = None
        self.jac = jacobian(self.err_anp)

    def err_anp(self, F):
        return self.err(F, np=anp)

    def err(self, F, np=np):
        pt_a = self.pt_a_
        pt_b = self.pt_b_
        e = np.einsum('...a,ab,...b', M.to_h(pt_b), F.reshape(3,3), M.to_h(pt_a))
        return np.ravel(e)

    @staticmethod
    def jac_to_cov(J):
        """ from scipy/optimize/minpack.py#L739 """
        _, s, VT = np.linalg.svd(J, full_matrices=False)
        thresh = np.finfo(np.float32).eps * max(J.shape) * s[0]
        s = s[s > thresh]
        VT = VT[:s.size]
        cov = np.dot(VT.T / s**2, VT)
        return cov

    def __call__(self, F, pt_a, pt_b):
        self.pt_a_ = pt_a
        self.pt_b_ = pt_b
        J = self.jac(F.ravel())
        return self.jac_to_cov(J)

def jac_F(F, pt_a, pt_b, np=np):
    jacobian(err_F,self.err_anp)#, np=anp)
    return

def main():
    w, h = 640, 480
    foc_solver = FocalSolverStrum(w, h)
    cov_solver = FundCov()
    #solver = KruppaSolverMC()
    #solver = KruppaSolver()
    solver = KruppaSolverRANSAC(w,h,method='esv')
    K0 = np.float32([
        (w+h), 0.0, w/2.0,
        0.0, (w+h), h/2.0,
        0, 0, 1]).reshape(3,3)
    #K0 = np.float32([200,0,200,0,200,200,0,0,1]).reshape(3,3)
    #K0 = np.float32([
    #    1260,0,280,0,1260,230,0,0,1]).reshape(3,3)
    #feat = cv2.ORB_create(nfeatures=1024)
    feat = cv2.xfeatures2d.SURF_create()#nfeatures=1024)
    matcher = Matcher(des=feat)

    mcfg = Matcher.PRESET_HARD.copy()
    mcfg['maxd'] = 256.0

    imgs = np.load('/tmp/db_imgs.npy')
    kpts = np.load('/tmp/db_kpts.npy')
    dess = np.load('/tmp/db_dess.npy')
    Fs = []
    Ws = []

    for i in range(0, len(imgs)):
        print '{}/{}'.format(i, len(imgs))
        for j in range(i-8, i+8):
            if i==j: continue
            if j<0 or j>=len(imgs): continue

            img0, kpt0, des0 = [e[i] for e in [imgs,kpts,dess]]
            img1, kpt1, des1 = [e[j] for e in [imgs,kpts,dess]]

            i_m_h0, i_m_h1 = matcher(des0, des1,
                    **mcfg
                    )

            if len(i_m_h0) < 32:
                continue

            pt_a = kpt0[i_m_h0]
            pt_b = kpt1[i_m_h1]

            F, msk = W.F(pt_a, pt_b,
                    method=cv2.FM_RANSAC,
                    ransacReprojThreshold=2.0,
                    confidence=0.999
                    )
            #print 'rank', np.linalg.matrix_rank(F)
            n_in = np.count_nonzero(msk)
            r_in = n_in / float(msk.size)
            print 'n_in / r_in : {}, {}'.format(n_in, r_in)

            r_err = np.einsum('...a,ab,...b', M.to_h(pt_b), F, M.to_h(pt_a))
            r_err = np.sqrt(np.square(r_err).mean())
            if n_in > 64 and r_in > 0.8 and r_err < 1.0:
                mim = V.draw_matches(img0, img1, pt_a, pt_b)
                cv2.imshow('mim', mim)
                k = cv2.waitKey(1)
                if k == ord('q'):
                    break
                if k == 27:
                    return
                Fs.append( F )
                #print err_F(F, pt_a[msk], pt_b[msk])
                cov = cov_solver(F, pt_a[msk], pt_b[msk])
                Ws.append( 1.0 / np.diag(cov).sum() )

    print 'fs-len', len(Fs)

    Fs = np.asarray(Fs)
    Ws = np.asarray(Ws)
    Ws *= Ws.size / Ws.sum()

    # solve focal length first
    f = foc_solver(Fs)
    K0[0,0] = K0[1,1] = f
    print('updated Kmat through focal length initialization : {}'.format(K0))

    #Fs = Fs[np.random.choice(len(Fs), size=128)]

    #K1 = solver0(K0, Fs)
    ##if mcheck(K1):
    ##    K0 = K1
    #print 'K1', K1

    n_it, K1, inl = solver(K0, Fs, Ws)
    print 'n_it', n_it
    print 'K1', K1
    print 'inl', inl.sum() / float(inl.size)
    #if mcheck(K1):
    #    K0 = K1

if __name__ == "__main__":
    main()


text:
from Stats import ticks, tags, video_info
from plotly import graph_objs
from plotly.offline import plot
import matplotlib.pyplot as plt
plt.rcdefaults()
import numpy as np

def draw():

    #PIE CHART OF 30 MOST POPULAR TAGS
    pie_indeces = np.argpartition(ticks, -50)[-50:]
    pie_ticks = [int(ticks[x]) for x in pie_indeces]
    pie_tags = [str(tags[x]) for x in pie_indeces]

    frame = graph_objs.Pie(labels=pie_tags, values=pie_ticks)

    plot([frame], filename='html/YouTubeTags_PieChart.html')

    #GET THE 15 MOST POPULAR TAGS
    #GET THEIR INDECES TO ALSO GET CORRESPONDING TAG NAME
    indeces = np.argpartition(ticks, -20)[-20:]
    ticks_sub = [int(ticks[x]) for x in indeces]
    tags_sub = [str(tags[x]) for x in indeces]

    x = tags_sub
    y = ticks_sub

    data = [graph_objs.Bar(
        x=x, y=y,
        marker=dict(
            color='rgb(158,202,225)',
            line=dict(
                color='rgb(8,48,107)',
                width=1.5),
        ),
        opacity=0.6
    )]

    layout = graph_objs.Layout(
        annotations=[
            dict(x=xi, y=yi,
                 text=str(yi),
                 xanchor='center',
                 yanchor='bottom',
                 showarrow=False,
                 ) for xi, yi in zip(x, y)]
    )

    fig = graph_objs.Figure(data=data, layout=layout)
    plot(fig, filename='html/YouTubeTags_BarChart.html')

    #MAKE A VISUAL BASED ON MOST VIEWED VIDEOS
    #GRAB 15 MOST VIEWED VIDEOS
    #VIDEO_INFO: (TITLE, CREATOR, VIEWS, TAGS, DATE PUBLISHED)
    video_title, creator, views, tag, date = zip(*video_info)
    video_indeces = np.argpartition(views, -20)[-20:]
    total_views = list()
    total_tags = list()
    total_video_titles = list()

    for x in video_indeces:
        total_views.append(views[x])
        total_tags.append(tag[x])
        total_video_titles.append(video_title[x])

    c = 0
    for x in total_tags:
        total_tags[c] = ' '.join(x)
        c += 1
    trace0 = graph_objs.Bar(
        x=total_video_titles,
        y=total_views,
        text=total_tags,
        marker=dict(
            color='rgb(158,202,225)',
            line=dict(
                color='rgb(8,48,107)',
                width=1.5,
            )
        ),
        opacity=0.6
    )

    data = [trace0]
    layout = graph_objs.Layout(
        title='Most Popular Videos',
    )

    fig = graph_objs.Figure(data=data, layout=layout)
    plot(fig, filename='html/YouTubeTags_PopularVideos.html')

max_stars_repo_name: nashahzad/YoutubeTagsCrawler

max_issues_repo_name: nashahzad/YoutubeTagsCrawler

max_forks_repo_name: nashahzad/YoutubeTagsCrawler


text:
import numpy as np
import unittest
import ps6
import os
import cv2

INPUT_DIR = "input_images/input_test"

class PCA(unittest.TestCase):

    def test_mean_face(self):

        for i in range(1, 4):
            file_name = "x_data_mean_{}.npy".format(i)
            file_path = os.path.join(INPUT_DIR, file_name)
            x_data = np.load(file_path)

            file_name = "correct_mean_{}.npy".format(i)
            file_path = os.path.join(INPUT_DIR, file_name)
            x_mean = np.load(file_path)

            result = ps6.get_mean_face(x_data)

            correct = np.allclose(result, x_mean, atol=1)
            message = "Values do not match the reference solution. " \
                      "This function should only compute the mean of each " \
                      "column."
            self.assertTrue(correct, message)

    def test_pca(self):

        k_list = [3, 3, 4]

        for i in range(1, 4):

            k = k_list[i-1]

            file_name = "x_data_pca_k{}_{}.npy".format(k, i)
            file_path = os.path.join(INPUT_DIR, file_name)
            x_data = np.load(file_path)

            file_name = "ref_val_pca_k{}_{}.npy".format(k, i)
            file_path = os.path.join(INPUT_DIR, file_name)
            ref_val = np.load(file_path)

            eig_vecs, eig_vals = ps6.pca(x_data, k)

            correct = eig_vals.size == k
            message = "Wrong number of eigenvalues. K used: {}. " \
                      "Returned eigenvalues: {}".format(k, eig_vals)
            self.assertTrue(correct, message)

            for j in range(k):
                correct = np.allclose(ref_val[j], eig_vals[j])
                message = "Wrong eigenvalue at position {}".format(i)
                self.assertTrue(correct, message)

class Boosting(unittest.TestCase):

    def setUp(self):
        self.threshold_list = [[87, 62, 122, 59, 125],
                               [138, 69, 54, 117, 61],
                               [65, 108,  95,  82, 132]]

        self.y_pred_ref = ["ypred_boosting_1.npy", "ypred_boosting_2.npy",
                           "ypred_boosting_3.npy"]

        self.x_data = np.load(os.path.join(INPUT_DIR, "x_boosting_1.npy"))
        self.alphas = np.load(os.path.join(INPUT_DIR, "alpha_1.npy"))

    def test_predict(self):
        thresh_vals_list = self.threshold_list
        y_pred_ref_list = self.y_pred_ref
        x_data = self.x_data
        alphas = self.alphas

        for thresh_vals, y_ref_file in zip(thresh_vals_list, y_pred_ref_list):
            y_pred_ref = np.load(os.path.join(INPUT_DIR, y_ref_file))

            # Create a list of untrained classifiers to verify
            # predict operation.
            classifiers = [ps6.WeakClassifier([], [], [], thresh=t) for t in
                           thresh_vals]

            ytrain = np.zeros((x_data.shape[0],))  # Doesn't matter for this test
            numIte = 0  # Doesn't matter for this test

            test = ps6.Boosting(x_data, ytrain, numIte)

            # Modifying object variables for testing
            test.weakClassifiers = classifiers
            test.alphas = alphas

            y_pred = test.predict(x_data)

            correct = np.allclose(y_pred, y_pred_ref)
            message = "Predictions do not match reference predictions."
            self.assertTrue(correct, message)

    def test_evaluate(self):
        thresh_vals_list = self.threshold_list

        x_data = self.x_data
        alphas = self.alphas

        corr_inc_list = [(277, 283), (279, 281), (274, 286)]
        y_data = np.load(os.path.join(INPUT_DIR, "y_boosting_1.npy"))

        iter_lists = zip(thresh_vals_list, corr_inc_list)
        for thresh_vals, corr_inc in iter_lists:

            # Create a list of untrained classifiers to verify
            # predict operation.
            classifiers = [ps6.WeakClassifier([], [], [], thresh=t) for t in
                           thresh_vals]

            ytrain = np.zeros((x_data.shape[0],))  # Doesn't matter for this test
            numIte = 0  # Doesn't matter for this test

            test = ps6.Boosting(x_data, ytrain, numIte)

            # Modifying object variables for testing
            test.weakClassifiers = classifiers
            test.alphas = alphas
            test.Xtrain = x_data
            test.ytrain = y_data

            correct, incorrect = test.evaluate()

            correct_test = np.allclose(correct, corr_inc[0])
            message = "Correct values do not match the reference."
            self.assertTrue(correct_test, message)

            incorrect_test = np.allclose(incorrect, corr_inc[1])
            message = "Incorrect values do not match the reference."
            self.assertTrue(incorrect_test, message)

class HaarFeature(unittest.TestCase):

    def setUp(self):
        self.input_dir = os.path.join("input_images", "input_test")

    def test_preview(self):
        feat_types = 2 * [(1, 2)] + 2 * [(2, 1)] + 2 * [(1, 3)] \
                     + 2 * [(3, 1)] + 2 * [(2, 2)]
        positions = 5 * [(6, 9), (8, 5)]
        sizes = 5 * [(14, 18), (11, 15)]
        samples = 5 * [0, 1]

        for j in range(10):
            feat_type = feat_types[j]
            pos = positions[j]
            size = sizes[j]
            s = samples[j]
            file_name = "haar_preview_ft{}_pos{}_sz{}_{}.npy".format(feat_type, pos, size, s)

            ref_img = np.load(os.path.join(self.input_dir, file_name))

            # Uncomment if you want to see the reference image
            # cv2.imshow("ref_img", ref_img.astype("uint8"))
            # cv2.waitKey(0)

            hf = ps6.HaarFeature(feat_type, pos, size)
            hf_img = hf.preview((50, 50))

            correct = np.allclose(hf_img, ref_img)
            message = "Output image does not match the reference solution."
            self.assertTrue(correct, message)

    def test_integral_images(self):
        ti_path = os.path.join(self.input_dir, "test_image_ii.npy")
        ii_path = os.path.join(self.input_dir, "integral_image_ii.npy")

        test_image = np.load(ti_path)
        integral_image = np.load(ii_path)

        r = np.random.randint(1, integral_image.shape[0])
        c = np.random.randint(1, integral_image.shape[1])

        result = ps6.convert_images_to_integral_images([test_image])

        ref_val = integral_image[r, c]

        correct = np.allclose(ref_val, result[0][r, c])
        message = "Value at row:{}, col:{} is not correct. "
        self.assertTrue(correct, message)

    def test_HaarFeature_evaluate(self):
        ti_path = os.path.join(self.input_dir, "test_image_ii.npy")
        ii_path = os.path.join(self.input_dir, "integral_image_ii.npy")

        test_image = np.load(ti_path)
        integral_image = np.load(ii_path)

        feat_type = (2, 2)  # Change feature type
        pos = (5, 5)
        size = (30, 30)

        if feat_type == (2, 1):
            A = np.sum(test_image[pos[0]:pos[0] + size[0] / 2,
                                  pos[1]:pos[1] + size[1]])
            B = np.sum(test_image[pos[0] + size[0] / 2:pos[0] + size[0],
                                  pos[1]:pos[1] + size[1]])
            ref = A - B

        if feat_type == (1, 2):
            A = np.sum(test_image[pos[0]:pos[0] + size[0],
                                  pos[1]:pos[1] + size[1] / 2])
            B = np.sum(test_image[pos[0]:pos[0] + size[0],
                                  pos[1] + size[1] / 2:pos[1] + size[1]])
            ref = A - B

        if feat_type == (3, 1):
            A = np.sum(test_image[pos[0]:pos[0] + size[0] / 3,
                                  pos[1]:pos[1] + size[1]])
            B = np.sum(test_image[pos[0] + size[0] / 3:pos[0] + 2 * size[0] / 3,
                                  pos[1]:pos[1] + size[1]])
            C = np.sum(test_image[pos[0] + 2 * size[0] / 3:pos[0] + size[0],
                                  pos[1]:pos[1] + size[1]])
            ref = A - B + C

        if feat_type == (1, 3):
            A = np.sum(test_image[pos[0]:pos[0] + size[0],
                                  pos[1]:pos[1] + size[1] / 3])
            B = np.sum(test_image[pos[0]:pos[0] + size[0],
                                  pos[1] + size[1] / 3:pos[1] + 2 * size[1] / 3])
            C = np.sum(test_image[pos[0]:pos[0] + size[0],
                                  pos[1] + 2 * size[1] / 3:pos[1] + size[1]])
            ref = A - B + C

        if feat_type == (2, 2):
            A = np.sum(test_image[pos[0]:pos[0] + size[0] / 2,
                                  pos[1]:pos[1] + size[1] / 2])
            B = np.sum(test_image[pos[0]:pos[0] + size[0] / 2,
                                  pos[1] + size[1] / 2:pos[1] + size[1]])
            C = np.sum(test_image[pos[0] + size[0] / 2:pos[0] + size[0],
                                  pos[1]:pos[1] + size[1] / 2])
            D = np.sum(test_image[pos[0] + size[0] / 2:pos[0] + size[0],
                                  pos[1] + size[1] / 2:pos[1] + size[1]])
            ref = -A + B + C - D

        hf = ps6.HaarFeature(feat_type, pos, size)
        score = hf.evaluate(integral_image)

        correct = np.allclose(score, ref)

        message = "Wrong score returned."
        self.assertTrue(correct, message)

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

max_stars_repo_name: jperuggia/ComputerVision

max_issues_repo_name: jperuggia/ComputerVision

max_forks_repo_name: jperuggia/ComputerVision


text:
import numpy as np

class CrossEntropyLoss:
    def __init__(self, weight=None, size_average=True):
        """
        初始化参数，因为要实现 torch.nn.CrossEntropyLoss 的两个比较重要的参数

        :param weight: 给予每个类别不同的权重
        :param size_average: 是否要对 loss 求平均
        """

        self.weight = weight
        self.size_average = size_average

    def __call__(self, input, target):
        """
        计算损失
        这个方法让类的实例表现的像函数一样，像函数一样可以调用

        :param input: (batch_size, C)，C是类别的总数
        :param target: (batch_size, 1)
        :return: 损失
        """

        batch_loss = 0.0
        input = input.detach().numpy()
        target = target.detach().numpy()
        print(input.shape[0], "input.shape[0]")
        for i in range(input.shape[0]):

            numerator = np.exp(input[i, target[i]])  # 分子
            denominator = np.sum(np.exp(input[i, :]))  # 分母

            # 计算单个损失
            loss = -np.log(numerator / denominator)
            if self.weight:
                loss = self.weight[target[i]] * loss
            print("单个损失： ", loss)

            # 损失累加
            batch_loss += loss

        # 整个 batch 的总损失是否要求平均
        if self.size_average == True:
            batch_loss /= input.shape[0]

        return batch_loss

    def forward(self, input, target):
        return CrossEntropyLoss(input, target)


text:
import numpy as np

import matplotlib.pyplot as plt
from matplotlib import gridspec
import sys
import os
path = os.getcwd()
parent = os.path.join(path, os.pardir)
sys.path.append(os.path.abspath(parent))
from bayes_opt1 import BayesianOptimization
from bayes_opt1 import UtilityFunction
sys.path.append(os.path.abspath(os.path.join(parent, 'bayes_opt1')))
from util import _kg2

def normalization(utility):
    return (utility- min(utility)) / (max(utility) - min(utility))

def plot_gp(optimizers, x, target, params, n_init=2):

    x = x.reshape(-1,1)
    y = target(x)

    def posterior(optimizer, x_obs, y_obs, grid):
        optimizer._gp.fit(x_obs, y_obs)
        mu, sigma = optimizer._gp.predict(grid, return_std=True)
        return mu, sigma

    fig = plt.figure(figsize=(16, 10))
    steps = len(list(optimizers.values())[0].space)

    gs = gridspec.GridSpec(2, 1, height_ratios=[2, 2])
    axis = plt.subplot(gs[0])
    af = plt.subplot(gs[1])
    n = len(optimizers)
    colmap = plt.cm.get_cmap(name='rainbow', lut=n)

    x_obs = {}
    y_obs={}
    mu={}
    sigma={}
    utility={}
    axis.plot(x, y, linewidth=3, label='Target')

    i=0
    for acq, optimizer in optimizers.items():
        x_obs_i = np.array([[res["params"]["x"]] for res in optimizer.res])
        x_obs[acq] = x_obs_i
        y_obs_i = np.array([res["target"] for res in optimizer.res])
        y_obs[acq] =  y_obs_i
        mu[acq], sigma[acq] = posterior(optimizer, x_obs_i, y_obs_i, x)
        utility_function = UtilityFunction(kind=acq, kappa = params[acq]['kappa'], xi = params[acq]['xi'])
        if acq != 'kg':
            utility[acq] = utility_function.utility(x, optimizer, optimizer._gp, optimizer._space.target.max())
        else:
            utility[acq] = _kg2(x, optimizer, optimizer._gp, n_grid = 100, J = 30)
        axis.plot(x, mu[acq], '--', color=colmap(i), label='Prediction')
        i=i+1

    fig.suptitle(
        'Utility Functions After {} Iterations And {} Initial Points'.format(steps-n_init, n_init),
        fontdict={'size':50}
    )

    axis.set_xlim((min(x), max(x)))
    axis.set_ylim((None, None))
    axis.set_ylabel('f(x)', fontdict={'size':20})
    axis.set_xlabel('x', fontdict={'size':20})

    for i, acq in enumerate(utility):
        af.plot(x, normalization(utility[acq]), label=acq, color= colmap(i))
        af.plot(x[np.argmax(utility[acq])], np.max(normalization(utility[acq])), '*', markersize=15, markerfacecolor=colmap(i), markeredgecolor='k', markeredgewidth=1)

    af.set_xlim((min(x), max(x)))
    af.set_ylim((-0.2,1.2))
    af.set_ylabel('Utility', fontdict={'size':20})
    af.set_xlabel('x', fontdict={'size':20})

    axis.legend(loc=2, bbox_to_anchor=(1.01, 1), borderaxespad=0.)
    af.legend(loc=2, bbox_to_anchor=(1.01, 1), borderaxespad=0.)

    plt.show()

def plot_convergence(optimizers, x, target, params, it=20, init_point=2):
    x = x.reshape(-1,1)
    y = target(x)

    tar={}
    points={}

    utility_function={}
    for acq, optimizer in optimizers.items():
        utility_function[acq] = UtilityFunction(kind=acq, kappa = params[acq]['kappa'], xi = params[acq]['xi'])
        points[acq]=np.zeros(it)
        tar[acq]= np.zeros(it)
        for i in range(it):
            points[acq][i]=optimizer.suggest(utility_function[acq])['x']
            tar[acq][i]= target(points[acq][i])
            optimizer.register(params=points[acq][i], target=tar[acq][i])

    fig = plt.figure(figsize=(13, 6))

    af = plt.subplot()

    fig.suptitle(
        'Convergences To The Optimum After {} Iterations And {} Initial Points'.format(it, init_point),
        fontdict={'size':50}
    )

    steps = len(list(optimizers.values())[0].space)
    n = len(optimizers)
    colmap = plt.cm.get_cmap(name='rainbow', lut=n)

    num_iter=np.arange(1,it+1)

    for i, acq in enumerate(optimizers.keys()):
        af.plot(num_iter, points[acq], '*',markersize=15,markerfacecolor=colmap(i), markeredgecolor='k', markeredgewidth=1,label=acq,linestyle='solid',color=colmap(i))

    af.axhline(y=x[np.argmax(target(x))], linestyle=':', label='Optimum to be achieved')
    af.legend(loc=2, bbox_to_anchor=(1.01, 1), borderaxespad=0.)

    a = range(0,21)
    af.set_xticks(a)
    af.set_yticks(range(-2,11))
    plt.xlabel('Number of iterations')
    plt.ylabel('Suggested x')

    plt.show()

# Plot regret

def simple_regret(f_max, target_space):
    """
    The simple regret rT = max{x∈X} f(x) − max{t∈[1,T]} f(x_t) measures the value of
the best queried point so far.
    """
    return f_max-target_space.max()

def plot_simple_regret(optimizers, x, target, params, dim, it=20, init_point=2):

    tar={}
    points={}
    utility_function={}
    regrets={}

    if dim > 1:
        inputs = []
        for col in range(dim):
            inputs.append(x[:,col])
        y = target(*inputs)

        f_max=max(y)

        for acq, optimizer in optimizers.items():
            regrets[acq]=np.zeros(it)
            utility_function[acq] = UtilityFunction(kind=acq, kappa = params[acq]['kappa'], xi = params[acq]['xi'])

        for i in range(it):
            for acq, optimizer in optimizers.items():
                points[acq]=optimizer.suggest(utility_function[acq])
                tar[acq]= target(**points[acq])
                optimizer.register(params=points[acq], target=tar[acq])
                regrets[acq][i]= simple_regret(f_max, optimizer._space.target)
    else:
        x = x.reshape(-1,1)
        y = target(x)

        f_max=max(y)

        for acq, optimizer in optimizers.items():
            points[acq]=np.zeros(it)
            tar[acq]= np.zeros(it)
            regrets[acq]=np.zeros(it)
            utility_function[acq] = UtilityFunction(kind=acq, kappa = params[acq]['kappa'], xi = params[acq]['xi'])

        for i in range(it):
            for acq, optimizer in optimizers.items():
                points[acq][i]=optimizer.suggest(utility_function[acq])['x']
                tar[acq][i]= target(points[acq][i])
                optimizer.register(params=points[acq][i], target=tar[acq][i])
                regrets[acq][i]= simple_regret(f_max, optimizer._space.target)

    fig = plt.figure(figsize=(13, 6))

    af = plt.subplot()

    fig.suptitle(
        'Simple Regret After {} Iterations And {} Initial Points'.format(it, init_point),
        fontdict={'size':50}
    )

    steps = len(list(optimizers.values())[0].space)
    n = len(optimizers)
    colmap = plt.cm.get_cmap(name='rainbow', lut=n)

    num_iter=np.arange(1,it+1)

    for i, acq in enumerate(optimizers.keys()):
        af.plot(num_iter, regrets[acq], '*',markersize=15,markerfacecolor=colmap(i), markeredgecolor='k', markeredgewidth=1,label=acq,linestyle='solid',color=colmap(i))

    af.set_ylim((-0.2,0.4))
    af.axhline(y=0, linestyle=':', label='0.0')
    af.legend(loc=2, bbox_to_anchor=(1.01, 1), borderaxespad=0.)

    af.set_xticks(num_iter)

    plt.xlabel('Number of iterations')
    plt.ylabel('Simple regret')

    plt.show()

max_stars_repo_name: MartinaGaravaglia/BayesianOptimization

max_issues_repo_name: MartinaGaravaglia/BayesianOptimization

max_forks_repo_name: MartinaGaravaglia/BayesianOptimization


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# @Author: Brian Cherinka, José Sánchez-Gallego, and Brett Andrews
# @Date: 2017-11-03
# @Filename: spaxel.py
# @License: BSD 3-clause (http://www.opensource.org/licenses/BSD-3-Clause)
#
# @Last modified by: José Sánchez-Gallego (gallegoj@uw.edu)
# @Last modified time: 2018-11-08 18:00:34

from __future__ import absolute_import, division, print_function

import inspect
import itertools
import warnings

import numpy as np

import marvin
import marvin.core.exceptions
import marvin.core.marvin_pickle
import marvin.tools.cube
import marvin.tools.maps
import marvin.tools.modelcube
import marvin.utils.general.general
from marvin.core.exceptions import MarvinBreadCrumb, MarvinError, MarvinUserWarning
from marvin.utils.datamodel.dap import datamodel as dap_datamodel
from marvin.utils.datamodel.drp import datamodel as drp_datamodel
from marvin.utils.general.structs import FuzzyDict

breadcrumb = MarvinBreadCrumb()

class DataModel(object):
    """A single object that holds the DRP and DAP datamodel."""

    def __init__(self, release):

        self.drp = drp_datamodel[release]
        self.dap = dap_datamodel[release]

class Spaxel(object):
    """A base class that contains information about a spaxel.

    This class represents an spaxel with information from the reduced DRP
    spectrum, the DAP maps properties, and the model spectrum from the DAP
    logcube. A `.SpaxelBase` can be initialised with all or only part of that
    information, and either from a file, a database, or remotely via the
    Marvin API.

    The `~marvin.tools.cube.Cube`, `~marvin.tools.maps.Maps` , and
    `~marvin.tools.modelcube.ModelCube` quantities for the spaxel are available
    in ``cube_quantities``, ``maps_quantities``, and ``modelcube_quantities``,
    respectively. For convenience, the quantities can also be accessed directly
    from the `.SpaxelBase` itself (e.g., ``spaxel.emline_gflux_ha_6465``).

    Parameters:
        x,y (int):
            The `x` and `y` coordinates of the spaxel in the cube (0-indexed).
        cube (`~marvin.tools.cube.Cube` object or path or bool):
            If ``cube`` is a `~marvin.tools.cube.Cube` object, that
            cube will be used for the `.SpaxelBase` instantiation. This mode
            is mostly intended for `~marvin.utils.general.general.getSpaxel`
            as it significantly improves loading time. Otherwise, ``cube`` can
            be ``True`` (default), in which case a cube will be instantiated
            using the input ``filename``, ``mangaid``, or ``plateifu``. If
            ``cube=False``, no cube will be used and the cube associated
            quantities will not be available. ``cube`` can also be the
            path to the DRP cube to use.
        maps (`~marvin.tools.maps.Maps` object or path or bool)
            As ``cube`` but for the DAP measurements corresponding to the
            spaxel in the `.Maps`.
        modelcube (`marvin.tools.modelcube.ModelCube` object or path or bool)
            As ``maps`` but for the DAP measurements corresponding to the
            spaxel in the `.ModelCube`.
        lazy (bool):
            If ``False``, the spaxel data is loaded on instantiation.
            Otherwise, only the metadata is created. The associated quantities
            can be then loaded by calling `.SpaxelBase.load()`.
        kwargs (dict):
            Arguments to be passed to `.Cube`, `.Maps`, and `.ModelCube`
            when (and if) they are initialised.

    Attributes:
        cube_quantities (`~marvin.utils.general.structs.FuzzyDict`):
            A querable dictionary with the `.Spectrum` quantities
            derived from `.Cube` and matching ``x, y``.
        datamodel (object):
            An object containing the DRP and DAP datamodels.
        maps_quantities (`~marvin.utils.general.structs.FuzzyDict`):
            A querable dictionary with the `.AnalysisProperty` quantities
            derived from `.Maps` and matching ``x, y``.
        model_quantities (`~marvin.utils.general.structs.FuzzyDict`):
            A querable dictionary with the `.Spectrum` quantities
            derived from `.ModelCube` and matching ``x, y``.
        ra,dec (float):
            Right ascension and declination of the spaxel. Not available until
            the spaxel has been `loaded <.SpaxelBase.load>`.

    """

    def __init__(self, x, y, cube=True, maps=True, modelcube=True, lazy=False, **kwargs):

        if not cube and not maps and not modelcube:
            raise MarvinError('no inputs defined.')

        self.cube_quantities = FuzzyDict({})
        self.maps_quantities = FuzzyDict({})
        self.modelcube_quantities = FuzzyDict({})

        self._cube = cube
        self._maps = maps
        self._modelcube = modelcube

        for attr in ['mangaid', 'plateifu', 'release', 'bintype', 'template']:

            value = kwargs.pop(attr, None) or \
                getattr(cube, attr, None) or \
                getattr(maps, attr, None) or \
                getattr(modelcube, attr, None)

            setattr(self, attr, value)

        self._kwargs = kwargs
        self._parent_shape = None

        # drop breadcrumb
        breadcrumb.drop(message='Initializing MarvinSpaxel {0}'.format(self.__class__),
                        category=self.__class__)

        self.x = int(x)
        self.y = int(y)

        self.loaded = False
        self.datamodel = None

        if lazy is False:
            self.load()

        # Load VACs
        from marvin.contrib.vacs.base import VACMixIn
        self.vacs = VACMixIn.get_vacs(self)

    def __dir__(self):

        class_members = list(list(zip(*inspect.getmembers(self.__class__)))[0])
        instance_attr = list(self.__dict__.keys())

        items = self.cube_quantities.__dir__()
        items += self.maps_quantities.__dir__()
        items += self.modelcube_quantities.__dir__()
        items += class_members + instance_attr

        return sorted(items)

    def __getattr__(self, value):

        _getattr = super(Spaxel, self).__getattribute__

        for tool_quantity_dict in ['cube_quantities', 'maps_quantities', 'modelcube_quantities']:
            if value in _getattr(tool_quantity_dict):
                return _getattr(tool_quantity_dict)[value]

        return super(Spaxel, self).__getattribute__(value)

    def __repr__(self):
        """Spaxel representation."""

        if not self.loaded:
            return '<Marvin Spaxel (x={0.x:d}, y={0.y:d}, loaded=False)'.format(self)

        # Gets the coordinates relative to the centre of the cube/maps.
        y_mid, x_mid = np.array(self._parent_shape) / 2.
        x_centre = int(self.x - x_mid)
        y_centre = int(self.y - y_mid)

        # Determine what tools are loaded.
        tools = np.array(['cube', 'maps', 'modelcube'])
        load_idx = np.where([self._cube, self._maps, self._modelcube])[0]
        flags = '/'.join(tools[load_idx])

        return ('<Marvin Spaxel (plateifu={0.plateifu}, x={0.x:d}, y={0.y:d}; '
                'x_cen={1:d}, y_cen={2:d}, loaded={3})>'.format(self, x_centre, y_centre, flags))

    def _check_versions(self, attr):
        """Checks that all input object have the same versions.

        Runs sanity checks to make sure that ``attr`` has the same value
        in the input `.Cube`, `.Maps`, and `.ModelCube`. Returns the value
        for the attribute or ``None`` if the attribute does not exist.

        """

        out_value = None

        inputs = []
        for obj in [self._cube, self._maps, self._modelcube]:
            if obj is not None and not isinstance(obj, bool):
                inputs.append(obj)

        if len(inputs) == 1:
            return getattr(inputs[0], attr, None)

        for obj_a, obj_b in itertools.combinations(inputs, 2):
            if hasattr(obj_a, attr) and hasattr(obj_b, attr):
                assert getattr(obj_a, attr) == getattr(obj_b, attr), \
                    'inconsistent {!r} between {!r} and {!r}'.format(attr, obj_a, obj_b)
            out_value = getattr(obj_a, attr, None) or getattr(obj_b, attr, None)

        return out_value

    def _set_radec(self):
        """Calculates ra and dec for this spaxel."""

        self.ra = None
        self.dec = None

        for obj in [self._cube, self._maps, self._modelcube]:
            if hasattr(obj, 'wcs'):
                if obj.wcs.naxis == 2:
                    self.ra, self.dec = obj.wcs.wcs_pix2world([[self.x, self.y]], 0)[0]
                elif obj.wcs.naxis == 3:
                    self.ra, self.dec, __ = obj.wcs.wcs_pix2world([[self.x, self.y, 0]], 0)[0]

    def save(self, path, overwrite=False):
        """Pickles the spaxel to a file.

        Parameters:
            path (str):
                The path of the file to which the `.Spaxel` will be saved.
                Unlike for other Marvin Tools that derive from
                `~marvin.tools.core.MarvinToolsClass`, ``path`` is
                mandatory for `.Spaxel.save` as there is no default path for a
                given spaxel.
            overwrite (bool):
                If True, and the ``path`` already exists, overwrites it.
                Otherwise it will fail.

        Returns:
            path (str):
                The realpath to which the file has been saved.

        """

        return marvin.core.marvin_pickle.save(self, path=path, overwrite=overwrite)

    @classmethod
    def restore(cls, path, delete=False):
        """Restores a Spaxel object from a pickled file.

        If ``delete=True``, the pickled file will be removed after it has been
        unplickled. Note that, for objects with ``data_origin='file'``, the
        original file must exists and be in the same path as when the object
        was first created.

        """

        return marvin.core.marvin_pickle.restore(path, delete=delete)

    def load(self, force=None):
        """Loads the spaxel data.

        Loads the spaxel data for cubes/maps/modelcube. By default attempts
        to load whatever is specified when spaxels are instantianted from other
        Marvin Tools. Can manually force load a data type with the force
        keyword.

        Parameters:
        -----------
        force : {cube|maps|modelcube}
            Datatype to force load.

        """

        if self.loaded and force is None:
            warnings.warn('already loaded', MarvinUserWarning)
            return

        assert force in [None, 'cube', 'maps', 'modelcube'], \
            'force can only be cube, maps, or modelcube'

        for tool in ['cube', 'maps', 'modelcube']:
            self._load_tool(tool, force=(force is not None and force == tool))

        self._set_radec()
        self.loaded = True

        for attr in ['mangaid', 'plateifu', 'release', 'bintype', 'template']:
            setattr(self, attr, self._check_versions(attr))

        self.datamodel = DataModel(self.release)

    def _load_tool(self, tool, force=False):
        """Loads the tool and the associated quantities."""

        if tool == 'cube':
            class_name = marvin.tools.cube.Cube
            method = self.getCube
            quantities_dict = 'cube_quantities'
        elif tool == 'maps':
            class_name = marvin.tools.maps.Maps
            method = self.getMaps
            quantities_dict = 'maps_quantities'
        elif tool == 'modelcube':
            class_name = marvin.tools.modelcube.ModelCube
            method = self.getModelCube
            quantities_dict = 'modelcube_quantities'

        attr_value = getattr(self, '_' + tool)
        if (attr_value is False or attr_value is None) and force is False:
            setattr(self, '_' + tool, None)
            return

        if not isinstance(attr_value, class_name):
            if tool == 'modelcube' and self.release == 'MPL-4':
                warnings.warn('ModelCube cannot be instantiated for MPL-4.', MarvinUserWarning)
                self._modelcube = None
                return
            else:
                setattr(self, '_' + tool, method())
        else:
            if force is True:
                warnings.warn('{0} is already loaded'.format(tool), MarvinUserWarning)

        self._parent_shape = getattr(getattr(self, '_' + tool), '_shape')

        setattr(self, quantities_dict,
                getattr(getattr(self, '_' + tool), '_get_spaxel_quantities')(self.x, self.y,
                                                                             spaxel=self))

    def getCube(self):
        """Returns the associated `~marvin.tools.cube.Cube`"""

        if isinstance(self._cube, marvin.tools.cube.Cube):
            return self._cube

        cube_input = (self._cube if self._cube is not True else None) \
            or self.plateifu or self.mangaid

        return marvin.tools.cube.Cube(cube_input, release=self.release, **self._kwargs)

    def getMaps(self):
        """Returns the associated `~marvin.tools.maps.Maps`"""

        if isinstance(self._maps, marvin.tools.maps.Maps):
            return self._maps

        maps_input = (self._maps if self._maps is not True else None) \
            or self.plateifu or self.mangaid

        return marvin.tools.maps.Maps(maps_input, bintype=self.bintype,
                                      template=self.template, release=self.release,
                                      **self._kwargs)

    def getModelCube(self):
        """Returns the associated `~marvin.tools.modelcube.ModelCube`"""

        if isinstance(self._modelcube, marvin.tools.modelcube.ModelCube):
            return self._modelcube

        modelcube_input = (self._modelcube if self._modelcube is not True else None) \
            or self.plateifu or self.mangaid

        return marvin.tools.modelcube.ModelCube(modelcube_input, bintype=self.bintype,
                                                template=self.template, release=self.release,
                                                **self._kwargs)

    @property
    def quality_flags(self):
        """Bundle Cube DRP3QUAL and Maps DAPQUAL flags."""

        drp3qual = self.datamodel.drp.bitmasks['MANGA_DRP3QUAL']
        cube = self.getCube()
        drp3qual.mask = int(cube.header['DRP3QUAL'])
        qual_flags = [drp3qual]

        if self.release != 'MPL-4':
            qual_flags.append(self.datamodel.dap.bitmasks['MANGA_DAPQUAL'])

        return qual_flags

max_stars_repo_name: margudo/marvin

max_issues_repo_name: margudo/marvin

max_forks_repo_name: margudo/marvin