text:
#!/usr/bin/env python3

r'''
Backend script to plot drifting, narrowband events in a generalized cadence of
ON-OFF radio SETI observations. The main function contained in this file is
:func:`~.plot_candidate_events` uses the other helper functions
in this file (described below) to plot events from a turboSETI event .csv file.
'''

from os import mkdir
from os.path import dirname, abspath, isdir
import gc
import logging
logger_plot_event_name = 'plot_event'
logger_plot_event = logging.getLogger(logger_plot_event_name)
logger_plot_event.setLevel(logging.INFO)

# Plotting packages import
import matplotlib
import matplotlib.pyplot as plt
matplotlib.use('agg')

# Math/Science package imports
import numpy as np
from astropy.time import Time

# BL imports
import blimpy as bl
from blimpy.utils import rebin

# preliminary plot arguments
fontsize=16
font = {'family' : 'DejaVu Sans',
'size' : fontsize}
MAX_IMSHOW_POINTS = (4096, 1268)

def overlay_drift(f_event, f_start, f_stop, drift_rate, t_duration, offset=0, alpha=1, color='#cc0000'):
    r'''
    Creates a dashed red line at the recorded frequency and drift rate of
    the plotted event - can overlay the signal exactly or be offset by
    some amount (offset can be 0 or 'auto').
    '''
    # determines automatic offset and plots offset lines
    if offset == 'auto':
        offset = ((f_start - f_stop) / 10)
        plt.plot((f_event - offset, f_event),
                 (10, 10),
                 "o-",
                 c=color,
                 lw=2,
                 alpha=alpha)

    # plots drift overlay line, with offset if desired
    plt.plot((f_event + offset, f_event + drift_rate/1e6 * t_duration + offset),
             (0, t_duration),
             c=color,
             ls='dashed', lw=2,
             alpha=alpha)

def plot_waterfall(wf, source_name, f_start=None, f_stop=None, **kwargs):
    r"""
    Plot waterfall of data in a .fil or .h5 file.

    Parameters
    ----------
    wf : blimpy.Waterfall object
        Waterfall object of an H5 or Filterbank file containing the dynamic spectrum data.
    source_name : str
        Name of the target.
    f_start : float
        Start frequency, in MHz.
    f_stop : float
        Stop frequency, in MHz.
    kwargs : dict
        Keyword args to be passed to matplotlib imshow().

    Notes
    -----
    Plot a single-panel waterfall plot (frequency vs. time vs. intensity)
    for one of the on or off observations in the cadence of interest, at the
    frequency of the expected event. Calls :func:`~overlay_drift`

    """
    # prepare font
    matplotlib.rc('font', **font)

    # Load in the data from fil
    plot_f, plot_data = wf.grab_data(f_start=f_start, f_stop=f_stop)

    # Make sure waterfall plot is under 4k*4k
    dec_fac_x, dec_fac_y = 1, 1

    # rebinning data to plot correctly with fewer points
    try:
        if plot_data.shape[0] > MAX_IMSHOW_POINTS[0]:
            dec_fac_x = plot_data.shape[0] / MAX_IMSHOW_POINTS[0]
        if plot_data.shape[1] > MAX_IMSHOW_POINTS[1]:
            dec_fac_y =  int(np.ceil(plot_data.shape[1] /  MAX_IMSHOW_POINTS[1]))
        plot_data = rebin(plot_data, dec_fac_x, dec_fac_y)
    except Exception as ex:
        print('\n*** Oops, grab_data returned plot_data.shape={}, plot_f.shape={}'
              .format(plot_data.shape, plot_f.shape))
        print('Waterfall info for {}:'.format(wf.filename))
        wf.info()
        raise ValueError('*** Something is wrong with the grab_data output!') from ex

    # Rolled back PR #82

    # determine extent of the plotting panel for imshow
    extent=(plot_f[0], plot_f[-1], (wf.timestamps[-1]-wf.timestamps[0])*24.*60.*60, 0.0)

    # plot and scale intensity (log vs. linear)
    kwargs['cmap'] = kwargs.get('cmap', 'viridis')
    plot_data = 10.0 * np.log10(plot_data)

    # get normalization parameters
    vmin = plot_data.min()
    vmax = plot_data.max()
    normalized_plot_data = (plot_data - vmin) / (vmax - vmin)

    # display the waterfall plot
    this_plot = plt.imshow(normalized_plot_data,
        aspect='auto',
        rasterized=True,
        interpolation='nearest',
        extent=extent,
        **kwargs
    )

    # add plot labels
    plt.xlabel("Frequency [Hz]",fontdict=font)
    plt.ylabel("Time [s]",fontdict=font)

    # add source name
    ax = plt.gca()
    plt.text(0.03, 0.8, source_name, transform=ax.transAxes, bbox=dict(facecolor='white'))
    # if plot_snr != False:
    #     plt.text(0.03, 0.6, plot_snr, transform=ax.transAxes, bbox=dict(facecolor='white'))
    # return plot

    del plot_f, plot_data
    gc.collect()

    return this_plot

def make_waterfall_plots(fil_file_list, on_source_name, f_start, f_stop, drift_rate, f_mid,
                         filter_level, source_name_list, offset=0, plot_dir=None, **kwargs):
    r'''
    Makes waterfall plots of an event for an entire on-off cadence.

    Parameters
    ----------
    fil_file_list : str
        List of filterbank files in the cadence.
    on_source_name : str
        Name of the on_source target.
    f_start : float
        Start frequency, in MHz.
    f_stop : float
        Stop frequency, in MHz.
    drift_rate : float
        Drift rate in Hz/s.
    f_mid : float
        <iddle frequency of the event, in MHz.
    filter_level : int
        Filter level (1, 2, or 3) that produced the event.
    source_name_list : list
        List of source names in the cadence, in order.
    bandwidth : int
        Width of the plot, incorporating drift info.
    kwargs : dict
        Keyword args to be passed to matplotlib imshow().

    Notes
    -----
    Makes a series of waterfall plots, to be read from top to bottom, displaying a full cadence
    at the frequency of a recorded event from find_event. Calls :func:`~plot_waterfall`

    '''
    global logger_plot_event

    # prepare for plotting
    matplotlib.rc('font', **font)

    # set up the sub-plots
    n_plots = len(fil_file_list)
    fig = plt.subplots(n_plots, sharex=True, sharey=True,figsize=(10, 2*n_plots))

    # get directory path for storing PNG files
    if plot_dir is None:
        dirpath = dirname(abspath(fil_file_list[0])) + '/'
    else:
        if not isdir(plot_dir):
            mkdir(plot_dir)
        dirpath = plot_dir

    # read in data for the first panel
    max_load = bl.calcload.calc_max_load(fil_file_list[0])
    #print('plot_event make_waterfall_plots: max_load={} is required for {}'.format(max_load, fil_file_list[0]))
    wf1 = bl.Waterfall(fil_file_list[0], f_start=f_start, f_stop=f_stop, max_load=max_load)
    t0 = wf1.header['tstart']
    plot_f1, plot_data1 = wf1.grab_data()

    # rebin data to plot correctly with fewer points
    dec_fac_x, dec_fac_y = 1, 1
    if plot_data1.shape[0] > MAX_IMSHOW_POINTS[0]:
        dec_fac_x = plot_data1.shape[0] / MAX_IMSHOW_POINTS[0]
    if plot_data1.shape[1] > MAX_IMSHOW_POINTS[1]:
        dec_fac_y =  int(np.ceil(plot_data1.shape[1] /  MAX_IMSHOW_POINTS[1]))
    plot_data1 = rebin(plot_data1, dec_fac_x, dec_fac_y)

    # define more plot parameters
    # never used: delta_f = 0.000250
    mid_f = np.abs(f_start+f_stop)/2.

    subplots = []
    del wf1, plot_f1, plot_data1
    gc.collect()

    # Fill in each subplot for the full plot
    for ii, filename in enumerate(fil_file_list):
        logger_plot_event.debug('make_waterfall_plots: file {} in list: {}'.format(ii, filename))
        # identify panel
        subplot = plt.subplot(n_plots, 1, ii + 1)
        subplots.append(subplot)

        # read in data
        max_load = bl.calcload.calc_max_load(filename)
        #print('plot_event make_waterfall_plots: max_load={} is required for {}'.format(max_load, filename))
        wf = bl.Waterfall(filename, f_start=f_start, f_stop=f_stop, max_load=max_load)
        # make plot with plot_waterfall
        source_name = source_name_list[ii]
        this_plot = plot_waterfall(wf,
                                   source_name,
                                   f_start=f_start,
                                   f_stop=f_stop,
                                   **kwargs)

        # calculate parameters for estimated drift line
        t_elapsed = Time(wf.header['tstart'], format='mjd').unix - Time(t0, format='mjd').unix
        t_duration = (wf.n_ints_in_file - 1) * wf.header['tsamp']
        f_event = f_mid + drift_rate / 1e6 * t_elapsed

        # plot estimated drift line
        overlay_drift(f_event, f_start, f_stop, drift_rate, t_duration, offset)

        # Title the full plot
        if ii == 0:
            plot_title = "%s \n $\\dot{\\nu}$ = %2.3f Hz/s, MJD:%5.5f" % (on_source_name, drift_rate, t0)
            plt.title(plot_title)
        # Format full plot
        if ii < len(fil_file_list)-1:
            plt.xticks(np.linspace(f_start, f_stop, num=4), ['','','',''])

        del wf
        gc.collect()

    # More overall plot formatting, axis labelling
    factor = 1e6
    units = 'Hz'

    #ax = plt.gca()
    #ax.get_xaxis().get_major_formatter().set_useOffset(False)
    xloc = np.linspace(f_start, f_stop, 5)
    xticks = [round(loc_freq) for loc_freq in (xloc - mid_f)*factor]
    if np.max(xticks) > 1000:
        xticks = [xt/1000 for xt in xticks]
        units = 'kHz'
    plt.xticks(xloc, xticks)
    plt.xlabel("Relative Frequency [%s] from %f MHz"%(units,mid_f),fontdict=font)

    # Add colorbar
    cax = fig[0].add_axes([0.94, 0.11, 0.03, 0.77])
    fig[0].colorbar(this_plot,cax=cax,label='Normalized Power (Arbitrary Units)')

    # Adjust plots
    plt.subplots_adjust(hspace=0,wspace=0)

    # save the figures
    path_png = dirpath + str(filter_level) + '_' + on_source_name + '_dr_' + "{:0.2f}".format(drift_rate) + '_freq_' "{:0.6f}".format(f_start) + ".png"
    plt.savefig(path_png, bbox_inches='tight')
    logger_plot_event.debug('make_waterfall_plots: Saved file {}'.format(path_png))

    # show figure before closing if this is an interactive context
    mplbe = matplotlib.get_backend()
    logger_plot_event.debug('make_waterfall_plots: backend = {}'.format(mplbe))
    if mplbe != 'agg':
        plt.show()

    # close all figure windows
    plt.close('all')

    return subplots

