text:
import logging

from deepthought.experiments.encoding.experiment_templates.base import GenericNNEncoderExperiment, NestedCVExperimentTemplate

log = logging.getLogger('deepthought.experiments')

class SiameseEncoderExperiment(GenericNNEncoderExperiment):

    _default_params = NestedCVExperimentTemplate._default_params.copy()
    _default_params['group_attribute'] = 'subject'

    def __init__(self,
                 job_id,
                 hdf5name,
                 fold_generator,
                 encoder_pipeline_factory,
                 **kwargs):

        super(SiameseEncoderExperiment, self).__init__(
            job_id, hdf5name, fold_generator, encoder_pipeline_factory, **kwargs)

    def build_pretrain_model(self, data_dict, hyper_params):
        """
        pretrain-method specific;
        constucts a siamese net;
        works with any network structure of the pipeline
        :param data_dict:
        :param hyper_params:
        :return:
        """
        from theano import tensor
        from blocks.model import Model

        # Note: this has to match the sources defined in the dataset
        indices = [tensor.ivector('{}_indices'.format(i)) for i in range(2)]

        pipeline = self.encoder_pipeline_factory.build_pipeline(
            input_shape=data_dict.get_value().shape, params=hyper_params)

        # compute feature represenation
        rep = [pipeline.apply(data_dict[indices[i]]) for i in range(2)]
        # for r in rep: print r.type

        # flatten representations
        rep = [r.flatten(ndim=2) for r in rep]
        # for r in rep: print r.type

        # TODO: choose based on hyper-params
        if hyper_params['siamese_distance'] == 'L1-norm':
            rval = abs(rep[0] - rep[1]).sum(axis=1)
        elif hyper_params['siamese_distance'] == 'euclidean':  # L2-norm
            # compute Euclidean distance
            rval = tensor.sqrt(((rep[0] - rep[1])**2).sum(axis=1))
        elif hyper_params['siamese_distance'] == 'dot-product':
            # alternative: negative dot-product distance
            rval = -(rep[0] * rep[1]).sum(axis=1)  # element-wise multiplication and row sum
        else:
            raise ValueError('Unsupported value for siamese_distance: {}'.format(hyper_params['siamese_distance']))
        rval = tensor.reshape(rval, (rval.shape[0], 1))

        return Model(rval)

    @staticmethod
    def pretrain(model, hyper_params, full_hdf5, full_meta, train_selectors, valid_selectors=None):
        """
        generic training method for siamese networks;
        works with any network structure
        :return:
        """
        from theano import tensor
        from deepthought.datasets.siamese import PairsIndexDataset

        train_pairs = PairsIndexDataset(full_hdf5, full_meta,
                                        train_selectors,
                                        targets_source=hyper_params['pretrain_target_source'],
                                        group_attribute=hyper_params['group_attribute'])

        if valid_selectors is not None:
            if hyper_params['use_ext_dataset_for_validation']:
                ext_selectors = train_selectors
            else:
                ext_selectors = None

            valid_pairs = PairsIndexDataset(full_hdf5, full_meta,
                                            valid_selectors,
                                            ext_selectors=ext_selectors,
                                            targets_source=hyper_params['pretrain_target_source'],
                                            group_attribute=hyper_params['group_attribute'])
        else:
            valid_pairs = None

        # Note: this has to match the sources defined in the dataset
        y = tensor.lvector('targets')

        margin = hyper_params['siamese_margin']  # 0  # -0.5 # FIXME: hyper-param
        d = model.outputs[0]
        cost = 0.5 * (y * d + (1 - y) * tensor.maximum(margin - d, 0))
        cost = tensor.mean(cost)

        return GenericNNEncoderExperiment.run_pretrain(model, hyper_params, cost, train_pairs, valid_pairs)

reason: from theano


text:
"""
"""

# add directories in src/ to path
# sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__),'..')))

import argparse
import os

import keras.backend.tensorflow_backend as ktf
import numpy as np
import tensorflow as tf
from keras.layers import Input
from sklearn.externals import joblib

from applications.config import get_spectralnet_config, get_siamese_config
from core import networks, costs
from core.data import build_siamese_data, load_spectral_data, decode_data, load_base_data, build_spectral_data, \
    load_data, load_siamese_data, embed_if_needed

# PARSE ARGUMENTS
from core.util import print_accuracy, get_session
from keras import backend as K

parser = argparse.ArgumentParser()
parser.add_argument('--gpu', type=str, help='gpu number to use', default='')
parser.add_argument('--gpu_memory_fraction', type=float, help='gpu percentage to use', default='0.8')
parser.add_argument('--dset', type=str, help='dataset to use', default='mnist')
args = parser.parse_args()

params = get_spectralnet_config(args)
params['train_set_fraction'] = 0.8
data = load_spectral_data(params['data_path'], args.dset)

#
# params_noencode = params.copy()
# params_noencode['use_code_space'] = False
# # base_data = load_data(params_noencode)
# origin_data = build_spectral_data(params_noencode, base_data)

ktf.set_session(get_session(args.gpu_memory_fraction))

x_train, y_train, x_val, y_val, x_test, y_test = data['spectral']['train_and_test']
x_train_unlabeled, y_train_unlabeled, x_train_labeled, y_train_labeled = data['spectral'][
    'train_unlabeled_and_labeled']
x_val_unlabeled, y_val_unlabeled, x_val_labeled, y_val_labeled = data['spectral']['val_unlabeled_and_labeled']

if 'siamese' in params['affinity']:
    siamese_data = load_siamese_data(get_siamese_config(args)['data_path'], args.dset)
    pairs_train, dist_train, pairs_val, dist_val = siamese_data['train_and_test']

# x = np.ones((10, 2))
# y = np.ones((10, 1))

# x = np.concatenate([x_train, np.ones((10, 2))], axis=0)
x = x_val[:100]
y = y_val[:100]
batch_sizes = {
    'Unlabeled': x.shape[0],
    'Labeled': x.shape[0],
    'Orthonorm': x.shape[0],
}

def run_predict(params):
    K.set_learning_phase(0)
    input_shape = x.shape[1:]

    y_labeled_onehot = np.empty((0, params['n_clusters']))

    # spectralnet has three inputs -- they are defined here
    inputs = {
        'Unlabeled': Input(shape=input_shape, name='UnlabeledInput'),
        'Labeled': Input(shape=input_shape, name='LabeledInput'),
        'Orthonorm': Input(shape=input_shape, name='OrthonormInput'),
    }
    y_true = tf.placeholder(tf.float32, shape=(None, params['n_clusters']), name='y_true')

    # Load Siamese network
    if params['affinity'] == 'siamese':
        siamese_input_shape = [params['n_clusters']]
        siamese_inputs = {
            'Unlabeled': Input(shape=siamese_input_shape, name='UnlabeledInput'),
            'Labeled': Input(shape=siamese_input_shape, name='LabeledInput'),
        }
        siamese_net = networks.SiameseNet(siamese_inputs, params['arch'], params.get('siam_reg'), y_true, params['siamese_model_path'])

    else:
        siamese_net = None

    # Load Spectral net
    spectralnet_model_path = os.path.join(params['model_path'], 'spectral_net')
    spectral_net = networks.SpectralNet(inputs, params['arch'],
                                        params.get('spec_reg'),
                                        y_true, y_labeled_onehot,
                                        params['n_clusters'], params['affinity'], params['scale_nbr'],
                                        params['n_nbrs'], batch_sizes,
                                        spectralnet_model_path,
                                        siamese_net, train=False
                                        )
    # get final embeddings
    W_tensor = costs.knn_affinity(siamese_net.outputs['A'], params['n_nbrs'], scale=None, scale_nbr=params['scale_nbr'])

    x_spectralnet = spectral_net.predict_unlabelled(x)
    W = spectral_net.run_tensor(x, W_tensor)
    print('x_spectralnet', x_spectralnet.shape)
    clustering_algo = joblib.load(os.path.join(params['model_path'], 'spectral_net', 'clustering_aglo.sav'))

    kmeans_assignments = clustering_algo.predict_cluster_assignments(x_spectralnet)
    y_spectralnet = clustering_algo.predict(x_spectralnet)
    print_accuracy(kmeans_assignments, y, params['n_clusters'])
    # x_dec = decode_data(x, params, params['dset'])
    return x_spectralnet, y_spectralnet, x_spectralnet, W

# RUN EXPERIMENT
x_spectralnet, y_spectralnet, x, W = run_predict(params)
if args.dset in ['cc', 'cc_semisup']:
    # run plotting script
    import plot_2d
    plot_2d.process_plot(x_spectralnet, y_spectralnet, x, y=None, params=params, W=W)

ext: py

max_stars_repo_name: myelinio/SpectralNet

max_issues_repo_name: myelinio/SpectralNet

max_forks_repo_name: myelinio/SpectralNet

reason: import numpy


text:
import logging
import os
from typing import List

import matplotlib.pyplot as plt
import numpy as np

from ..features.speed_profile_handling import compute_features
from ..utils.io import makedirs, read_json, write_json

def plot_profiles(df, metadata=None, save_dir='/tmp'):
    plt.figure(figsize=(15, 5))
    if metadata is not None:
        vid = metadata['VehicleID']
        rel = metadata['Relative Error (%)']
        com = metadata["Comment"]
        plt.title(
            f"Speed profiles - VehicleID {vid}, Relative Error = {rel}, reason: {com}")
    plt.plot(df["axle1RawSpeed"], label="axle1RawSpeed")
    plt.plot(df["axle2RawSpeed"], label="axle2RawSpeed")
    plt.plot(df["trainSpeed"], label="trainSpeed")
    plt.legend()
    plt.savefig(f'{save_dir}/profiles.jpg')
    return f'{save_dir}/profiles.jpg'

def plot_features(df, features=[], nominal_feature_name='is_nominal', metadata=None, save_dir='/tmp'):
    for f in features:
        plt.figure(figsize=(10, 7))
        if metadata is not None:
            vid = metadata['VehicleID']
            rel = metadata['Relative Error (%)']
            com = metadata["Comment"]
            plt.title(
                f"Speed profiles - VehicleID {vid}, Relative Error = {rel}, reason: {com}")
        else:
            plt.title(f)

        plt.hist(df[df[nominal_feature_name] == True][f],  bins=40, density=True, range=(
            df[f].min()-1e-6, df[f].max()+1e-6), alpha=0.5, label="nominal")
        plt.hist(df[df[nominal_feature_name] == False][f], bins=40, density=True, range=(
            df[f].min()-1e-6, df[f].max()+1e-6), alpha=0.5, label="anomalous")
        plt.legend()
        plt.savefig(f'{save_dir}/features.jpg')
    return f'{save_dir}/features.jpg'

def plot_features_grid(df, features=[], vehicle_id: str = "2048", nominal_feature_name='is_nominal', metadata=None, log=False, save_dir='/tmp'):
    if isinstance(vehicle_id, int):
        vehicle_id = str(vehicle_id)

    one_vehicle_df = df[[vehicle_id in f for f in df.index]]
    df = one_vehicle_df

    features = df.columns.drop(nominal_feature_name)
    n_features = len(features)

    n_cols = 5
    nrows = n_features//n_cols if (n_features %
                                   n_cols == 0) else n_features//n_cols + 1
    fig, axs = plt.subplots(nrows, n_cols, figsize=(n_cols*5, nrows*4))
    axs = axs.ravel()

    for i, f in enumerate(features):
        if metadata is not None:
            vid = metadata['VehicleID']
            rel = metadata['Relative Error (%)']
            com = metadata["Comment"]
            axs[i].set_title(
                f"Speed profiles - VehicleID {vid}, Relative Error = {rel}, reason: {com}")
        else:
            axs[i].set_title(f)

        xn, binsn, pn = axs[i].hist(df[df[nominal_feature_name] == True][f],  bins=40, density=True, range=(
            df[f].min()-1e-6, df[f].max()+1e-6), alpha=0.5, label="nominal")
        xa, binsa, pa = axs[i].hist(df[df[nominal_feature_name] == False][f], bins=40, density=True, range=(
            df[f].min()-1e-6, df[f].max()+1e-6), alpha=0.5, label="anomalous")

        for item in pn:
            item.set_height(item.get_height()/np.max(xn))
        for item in pa:
            item.set_height(item.get_height()/np.max(xa))
        if log:
            axs[i].set_yscale('log')
            axs[i].set_ylim(0, 10)
        else:
            axs[i].set_ylim(0, 1.2)

    plt.legend()
    plt.savefig(f'{save_dir}/features_grid.jpg')
    return f'{save_dir}/features_grid.jpg'

