text:
import logging
import os
import random
import time

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import tensorflow as tf
import helper as util

matplotlib.use('agg')

INPUT_IMAGE_DIR = "input"
INTERPOLATED_IMAGE_DIR = "interpolated"
TRUE_IMAGE_DIR = "true"

class SuperResolution:
    def __init__(self, flags, model_name=""):

        # Model Parameters
        self.filters = flags.filters
        self.min_filters = flags.min_filters
        self.nin_filters = flags.nin_filters
        self.nin_filters2 = flags.nin_filters2 if flags.nin_filters2 != 0 else flags.nin_filters // 2
        self.cnn_size = flags.cnn_size
        self.last_cnn_size = flags.last_cnn_size
        self.cnn_stride = 1
        self.layers = flags.layers
        self.nin = flags.nin
        self.bicubic_init = flags.bicubic_init
        self.dropout = flags.dropout
        self.activator = flags.activator
        self.filters_decay_gamma = flags.filters_decay_gamma

        # Training Parameters
        self.initializer = flags.initializer
        self.weight_dev = flags.weight_dev
        self.l2_decay = flags.l2_decay
        self.optimizer = flags.optimizer
        self.beta1 = flags.beta1
        self.beta2 = flags.beta2
        self.momentum = flags.momentum
        self.batch_num = flags.batch_num
        self.batch_image_size = flags.batch_image_size
        if flags.stride_size == 0:
            self.stride_size = flags.batch_image_size // 2
        else:
            self.stride_size = flags.stride_size

        # Learning Rate Control for Training
        self.initial_lr = flags.initial_lr
        self.lr_decay = flags.lr_decay
        self.lr_decay_epoch = flags.lr_decay_epoch

        # Dataset or Others
        self.dataset = flags.dataset
        self.test_dataset = flags.test_dataset

        # Image Processing Parameters
        self.scale = flags.scale
        self.max_value = flags.max_value
        self.channels = flags.channels
        self.jpeg_mode = flags.jpeg_mode
        self.output_channels = self.scale * self.scale

        # Environment (all directory name should not contain '/' after )
        self.checkpoint_dir = flags.checkpoint_dir
        self.tf_log_dir = flags.tf_log_dir

        # Debugging or Logging
        self.debug = flags.debug
        self.save_loss = flags.save_loss
        self.save_weights = flags.save_weights
        self.save_images = flags.save_images
        self.save_images_num = flags.save_images_num
        self.log_weight_image_num = 32

        # initialize variables
        self.name = self.get_model_name(model_name)
        self.batch_input = self.batch_num * [None]
        self.batch_input_quad = np.zeros(
            shape=[self.batch_num, self.batch_image_size, self.batch_image_size, self.scale * self.scale])
        self.batch_true_quad = np.zeros(
            shape=[self.batch_num, self.batch_image_size, self.batch_image_size, self.scale * self.scale])
        self.receptive_fields = 2 * self.layers + self.cnn_size - 2
        self.complexity = 0

        # initialize environment
        util.make_dir(self.checkpoint_dir)
        util.make_dir(flags.graph_dir)
        util.make_dir(self.tf_log_dir)
        if flags.initialise_tf_log:
            util.clean_dir(self.tf_log_dir)
        util.set_logging(flags.log_filename, stream_log_level=logging.INFO, file_log_level=logging.INFO,
                         tf_log_level=tf.logging.WARN)
        config = tf.ConfigProto()
        config.gpu_options.allow_growth = True
        self.sess = tf.InteractiveSession(config=config)
        self.init_train_step()

        logging.info("\nDCSCN -------------------------------------")
        logging.info("%s [%s]" % (util.get_now_date(), self.name))

    def get_model_name(self, model_name):
        if model_name is "":
            name = "dcscn_L%d_F%d" % (self.layers, self.filters)
            if self.min_filters != 0:
                name += "to%d" % self.min_filters
            if self.filters_decay_gamma != 1.0:
                name += "_G%2.2f" % self.filters_decay_gamma
            if self.cnn_size != 3:
                name += "_C%d" % self.cnn_size
            if self.scale != 2:
                name += "_Sc%d" % self.scale
            if self.nin:
                name += "_NIN"
            if self.nin_filters != 0:
                name += "_A%d" % self.nin_filters
                if self.nin_filters2 != self.nin_filters // 2:
                    name += "_B%d" % self.nin_filters2
            if self.bicubic_init:
                name += "_BI"
            if self.dropout != 1.0:
                name += "_D%0.2f" % self.dropout
            if self.max_value != 255.0:
                name += "_M%2.1f" % self.max_value
            if self.activator != "relu":
                name += "_%s" % self.activator
            if self.dataset != "yang91":
                name += "_" + self.dataset
            if self.batch_image_size != 32:
                name += "_B%d" % self.batch_image_size
            if self.last_cnn_size != 1:
                name += "_L%d" % self.last_cnn_size
        else:
            name = "dcscn_%s" % model_name

        return name

    def load_datasets(self, target, data_dir, batch_dir, batch_image_size, stride_size=0):

        print("Loading datasets for [%s]..." % target)
        util.make_dir(batch_dir)

        if stride_size == 0:
            stride_size = batch_image_size // 2

        if self.bicubic_init:
            resampling_method = "bicubic"
        else:
            resampling_method = "nearest"

        datasets = util.DataSets(self.scale, batch_image_size, stride_size, channels=self.channels,
                                   jpeg_mode=self.jpeg_mode, max_value=self.max_value,
                                   resampling_method=resampling_method)

        if not datasets.is_batch_exist(batch_dir):
            datasets.build_batch(data_dir, batch_dir)

        if target == "training":
            datasets.load_batch_train(batch_dir)
            self.train = datasets
        else:
            datasets.load_batch_test(batch_dir)
            self.test = datasets

    def init_epoch_index(self):

        self.batch_index = random.sample(range(0, self.train.input.count), self.train.input.count)
        self.index_in_epoch = 0
        self.training_psnr_sum = 0
        self.training_step = 0

    def build_conv_and_bias(self, name, input_tensor, cnn_size, input_feature_num, output_feature_num,
                            use_activator=True,
                            use_dropout=True):
        with tf.variable_scope(name):
            w = util.weight([cnn_size, cnn_size, input_feature_num, output_feature_num],
                            stddev=self.weight_dev, name="conv_W", initializer=self.initializer)
            b = util.bias([output_feature_num], name="conv_B")
            h = self.conv2d(input_tensor, w, self.cnn_stride, bias=b,
                            activator=self.activator if use_activator else None,
                            name=name)

            if use_dropout and self.dropout != 1.0:
                h = tf.nn.dropout(h, self.dropout_input, name="dropout")

            if self.save_weights:
                util.add_summaries("weight", self.name, w, save_stddev=True, save_mean=True)
                util.add_summaries("bias", self.name, b, save_stddev=True, save_mean=True)

            if self.save_images and cnn_size > 1 and input_feature_num == 1:
                weight_transposed = tf.transpose(w, [3, 0, 1, 2])
                with tf.name_scope("image"):
                    tf.summary.image(self.name, weight_transposed, max_outputs=self.log_weight_image_num)

        return w, b, h

    def build_conv(self, name, input_tensor, cnn_size, input_feature_num, output_feature_num):
        with tf.variable_scope(name):
            w = util.weight([cnn_size, cnn_size, input_feature_num, output_feature_num],
                            stddev=self.weight_dev, name="conv_W", initializer=self.initializer)
            h = self.conv2d(input_tensor, w, self.cnn_stride, bias=None, activator=None, name=name)

            if self.save_weights:
                util.add_summaries("weight", self.name, w, save_stddev=True, save_mean=True)

            if self.save_images and cnn_size > 1 and input_feature_num == 1:
                weight_transposed = tf.transpose(w, [3, 0, 1, 2])
                with tf.name_scope("image"):
                    tf.summary.image(self.name, weight_transposed, max_outputs=self.log_weight_image_num)

        return w, h

    def build_input_batch(self, batch_dir):

        for i in range(self.batch_num):
            if self.index_in_epoch >= self.train.input.count:
                self.init_epoch_index()
                self.epochs_completed += 1

            image_no = self.batch_index[self.index_in_epoch]
            self.batch_input[i] = util.load_image(batch_dir + "/" + INPUT_IMAGE_DIR + "/%06d.bmp" % image_no,
                                                  print_console=False)
            batch_input_quad = util.load_image(batch_dir + "/" + INTERPOLATED_IMAGE_DIR + "/%06d.bmp" % image_no,
                                               print_console=False)
            util.convert_to_multi_channel_image(self.batch_input_quad[i], batch_input_quad, self.scale)
            batch_true_quad = util.load_image(batch_dir + "/" + TRUE_IMAGE_DIR + "/%06d.bmp" % image_no,
                                              print_console=False)
            util.convert_to_multi_channel_image(self.batch_true_quad[i], batch_true_quad, self.scale)
            self.index_in_epoch += 1

    def build_graph(self):

        input_feature_num = self.channels

        self.x = tf.placeholder(tf.float32, shape=[None, None, None, input_feature_num], name="X")
        self.y = tf.placeholder(tf.float32, shape=[None, None, None, self.output_channels], name="Y")
        self.x2 = tf.placeholder(tf.float32, shape=[None, None, None, self.output_channels], name="X2")
        self.dropout_input = tf.placeholder(tf.float32, shape=[], name="dropout_keep_rate")

        # building feature extraction layers

        if self.nin_filters == 0:
            self.W_conv = (self.layers + 1) * [None]
            self.B_conv = self.layers * [None]
        else:
            if self.nin:
                self.W_conv = (self.layers + 4) * [None]
                self.B_conv = (self.layers + 3) * [None]
                self.H2 = 2 * [None]
            else:
                self.W_conv = (self.layers + 2) * [None]
                self.B_conv = (self.layers + 1) * [None]

        self.H = self.layers * [None]

        output_feature_num = self.layers * [0]
        total_output_feature_num = 0
        features = ""
        input_feature_num = self.channels
        input_tensor = self.x

        for i in range(self.layers):
            if self.min_filters != 0:
                if i == 0:
                    output_feature_num[i] = self.filters
                else:
                    x1 = i / float(self.layers - 1)
                    y1 = pow(x1, 1.0 / self.filters_decay_gamma)
                    output_feature_num[i] = int((self.filters - self.min_filters) * (1 - y1) + self.min_filters)
            else:
                output_feature_num[i] = self.filters
            total_output_feature_num += output_feature_num[i]
            features += "%d " % output_feature_num[i]

            self.W_conv[i], self.B_conv[i], self.H[i] = self.build_conv_and_bias("conv%d" % i,
                                                                                input_tensor, self.cnn_size,
                                                                                input_feature_num,
                                                                                output_feature_num[i])
            input_feature_num = output_feature_num[i]
            input_tensor = self.H[i]

        with tf.variable_scope("concat"):
            self.H_concat = tf.concat(self.H, 3, name="H_concat")
        features += " Total: (%d)" % total_output_feature_num

        # building reconstruction layers ---

        if self.nin_filters == 0:
            self.W_conv[self.layers], self.H_out = \
                self.build_conv("L", self.H_concat, 1, total_output_feature_num, self.output_channels)
        else:
            if self.nin:
                self.W_conv[self.layers], self.B_conv[self.layers], self.H2[0] = \
                    self.build_conv_and_bias("A1", self.H_concat, 1, total_output_feature_num, self.nin_filters)

                self.W_conv[self.layers + 1], self.B_conv[self.layers + 1], self.H_B1 = \
                    self.build_conv_and_bias("B1", self.H_concat, 1, total_output_feature_num, self.nin_filters2)

                self.W_conv[self.layers + 2], self.B_conv[self.layers + 2], self.H2[1] = \
                    self.build_conv_and_bias("B2", self.H_B1, 3, self.nin_filters2, self.nin_filters2)

                self.H_concat2 = tf.concat(self.H2, 3, name="H_concat2")

                self.W_conv[self.layers + 3], self.H_out = \
                    self.build_conv("L", self.H_concat2, self.last_cnn_size, self.nin_filters + self.nin_filters2,
                                    self.output_channels)
            else:
                self.W_conv[self.layers], self.B_conv[self.layers], self.H_node = \
                    self.build_conv_and_bias("A1", self.H_concat, 1, total_output_feature_num, self.nin_filters)

                self.W_conv[self.layers + 1], self.H_out = \
                    self.build_conv("L", self.H_node, 1, self.nin_filters, self.output_channels)

        self.y_ = self.H_out + self.x2
        self.weights = self.W_conv

        logging.info("Feature:%s Complexity:%s Receptive Fields:%d" % (
            features, "{:,}".format(self.complexity), self.receptive_fields))

    def conv2d(self, x, w, stride, bias=None, activator=None, leaky_relu_alpha=0.1, name=""):
        conv = tf.nn.conv2d(x, w, strides=[stride, stride, 1, 1], padding="SAME", name=name + "_conv")

        self.complexity += int(w.shape[0] * w.shape[1] * w.shape[2] * w.shape[3])

        if bias is not None:
            conv = tf.add(conv, bias, name=name + "_add")
            self.complexity += int(bias.shape[0])

        if activator is not None:
            if activator == "relu":
                conv = tf.nn.relu(conv, name=name + "_relu")
            elif activator == "sigmoid":
                conv = tf.nn.sigmoid(conv, name=name + "_sigmoid")
            elif activator == "tanh":
                conv = tf.nn.tanh(conv, name=name + "_tanh")
            elif activator == "leaky_relu":
                conv = tf.maximum(conv, leaky_relu_alpha * conv, name=name + "_leaky")
            elif activator == "prelu":
                with tf.variable_scope("prelu"):
                    alphas = tf.Variable(tf.constant(0.1, shape=[w.get_shape()[3]]), name=name + "_prelu")
                    if self.save_weights:
                        util.add_summaries("prelu_alpha", self.name, alphas, save_stddev=False, save_mean=False)
                    conv = tf.nn.relu(conv) + tf.multiply(alphas, (conv - tf.abs(conv))) * 0.5
            else:
                raise NameError('Not implemented activator:%s' % activator)
            self.complexity += int(bias.shape[0])

        return conv

    def build_optimizer(self):

        self.lr_input = tf.placeholder(tf.float32, shape=[], name="LearningRate")

        diff = self.y_ - self.y

        mse = tf.reduce_mean(tf.square(diff), name="mse")

        if self.debug:
            mse = tf.Print(mse, [mse], message="MSE: ")

        loss = mse

        if self.l2_decay > 0:
            l2_losses = [tf.nn.l2_loss(w) for w in self.weights]
            l2_loss = self.l2_decay * tf.add_n(l2_losses)
            loss += l2_loss
            if self.save_loss:
                tf.summary.scalar("l2_loss/" + self.name, l2_loss)

        if self.save_loss:
            tf.summary.scalar("test_PSNR/" + self.name, self.get_psnr_tensor(mse))
            tf.summary.scalar("test_loss/" + self.name, loss)

        self.loss = loss
        self.mse = mse
        self.training_optimizer = self.add_optimizer_op(loss, self.lr_input)

        util.print_num_of_total_parameters(output_detail=True)

    def get_psnr_tensor(self, mse):

        with tf.variable_scope('get_PSNR'):
            value = tf.constant(self.max_value, dtype=mse.dtype) / tf.sqrt(mse)
            numerator = tf.log(value)
            denominator = tf.log(tf.constant(10, dtype=mse.dtype))
            return tf.constant(20, dtype=mse.dtype) * numerator / denominator

    def add_optimizer_op(self, loss, lr_input):

        if self.optimizer == "gd":
            optimizer = tf.train.GradientDescentOptimizer(lr_input)
        elif self.optimizer == "adadelta":
            optimizer = tf.train.AdadeltaOptimizer(lr_input)
        elif self.optimizer == "adagrad":
            optimizer = tf.train.AdagradOptimizer(lr_input)
        elif self.optimizer == "adam":
            optimizer = tf.train.AdamOptimizer(lr_input, beta1=self.beta1, beta2=self.beta2)
        elif self.optimizer == "momentum":
            optimizer = tf.train.MomentumOptimizer(lr_input, self.momentum)
        elif self.optimizer == "rmsprop":
            optimizer = tf.train.RMSPropOptimizer(lr_input, momentum=self.momentum)
        else:
            print("Optimizer arg should be one of [gd, adadelta, adagrad, adam, momentum, rmsprop].")
            return None

        training_optimizer = optimizer.minimize(loss)

        return training_optimizer

    def build_summary_saver(self):
        if self.save_loss or self.save_weights:
            self.summary_op = tf.summary.merge_all()
            self.summary_writer = tf.summary.FileWriter(self.tf_log_dir, graph=self.sess.graph)

        self.saver = tf.train.Saver(max_to_keep=None)

    def init_all_variables(self, load_model_name="", trial=0):

        if load_model_name is "":
            self.sess.run(tf.global_variables_initializer())
            print("Model Initialised.")
        else:
            self.load_model(load_model_name, trial)

    def load_model(self, load_model_name="default", trial=0, output_log=True):

        if load_model_name == "default":
            load_model_name = self.name + "_" + str(trial)

        filename = self.checkpoint_dir + "/" + load_model_name + ".ckpt"

        if not os.path.isfile(filename + ".index"):
            print("Error. Model[%s] is not exist!" % filename)
            exit(-1)

        self.saver.restore(self.sess, filename)
        if output_log:
            logging.info("Model restored [ %s ]." % filename)

    def train_batch(self):

        _, mse = self.sess.run([self.training_optimizer, self.mse], feed_dict={self.x: self.batch_input,
                                                                               self.x2: self.batch_input_quad,
                                                                               self.y: self.batch_true_quad,
                                                                               self.lr_input: self.lr,
                                                                               self.dropout_input: self.dropout})
        self.training_psnr_sum += util.get_psnr(mse, max_value=self.max_value)
        self.training_step += 1
        self.step += 1

    def evaluate(self):

        if self.save_loss or self.save_weights:
            summary_str, mse = self.sess.run([self.summary_op, self.mse],
                                             feed_dict={self.x: self.test.input.images,
                                                        self.x2: self.test.input.quad_images,
                                                        self.y: self.test.true.quad_images,
                                                        self.dropout_input: 1.0})

            self.summary_writer.add_summary(summary_str, self.step)
            self.summary_writer.flush()
        else:
            mse = self.sess.run(self.mse,
                                feed_dict={self.x: self.test.input.images,
                                           self.x2: self.test.input.quad_images,
                                           self.y: self.test.true.quad_images,
                                           self.dropout_input: 1.0})
        return mse

    def update_epoch_and_lr(self, mse):
        lr_updated = False

        if self.min_validation_mse < 0 or self.min_validation_mse > mse:
            # update new mse
            self.min_validation_epoch = self.epochs_completed
            self.min_validation_mse = mse
        else:
            if self.epochs_completed > self.min_validation_epoch + self.lr_decay_epoch:
                # set new learning rate
                self.min_validation_epoch = self.epochs_completed
                self.lr *= self.lr_decay
                lr_updated = True

        psnr = util.get_psnr(mse, max_value=self.max_value)
        self.csv_epochs.append(self.epochs_completed)
        self.csv_psnr.append(psnr)
        self.csv_training_psnr.append(self.training_psnr_sum / self.training_step)

        return lr_updated

    def save_summary(self):

        summary_str = self.sess.run(self.summary_op,
                                    feed_dict={self.x: self.test.input.images,
                                               self.x2: self.test.input.quad_images,
                                               self.y: self.test.true.quad_images})

        self.summary_writer.add_summary(summary_str, 0)
        self.summary_writer.flush()

    def print_status(self, mse):

        psnr = util.get_psnr(mse, max_value=self.max_value)

        if self.step == 0:
            print("Initial MSE:%f PSNR:%f" % (mse, psnr))
        else:
            processing_time = (time.time() - self.start_time) / self.step
            print("%s Step:%d MSE:%f PSNR:%f (Training PSNR:%0.3f)" % (
                util.get_now_date(), self.step, mse, psnr, self.training_psnr_sum / self.training_step))
            print("Epoch:%d (Step:%s) LR:%f (%2.3fsec/step) MinPSNR:%0.3f" % (
                self.epochs_completed, "{:,}".format(self.step), self.lr, processing_time,
                util.get_psnr(self.min_validation_mse)))

    def print_weight_variables(self):

        for bias in self.B_conv:
            util.print_filter_biases(bias)

        for weight in self.W_conv:
            util.print_filter_weights(weight)

    def save_model(self, model_dir, trial=0, name=None):

        if name is None:
            name = self.name

        filename = model_dir + "/" + name + "_" + str(trial) + ".ckpt"
        self.saver.save(self.sess, filename)
        print("Model saved [%s]." % filename)

    def save_graphs(self, checkpoint_dir, trial):

        psnr_graph = np.column_stack((self.csv_epochs, self.csv_psnr, self.csv_training_psnr))
        filename = checkpoint_dir + "/" + self.name + ("_%d.csv" % trial)
        np.savetxt(filename, psnr_graph, delimiter=",")

        filename2 = checkpoint_dir + "/" + self.name + ("_%d.png" % trial)
        plt.plot(self.csv_epochs, self.csv_training_psnr, "b", label='Training PSNR')
        plt.plot(self.csv_epochs, self.csv_psnr, "r", label='Test PSNR')
        plt.vlines(self.csv_epochs[-1], 0, self.csv_psnr[-1], color='0.75')
        plt.hlines(self.csv_psnr[-1], 0, self.csv_epochs[-1], color='0.75')
        plt.ylim(ymin=30)

        if trial == 0:
            plt.legend(loc='lower right')
        plt.savefig(filename2)
        plt.savefig("PSNR.png")

        print("Graph saved [%s / %s]." % (filename, filename2))

    def do(self, input_image, bicubic_input_image=None):

        h, w = input_image.shape[:2]
        ch = input_image.shape[2] if len(input_image.shape) > 2 else 1

        if self.max_value != 255.0:
            input_image = np.multiply(input_image, self.max_value / 255.0)  # type: np.ndarray

        if bicubic_input_image is not None:
            bicubic_image = bicubic_input_image
        elif self.bicubic_init:
            bicubic_image = util.resize_image_by_pil(input_image, self.scale, resampling_method="bicubic")
        else:
            bicubic_image = util.resize_image_by_pil(input_image, self.scale, resampling_method="nearest")

        quad_image = np.zeros([1, h, w, self.output_channels])  # type: np.ndarray
        util.convert_to_multi_channel_image(quad_image[0], bicubic_image, self.scale)

        y = self.sess.run(self.y_, feed_dict={self.x: input_image.reshape(1, h, w, ch), self.x2: quad_image,
                                              self.dropout_input: 1.0})

        if self.max_value != 255.0:
            quad_image = np.multiply(y[0], 255.0 / self.max_value)
        else:
            quad_image = y[0]
        image = np.zeros(shape=[h * self.scale, w * self.scale, 1])  # type: np.ndarray
        util.convert_from_multi_channel_image(image, quad_image, self.scale)

        return image

    def do_for_file(self, file_path, output_folder="output"):

        filename, extension = os.path.splitext(file_path)
        output_folder += "/"
        org_image = util.load_image(file_path)
        util.save_image(output_folder + file_path, org_image)

        if len(org_image.shape) >= 3 and org_image.shape[2] == 3 and self.channels == 1:
            input_y_image = util.convert_rgb_to_y(org_image, jpeg_mode=self.jpeg_mode)
            scaled_image = util.resize_image_by_pil(input_y_image, self.scale)
            util.save_image(output_folder + filename + "_bicubic_y" + extension, scaled_image)
            output_y_image = self.do(input_y_image)
            util.save_image(output_folder + filename + "_result_y" + extension, output_y_image)

            scaled_ycbcr_image = util.convert_rgb_to_ycbcr(util.resize_image_by_pil(org_image, self.scale),
                                                           jpeg_mode=self.jpeg_mode)
            image = util.convert_y_and_cbcr_to_rgb(output_y_image, scaled_ycbcr_image[:, :, 1:3],
                                                   jpeg_mode=self.jpeg_mode)
        else:
            scaled_image = util.resize_image_by_pil(org_image, self.scale)
            util.save_image(output_folder + filename + "_bicubic_y" + extension, scaled_image)
            image = self.do(org_image)

        util.save_image(output_folder + filename + "_result" + extension, image)
        return 0

    def do_for_evaluate(self, file_path, output_directory="output", output=True, print_console=True):

        filename, extension = os.path.splitext(file_path)
        output_directory += "/"
        true_image = util.set_image_alignment(util.load_image(file_path), self.scale)

        if true_image.shape[2] == 3 and self.channels == 1:
            input_y_image = util.build_input_image(true_image, channels=self.channels, scale=self.scale,
                                                     alignment=self.scale, convert_ycbcr=True, jpeg_mode=self.jpeg_mode)
            # for color images
            if output:
                input_bicubic_y_image = util.resize_image_by_pil(input_y_image, self.scale)
                true_ycbcr_image = util.convert_rgb_to_ycbcr(true_image, jpeg_mode=self.jpeg_mode)

                output_y_image = self.do(input_y_image, input_bicubic_y_image)
                psnr, ssim = util.compute_psnr_and_ssim(true_ycbcr_image[:, :, 0:1], output_y_image, border_size=self.scale)
                loss_image = util.get_loss_image(true_ycbcr_image[:, :, 0:1], output_y_image, border_size=self.scale)

                output_color_image = util.convert_y_and_cbcr_to_rgb(output_y_image, true_ycbcr_image[:, :, 1:3],
                                                                    jpeg_mode=self.jpeg_mode)

                util.save_image(output_directory + file_path, true_image)
                util.save_image(output_directory + filename + "_input" + extension, input_y_image)
                util.save_image(output_directory + filename + "_input_bicubic" + extension, input_bicubic_y_image)
                util.save_image(output_directory + filename + "_true_y" + extension, true_ycbcr_image[:, :, 0:1])
                util.save_image(output_directory + filename + "_result" + extension, output_y_image)
                util.save_image(output_directory + filename + "_result_c" + extension, output_color_image)
                util.save_image(output_directory + filename + "_loss" + extension, loss_image)
            else:
                true_y_image = util.convert_rgb_to_y(true_image, jpeg_mode=self.jpeg_mode)
                output_y_image = self.do(input_y_image)
                psnr, ssim = util.compute_psnr_and_ssim(true_y_image, output_y_image, border_size=self.scale)

        elif true_image.shape[2] == 1 and self.channels == 1:

            # for monochrome images
            input_image = util.build_input_image(true_image, channels=self.channels, scale=self.scale,
                                                   alignment=self.scale)
            output_image = self.do(input_image)
            psnr, ssim = util.compute_psnr_and_ssim(true_image, output_image, border_size=self.scale)
            if output:
                util.save_image(output_directory + file_path, true_image)
                util.save_image(output_directory + filename + "_result" + extension, output_image)
        else:
            psnr = ssim = None

        if print_console:
            print("PSNR:%f SSIM:%f" % (psnr, ssim))
        return psnr, ssim

    def init_train_step(self):
        self.lr = self.initial_lr
        self.csv_epochs = []
        self.csv_psnr = []
        self.csv_training_psnr = []
        self.epochs_completed = 0
        self.min_validation_mse = -1
        self.min_validation_epoch = -1
        self.step = 0

        self.start_time = time.time()

    def end_train_step(self):
        self.total_time = time.time() - self.start_time

    def print_steps_completed(self, output_to_logging=False):

        if self.step == 0:
            return

        processing_time = self.total_time / self.step
        h = self.total_time // (60 * 60)
        m = (self.total_time - h * 60 * 60) // 60
        s = (self.total_time - h * 60 * 60 - m * 60)

        status = "Finished at Total Epoch:%d Steps:%s Time:%02d:%02d:%02d (%2.3fsec/step)" % (
            self.epochs_completed, "{:,}".format(self.step), h, m, s, processing_time)

        if output_to_logging:
            logging.info(status)
        else:
            print(status)