def plot_candidate_events(candidate_event_dataframe, fil_file_list, filter_level, source_name_list,
                          offset=0, plot_dir=None, **kwargs):
    r'''
    Calls :func:`~make_waterfall_plots` on each event in the input .csv file.

    Arguments
    ---------
    candidate_event_dataframe : dict
        A pandas dataframe containing information
        about a candidate event. The necessary data
        includes the start and stop frequencies, the
        drift rate, and the source name. To determine
        the required variable names and formatting
        conventions, see the output of
        find_event_pipeline.
    fil_file_list : list
        A Python list that contains a series of
        strings corresponding to the filenames of .fil
        files, each on a new line, that corresponds to
        the cadence used to create the .csv file used
        for event_csv_string.
    filter_level : int
        A string indicating the filter level of the
        cadence used to generate the
        candidate_event_dataframe. Used only for
        output file naming, convention is "f1", "f2",
        or "f3". Descriptions for the three levels of
        filtering can be found in the documentation
        for find_event.py
    source_name_list : list
        A Python list that contains a series of strings
        corresponding to the source names of the
        cadence in chronological (descending through
        the plot panels) cadence.
    offset : int, optional
        The amount that the overdrawn "best guess"
        line from the event parameters in the csv
        should be shifted from its original position
        to enhance readability. Can be set to 0
        (default; draws line on top of estimated
        event) or 'auto' (shifts line to the left by
        an auto-calculated amount, with addition lines
        showing original position).
    kwargs : dict

    Examples
    --------
    It is highly recommended that users interact with this program via the
    front-facing plot_event_pipeline.py script. See the usage of that file in
    its own documentation.

    If you would like to run plot_candidate_events without calling
    plot_event_pipeline.py, the usage is as follows:

    >>> plot_event.plot_candidate_events(candidate_event_dataframe, fil_file_list,
    ...                                  filter_level, source_name_list, offset=0)

    '''
    global logger_plot_event

    # load in the data for each individual hit
    if candidate_event_dataframe is None:
        print('*** plot_candidate_events: candidate_event_dataframe is None, nothing to do.')
        return
    len_df = len(candidate_event_dataframe)
    if len_df < 1:
        print('*** plot_candidate_events: len(candidate_event_dataframe) = 0, nothing to do.')
        return
    for i in range(0, len_df):
        candidate = candidate_event_dataframe.iloc[i]
        on_source_name = candidate['Source']
        f_mid = candidate['Freq']
        drift_rate = candidate['DriftRate']

        # calculate the length of the total cadence from the fil files' headers
        first_fil = bl.Waterfall(fil_file_list[0], load_data=False)
        tfirst = first_fil.header['tstart']
        last_fil = bl.Waterfall(fil_file_list[-1], load_data=False)
        tlast = last_fil.header['tstart']
        t_elapsed = Time(tlast, format='mjd').unix - Time(tfirst, format='mjd').unix + (last_fil.n_ints_in_file -1) * last_fil.header['tsamp']

        # calculate the width of the plot based on making sure the full drift is visible
        bandwidth = 2.4 * abs(drift_rate)/1e6 * t_elapsed
        bandwidth = np.max((bandwidth, 500./1e6))

        # Get start and stop frequencies based on midpoint and bandwidth
        f_start, f_stop = np.sort((f_mid - (bandwidth/2),  f_mid + (bandwidth/2)))

        # logger_plot_event.debug useful values
        logger_plot_event.debug('*************************************************')
        logger_plot_event.debug('***     The Parameters for This Plot Are:    ****')
        logger_plot_event.debug('Target = {}'.format(on_source_name))
        logger_plot_event.debug('Bandwidth = {} MHz'.format(round(bandwidth, 5)))
        logger_plot_event.debug('Time Elapsed (inc. Slew) = {} s'.format(round(t_elapsed)))
        logger_plot_event.debug('Middle Frequency = {} MHz'.format(round(f_mid, 4)))
        logger_plot_event.debug('Expected Drift = {} Hz/s'.format(round(drift_rate, 4)))
        logger_plot_event.debug('*************************************************')

        # Pass info to make_waterfall_plots() function
        make_waterfall_plots(fil_file_list,
                             on_source_name,
                             f_start,
                             f_stop,
                             drift_rate,
                             f_mid,
                             filter_level,
                             source_name_list,
                             offset=offset,
                             plot_dir=plot_dir,
                             **kwargs)

ext: py

reason: import numpy,from astropy


text:
# coding:utf-8
from captcha.image import ImageCaptcha
import random
from PIL import Image
import numpy as np
import tensorflow as tf

number = ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9']
alphabet = ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', 'p', 'q', 'r', 's', 't', 'u',
            'v', 'w', 'x', 'y', 'z']
ALPHABET = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q', 'R', 'S', 'T', 'U',
            'V', 'W', 'X', 'Y', 'Z']
SAVE_PATH = "./"
CHAR_SET = number + alphabet + ALPHABET
CHAR_SET_LEN = len(CHAR_SET)
IMAGE_HEIGHT = 60
IMAGE_WIDTH = 160

def random_captcha_text(char_set=None, captcha_size=4):
    if char_set is None:
        char_set = number + alphabet + ALPHABET

    captcha_text = []
    for i in range(captcha_size):
        c = random.choice(char_set)
        captcha_text.append(c)
    return captcha_text

def gen_captcha_text_and_image(width=160, height=60, char_set=CHAR_SET):
    image = ImageCaptcha(width=width, height=height)

    captcha_text = random_captcha_text(char_set)
    captcha_text = ''.join(captcha_text)

    captcha = image.generate(captcha_text)

    captcha_image = Image.open(captcha)
    captcha_image = np.array(captcha_image)
    return captcha_text, captcha_image

text, image = gen_captcha_text_and_image(char_set=CHAR_SET)
MAX_CAPTCHA = len(text)
print('CHAR_SET_LEN=', CHAR_SET_LEN, ' MAX_CAPTCHA=', MAX_CAPTCHA)

def convert2gray(img):
    if len(img.shape) > 2:
        gray = np.mean(img, -1)
        return gray
    else:
        return img

def text2vec(text):
    vector = np.zeros([MAX_CAPTCHA, CHAR_SET_LEN])
    for i, c in enumerate(text):
        idx = CHAR_SET.index(c)
        vector[i][idx] = 1.0
    return vector

def vec2text(vec):
    text = []
    for i, c in enumerate(vec):
        text.append(CHAR_SET[c])
    return "".join(text)

def get_next_batch(batch_size=128):
    batch_x = np.zeros([batch_size, IMAGE_HEIGHT, IMAGE_WIDTH, 1])
    batch_y = np.zeros([batch_size, MAX_CAPTCHA, CHAR_SET_LEN])

    def wrap_gen_captcha_text_and_image():
        while True:
            text, image = gen_captcha_text_and_image(char_set=CHAR_SET)
            if image.shape == (60, 160, 3):
                return text, image

    for i in range(batch_size):
        text, image = wrap_gen_captcha_text_and_image()
        image = tf.reshape(convert2gray(image), (IMAGE_HEIGHT, IMAGE_WIDTH, 1))
        batch_x[i, :] = image
        batch_y[i, :] = text2vec(text)

    return batch_x, batch_y

def crack_captcha_cnn():
    model = tf.keras.Sequential()

    model.add(tf.keras.layers.Conv2D(32, (3, 3)))
    model.add(tf.keras.layers.PReLU())
    model.add(tf.keras.layers.MaxPool2D((2, 2), strides=2))

    model.add(tf.keras.layers.Conv2D(64, (5, 5)))
    model.add(tf.keras.layers.PReLU())
    model.add(tf.keras.layers.MaxPool2D((2, 2), strides=2))

    model.add(tf.keras.layers.Conv2D(128, (5, 5)))
    model.add(tf.keras.layers.PReLU())
    model.add(tf.keras.layers.MaxPool2D((2, 2), strides=2))

    model.add(tf.keras.layers.Flatten())
    model.add(tf.keras.layers.Dense(MAX_CAPTCHA * CHAR_SET_LEN))
    model.add(tf.keras.layers.Reshape([MAX_CAPTCHA, CHAR_SET_LEN]))
    model.add(tf.keras.layers.Softmax())

    return model

def train():
    try:
        model = tf.keras.models.load_model(SAVE_PATH + 'model')
    except Exception as e:
        print('#######Exception', e)
        model = crack_captcha_cnn()

    model.compile(optimizer='Adam',
                  metrics=['accuracy'],
                  loss='categorical_crossentropy')

    for times in range(5):
        batch_x, batch_y = get_next_batch(512)
        print('times=', times, ' batch_x.shape=', batch_x.shape, ' batch_y.shape=', batch_y.shape)
        model.fit(batch_x, batch_y, epochs=4)
        print("y预测=\n", np.argmax(model.predict(batch_x), axis=2))
        print("y实际=\n", np.argmax(batch_y, axis=2))

        if 0 == times % 10:
            print("save model at times=", times)
            model.save(SAVE_PATH + 'model')

def predict():
    model = tf.keras.models.load_model(SAVE_PATH + 'model')
    success = 0
    count = 100
    for _ in range(count):
        data_x, data_y = get_next_batch(1)
        prediction_value = model.predict(data_x)
        data_y = vec2text(np.argmax(data_y, axis=2)[0])
        prediction_value = vec2text(np.argmax(prediction_value, axis=2)[0])

        if data_y.upper() == prediction_value.upper():
            print("y预测=", prediction_value, "y实际=", data_y, "预测成功。")
            success += 1
        else:
            print("y预测=", prediction_value, "y实际=", data_y, "预测失败。")

    print("预测", count, "次", "成功率 =", success / count)

    pass

if __name__ == "__main__":
    train()
    predict()

reason: import numpy


text:
import bs4 as bs
import datetime as dt
import matplotlib.pyplot as plt
from matplotlib import style
import numpy as np
import os
import pandas as pd
import pandas_datareader.data as web
import pickle
import requests
from collections import Counter
style.use('ggplot')

################################################
################################################get stock data
def save_100_tickers():
    resp = requests.get('https://finance.yahoo.com/most-active?offset=0&count=100')
    soup = bs.BeautifulSoup(resp.text, 'lxml')
    table = soup.find('table', {'class': 'W(100%)'})
    tickers=[]
    for i in table.findAll('tr')[1:]:
        deta=i.findAll('td')[0].text
        tickers.append(deta)
    with open("100tickers.pickle", "wb") as f:
        pickle.dump(tickers, f)
    return tickers

def get_stock_from_yahoo():

    with open("100tickers.pickle", "rb") as f:
        tickers = pickle.load(f)
    if not os.path.exists('100stock'):
        os.makedirs('100stock')

    start = dt.datetime(2019, 1, 1)
    end = dt.datetime.now()
    for ticker in tickers:
        # just in case your connection breaks, we'd like to save our progress!
        if not os.path.exists('100stock/{}.csv'.format(ticker)):
            df = web.DataReader(ticker, 'yahoo', start, end)
            df.reset_index(inplace=True)
            df.set_index("Date", inplace=True)
            #df = df.drop("Symbol", axis=1)
            df.to_csv('100stock/{}.csv'.format(ticker))
        else:
            print('Already have {}'.format(ticker))

def stock_close_data():
    with open("100tickers.pickle", "rb") as f:
        tickers = pickle.load(f)

    main_df = pd.DataFrame()
    for count, ticker in enumerate(tickers):
        df = pd.read_csv('100stock/{}.csv'.format(ticker))
        df.set_index('Date', inplace=True)
        df.rename(columns={'Adj Close': ticker}, inplace=True)
        df.drop(['Open', 'High', 'Low', 'Close', 'Volume'], 1, inplace=True)

        if main_df.empty:
            main_df = df
        else:
            main_df = main_df.join(df, how='outer')
        main_df.fillna(0,inplace=True)
    main_df.to_csv('100stock_closes.csv')
################################################
################################################

def visualize_data():
    df = pd.read_csv('100stock_closes.csv')
    df_corr = df.corr()
    #print(df_corr.head())
    df_corr.to_csv('100corr.csv')
    data1 = df_corr.values
    fig1 = plt.figure()
    ax1 = fig1.add_subplot(111)

    heatmap1 = ax1.pcolor(data1, cmap=plt.cm.RdYlGn)
    fig1.colorbar(heatmap1)

    ax1.set_xticks(np.arange(data1.shape[1]) + 0.5, minor=False)
    ax1.set_yticks(np.arange(data1.shape[0]) + 0.5, minor=False)
    # ax1.invert_yaxis()
    # ax1.xaxis.tick_top()
    # column_labels = df_corr.columns
    # row_labels = df_corr.index
    # ax1.set_xticklabels(column_labels)
    # ax1.set_yticklabels(row_labels)
    # plt.xticks(rotation=90)
    # heatmap1.set_clim(-1, 1)
    plt.tight_layout()
    plt.show()

def next7days_stock():
    hm_days = 7
    df = pd.read_csv('100stock_closes.csv', index_col=0)
    tickers = df.columns.values.tolist()
    for m in tickers:
        for i in range(1, hm_days+1):
            df['{}_{}d'.format(m, i)] = (df[m].shift(-i) - df[m]) / df[m]
            df.fillna(0, inplace=True)
            daf=df.iloc[:,100:]
    daf.to_csv('next7days_stock.csv')
    return tickers, daf

def buy_sell_hold(*args):
    cols = [c for c in args]
    requirement = 0.05
    for col in cols:
        if col > requirement:
            return 1
        if col < -requirement:
            return -1
    return 0

'''This will let us see the distributions of classes both in our dataset and in our algorithm's predictions'''