def plot_event(event_filename: str,
                  bucket_name: str,
                  minio_client,
                  features=[],
                  sample_value: str = "500L",
                  speed_vars: List[str] = None,
                  save_dir='/tmp'):
    if speed_vars is None:
        speed_vars = ["axle1RawSpeed", "axle2RawSpeed", "trainSpeed"]
    if isinstance(sample_value, int):
        sample_value = str(str(sample_value) + 'L')

    df_event = minio_client.get_pickle(event_filename, bucket_name)
    df_event = df_event.resample(sample_value).mean()

    for v in speed_vars:
        df_event[v] = df_event[v].astype(np.int16)

    plot_profiles(df_event)
    print(df_event["axle1RawSpeed"].diff().dropna().mean())
#     df_event["axle1RawSpeed"].diff().dropna().plot()
#     df_event["axle1RawSpeed"].plot()

    features = compute_features(df_event, *speed_vars)
    # pprint(features)

def plot_all_metrics(recalls, precisions, f1s, thresholds, bt, save_dir='/tmp'):
    plt.plot(thresholds, precisions, alpha=0.5, label="precision")
    plt.plot(thresholds, recalls,    alpha=0.5, label="recall")
    plt.plot(thresholds, f1s, label="f1")
    plt.axvline(bt, color="red", linewidth=3)

    plt.title(f"Best threshold for anomaly detection: {bt}")

    plt.legend()
    plt.savefig(f'{save_dir}/metrics.jpg')
    logging.info(f"prf metrics saved !")
    return f'{save_dir}/metrics.jpg'

def plot_model_acc(model, save_dir='/tmp'):
    plt.figure()
    plt.plot(model["val_axle_output_accuracy"], label="accuracy axle event")
    plt.plot(model["val_odo_output_accuracy"], label="accuracy odometry algo issue")
    plt.plot(model["val_speed_output_accuracy"], label="accuracy speed estimation")
    # plt.ylim(0.9,1.01)
    plt.legend()
    plt.savefig(f'{save_dir}/accuracy.jpg')
    logging.info(f"model accuracy figure saved !")
    return f'{save_dir}/accuracy.jpg'

def plot_val_loss(model, save_dir='/tmp'):
    plt.figure()
    plt.plot(model["val_axle_output_loss"], label="validation loss axle event")
    plt.plot(model["val_odo_output_loss"], label="validation loss odometry algo issue")
    plt.plot(model["val_speed_output_loss"], label="validation loss speed estimation")
    plt.legend()
    plt.semilogy()
    plt.savefig(f'{save_dir}/val_loss.jpg')
    logging.info(f"validation loss figure saved !")
    return f'{save_dir}/val_loss.jpg'

def plot_train_val_loss(model, save_dir='/tmp'):
    plt.figure()
    plt.plot(model["axle_output_loss"], label="train loss axle event")
    plt.plot(model["val_axle_output_loss"], label="validation loss axle event")
    plt.legend()
    plt.semilogy()
    plt.savefig(f'{save_dir}/axle_train_val_loss.jpg')
    logging.info(f"axle train loss figure saved !")

    plt.figure()
    plt.plot(model["odo_output_loss"], label="train loss odometry algo issue")
    plt.plot(model["val_odo_output_loss"], label="validation loss odometry algo issue")
    plt.legend()
    plt.semilogy()
    plt.savefig(f'{save_dir}/odo_train_val_loss.jpg')
    logging.info(f"odo train loss figure saved !")

    plt.figure()
    plt.plot(model["speed_output_loss"], label="train loss speed estimation")
    plt.plot(model["val_speed_output_loss"], label="validation loss speed estimation")
    plt.legend()
    plt.semilogy()
    plt.savefig(f'{save_dir}/speed_train_val_loss.jpg')
    logging.info(f"speed train loss figure saved !")
    return f'{save_dir}/axle_train_val_loss.jpg', f'{save_dir}/odo_train_val_loss.jpg', f'{save_dir}/speed_train_val_loss.jpg'

def plot_train_valid_global_loss(model, save_dir='/tmp'):
    plt.figure()
    plt.plot(model["loss"], label="train loss global")
    plt.plot(model["val_loss"], label="validation loss global")
    plt.legend()
    plt.semilogy()
    plt.savefig(f'{save_dir}/train_valid_global_loss.jpg')
    logging.info(f"global train loss figure saved !")
    return f'{save_dir}/train_valid_global_loss.jpg'

def get_all_plots(history, save_dir):
    save_figs_dir = os.path.join(save_dir, 'figures')
    makedirs(save_figs_dir)
    plt_model_acc = plot_model_acc(history, save_dir=save_figs_dir)
    plt_val_loss = plot_val_loss(history, save_dir=save_figs_dir)
    plt_axle_loss, plt_odo_loss, plt_speed_loss = plot_train_val_loss(
        history, save_dir=save_figs_dir)
    plt_train_valid_global_loss = plot_train_valid_global_loss(
        history, save_dir=save_figs_dir)
    figures = [plt_model_acc, plt_val_loss, plt_axle_loss,
            plt_odo_loss, plt_speed_loss, plt_train_valid_global_loss]
    return figures


text:
# -*- coding: utf8 -*-
import os

import tensorflow as tf
from tensorflow.examples.tutorials.mnist import input_data
import numpy as np
import matplotlib.pyplot as plt
import mnist_inference

# define input, output, batch and training params

BATCH_SIZE = 50
LEARNING_RATE_BASE = 0.8
LEARNING_RATE_DECAY = 0.99
REGULARIZATION_RATE = 0.0001
TRAINING_STEPS = 10000
MOVING_AVERAGE_DECAY = 0.99

MODEL_SAVE_PATH = "model/"
MODEL_NAME = "model.ckpt"
score_filename = "accuracy_score_cnn.txt"

# train a convolutional neural network
def train(mnist, continue_train=False):
    x = tf.placeholder(tf.float32, [BATCH_SIZE,
                                    mnist_inference.IMAGE_SIZE,
                                    mnist_inference.IMAGE_SIZE,
                                    mnist_inference.NUM_CHANNELS], name='x-input')
    y_ = tf.placeholder(tf.float32, [None, mnist_inference.OUTPUT_NODE], name='y-input')

    regularizer = tf.contrib.layers.l2_regularizer(REGULARIZATION_RATE)
    y = mnist_inference.inference(x, True, regularizer)
    global_step = tf.Variable(0, trainable=False)

    # moving average, cross entropy, loss function with regularization and learning rate
    variable_average = tf.train.ExponentialMovingAverage(MOVING_AVERAGE_DECAY, global_step)
    variable_average_op = variable_average.apply(tf.trainable_variables())
    cross_entropy = tf.nn.sparse_softmax_cross_entropy_with_logits(logits=y, labels=tf.argmax(y_, 1))
    cross_entropy_mean = tf.reduce_mean(cross_entropy)
    loss = cross_entropy_mean + tf.add_n(tf.get_collection('losses'))
    learning_rate = tf.train.exponential_decay(
        LEARNING_RATE_BASE,
        global_step,
        mnist.train.num_examples / BATCH_SIZE,
        LEARNING_RATE_DECAY
    )

    train_step = tf.train.GradientDescentOptimizer(learning_rate).minimize(loss, global_step=global_step)
    with tf.control_dependencies([train_step, variable_average_op]):
        train_op = tf.no_op(name='train')

    # initialize persistence class
    saver = tf.train.Saver()

    config = tf.ConfigProto(allow_soft_placement=True)
    gpu_options = tf.GPUOptions(per_process_gpu_memory_fraction=0.4)
    config.gpu_options.allow_growth = True
    with tf.Session(config=config) as sess:
        if continue_train:
            ckpt = tf.train.get_checkpoint_state(
                MODEL_SAVE_PATH
            )
            if ckpt and ckpt.model_checkpoint_path:
                saver.restore(sess, ckpt.model_checkpoint_path)
        else:
            sess.run(tf.global_variables_initializer())
        # create directory
        try:
            os.mkdir(MODEL_SAVE_PATH)
        except:
            print("directory already exist")

        # define accuracy
        correct_prediction = tf.equal(tf.arg_max(y, 1), tf.arg_max(y_, 1))
        accuracy = tf.reduce_mean(tf.cast(correct_prediction, tf.float32))
        test_result = list(range(int(TRAINING_STEPS / 1000)))

        for i in range(TRAINING_STEPS):
            xs, ys = mnist.train.next_batch(BATCH_SIZE)
            reshaped_xs = np.reshape(xs, (
                BATCH_SIZE,
                mnist_inference.IMAGE_SIZE,
                mnist_inference.IMAGE_SIZE,
                mnist_inference.NUM_CHANNELS))

            _, loss_value, step = sess.run([train_op, loss, global_step], feed_dict={x: reshaped_xs, y_: ys})

            txs = mnist.test.images[0:BATCH_SIZE]
            test_feed = {
                x: np.reshape(txs, (BATCH_SIZE,
                                    mnist_inference.IMAGE_SIZE,
                                    mnist_inference.IMAGE_SIZE,
                                    mnist_inference.NUM_CHANNELS)),
                y_: mnist.test.labels[0:BATCH_SIZE]}

            accuracy_score = sess.run(accuracy, feed_dict=test_feed)
            test_result[int(i / 1000)] = accuracy_score

            if i % 1000 == 0:
                print("after %d training step(s), loss on training batch is %g , validation accuracy = %g" % (
                    step, loss_value, accuracy_score))
                saver.save(
                    sess, os.path.join(MODEL_SAVE_PATH, MODEL_NAME), global_step=global_step
                )

    # np.savetxt(score_filename, test_result, fmt="%0.4f")
    #
    # dispImg(test_result, 'accuracy_score')
    # plt.show()

def dispImg(test_result, filename):
    # draw a graph of accuracy using matplotlib
    iteration_count = range(0, TRAINING_STEPS, 1000)
    plt.figure(num=1, figsize=(15, 8))
    plt.title("Plot accuracy", size=20)
    plt.xlabel("iteration count", size=14)
    plt.ylabel("accuracy/%", size=14)
    test_note = [TRAINING_STEPS - 1000, test_result[TRAINING_STEPS / 1000 - 1]]
    plt.annotate('test-' + str(test_note), xy=(test_note[0], test_note[1]),
                 xytext=(test_note[0] + 1000, test_note[1] - 0.07), arrowprops=dict(facecolor='black', shrink=0.05))
    plt.grid(True)
    plt.plot(iteration_count, test_result, linestyle='-.', marker='X', label='test data')
    plt.legend(loc="upper left")
    try:
        os.mkdir('images/')
    except:
        print("directory already exist")
    plt.savefig('images/%s.png' % filename, format='png')

def main(argv=None):
    mnist = input_data.read_data_sets("../MNIST_data", one_hot=True)
    print("start")
    train(mnist, True)

if __name__ == '__main__':
    tf.app.run()

max_stars_repo_name: carbo-T/TF

max_issues_repo_name: carbo-T/TF

max_forks_repo_name: carbo-T/TF


text:
from PyQt5 import QtWidgets
import numpy as np