text:
#!/usr/bin/env python
# coding: utf-8

# # BCG Gamma Challenge

# # Libraries

# In[1]:

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from scipy import stats

# In[2]:

pd.set_option('display.max_rows', 500)
pd.set_option('display.max_columns', 500)

# # Dataset

# In[3]:

df_municipios_2015 = pd.read_csv('bcggammachallenge/municipios/municipios20150101.csv')

# In[4]:

df_municipios_2016 = pd.read_csv('bcggammachallenge/municipios/municipios20160101.csv')

# In[5]:

df_municipios_2017 = pd.read_csv('bcggammachallenge/municipios/municipios20170101.csv')

# In[6]:

df_municipios_2015.shape

# In[7]:

df_municipios_2016.shape

# In[8]:

df_municipios_2017.shape

# In[9]:

df = pd.concat([df_municipios_2015, df_municipios_2016, df_municipios_2017])

# In[10]:

df.shape

# In[11]:

df.head(5)

# In[12]:

df_ideb_ini = pd.read_csv('bcggammachallenge/ideb/ideb_municipios_anosiniciais2005_2017.csv',sep = ',',encoding='latin-1')

# In[13]:

df_ideb_ini.columns

# In[14]:

df_ideb_ini[['Cod_Municipio_Completo', 'Ideb2017']].head()