##Couter:dict subclass for counting hashable objects
def extract_featuresets(tickers):
    #tickers, daf = next7days_stock()
    #vals=[]
    for n in tickers:
        daf['{}_target'.format(n)] = list(map( buy_sell_hold,
                                               daf['{}_1d'.format(n)],
                                               daf['{}_2d'.format(n)],
                                               daf['{}_3d'.format(n)],
                                               daf['{}_4d'.format(n)],
                                               daf['{}_5d'.format(n)],
                                               daf['{}_6d'.format(n)],
                                               daf['{}_7d'.format(n)] ))
        vals = df['{}_target'.format(ticker)].values.tolist()
        str_vals = [str(i) for i in vals]
        #print('Data spread:', Counter(str_vals))

        df.fillna(0, inplace=True)
        df = df.replace([np.inf, -np.inf], np.nan)
        df.dropna(inplace=True)

        df_vals = df[[ticker for ticker in tickers]].pct_change()
        df_vals = df_vals.replace([np.inf, -np.inf], 0)
        df_vals.fillna(0, inplace=True)

        X = df_vals.values
        y = df['{}_target'.format(ticker)].values
    return X, y, df
    target_data=[daf['{}_target'.format(x)] for x in tickers]

from sklearn import svm, neighbors
from sklearn.model_selection import cross_validate
from sklearn.ensemble import VotingClassifier, RandomForestClassifier

def do_ml(ticker):

    X,y,df = extract_featuresets(tickers)
    X_train, X_test, y_train, y_test = cross_validate.train_test_split(X,
                                                        y,
                                                        test_size=0.25)


text:
import numpy as np
import json
import random
import gc

from keras.models import load_model
from keras import backend as K
from model import fm_model
from data_stuff import feature_scaling, normalize_output, split_data
from csv_stuff import get_permutation_csv, write_csv

def round_1(x_train, y_train, param, write = 'result/round_1.csv'):
    rand_p = rand_params(param)

    history, model = fm_model(x_train, y_train, rand_p)
    train_result = model.evaluate(x_train, y_train,verbose=0)
    validation_result = [history.history['val_loss'][-1] ,history.history['val_mean_absolute_error'][-1]]

    jsonparam = json.dumps(rand_p)
    row = [train_result[0], train_result[1], validation_result[0], validation_result[1], jsonparam]
    write_csv(write, row)

def round_2(x_train, y_train, x_test, y_test, param_no, fold=10, read = 'result/round_1.csv', write = 'result/round_2.csv'):
    minp={
        'hidden_layers':1,
        'l1':1,
        'l2':1,
        'l3':1,
        'l4':1,
        'l5':1,
        'ac1':'relu',
        'ac2':'relu',
        'ac3':'relu',
        'ac4':'relu',
        'ac5':'relu',
        'd1':0,
        'd2':0,
        'd3':0,
        'd4':0,
        'd5':0,
        'kernel_init':'random_uniform',
        'loss_func':'mse',
        'metrics':'mae',
        'val_split':0.2,
        'epoch_size':2,
        'batch_size':1000,
        'lr':1
    }
    permutation = get_permutation_csv(read)
    param = permutation[param_no]['param']
    vloss = []
    vmetric = []
    loss = []
    metric = []
    highscore = 999
    _,highmodel = fm_model(x_train, y_train, minp)

    for i in range(fold):
        print('==================== round_2: validating {} ===================='.format(i+1))
        history, model = fm_model(x_train, y_train, param)
        train_result = model.evaluate(x_train, y_train,verbose=0)

        loss.append(train_result[0])
        metric.append(train_result[1])
        vloss.append(history.history['val_loss'][-1])
        vmetric.append(history.history['val_mean_absolute_error'][-1])

        if vmetric[-1] < highscore:
            highscore = vmetric[-1]
            high_model = model

    vloss_avg = sum(vloss)/len(vloss)
    vmetric_avg = sum(vmetric)/len(vmetric)
    loss_avg = sum(loss)/len(loss)
    metric_avg = sum(metric)/len(metric)

    jsonparam = json.dumps(param)
    savedir = 'result/model_{}.h5'.format(param_no)
    highmodel.save(savedir)
    # jsonmodel = highmodel.to_json()
    # jsonweight = highmodel.get_weights()
    test_score = model.evaluate(x_test, y_test)
    row = [loss_avg, metric_avg, vloss_avg, vmetric_avg ,jsonparam, test_score[1]]
    write_csv(write, row)

    gc.collect()
    del model
    K.clear_session()

def rand_params(p):
    rand_p = p
    for items in rand_p:
        next = rand_p[items]
        if isinstance(next, tuple):
            if items not in ['d1', 'd2', 'd3', 'd4', 'd5', 'lr']:
                next = int(random.randrange(next[0], next[1] ,next[2])/next[3])
            else:
                next = random.randrange(next[0], next[1] ,next[2])/next[3]
        if isinstance(next, list):
            next = random.choice(next)
        rand_p[items] = next
    return rand_p

max_stars_repo_path: Neural Network/Task2/param.py

max_issues_repo_path: Neural Network/Task2/param.py

max_forks_repo_path: Neural Network/Task2/param.py


text:
"""

                            mfp5io

    implements a hdf5/h5py based class to handle storage of
    molecular systems (used for restart and trajectory data)

    This is a rewrite of the V1.0 that ships with molsys and is based on a
    molsys type mol object. It is meant to work with any MM engine that provides
    a molsys mol object and provides a certain API. It can be used with pydlpoly
    (using the new startup mode with molsys) or pylmps currently. Note that in pylmps
    trajectory info is written in C++ in lammps using dump_mfp5.cpp (which needs to be
    compiled in from the package PYLMPS)

    one part of the data is "fixed"
    - number of atoms
    - elements
    - atom types
    - connectivity
    - boundary conditions
    - forcefield parameters
    in a restart this data is used as a replacement to the usual file input

    the other part is "staged"
    in other words if needed multiple stages can be defined like "equilibration",
    "sample1", "sample2" etc.
    the initial and default stage is called "default" by default ... yes!
    for each stage current restart information (position/velocitites/cell) is
    automatically stored every N steps. restart data is stored in double precision
    in addition, trajectory info is stored every M steps. (both N and M can be changed)
    trajectory data (pos or pos/vel or pos/vel/force) can be stored both in single and
    double precision on request.

    Versions (in file attribute "version"
    - no version info == 1.0 : original
    - version 1.1            : nstep of trajectories is stored per dataset and not for the stage
    - version 1.2            : mfp5io version to work with new mol objects (can write par info)

This has been renamed in 2021 from pdlpio2 to mfp5io (according with the name change in lammps)
"""

from h5py import File, special_dtype
import numpy as np

import types
import os
import string

from molsys import mpiobject
import molsys