def read_measurements(self):
    from f_read_measurements_excel import f_read_measurements_excel
    filename = QtWidgets.QFileDialog.getOpenFileName()
    try:
        pathnamemeasurement = str(filename[0])
        PathForMeasurements.setText(self.pathnamemeasurement)
        measurements = f_read_measurements_excel(self.pathnamemeasurement)
        size_of_var = np.shape(self.measurements)
        self.number_measurements = int(size_of_var[1])
        content_to_print = ' Number of measurements uploaded is: {}.'.format(self.number_measurements)
        MeasFileContents.setText(str(content_to_print))
    except:
        PathForMeasurements.setText("Enter path for file containing measurement data")
        MeasFileContents.clear()
    from f_check_inputs import f_check_inputs
    checkinputs = f_check_inputs(measurements, ims_predictions,
                                      combined_uncertainties, checkboxvariable, numberofstage,
                                      numberofmeasurementperstage)
    print(checkinputs)
    return measurements

max_stars_repo_name: MeDIUM-FCL/MeDIUM

max_issues_repo_name: MeDIUM-FCL/MeDIUM

max_forks_repo_name: MeDIUM-FCL/MeDIUM


text:
""" Session defined for the IDGAN system. Includes the training process, a quick evaluation method that
uses sampled test items, and a more comprehensive evaluation on the entirety of the test set. Significant
cross-over with SAE session, as the evaluated system component is the IDGAN generator, instantiated as an SAE. """

import os
import time
import codecs
import logging
import numpy as np
import tensorflow as tf

from shared.util import make_pickle, load_pickle, load_model
from autoencoder.src.codebase.preprocessing import prepare_data
from autoencoder.src.codebase.batching import DataServer

from IDGAN.src.codebase.seq_gan import IDGAN
from IDGAN.src.codebase.gan_trainer import IDGANTrainer
from IDGAN.src.codebase.interface import IDGANInterface

# Options
from cognitive_language_model.src.options.test_options import TestOptions as LMTestOpt
from sentence_similarity_classifier.src.options.test_options import TestOptions as SSCTestOpt
from autoencoder.src.options.train_options import TrainOptions as AETrainOpt
from IDGAN.src.options.train_options import TrainOptions as GANTrainOpt
# from IDGAN.src.options.small_options import TrainOptionsSmall as GANTrainOptSmall
from IDGAN.src.options.test_options import TestOptions as GANTestOpt

# Suppress TensorFlow info-level messages
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '1'

# Initialize options (uncomment the preferred set of IDGAN training options)
# Component options are initialized from component-specific option classes
lm_gan_opt = LMTestOpt()
ssc_gan_opt = SSCTestOpt()
ae_gan_opt = AETrainOpt()
train_opt = GANTrainOpt()
# train_opt = GANTrainOptSmall()
test_opt = GANTestOpt()
opts = [train_opt, ae_gan_opt, lm_gan_opt, ssc_gan_opt]

# Training is logged to allow for post-hoc examination of possible error sources
logging_path = os.path.join(train_opt.out_dir, 'console_log.txt')
logging.basicConfig(filename=logging_path, level=logging.INFO)
# The logs are also printed to console for immediate evaluation
console = logging.StreamHandler()
console.setLevel(logging.INFO)
logging.getLogger('').addHandler(console)

# Specify training, validation, and test sets
if train_opt.use_toy:
    # Toy sets for quick identification of syntactic and semantic errors
    data_dir = os.path.join(train_opt.data_dir, 'toy')
    source_train_name, source_valid_name, source_test_name = 'toy_train', 'toy_valid', 'toy_test'
    target_train_name, target_valid_name, target_test_name = 'toy_train', 'toy_valid', 'toy_test'
else:
    # ID-variant Europarl corpora provide the input data for IDGAN training
    data_dir = os.path.join(train_opt.data_dir, 'europarl')
    source_train_name, source_valid_name, source_test_name = 'europarl_v7_high_train', 'europarl_v7_high_valid', \
                                                             'europarl_v7_high_test'
    target_train_name, target_valid_name, target_test_name = 'europarl_v7_low_train', 'europarl_v7_low_valid', \
                                                             'europarl_v7_low_test'

corpus_names = [source_train_name, source_valid_name, source_test_name, target_train_name, target_valid_name,
                target_test_name]
# Construct paths pointing to source files
source_paths = [os.path.join(data_dir, '{:s}.txt'.format(name)) for name in corpus_names]
# Construct paths pointing to pickles containing corpus items and vocabulary
pickle_paths = [os.path.join(data_dir, '{:s}_data.pkl'.format(name)) for name in corpus_names]
# Create pickles - this way, data pre-processing has only to be done once
if not os.path.isfile(pickle_paths[0]):
    for i in range(len(corpus_names)):
        make_pickle(train_opt, prepare_data, corpus_names[i], source_paths[i], pickle_paths[i], vocab_path=None,
                    is_train='train' in corpus_names[i], is_valid='valid' in corpus_names[i],
                    is_test='test' in corpus_names[i])

# Load the vocabulary pickle shared among all of IDGAN's components
if train_opt.use_toy:
    vocab_pickle = os.path.join(train_opt.root_dir, 'data/toy/toy_train_vocab.pkl')
else:
    vocab_pickle = os.path.join(train_opt.root_dir, 'data/europarl/europarl_v7_train_vocab.pkl')
vocab = load_pickle(vocab_pickle)

# Write vocabulary contents to file for manual inspection
vocab_log_path = os.path.join(train_opt.out_dir, '{:s}_vocab_log.txt'.format('shared'))
with codecs.open(vocab_log_path, 'w', encoding='utf8') as in_file:
    for key, value in vocab.index_to_word.items():
        in_file.write('{:s}, {:s}\n'.format(str(key), value))
print('Vocab log written.')

# Define TensorFlow session configuration - active during all calls to session.py (uncomment desired settings)
config = tf.ConfigProto(allow_soft_placement=True)

# config.gpu_options.allow_growth = True
# config.gpu_options.per_process_gpu_memory_fraction = 0.5

# ======================================================================================================

def get_ckpt_vars(ckpt_dir):
    """ Returns the names of variables stored within a TensorFlow checkpoint; used in the initialization
    of IDGAN's components with pre-trained parameters. """
    ckpt_vars = list()
    with tf.Session():
        for var_name, _ in tf.contrib.framework.list_variables(ckpt_dir):
            ckpt_vars.append(var_name)
    return ckpt_vars

# ======================================================================================================