# In[15]:

df_ideb_ini = df_ideb_ini.rename(columns={'Cod_Municipio_Completo': 'cod_municipio'})

# In[16]:

df_ideb_ini_2015 = df_ideb_ini.copy()

# In[17]:

df_ideb_ini = df_ideb_ini[['cod_municipio', 'Ideb2017']]

# In[18]:

df_ideb_ini_2015 = df_ideb_ini_2015[['cod_municipio', 'Ideb2015']]

# In[19]:

df_ideb_ini.head()

# In[20]:

df_ideb_ini_2015.head()

# In[21]:

df_ideb_ini['cod_municipio'] = df_ideb_ini.cod_municipio.astype(float)
df_ideb_ini_2015['cod_municipio'] = df_ideb_ini_2015.cod_municipio.astype(float)

# In[22]:

df_result_2017 = pd.merge(df_municipios_2017, df_ideb_ini, how='inner', on='cod_municipio')
df_result_2015 = pd.merge(df_municipios_2015, df_ideb_ini_2015, how='inner', on='cod_municipio')

# In[23]:

df_result_2017.head()

# In[ ]:

# In[24]:

df_result_2015 = df_result_2015.rename(columns={'Ideb2015': 'ideb'})
df_result_2017 = df_result_2017.rename(columns={'Ideb2017': 'ideb'})

# In[25]:

df_result_2015.shape

# In[26]:

df_result_2017.shape

# In[27]:

df_result_2015.head()

# In[28]:

df_result_2017.head()

# In[30]:

df = pd.concat([df_result_2015, df_result_2017])

# In[31]:

df.shape

# In[56]:

df['ano_censo'].head()

# In[80]:

df['ideb'].head()

# In[97]:

df['ideb'] = df.ideb.apply(lambda x: x.replace('-',''))

# In[98]:

df['ideb'] = df['ideb'].replace('', 0)

# In[84]:

df_ideb[df_ideb['Ideb2017'] != '-']['Ideb2017'].order

# In[ ]:

df_ideb[df_ideb['Ideb2017'] != '-']['Ideb2017'].min()

# In[ ]:

df_ideb.sort_values(by=['Ideb2017'], ascending=False).head(10)

# In[ ]:

ext: py

max_stars_repo_name: alexandrerays/E-xray

max_issues_repo_name: alexandrerays/E-xray

max_forks_repo_name: alexandrerays/E-xray

reason: import numpy,from scipy


text:
import abc
from math import pi

import numpy as np
import tsplib95

reason: import numpy


text:
import tkinter as tk
from matplotlib import pyplot as plt
import numpy as np
import math
import random
from PIL import ImageTk, Image
import load_model_v2 as load_model
import time
from threading import Timer

# 可以在這邊設定繪圖畫板長寬，範例畫板是16:9
width = 960
height = 540
topic_class = ['飛機', '蝙蝠', '熊', '鳥', '蠟燭', '車', '大象', '眼睛', '花朵', '冰淇淋', '鑰匙', '棒棒糖', '老鼠', '嘴巴', '鼻子', '豬', '鯊魚', '草莓', '太陽', '電視']
topic_class_eng = ['airplane', 'bat', 'bear', 'bird', 'candle', 'car', 'elephant', 'eye', 'flower', 'ice_cream', 'key', 'lollipop', 'mouse', 'mouth', 'nose', 'pig', 'shark', 'strawberry', 'sun', 'television']
dictionary = dict(zip(topic_class_eng, topic_class))

# #########主畫面TK區
class SampleApp(tk.Tk):
    def __init__(self):
        tk.Tk.__init__(self)
        # super().__init__(self) # 太詭異了吧 為什麼這個不行上面可以
        self._frame = None
        self.resizable(False, False)
        self.switch_frame(StartPage)
        self.title('造畫弄人!!')

    def switch_frame(self, frame_class):
        # frame_class 要繼承自tk.Frame，傳入自己(i.e. SampleApp， 也就是tk.Tk())建構
        new_frame = frame_class(self)
        if self._frame is not None:
            self._frame.destroy()
        self._frame = new_frame
        self._frame.pack()
# ##########主畫面TK區

# ##########首頁Frame

class StartPage(tk.Frame):
    def __init__(self, master):
        global width, height
        tk.Frame.__init__(self, master)
        self.ref = Image.open('image/main.png')
        self.ref = self.ref.resize((width, height), Image.BILINEAR)
        self.bg_image = ImageTk.PhotoImage(self.ref)

        # tk.Label(self,image = bg_image,width = width, height = height*1.2).pack()
        tk.Label(self, image=self.bg_image).pack()

        tk.Button(self, text='Game Start', font=('Arial', 25), fg='Blue',
                  command=lambda: master.switch_frame(MainPage)).pack()
# ##########首頁Frame

# ##########遊戲頁Frame