class mfp5io(mpiobject):

    def __init__(self, fname, ffe=None, restart=None, filemode="a", mpi_comm=None, out=None):
        """generates a mfp5 file object

        This object knows three states: conencted to a ffe (force field engine) (with a mol object),
        connected to a ffe but without mol object, which implies a restart (will generate mol from mfp5 file)
        or the analysis mode. In the analysis mode writing to trajectory data is forbidden.
        In the connected mode one can request appending to a stage.
        Otherwise a new stage is required.

        NOTE: always use self.open() as a first command in a method that will access the file. it will check if the file
              is already open and only open if it is not open. in certain cases (e.g. dumping with lammps) the file needs to be
              closed before passing to lammps. To make sure the file is connected call open().

        Args:
            fname (string): filename (including .mfp5 extension)
            ffe (object, optional): Defaults to None. Force Field engine (must contain a ffe.mol instance of type molsys)
            restart (string, optional): Defaults to None. Name of the stage from which a restart should be read
        """
        super(mfp5io, self).__init__(mpi_comm,out)
        self.verbose = 0
        self.filemode = filemode
        # helper object for hdf5 variable length strings
        self.str_dt = special_dtype(vlen=str)
        #
        self.fname = fname
        # check if this file exists or should be generated
        self.fexists = os.path.isfile(self.fname)
        self.h5file_isopen = False
        # check if we run with an FFE or not
        self.ffe  = ffe
        if ffe is not None:
            self.ffe_name = type(ffe).__name__
            if restart is None:
                # a force field engine is connected => check if mol object is available and of the right type
                try:
                    self.mol = ffe.mol
                except AttributeError:
                    self.pprint("PDLP ERROR: Force Field Engine has no mol object!")
                    raise
                # ok, now let us check if mol is actually a molsys object (require it to have a is_topo attribute)
                try:
                    is_topo = self.mol.is_topo
                except AttributeError:
                    self.pprint("PDLP ERROR: mol object seems not be a molsys object! This is required for mfp5io")
                    raise
                # last but not least check that mol is not a topo ... we do not want that
                if is_topo is True:
                    self.pprint("PDLP ERROR: mol object is a topo object .. that should not happen!")
                    raise Exception("mfp5 error")
                self.mode = "ffe"
                # now all is fine .. set defaults for this runmode
                self.stage = "default" # default stage which can be overwritten (contains no trajectory data, only restarts)
            else:
                # this is a restart, which means the file must exist and also the stage name
                assert self.fexists, "PDLP ERROR: For restart the file must exist. Check filename and path!"
                self.restart_stage = restart
                self.mode = "restart"
        else:
            # if ffe is None we are in analysis mode (means the file must exist)
            assert self.fexists, "PDLP ERROR: For analysis the file must exist!"
            self.mode = "analysis"
            self.restart_stage = "default" # use default stage to start up (should be always there, right?)
        # now open the file
        self.open()
        #
        self.file_version = 1.2
        if self.fexists is True:
            if self.mode == "ffe" or self.mode == "restart":
                # this is an exisiting file and we are in ffe/restart mode and want to add data => make sure all is matching
                if self.is_master:
                    file_version = self.h5file.attrs["version"]
                else:
                    file_version = None
                file_version = self.mpi_comm.bcast(file_version)
                assert self.file_version == file_version, "PDLP ERROR: Exisiting file has a different version! Can not add data"
                if self.mode == "ffe":
                    # check the system if it contains the same molecule (we could be more clever here, but if the list of atomtypes is identical we should be safe)
                    # TBI mor consistent tests like paramters etc .. we could read the complete mol out and compare
                    self.compare_system()
        else:
            # file is new .. set version
            if self.is_master:
                self.h5file.attrs["version"] = self.file_version
            # this must be ffe mode so we can start generating the system group
            self.write_system()
            # new file so we write the initital structure to the restart
            self.add_stage(self.stage)
            self.write_restart()
        return

    # def __del__(self):
    #     self.close()
    #     return

    def open(self):
        """ method to make sure that h5file is open """
        if self.h5file_isopen:
            return
        self.h5file_isopen = True
        if self.is_master:
            self.h5file = File(self.fname, self.filemode)
        else:
            self.h5file = None
        return

    def close(self):
        """ method to make sure that the file is closed """
        if not self.h5file_isopen:
            return
        self.h5file_isopen = False
        if self.is_master:
            self.h5file.close()
        return

    # system IO
    #
    # The system group contains all the time independent information of the system (connectivity, atomtypes, parameter)
    # Note: in the old pdlpio system this data was prepared by the mol (assign_FF) class and just written here
    #       Now in mfpfio we expect a certain interface, namely the mol object from molsys and collect all the
    #       data in an active way

    def write_system(self):
        """ Writes all the info from a mol object to the system group in the mfp5 file

        The following info is written:
            - elems       (strings)
            - atypes      (strings)
            - fragtypes   (strings)
            - fragnumbers (int)
            - bcd         (int)

        Further info is stored in subgroups for different addons
            - ff addon (parameter and rics)
            - molecules addon (TBI)

        """
        self.open()
        if self.is_master:
            system = self.h5file.require_group("system")
            # elems
            na = self.ffe.mol.get_natoms()
            elems = self.ffe.mol.get_elems()
            mfp5_elems = system.require_dataset("elems", shape=(na,), dtype=self.str_dt)
            mfp5_elems[...] = elems
            # atypes
            atypes = self.ffe.mol.get_atypes()
            mfp5_atypes = system.require_dataset("atypes", shape=(na,), dtype=self.str_dt)
            mfp5_atypes[...] = atypes
            # fragtypes
            fragtypes = self.ffe.mol.get_fragtypes()
            mfp5_fragtypes = system.require_dataset("fragtypes", shape=(na,), dtype=self.str_dt)
            mfp5_fragtypes[...] = fragtypes
            # fragnumbers
            fragnumbers = self.ffe.mol.get_fragnumbers()
            mfp5_fragnumbers = system.require_dataset("fragnumbers", shape=(na,), dtype="i")
            mfp5_fragnumbers[...] = fragnumbers
            # get connectivity table
            cnc_table = self.ffe.mol.get_ctab()
            if len(cnc_table) > 0:
                # catch if there are no bonds at all: h5py does not like zero size selections
                cnc_table = np.array(cnc_table, dtype="i")
                mfp5_cnc_table = system.require_dataset("cnc_table", shape=cnc_table.shape, dtype=cnc_table.dtype)
                mfp5_cnc_table[...] = cnc_table
            system.attrs["bcd"] = self.ffe.mol.get_bcond()
            if "ff" in self.ffe.mol.loaded_addons:
                # a force field is loaded: we assume it is also initialized. TBI: a switch in the ff addon to verify this
                data = self.ffe.mol.ff.pack()
                ff = system.require_group("ff")
                ric = ff.require_group("ric")
                par = ff.require_group("par")
                if data["FF"] is not None:
                    par.attrs["FF"] = data["FF"]
                else:
                    par.attrs["FF"] = ""
                for r in ["bnd", "ang", "dih", "oop", "cha", "vdw"]:
                    # write ric arrays to mfp5 file
                    if r in data:
                        d = data[r]
                        fd = ric.require_dataset(r, shape=d.shape, dtype="i")
                        fd[...] = d
                    # write par data to mfp5 file
                    #  order is ptype, names, npars, pars
                    if r+"_par" in data:
                        p = data[r+"_par"]
                        fp = par.require_group(r)
                        n = len(p[0]) # number of paramters for this ric type
                        fptypes = fp.require_dataset("ptypes", shape=(n,), dtype=self.str_dt)
                        fptypes[...] = p[0]
                        fnames = fp.require_dataset("names", shape=(n,), dtype=self.str_dt)
                        fnames[...] = p[1]
                        fnpars = fp.require_dataset("npars", shape=p[2].shape, dtype="i")
                        fnpars[...] = p[2]
                        fpars = fp.require_dataset("pars", shape=p[3].shape, dtype="float64")
                        fpars[...] = p[3]
        return

    def compare_system(self):
        """ method to verify that the mol object in the ffe is consistent with the stuff on file
            NOTE: this can get arbitrary complex .. do a simple test of atypes here """
        self.open()
        OK = None
        if self.is_master:
            system = self.h5file["system"]
            mfp5_atypes = system["atypes"]
            atypes = self.ffe.mol.get_atypes()
            OK = True
            if len(mfp5_atypes)==len(atypes):
                for a1,a2 in zip(mfp5_atypes, atypes):
                    if a1 != a2:
                        OK = False
            else:
                OK = False
        OK = self.mpi_comm.bcast(OK)
        assert OK, "PDLP ERROR: The system in the mfp5 file is not equivalent to your actual system. Aborting!"
        return

    def get_mol_from_system(self, vel=False, restart_ff=True):
        """ read mol info from system group and generate a mol object

        in parallel this is done on the master only and the data is broadcasted to the other nodes

        Args:
            - vel (bool, otional): Defaults to False. If True: return velocity array in addtion to mol object
        """
        self.open()
        # read the data from file on the master
        if self.is_master:
            system = self.h5file["system"]
            # for python3: now all strings coming from hdf5 are byte class types .. need to convert to strings
            elems  = list(np.array(system["elems"]).astype("str"))
            atypes = list(np.array(system["atypes"]).astype("str"))
            fragtypes = list(np.array(system["fragtypes"]).astype("str"))
            fragnumbers = list(system["fragnumbers"])
            if "cnc_table" in list(system.keys()):
                cnc_table = list(system["cnc_table"])
            else:
                cnc_table = []
            bcd = system.attrs["bcd"]
        # broadcast if we run parallel
        if self.mpi_size > 1:
            if self.is_master:
                data = (elems, atypes, fragtypes, fragnumbers, cnc_table, bcd)
            else:
                data = None
            elems, atypes, fragtypes, fragnumbers, cnc_table, bcd = self.mpi_comm.bcast(data)
        # generate mol object
        na = len(elems)
        mol = molsys.mol()
        mol.set_natoms(na)
        mol.set_elems(elems)
        mol.set_atypes(atypes)
        mol.set_fragnumbers(fragnumbers)
        mol.set_fragtypes(fragtypes)
        mol.set_ctab(cnc_table, conn_flag=True)
        mol.bcond = bcd
        # now read the restart info that needs to be passed to the mol instance (Note: no velocities etc are read here)
        if self.is_master:
            try:
                rstage = self.h5file[self.restart_stage]
            except KeyError:
                self.pprint("PDLP ERROR: The requested restart stage %s does not exist in the file!" % self.restart_stage)
                raise
            restart = rstage["restart"]
            xyz = np.array(restart["xyz"], dtype="float64")
            cell = np.array(restart["cell"], dtype="float64")
            if vel:
                velocities = np.array(restart["vel"], dtype="float64")
        if self.mpi_size>1:
            if self.is_master:
                self.mpi_comm.Bcast(xyz)
                self.mpi_comm.Bcast(cell)
                if vel:
                    self.mpi_comm.Bcast(velocities)
            else:
                xyz = np.empty([na, 3], dtype="float64")
                cell = np.empty([3, 3], dtype="float64")
                self.mpi_comm.Bcast(xyz)
                self.mpi_comm.Bcast(cell)
                if vel:
                    velocities = np.empty([na, 3], dtype="float64")
                    self.mpi_comm.Bcast(velocities)
        mol.set_xyz(xyz)
        mol.set_cell(cell)
        # new check if addon data is present in the system group
        # start with ff
        ff_data = None
        if self.is_master:
            if ("ff" in list(system.keys())) and restart_ff:
                # ok, there is force field data in this file lets read it in as a packed directory
                ff_data = {}
                ff = system["ff"]
                ric = ff["ric"]
                par = ff["par"]
                ff_data["FF"] = par.attrs["FF"]
                for r in ["bnd", "ang", "dih", "oop", "cha", "vdw"]:
                    # read in ric integer arrays
                    if r in ric:
                        ff_data[r] = np.array(ric[r])
                    # now get the parameter
                    if r in par:
                        rpar = par[r]
                        ff_data[r+"_par"] = (\
                                        list(rpar["ptypes"]),\
                                        list(rpar["names"]),\
                                        np.array(rpar["npars"]),\
                                        np.array(rpar["pars"]))
        if self.mpi_size>1:
            ff_data = self.mpi_comm.bcast(ff_data)
        if ff_data is not None:
            # yes there was ff data .. set up mol addon ff
            mol.addon("ff")
            mol.ff.unpack(ff_data)
        if vel:
            return (mol, velocities)
        else:
            return mol

    def setup_molecules(self,mol):
        mol.addon('molecules')
        mm = mol.molecules
        molnames = ['xyz']
        for i,wm in enumerate(mm.whichmol):
            if wm != 0:
                if mol.fragtypes[i] not in molnames:
                    molnames.append(mol.fragtypes[i])
        mm.molnames = molnames
        mm.moltypes = [0] + [i+1 for i in mm.moltypes[1:]]
        return mol

    def add_stage(self, stage):
        """add a new stage

        Args:
            stage (string): stage name to add

        Returns:
            (bool) True if it worked and False if the stage already existed
        """
        if stage == "system":
            self.pprint("PDLP ERROR: Stage name system is not allowed!")
            raise IOError
        self.open()
        OK = True
        if self.is_master:
            if stage in list(self.h5file.keys()):
                OK = False
            else:
                sgroup = self.h5file.create_group(stage)
                rgroup = sgroup.create_group("restart")
                # generate restart arrays for xyz, cell and vel
                na = self.ffe.get_natoms()
                rgroup.require_dataset("xyz",shape=(na,3), dtype="float64")
                rgroup.require_dataset("vel",shape=(na,3), dtype="float64")
                rgroup.require_dataset("cell",shape=(3,3), dtype="float64")
        OK = self.mpi_comm.bcast(OK)
        return OK

    def has_stage(self,stagename):
        return stagename in self.get_stages()

    def get_stages(self):
        stagelist = list(self.h5file.keys())
        stagelist.remove("system")
        return stagelist

    def get_traj_from_stage(self, sname):
        if sname not in self.h5file.keys():
            return False
        stage = self.h5file[sname]
        if "traj" in stage.keys():
            traj = stage["traj"]
        else:
            traj = {}
        return traj

    def write_restart(self, stage=None, velocities=False):
        """write restart info to file

        Args:
            stage (string, optional): Defaults to None. Name of the stage to write restart (must exist). Uses current stage by default
            velocities (bool, optional): Defaults to False. When True writes also velocities
        """
        self.open()
        if stage is None:
            stage = self.stage
        OK = True
        # we need to do this in parallel on nodes (lammps gather is a parallel step)
        xyz = self.ffe.get_xyz()
        if velocities:
            vel = self.ffe.get_vel()
        if self.is_master:
            if stage not in list(self.h5file.keys()):
                OK = False
            else:
                restart = self.h5file[stage+"/restart"]
                cell = self.ffe.get_cell()
                rest_xyz = restart.require_dataset("xyz",shape=xyz.shape, dtype=xyz.dtype)
                rest_xyz[...] = xyz
                rest_cell = restart.require_dataset("cell", shape=cell.shape, dtype=cell.dtype)
                rest_cell[...] = cell
                if velocities:
                    rest_vel = restart.require_dataset("vel", shape=vel.shape, dtype=vel.dtype)
                    rest_vel[...] = vel
        OK = self.mpi_comm.bcast(OK)
        if not OK:
            self.pprint("PDLP ERROR: writing restart to stage %s failed. Stage does not exist" % stage)
            raise IOError
        return

    def prepare_stage(self, stage, traj_data, traj_nstep, data_nstep=1, thermo_values=[], prec="float32", tstep=0.001):
        """prepare a stage for trajectory writing

        Args:
            stage (string): name of the stage (must exist)
            traj_data (list of strings): name of the data to be written
            traj_nstep (int): frequency in steps to write
            data_nstep (int or list of ints), optional): Defaults to 1. freq to write each datatype
            thermo_values (list): Defaults to [], list of strings of thermodynamic data ... each code needs to know what is meant here
            prec (str, optional): Defaults to "float64". precision to use in writing traj data
            tstep (float, optional): Defaults to 0.001. Timestep in ps
        """
        # RS function rewritten for parallel runs (in lammps soem of the data_funcs need to be exectued by all nodes)
        self.open()
        if type(data_nstep)== type([]):
            assert len(data_nstep)==len(traj_data)
        else:
            data_nstep = len(traj_data)*[data_nstep]
        OK = True
        if self.is_master:
            if stage not in list(self.h5file.keys()):
                OK = False
        OK = self.mpi_comm.bcast(OK)
        if not OK:
            raise IOError("PDLP ERROR: preparing stge %s failed. Stage does not exist!" % stage)
        if self.is_master:
            traj = self.h5file.require_group(stage+"/traj")
            traj.attrs["nstep"] = traj_nstep
            traj.attrs["tstep"] = tstep
        for dname,dnstep in zip(traj_data, data_nstep):
            assert dname in list(self.ffe.data_funcs.keys())
            data = self.ffe.data_funcs[dname]()
            dshape = list(data.shape)
            if self.is_master:
                mfp5_data = traj.require_dataset(dname,
                                shape=tuple([1]+dshape),
                                maxshape=tuple([None]+dshape),
                                chunks=tuple([1]+dshape),
                                dtype=prec)
                mfp5_data[...] = data
                mfp5_data.attrs["nstep"] = dnstep
        if len(thermo_values)>0:
            # write thermodynamic data .. size of array is length of list
            nthermo = len(thermo_values)
            if self.is_master:
                mfp5_data = traj.require_dataset("thermo",
                                shape=(1, nthermo),
                                maxshape=(None, nthermo),
                                chunks=(1, nthermo),
                                dtype=prec)
                # TBI .. get thermo data in ffe on python level
                mfp5_data[...] = 0.0
                mfp5_data.attrs["labels"] = " ".join(thermo_values)
        self.mpi_comm.barrier()
        return

    def add_bondtab(self, stage, nbondsmax, bond_tab = 0, bond_order = 0.0):
        """generates a bondtab/bondorder entry in the trajectory group of the current stage for ReaxFF

        Args:
            stage (string): name of the current stage
            nbondsmax (int): size of the tables
            bond_tab (list) : bond table
            bond_order (list) : bond order
        """
        if self.is_master:
            st = self.h5file[stage]
            traj = st["traj"]
            bondtab = traj.require_dataset("bondtab",
                            shape=(1,nbondsmax,2),
                            maxshape=(None, nbondsmax,2),
                            chunks=(1,nbondsmax,2),
                            dtype = "int32")
            bondord = traj.require_dataset("bondord",
                            shape=(1,nbondsmax),
                            maxshape=(None, nbondsmax),
                            chunks=(1,nbondsmax),
                            dtype = "float32")
            bondtab[...] = bond_tab
            bondord[...] = bond_order
        return

    def __call__(self, force_wrest=False):
        """ this generic routine is called to save restart info to the hdf5 file
            if force_wrest is True then restart is written in any case """
        if not self.pd: raise IOError("No pydlpoly instance")
        data_written = False
        if (self.counter%self.rest_nstep == 0) or force_wrest:
            # restart is written .. fetch all data
            for d in self.rest_data:
                if self.verbose:
                    self.pprint("Writing restart data %s to mfp5 file" % d)
                data = self.data_funcs[d]()
                self.rest_datasets[d][...] = data
            data_written = True
        if self.track_data != None:
            for i,d in enumerate(self.track_data):
                if (self.counter%self.traj_nstep[i] == 0):
                    self.traj_frame[i] += 1
                    data = self.data_funcs[d]()
                    tds = self.traj_datasets[d]
                    tds.resize(self.traj_frame[i], axis=0)
                    tds[self.traj_frame[i]-1,...] = data
                    data_written = True
        # now we are done
        if data_written : self.h5file.flush()
        self.counter += 1
        return