def train_session():
    """ Executes a training session on the IDGAN system. """
    # Clear the default graph within which the model graph is constructed
    tf.reset_default_graph()

    # Load data
    source_train_data = load_pickle(pickle_paths[0])
    source_valid_data = load_pickle(pickle_paths[1])
    target_train_data = load_pickle(pickle_paths[3])
    target_valid_data = load_pickle(pickle_paths[4])

    # Construct the system graph (component-specific graphs are constructed within the IDGAN graph)
    seq_gan = IDGAN(opts, vocab, 'IDGAN')

    # Initialize IDGAN's component models with pre-trained parameters
    # Declare paths pointing to checkpoints containing desired parameter values
    component_ckpt_dir = os.path.join(train_opt.local_dir, 'checkpoints/components')
    lm_dir = os.path.join(component_ckpt_dir, 'lm')
    source_encoder_dir = os.path.join(component_ckpt_dir, 'source')
    if train_opt.cross_dec:
        source_decoder_dir = os.path.join(component_ckpt_dir, 'source_decoder')
    else:
        source_decoder_dir = os.path.join(component_ckpt_dir, 'source')  # NO crossing
    target_dir = os.path.join(component_ckpt_dir, 'target')
    chosen_epoch = 'best'

    # Isolate parameters to be loaded into the IDGAN's system
    # Excludes optimization variables as well as variables connected to the training of 'frozen' IDGAN components
    # Get lists of variables contained within component checkpoint files
    lm_vars_plus_optimization = get_ckpt_vars(lm_dir)
    source_encoder_vars_plus_optimization = get_ckpt_vars(source_encoder_dir)
    source_decoder_vars_plus_optimization = get_ckpt_vars(source_decoder_dir)
    target_vars_plus_optimization = get_ckpt_vars(target_dir)
    # Exclude training-specific variables from IDGAN initialization
    lm_vars = [var_name for var_name in lm_vars_plus_optimization if 'optimization' not in var_name]
    # To enable the 'crossed decoder' training condition, separate the encoder and decoder variables
    # of the SAE pre-trained on the high-ID corpus ('translator SAE' within IDGAN)
    source_encoder_vars = \
        [var_name for var_name in source_encoder_vars_plus_optimization if 'optimization' not in var_name]
    source_encoder_vars = [var_name for var_name in source_encoder_vars if 'encoder' in var_name]
    source_decoder_vars = \
        [var_name for var_name in source_decoder_vars_plus_optimization if 'optimization' not in var_name]
    source_decoder_vars = [var_name for var_name in source_decoder_vars if 'decoder' in var_name]
    target_vars = [var_name for var_name in target_vars_plus_optimization if 'optimization' not in var_name]
    # Obtain list of all variables which have to be initialized (either randomly or from pre-trained values)
    # within the IDGAN system
    all_vars = tf.get_collection(tf.GraphKeys.GLOBAL_VARIABLES)

    # Check for matches between variables found within the pre-trained checkpoints and IDGAN variables
    lm_parameters = [var for var in all_vars if var.name.split(':')[0] in lm_vars]
    source_encoder_parameters = [var for var in all_vars if var.name.split(':')[0] in source_encoder_vars]
    source_decoder_parameters = [var for var in all_vars if var.name.split(':')[0] in source_decoder_vars]
    target_parameters = [var for var in all_vars if var.name.split(':')[0] in target_vars]
    # Load matching variables from corresponding checkpoints
    loaded_parameters = lm_parameters + source_encoder_parameters + source_decoder_parameters + target_parameters
    # Rest is initialized randomly
    initialized_parameters = [var for var in all_vars if var not in loaded_parameters]

    # Initialize saver objects tasked with loading in the pre-trained parameters
    lm_saver = tf.train.Saver(lm_parameters)
    source_encoder_saver = tf.train.Saver(source_encoder_parameters)
    source_decoder_saver = tf.train.Saver(source_decoder_parameters)
    target_saver = tf.train.Saver(target_parameters)
    # Declare random initialization OP
    init_op = tf.variables_initializer(initialized_parameters)

    # Time training duration
    starting_time = time.time()

    with tf.Session(config=config) as train_sess:
        # Load pre-trained parameters into the IDGAN graph
        load_model(train_sess, lm_saver, lm_dir, chosen_epoch)
        load_model(train_sess, source_encoder_saver, source_encoder_dir, chosen_epoch)
        load_model(train_sess, source_decoder_saver, source_decoder_dir, chosen_epoch)
        load_model(train_sess, target_saver, target_dir, chosen_epoch)
        # Initialize the rest
        train_sess.run(init_op)

        # Initialize IDGAN interface and trainer, used for inference and training steps, respectively
        interface = IDGANInterface(seq_gan, vocab, train_sess, test_opt)
        trainer = IDGANTrainer(vocab, train_opt, seq_gan, train_sess, source_train_data, source_valid_data,
                               target_train_data, target_valid_data, test_opt, interface, verbose=True)
        # Train system (either for a predefined number of epochs or until early stopping)
        print('+++TRAINING+++')
        trainer.train_gan()

    # Report training duration
    elapsed = time.time() - starting_time
    logging.info('Training took {:d} hours, {:d} minutes, and {:.4f} seconds.'.format(
        int(elapsed // 3600), int((elapsed % 3600)) // 60, elapsed % 60))

# ======================================================================================================

def test_session(batch_size=1, target_epoch='best', beam_decoding=False):
    """ Executes a quick test session on the IDGAN system by sampling a small quantity of items from the test set
    and using the model to first compress them into a meaningful representation and subsequently reconstruct them;
    the evaluation process focuses exclusively on the translator SAE. """
    # Clear the default graph
    tf.reset_default_graph()
    # Declare the batch size, if left unspecified in test options
    if train_opt.batch_size is None:
        train_opt.batch_size = batch_size
    # Load data
    source_test_data = load_pickle(pickle_paths[2])
    # Build system graph
    seq_gan = IDGAN(opts, vocab, 'IDGAN')
    # Declare saver object for restoring learned IDGAN parameters
    test_saver = tf.train.Saver()

    # Initiate testing session
    with tf.Session(config=config) as test_sess:
        # Load learned model parameters
        load_model(test_sess, test_saver, test_opt.save_dir, target_epoch)
        # Initialize system interface containing inference methods
        interface = IDGANInterface(seq_gan, vocab, test_sess, test_opt)
        # Sample candidate sentences from the test set
        samples = np.random.choice(source_test_data, test_opt.num_samples).tolist()
        while max([len(sample.split()) for sample in samples]) > 10:
            samples = np.random.choice(source_test_data, test_opt.num_samples).tolist()
        # Initialize a loader object to pre-process the sampled sentences
        sample_loader = DataServer(samples, vocab, test_opt)
        samples_read = 0

        print('Sampled sentences:')
        for s_id, s in enumerate(samples):
            print('{:d}: {:s}'.format(s_id, s))
        print('-' * 10 + '\n')

        if not beam_decoding:
            # Perform greedy ID-reduction on the sampled sentences
            print('Greedy decoding:')
            for _, sample_data in enumerate(sample_loader):
                enc_labels, enc_inputs, dec_inputs = sample_data
                generated = interface.greedy_generation(enc_labels, enc_inputs, dec_inputs)
                for j in range(test_opt.batch_size):
                    print('Encoded: {:s}\nDecoded: {:s}\n'.format(samples[samples_read + j], generated[j]))
                samples_read += test_opt.batch_size
        else:
            # Perform greedy ID-reduction with beam-search on the sampled sentences
            assert (test_opt.batch_size == 1), 'Beam search not defined for batches with more than one element.'
            print('Beam search decoding:')
            for _, sample_data in enumerate(sample_loader):
                enc_labels, enc_input, _ = sample_data
                print('Encoded: {:s}'.format(samples[i]))
                interface.beam_generation(enc_labels, enc_input, print_results=True)

    print('-' * 10 + '\n')
    print('=' * 10 + '\n')
    print('IDGAN evaluation completed!')

# ======================================================================================================

def test_to_file(batch_size=1, target_epoch='best', beam_decoding=True):
    """ Executes a comprehensive test session on the entire test corpus;
    output is written to file for the calculation of the achieved corpus-wide ID reduction and
    the BLEU score between source sentences and their ID-reduced translations. """
    def _reconstruct_input(input_array):
        """ Reconstructs input sentences from numpy arrays; used to derive an accurate representation of the
        pre-processed, encoded sequences. """
        # Convert input array to list of lists of word indices
        input_idx = [np.squeeze(array).tolist() for array in np.split(input_array, input_array.shape[0], axis=0)]
        # Translate indices into corresponding word tokens; truncated after sentence-final <EOS>
        input_boundaries = [idx_list.index(vocab.eos_id) if vocab.eos_id in idx_list else len(idx_list)
                            for idx_list in input_idx]
        input_sentences = [[vocab.index_to_word[idx] for idx in input_idx[j][:input_boundaries[j]]]
                           for j in range(len(input_idx))]
        input_sentences = [' '.join(word_list) + '.' for word_list in input_sentences]
        return input_sentences

    assert (test_opt.batch_size == 1), \
        'Function is defined for a batch size of 1 due to the nature of beam search implementation.'

    # Clear the default graph
    tf.reset_default_graph()
    # Declare the batch size
    if train_opt.batch_size is None:
        train_opt.batch_size = batch_size
    # Load test data from the high-ID corpus
    source_test_data = load_pickle(pickle_paths[2])
    # Build model graph
    seq_gan = IDGAN(opts, vocab, 'IDGAN')
    # Declare saver object for restoring learned model parameters
    test_saver = tf.train.Saver()

    # Declare paths pointing to locations of output files (i.e. reference and translations sets for BLEU)
    encoded_path = os.path.join(test_opt.out_dir, 'source_encoded_test_corpus_beam_{:s}.txt'.format(str(beam_decoding)))
    decoded_path = os.path.join(test_opt.out_dir, 'source_decoded_test_corpus_beam_{:s}.txt'.format(str(beam_decoding)))

    # Initiate testing session
    with tf.Session(config=config) as test_sess:
        # Load learned model parameters
        load_model(test_sess, test_saver, test_opt.save_dir, target_epoch)
        # Initialize the model interface containing inference methods
        interface = IDGANInterface(seq_gan, vocab, test_sess, test_opt)
        # Initialize a loader object to pre-process the test corpus
        test_loader = DataServer(source_test_data, vocab, test_opt)

        with open(encoded_path, 'w') as enc_file:
            with open(decoded_path, 'w') as dec_file:
                if not beam_decoding:
                    # Perform greedy ID-reduction on the sampled sentences
                    print('Greedy decoding:')
                    for s_id, test_items in enumerate(test_loader):
                        enc_labels, enc_inputs, dec_inputs = test_items
                        generated = interface.greedy_generation(enc_labels, enc_inputs, dec_inputs)
                        enc_file.write(_reconstruct_input(enc_labels)[0] + '\n')
                        dec_file.write(generated[0] + '\n')
                        if s_id % 10 == 0 and s_id > 0:
                            print('{:d} sentences written to file.'.format(s_id))

                else:
                    # Perform greedy ID-reduction with beam-search on the sampled sentences
                    assert (test_opt.batch_size == 1), 'Beam search not defined for batches with more than one element.'
                    print('Beam search decoding:')
                    for s_id, test_items in enumerate(test_loader):
                        enc_labels, enc_input, _ = test_items
                        generated = interface.beam_generation(enc_labels, enc_input, print_results=False)
                        # Write best beam result only
                        enc_file.write(_reconstruct_input(enc_labels)[0] + '\n')
                        dec_file.write(generated[0][0] + '\n')
                        if s_id % 10 == 0 and s_id > 0:
                            print('{:d} sentences written to file.'.format(s_id))

    print('-' * 10 + '\n')
    print('=' * 10 + '\n')
    print('IDGAN documented evaluation completed!')

max_stars_repo_name: demelin/IDGAN

max_issues_repo_name: demelin/IDGAN

max_forks_repo_name: demelin/IDGAN

reason: import numpy,from numpy


text:
import joblib
import numpy as np

from flask import Flask
from flask import jsonify

app = Flask(__name__)

#POSTMAN TESTING
@app.route('/predict', methods=['GET'])
def predict():
    X_test = np.array([7.594444821,7.479555538,1.616463184,1.53352356,0.796666503,0.635422587,0.362012237,0.315963835,2.277026653])
    prediction = model.predict(X_test.reshape(1,-1))
    return jsonify({'prediccion' : list(prediction)})

if __name__ == "__main__":
    model = joblib.load('./models/best_model.pkl')
    app.run(port=8080)


text:
import numpy as np
import pandas as pd

def create_mesh(deck):
    edges = [np.linspace(x.start, x.end, deck.ctrl.cells_per_region + 1,
                         dtype=np.float32) for x in deck.reg.values()]
    edges = np.unique(np.concatenate(edges))
    mats = [None] * len(edges)
    sigma_a = np.zeros(len(edges), dtype=np.float32)
    sigma_s0 = np.zeros(len(edges), dtype=np.float32)
    sigma_s1 = np.zeros(len(edges), dtype=np.float32)
    for edge_index, edge in enumerate(edges):
        for reg in deck.reg.values():
            if (reg.start <= edge) and (edge <= reg.end):
                mats[edge_index] = reg.mat
                sigma_a[edge_index] = deck.mat[reg.mat].sigma_a
                sigma_s0[edge_index] = deck.mat[reg.mat].sigma_s0
                sigma_s1[edge_index] = deck.mat[reg.mat].sigma_s1
                break

    source = np.zeros(len(edges), dtype=np.float32)
    for edge_index, edge in enumerate(edges):
        for src in deck.src.values():
            if (src.start <= edge) and (edge <= src.end):
                source[edge_index] = src.magnitude
                break

    sigma_t = sigma_a + sigma_s0
    df = pd.DataFrame({'edge': edges,
                       'mat': mats,
                       'sigma_t': sigma_t,
                       'sigma_a': sigma_a,
                       'sigma_s0': sigma_s0,
                       'sigma_s1': sigma_s1,
                       'source': source})
    return df

max_stars_repo_name: pozulp/narrows

max_issues_repo_name: LLNL/narrows

max_forks_repo_name: LLNL/narrows


text:
from __future__ import absolute_import, division #makes KratosMultiphysics backward compatible with python 2.6 and 2.7

# Importing the Kratos Library
import KratosMultiphysics

# Import applications
import KratosMultiphysics.ConvectionDiffusionApplication as KratosConvDiff
import KratosMultiphysics.MultilevelMonteCarloApplication as KratosMLMC

# Avoid printing of Kratos informations
KratosMultiphysics.Logger.GetDefaultOutput().SetSeverity(KratosMultiphysics.Logger.Severity.WARNING) # avoid printing of Kratos things

# Importing the base class
from analysis_stage import AnalysisStage

# Import packages
import numpy as np

# Import Monte Carlo library
import mc_utilities as mc

# Import cpickle to pickle the serializer
try:
    import cpickle as pickle  # Use cPickle on Python 2.7
except ImportError:
    import pickle

# Import exaqute
from exaqute.ExaquteTaskPyCOMPSs import *   # to execute with pycompss
# from exaqute.ExaquteTaskHyperLoom import *  # to execute with the IT4 scheduler
# from exaqute.ExaquteTaskLocal import *      # to execute with python3
'''
get_value_from_remote is the equivalent of compss_wait_on: a synchronization point
in future, when everything is integrated with the it4i team, importing exaqute.ExaquteTaskHyperLoom you can launch your code with their scheduler instead of BSC
'''

'''Adapt the following class depending on the problem, deriving the MonteCarloAnalysis class from the problem of interest'''

'''
This Analysis Stage implementation solves the elliptic PDE in (0,1)^2 with zero Dirichlet boundary conditions
-lapl(u) = xi*f,    f= -432*x*(x-1)*y*(y-1)
                    f= -432*(x**2+y**2-x-y)
where xi is a Beta(2,6) random variable, and computes statistic of the QoI
Q = int_(0,1)^2 u(x,y)dxdy
'''
class MonteCarloAnalysis(AnalysisStage):
    '''Main analysis stage for Monte Carlo simulations'''
    def __init__(self,input_model,input_parameters,sample):
        self.sample = sample
        super(MonteCarloAnalysis,self).__init__(input_model,input_parameters)
        self._GetSolver().main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.NODAL_AREA)

    def _CreateSolver(self):
        import convection_diffusion_stationary_solver
        return convection_diffusion_stationary_solver.CreateSolver(self.model,self.project_parameters["solver_settings"])

    def _GetSimulationName(self):
        return "Monte Carlo Analysis"

    '''Introduce here the stochasticity in the right hand side defining the forcing function and apply the stochastic contribute'''
    def ModifyInitialProperties(self):
        for node in self.model.GetModelPart("MLMCLaplacianModelPart").Nodes:
            coord_x = node.X
            coord_y = node.Y
            # forcing = -432.0 * coord_x * (coord_x - 1) * coord_y * (coord_y - 1)
            forcing = -432.0 * (coord_x**2 + coord_y**2 - coord_x - coord_y) # this forcing presents the below commented analytical solution
            node.SetSolutionStepValue(KratosMultiphysics.HEAT_FLUX,forcing*self.sample)

##################################################
######## END OF CLASS MONTECARLOANALYSIS #########
##################################################

'''
function generating the random sample
here the sample has a beta distribution with parameters alpha = 2.0 and beta = 6.0
'''
def GenerateSample():
    alpha = 2.0
    beta = 6.0
    number_samples = 1
    sample = np.random.beta(alpha,beta,number_samples)
    return sample

'''
function evaluating the QoI of the problem: int_{domain} TEMPERATURE(x,y) dx dy
right now we are using the midpoint rule to evaluate the integral: improve!
'''
def EvaluateQuantityOfInterest(simulation):
    """here we evaluate the QoI of the problem: int_{domain} SOLUTION(x,y) dx dy
    we use the midpoint rule to evaluate the integral"""
    KratosMultiphysics.CalculateNodalAreaProcess(simulation._GetSolver().main_model_part,2).Execute()
    Q = 0.0
    for node in simulation._GetSolver().main_model_part.Nodes:
        Q = Q + (node.GetSolutionStepValue(KratosMultiphysics.NODAL_AREA)*node.GetSolutionStepValue(KratosMultiphysics.TEMPERATURE))
    return Q

'''
function called in the main returning a future object (the result class) and an integer (the finer level)
input:
        pickled_coarse_model      : pickled model
        pickled_coarse_parameters : pickled parameters
output:
        MonteCarloResults class   : class of the simulation results
        current_MC_level          : level of the current MLMC simulation
'''
def ExecuteMonteCarloAnalysis(pickled_model, pickled_parameters):
    current_MC_level = 0 # MC has only level 0
    return (ExecuteMonteCarloAnalysis_Task(pickled_model, pickled_parameters),current_MC_level)

'''
function executing the problem
input:
        model       : serialization of the model
        parameters  : serialization of the Project Parameters
output:
        QoI         : Quantity of Interest
'''
@ExaquteTask(returns=1)
def ExecuteMonteCarloAnalysis_Task(pickled_model, pickled_parameters):
    '''overwrite the old model serializer with the unpickled one'''
    model_serializer = pickle.loads(pickled_model)
    current_model = KratosMultiphysics.Model()
    model_serializer.Load("ModelSerialization",current_model)
    del(model_serializer)
    '''overwrite the old parameters serializer with the unpickled one'''
    serialized_parameters = pickle.loads(pickled_parameters)
    current_parameters = KratosMultiphysics.Parameters()
    serialized_parameters.Load("ParametersSerialization",current_parameters)
    del(serialized_parameters)
    '''initialize the MonteCarloResults class'''
    current_level = 0 # always 0 for MC
    mc_results_class = mc.MonteCarloResults(current_level)
    sample = GenerateSample()
    simulation = MonteCarloAnalysis(current_model,current_parameters,sample)
    simulation.Run()
    QoI = EvaluateQuantityOfInterest(simulation)
    mc_results_class.QoI[current_level].append(QoI) # saving results in the corresponding list, for MC only list of level 0
    return mc_results_class

'''
function serializing and pickling the model and the parameters of the problem
the idea is the following:
i)   from Model/Parameters Kratos object to StreamSerializer Kratos object
ii)  from StreamSerializer Kratos object to pickle string
iii) from pickle string to StreamSerializer Kratos object
iv)  from StreamSerializer Kratos object to Model/Parameters Kratos object
input:
        parameter_file_name   : path of the Project Parameters file
output:
        pickled_model      : model serializaton
        pickled_parameters : project parameters serialization
'''
@ExaquteTask(parameter_file_name=FILE_IN,returns=2)
def SerializeModelParameters_Task(parameter_file_name):
    with open(parameter_file_name,'r') as parameter_file:
        parameters = KratosMultiphysics.Parameters(parameter_file.read())
    local_parameters = parameters
    model = KratosMultiphysics.Model()
    # local_parameters["solver_settings"]["model_import_settings"]["input_filename"].SetString(model_part_file_name[:-5])
    fake_sample = GenerateSample()
    simulation = MonteCarloAnalysis(model,local_parameters,fake_sample)
    simulation.Initialize()
    serialized_model = KratosMultiphysics.StreamSerializer()
    serialized_model.Save("ModelSerialization",simulation.model)
    serialized_parameters = KratosMultiphysics.StreamSerializer()
    serialized_parameters.Save("ParametersSerialization",simulation.project_parameters)
    # pickle dataserialized_data
    pickled_model = pickle.dumps(serialized_model, 2) # second argument is the protocol and is NECESSARY (according to pybind11 docs)
    pickled_parameters = pickle.dumps(serialized_parameters, 2)
    print("\n","#"*50," SERIALIZATION COMPLETED ","#"*50,"\n")
    return pickled_model,pickled_parameters

if __name__ == '__main__':

    '''set the ProjectParameters.json path'''
    parameter_file_name = "../tests/PoissonSquareTest/parameters_poisson_finer.json"
    '''create a serialization of the model and of the project parameters'''
    pickled_model,pickled_parameters = SerializeModelParameters_Task(parameter_file_name)
    '''customize setting parameters of the ML simulation'''
    settings_MC_simulation = KratosMultiphysics.Parameters("""
    {
        "tolerance" : 0.1,
        "cphi" : 5e-1,
        "batch_size" : 20,
        "convergence_criteria" : "MC_higher_moments_sequential_stopping_rule"
    }
    """)
    '''contruct MonteCarlo class'''
    mc_class = mc.MonteCarlo(settings_MC_simulation)
    '''start MC algorithm'''
    while mc_class.convergence is not True:
        mc_class.InitializeMCPhase()
        mc_class.ScreeningInfoInitializeMCPhase()
        for instance in range (mc_class.difference_number_samples[0]):
            mc_class.AddResults(ExecuteMonteCarloAnalysis(pickled_model,pickled_parameters))
        mc_class.FinalizeMCPhase()
        mc_class.ScreeningInfoFinalizeMCPhase()

    mc_class.QoI.mean = get_value_from_remote(mc_class.QoI.mean)
    print("\nMC mean = ",mc_class.QoI.mean)

    ''' The below part evaluates the relative L2 error between the numerical solution SOLUTION(x,y,sample) and the analytical solution, also dependent on sample.
    Analytical solution available in case FORCING = sample * -432.0 * (coord_x**2 + coord_y**2 - coord_x - coord_y)'''
    # model_serializer = pickle.loads(pickled_model)
    # current_model = KratosMultiphysics.Model()
    # model_serializer.Load("ModelSerialization",current_model)
    # del(model_serializer)
    # serialized_parameters = pickle.loads(pickled_parameters)
    # current_parameters = KratosMultiphysics.Parameters()
    # serialized_parameters.Load("ParametersSerialization",current_parameters)
    # del(serialized_parameters)
    # sample = 1.0
    # simulation = MonteCarloAnalysis(current_model,current_parameters,sample)
    # simulation.Run()
    # KratosMultiphysics.CalculateNodalAreaProcess(simulation._GetSolver().main_model_part,2).Execute()
    # error = 0.0
    # L2norm_analyticalsolution = 0.0
    # for node in simulation._GetSolver().main_model_part.Nodes:
    #     local_error = ((node.GetSolutionStepValue(KratosMultiphysics.TEMPERATURE) - (432.0*simulation.sample*node.X*node.Y*(1-node.X)*(1-node.Y)*0.5))**2) * node.GetSolutionStepValue(KratosMultiphysics.NODAL_AREA)
    #     error = error + local_error
    #     local_analyticalsolution = (432.0*simulation.sample*node.X*node.Y*(1-node.X)*(1-node.Y)*0.5)**2 * node.GetSolutionStepValue(KratosMultiphysics.NODAL_AREA)
    #     L2norm_analyticalsolution = L2norm_analyticalsolution + local_analyticalsolution
    # error = np.sqrt(error)
    # L2norm_analyticalsolution = np.sqrt(L2norm_analyticalsolution)
    # print("L2 relative error = ", error/L2norm_analyticalsolution)


text:
import numpy as np

class Variance():
    @staticmethod
    def variance(data):
        return np.var(data)


text:
import os
from os.path import join as opj
import cv2
import numpy as np

BBOX_B = 10/15

class zone_tracker():
    def __init__(self):
        self.zone_list = list()
        self.id2id = dict()

    def update(self,det):
        have_up = False
        det_time = int(det['frame'][3:])
        det_x = 1280 - (det['bbox'][2]+det['bbox'][0])/2
        det_y = (det['bbox'][3]+det['bbox'][1])/2
        max_sim = 0
        max_zi = 0
        for zi,zl in enumerate(self.zone_list):
            zl_time = int(zl['frame'][3:])
            if det_time==zl_time:
                continue
            if det_time>zl_time+4:
                self.zone_list[zi] = -1
                continue
            if zl['id']==det['id']:
                self.zone_list[zi] = det
                have_up = True
                break
            if det['id'] in self.id2id:
                if self.id2id[det['id']]==zl['id']:
                    if det_time==zl_time:
                        have_up = True
                        break
                    else:
                        det['id'] = zl['id']
                        self.id2id[det['id']] = zl['id']
                        self.zone_list[zi] = det
                        have_up = True
                        break
            zl_x = 1280 - (zl['bbox'][2] + zl['bbox'][0]) / 2
            zl_y = (zl['bbox'][3] + zl['bbox'][1]) / 2
            if det_x<zl_x+10 and det_y<zl_y+10:
                continue
            if abs(det_x-zl_x)<15 and abs(det_y-zl_y)<15:
                # print("det:{}  zl:{}".format([det['id'],det_x,det_y],[zl['id'],zl_x,zl_y]))
                if det['id'] in self.id2id:
                    continue
                det['id'] = zl['id']
                self.id2id[det['id']]=zl['id']
                self.zone_list[zi] = det
                have_up = True
                break
        self.zone_list = [zl for zl in self.zone_list if zl!=-1]
        if not have_up:
            self.zone_list.append(det)
        return det

class zone():
    def __init__(self):
        # 0: b 1: g 3: r 123:w
        # w r 非高速
        # b g 高速
        zones = {}
        zone_path = "./zone"
        for img_name in os.listdir(zone_path):
            camnum = int(img_name.split('.')[0][-3:])
            zone_img = cv2.imread(opj(zone_path, img_name))
            zones[camnum] = zone_img

        self.zones = zones
        self.current_cam = 0

    def set_cam(self,cam):
        self.current_cam = cam

    def get_zone(self,bbox):
        cx = int((bbox[0] + bbox[2]) / 2)
        cy = int((bbox[1] + bbox[3]) / 2)
        pix = self.zones[self.current_cam][cy, cx, :]
        zone_num = 0
        if pix[0] > 50 and pix[1] > 50 and pix[2] > 50:  # w
            zone_num = 1
        if pix[0] < 50 and pix[1] < 50 and pix[2] > 50:  # r
            zone_num = 2
        if pix[0] < 50 and pix[1] > 50 and pix[2] < 50:  # g
            zone_num = 3
        if pix[0] > 50 and pix[1] < 50 and pix[2] < 50:  # b
            zone_num = 4
        return zone_num

    def is_ignore(self,zone_list,frame_list, cid):
        # 0 不在任何路口 1 白色 2 红色 3 绿色 4 蓝色
        zs, ze = zone_list[0], zone_list[-1]
        fs, fe = frame_list[0],frame_list[-1]
        if zs == ze:
            # 如果一直在一个区域里，排除
            if ze in [1,2]:
                return True
            if zs!=0 and 0 in zone_list:
                return False
            if fe-fs>1500:
                return True
            if fs<2:
                if cid in [45]:
                    return True
            if fe > 1999:
                if cid in [41]:
                    return True
            if fs<2 or fe>1999:
              if ze in [3,4]:
                return False
            return True
        else:
            # 如果区域发生变化
            if cid in [41, 42, 43, 44, 45, 46]:
                # 如果从支路进支路出，排除
                if zs == 1 and ze == 2:
                    return True
                if zs == 2 and ze == 1:
                    return True
            if cid in [41]:
                # 在41相机，车辆没有进出42相机
                if (zs in [1, 2]) and ze == 4:
                    return True
                if zs == 4 and (ze in [1, 2]):
                    return True
            if cid in [46]:
                # 在46相机，车辆没有进出45相机
                if (zs in [1, 2]) and ze == 3:
                    return True
                if zs == 3 and (ze in [1, 2]):
                    return True
            return False

    def filter_mot(self,mot_list, cid):
        new_mot_list = dict()
        for tracklet in mot_list:
            tracklet_dict = mot_list[tracklet]
            frame_list = list(tracklet_dict.keys())
            frame_list.sort()
            zone_list = []
            for f in frame_list:
                zone_list.append(tracklet_dict[f]['zone'])
            if not self.is_ignore(zone_list,frame_list, cid):
                new_mot_list[tracklet] = tracklet_dict
        return new_mot_list

    def filter_bbox(self,mot_list,cid):
        new_mot_list = dict()
        yh = self.zones[cid].shape[0]
        for tracklet in mot_list:
            # if tracklet==1812:
            #     print(tracklet)
            tracklet_dict = mot_list[tracklet]
            frame_list = list(tracklet_dict.keys())
            frame_list.sort()
            bbox_list = []
            for f in frame_list:
                bbox_list.append(tracklet_dict[f]['bbox'])
            bbox_x = [b[0] for b in bbox_list]
            bbox_y = [b[1] for b in bbox_list]
            bbox_w = [b[2]-b[0] for b in bbox_list]
            bbox_h = [b[3]-b[1] for b in bbox_list]
            new_frame_list = list()
            if 0 in bbox_x or 0 in bbox_y:
                b0 = [i for i, f in enumerate(frame_list) if bbox_x[i]<5 or bbox_y[i]+bbox_h[i]>yh-5]
                if len(b0)==len(frame_list):
                    if cid in [41,42,44,45,46]:
                        continue
                    max_w = max(bbox_w)
                    max_h = max(bbox_h)
                    for i,f in enumerate(frame_list):
                        if bbox_w[i] > max_w * BBOX_B and bbox_h[i] > max_h * BBOX_B:
                            new_frame_list.append(f)
                else:
                    l_i,r_i = 0,len(frame_list)-1
                    if b0[0]==0:
                        for i in range(len(b0)-1):
                            if b0[i]+1==b0[i+1]:
                                l_i=b0[i+1]
                            else:
                                break
                    if b0[-1]==len(frame_list)-1:
                        for i in range(len(b0) - 1):
                            i = len(b0)-1-i
                            if b0[i]-1==b0[i-1]:
                                r_i=b0[i-1]
                            else:
                                break

                    max_lw, max_lh = bbox_w[l_i], bbox_h[l_i]
                    max_rw, max_rh = bbox_w[r_i], bbox_h[r_i]
                    for i,f in enumerate(frame_list):
                        if i < l_i:
                            if bbox_w[i] > max_lw * BBOX_B and bbox_h[i] > max_lh * BBOX_B:
                                new_frame_list.append(f)
                        elif i > r_i:
                            if bbox_w[i] > max_rw * BBOX_B and bbox_h[i] > max_rh * BBOX_B:
                                new_frame_list.append(f)
                        else:
                            new_frame_list.append(f)
                new_tracklet_dict=dict()
                for f in new_frame_list:
                    new_tracklet_dict[f]=tracklet_dict[f]
                new_mot_list[tracklet] = new_tracklet_dict
            else:
                new_mot_list[tracklet] = tracklet_dict
        return new_mot_list

    def break_mot(self,mot_list,cid):
        # if not cid in [41,44,45,46]:
        #     return mot_list
        new_mot_list = dict()
        new_num_tracklets = max(mot_list)+1
        for tracklet in mot_list:
            tracklet_dict = mot_list[tracklet]
            frame_list = list(tracklet_dict.keys())
            frame_list.sort()
            zone_list = []
            back_tracklet = False
            new_zone_f = 0
            pre_frame = frame_list[0]
            time_break = False
            # if tracklet == 714:
            #     print(tracklet)
            for f in frame_list:
                if f-pre_frame>100:
                    if cid in [44,45]:
                        time_break = True
                        break
                if not cid in [41, 44, 45, 46]:
                    break
                pre_frame=f
                new_zone = tracklet_dict[f]['zone']
                if len(zone_list)>0 and zone_list[-1] == new_zone:
                    continue
                if new_zone_f>1:
                    if len(zone_list) > 1 and new_zone in zone_list:
                        back_tracklet = True
                    zone_list.append(new_zone)
                    new_zone_f=0
                else:
                    new_zone_f+=1
            if back_tracklet:
                new_tracklet_dict = dict()
                pre_bbox = -1
                pre_arrow = 0
                have_break = False
                for f in frame_list:
                    now_bbox = tracklet_dict[f]['bbox']
                    if pre_bbox == -1:
                        pre_bbox = now_bbox
                    now_arrow = now_bbox[0] - pre_bbox[0]
                    if pre_arrow*now_arrow<0 and len(new_tracklet_dict)>15 and not have_break:
                        new_mot_list[tracklet] = new_tracklet_dict
                        new_tracklet_dict = dict()
                        have_break = True
                    if have_break:
                        tracklet_dict[f]['id'] = new_num_tracklets
                    new_tracklet_dict[f] = tracklet_dict[f]
                    pre_bbox,pre_arrow = now_bbox,now_arrow
                if have_break:
                    new_mot_list[new_num_tracklets] = new_tracklet_dict
                    new_num_tracklets += 1
                else:
                    new_mot_list[tracklet] = new_tracklet_dict
            elif time_break:
                new_tracklet_dict = dict()
                have_break = False
                pre_frame = frame_list[0]
                for f in frame_list:
                    if f-pre_frame>100:
                        new_mot_list[tracklet] = new_tracklet_dict
                        new_tracklet_dict = dict()
                        have_break = True
                    new_tracklet_dict[f] = tracklet_dict[f]
                    pre_frame=f
                if have_break:
                    new_mot_list[new_num_tracklets] = new_tracklet_dict
                    new_num_tracklets += 1
                else:
                    new_mot_list[tracklet] = new_tracklet_dict
            else:
                new_mot_list[tracklet] = tracklet_dict
        print("old:{} new:{}".format(len(mot_list),len(new_mot_list)))
        return new_mot_list

    def comb_mot(self,mot_list,cid):
        time_mot = dict()
        for tid in mot_list:
            tracklet = mot_list[tid]
            for tf in tracklet:
                if tf in time_mot:
                    time_mot[tf].append(tracklet[tf])
                else:
                    time_mot[tf]=[tracklet[tf]]
        time_list = list(time_mot)
        time_list.sort()
        zone_tack = zone_tracker()
        new_mot_list = dict()
        for t in time_list:
            det_list = time_mot[t]
            for di,dt in enumerate(det_list):
                if dt['zone'] in [4]:
                    dt = zone_tack.update(dt)
                else:
                    if dt['id'] in zone_tack.id2id:
                        dt['id']=zone_tack.id2id[dt['id']]
                time_mot[t][di] = dt
                if dt['id'] in new_mot_list:
                    new_mot_list[dt['id']][int(dt['frame'][3:])]=dt
                else:
                    new_mot_list[dt['id']] = {int(dt['frame'][3:]): dt}
        print("old:{} new:{}".format(len(mot_list),len(new_mot_list)))
        return new_mot_list


text:
"""
The sequence execution part of

  Prefrontal neural correlates of memory for sequences.
  B. B. Averbeck & D. Lee, JNS 2007.

  http://www.jneurosci.org/content/27/9/2204.abstract

"""
from __future__ import division

import numpy as np

from pycog import tasktools

#-----------------------------------------------------------------------------------------
# Sequence-related definitions
#-----------------------------------------------------------------------------------------

# Sequence definitions
sequences = {
    1: [4] + [5, 1, 0],
    2: [4] + [3, 1, 2],
    3: [4] + [3, 7, 8],
    4: [4] + [5, 7, 6],
    5: [4] + [5, 1, 2],
    6: [4] + [3, 1, 0],
    7: [4] + [3, 7, 6],
    8: [4] + [5, 7, 8]
}
nseq = len(sequences)

# Possible targets from each position
#
#   0 1 2
#   3 4 5
#   6 7 8
#
options = {
    1: [0, 2],
    3: [1, 7],
    4: [3, 5],
    5: [1, 7],
    7: [6, 8]
}

#-----------------------------------------------------------------------------------------
# Network structure
#-----------------------------------------------------------------------------------------

Nin  = 9 + nseq
N    = 200
Nout = 2

# For addressing inputs
DOTS     = range(9)
SEQUENCE = range(9, 9+nseq)

# E/I
ei, EXC, INH = tasktools.generate_ei(N)

# Time constant
tau = 50

#-----------------------------------------------------------------------------------------
# Noise
#-----------------------------------------------------------------------------------------

var_rec = 0.01**2

#-----------------------------------------------------------------------------------------
# Task structure
#-----------------------------------------------------------------------------------------

# Screen coordinates
x0, y0 = -0.5, -0.5
dx, dy = +0.5, +0.5

# Convert dots to screen coordinates
def target_position(k):
    j = 2 - k//3
    i = k % 3

    return x0+i*dx, y0+j*dy

def generate_trial(rng, dt, params):
    #---------------------------------------------------------------------------------
    # Select task condition
    #---------------------------------------------------------------------------------

    if params['name'] in ['gradient', 'test']:
        seq = params.get('seq', rng.choice(sequences.keys()))
    elif params['name'] == 'validation':
        b = params['minibatch_index'] % nseq
        if b == 0:
            generate_trial.seqs = rng.permutation(nseq)

        seq = generate_trial.seqs[b] + 1
    else:
        raise ValueError("Unknown trial type.")

    #-------------------------------------------------------------------------------------
    # Epochs
    #-------------------------------------------------------------------------------------

    iti      = 1000
    fixation = 1000
    M1       = 500
    M2       = 500
    M3       = 500
    T        = iti + fixation + M1 + M2 + M3

    epochs = {
        'iti':      (0, iti),
        'fixation': (iti, iti + fixation),
        'M1':       (iti + fixation, iti + fixation + M1),
        'M2':       (iti + fixation + M1, iti + fixation + M1 + M2),
        'M3':       (iti + fixation + M1 + M2, iti + fixation + M1 + M2 + M3),
        }
    epochs['T'] = T

    #---------------------------------------------------------------------------------
    # Trial info
    #---------------------------------------------------------------------------------

    t, e  = tasktools.get_epochs_idx(dt, epochs) # Time, task epochs in discrete time
    trial = {'t': t, 'epochs': epochs}           # Trial

    trial['info'] = {'seq': seq}

    #---------------------------------------------------------------------------------
    # Inputs
    #---------------------------------------------------------------------------------

    # Input matrix
    X = np.zeros((len(t), Nin))

    # Which sequence?
    X[:,SEQUENCE[seq-1]] = 1

    # Sequence
    sequence = sequences[seq]

    # Options
    X[e['fixation'],sequence[0]] = 1
    for I, J in zip([e['M1'], e['M2'], e['M3']],
                    [[sequence[0]] + options[sequence[0]],
                     [sequence[1]] + options[sequence[1]],
                     [sequence[2]] + options[sequence[2]]]):
        for j in J:
            X[I,j] = 1

    # Inputs
    trial['inputs'] = X

    #---------------------------------------------------------------------------------
    # Target output
    #---------------------------------------------------------------------------------

    if params.get('target_output', False):
        Y = np.zeros((len(t), Nout)) # Output matrix
        M = np.zeros((len(t), Nout)) # Mask matrix

        # Hold gaze
        Y[e['fixation'],:] = target_position(sequence[0])
        Y[e['M1'],:]       = target_position(sequence[1])
        Y[e['M2'],:]       = target_position(sequence[2])
        Y[e['M3'],:]       = target_position(sequence[3])

        # We don't constrain the intertrial interval
        M[e['fixation']+e['M1']+e['M2']+e['M3'],:] = 1

        # Output and mask
        trial['outputs'] = Y
        trial['mask']    = M

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

    return trial

min_error = 0.05

mode         = 'continuous'
n_validation = 100*nseq

max_stars_repo_name: oesst/pycog

max_issues_repo_name: oesst/pycog

max_forks_repo_name: oesst/pycog


text:
import numpy as np

from berenet.base_layer import BaseLayer
from copy import deepcopy

class RecurrentLayer(BaseLayer):
	def __init__(self, inputs, outputs, recurrent=True, is_output=False, bias=1, function=BaseLayer.FUNCTIONS[1]):
		assert function in super().FUNCTIONS, 'Invalid function, please make sure it\'s one of: ' + str(super().FUNCTIONS)

		self.Fp = list()
		self.Y = list()

		self.D = list()

		self.inputs = inputs
		self.outputs = outputs
		self.is_output = is_output
		self.recurrent = recurrent
		self.function = function
		self.bias = bias

		if not is_output:
			if bias:
				self.W = np.random.normal(size=(inputs + bias, outputs), scale=1E-4)
			else:
				self.W = np.random.normal(size=(inputs, outputs), scale=1E-4)
		else:
			self.bias = False

		if recurrent:
			self.RW = np.random.normal(size=(1, inputs), scale=1E-4)
			self.recurrentD = list()

	def forward_pass(self, X):
		if len(self.Y) > 0 and self.recurrent:
			if self.bias:
				X += self.RW * self.Y[-1][:,:-self.bias]
			else:
				X += self.RW * self.Y[-1]

		activation = self._activate(X)
		activated_output = activation[0]

		if self.bias:
			activated_output = np.ones(shape=(activation[0].shape[0], activation[0].shape[1] + self.bias))
			activated_output[:,:-self.bias] = activation[0]

		self.Y.append(activated_output)
		self.Fp.append(activation[1])

		if self.is_output:
			return activation[0]

		return np.dot(activated_output, self.W)

	def reset(self):
		del self.Fp[:]
		del self.Y[:]

		self.clear_deltas()

	def clear_deltas(self):
		del self.D[:]

		if self.recurrent:
			del self.recurrentD[:]

	@staticmethod
	def softmax(X):
		for i in range(X.shape[0]):
			e_x = np.exp(X[i] - np.max(X[i]))
			X[i] = e_x / e_x.sum()


text:
"""
Simple thresholding applied to a real image using np.arange() to create the different threshold values
"""

# Import required packages:
import numpy as np
import cv2
from matplotlib import pyplot as plt

def show_img_with_matplotlib(color_img, title, pos):
    """Shows an image using matplotlib capabilities"""

    # Convert BGR image to RGB
    img_RGB = color_img[:, :, ::-1]

    ax = plt.subplot(3, 3, pos)
    plt.imshow(img_RGB)
    plt.title(title)
    plt.axis('off')

# Create the dimensions of the figure and set title and color:
fig = plt.figure(figsize=(9, 9))
plt.suptitle("Thresholding using np.arange() to create the different threshold values", fontsize=14, fontweight='bold')
fig.patch.set_facecolor('silver')

# Load the image and convert it to grayscale:
image = cv2.imread('sudoku.png')
gray_image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

# Plot the grayscale images and the histograms:
show_img_with_matplotlib(cv2.cvtColor(gray_image, cv2.COLOR_GRAY2BGR), "img", 1)

# Get the array with the values for thresholding in the range [60-130] with step 10:
# This function returns an array of evenly spaced values:
threshold_values = np.arange(start=60, stop=140, step=10)
# print(threshold_values)

# Apply cv2.threshold() with the different threshold values defined in 'threshold_values'
thresholded_images = []
for threshold in threshold_values:
    ret, thresh = cv2.threshold(gray_image, threshold, 255, cv2.THRESH_BINARY)
    thresholded_images.append(thresh)

# Show the thresholded images:
for index, (thresholded_image, threshold_value) in enumerate(zip(thresholded_images, threshold_values)):
    show_img_with_matplotlib(cv2.cvtColor(thresholded_image, cv2.COLOR_GRAY2BGR), "threshold = " + str(threshold_value),
                             index + 2)

# Show the Figure:
plt.show()


text:
"""
Stack diffuse X-ray emission in galaxy groups from ROSAT All-Sky Survey Data.
Author: Zack Hutchens
"""
from scipy import ndimage
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import colors, cm
from astroquery.skyview import SkyView as sv
from astropy.io import fits
from astropy.stats import sigma_clipped_stats
from photutils import DAOStarFinder, CircularAperture
from scipy.ndimage import gaussian_filter
import os
import sys
from skyview_downloader import download_images_java
import pickle

def scale_image(output_coords,scale):
    return (output_coords[0]/scale+150-150/scale, output_coords[1]/scale+150-150/scale)

def get_circle(R):
    theta=np.linspace(0,2*np.pi,1000)
    x = R*np.cos(theta)
    y = R*np.sin(theta)
    return x,y

class rosat_xray_stacker:
    """
    A class for stacking ROSAT All-Sky Survey (RASS) X-ray images of galaxy
    groups, including for image retreival (from NASA SkyView), image processing,
    and stacking.

    Attributes
    -----------------
    grpid : iterable
        List of group ID numbers. Must be unique, no duplicates.
        (Must pass by group and not by galaxy.)
    grpra : iterable
        Right-ascensions of galaxy groups in decimal degrees.
        This coordinate determines the position for image extraction/stacking.
    grpdec : iterable
        Declinations of galaxy groups in decimal degrees.
        This coordinate determines the position for image extraction/stacking.
    grpcz : iterable
        Recessional velocities of groups in km/s.
    centralname : str or iterable, default None
        Name of central galaxy in each group, or other group identifier
        beyond grpid.
    surveys : iterable, default None
        Surveys from which to extract/stack images, corresponding to the
        documentation from NASA SkyView (e.g. 'RASS-Int Hard').
    """
    def __init__(self, grpid, grpra, grpdec, grpcz, centralname='', surveys=None):
        self.grpid = np.array(grpid)
        self.grpra = np.array(grpra)
        self.grpdec = np.array(grpdec)
        self.grpcz = np.array(grpcz)
        if isinstance(surveys, str):
            surveys = [surveys]
        if not isinstance(centralname, str):
            centralname=list(centralname)
        self.surveys = surveys
        self.centralname = centralname
        self.goodflag = np.full(len(self.grpid),1)
        self.detection = np.zeros(len(self.grpid))

    def __repr__(self):
        return "rosat_xray_stacker instance at "+str(hex(id(self)))+" (N = {a} groups)".format(a=len(self.grpid))

    def __str__(self):
        return self.__repr__()

    def to_pickle(self, savename):
        """
        Save the rosat_xray_stacker instance to a serialized
        Pickle file, for later access. e.g. preparing images
        in one code, and stacking in another.

        Parameters
        ----------------
        savename : str
            Location on disk where serial file should be saved.
        """
        filep = open(savename, 'wb')
        pickle.dump(self, filep, protocol=pickle.HIGHEST_PROTOCOL)
        filep.close()

    def download_images(self, imgfilepath):
        """
        Download RASS images corresponding to the given group catalog.
        The number of download images is N*M, where N is the number of
        galaxy groups (i.e. len(grpra)) and M is len(surveys). If the
        provided path already contains <= N*M FITS files, the program
        will ask permission to re-download the images. This will over-
        write any existing data in imgfilepath.

        Parameters
        ------------------
        imgfilepath : str
            System file path where download images should be stored.

        Returns
        -------------------
        None
        """
        dirfiles = os.listdir(imgfilepath)
        nimagesneeded = len(self.surveys)*len(self.grpra)
        nimagesindir = np.sum([int('.fits' in fname) for fname in dirfiles])
        if nimagesindir==0:
            print('Provided directory is empty... downloading RASS images.')
            download_images_java(imgfilepath, self.grpra, self.grpdec, self.grpid, self.surveys, self.centralname)
        elif nimages>0 and nimagesindir<=nimagesneeded:
            print('The provided directory already contains some FITS images.')
            answ=input("Enter 'c' to re-download RASS images into the provided directory, or enter any other key to exit the program: ")
            if answ=='c':download_images_java(imgfilepath, self.grpra, self.grpdec, self.grpid, self.surveys, self.centralname)
            else:sys.exit()
        else:
            print('Error in download_images: please clear the directory and try re-downloading.')
            sys.exit()

    def sort_images(self, imagefiledir, snrmin=0.2, maxzero=0.98):
        """
        Sort raw RASS images into good and poor coverage based on a SNR
        and zero-pixel threshold. This function creates a new attribute
        of the stacker object called `rosat_xray_stacker.goodimage`,
        a boolean array that indicates whether the raw RASS image had good
        or poor coverage.

        Parameters
        ----------------------------
        rawpath : str
            Path location where raw images are stored.
        snrmin : scalar
            Minimum signal-to-noise ratio used to make cut. Default
            is 0.2 (anything below 0.2 marked 'poor'.
        maxzero : scalar
            Maximum fraction of zero pixels for a good image, default 0.98.
        outpath : str, default None
            Path where good images should be written. Default none.

        Returns
        ----------------------------
        None. The good image flag attribute is modified in place.

        """
        # sort directory files with order of group catalog
        imagefiles = os.listdir(imagefiledir)
        imagenames = np.array(os.listdir(imagefiledir))
        assert len(self.grpid)==len(imagenames), "Number of files in directory must match number of groups."
        imageIDs = np.array([float(imgnm.split('_')[2][3:]) for imgnm in imagenames])
        _, order = np.where(self.grpid[:,None]==imageIDs)
        imageIDs = imageIDs[order]
        imagenames = imagenames[order]
        assert (imageIDs==self.grpid).all(), "ID numbers are not sorted properly."

        # make SNR and zero-pixel assessments of images
        snr = np.zeros(len(imagefiles))
        zcount = np.zeros(len(imagefiles))
        goodflag = np.zeros(len(imagefiles))
        for i,imgpath in enumerate(imagenames):
            if imgpath.endswith('.fits'):
                hdulist = fits.open(imagefiledir+imgpath,memap=False)
                image = hdulist[0].data.flatten()
                hdulist.close()
                snr[i] = np.mean(image)/np.std(image)
                zcount[i] = len(image[image==0])/len(image)
        sel = np.where(np.logical_and(snr>=snrmin, zcount<maxzero))
        goodflag[sel]=1
        self.goodflag = goodflag

    @staticmethod
    def calculate_rms(nparray):
        return np.sqrt(np.mean(nparray**2, axis=0))

    def measure_SNR_1Mpc(self,imagefiledir,snrthreshold=3):
        """
        Measure the signal-to-noise ratio within the central 1Mpc sky area of the
        group image. This SNR measurement is used to determine whether the good
        image is marked as a detection or nondetection, and thus whether point
        source removal should be applied.

        Parameters
        -------------------
        imagepath : str
            Directory where raw images were stored. Only images with goodflag=1
            will have their SNR measured.
        snrthreshold : float, default 3
            Confidence threshold for which signals are considered detections. Default 3-sigma.

        Returns
        -------------------
        None. Modifies rosat_xray_stacker.detection in place and creates rosat_xray_stacker.centralSNR.
        """
        imagefiles = os.listdir(imagefiledir)
        imagenames = np.array(os.listdir(imagefiledir))
        assert len(self.grpid)==len(imagenames), "Number of files in directory must match number of groups."
        imageIDs = np.array([float(imgnm.split('_')[2][3:]) for imgnm in imagenames])
        _, order = np.where(self.grpid[:,None]==imageIDs)
        imageIDs = imageIDs[order]
        imagenames = imagenames[order]
        assert (imageIDs==self.grpid).all(), "ID numbers are not sorted properly."

        # make SNR measurement
        snr = np.zeros(len(imagenames))
        X,Y = np.meshgrid(np.arange(0,300,1), np.arange(0,300,1))
        for i,imgpath in enumerate(imagenames):
            if imgpath.endswith('.fits'):
                hdulist = fits.open(imagefiledir+imgpath,memap=False)
                image = hdulist[0].data
                hdulist.close()
                radius = (1/(self.grpcz[i]/70.))*206265/45. # in px
                dist_from_center = np.sqrt((X-150.)**2. + (Y-150.)**2.)
                measuresel = np.where(np.logical_and(dist_from_center<radius, image>0))
                #snr[i] = np.mean(image[measuresel])/np.std(image[measuresel])
                snr[i] = np.mean(image[measuresel])/np.std(image)
                #snr[i] = np.sum(image[measuresel])/np.sqrt(np.mean(image[image>0]**2.))
                #snr[i] = np.sum(image[measuresel])/self.calculate_rms(image[image>0].flatten())
                print(np.sum(image[measuresel])/np.sqrt(np.sum(image[measuresel])))
        self.centralSNR = snr
        self.detection = (snr>snrthreshold) & (self.goodflag==1)

    def mask_point_sources(self, imgfiledir, outfiledir, scs_cenfunc=np.mean, scs_sigma=3, scs_maxiters=2, smoothsigma=1.0,\
                        starfinder_fwhm=3, starfinder_threshold=8, mask_aperture_radius=5, imagewidth=300,\
                        imageheight=300, examine_result=False):
        """
        Mask point sources in a series of X-ray FITS images.
        This code will read the raw images, find sources and
        apply masks, and write-out masked images.

        Parameters
        ------------------
        imgfiledir : str
            File path where raw images are stored (use a trailing /; e.g. "/home/images/").
        outfiledir : str
            Directory where masked images should be written (use trailing /).
        scs_cenfunc : callable, default np.mean
            Function to set the center value for sigma-clipping statistics, use np.mean or np.median.
        scs_sigma : int, default 3
            Number of standard deviations to use as the upper/lower clipping limit.
        scs_maxiters : int, default 2
            Maximum number of sigma-clipping iterations to perform.
        smoothsigma : float, default 1.0
            Standard deviation of Gaussian kernel used for smoothing, prior to point source detection.
            If this parameter is None, then no smoothing is performed.
        starfinder_fwhm : int, default 3
            The full width-half maximum (FWHM) of the major axis of Gaussian kernel (in pixels) used by DAOStarFinder.
        starfinder_threshold : int, default 8
            From DAOStarFinder, "the absolute image value above which to select sources".
        mask_aperture_radius : int, default 5
            Radius of applied masks (pixels).
        imagewidth : int, default 300
            Width of masked image for DAOStarFinder; should match size of input image.
        imageheight : int, default 300
            Height of masked image for DAOStarFinder; should match size of input image.
        examine_result : bool, default False
            If True, display the raw image, masks, and masked image in each iteration.

        Returns
        ------------------
        None. All output images (with point sources masked) are written to `outfilepath`.

        """
        assert callable(scs_cenfunc), "Argument `cenfunc` for sigma_clipped_stats must be callable."
        imagefiles = os.listdir(imgfiledir)
        for grp in self.grpid:
            pass
        for i,imgpath in enumerate(imagefiles): ### iterate by group and not by image!!
            # get image
            hdulist = fits.open(imgfiledir+imgpath)
            image = hdulist[0].data

            # get image stats
            mean, median, std = sigma_clipped_stats(image[image!=0],sigma=scs_sigma, maxiters=scs_maxiters, cenfunc=scs_cenfunc)
            mean2, median2, std2 = sigma_clipped_stats(image,sigma=scs_sigma, maxiters=np.max([1,scs_maxiters-1]), cenfunc=scs_cenfunc)

            # smooth image before finding sources
            if smoothsigma is not None:
                smoothimg = gaussian_filter(image, sigma=smoothsigma)
            else:
                smoothimg = np.copy(image)

            # find point sources using DAOStarFinder (photutils)
            daofind = DAOStarFinder(fwhm=starfinder_fwhm, threshold=mean+starfinder_threshold*std)
            table = daofind.find_stars(smoothimg)

            if table is not None:
                # create and apply masks (unless it is a diffuse bright source?)
                positions=np.transpose(np.array([table['xcentroid'],table['ycentroid']]))
                apertures=CircularAperture(positions,r=mask_aperture_radius)
                masks=apertures.to_mask(method='center')

                # Create new image
                newimage = np.zeros_like(image)
                newmask = np.zeros_like(image)
                newmask = np.sum(np.array([msk.to_image(shape=((imagewidth,imageheight))) for msk in masks]), axis=0)

                replacesel = np.logical_and(newmask>0,image>mean+std)
                newimage[replacesel] = mean2
                newimage[~replacesel] = image[~replacesel]

                # examine, if requested
                if examine_result:
                    print(mean2)
                    fig, ax = plt.subplots(ncols=4, figsize=(16,7))
                    ax[0].imshow(image, vmin=0, vmax=np.max(image))
                    #ax[0].plot(table['xcentroid'], table['ycentroid'], 'x', color='yellow')
                    ax[0].set_title("Raw Image")
                    ax[1].imshow(smoothimg, vmin=0, vmax=np.max(image))
                    ax[1].set_title("Gaussian-smoothed Image\n Smoothing Sigma = "+str(smoothsigma))
                    ax[2].imshow(newmask, cmap='binary')
                    ax[2].set_title("Masks")
                    ax[3].imshow(newimage, vmin=0, vmax=np.max(image))
                    ax[3].set_title("Masked Image")
                    plt.show()
            else:
                print('skipping '+imgpath+': no point sources found')
                newimage=np.copy(image)
            # write to file and continue
            hdulist[0].data=newimage
            savepath=outfiledir+imgpath#[:-5]+"_pntsourcesremoved.fits"
            hdulist.writeto(savepath)
            hdulist.close()

    def scale_subtract_images(self, imagefiledir, outfiledir, progressConf=False):
        """
        Subtract >5*sigma pixels from images and scale images to
        a common redshift.

        Parameters
        -------------
        imgfiledir : str
            Path to directory containing input images for scaling.
            Each FITS file in this directory must be named consistently
            with the rest of this program (e.g. RASS-Int_Broad_grp13_ECO03822.fits).
        outfiledir : str
            Path where scaled images should be written.

        Returns
        -------------
        Scaled/subtracted images are written to the specified path.
        """
        czmax = np.max(self.grpcz)
        imagenames = np.array(os.listdir(imagefiledir))
        imageIDs = np.array([float(imgnm.split('_')[2][3:]) for imgnm in imagenames])
        for k in range(0,len(imagenames)):
            hdulist = fits.open(imagefiledir+imagenames[k], memap=False)
            img = hdulist[0].data
            im2=img*1
            mean=np.mean(im2[np.where(im2!=0)])
            std = np.std(im2[np.where(im2!=0)])
            im2[im2 > mean+5*std]=mean
            im2[130:170,130:170]=img[130:170,130:170] # preserve inner portion
            img = np.copy(im2)
            czsf = self.grpcz[self.grpid==imageIDs[k]]/czmax
            img = ndimage.geometric_transform(img, scale_image, cval=0, extra_keywords={'scale':czsf})
            hdulist[0].data = img
            hdulist.writeto(outfiledir+imagenames[k], overwrite=True)
            hdulist.close()
            if progressConf: print("Image {} complete.".format(k))

    def stack_images(self, imagefiledir, stackproperty, binedges):
        return self.stack_images_func(self.grpid, self.grpcz, imagefiledir, stackproperty, binedges)

    @staticmethod
    def stack_images_func(grpid, grpcz, imagefiledir, stackproperty, binedges):
        """
        Stack X-ray images of galaxy groups in bins of group properties
        (e.g. richness or halo mass).

        Parameters
        --------------------
        imgfiledir : str
            Path to directory containing input images for stacking.
            Each FITS file in this directory must be named consistently
            with the rest of this program (e.g. RASS-Int_Broad_grp13_ECO03822.fits).
        stackproperty : iterable
            Group property to be used for binning (e.g. halo mass). This
            list should include an entry for *every* group, as to match
            the length of self.grpid.
        bins : iterable
            Array of bins for stacking. It should represent the bin *edges*.
            Example: if bins=[11,12,13,14,15,16], then the resulting bins
            are [11,12], [12,13], [13,14], [14,15], [15,16].

        Returns
        --------------------
        groupstackID : np.array
            The ID of stack to which each galaxy group in self.grpid belongs.
        n_in_bin : list
            Number of images contributing to each stack.
        bincenters : np.array
            Center of each bin used in stacking.
        finalimagelist : list
            List of final stacked images, i.e. finalimagelist[i] is a sigma-clipped
            average of all images whose groups satisifed the boudnary conditions of
            bin[i]. Each image is a 2D numpy array.
        """
        imagenames = np.array(os.listdir(imagefiledir))
        assert len(grpid)==len(imagenames), "Number of files in directory must match number of groups."
        imageIDs = np.array([float(imgnm.split('_')[2][3:]) for imgnm in imagenames])
        _, order = np.where(grpid[:,None]==imageIDs)
        imageIDs = imageIDs[order]
        imagenames = imagenames[order]
        assert (imageIDs==grpid).all(), "ID numbers are not sorted properly."

        stackproperty = np.asarray(stackproperty)
        groupstackID = np.zeros_like(stackproperty)
        binedges = np.array(binedges)
        leftedges = binedges[:-1]
        rightedges = binedges[1:]
        bincenters = (leftedges+rightedges)/2.
        finalimagelist = []
        n_in_bin=[]
        for i in range(0,len(bincenters)):
            stacksel = np.where(np.logical_and(stackproperty>=leftedges[i], stackproperty<rightedges[i]))
            imagenamesneeded = imagenames[stacksel]
            imageIDsneeded = imageIDs[stacksel]
            groupstackID[stacksel]=i+1
            images_to_stack = []
            for j in range(0,len(imagenamesneeded)):
                img = imagenamesneeded[j]
                hdulist = fits.open(imagefiledir+img, memmap=False)
                img = hdulist[0].data
                hdulist.close()
                images_to_stack.append(img)
            avg, median, std = sigma_clipped_stats(np.array(images_to_stack), sigma=10., maxiters=1, axis=0)
            n_in_bin.append(len(images_to_stack))
            finalimagelist.append(avg)
            print("Bin {} done.".format(i))
        return groupstackID, n_in_bin, bincenters, finalimagelist

if __name__=='__main__':
    ecocsv = pd.read_csv("../g3groups/ECOdata_G3catalog_luminosity.csv")
    ecocsv = ecocsv[ecocsv.g3fc_l==1] # centrals only
    eco = rosat_xray_stacker(ecocsv.g3grp_l, ecocsv.g3grpradeg_l, ecocsv.g3grpdedeg_l, ecocsv.g3grpcz_l, centralname=ecocsv.name, surveys=['RASS-Int Hard'])
    eco.sort_images('./g3rassimages/eco/')
    eco.measure_SNR_1Mpc('./g3rassimages/eco/')
    #eco.download_images('./g3rassimages/eco/')
    #eco.mask_point_sources('/srv/two/zhutchen/rosat_xray_stacker/g3rassimages/eco/', 'whatever/', examine_result=True, smoothsigma=3, starfinder_threshold=5)
    #eco.mask_point_sources('/srv/scratch/zhutchen/eco03822files/', 'whatever/', examine_result=True, smoothsigma=3, starfinder_threshold=5)
    #eco.scale_subtract_images("./g3rassimages/eco/", "./g3rassimages/eco_scaled/", True)
    nbin, bincenters, images = eco.stack_images("./g3rassimages/eco_scaled/", "whatever", np.asarray(ecocsv.g3logmh_l), binedges=np.arange(12,16,1))
    Rvirs = ((3*10**bincenters) / (4*np.pi*337*0.3*1.36e11) )**(1/3)
    rvirscales = 0.5 * Rvirs/(7000/70.) * 206265 / 45. / 3.
    print(rvirscales)
    images = [gaussian_filter(images[i],2) for i in range(0,len(images))]
    maxes = np.asarray([np.max(im) for im in images])
    scaleto = np.mean(maxes)-0.5*np.std(maxes)
    czmax = np.max(ecocsv.g3grpcz_l)
    print(czmax)
    for index, image in enumerate(images):
        plt.figure()
        plt.imshow(image, extent=[-150,150,-150,150],vmax=scaleto,vmin=0)
        Rvir = ((3*10**bincenters[index]) / (4*np.pi*337*0.3*1.36e11) )**(1/3)
        print(Rvir)
        Rvirx, Rviry = get_circle(Rvir/(czmax/70.) * 206265 * (1/45))
        plt.plot(Rvirx, Rviry, color='orange', linewidth=2)
        plt.title(r'$<\log M_{\rm vir}>=$ '+'{:0.2f}'.format(bincenters[index])+' (N={})'.format(nbin[index]))
        plt.show()

reason: import numpy,from scipy,from astropy