class MainPage(tk.Frame):
    # constructor
    def __init__(self, master):
        tk.Frame.__init__(self, master)
        self.stroke = 0
        self.draw = [[[], []]]
        self.predictions_stack = []
        self.guessed = []
        self.guessed_index = 0
        self.construct()
        self.t = Timer(20, self.times_up)
        self.t.start()
    # constructor

    # 計時器(背景執行緒)
    def times_up(self):
        self.master.switch_frame(FailureEndingPage)
    # 計時器(背景執行緒)

    # 畫畫中
    def drawing(self, event):
        global width, height
        color = 'Black'
        x, y = event.x, event.y
        self.canvas.create_oval((x-12), (y-12), (x+12), (y+12), fill=color)
        self.draw[self.stroke][0].append(x*256//(width-100))
        self.draw[self.stroke][1].append(y*256//(height-30))
    # 畫畫中

    # 提起畫筆
    def end_stroke(self, event):
        global dictionary
        if self.draw[self.stroke] != [[], []]:
            self.stroke += 1
            self.draw.append([[], []])
        if (self.stroke > 0) and (self.stroke % 3 == 0):
            self.auto_predict()

        if self.predictions_stack:
            guess = dictionary[self.predictions_stack[self.guessed_index]]
            if guess not in self.guessed:
                self.answer_lbl.config(
                    text=f'我猜你畫的是...{guess}', font=('微軟正黑體', 20))
                self.guessed_index += 1
                self.guessed_index %= 3
            if self.answer == guess:
                print('猜中了')
                self.t.cancel()
                time.sleep(0.5)
                self.master.switch_frame(EndingPage)
            else:
                self.guessed.append(guess)
                print(f'沒中, guessed = {self.guessed}')

    # 提起畫筆

    # 清除畫布
    def clear(self):
        self.canvas.delete('all')
        self.stroke = 0
        self.draw = [[[], []]]
        self.answer_lbl.config(text='...')
    # 清除畫布

    # 換題目
    def change_topic(self):
        ind = random.randint(1, len(topic_class))
        self.answer = topic_class[ind-1]
        self.question.config(text=f'題目: {self.answer}')
        self.clear()
        self.t.cancel()
        self.t = Timer(20, self.times_up)
        self.t.start()
    # 換題目

    # 跑預測
    def auto_predict(self):
        global predictions
        if self.isListEmpty(self.draw):
            return
        # 可能
        self.draw.pop()

        self.predictions_stack = load_model.predict(self.draw)

        self.draw.append([[], []])
        # self.master.switch_frame(PredictionPage)
    # 跑預測

    # 檢查陣列是不是空空的
    def isListEmpty(self, inList):
        if isinstance(inList, list):  # Is a list
            return all(map(self.isListEmpty, inList))
        return False
    # 檢查陣列是不是空空的

    def construct(self):
        global width, height
        # 上層區塊
        self.frame_top = tk.Frame(
            self, width=width-100, height=30, bg='#DCDCDC')
        self.frame_top.pack()

        # 清除畫布按鈕
        self.clear_btn = tk.Button(self.frame_top)
        self.clear_btn.config(text='神奇橡皮擦', font=(
            '微軟正黑體', 14), command=self.clear, fg="Blue")
        self.clear_btn.place(relx=0.16, rely=0)
        # 清除畫布按鈕

        # 題目標示
        self.question = tk.Label(self.frame_top)
        # #隨機出題
        ind = random.randint(1, len(topic_class))
        self.answer = topic_class[ind-1]
        # #隨機出題
        self.question.config(
            text=f'題目: {self.answer}', font=('微軟正黑體', 16), fg='Blue')
        self.question.place(relx=0.4, rely=0)
        # 題目標示

        # 切換題目按鈕
        self.submit_btn = tk.Button(self.frame_top)
        self.submit_btn.config(text='重抽題目', font=(
            '微軟正黑體', 14), command=self.change_topic, fg="Blue")
        self.submit_btn.place(relx=0.6, rely=0)
        # 切換題目按鈕
        # 上層區塊

        # 中間繪圖畫布
        self.frame_canvas = tk.Frame(self, width=width-100, height=height-30)
        self.frame_canvas.pack()

        self.canvas = tk.Canvas(
            self.frame_canvas, width=width-100, height=height-30, bg='White')
        self.canvas.place(relx=0, rely=0)
        self.canvas.bind('<Button1-Motion>', self.drawing)  # 按著滑鼠左鍵滑動
        self.canvas.bind('<Button1-ButtonRelease>', self.end_stroke)  # 放開滑鼠左鍵
        # 中間繪圖畫布

        # 下層顯示作答區塊
        self.frame_bottom = tk.Frame(
            self, width=width-100, height=50, bg='#DCDCDC')
        self.frame_bottom.pack()

        self.answer_lbl = tk.Label(
            self.frame_bottom, text='...', font=('Arial', 20), fg='Red')
        self.answer_lbl.pack()
        # 下層顯示作答區塊

# ##########遊戲頁Frame

# ##########看之後要不要做個統計頁面這樣
# class PredictionPage(tk.Frame):
#   def __init__(self, master):
#     global predictions, dictionary
#     tk.Frame.__init__(self, master)
#     self.ref = Image.open('image/ending.png')
#     self.ref = self.ref.resize((400,400), Image.BILINEAR)
#     self.bg_image = ImageTk.PhotoImage(self.ref)
#     tk.Label(self, image = self.bg_image).pack()

#     tk.Label(self, text = f'我們的神經網路判斷你畫的是...{dictionary[predictions[0]]}', font = ('微軟正黑體', 30), fg = 'Red').pack()
#     tk.Label(self, text = f'我們的神經網路覺得你畫的還可能是...{dictionary[predictions[1]]}', font = ('微軟正黑體', 16), fg = 'Green').pack()
#     tk.Label(self, text = f'我們的神經網路覺得你畫的還可能是...{dictionary[predictions[2]]}', font = ('微軟正黑體', 16), fg = 'Green').pack()
#     tk.Label(self, text = f'我們的神經網路覺得你畫的還可能是...{dictionary[predictions[3]]}', font = ('微軟正黑體', 16), fg = 'Green').pack()
#     tk.Label(self, text = f'我們的神經網路覺得你畫的還可能是...{dictionary[predictions[4]]}', font = ('微軟正黑體', 16), fg = 'Green').pack()

#     tk.Button(self, text = '想回到最初的起點', font = ('微軟正黑體', 25), fg = 'Blue', command = lambda:master.switch_frame(StartPage)).pack()
# ##########看之後要不要做個統計頁面這樣

# ##########結束頁Frame
class EndingPage(tk.Frame):
    def __init__(self, master):
        tk.Frame.__init__(self, master)
        tk.Label(self, text='恭喜！你是靈魂繪師', font=('微軟正黑體', 25), fg='Red').pack()
        self.ref = Image.open('image/ending.png')
        self.ref = self.ref.resize((960, 485), Image.BILINEAR)
        self.bg_image = ImageTk.PhotoImage(self.ref)
        tk.Label(self, image=self.bg_image).pack()
        tk.Button(self, text='想回到最初的起點', font=('微軟正黑體', 25), fg='Blue',
                  command=lambda: master.switch_frame(StartPage)).pack()
# ##########結束頁Frame

# ##########失敗頁Frame

class FailureEndingPage(tk.Frame):
    def __init__(self, master):
        tk.Frame.__init__(self, master)
        tk.Label(self, text='施主,您還是多修練修練畫工吧',
                 font=('微軟正黑體', 25), fg='Red').pack()
        self.ref = Image.open('image/ending.png').convert('LA')
        self.ref = self.ref.resize((960, 485), Image.BILINEAR)
        self.bg_image = ImageTk.PhotoImage(self.ref)
        tk.Label(self, image=self.bg_image).pack()
        tk.Button(self, text='想回到最初的起點', font=('微軟正黑體', 25), fg='Blue',
                  command=lambda: master.switch_frame(StartPage)).pack()
# ##########失敗頁Frame

if __name__ == "__main__":
    app = SampleApp()
    app.mainloop()


text:
import TccLib

import numpy as np

# Define boundary conditions and parameters
vel = 1.
dt = 0.1
total_time = 50.

# Set liquid parameters or declare liquid
# density = 1e3
# viscosity = 0.89e-3
liquid = "super_oil"

# Import gmsh created mesh
mesh = TccLib.Mesh("Poiseuille_xref2", liquid=liquid)

# Show mesh geometry
# mesh.show_geometry(names=True, save=True)

# Define analytic comparison expression
analytic_expression = lambda y: 4 * 1.5 * vel / (mesh.length_y**2) * (mesh.length_y - y) * y

# -------------------------------------------------- PSI ---------------------------------------------------------------
# Define boundary conditions for psi
boundary_conditions_values_psi = {
    "north": lambda y: mesh.y[y],
    "east": lambda y: mesh.y[y],
    "south": lambda y: mesh.y[y],
    "west": lambda y: mesh.y[y],
}

boundary_conditions_types_psi = {
    "north": True,
    "east": True,
    "south": True,
    "west": True,
}

# Get vectors of nodes values
xy_indices, xy_values, xy_types = TccLib.build_boundary_conditions(mesh, boundary_conditions_values_psi,
                                                                   boundary_conditions_types_psi)

# Set Boundary Conditions for psi
mesh.new_boundary_condition("psi", point_index=xy_indices, values=xy_values,
                            type_of_boundary=xy_types)

# -------------------------------------------------- VEL X -------------------------------------------------------------
# Define boundary conditions for vel_x
boundary_conditions_values_x = {
    "north": 0.,
    "east": 0.,
    "south": 0.,
    "west": vel,
}

boundary_conditions_types_x = {
    "north": True,
    "east": False,
    "south": True,
    "west": True,
}

# Get vectors of nodes values
xy_indices, xy_values, xy_types = TccLib.build_boundary_conditions(mesh, boundary_conditions_values_x,
                                                                   boundary_conditions_types_x)

# Set Boundary Conditions for vel_x
mesh.new_boundary_condition("vel_x", point_index=xy_indices, values=xy_values,
                            type_of_boundary=xy_types)

# -------------------------------------------------- VEL Y -------------------------------------------------------------
# Define boundary conditions for vel_y
boundary_conditions_values_y = {
    "north": 0.,
    "east": 0.,
    "south": 0.,
    "west": 0.,
}

boundary_conditions_types_y = {
    "north": True,
    "east": False,
    "south": True,
    "west": True,
}

# Get vectors of nodes values
xy_indices, xy_values, xy_types = TccLib.build_boundary_conditions(mesh, boundary_conditions_values_y,
                                                                   boundary_conditions_types_y)

# Set Boundary Conditions for vel_y
mesh.new_boundary_condition("vel_y", point_index=xy_indices, values=xy_values,
                            type_of_boundary=xy_types)

# Solve for FEM velocity field solution
velocity_x, velocity_y = TccLib.solve_velocity_field(mesh, dt=dt, total_time=total_time, save_each_frame=True,
                                                     stop_rule=1e-5)
TccLib.util.save(velocity_x, "vel_x")
velocity_x = TccLib.util.load("vel_x")

# Show results in quiver plot
# mesh.show_velocity_quiver(velocity_x, velocity_y)

# Define x vector of positions
x_vector = np.linspace(min(mesh.y), max(mesh.y), 100)
small_dict = velocity_x.reduced_dict_log(5)

# Find values of property in the mesh at a determined set position for every value in the vector of x
x_position = mesh.length_x*0.6
y_vector = {key: [mesh.get_interpolated_value([x_position, x], small_dict[key]) for x in x_vector]
            for key in small_dict}

# Show comparison graph
TccLib.util.show_comparison(x_vector, analytic_expression, y_vector, numeric_label="Solução Numérica",
                            analytic_label="Solução Analítica", title="Equação de Hagen-Poiseuille",
                            x_label="Posição no Eixo Y(m)", y_label="Valor da Velocidade (m/s)",
                            save_file_as="{0}_validation".format(mesh.name))

max_stars_repo_name: encarvlucas/LucasCarvalhoTCC

max_issues_repo_name: encarvlucas/LucasCarvalhoTCC

max_forks_repo_name: encarvlucas/LucasCarvalhoTCC


text:
import os
import cv2
import torch
import segmentation_models_pytorch as smp
from config import configs
from IPython import embed
import numpy as np
from glob import glob
from tqdm import tqdm
import matplotlib.pyplot as plt
from albumentations import pytorch as AT
import albumentations as albu

def get_valid_augmentation():
    # for valid
    valid_transform = [
        albu.Resize(height=configs.input_size,width=configs.input_size),
        albu.Normalize(mean=(0.485, 0.456, 0.406), std=(0.229, 0.224, 0.225), p=1),
        AT.ToTensor(),
    ]
    return albu.Compose(valid_transform)

model = smp.Unet(
    encoder_name=configs.encoder,
    encoder_weights=configs.encoder_weights,
    classes=configs.num_classes,
    activation=configs.activation)
model.cuda()
model.eval()
state = torch.load("checkpoints/cloud/resnet50-best_dice.pth", map_location=lambda storage, loc: storage)
model.load_state_dict(state["state_dict"])

images = glob("/data/zcj/dataset/detection/cloud/dataset/test/images/*")
#images = glob("./webank/images/*")
for filename in tqdm(images[:20]):
    fname = os.path.basename(filename)
    image = cv2.imread(filename)
    img = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)
    augmented = get_valid_augmentation()(image=img)
    img = augmented['image'].unsqueeze(0)
    preds_batch = model(img.cuda())
    for preds in preds_batch:
        pred_mask = np.argmax(preds.detach().cpu().numpy(),0).astype(np.uint8)
        #cv2.imwrite(configs.pred_mask+"pred_"+fname,pred_mask*255)
        plt.figure(figsize=(18,18))
        plt.subplot(1,2,1)
        plt.imshow(image)
        plt.subplot(1,2,2)
        plt.imshow(pred_mask)
        plt.savefig(configs.pred_mask+"pred_"+fname)
        plt.close()


text:
# -*- coding: utf-8 -*-
"""
This is a function to test the day2week2month functions

This module is part of the iRONS toolbox by A. Peñuela and F. Pianosi and at
Bristol University (2020).
"""
import pandas as pd
import numpy as np
from numpy.testing import assert_array_equal

if __name__ == '__main__':
    import sys
    sys.path.append("..") # Adds higher directory to python modules path.
    ### Function to test ###
    from day2week2month import day2week
else:
    ### Function to test ###
    from irons.Software.day2week2month import day2week

# Test inputs
dates = pd.date_range(start = '2020-06-09', end = '2020-07-09', freq = 'D')
data  = np.ones(dates.size+1)
# Run the function to test
dates_week,data_week,data_cum_week = day2week(dates,data)

### Testing functions ###
def test_dates_week():
    # Expected output
    dates_week_expect = pd.to_datetime(np.array(['2020-06-09', '2020-06-16',
                                                 '2020-06-23', '2020-06-30',
                                                 '2020-07-07']))
    # Test
    assert_array_equal(dates_week,dates_week_expect)

### Testing functions ###
def test_data_week():
    # Expected output
    data_week_expect = np.array([[0],[7],[7],[7],[7]])
    # Test
    assert_array_equal(data_week,data_week_expect)

### Testing functions ###
def test_data_cum_week():
    # Expected output
    data_cum_week_expect = np.array([[0],[7],[14],[21],[28]])
    # Test
    assert_array_equal(data_cum_week,data_cum_week_expect)

max_stars_repo_name: csiro-hydrology/iRONS

max_issues_repo_name: csiro-hydrology/iRONS

max_forks_repo_name: csiro-hydrology/iRONS

reason: import numpy,from numpy


text:
import cntk as C
from cntk.ops.functions import load_model

import numpy as np
import imp
from scipy.spatial.distance import cosine

def lstm(qry, ans, ans2):
  # Get the model
  query_vector = load_model("models/query_vector.model")
  answer_vector = load_model("models/answer_vector.model")
  # Get the data file
  value = getVarFromFile("variables.txt")
  queryf = open(value.query_wf)
  ansf = open(value.answer_wf)

  query_wl = [line.rstrip('\n') for line in queryf]
  answers_wl = [line.rstrip('\n') for line in ansf]
  query_dict = {query_wl[i]:i for i in range(len(query_wl))}
  answers_dict = {answers_wl[i]:i for i in range(len(answers_wl))}
  queryf.close()
  ansf.close()

  # let's run a sequence through
  qry_idx = [query_dict[w+' '] for w in qry.split()] # convert to query word indices
  ans_idx = [answers_dict[w+' '] for w in ans.split()] # convert to answer word indices
  ans2_idx = [answers_dict[w+' '] for w in ans2.split()] # convert to fake answer word indices

  # Create the one hot representations
  qry_onehot = np.zeros([len(qry_idx),len(query_dict)], np.float32)
  for t in range(len(qry_idx)):
    qry_onehot[t,qry_idx[t]] = 1

  ans_onehot = np.zeros([len(ans_idx),len(answers_dict)], np.float32)
  for t in range(len(ans_idx)):
    ans_onehot[t,ans_idx[t]] = 1

  ans2_onehot = np.zeros([len(ans2_idx),len(answers_dict)], np.float32)
  for t in range(len(ans2_idx)):
    ans2_onehot[t, ans2_idx[t]] = 1

  qry_embedding = query_vector.eval([qry_onehot])
  ans_embedding = answer_vector.eval([ans_onehot])
  ans2_embedding = answer_vector.eval([ans2_onehot])
  qry_ans_similarity = 1-cosine(qry_embedding, ans_embedding)
  qry_ans2_similarity = 1-cosine(qry_embedding, ans2_embedding)
  return qry_ans_similarity, qry_ans2_similarity

# Reading from file
def getVarFromFile(filename):
  f = open(filename)
  global data
  data = imp.load_source('data', '', f)
  f.close()
  return data

max_stars_repo_name: EskenderTamrat/DSSM-LSTM


text:
"""
This the main estimator module
"""
import numpy as np
import pandas as pd

import warnings
from time import time

from sklearn.base import BaseEstimator, DensityMixin
from sklearn.utils.validation import check_is_fitted
from sklearn.utils import check_random_state, check_array
from sklearn.exceptions import ConvergenceWarning

from .random_matrix import generateRandomDefPosMat
from ._euler_model import EulerForceVisibleModel

from ._gle_basis_projection import GLE_BasisTransform

# In case the fortran module is not available, there is the python fallback
try:
    from ._kalman_python import filtersmoother
except ImportError as err:
    print(err)
    warnings.warn("Python fallback will been used for filtersmoother module. Consider compiling the fortran module")
    from ._kalman_python import filtersmoother

import multiprocessing

model_class = {"euler": EulerForceVisibleModel}

def sufficient_stats(traj, dim_x):
    """
    Given a sample of trajectory, compute the averaged values of the sufficient statistics
    Datas are stacked as (xv_plus_proj, xv_proj, v, bk)
    """

    xval = traj[:-1, 2 * dim_x : 3 * dim_x]
    dx = traj[:-1, :dim_x] - traj[:-1, dim_x : 2 * dim_x]
    bk = traj[:-1, 3 * dim_x :]
    xx = np.mean(xval[:, :, np.newaxis] * xval[:, np.newaxis, :], axis=0)
    xdx = np.mean(xval[:, :, np.newaxis] * dx[:, np.newaxis, :], axis=0)
    dxdx = np.mean(dx[:, :, np.newaxis] * dx[:, np.newaxis, :], axis=0)
    bkx = np.mean(bk[:, :, np.newaxis] * xval[:, np.newaxis, :], axis=0)
    bkdx = np.mean(bk[:, :, np.newaxis] * dx[:, np.newaxis, :], axis=0)
    bkbk = np.mean(bk[:, :, np.newaxis] * bk[:, np.newaxis, :], axis=0)

    return pd.Series({"dxdx": dxdx, "xdx": xdx, "xx": xx, "bkx": bkx, "bkdx": bkdx, "bkbk": bkbk, "µ_0": 0, "Σ_0": 1, "hS": 0})

def sufficient_stats_hidden(muh, Sigh, traj, old_stats, dim_x, dim_h, dim_force, model="aboba"):
    """
    Compute the sufficient statistics averaged over the hidden variable distribution
    Datas are stacked as (xv_plus_proj, xv_proj, v, bk)
    """
    # print("Suff_stats")
    xx = np.zeros((dim_x + dim_h, dim_x + dim_h))
    xx[:dim_x, :dim_x] = old_stats["xx"]
    xdx = np.zeros_like(xx)
    xdx[:dim_x, :dim_x] = old_stats["xdx"]
    dxdx = np.zeros_like(xx)
    dxdx[:dim_x, :dim_x] = old_stats["dxdx"]
    bkx = np.zeros((dim_force, dim_x + dim_h))
    bkx[:, :dim_x] = old_stats["bkx"]
    bkdx = np.zeros_like(bkx)
    bkdx[:, :dim_x] = old_stats["bkdx"]

    xval = traj[:-1, 2 * dim_x : 3 * dim_x]
    dx = traj[:-1, :dim_x] - traj[:-1, dim_x : 2 * dim_x]
    bk = traj[:-1, 3 * dim_x :]

    dh = muh[:-1, :dim_h] - muh[:-1, dim_h:]

    Sigh_tptp = np.mean(Sigh[:-1, :dim_h, :dim_h], axis=0)
    Sigh_ttp = np.mean(Sigh[:-1, dim_h:, :dim_h], axis=0)
    Sigh_tpt = np.mean(Sigh[:-1, :dim_h, dim_h:], axis=0)
    Sigh_tt = np.mean(Sigh[:-1, dim_h:, dim_h:], axis=0)

    muh_tptp = np.mean(muh[:-1, :dim_h, np.newaxis] * muh[:-1, np.newaxis, :dim_h], axis=0)
    muh_ttp = np.mean(muh[:-1, dim_h:, np.newaxis] * muh[:-1, np.newaxis, :dim_h], axis=0)
    muh_tpt = np.mean(muh[:-1, :dim_h, np.newaxis] * muh[:-1, np.newaxis, dim_h:], axis=0)
    muh_tt = np.mean(muh[:-1, dim_h:, np.newaxis] * muh[:-1, np.newaxis, dim_h:], axis=0)

    xx[dim_x:, dim_x:] = Sigh_tt + muh_tt
    xx[dim_x:, :dim_x] = np.mean(muh[:-1, dim_h:, np.newaxis] * xval[:, np.newaxis, :], axis=0)

    xdx[dim_x:, dim_x:] = Sigh_ttp + muh_ttp - Sigh_tt - muh_tt
    xdx[dim_x:, :dim_x] = np.mean(muh[:-1, dim_h:, np.newaxis] * dx[:, np.newaxis, :], axis=0)
    xdx[:dim_x, dim_x:] = np.mean(xval[:, :, np.newaxis] * dh[:, np.newaxis, :], axis=0)

    dxdx[dim_x:, dim_x:] = Sigh_tptp + muh_tptp - Sigh_ttp - Sigh_tpt - muh_ttp - muh_tpt + Sigh_tt + muh_tt
    dxdx[dim_x:, :dim_x] = np.mean(dh[:, :, np.newaxis] * dx[:, np.newaxis, :], axis=0)

    bkx[:, dim_x:] = np.mean(bk[:, :, np.newaxis] * muh[:-1, np.newaxis, dim_h:], axis=0)
    bkdx[:, dim_x:] = np.mean(bk[:, :, np.newaxis] * dh[:, np.newaxis, :], axis=0)

    xx[:dim_x, dim_x:] = xx[dim_x:, :dim_x].T
    dxdx[:dim_x, dim_x:] = dxdx[dim_x:, :dim_x].T

    detd = np.linalg.det(Sigh[:-1, :, :])
    dets = np.linalg.det(Sigh[:-1, dim_h:, dim_h:])
    hSdouble = 0.5 * np.log(detd[detd > 0.0]).mean()
    hSsimple = 0.5 * np.log(dets[dets > 0.0]).mean()
    # TODO take care of initial value that is missing
    return pd.Series({"dxdx": dxdx, "xdx": xdx, "xx": xx, "bkx": bkx, "bkdx": bkdx, "bkbk": old_stats["bkbk"], "µ_0": muh[0, dim_h:], "Σ_0": Sigh[0, dim_h:, dim_h:], "hS": 0.5 * dim_h * (1 + np.log(2 * np.pi)) + hSdouble - hSsimple})

def e_step_worker_pool(est, traj, datas_visible, N):
    muh, Sigh = est._e_step(traj)  # Compute hidden variable distribution
    return sufficient_stats_hidden(muh, Sigh, traj, datas_visible, est.dim_x, est.dim_h, est.dim_coeffs_force) / N

class GLE_Estimator(DensityMixin, BaseEstimator):
    """A GLE estimator based on Expectation-Maximation algorithm.

    Parameters
    ----------
    dim_x : int, default=1
        The number of visible dimensions

    dim_h : int, default=1
        The number of hidden dimensions

    tol : float, defaults to 1e-5.
        The convergence threshold. EM iterations will stop when the lower bound average gain is below this threshold.

    max_iter: int, default=100
        The maximum number of EM iterations

    OptimizeForce: bool, default=True
        Optimize or not the force coefficients, to be set to False if the force or the potential have been externally determined

    OptimizeDiffusion: bool, default=True
        Optimize or not the diffusion coefficients

    init_params : {'user','random','markov'}, defaults to 'random'.

        The method used to initialize the fitting coefficients.
        Must be one of::

            'user' : coefficients are initialized at values provided by the user
            'random' : coefficients are initialized randomly.
            'markov' : coefficients are initialized with Markovian estimation of the visible part

    model : {}, default to 'euler'.
        Choice of time discretized model to be fitted. For now only euler model is implemented

    basis: a scikit-learn Transformer class, default to linear basis.
        Transformer to get value of the basis function

    A_init, C_init, force_init, mu_init, sig_init: array, optional
        The user-provided initial coefficients, defaults to None.
        If it None, coefficients are initialized using the `init_params` method.

    n_init : int, defaults to 1.
        The number of initializations to perform. The best results are kept.

    no_stop: bool, default to False
        Does not stop the iterations if the algorithm have converged

    warm_start : bool, default to False.
        If 'warm_start' is True, the solution of the last fitting is used as
        initialization for the next call of fit(). This can speed up
        convergence when fit is called several times on similar problems.

    random_state : int, RandomState instance or None, optional (default=None)
        Controls the random seed given to the method chosen to initialize the
        parameters (see `init_params`).
        In addition, it controls the generation of random samples from the
        fitted distribution (see the method `sample`).
        Pass an int for reproducible output across multiple function calls.
        See :term:`Glossary <random_state>`.

    verbose : int, default to 0.
        Enable verbose output. If 1 then it prints the current
        initialization and each iteration step. If greater than 1 then
        it prints also the log probability and the time needed
        for each step.

    verbose_interval : int, default to 10.
        Number of iteration done before the next print.

    multiprocessing: int, default to 1
        Number of process to use for E step
    """

    def __init__(
        self,
        dim_x=1,
        dim_h=1,
        tol=1e-5,
        max_iter=100,
        OptimizeForce=True,
        OptimizeDiffusion=True,
        init_params="random",
        model="euler",
        basis=GLE_BasisTransform(basis_type="linear"),
        A_init=None,
        C_init=None,
        force_init=None,
        mu_init=None,
        sig_init=None,
        n_init=1,
        random_state=None,
        warm_start=False,
        no_stop=False,
        verbose=0,
        verbose_interval=10,
        multiprocessing=1,
    ):
        self.dim_x = dim_x
        self.dim_h = dim_h

        self.OptimizeForce = OptimizeForce
        self.OptimizeDiffusion = OptimizeDiffusion

        self.model = model
        self.basis = basis

        self.A_init = A_init
        self.C_init = C_init
        self.force_init = force_init
        self.mu_init = mu_init
        self.sig_init = sig_init

        self.tol = tol
        self.max_iter = max_iter
        self.n_init = n_init
        self.warm_start = warm_start

        self.no_stop = no_stop

        self.init_params = init_params
        self.random_state = random_state
        self.verbose = verbose
        self.verbose_interval = verbose_interval
        self.multiprocessing = multiprocessing

    def _more_tags(self):
        return {"X_types": "2darray"}

    def set_init_coeffs(self, coeffs):
        """Set the initial values of the coefficients via a dict

        Parameters
        ----------
        coeffs : dict
        Contains the wanted values of the initial coefficients,
        if a key is absent the coefficients is set to None
        """
        keys_coeffs = ["A", "C", "force", "µ_0", "Σ_0"]
        init_coeffs = [None] * len(keys_coeffs)
        for n, key in enumerate(keys_coeffs):
            if key in coeffs:
                init_coeffs[n] = coeffs[key]
        self.A_init, self.C_init, self.force_init, self.mu_init, self.sig_init = init_coeffs

    def _check_initial_parameters(self):
        """Check values of the basic parameters."""
        if self.dt <= 0.0:
            raise ValueError("Invalid value for 'dt': %d " "Timestep should be positive" % self.dt)
        if self.dim_h < 0:
            raise ValueError("Invalid value for 'dim_h': %d " "Estimator requires non-negative hidden dimension" % self.dim_h)
        if self.tol < 0.0:
            raise ValueError("Invalid value for 'tol': %.5f " "Tolerance used by the EM must be non-negative" % self.tol)
        if self.max_iter < 1:
            raise ValueError("Invalid value for 'max_iter': %d " "Estimation requires at least one iteration" % self.max_iter)

        self.model = self.model.casefold()

        if self.model not in model_class.keys():
            raise ValueError("Model {} not implemented".format(self.model))
        self.model_class = model_class[self.model](self.dim_x)

        if self.init_params == "user":
            if self.n_init != 1:
                self.n_init = 1
                warnings.warn("The number of initialization have been put to 1 as the coefficients are user initialized.")

            if self.A_init is None:
                raise ValueError("No initial values for A is provided and init_params is set to user defined")
            if self.force_init is None:
                raise ValueError("No initial values for the force is provided and init_params is set to user defined")
            # if self.mu_init is None or self.sig_init is None:
            #     raise ValueError("No initial values for initial conditions are provided and init_params is set to user defined")

        if self.A_init is not None:
            if np.asarray(self.A_init).shape != (self.dim_x + self.dim_h, self.dim_x + self.dim_h):
                raise ValueError("Wrong dimensions for A_init")
            if self.C_init is None:
                self.C_init = np.identity(self.dim_x + self.dim_h)
            expA, SST = self.model_class._convert_user_coefficients(np.asarray(self.A_init), np.asarray(self.C_init), self.dt)
            if not np.all(np.linalg.eigvals(SST) > 0):
                raise ValueError("Provided user values does not lead to definite positive diffusion matrix")

        if self.C_init is not None:
            if np.asarray(self.C_init).shape != (self.dim_x + self.dim_h, self.dim_x + self.dim_h):
                raise ValueError("Wrong dimensions for C_init")

        if self.mu_init is not None:
            if np.asarray(self.mu_init).shape != (self.dim_h,):
                raise ValueError("Provided user values for initial mean of hidden variables have wrong shape, provided {}, wanted {}".format(np.asarray(self.mu_init).shape, (self.dim_h,)))

        if self.sig_init is not None:
            if np.asarray(self.sig_init).shape != (self.dim_h, self.dim_h):
                raise ValueError("Provided user values for initial variance of hidden variables have wrong shape, provided {}, wanted {}".format(np.asarray(self.sig_init).shape, (self.dim_h, self.dim_h)))
            if not np.all(np.linalg.eigvals(self.sig_init) >= 0):
                raise ValueError("Provided user values for initial variance of hidden variables is not a definite positive diffusion matrix")

        # We initialize the coefficients value to dummy values to ensure existence of the variables
        if not hasattr(self, "friction_coeffs"):
            self.friction_coeffs = np.identity(self.dim_x + self.dim_h)
        if not hasattr(self, "diffusion_coeffs"):
            self.diffusion_coeffs = np.identity(self.dim_x + self.dim_h)

    def _initialize_parameters(self, random_state, traj_len=50):
        """Initialize the model parameters.

        Parameters
        ----------
        X : array-like, shape  (n_samples, n_features)
        random_state : RandomState
            A random number generator instance that controls the random seed
            used for the method chosen to initialize the parameters.
        """
        self.random_state = check_random_state(random_state)
        if self.init_params == "random" or self.init_params == "markov":
            A = generateRandomDefPosMat(dim_x=self.dim_x, dim_h=self.dim_h, rng=self.random_state, max_ev=(1.0 / 50) / self.dt, min_re_ev=(0.5 / traj_len) / self.dt)  # We ask the typical time scales to be correct with minimum and maximum timescale of the trajectory
            if self.C_init is None:
                # temp_mat = generateRandomDefPosMat(self.dim_h + self.dim_x, random_state)
                # C = temp_mat + temp_mat.T
                C = np.identity(self.dim_x + self.dim_h)
            else:
                C = np.asarray(self.C_init)
            (self.friction_coeffs, self.diffusion_coeffs) = self.model_class._convert_user_coefficients(A, C, self.dt)
        elif self.init_params == "user":
            (self.friction_coeffs, self.diffusion_coeffs) = self.model_class._convert_user_coefficients(np.asarray(self.A_init), np.asarray(self.C_init), self.dt)
            self.force_coeffs = np.asarray(self.force_init).reshape(self.dim_x, -1)
        else:
            raise ValueError("Unimplemented initialization method '%s'" % self.init_params)

        if not self.OptimizeDiffusion and self.A_init is not None and self.C_init is not None:
            _, self.diffusion_coeffs = self.model_class._convert_user_coefficients(np.asarray(self.A_init), np.asarray(self.C_init), self.dt)

        if not hasattr(self.basis, "fitted_"):  # Setup the basis if needed
            dummytraj = np.zeros((1, self.dim_x))
            self.basis.fit(dummytraj)
        self.dim_coeffs_force = self.basis.nb_basis_elt_

        if self.force_init is not None:
            self.force_coeffs = np.asarray(self.force_init).reshape(self.dim_x, -1)
        else:
            self.force_coeffs = -self.random_state.random(size=(self.dim_x, self.dim_coeffs_force))  # -np.ones((self.dim_x, self.dim_coeffs_force))

        # Initial conditions for hidden variables, either user provided or chosen from stationnary state probability fo the hidden variables
        if self.mu_init is not None:
            self.mu0 = np.asarray(self.mu_init)
        else:
            self.mu0 = np.zeros((self.dim_h))

        if self.sig_init is not None:
            self.sig0 = np.asarray(self.sig_init)
        elif self.C_init is not None:
            self.sig0 = self.C_init[self.dim_x :, self.dim_x :]
        else:
            self.sig0 = np.identity(self.dim_h)
        self.initialized_ = True

    def fit(self, X, y=None, idx_trajs=[]):
        """Estimate model parameters with the EM algorithm.
        The method iterates between E-step and M-step for ``max_iter``
        times until the change of likelihood or lower bound is less than
        ``tol``, otherwise, a ``ConvergenceWarning`` is raised.
        Upon consecutive calls, training starts where it left off.

        Parameters
        ----------
        X : array-like, shape (n_samples, dim_x)
            List of positions data.
        idx_trajs: array, default []
            Location of split if multiple trajectory are inputed

        Returns
        -------
        self
        """
        X = check_array(X, ensure_min_samples=4, ensure_min_features=self.dim_x)
        self.dt = X[1, 0] - X[0, 0]
        self._check_initial_parameters()

        Xproc, idx_trajs = self.model_class.preprocessingTraj(self.basis, X, idx_trajs=idx_trajs)
        traj_list = np.split(Xproc, idx_trajs)
        _min_traj_len = np.min([trj.shape[0] for trj in traj_list])
        # print(traj_list)
        # if we enable warm_start, we will have a unique initialisation
        do_init = not (self.warm_start and hasattr(self, "converged_"))
        n_init = self.n_init if do_init else 1

        max_lower_bound = -np.infty
        self.converged_ = False

        self.random_state = check_random_state(self.random_state)

        self.logL = np.empty((n_init, self.max_iter))
        self.logL[:] = np.nan

        # Initial evalution of the sufficient statistics for observables
        datas_visible = 0.0
        for traj in traj_list:
            datas_visible += sufficient_stats(traj, self.dim_x) / len(traj_list)

        self.coeffs_list_all = []

        # For referencement
        best_coeffs = None
        best_n_iter = -1
        best_n_init = -1

        for init in range(n_init):
            coeff_list_init = []
            if do_init:
                self._initialize_parameters(self.random_state, traj_len=_min_traj_len)
                if self.init_params == "markov":
                    self._m_step_markov(datas_visible)  # Initialize the visibles coefficients from markovian approx

            self._print_verbose_msg_init_beg(init)
            lower_bound = -np.infty if do_init else self.lower_bound_
            lower_bound_m_step = -np.infty
            # Algorithm loop
            for n_iter in range(1, self.max_iter + 1):
                prev_lower_bound = lower_bound
                new_stat = self._e_step_stats(traj_list, datas_visible)

                # new_stat_rs = self._rescale_hidden(new_stat)
                # new_stat_rs=new_stat

                lower_bound = self.loglikelihood(new_stat)
                if self.verbose >= 2:
                    if lower_bound - lower_bound_m_step < 0:
                        print("Delta ll after E step:", lower_bound - lower_bound_m_step)
                curr_coeffs = self.get_coefficients()
                curr_coeffs["ll"] = lower_bound
                coeff_list_init.append(curr_coeffs)

                self._m_step(new_stat)
                lower_bound_m_step = self.loglikelihood(new_stat)
                if self.verbose >= 2 and lower_bound_m_step - lower_bound < 0:
                    print("Delta ll after M step:", lower_bound_m_step - lower_bound)
                if np.isnan(lower_bound_m_step) or not self._check_finiteness():  # If we have nan value we simply restart the iteration
                    warnings.warn("Initialization %d has NaN values. Ends iteration" % (init), ConvergenceWarning)
                    if self.verbose >= 2:
                        print("Friction:\n{} \n Diffusion:\n{} \n Force :\n{} \n µ0 :\n{} \n Σ0:\n{} \n".format(self.friction_coeffs, self.diffusion_coeffs, self.force_coeffs, self.mu0, self.sig0))
                        print("ll: {}".format(lower_bound))
                    break

                self.logL[init, n_iter - 1] = lower_bound
                change = lower_bound - prev_lower_bound
                self._print_verbose_msg_iter_end(n_iter, change, lower_bound)

                if lower_bound > max_lower_bound:
                    max_lower_bound = lower_bound
                    best_coeffs = self.get_coefficients()
                    best_n_iter = n_iter
                    best_n_init = init

                if abs(change) < self.tol:
                    self.converged_ = True
                    if not self.no_stop:
                        break

            self._print_verbose_msg_init_end(lower_bound, n_iter)
            self.coeffs_list_all.append(coeff_list_init)
            if not self.converged_:
                warnings.warn("Initialization %d did not converge. " "Try different init parameters, " "or increase max_iter, tol " "or check for degenerate data." % (init + 1), ConvergenceWarning)
        if best_coeffs is not None:
            self.set_coefficients(best_coeffs, with_basis=False)  # Don't set basis coefficients
        self.n_iter_ = best_n_iter
        self.n_best_init_ = best_n_init
        self.lower_bound_ = max_lower_bound
        self._print_verbose_msg_fit_end(max_lower_bound, best_n_init, best_n_iter)

        return self

    def _e_step(self, traj):
        """E step.
        Parameters
        ----------
        traj : array-like, shape (n_timstep, dim_x) One trajectory

        Returns
        -------
        muh : array-like, shape (n_timstep, 2*dim_h)
            Mean values of the pair of the hidden variables
        Sigh : array-like, shape (n_timstep, 2*dim_h,2*dim_h)
            Covariances of the pair of the hidden variables
        """
        Xtplus, mutilde, R = self.model_class.compute_expectation_estep(traj, self.friction_coeffs, self.force_coeffs, self.dim_h, self.dt)

        return filtersmoother(Xtplus, mutilde, R, self.diffusion_coeffs, self.mu0, self.sig0)

    def _e_step_stats(self, traj_list, datas_visible):
        new_stat = 0.0
        if self.multiprocessing > 1:  # If we ask for more than one process
            with multiprocessing.Pool(processes=self.multiprocessing) as pool:
                proc = [pool.apply_async(e_step_worker_pool, args=(self, traj, datas_visible, len(traj_list))) for traj in traj_list]
                for p in proc:
                    ret = p.get()  # will block
                    new_stat += ret
        else:
            for traj in traj_list:
                muh, Sigh = self._e_step(traj)  # Compute hidden variable distribution
                new_stat += sufficient_stats_hidden(muh, Sigh, traj, datas_visible, self.dim_x, self.dim_h, self.dim_coeffs_force) / len(traj_list)
        return new_stat

    def _m_step(self, sufficient_stat):
        """M step.
        .. todo::   -Select dimension of fitted parameters from the sufficient stats (To deal with markovian initialization)
        """
        friction, force, diffusion = self.model_class.m_step(self.friction_coeffs, self.diffusion_coeffs, self.force_coeffs, sufficient_stat, self.dim_h, self.dt, self.OptimizeDiffusion, self.OptimizeForce)

        self.friction_coeffs = friction
        if self.OptimizeForce:
            self.force_coeffs = force
        self.mu0 = sufficient_stat["µ_0"]
        # self.sig0 = sufficient_stat["Σ_0"]
        # A, C = self.model_class._convert_local_coefficients(self.friction_coeffs, self.diffusion_coeffs, self.dt)
        # self.sig0 = C[self.dim_x :, self.dim_x :]
        if self.OptimizeDiffusion:
            self.diffusion_coeffs = diffusion

    def _em_step(self, vect_coeff, traj_list, datas_visible):
        """
        Wrapper of the E and M for gradient descent algorithm
        """
        self.unvectorization_coefficient(vect_coeff)
        new_stat = self._e_step_stats(traj_list, datas_visible)
        return self.loglikelihood(new_stat)

    def _check_finiteness(self):
        """
        Check that all quantities are finite
        """
        return np.isfinite(np.sum(self.friction_coeffs)) and np.isfinite(np.sum(self.diffusion_coeffs)) and np.isfinite(np.sum(self.force_coeffs)) and np.isfinite(np.sum(self.mu0)) and np.isfinite(np.sum(self.sig0))

    def _m_step_markov(self, sufficient_stat_vis):
        """Compute coefficients estimate via Markovian approximation to provide initialization"""
        A_full, C_full = self.model_class._convert_local_coefficients(self.friction_coeffs, self.diffusion_coeffs, self.dt)
        friction, force, diffusion = self.model_class.m_step(self.friction_coeffs, self.diffusion_coeffs, self.force_coeffs, sufficient_stat_vis, 0, self.dt, self.OptimizeDiffusion, self.OptimizeForce)

        A, C = self.model_class._convert_local_coefficients(friction, diffusion, self.dt)
        A_full[: self.dim_x, : self.dim_x] = A

        if self.OptimizeForce:
            self.force_coeffs = force
        if self.OptimizeDiffusion:
            C_full[: self.dim_x, : self.dim_x] = C
        (self.friction_coeffs, self.diffusion_coeffs) = self.model_class._convert_user_coefficients(A_full, C_full, self.dt)

    def loglikelihood(self, suff_datas, dim_h=None):
        """
        Return the current value of the negative log-likelihood
        """
        if dim_h is None:
            dim_h = self.dim_h
        ll = self.model_class.loglikelihood(suff_datas, self.friction_coeffs, self.diffusion_coeffs, self.force_coeffs, dim_h, self.dt)
        if dim_h > 0 and not np.isnan(suff_datas["hS"]):
            return ll + suff_datas["hS"]
        else:
            if dim_h > 0:
                warnings.warn("NaN value in hidden entropy")
            return ll

    def score(self, X, y=None, idx_trajs=[], Xh=None):
        """Compute the per-sample average log-likelihood of the given data X.

        Parameters
        ----------
        X : array-like, shape (n_samples, dim_x)
            List of positions data.
        idx_trajs: array, default []
            Location of split if multiple trajectory are inputed

        Returns
        -------
        log_likelihood : float
            Log likelihood of the generalized Langevin model given X.
        """
        check_is_fitted(self, "initialized_")
        X = check_array(X, ensure_min_samples=4, ensure_min_features=self.dim_x)
        self.dt = X[1, 0] - X[0, 0]

        Xproc, idx_trajs = self.model_class.preprocessingTraj(self.basis, X, idx_trajs=idx_trajs)
        traj_list = np.split(Xproc, idx_trajs)
        # Initial evalution of the sufficient statistics for observables
        new_stat = 0.0
        if Xh is None:
            datas_visible = 0.0
            for traj in traj_list:
                datas_visible += sufficient_stats(traj, self.dim_x) / len(traj_list)
            new_stat = self._e_step_stats(traj_list, datas_visible)
            # muh, Sigh = self._e_step(traj)  # Compute hidden variable distribution
            # new_stat += sufficient_stats_hidden(muh, Sigh, traj, datas, self.dim_x, self.dim_h, self.dim_coeffs_force) / len(traj_list)
        else:
            traj_list_h = np.split(Xh, idx_trajs)
            for n, traj in enumerate(traj_list):
                datas_visible = sufficient_stats(traj, self.dim_x)
                zero_sig = np.zeros((len(traj), 2 * self.dim_h, 2 * self.dim_h))
                muh = np.hstack((np.roll(traj_list_h[n], -1, axis=0), traj_list_h[n]))
                new_stat += sufficient_stats_hidden(muh, zero_sig, traj, datas_visible, self.dim_x, self.dim_h, self.dim_coeffs_force) / len(traj_list)
        lower_bound = self.loglikelihood(new_stat)
        return lower_bound

    def predict(self, X, idx_trajs=[]):
        """Predict the hidden variables for the data samples in X using trained model.

        Parameters
        ----------
        X : array-like, shape (n_samples, dim_x)
            List of positions data.
        idx_trajs: array, default []
            Location of split if multiple trajectory are inputed

        Returns
        -------
        muh : array, shape (n_samples,dim_h)
            Average value of the hidden variables given current estimation of the parameters and provided trajectories.
        """
        check_is_fitted(self, "converged_")
        X = check_array(X, ensure_min_samples=4, ensure_min_features=self.dim_x)
        Xproc, idx_trajs = self.model_class.preprocessingTraj(self.basis, X, idx_trajs=idx_trajs)
        traj_list = np.split(Xproc, idx_trajs)
        muh_out = None
        for traj in traj_list:
            muh, Sigh = self._e_step(traj)  # Compute hidden variable distribution
            if muh_out is None:
                muh_out = muh[:, self.dim_h :]
            else:
                muh_out = np.hstack((muh_out, muh[:, self.dim_h :]))
        return muh_out

    def sample(self, n_samples=50, n_trajs=1, x0=None, v0=None, dt=5e-3, burnout=0, rng=None):
        """Generate random samples from the fitted GLE model.

        Use the provided basis to compute the force term. The basis should be fitted first, if not will be fitted using a dummy trajectory.

        Parameters
        ----------
        n_samples : int, default=50
            Number of timestep per trajectory to generate. Defaults to 50.
        n_trajs : int, default=1
            Number of trajectory to generate
        x0,v0 : array-like, optionnal
            Initial value of the trajectory
        dt : float, default = 5e-3
            Timestep to use for the integration
        burnout : int, default to 0
            Remove [burnout] step from the start of the trajectories
        rng: random number generator
            Provide if wanted a random number generator. See numpy documention about RNG

        Returns
        -------
        X : array, shape (n_samples, n_features)
            Randomly generated trajectory
        idx : array
            Index of the trajectories
        y : array, shape (nsamples,)
            Hidden variables values
        """
        if rng is None:
            self.random_state = check_random_state(self.random_state)
        else:
            self.random_state = check_random_state(rng)
        self.dt = dt
        self._check_initial_parameters()
        if not (self.warm_start or hasattr(self, "converged_")):
            self._initialize_parameters(self.random_state, traj_len=n_samples)

        if x0 is None:
            x0 = np.zeros((self.dim_x))
        if v0 is None:
            v0 = np.zeros((self.dim_x))

        X = None
        idx_trajs = []
        X_h = None

        if self.multiprocessing > 1:  # If we ask for more than one process
            child_seeds = np.random.SeedSequence(self.random_state.get_state()[1]).spawn(n_trajs)
            with multiprocessing.Pool(processes=self.multiprocessing) as pool:
                proc = [
                    pool.apply_async(
                        self.model_class.generator,
                        args=(),
                        kwds={
                            "nsteps": n_samples,
                            "rng": np.random.default_rng(child_seeds[n]),
                            "dt": self.dt,
                            "dim_h": self.dim_h,
                            "x0": x0,
                            "v0": v0,
                            "friction": self.friction_coeffs,
                            "SST": self.diffusion_coeffs,
                            "force_coeffs": self.force_coeffs,
                            "muh0": self.mu0,
                            "sigh0": self.sig0,
                            "basis": self.basis,
                        },
                    )
                    for n in range(n_trajs)
                ]
                for p in proc:
                    txv, h = p.get()  # will block
                    if X is None:
                        X = txv[burnout:, :]
                    else:
                        idx_trajs.append(len(X))
                        X = np.vstack((X, txv[burnout:, :]))

                    if X_h is None:
                        X_h = h[burnout:, :]
                    else:
                        X_h = np.vstack((X_h, h[burnout:, :]))
        else:
            for n in range(n_trajs):
                txv, h = self.model_class.generator(nsteps=n_samples, dt=self.dt, dim_h=self.dim_h, x0=x0, v0=v0, friction=self.friction_coeffs, SST=self.diffusion_coeffs, force_coeffs=self.force_coeffs, muh0=self.mu0, sigh0=self.sig0, basis=self.basis, rng=self.random_state)

                if X is None:
                    X = txv[burnout:, :]
                else:
                    idx_trajs.append(len(X))
                    X = np.vstack((X, txv[burnout:, :]))

                if X_h is None:
                    X_h = h[burnout:, :]
                else:
                    X_h = np.vstack((X_h, h[burnout:, :]))
        return X, idx_trajs, X_h

    def get_coefficients(self):
        """Return the actual values of the fitted coefficients."""
        A, C = self.model_class._convert_local_coefficients(self.friction_coeffs, self.diffusion_coeffs, self.dt)
        return {"A": A, "C": C, "force": self.force_coeffs, "µ_0": self.mu0, "Σ_0": self.sig0, "SST": self.diffusion_coeffs, "dt": self.dt, "basis": self.basis.get_coefficients()}

    def set_coefficients(self, coeffs, with_basis=True):
        """Set the value of the coefficients

        Parameters
        ----------
        coeffs : dict
            Contains the value of the coefficients to set.
        """
        (self.friction_coeffs, self.diffusion_coeffs) = self.model_class._convert_user_coefficients(np.asarray(coeffs["A"]), np.asarray(coeffs["C"]), self.dt)
        (self.force_coeffs, self.mu0, self.sig0) = (coeffs["force"], coeffs["µ_0"], coeffs["Σ_0"])
        if with_basis:
            self.basis.set_coefficients(coeffs["basis"])

    def _n_parameters(self):
        """Return the number of free parameters in the model."""
        return 2 * (self.dim_x + self.dim_h) ** 2 + self.dim_h + self.dim_h ** 2

    def bic(self, X):
        """Bayesian information criterion for the current model on the input X.

        Parameters
        ----------
        X : array of shape (n_samples, n_features)

        Returns
        -------
        bic : float
            The lower the better.
        """
        return -2 * self.score(X) * X.shape[0] + self._n_parameters() * np.log(X.shape[0])

    def aic(self, X):
        """Akaike information criterion for the current model on the input X.

        Parameters
        ----------
        X : array of shape (n_samples, n_features)

        Returns
        -------
        aic : float
            The lower the better.
        """
        return -2 * self.score(X) * X.shape[0] + 2 * self._n_parameters()

    def _print_verbose_msg_init_beg(self, n_init):
        """Print verbose message on initialization."""
        if self.verbose == 1:
            print("Initialization %d" % n_init)
        elif self.verbose >= 2:
            print("Initialization %d" % n_init)
            self._init_prev_time = time()
            self._iter_prev_time = self._init_prev_time
        if self.verbose >= 3:
            print("----------------Current parameters values------------------")
            print(self.get_coefficients())

    def _print_verbose_msg_iter_end(self, n_iter, diff_ll, log_likelihood):
        """Print verbose message on initialization."""
        if n_iter % self.verbose_interval == 0:
            if self.verbose == 1:
                print("***Iteration EM*** : {} / {} --- Current loglikelihood {}".format(n_iter, self.max_iter, log_likelihood))
            elif self.verbose >= 2:
                cur_time = time()
                print("***Iteration EM*** :%d / %d\t time lapse %.5fs\t Current loglikelihood %.5f loglikelihood change %.5f" % (n_iter, self.max_iter, cur_time - self._iter_prev_time, log_likelihood, diff_ll))
                self._iter_prev_time = cur_time
            if self.verbose >= 3:
                print("----------------Current parameters values------------------")
                print(self.get_coefficients())

    def _print_verbose_msg_init_end(self, ll, best_iter):
        """Print verbose message on the end of iteration."""
        if self.verbose == 1:
            print("Initialization converged: %s at step %i \t ll %.5f" % (self.converged_, best_iter, ll))
        elif self.verbose >= 2:
            print("Initialization converged: %s at step %i \t time lapse %.5fs\t ll %.5f" % (self.converged_, best_iter, time() - self._init_prev_time, ll))
            print("----------------Current parameters values------------------")
            print(self.get_coefficients())

    def _print_verbose_msg_fit_end(self, ll, best_init, best_iter):
        """Print verbose message on the end of iteration."""
        if self.verbose == 1:
            print("Fit converged: %s Init: %s at step %i \t ll %.5f" % (self.converged_, best_init, best_iter, ll))
        elif self.verbose >= 2:
            print("Fit converged: %s Init: %s at step %i \t time lapse %.5fs\t ll %.5f" % (self.converged_, best_init, best_iter, time() - self._init_prev_time, ll))
            print("----------------Fitted parameters values------------------")
            print(self.get_coefficients())


text:
#!/usr/bin/env python

import argparse
import json
import os
import random

import numpy as np

import ray
from ray.tune import Trainable, run, sample_from
from ray.tune.schedulers import HyperBandScheduler

class MyTrainableClass(Trainable):
    """Example agent whose learning curve is a random sigmoid.

    The dummy hyperparameters "width" and "height" determine the slope and
    maximum reward value reached.
    """

    def _setup(self, config):
        self.timestep = 0

    def _train(self):
        self.timestep += 1
        v = np.tanh(float(self.timestep) / self.config.get("width", 1))
        v *= self.config.get("height", 1)

        # Here we use `episode_reward_mean`, but you can also report other
        # objectives such as loss or accuracy.
        return {"episode_reward_mean": v}

    def _save(self, checkpoint_dir):
        path = os.path.join(checkpoint_dir, "checkpoint")
        with open(path, "w") as f:
            f.write(json.dumps({"timestep": self.timestep}))
        return path

    def _restore(self, checkpoint_path):
        with open(checkpoint_path) as f:
            self.timestep = json.loads(f.read())["timestep"]

if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument(
        "--smoke-test", action="store_true", help="Finish quickly for testing")
    args, _ = parser.parse_known_args()
    ray.init(num_cpus=4 if args.smoke_test else None)

    # Hyperband early stopping, configured with `episode_reward_mean` as the
    # objective and `training_iteration` as the time unit,
    # which is automatically filled by Tune.
    hyperband = HyperBandScheduler(
        time_attr="training_iteration",
        metric="episode_reward_mean",
        mode="max",
        max_t=200)

    run(MyTrainableClass,
        name="hyperband_test",
        num_samples=20,
        stop={"training_iteration": 1 if args.smoke_test else 99999},
        config={
            "width": sample_from(lambda spec: 10 + int(90 * random.random())),
            "height": sample_from(lambda spec: int(100 * random.random()))
        },
        scheduler=hyperband,
        fail_fast=True)

max_stars_repo_name: acmore/ray

max_issues_repo_name: acmore/ray

max_forks_repo_name: acmore/ray


text:
#!/usr/bin/env python
from __future__ import print_function, division

import os
import sys
sys.path.append(os.path.dirname(sys.path[0]))
import argparse
from random import randint

from matplotlib import pyplot as plt
import numpy as np
from scipy.interpolate import RegularGridInterpolator

from cryoio import mrc
import density
import geometry

def gen_slices(model_files, fspace=False, log_scale=True):
    for model in model_files:
        M = mrc.readMRC(model)
        N = M.shape[0]
        print('model size: {0}x{0}x{0}'.format(N))

        oversampling_factor = 3
        zeropad = oversampling_factor - 1  # oversampling factor = zeropad + 1
        psize = 3.0 * oversampling_factor
        beamstop_freq = 0.003
        mask = geometry.gen_dense_beamstop_mask(N, 2, beamstop_freq, psize=psize)
        # mask = None

        if fspace:
            fM = M
        else:
            M_totalmass = 1000000
            M *= M_totalmass / M.sum()

            V = density.real_to_fspace_with_oversampling(M, oversampling_factor)
            fM = V.real ** 2 + V.imag ** 2

            mask_3D = geometry.gen_dense_beamstop_mask(N, 3, beamstop_freq, psize=psize)
            fM *= mask_3D
            mrc.writeMRC('particle/{}_fM_totalmass_{}_oversampling_{}.mrc'.format(
                os.path.splitext(os.path.basename(model))[0], str(int(M_totalmass)).zfill(5), oversampling_factor
                ), fM, psz=psize)

        slicing_func = RegularGridInterpolator([np.arange(N),]*3, fM, bounds_error=False, fill_value=0.0)
        coords = geometry.gencoords_base(N, 2)

        fig, axes = plt.subplots(3, 3, figsize=(12.9, 9.6))
        for i, ax in enumerate(axes.flat):
            row, col = np.unravel_index(i, (3, 3))

            # Randomly generate the viewing direction/shift
            pt = np.random.randn(3)
            pt /= np.linalg.norm(pt)
            psi = 2 * np.pi * np.random.rand()
            EA = geometry.genEA(pt)[0]
            EA[2] = psi
            R = geometry.rotmat3D_EA(*EA)[:, 0:2]
            rotated_coords = R.dot(coords.T).T + int(N/2)
            img = slicing_func(rotated_coords).reshape(N, N)
            img = np.require(np.random.poisson(img), dtype=np.float32)

            if log_scale:
                img = np.log(np.maximum(img, 0))
            if mask is not None:
                img *= mask

            im = ax.imshow(img, origin='lower')  # cmap='Greys'
            ticks = [0, int(N/4.0), int(N/2.0), int(N*3.0/4.0), int(N-1)]
            if row == 2:
                ax.set_xticks(ticks)
            else:
                ax.set_xticks([])
            if col == 0:
                ax.set_yticks(ticks)
            else:
                ax.set_yticks([])
            fig.colorbar(im, ax=ax)

        # fig.subplots_adjust(right=0.8)
        # cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
        # fig.colorbar(im, cax=cbar_ax)
        fig.suptitle('simulated experimental data of XFEL for {}'.format(model))
        # fig.tight_layout()
    plt.show()

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument("model_files", help="list of mrcs files.", nargs='+')
    parser.add_argument("-l", "--log_scale", help="show image in log scale.",
                        action="store_true")
    parser.add_argument('-f', '--fspace', help='input model has been already in Fourier space.', action='store_true')
    args = parser.parse_args()

    log_scale = args.log_scale
    model_files = args.model_files
    fspace = args.fspace
    print('model_files:', model_files)
    print('fspace:', fspace)
    print('log_scale', log_scale)
    gen_slices(model_files, fspace=fspace, log_scale=log_scale)

max_stars_repo_name: lqhuang/SOD-xfel

max_issues_repo_name: lqhuang/SOD-xfel

max_forks_repo_name: lqhuang/SOD-xfel


text:
# -*- coding: utf-8 -*-
"""
Jeu de la bataille

@author: JH218595
"""

import matplotlib.pyplot as plt
import numpy as np

def new_card_set():
    '''
    Creates an returns a brand new card set.
    '''
    cards = []
    cards_nb = list(range(9,15)) #
    cards_sym= list(['P', 'C', 'Q', 'T'])
    for nb in cards_nb:
        for sym in cards_sym:
            cards.append((nb, sym))
    return cards

def sort_cards(cards):
        '''
        Creates two card sets A and B from a new card set
        '''
        from numpy.random import randint
        A = []
        B = []
        while len(cards) != 0: # tant que le paquet n'est pas vide on distribue
            A.append(cards.pop(randint(low=0, high=len(cards))))
            B.append(cards.pop(randint(low=0, high=len(cards))))
        return A, B

def play_bataille(A,B):
    '''
    Joue à la bataille et retourne l'évolution du nombre de cartes par joueur
    au cours d'une partie.
    '''
    manche_A = []
    manche_B = []

    taille_A = []
    taille_B = []

    continuer = True

    while continuer:
        taille_A.append(len(A))
        taille_B.append(len(B))

#        if len(taille_B)>100: # trop de manche, on arrête la partie
#            continuer = False
        if len(A) == 0:
            #print('### B gagne la partie ###')
            continuer = False
        elif len(B) == 0:
            #print('### A gagne la partie ###')
            continuer = False
        else:
            manche_A.append(A.pop())
            manche_B.append(B.pop())

            if manche_A[-1][0] == manche_B[-1][0]:
                #print('Bataille !!')
                pass
            elif manche_A[-1][0] > manche_B[-1][0]:
                #print('A remporte la manche')
                A = manche_A + manche_B + A
                manche_A.clear()
                manche_B.clear()
            else:
                #print('B remporte la manche')
                B = manche_B + manche_A + B
                manche_A.clear()
                manche_B.clear()

    return taille_A, taille_B

cards = new_card_set()
A, B = sort_cards(cards)
taille_A, taille_B = play_bataille(A, B)

plt.figure(1)
plt.clf()
plt.plot(taille_A, lw=2)
plt.plot(taille_B, lw=2)
plt.xlabel('Nombre de manches')
plt.grid(True)

# fait un nombre N de partie et affiche les statistiques
def play_N_batailles(N):
    nb_manches = []
    for idx in range(N):
        cards = new_card_set()
        A, B = sort_cards(cards)
        taille_A, taille_B = play_bataille(A, B)
        nb_manches.append(len(taille_A))
    return nb_manches

N = 10000
nb_manches = np.array(play_N_batailles(N))

plt.figure(2)
plt.hist(nb_manches, bins=500, color='b', alpha=0.5)

from scipy.stats import gamma

params = gamma.fit(nb_manches)

x = np.arange(1, N)
#plt.figure(3)
plt.plot(x, 2*N*gamma.pdf(x, *params[:-2], loc=params[-2], scale=params[-1]),
                          lw=2, color='r')
plt.xlim(0, 300)

max_stars_repo_name: Hash--/documents

max_issues_repo_name: Hash--/documents

max_forks_repo_name: Hash--/documents

reason: import numpy,from numpy,from scipy


text:
import numpy as np

def _summarize_peaks(peaks):
    """
    merge peaks position if closer than 10
    """
    previous = peaks[0]
    new_peaks = [previous]
    for pos in peaks:
        if pos > previous + 10:
            new_peaks.add(pos)
        previous = pos
    return new_peaks

def find_mature(x, y, win=10):
    """
    Window apprach to find hills in the expression profile
    """
    previous = min(y)
    peaks = []
    intervals = range(x, y, win)
    for pos in intervals:
        if y[pos] > previous * 10:
            previous = y[pos]
            peaks.add(pos)
    peaks = _summarize_peaks(peaks)

max_stars_repo_name: kkarolis/seqcluster

max_issues_repo_name: kkarolis/seqcluster

max_forks_repo_name: kkarolis/seqcluster