text:
#!/usr/bin/env python
import rospy
from sensor_msgs.msg import Image
from cv_bridge import CvBridge, CvBridgeError
import cv2
import numpy as np
import matplotlib.pyplot as plt

class DetectPPL(object):

    def __init__(self):

        self.image_sub = rospy.Subscriber("/camera/rgb/image_raw",Image,self.camera_callback)
        self.bridge_object = CvBridge()

    def camera_callback(self,data):
        try:
            # We select bgr8 because its the OpenCV encoding by default
            cv_image = self.bridge_object.imgmsg_to_cv2(data, desired_encoding="bgr8")
        except CvBridgeError as e:
            print(e)

        photo = cv_image

        #Reduce photo size to reduce sensitivity
        photo = cv2.resize(photo,(600, 400))
        gray_photo = cv2.cvtColor(photo, cv2.COLOR_RGB2GRAY)
        # ret,thresh = cv2.adaptiveThreshold(gray_photo,255,cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY,11,2)

        blur = cv2.GaussianBlur(gray_photo,(5,5),0)
        ret,thresh = cv2.threshold(blur,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
        cv2.imshow("Thresh", thresh)
        cv2.imshow('Detection',photo)

        cv2.waitKey(1)

def main():
    detect_wanted_object = DetectPPL()
    rospy.init_node('detect_wanted_node', anonymous=True)

    try:
        rospy.spin()

    except KeyboardInterrupt:
        print("Shutting down")
    cv2.destroyAllWindows()

if __name__ == '__main__':
    main()


text:
# -*- coding: utf-8 -*-
"""
Created on Fri Jun 30 17:56:30 2017

@author: Paul
"""
import os,sys
sys.path.append('C:\Users\Paul\Documents\These\codes')

from scipy import ndimage as nd
from skimage import io,exposure,measure
import PIL
import numpy as np

from Savitzjy_goulay import savitzky_golay
from remove_small_objects import remove_small_objects
from my_tools_img_proc import finddata

class myImage(object):
    def __init__(self,filename,dimension):
        self.filename=filename
        self.ini = io.imread(filename)
        self.img = self.ini[:dimension]
        self.isPrepared = False
        self.segparam = 0
        self.pR = 1
        self.size = np.size(self.img)
        self.char = {}

    def set_segparam(self,segparam):
        '''
        Set the segmentation paramter to an arbitrary value. If given, the code
        will not try to find it and take this one instead.
        '''
        self.segparam = segparam

    def preparerMedian(self,param = 6,
                       size = (2,2),
                       method = 'find',
                       obj_size = 80,
                       gamma=False):

        self.median = nd.filters.median_filter(self.img, size = size)

        if self.segparam == 0 :
            self.segparam = find_segparam(self.median,gamma = gamma)

        tmp = remove_small_objects(self.median>self.segparam,obj_size)
        self.ready = np.logical_not(remove_small_objects(np.logical_not(tmp),obj_size))
        self.isPrepared = True

    def getChar(self):
        if not self.isPrepared:
            self.preparerMedian()
            print(u'Image en préparation automatique, attention aux paramètres')
        self.char['taux'] = 100*np.sum(self.ready)/float(self.size)
        self.char['temp'],self.char['time'] = finddata(self.filename)
        self.char['perimeter']=measure.perimeter(self.ready)

max_stars_repo_name: Samarargon/sem-proc

max_issues_repo_name: Samarargon/sem-proc

max_forks_repo_name: Samarargon/sem-proc

reason: import numpy,from scipy


text:
"""Implementations of algorithms for continuous control."""

import pdb
from typing import Optional, Sequence, Tuple

import flax
import jax
import jax.numpy as jnp
import numpy as np

from jax_rl.agents.actor_critic_temp import ActorCriticTemp
from jax_rl.agents.sac import actor, critic, temperature
from jax_rl.datasets import Batch
from jax_rl.networks import critic_net, policies
from jax_rl.networks.common import InfoDict, create_model

@jax.partial(jax.jit, static_argnums=(2, 3, 4, 5))
def _update_jit(sac: ActorCriticTemp, batch: Batch, discount: float,
                tau: float, target_update_period: int,
                target_entropy: float) -> Tuple[ActorCriticTemp, InfoDict]:

    sac, critic_info = critic.update(sac, batch, discount, soft_critic=True)
    sac = critic.target_update(sac, tau, target_update_period)

    sac, actor_info = actor.update(sac, batch)
    sac, alpha_info = temperature.update(sac, actor_info['entropy'],
                                         target_entropy)

    return sac, {**critic_info, **actor_info, **alpha_info}

class SACLearner(object):
    def __init__(self,
                 seed: int,
                 observations: jnp.ndarray,
                 actions: jnp.ndarray,
                 actor_lr: float = 3e-4,
                 critic_lr: float = 3e-4,
                 temp_lr: float = 3e-4,
                 hidden_dims: Sequence[int] = (256, 256),
                 discount: float = 0.99,
                 tau: float = 0.005,
                 target_update_period: int = 1,
                 target_entropy: Optional[float] = None,
                 init_temperature: float = 1.0):

        action_dim = actions.shape[-1]

        if target_entropy is None:
            self.target_entropy = -action_dim / 2
        else:
            self.target_entropy = target_entropy

        self.tau = tau
        self.target_update_period = target_update_period
        self.discount = discount

        # pdb.set_trace()
        rng = jax.random.PRNGKey(seed)
        rng, actor_key, critic_key, temp_key = jax.random.split(rng, 4)

        # One actor model and two critic model?
        # [creat_model] means model + params?
        actor = create_model(
            policies.NormalTanhPolicy(hidden_dims, action_dim),
            [actor_key, observations])
        actor = actor.with_optimizer(flax.optim.Adam(learning_rate=actor_lr))

        critic = create_model(critic_net.DoubleCritic(hidden_dims),
                              [critic_key, observations, actions])
        critic = critic.with_optimizer(
            flax.optim.Adam(learning_rate=critic_lr))
        target_critic = create_model(critic_net.DoubleCritic(hidden_dims),
                                     [critic_key, observations, actions])

        temp = create_model(temperature.Temperature(init_temperature),
                            [temp_key])
        temp = temp.with_optimizer(flax.optim.Adam(learning_rate=temp_lr))

        # sac have 5 params: actor, 2 * critic, temp, rng?
        self.sac = ActorCriticTemp(actor=actor,
                                   critic=critic,
                                   target_critic=target_critic,
                                   temp=temp,
                                   rng=rng)

    def sample_actions(self,
                       observations: np.ndarray,
                       temperature: float = 1.0) -> jnp.ndarray:
        rng, actions = policies.sample_actions(self.sac.rng, self.sac.actor.fn,
                                               self.sac.actor.optimizer.target,
                                               observations, temperature)

        self.sac = self.sac.replace(rng=rng)

        actions = np.asarray(actions)
        return np.clip(actions, -1, 1)

    def update(self, batch: Batch) -> InfoDict:
        self.sac, info = _update_jit(self.sac, batch, self.discount, self.tau,
                                     self.target_update_period,
                                     self.target_entropy)
        return info

reason: import numpy,import jax,from jax


text:
"""Initializer for bootstraping initial states."""
from typing import Optional

from absl import logging
import jax
from jax import random
import jax.numpy as jnp
import jax.scipy.special as jsp

from sympais import constraint
from sympais import profiles
from sympais import realpaver_solver
from sympais import types
from sympais import z3_solver

class Initializer:
  """Initializer for finding feasible initial solutions for bootstraping MCMC chains.

  This implements the solver used in SYMPAIS to find initial feasible solution.
  """

  def __init__(
      self,
      profile: profiles.Profile,
      constraints: types.Constraints,
      domains: types.DomainMap,
      init_strategy: str,
      resample: bool,
  ):
    """Create a initializer for finding initial feasible solutions.
      Args:
        profile: input `Profile` for the probabilistic program.
        constraints: input `Constraint` for the program analysis query.
        domains: input `Domain` for the input profile.
        init_strategy: "z3" or "realpaver", strategy used for finding initial solution.
        resample: if True, then the `num_solutions` of initial solution are found by
          re-sampling the set of candidate solutions via importance sampling (IS).
          Otherwise, the candidate solutions will be chosen uniformly.
      Raises:
        ValueError: if init_strategy is not in ("z3", "realpaver").
    """
    if init_strategy not in ("z3", "realpaver"):
      raise ValueError(f"Unknown initialization method {init_strategy}")
    self._profile = profile
    self._constraints = constraints
    self._domains = domains
    self._init = init_strategy
    self._resample = resample
    constraint_fn = jax.jit(constraint.build_constraint_fn(constraints, domains))
    self._target_log_prob_fn = jax.jit(
        constraint.build_target_log_prob_fn(profile, domains, constraint_fn))

  @staticmethod
  def _find_feasible_solution_realpaver(constraints: types.Constraints,
                                        domains: Optional[types.DomainMap] = None):
    constraint_fn = constraint.build_constraint_fn(constraints, domains)
    # First find interval boxes using RealPaver
    boxes = realpaver_solver.RealPaverSolver().solve(
        constraints, domains=domains, precision=0.01)
    solutions = []
    # All centers of the inner box would satisfy the constraints
    for box in list(boxes.inner_boxes):
      values = {
          site: interval.lower + (interval.upper - interval.lower) / 2
          for site, interval in box.intervals.items()
      }
      solutions.append(values)
    # Here we use the heuristic to find feasible solution
    # by computing the geometric center of each interval box and test for
    # constraint satisfaction.
    for box in list(boxes.outer_boxes):
      values = {
          site: interval.lower + (interval.upper - interval.lower) / 2
          for site, interval in box.intervals.items()
      }
      if constraint_fn(values):
        solutions.append(values)
    return solutions

  def __call__(self, num_solutions: int, rng: types.PRNGKey):
    """Find initial feasible solutions for the MCMC chains.
    Args:
      num_solutions: the number of solutions to be returned.
      rng: random seed to use
    """
    if self._init == "z3":
      states = [z3_solver.Z3Solver().solve(self._constraints, self._domains)]
      proposal_indices = random.randint(rng, (num_solutions,), minval=0, maxval=1)
      states_selected = []
      for pi in proposal_indices:
        states_selected.append(states[pi])
      proposal_states = (
          jax.tree_multimap(  # pylint: disable=no-value-for-parameter
              lambda *xs: jnp.stack(xs), *states_selected))
    elif self._init == "realpaver":
      feasible_states = self._find_feasible_solution_realpaver(
          self._constraints, self._domains)
      # Also add z3 feasible states. This guards against failures of not
      # getting a single feasible solution from RP.
      z3_feasible_states = [
          z3_solver.Z3Solver().solve(
              self._constraints,
              self._domains,
          )
      ]
      feasible_states.extend(z3_feasible_states)
      logging.info("Number of feasible solutions found is %d.", len(feasible_states))
      if len(feasible_states) < num_solutions:
        # First ensure that all feasible states are selected at least once.
        states_selected = []
        for state in feasible_states:
          states_selected.append(state)
        # Populate the rest of the states by sampling from feasible solutions
        if self._resample:
          logp = jnp.array(
              [self._target_log_prob_fn(state) for state in feasible_states])
          logp = logp - jsp.logsumexp(logp)
        else:
          logp = jnp.ones((len(feasible_states),)) / len(feasible_states)
        logging.debug("Initial state probs %s", jnp.exp(logp))
        proposal_indices = random.choice(
            rng,
            len(feasible_states),
            shape=(num_solutions - len(feasible_states),),
            replace=True,
            p=jnp.exp(logp),
        )
        for pi in proposal_indices:
          states_selected.append(feasible_states[pi])
      else:
        if self._resample:
          logp = jnp.array(
              [self._target_log_prob_fn(state) for state in feasible_states])
          logp = logp - jsp.logsumexp(logp)
        else:
          logp = jnp.ones((len(feasible_states),)) / len(feasible_states)
        logging.debug("Initial state probs %s", jnp.exp(logp))
        proposal_indices = random.choice(
            rng,
            len(feasible_states),
            shape=(num_solutions,),
            replace=True,
            p=jnp.exp(logp),
        )
        states_selected = []
        for pi in proposal_indices:
          states_selected.append(feasible_states[pi])
      proposal_states = (
          jax.tree_multimap(  # pylint: disable=no-value-for-parameter
              lambda *xs: jnp.stack(xs), *states_selected))
    else:
      raise ValueError("Unknown initialization strategy: {}".format(self._init))
    return proposal_states

max_stars_repo_name: ethanluoyc/sympais

max_issues_repo_name: ethanluoyc/sympais

max_forks_repo_name: ethanluoyc/sympais

reason: import jax,from jax


text:
import unittest
import tensorflow as tf
import numpy as np
import spn.experiments.RandomSPNs.RAT_SPN as RAT_SPN
import spn.experiments.RandomSPNs.region_graph as region_graph

import spn.algorithms.Inference as Inference

class TestRatSpn(unittest.TestCase):
    def test_inference_results(self):
        np.random.seed(123)
        tf.set_random_seed(123)

        num_dims = 20

        rg = region_graph.RegionGraph(range(num_dims))
        for _ in range(0, 10):
            rg.random_split(2, 3)

        args = RAT_SPN.SpnArgs()
        args.normalized_sums = True
        spn = RAT_SPN.RatSpn(10, region_graph=rg, name="obj-spn", args=args)

        sess = tf.Session()
        sess.run(tf.global_variables_initializer())
        dummy_input = np.random.normal(0.0, 1.2, [10, num_dims])
        input_ph = tf.placeholder(tf.float32, [10, num_dims])
        output_tensor = spn.forward(input_ph)
        tf_output = sess.run(output_tensor, feed_dict={input_ph: dummy_input})

        output_nodes = spn.get_simple_spn(sess)
        simple_output = []
        for node in output_nodes:
            simple_output.append(Inference.likelihood(node, dummy_input))
        simple_output = np.stack(simple_output)
        deviation = simple_output / np.exp(tf_output)
        rel_error = np.abs(deviation - 1.0)
        # print(rel_error)

        self.assertTrue(np.all(rel_error < 1e-2))

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

max_issues_repo_name: minimrbanana/SPFlow

max_forks_repo_name: minimrbanana/SPFlow


text:
from enum import Enum
from PuzzleLib import Config

ConvFwdAlgo = None
ConvBwdFilterAlgo = None
ConvBwdDataAlgo = None
convNd = None
convNdBackwardData = None
convNdBackwardParams = None

convNdbenchmark = None

deconvNd = None
deconvNdBackwardData = None
deconvNdBackwardParams = None

PoolMode = None
poolNd = None
poolNdBackward = None

BatchNormMode = None
batchNormNd = None
batchNormNdBackward = None

SoftMaxMode = None
softmaxNd = None
softmaxNdBackward = None

mapLRN = None
mapLRNBackward = None

crossMapLRN = None
crossMapLRNBackward = None

instanceNorm2d = None
instanceNorm2dBackward = None

spatialTf = None
spatialTfBackward = None

RNNMode = None
DirectionMode = None

createRnn = None

acquireRnnParams = None
updateRnnParams = None

forwardRnn = None
backwardDataRnn = None
backwardParamsRnn = None

deviceSupportsBatchHint = None

def autoinit():
	if not Config.shouldInit():
		return

	if Config.backend == Config.Backend.cuda:
		initCuda()
	elif Config.backend == Config.Backend.hip:
		initHip()
	elif Config.backend == Config.Backend.cpu:
		initCPU()
	elif Config.backend == Config.Backend.intel:
		initIntel()
	else:
		raise Config.ConfigError(Config.backend)

def initCuda():
	from PuzzleLib.Cuda import Backend

	backend = initGPU(Backend)
	memoryPool, dnn = backend.memoryPool, backend.dnn

	def wrapPoolNd(data, size, stride, pad, mode, test):
		return dnn.poolNd(data, size, stride, pad, mode.value, None, memoryPool), None

	def wrapPoolNdBackward(indata, outdata, grad, _, size, stride, pad, mode):
		return dnn.poolNdBackward(grad, indata, outdata, size, stride, pad, mode.value, None, memoryPool)

	global PoolMode, poolNd, poolNdBackward
	PoolMode = backend.PoolMode
	poolNd = wrapPoolNd
	poolNdBackward = wrapPoolNdBackward

	def wrapMapLRN(data, means, N, alpha, beta, K, test):
		return dnn.mapLRN(data, means, N, alpha, beta, K, allocator=memoryPool), None

	def wrapMapLRNBackward(data, _, grad, means, __, N, alpha, beta, K):
		return dnn.mapLRNBackward(data, grad, means, N, alpha, beta, K, allocator=memoryPool)

	global mapLRN, mapLRNBackward
	mapLRN = wrapMapLRN
	mapLRNBackward = wrapMapLRNBackward

	def wrapCrossMapLRN(data, N, alpha, beta, K, test):
		return dnn.crossMapLRN(data, N, alpha, beta, K, allocator=memoryPool), None

	def wrapCrossMapLRNBackward(data, outdata, grad, _, N, alpha, beta, K):
		return dnn.crossMapLRNBackward(data, outdata, grad, N, alpha, beta, K, allocator=memoryPool)

	global crossMapLRN, crossMapLRNBackward
	crossMapLRN = wrapCrossMapLRN
	crossMapLRNBackward = wrapCrossMapLRNBackward

	def wrapSpatialTf(data, transform, outshape, getGrid):
		return dnn.spatialTf(data, transform, outshape, getGrid, allocator=memoryPool)

	def wrapSpatialTfBackward(grad, data, grid):
		return dnn.spatialTfBackward(grad, data, grid, allocator=memoryPool)

	global spatialTf, spatialTfBackward
	spatialTf = wrapSpatialTf
	spatialTfBackward = wrapSpatialTfBackward

def initHip():
	from PuzzleLib.Hip import Backend

	backend = initGPU(Backend)
	memoryPool, dnn = backend.memoryPool, backend.dnn

	def wrapPoolNd(data, size, stride, pad, mode, test):
		result = dnn.poolNd(data, size, stride, pad, mode.value, test, None, memoryPool)
		return result if not test else (result, None)

	def wrapPoolNdBackward(indata, outdata, grad, workspace, size, stride, pad, mode):
		return dnn.poolNdBackward(grad, indata, outdata, workspace, size, stride, pad, mode.value, None, memoryPool)

	global PoolMode, poolNd, poolNdBackward
	PoolMode = backend.PoolMode
	poolNd = wrapPoolNd
	poolNdBackward = wrapPoolNdBackward

	def wrapMapLRN(data, means, N, alpha, beta, K, test):
		assert means is None

		result = dnn.lrn(data, N, alpha, beta, K, backend.LRNMode.map.value, test, allocator=memoryPool)
		return result if not test else (result, None)

	def wrapMapLRNBackward(data, outdata, grad, means, workspace, N, alpha, beta, K):
		assert means is None
		return dnn.lrnBackward(
			grad, data, outdata, workspace, N, alpha, beta, K, backend.LRNMode.map.value, allocator=memoryPool
		)

	global mapLRN, mapLRNBackward
	mapLRN = wrapMapLRN
	mapLRNBackward = wrapMapLRNBackward

def initGPU(Backend):
	backend = Backend.getBackend(Config.deviceIdx, initmode=1, logger=Config.getLogger())

	initBaseGPU(backend)
	initInstanceNormGPU(backend)
	initRnnGPU(backend)

	return backend

def initBaseGPU(backend):
	import numpy as np
	memoryPool, dnn = backend.memoryPool, backend.dnn

	global ConvFwdAlgo, ConvBwdDataAlgo, ConvBwdFilterAlgo
	ConvFwdAlgo = backend.ConvFwdAlgo
	ConvBwdDataAlgo = backend.ConvBwdDataAlgo
	ConvBwdFilterAlgo = backend.ConvBwdFilterAlgo

	def wrapConvNd(data, W, bias, stride, pad, dilation, groups, algo):
		return dnn.convNd(
			data, W, bias.ravel() if bias is not None else None, stride, pad, dilation, groups,
			algo.value, None, memoryPool
		)

	def wrapConvNdBackwardData(grad, W, data, stride, pad, dilation, groups, algo):
		return dnn.convNdBackwardData(
			grad, W, None, data, stride, pad, dilation, None, groups, algo.value, None, memoryPool
		)

	def wrapConvNdBackwardParams(data, grad, W, bias, stride, pad, dilation, groups,
								 wgrad, bgrad, scale, momentum, algo):
		return dnn.convNdBackwardParams(
			data, grad, W, stride, pad, dilation, groups, bias is not None, False,
			wgrad, bgrad.ravel() if bgrad is not None else None, scale, momentum, algo.value, memoryPool
		)

	global convNd, convNdBackwardData, convNdBackwardParams
	convNd = wrapConvNd
	convNdBackwardData = wrapConvNdBackwardData
	convNdBackwardParams = wrapConvNdBackwardParams

	def wrapConvNdbenchmark(datashape, Wshape, stride, pad, dilation, groups, transpose):
		fwdResults, bwdDataResults, bwdParamResults = backend.convNdbenchmark(
			datashape, Wshape, np.float32, stride, pad, dilation, groups
		)

		return fwdResults, bwdParamResults, bwdDataResults

	global convNdbenchmark
	convNdbenchmark = wrapConvNdbenchmark

	def wrapDeconvNd(data, W, bias, stride, pad, dilation, postpad, groups, algo):
		return dnn.convNdBackwardData(
			data, W, bias.ravel() if bias is not None else None, None, stride, pad, dilation, postpad, groups,
			algo.value, None, memoryPool
		)

	def wrapDeconvNdBackwardData(grad, W, data, stride, pad, dilation, groups, algo):
		assert data is not None
		return dnn.convNd(grad, W, None, stride, pad, dilation, groups, algo.value, None, memoryPool)

	def wrapDeconvNdBackwardParams(data, grad, W, bias, stride, pad, dilation, groups,
								   wgrad, bgrad, scale, momentum, algo):
		return dnn.convNdBackwardParams(
			grad, data, W, stride, pad, dilation, groups, bias is not None, True,
			wgrad, bgrad.ravel() if bgrad is not None else None, scale, momentum, algo.value, memoryPool
		)

	global deconvNd, deconvNdBackwardData, deconvNdBackwardParams
	deconvNd = wrapDeconvNd
	deconvNdBackwardData = wrapDeconvNdBackwardData
	deconvNdBackwardParams = wrapDeconvNdBackwardParams

	global BatchNormMode
	BatchNormMode = backend.BatchNormMode

	def wrapBatchNormNd(data, scale, bias, mean, var, epsilon, factor, test, mode=BatchNormMode.spatial, out=None):
		shape = scale.shape
		result = dnn.batchNormNd(
			data, mean.ravel(), var.ravel(), scale.ravel(), bias.ravel(), epsilon, factor, test, mode.value, out=out,
			allocator=memoryPool
		)

		if test:
			return result

		outdata, savemean, saveinvvar = result
		return outdata, savemean.reshape(shape), saveinvvar.reshape(shape)

	def wrapBatchNormNdBackward(data, grad, scale, savemean, saveinvvar, epsilon, mode=BatchNormMode.spatial):
		shape = scale.shape
		ingrad, scalegrad, bgrad = dnn.batchNormNdBackward(
			grad, data, scale.ravel(), savemean.ravel(), saveinvvar.ravel(), epsilon, mode.value, allocator=memoryPool
		)

		return ingrad, scalegrad.reshape(shape), bgrad.reshape(shape)

	global batchNormNd, batchNormNdBackward
	batchNormNd = wrapBatchNormNd
	batchNormNdBackward = wrapBatchNormNdBackward

	global SoftMaxMode
	SoftMaxMode = backend.SoftMaxMode

	def wrapSoftmaxNd(data, mode=SoftMaxMode.spatial):
		return dnn.softmaxNd(data, mode.value, allocator=memoryPool)

	def wrapSoftmaxNdBackward(outdata, grad):
		return dnn.softmaxNdBackward(grad, outdata, allocator=memoryPool)

	global softmaxNd, softmaxNdBackward
	softmaxNd = wrapSoftmaxNd
	softmaxNdBackward = wrapSoftmaxNdBackward

def initInstanceNormGPU(backend):
	memoryPool = backend.memoryPool

	def wrapInstanceNorm2d(data, scale, bias, epsilon=1e-5):
		return backend.instanceNorm2d(data, scale.ravel(), bias.ravel(), epsilon, allocator=memoryPool)

	def wrapInstanceNorm2dBackward(grad, data, extscale, savemean, saveinvvar, epsilon, affine):
		return backend.instanceNorm2dBackward(
			grad, data, extscale, savemean, saveinvvar, epsilon, affine, allocator=memoryPool
		)

	global instanceNorm2d, instanceNorm2dBackward
	instanceNorm2d = wrapInstanceNorm2d
	instanceNorm2dBackward = wrapInstanceNorm2dBackward

def initRnnGPU(backend):
	import numpy as np
	memoryPool = backend.memoryPool

	global RNNMode, DirectionMode
	RNNMode = backend.RNNMode
	DirectionMode = backend.DirectionMode

	def wrapCreateRnn(insize, hsize, layers, mode, direction, dropout, seed, batchsize):
		rnn, W, params = backend.createRnn(
			insize, hsize, np.float32, layers, mode=mode, direction=direction, dropout=dropout,
			seed=seed, batchsize=0 if batchsize is None else batchsize
		)

		return rnn, W, {i: layer for i, layer in enumerate(params)}

	def wrapAcquireRnnParams(descRnn, w):
		params = backend.acquireRnnParams(descRnn, w)
		return w, params

	def wrapUpdateRnnParams(descRnn, w, params):
		params = [params[layer] for layer in sorted(params.keys())]
		backend.updateRnnParams(descRnn, w, params)

	global createRnn, acquireRnnParams, updateRnnParams
	createRnn = wrapCreateRnn
	acquireRnnParams = wrapAcquireRnnParams
	updateRnnParams = wrapUpdateRnnParams

	def wrapForwardRnn(data, W, descRnn, test=False):
		return descRnn.forward(data, W, test=test, allocator=memoryPool)

	def wrapBackwardDataRnn(grad, outdata, W, reserve, descRnn):
		ingrad, _, _ = descRnn.backwardData(grad, outdata, W, reserve, allocator=memoryPool)
		return ingrad, reserve

	def wrapBackwardParamsRnn(data, outdata, _, reserve, descRnn):
		return descRnn.backwardParams(data, outdata, reserve, allocator=memoryPool)

	global forwardRnn, backwardDataRnn, backwardParamsRnn
	forwardRnn = wrapForwardRnn
	backwardDataRnn = wrapBackwardDataRnn
	backwardParamsRnn = wrapBackwardParamsRnn

	global deviceSupportsBatchHint
	deviceSupportsBatchHint = backend.deviceSupportsBatchHint

	return backend

def initCPU():
	from PuzzleLib.CPU.Wrappers import NumpyDnn

	def wrapConvNd(data, W, bias, stride, pad, dilation, groups, algo):
		assert groups == 1
		return NumpyDnn.conv2d(data, W, bias, stride, pad, dilation)

	global convNd, convNdBackwardData, convNdBackwardParams
	convNd = wrapConvNd

	def wrapPoolNd(data, size, stride, pad, mode, test):
		return NumpyDnn.pool2d(data, size, stride, pad, mode), None

	global PoolMode, poolNd, poolNdBackward
	PoolMode = NumpyDnn.PoolMode
	poolNd = wrapPoolNd

	class ProxyBatchNormMode(Enum):
		perActivation = 0
		spatial = 1

	def wrapBatchNormNd(data, scale, bias, mean, var, epsilon, factor, test, mode=None, out=None):
		outdata = NumpyDnn.batchNorm2d(data, scale, bias, mean, var, epsilon, test, out=out)
		return outdata if test else (outdata, mean, var)

	global BatchNormMode, batchNormNd
	BatchNormMode = ProxyBatchNormMode
	batchNormNd = wrapBatchNormNd

	global deviceSupportsBatchHint
	deviceSupportsBatchHint = lambda: False

def initIntel():
	from PuzzleLib.Intel.Wrappers import DNNL, DNNLInstanceNorm

	global ConvFwdAlgo, ConvBwdDataAlgo, ConvBwdFilterAlgo
	ConvFwdAlgo = DNNL.ConvAlgo
	ConvBwdDataAlgo = DNNL.ConvAlgo
	ConvBwdFilterAlgo = DNNL.ConvAlgo

	def wrapConvNd(data, W, bias, stride, pad, dilation, groups, algo):
		assert groups == 1
		return DNNL.convNd(data, W, bias, stride, pad, dilation, algo=algo)

	def wrapConvNdBackwardData(grad, W, data, stride, pad, dilation, groups, algo):
		assert groups == 1
		return DNNL.convNdBackwardData(grad, W, data, stride, pad, dilation, algo=algo)

	def wrapConvNdBackwardParams(data, grad, W, bias, stride, pad, dilation, groups, wgrad, bgrad, scale, momentum,
								 algo):
		assert groups == 1
		return DNNL.convNdBackwardParams(
			data, grad, W, bias, stride, pad, dilation, wgrad, bgrad, scale, momentum, algo=algo
		)

	global convNd, convNdBackwardData, convNdBackwardParams
	convNd = wrapConvNd
	convNdBackwardData = wrapConvNdBackwardData
	convNdBackwardParams = wrapConvNdBackwardParams

	def wrapConvNdbenchmark(datashape, Wshape, stride, pad, dilation, groups, transpose):
		assert groups == 1
		return DNNL.convNdbenchmark(datashape, Wshape, stride, pad, dilation, transpose)

	global convNdbenchmark
	convNdbenchmark = wrapConvNdbenchmark

	def wrapDeconvNd(data, W, bias, stride, pad, dilation, groups, postpad, algo):
		assert groups == 1
		return DNNL.convNd(data, W, bias, stride, pad, dilation, postpad, algo=algo, transpose=True)

	def wrapDeconvNdBackwardData(grad, W, data, stride, pad, dilation, groups, algo):
		assert groups == 1
		return DNNL.convNdBackwardData(grad, W, data, stride, pad, dilation, algo=algo, transpose=True)

	def wrapDeconvNdBackwardParams(data, grad, W, bias, stride, pad, dilation, groups, wgrad, bgrad, scale, momentum,
								   algo):
		assert groups == 1
		return DNNL.convNdBackwardParams(
			data, grad, W, bias, stride, pad, dilation, wgrad, bgrad, scale, momentum, algo=algo, transpose=True
		)

	global deconvNd, deconvNdBackwardData, deconvNdBackwardParams
	deconvNd = wrapDeconvNd
	deconvNdBackwardData = wrapDeconvNdBackwardData
	deconvNdBackwardParams = wrapDeconvNdBackwardParams

	def wrapPoolNd(data, size, stride, pad, mode, test):
		result = DNNL.poolNd(data, size, stride, pad, mode, test)
		return (result, None) if test else (result[0], result[1:])

	def wrapPoolNdBackward(indata, outdata, grad, workspace, size, stride, pad, mode):
		workspace, desc = workspace
		return DNNL.poolNdBackward(indata, grad, workspace, desc, size, stride, pad, mode)

	global PoolMode, poolNd, poolNdBackward
	PoolMode = DNNL.PoolMode
	poolNd = wrapPoolNd
	poolNdBackward = wrapPoolNdBackward

	class ProxyBatchNormMode(Enum):
		perActivation = 0
		spatial = 1

	def wrapBatchNormNd(data, scale, bias, mean, var, epsilon, factor, test, mode=None, out=None):
		outdata, mean, var, desc = DNNL.batchNormNd(data, scale, bias, mean, var, epsilon, test, out=out)
		return outdata if test else (outdata, mean, var)

	global BatchNormMode, batchNormNd
	BatchNormMode = ProxyBatchNormMode
	batchNormNd = wrapBatchNormNd

	global softmaxNd, softmaxNdBackward
	softmaxNd = DNNL.softmaxNd
	softmaxNdBackward = DNNL.softmaxNdBackward

	def wrapCrossMapLRN(data, N, alpha, beta, K, test):
		result = DNNL.lrn(data, DNNL.LRNMode.cross, N, alpha, beta, K, test)
		return (result[0], result[1:]) if not test else (result, None)

	def wrapCrossMapLRNBackward(data, outdata, grad, workspace, N, alpha, beta, K):
		workspace, desc = workspace
		return DNNL.lrnBackward(data, grad, workspace, desc, DNNL.LRNMode.cross, N, alpha, beta, K)

	global crossMapLRN, crossMapLRNBackward
	crossMapLRN = wrapCrossMapLRN
	crossMapLRNBackward = wrapCrossMapLRNBackward

	def wrapInstanceNorm2d(data, scale, bias, epsilon):
		result = DNNLInstanceNorm.instanceNorm2d(data, scale, bias, epsilon)

		outdata, savemean, savevar, extscale, extbias, desc = result
		return outdata, savemean, savevar, (extscale, extbias, desc)

	def wrapInstanceNorm2dBackward(grad, data, exts, savemean, savevar, epsilon, affine):
		extscale, extbias, desc = exts
		return DNNLInstanceNorm.instanceNorm2dBackward(
			grad, data, extscale, extbias, savemean, savevar, epsilon, desc, affine
		)

	global instanceNorm2d, instanceNorm2dBackward
	instanceNorm2d = wrapInstanceNorm2d
	instanceNorm2dBackward = wrapInstanceNorm2dBackward

	global deviceSupportsBatchHint
	deviceSupportsBatchHint = lambda: False

autoinit()

max_stars_repo_name: EmilPi/PuzzleLib

max_issues_repo_name: EmilPi/PuzzleLib

max_forks_repo_name: EmilPi/PuzzleLib


text:
import autograd.numpy as np
import scipy.sparse as sparse

# These work for sparse as well as CondensedA

def init(nx, nu, N, stateLim=True, inputLim=True, periodic=False, polyBlocks=None):
    '''Return scipy sparse.
    If polyBlocks is specified, additional polyhedron membership constraints are added at the bottom.
    The size of these needs to be fixed to maintain the sparsity structure (i.e. Nc-sized polyhedron).
    polyBlocks = [polyBlock1, ...], where
    polyBlock = [xstart, Nc, Ncx], where
    Apoly @ x[xstart : xstart+Ncx] <= upoly,
    with both sides having Nc rows (# constraints)
    '''
    Ad = np.ones((nx, nx))
    Bd = np.ones((nx, nu))
    # Ad, Bd = getLin(x0, ctrl, dt)
    Ax = sparse.kron(sparse.eye(N+1),-sparse.eye(nx)) + sparse.kron(sparse.eye(N+1, k=-1), Ad)
    Bu = sparse.kron(sparse.vstack([sparse.csc_matrix((1, N)), sparse.eye(N)]), Bd)
    # LTV dynamics constraint
    Aeq = sparse.hstack([Ax, Bu])

    if periodic:
        # Equality constraint x0 = xN
        Aeq = sparse.vstack((Aeq,
            sparse.hstack((-sparse.eye(nx), np.zeros((nx, (N-1)*nx)), sparse.eye(nx), np.zeros((nx, N * nu))))
        ))

    if stateLim or inputLim:
        nLim = 0
        if stateLim:
            nLim += (N+1)*nx
        if inputLim:
            nLim += N*nu
        Aineq = sparse.eye(nLim)
        # Put together eq and ineq
        A = sparse.vstack([Aeq, Aineq])
    else:
        A = Aeq

    # Add additional constraints for polyhedra
    if polyBlocks is not None:
        Cpolyx = np.zeros((0,nx))  # will stack rows below for each constraint
        for polyBlock in polyBlocks:
            # First get the LHS matrix corresponding to a single xk
            xstart, Nc, Ncx = polyBlock
            Cconstraint = np.tile(np.hstack([np.zeros(xstart), np.ones(Ncx), np.zeros(nx-(Ncx + xstart))]), (Nc,1))
            Cpolyx = sparse.vstack((Cpolyx, Cconstraint))
        A = sparse.vstack((A,
            # left block corresponds to x0,...,xN, right block to u0,...,u(N-1)
            sparse.hstack((sparse.kron(sparse.eye(N+1), Cpolyx), np.zeros(((N+1)*Cpolyx.shape[0], N*nu))))
        ))

    return A.tocsc()

def updateElems(matrix, row_index, col_index, vals):
    """Will update len(vals) elements; do nothing if that entry is zero. Returns 0 if nothing was updated (element was 0 in the sparse matrix), or 1 if updated successfully.
    Derived from https://rushter.com/blog/scipy-sparse-matrices/"""
    # Get col values
    col_start = matrix.indptr[col_index]
    col_end = matrix.indptr[col_index + 1]
    col_values = matrix.data[col_start:col_end]

    # Get row indices of occupied values
    index_start = matrix.indptr[col_index]
    index_end = matrix.indptr[col_index + 1]

    # contains indices of occupied cells at a specific col
    col_indices = list(matrix.indices[index_start:index_end])

    # Find a positional index for a specific row index
    try:
        i = col_indices.index(row_index)
        # Keep updating the col_values until the col ran out, or the user did not request any more
        Ntoupdate = min(len(vals), len(col_values)-i)
        for j in range(Ntoupdate):
            col_values[i + j] = vals[j]
        return Ntoupdate
    except ValueError:
        # Not found
        return 0

# def updateElem(obj, i, j, val):
#     '''
#     Returns 0 if nothing was updated (element was 0 in the sparse matrix), or 1 if updated successfully.
#     '''

#     # indptr has #cols elements, and the entry is the index into data/indices for the first element of that column
#     offs = obj.indptr[j]
#     #
#     if j < obj.shape[1] - 1:
#         maxOffsThisCol = obj.indptr[j+1]
#     else:
#         maxOffsThisCol = len(obj.indices)
#     while offs < maxOffsThisCol:
#         if obj.indices[offs] == i:
#             obj.data[offs] = val
#             return 1
#         offs += 1
#     return 0

def updateDynamics(obj, N, ti, Ad=None, Bd=None):
    '''Pass a block to update, and a matrix to go in that block, and it will update all the elements in that block.
    At ti=0, this updates the block equation x1 = A0 x0 + B0 u0 [+ c0],
    at ti=N-1, x[N] = A[N-1] x[N-1] + B[N-1] u[N-1] [+ c[N-1]]
    (recall X = (x0, ..., x[N], u0, ..., u[N-1]))
    Note that since this only deals with the constrain "A" matrix, the affine part ci must be updated separately.
    '''
    assert ti < N
    retval = 0
    if Ad is not None:
        nx = Ad.shape[0]
        for j in range(nx):
            retval += updateElems(obj, nx * (ti + 1), nx * ti + j, Ad[:,j])

    if Bd is not None:
        nx, nu = Bd.shape
        for j in range(nu):
            retval += updateElems(obj, nx * (ti + 1), (N + 1) * nx + nu * ti + j, Bd[:,j])

    return retval

def updatePolyBlock(obj, nx, nu, N, ti, polyBlocks, pbi, Cdi):
    '''Updates a block in the lower left (Aineq) with a matrix denoting a polyhedron.
    polyBlocks = Provide the same polyBlocks used in csc.init
    pbi = index into polyBlocks list
    Cdi = (Nc,Ncx)-shaped matrix to use to replace the existing block
    '''
    if polyBlocks is None:
        return
    assert ti <= N
    assert pbi < len(polyBlocks), "index too big for polyBlocks list"
    NcPerTi = 0
    NcBeforei = 0
    # Get information needed about all the constraints
    for ii in range(len(polyBlocks)):
        xstart, Nc, Ncx = polyBlocks[ii]
        NcPerTi += Nc
        if ii < pbi:
            NcBeforei += Nc
        elif ii == pbi:
            xstartj = xstart
            Nci = Nc
            Ncxi = Ncx

    # Lots of error checking
    assert Cdi.shape[0] == Nci, "Cdi shape is wrong"
    assert Cdi.shape[1] == Ncxi, "Cdi shape is wrong"

    ioffs = 2*(N+1)*nx + N*nu  # stanard Aeq,Aineq size of condensed A before polyBlocks
    ioffs += ti * NcPerTi + NcBeforei
    joffs = nx * ti + xstartj

    for j in range(Ncxi):
        updateElems(obj, ioffs, joffs + j, Cdi[:, j])

    return ioffs  # to be helpful

'''
These are all for the condensed A matrix
CSC matrix
three NumPy arrays: indices, indptr, data
In OSQP: sparse structure cannot be changed, i.e. can only change data
Ax_new: Vector of new elements in A->x
Ax_new_idx: Index mapping new elements to positions in A->x
A->x must be the data vector
Assuming that indices, indptr can't be changed
'''

class CondensedA:

    def __init__(self, nx, nu, N, polyBlocks=None, val=0):
        # Pass in nx=#states, nu=#inputs, and N=prediction horizon
        self.nx = nx
        self.nu = nu
        self.N = N
        if polyBlocks is not None:
            if np.sum(polyBlocks) != nx:
                raise ValueError('Sum of elements of polyBlocks must be = nx')
            # each element >= 1
            # all elements integers

        # create the indptr and index arrays for sparse CSC format
        self.data = np.zeros((self.nnz))
        self.indices = np.zeros_like(self.data, dtype=np.int32)
        numcols = self.shape[1]
        self.indptr = np.zeros((numcols + 1), dtype=np.int32)

        # NOTE: in the code below, i refers to row index, j to col index
        # similar for bj for block col index

        # populate indptr
        si = 0
        for j in range(numcols):
            if j < self.N * self.nx:
                si += nx + 2
            elif j < (self.N + 1) * self.nx:
                si += 2
            else:
                si += self.nx + 1
            self.indptr[j+1] = si

        # populate row indices
        self._offs = -1 # internal val to count through the data vector
        for j in range(self.shape[1]):
            # In each column, the Aeq blocks come first
            if j < (self.N + 1) * self.nx:
                # Leftmost cols
                bj = int(np.floor(j / self.nx)) # block col index - goes from [0,N+1)
                # the block diagonal -I
                self.addNZ(j, j, -1)
                    # The sub-block-diagonal Ad under the -I -- only in the left column block
                if j < self.N * self.nx:
                    for i in range(self.nx * (bj + 1), self.nx * (bj + 2)):
                        self.addNZ(j, i, val)
            else:
                # Right block column
                bj = int(np.floor((j - (self.N + 1) * self.nx) / self.nu))# block col index for the right - goes from [0,N)
                # The Bd blocks
                for i in range(self.nx * (bj + 1), self.nx * (bj + 2)):
                    self.addNZ(j, i, val)

            # The identity in Aineq
            if polyBlocks is not None and j < (self.N + 1) * self.nx:
                # we are at column j
                bj = int(np.floor(j / self.nx)) # block col index - goes from [0,N+1)
                iWithinX = j - bj * self.nx
                i0 = 0
                raise NotImplementedError
                for polyBlock in polyBlocks:
                    for ii in range(polyBlock):
                        # self.addNZ(j, i0, 1)
                        i0 += 1
                # 	if iWithinX >= polyBlock[0] and iWithinX < polyBlock[1]:
                # 		# rows
            else:
                self.addNZ(j, (self.N + 1) * self.nx + j, 1)

            # offs = self.matrixIdxToOffset(i, j)
            # if offs >= 0:
            # 	self.indices[offs] = j

        # sparse structure is now fixed, but data can be updated

    def addNZ(self, j, i, val):
        self._offs += 1
        self.indices[self._offs] = i # row index
        self.data[self._offs] = val

    # def blockToMatrixIdx(self, ib, jb, i, j):
    # 	# Return matrix indices from the block index (ib, jb) and element index within the block (i, j)
    # 	if ib < 2*(N+1):
    # 		rowi = self.nx * ib + i
    # 	else:
    # 		rowi = 2 * (self.N+1) * self.nx + (ib - 2 * (self.N+1)) * self.nu + i

    # 	if jb < (N+1):
    # 		rowi = self.nx * jb + j
    # 	else:
    # 		rowi = (self.N+1) * self.nx + (jb - (self.N+1)) * self.nu + j

    # 	return rowi, coli

    def get_nnz(self):
        # For testing
        return self.N * self.nx * (self.nx + 2) + self.nx * 2 + self.N * self.nu * (self.nx + 1)
    def get_shape(self):
        # For testing
        return (2 * (self.N + 1) * self.nx + self.N * self.nu), (self.N + 1) * self.nx + self.N * self.nu

    nnz = property(get_nnz, None, None, "Number of nonzero entries")
    shape = property(get_shape, None, None, "Number of nonzero entries")

max_stars_repo_name: avikde/controlutils

max_issues_repo_name: avikde/controlutils

max_forks_repo_name: avikde/controlutils

reason: import scipy


text:
import torch
import numpy as np

from utilities import conversions

def test_to_channels_last():
    batch_size, seq_len, height, width, channels = 5, 5, 32, 32, 3
    tensor = torch.randn((batch_size, seq_len, channels, height, width))
    converted = conversions.to_channels_last(tensor)

    for b in range(batch_size):
        for s in range(seq_len):
            for h in range(height):
                for w in range(width):
                    for c in range(channels):
                        assert tensor[b, s, c, h, w] == converted[b, s, h, w, c]

def test_to_channels_first():
    batch_size, seq_len, height, width, channels = 5, 5, 32, 32, 3
    tensor = torch.randn((batch_size, seq_len, height, width, channels))
    converted = conversions.to_channels_first(tensor)

    for b in range(batch_size):
        for s in range(seq_len):
            for h in range(height):
                for w in range(width):
                    for c in range(channels):
                        assert tensor[b, s, h, w, c] == converted[b, s, c, h, w]

def test_concat_rgb():
    """Test that concat_rgb correctly reshapes the tensor by concatenating sequences along
    channel dimensions
    """
    batch_len = 2
    seq_len = 5
    channels = 3
    img_size = 32
    batch = torch.randn((batch_len, seq_len, channels, img_size, img_size))
    concatenated = conversions.concat_rgb(batch)

    for concat_seq, batch_seq in zip(concatenated, batch):
        expected = torch.empty((seq_len * channels, 32, 32))
        for i in range(32):
            for j in range(32):
                rgb = torch.empty(channels * seq_len)
                for s in range(seq_len):
                    rgb[s * channels + 0] = batch_seq[s, 0, i, j]
                    rgb[s * channels + 1] = batch_seq[s, 1, i, j]
                    rgb[s * channels + 2] = batch_seq[s, 2, i, j]
                expected[:, i, j] = rgb
        assert torch.equal(expected, concat_seq)

def test_batch_to_sequence():
    batch_size, seq_len, height, width, channels = 15, 10, 32, 32, 3
    batch = np.random.normal(size=(batch_size, seq_len, height, width,
                                   channels))
    sequence = conversions.batch_to_sequence(batch)

    for b in range(batch_size):
        for s in range(seq_len):
            assert np.array_equal(batch[b],
                                  sequence[b * seq_len:(b + 1) * seq_len])