text:
#   Copyright (c) 2019 PaddlePaddle Authors. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""Mask, padding and batching."""
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
import numpy as np

def mask(batch_tokens, total_token_num, vocab_size, CLS=1, SEP=2, MASK=3):
    """
    Add mask for batch_tokens, return out, mask_label, mask_pos;
    Note: mask_pos responding the batch_tokens after padded;
    """
    max_len = max([len(sent) for sent in batch_tokens])
    mask_label = []
    mask_pos = []
    prob_mask = np.random.rand(total_token_num)
    # Note: the first token is [CLS], so [low=1]
    replace_ids = np.random.randint(1, high=vocab_size, size=total_token_num)
    pre_sent_len = 0
    prob_index = 0
    for sent_index, sent in enumerate(batch_tokens):
        mask_flag = False
        prob_index += pre_sent_len
        for token_index, token in enumerate(sent):
            prob = prob_mask[prob_index + token_index]
            if prob > 0.15:
                continue
            elif 0.03 < prob <= 0.15:
                # mask
                if token != SEP and token != CLS:
                    mask_label.append(sent[token_index])
                    sent[token_index] = MASK
                    mask_flag = True
                    mask_pos.append(sent_index * max_len + token_index)
            elif 0.015 < prob <= 0.03:
                # random replace
                if token != SEP and token != CLS:
                    mask_label.append(sent[token_index])
                    sent[token_index] = replace_ids[prob_index + token_index]
                    mask_flag = True
                    mask_pos.append(sent_index * max_len + token_index)
            else:
                # keep the original token
                if token != SEP and token != CLS:
                    mask_label.append(sent[token_index])
                    mask_pos.append(sent_index * max_len + token_index)
        pre_sent_len = len(sent)
        # ensure at least mask one word in a sentence
        while not mask_flag:
            token_index = int(np.random.randint(1, high=len(sent) - 1, size=1))
            if sent[token_index] != SEP and sent[token_index] != CLS:
                mask_label.append(sent[token_index])
                sent[token_index] = MASK
                mask_flag = True
                mask_pos.append(sent_index * max_len + token_index)
    mask_label = np.array(mask_label).astype("int64").reshape([-1, 1])
    mask_pos = np.array(mask_pos).astype("int64").reshape([-1, 1])
    return batch_tokens, mask_label, mask_pos

def prepare_batch_data(insts,
                       total_token_num,
                       max_len=None,
                       voc_size=0,
                       pad_id=None,
                       cls_id=None,
                       sep_id=None,
                       mask_id=None,
                       return_input_mask=True,
                       return_max_len=True,
                       return_num_token=False):
    """
    1. generate Tensor of data
    2. generate Tensor of position
    3. generate self attention mask, [shape: batch_size *  max_len * max_len]
    """
    batch_src_ids = [inst[0] for inst in insts]
    batch_sent_ids = [inst[1] for inst in insts]
    batch_pos_ids = [inst[2] for inst in insts]
    labels_list = []
    # compatible with mrqa, whose example includes start/end positions,
    # or unique id
    for i in range(3, len(insts[0]), 1):
        labels = [inst[i] for inst in insts]
        labels = np.array(labels).astype("int64").reshape([-1, 1])
        labels_list.append(labels)
    # First step: do mask without padding
    if mask_id >= 0:
        out, mask_label, mask_pos = mask(
            batch_src_ids,
            total_token_num,
            vocab_size=voc_size,
            CLS=cls_id,
            SEP=sep_id,
            MASK=mask_id)
    else:
        out = batch_src_ids
    # Second step: padding
    src_id, self_input_mask = pad_batch_data(
        out,
        max_len=max_len,
        pad_idx=pad_id, return_input_mask=True)
    pos_id = pad_batch_data(
        batch_pos_ids,
        max_len=max_len,
        pad_idx=pad_id,
        return_pos=False,
        return_input_mask=False)
    sent_id = pad_batch_data(
        batch_sent_ids,
        max_len=max_len,
        pad_idx=pad_id,
        return_pos=False,
        return_input_mask=False)
    if mask_id >= 0:
        return_list = [
            src_id, pos_id, sent_id, self_input_mask, mask_label, mask_pos
        ] + labels_list
    else:
        return_list = [src_id, pos_id, sent_id, self_input_mask] + labels_list
    return return_list if len(return_list) > 1 else return_list[0]

def pad_batch_data(insts,
                   max_len=None,
                   pad_idx=0,
                   return_pos=False,
                   return_input_mask=False,
                   return_max_len=False,
                   return_num_token=False):
    """
    Pad the instances to the max sequence length in batch, and generate the
    corresponding position data and input mask.
    """
    return_list = []
    if max_len is None:
        max_len = max(len(inst) for inst in insts)
    # Any token included in dict can be used to pad, since the paddings' loss
    # will be masked out by weights and make no effect on parameter gradients.
    inst_data = np.array([
        list(inst) + list([pad_idx] * (max_len - len(inst))) for inst in insts
    ])
    return_list += [inst_data.astype("int64").reshape([-1, max_len, 1])]
    # position data
    if return_pos:
        inst_pos = np.array([
            list(range(0, len(inst))) + [pad_idx] * (max_len - len(inst))
            for inst in insts
        ])
        return_list += [inst_pos.astype("int64").reshape([-1, max_len, 1])]
    if return_input_mask:
        # This is used to avoid attention on paddings.
        input_mask_data = np.array([[1] * len(inst) + [0] *
                                    (max_len - len(inst)) for inst in insts])
        input_mask_data = np.expand_dims(input_mask_data, axis=-1)
        return_list += [input_mask_data.astype("float32")]
    if return_max_len:
        return_list += [max_len]
    if return_num_token:
        num_token = 0
        for inst in insts:
            num_token += len(inst)
        return_list += [num_token]
    return return_list if len(return_list) > 1 else return_list[0]

if __name__ == "__main__":
    pass


text:
import torch
import pytorch_lightning as pl
from pytorch_lightning.logging import TestTubeLogger
from pytorch_lightning.callbacks import ModelCheckpoint

import sys
sys.path.append('..')

import argparse

from networks import RNNDecoder, RNNEncoder
from torch.utils.data import DataLoader
from datasets import MolecularDataset

from rdkit import Chem
import numpy as np

from rdkit import RDLogger
lg = RDLogger.logger()
lg.setLevel(RDLogger.CRITICAL)

class RNN_VAE(pl.LightningModule):
    def __init__(self):
        super(RNN_VAE, self).__init__()
        self.z_dim = 44

        self.enc = RNNEncoder(out_dim=2 * self.z_dim)
        self.dec = RNNDecoder(in_dim=self.z_dim)

        self.beta = 0.01

    @staticmethod
    def sample_repar_z(means, logvar):
        return means + torch.randn_like(means) * torch.exp(0.5 * logvar)

    @staticmethod
    def kl_div(means_q, logvar_q, means_p=None, logvar_p=None):
        if means_p is None:  # prior is N(0, I)
            return -0.5 * torch.mean(
                torch.sum(1 + logvar_q - means_q.pow(2) - logvar_q.exp(),
                          dim=-1))
        else:
            return -0.5 * torch.mean(
                torch.sum(1 - logvar_p + logvar_q -
                          (means_q.pow(2) + logvar_q.exp()) * (-logvar_p).exp(),
                          dim=-1))

    def training_step(self, batch, batch_nb):
        # pair of objects
        x, _ = batch

        # compute proposal distributions
        p_z_means, p_z_logvar = torch.split(self.enc(x), self.z_dim, -1)

        # sample z
        z_sample = self.sample_repar_z(p_z_means, p_z_logvar)

        # compute kl divergence
        kl = self.kl_div(p_z_means, p_z_logvar)

        # compute reconstrunction loss
        x_by_z_logprob = self.dec.get_log_prob(x=x, z=z_sample).mean()

        loss = (-x_by_z_logprob + self.beta * kl)

        return {'loss': loss,
                'log': {
                    'x_by_z_logprob': x_by_z_logprob,
                    'kl': kl}
                }

    def validation_step(self, batch, batch_nb):
        x, _ = batch

        # compute proposal distributions
        p_z_means, p_z_logvar = torch.split(self.enc(x), self.z_dim, -1)

        # sample z
        z_sample = self.sample_repar_z(p_z_means, p_z_logvar)

        sampled_sm = self.dec.sample(z_sample)

        valid_proc = len([s for s in sampled_sm if Chem.MolFromSmiles(s) is not None]) / len(sampled_sm)
        unique_proc = len(np.unique(sampled_sm)) / len(sampled_sm)
        eq_proc = len([s  for (x_ob, s) in zip(x, sampled_sm) if (Chem.MolFromSmiles(s) is not None) and (Chem.MolToSmiles(Chem.MolFromSmiles(s))) == x_ob]) / len(sampled_sm)

        return {'valid': valid_proc, 'unique': unique_proc, 'equal': eq_proc}

    def validation_end(self, outputs):
        val_stats = {}

        val_stats['val_valid'] = np.array([x['valid'] for x in outputs]).mean()
        val_stats['val_unique'] = np.array([x['unique'] for x in outputs]).mean()
        val_stats['val_equal'] = np.array([x['equal'] for x in outputs]).mean()

        return {'log': val_stats}

    def configure_optimizers(self):
        return torch.optim.Adam(self.parameters(), lr=3e-4)

    @pl.data_loader
    def train_dataloader(self):
        return DataLoader(MolecularDataset('../data/train.csv', train=True),
                          batch_size=512, shuffle=True, num_workers=4)

    @pl.data_loader
    def val_dataloader(self):
        return DataLoader(MolecularDataset('../data/train.csv', train=False),
                          batch_size=512, shuffle=True, num_workers=4)

if __name__ == '__main__':
    torch.manual_seed(777)
    torch.backends.cudnn.deterministic = True
    torch.backends.cudnn.benchmark = False

    parser = argparse.ArgumentParser(description=
                                     'Script to pretrain RNN encoder and decoder')
    parser.add_argument('--gpu', type=int, default=-1)

    args = parser.parse_args()

    model = RNN_VAE()

    model_checkpoint = ModelCheckpoint('../saved_models/rnn_vae',
                                    save_best_only=False, period=10)

    logger = TestTubeLogger(save_dir='../logs', name='rnn_vae')

    tr = pl.Trainer(gpus=([] if (args.gpu < 0) else [args.gpu]),
                    logger=logger,
                    early_stop_callback=False,
                    show_progress_bar=True,
                    max_nb_epochs=200,
                    default_save_path=None,
                    checkpoint_callback=model_checkpoint)
    tr.fit(model)
    torch.save(model.enc.state_dict(), '../saved_models/rnn_enc.ckpt')
    torch.save(model.dec.state_dict(), '../saved_models/rnn_dec.ckpt')

ext: py

max_stars_repo_name: magelead/BiAAE

max_issues_repo_name: magelead/BiAAE

max_forks_repo_name: magelead/BiAAE

reason: import numpy


text:
import torch
import torchvision.utils as vutils
import glob, os
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import numpy as np
import argparse

parser = argparse.ArgumentParser()
parser.add_argument("--save_image", type=bool, default=False)
parser.add_argument("--save_video", type=bool, default=False)

def get_files():
    files = glob.glob("*.pt")
    files = [(file, int(file.split('.')[0].split('-')[-1])) for file in files]
    files.sort(key=lambda f: f[1])
    files = [f[0] for f in files]
    return files

if __name__ == "__main__":
    args = parser.parse_args()
    img_list = []

    files = get_files()
    for file in files:
        module = torch.jit.load(file, map_location=torch.device("cpu"))
        images = list(module.parameters())[0].detach().cpu().transpose(2, 3)
        img_list.append(vutils.make_grid(images, padding=2, normalize=True))

        if args.save_image:
            # save fake images to directory
            directory = os.path.splitext(file)[0]
            if not os.path.exists(directory):
                os.mkdir(directory)
            images = [images[i].reshape(3, 64, 64) for i in range(64)]
            for i, image in enumerate(images):
                vutils.save_image(image,
                                  directory + '/' + str(i) + '.png',
                                  normalize=True)

    if args.save_video:
        fig = plt.figure(figsize=(8, 8))
        plt.axis("off")
        ims = [[plt.imshow(np.transpose(i, (1, 2, 0)), animated=True)]
               for i in img_list]
        ani = animation.ArtistAnimation(fig,
                                        ims,
                                        interval=1000,
                                        repeat_delay=1000,
                                        blit=True)

        writer = animation.writers['ffmpeg']
        writer = writer(fps=1, metadata=dict(artist='Me'), bitrate=1800)
        ani.save("dcgan.mp4", writer)

max_stars_repo_name: sneaxiy/gotorch

max_issues_repo_name: sneaxiy/gotorch

max_forks_repo_name: sneaxiy/gotorch


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

import numpy as np
import tensorflow as tf
import tensorflow.keras.backend as K
from tensorflow.keras.initializers import Constant
from tensorflow.keras.layers import InputSpec, Layer, Lambda, Dropout, Multiply

INFTY = -100

class Highway(Layer):

    def __init__(self, activation=None, bias_initializer=-1, **kwargs):
        super().__init__(**kwargs)
        self.activation = tf.keras.activations.get(activation)
        self.bias_initializer = bias_initializer
        if isinstance(self.bias_initializer, int):
            self.bias_initializer = Constant(self.bias_initializer)
        self.input_spec = [InputSpec(min_ndim=2)]

    def build(self, input_shape):
        assert len(input_shape) >= 2
        input_dim = input_shape[-1]

        self.gate_kernel = self.add_weight(
            shape=(input_dim, input_dim), initializer='uniform', name='gate_kernel')
        self.gate_bias = self.add_weight(
            shape=(input_dim,), initializer=self.bias_initializer, name='gate_bias')
        self.dense_kernel = self.add_weight(
            shape=(input_dim, input_dim), initializer='uniform', name='dense_kernel')
        self.dense_bias = self.add_weight(
            shape=(input_dim,), initializer=self.bias_initializer, name='dense_bias')
        self.input_spec = InputSpec(min_ndim=2, axes={-1: input_dim})
        self.built = True

    def call(self, inputs, **kwargs):
        gate = K.dot(inputs, self.gate_kernel)
        gate = K.bias_add(gate, self.gate_bias, data_format="channels_last")
        gate = self.activation(gate)
        new_value = K.dot(inputs, self.dense_kernel)
        new_value = K.bias_add(new_value, self.dense_bias, data_format="channels_last")
        return gate * new_value + (1.0 - gate) * inputs

    def compute_output_shape(self, input_shape):
        return input_shape

def weighted_sum(first, second, sigma, first_threshold=-np.inf, second_threshold=np.inf):
    logit_probs = first * sigma + second * (1.0 - sigma)
    infty_tensor = K.ones_like(logit_probs) * INFTY
    logit_probs = K.switch(K.greater(first, first_threshold), logit_probs, infty_tensor)
    logit_probs = K.switch(K.greater(second, second_threshold), logit_probs, infty_tensor)
    return logit_probs

class WeightedCombinationLayer(Layer):

    """
    A class for weighted combination of probability distributions
    """

    def __init__(self, first_threshold=None, second_threshold=None,
                 use_dimension_bias=False, use_intermediate_layer=False,
                 intermediate_dim=64, intermediate_activation=None,
                 from_logits=False, return_logits=False,
                 bias_initializer=1.0, **kwargs):
        # if 'input_shape' not in kwargs:
        #     kwargs['input_shape'] = [(None, input_dim,), (None, input_dim)]
        super(WeightedCombinationLayer, self).__init__(**kwargs)
        self.first_threshold = first_threshold if first_threshold is not None else INFTY
        self.second_threshold = second_threshold if second_threshold is not None else INFTY
        self.use_dimension_bias = use_dimension_bias
        self.use_intermediate_layer = use_intermediate_layer
        self.intermediate_dim = intermediate_dim
        self.intermediate_activation = tf.keras.activations.get(intermediate_activation)
        self.from_logits = from_logits
        self.return_logits = return_logits
        self.bias_initializer = bias_initializer
        self.input_spec = [InputSpec(), InputSpec(), InputSpec()]

    def build(self, input_shape):
        assert len(input_shape) == 3
        assert input_shape[0] == input_shape[1]
        assert input_shape[0][:-1] == input_shape[2][:-1]

        input_dim, features_dim = input_shape[0][-1], input_shape[2][-1]
        if self.use_intermediate_layer:
            self.first_kernel = self.add_weight(
                shape=(features_dim, self.intermediate_dim),
                initializer="random_uniform", name='first_kernel')
            self.first_bias = self.add_weight(
                shape=(self.intermediate_dim,),
                initializer="random_uniform", name='first_bias')
        self.features_kernel = self.add_weight(
            shape=(features_dim, 1), initializer="random_uniform", name='kernel')
        self.features_bias = self.add_weight(
            shape=(1,), initializer=Constant(self.bias_initializer), name='bias')
        if self.use_dimension_bias:
            self.dimensions_bias = self.add_weight(
                shape=(input_dim,), initializer="random_uniform", name='dimension_bias')
        super(WeightedCombinationLayer, self).build(input_shape)

    def call(self, inputs, **kwargs):
        assert isinstance(inputs, list) and len(inputs) == 3
        first, second, features = inputs[0], inputs[1], inputs[2]
        if not self.from_logits:
            first = K.clip(first, 1e-10, 1.0)
            second = K.clip(second, 1e-10, 1.0)
            first_, second_ = K.log(first), K.log(second)
        else:
            first_, second_ = first, second
        # embedded_features.shape = (M, T, 1)
        if self.use_intermediate_layer:
            features = K.dot(features, self.first_kernel)
            features = K.bias_add(features, self.first_bias, data_format="channels_last")
            features = self.intermediate_activation(features)
        embedded_features = K.dot(features, self.features_kernel)
        embedded_features = K.bias_add(
            embedded_features, self.features_bias, data_format="channels_last")
        if self.use_dimension_bias:
            tiling_shape = [1] * (K.ndim(first) - 1) + [K.shape(first)[-1]]
            embedded_features = K.tile(embedded_features, tiling_shape)
            embedded_features = K.bias_add(
                embedded_features, self.dimensions_bias, data_format="channels_last")
        sigma = K.sigmoid(embedded_features)

        result = weighted_sum(first_, second_, sigma,
                              self.first_threshold, self.second_threshold)
        probs = K.softmax(result)
        if self.return_logits:
            return [probs, result]
        return probs

    def compute_output_shape(self, input_shape):
        first_shape = input_shape[0]
        if self.return_logits:
            return [first_shape, first_shape]
        return first_shape

def TemporalDropout(inputs, dropout=0.0):
    """
    Drops with :dropout probability temporal steps of input 3D tensor
    """
    # TO DO: adapt for >3D tensors
    if dropout == 0.0:
        return inputs
    inputs_func = lambda x: K.ones_like(inputs[:, :, 0:1])
    inputs_mask = Lambda(inputs_func)(inputs)
    inputs_mask = Dropout(dropout)(inputs_mask)
    tiling_shape = [1, 1, K.shape(inputs)[2]] + [1] * (K.ndim(inputs) - 3)
    inputs_mask = Lambda(K.tile, arguments={"n": tiling_shape},
                         output_shape=inputs._keras_shape[1:])(inputs_mask)
    answer = Multiply()([inputs, inputs_mask])
    return answer

def positions_func(inputs, pad=0):
    """
    A layer filling i-th column of a 2D tensor with
    1+ln(1+i) when it contains a meaningful symbol
    and with 0 when it contains PAD
    """
    position_inputs = K.cumsum(K.ones_like(inputs, dtype="float32"), axis=1)
    position_inputs *= K.cast(K.not_equal(inputs, pad), "float32")
    return K.log(1.0 + position_inputs)

max_stars_repo_name: xbodx/DeepPavlov

max_issues_repo_name: xbodx/DeepPavlov

max_forks_repo_name: xbodx/DeepPavlov


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

import math
import os
import sys
import time

import numpy as np
from six.moves import xrange
import tensorflow as tf

from nlp.textsum import model as seq2seq_model
from nlp.textsum.dataset import data_utils, dataset

def read_data(buckets, source_path, target_path):

    data_set = [[] for _ in buckets]
    with tf.gfile.GFile(source_path, mode="r") as source_file:
        with tf.gfile.GFile(target_path, mode="r") as target_file:
            source, target = source_file.readline(), target_file.readline()
            counter = 0
            while source and target:
                counter += 1
                if counter % 10000 == 0:
                    print("  reading data line %d" % counter)
                    sys.stdout.flush()
                source_ids = [int(x) for x in source.split()]
                target_ids = [int(x) for x in target.split()]
                target_ids.append(data_utils.EOS_ID)
                for bucket_id, (source_size, target_size) in enumerate(buckets):
                    if len(source_ids) < source_size and len(target_ids) < target_size:
                        data_set[bucket_id].append([source_ids, target_ids])
                        break
                source, target = source_file.readline(), target_file.readline()
    return data_set

def train(args, steps_per_checkpoint=10):
    if not os.path.isdir(args.utils_dir):
        os.mkdir(args.utils_dir)

    if not os.path.isdir(args.train_dir):
        os.mkdir(args.train_dir)

    # Prepare Headline data.
    print("Preparing Headline data in %s" % args.data_dir)
    src_train, dest_train, src_dev, dest_dev, vocab_size =\
        data_utils.prepare_data(args.data_dir, args.utils_dir)

    args.vocab_size = vocab_size

    # device config for CPU usage
    config = tf.ConfigProto(device_count={"CPU": 4},  # limit to 4 CPU usage
                            inter_op_parallelism_threads=1,
                            intra_op_parallelism_threads=2)  # n threads parallel for ops

    with tf.Session(config=config) as sess:
        # Create model.
        print("Creating %d layers of %d units." % (args.num_layers, args.hidden_size))
        model = seq2seq_model.create_model(sess, args.train_dir, args, False)

        # Read data into buckets and compute their sizes.
        dev_set = read_data(args.buckets, src_dev, dest_dev)
        train_set = read_data(args.buckets, src_train, dest_train)

        train_bucket_sizes = [len(train_set[b]) for b in xrange(len(args.buckets))]
        train_total_size = float(sum(train_bucket_sizes))

        trainbuckets_scale = [sum(train_bucket_sizes[:i + 1]) / train_total_size
                              for i in xrange(len(train_bucket_sizes))]

        # This is the training loop.
        step_time, loss = 0.0, 0.0
        current_step = 0
        previous_losses = []
        while current_step<args.num_steps:

            random_number_01 = np.random.random_sample()
            bucket_id = min([i for i in xrange(len(trainbuckets_scale))
                             if trainbuckets_scale[i] > random_number_01])

            # Get a batch and make a step.
            start_time = time.time()
            encoder_inputs, decoder_inputs, target_weights =\
                dataset.DataSet(args, train_set).next_batch(bucket_id)

            _, step_loss, _ = model.step(sess, encoder_inputs, decoder_inputs,
                                         target_weights, bucket_id, False)
            step_time += (time.time() - start_time) / steps_per_checkpoint
            loss += step_loss / steps_per_checkpoint
            current_step += 1

            # Once in a while, we save checkpoint, print statistics, and run evals.
            if current_step % steps_per_checkpoint == 0:
                # Print statistics for the previous epoch.
                perplexity = math.exp(float(loss)) if loss < 300 else float("inf")
                print("global step %d learning rate %.4f step-time %.2f perplexity "
                      "%.2f" % (model.global_step.eval(), model.learning_rate.eval(),
                                step_time, perplexity))
                # Decrease learning rate if no improvement was seen over last 3 times.
                if len(previous_losses) > 2 and loss > max(previous_losses[-3:]):
                    sess.run(model.learning_rate_decay_op)
                previous_losses.append(loss)
                # Save checkpoint and zero timer and loss.
                checkpoint_path = os.path.join(args.train_dir, "headline_large.ckpt")
                model.saver.save(sess, checkpoint_path, global_step=model.global_step)
                step_time, loss = 0.0, 0.0
                # Run evals on development set and print their perplexity.
                for bucket_id in xrange(len(args.buckets)):
                    if len(dev_set[bucket_id]) == 0:
                        print("eval: empty bucket %d" % (bucket_id))
                        continue
                    encoder_inputs, decoder_inputs, target_weights =\
                        dataset.DataSet(args, dev_set).next_batch(bucket_id)
                    _, eval_loss, _ = model.step(sess, encoder_inputs, decoder_inputs,
                                                 target_weights, bucket_id, True)
                    eval_ppx = math.exp(float(eval_loss)) if eval_loss < 300 else float(
                        "inf")
                    print("  eval: bucket %d perplexity %.2f" % (bucket_id, eval_ppx))
                sys.stdout.flush()


text:
#!/usr/bin/env python

"""
@package ion_functions.test.opt_functions
@file ion_functions/test/opt_functions.py
@author Christopher Wingard, Russell Desiderio, Craig Risien
@brief Unit tests for opt_functions module
"""

from nose.plugins.attrib import attr
from ion_functions.test.base_test import BaseUnitTestCase

import numpy as np
from ion_functions.data import opt_functions as optfunc

@attr('UNIT', group='func')
class TestOptFunctionsUnit(BaseUnitTestCase):

    def test_opt_functions_OPTAA_sub_functions(self):
        """
        Test the OPTAA function subroutines in the opt_functions.py module.

        Values based on test data in DPSs and available on Alfresco.

        OOI (2013). Data Product Specification for Optical Beam Attenuation
            Coefficient. Document Control Number 1341-00690.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI
            >> Controlled >> 1000 System Level >>
            1341-00690_Data_Product_SPEC_OPTATTN_OOI.pdf)

        OOI (2013). Data Product Specification for Optical Absorption
            Coefficient. Document Control Number 1341-00700.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI
            >> Controlled >> 1000 System Level >>
            1341-00700_Data_Product_SPEC_OPTABSN_OOI.pdf)

        OOI (2014). OPTAA Unit Test. 1341-00700_OPTABSN Artifact.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI >>
            >> REFERENCE >> Data Product Specification Artifacts >> 1341-00700_OPTABSN >>
            OPTAA_unit_test.xlsx)

        Implemented by Christopher Wingard, April 2013

        Modified by Russell Desiderio, 20-Feb-2014
            The opt_scatter_corr function was modified to trap out commonly
            occurring instances for when the scattering correction to absorption
            should not be applied (in these cases the additive correction term is
            set to 0). This requires that the unit test values for the abs and c
            data be physically reasonable (c >= abs). Therefore the unit test c
            data were changed by adding 1.0 to the c clear water offset.
        """

        # test inputs
        tbins = np.array([14.5036, 15.5200, 16.4706, 17.4833, 18.4831, 19.5196, 20.5565])
        tarr = np.array([
            [-0.004929, -0.004611, 0.0, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.004611, -0.004418, 0.0, 0.0, 0.0, 0.0],
            [0.0, 0.0, -0.004418, -0.004355, 0.0, 0.0, 0.0],
            [0.0, 0.0, 0.0, -0.004355, -0.004131, 0.0, 0.0],
            [0.0, 0.0, 0.0, 0.0, -0.004131, -0.003422, 0.0],
            [0.0, 0.0, 0.0, 0.0, 0.0, -0.003422, -0.002442]
        ])

        sig = np.array([50, 150, 250, 350, 450, 495])
        ref = np.array([500, 500, 500, 500, 500, 500])
        traw = np.array([48750, 48355, 47950, 47535, 47115, 46684])
        wlngth = np.array([500., 550., 600., 650., 700., 715.])
        Tcal = 20.
        T = np.array([4., 8., 12., 16., 20., 24.])
        PS = np.array([10., 15., 20., 25., 30., 35])

        # clear water offsets
        c_off = 1.01  # previous tests used 0.01
        a_off = 0.1

        # also test case of user selected reference wavelength for abs scatter correction
        ref_wave = 695  # nearest wavelength (700nm) should become reference wavelength

        # expected outputs
        tint = np.array([15., 16., 17., 18., 19., 20.])
        deltaT = np.array([-0.0048, -0.0045, -0.0044, -0.0042, -0.0038, -0.0030])
        cpd = np.array([10.2251, 5.8304, 3.7870, 2.4409, 1.4352, 1.0532])
        apd = np.array([9.3151, 4.9204, 2.8770, 1.5309, 0.5252, 0.1432])
        cpd_ts = np.array([10.226025, 5.831245, 3.795494, 2.441203, 1.440651, 1.044652])
        apd_ts = np.array([9.3155, 4.9206, 2.8848, 1.5303, 0.5297, 0.1338])
        apd_ts_s_ref715 = np.array([9.181831, 4.786862, 2.751082, 1.396591, 0.396010, 0.000000])
        apd_ts_s_ref700 = np.array([8.785990, 4.390950, 2.355159, 1.000592, 0.000000, -0.396015])

        # compute beam attenuation and optical absorption values
        dgC = np.zeros(6)
        dT = np.zeros(6)
        c = np.zeros(6)
        c_ts = np.zeros(6)
        a = np.zeros(6)
        a_ts = np.zeros(6)
        a_ts_sdef = np.zeros(6)
        a_ts_s700 = np.zeros(6)
        for i in range(6):
            dgC[i] = optfunc.opt_internal_temp(traw[i])
            c[i], dT[i] = optfunc.opt_pd_calc(ref[i], sig[i], c_off, dgC[i], tbins, tarr[i])
            c_ts[i] = optfunc.opt_tempsal_corr('c', c[i], wlngth[i], Tcal, T[i], PS[i])
            a[i], foo = optfunc.opt_pd_calc(ref[i], sig[i], a_off, dgC[i], tbins, tarr[i])
            a_ts[i] = optfunc.opt_tempsal_corr('a', a[i], wlngth[i], Tcal, T[i], PS[i])

        # the scatter-correction-to-absorption check must be done outside the loop:
        #    for each iteration within the loop, abs and c for one wavelength are
        #    processed. because there is only one wavelength of data, that wavelength
        #    also becomes the reference wavelength, in which case the scatter-corrected
        #    abs value is calculated to be  a - (a/(c-a)) * (c-a) = identically 0 for
        #    each iteration wavelength.
        # case: default reference wavelength (715nm) used for scatter correction to abs.
        a_ts_sdef = optfunc.opt_scatter_corr(a_ts, wlngth, c_ts, wlngth)
        # case: user-selected reference wavelength (695nm --> 700nm is closest)
        a_ts_s700 = optfunc.opt_scatter_corr(a_ts, wlngth, c_ts, wlngth, ref_wave)

        # compare calculated results to expected
        np.testing.assert_allclose(dgC, tint, rtol=0.1, atol=0.1)
        np.testing.assert_allclose(dT, deltaT, rtol=1e-4, atol=1e-4)
        np.testing.assert_allclose(c, cpd, rtol=1e-4, atol=1e-4)
        np.testing.assert_allclose(a, apd, rtol=1e-4, atol=1e-4)
        np.testing.assert_allclose(c_ts, cpd_ts, rtol=1e-4, atol=1e-4)
        np.testing.assert_allclose(a_ts, apd_ts, rtol=1e-4, atol=1e-4)
        np.testing.assert_allclose(a_ts_sdef, apd_ts_s_ref715, rtol=1e-4, atol=1e-4)
        np.testing.assert_allclose(a_ts_s700, apd_ts_s_ref700, rtol=1e-4, atol=1e-4)

        # now test the traps set for the cases in which unphysical scatter correction
        # calculations would be applied to the absorption data if not for the traps.
        a_ts_save = np.copy(a_ts)
        c_ts_save = np.copy(c_ts)

        # case: if a_ref < 0, do not apply the scatter correction to abs.
        # subcase: default reference wavelength
        a_ts[-1] = -0.01
        a_ts_out = np.zeros(6)
        a_ts_out = optfunc.opt_scatter_corr(a_ts, wlngth, c_ts, wlngth)
        np.testing.assert_allclose(a_ts_out, a_ts, rtol=1e-8, atol=1e-8)
        # subcase: user selected reference wavelength
        a_ts = np.copy(a_ts_save)
        a_ts[-2] = -0.01
        a_ts_out = np.zeros(6)
        a_ts_out = optfunc.opt_scatter_corr(a_ts, wlngth, c_ts, wlngth, ref_wave)
        np.testing.assert_allclose(a_ts_out, a_ts, rtol=1e-8, atol=1e-8)

        # case: if a_ref > 0 but c_ref - a_ref < 0, do not apply the scatter correction to abs.
        # subcase: default reference wavelength
        a_ts = np.copy(a_ts_save)
        a_ts[-1] = 0.01
        c_ts[-1] = 0.005
        a_ts_out = np.zeros(6)
        a_ts_out = optfunc.opt_scatter_corr(a_ts, wlngth, c_ts, wlngth)
        np.testing.assert_allclose(a_ts_out, a_ts, rtol=1e-8, atol=1e-8)
        # subcase: user selected reference wavelength
        a_ts = np.copy(a_ts_save)
        c_ts = np.copy(c_ts_save)
        a_ts[-2] = 0.01
        c_ts[-2] = 0.005
        a_ts_out = np.zeros(6)
        a_ts_out = optfunc.opt_scatter_corr(a_ts, wlngth, c_ts, wlngth, ref_wave)
        np.testing.assert_allclose(a_ts_out, a_ts, rtol=1e-8, atol=1e-8)

        # case: when both a_ref < 0 and c_ref - a_ref < 0, the scatter ratio does have
        # the correct sign; however, the scatter correction to absorption should not be
        # applied. this test is included because it is possible to code for the two
        # traps above without catching this case.
        # subcase: default reference wavelength
        a_ts = np.copy(a_ts_save)
        c_ts = np.copy(c_ts_save)
        a_ts[-1] = -0.01
        c_ts[-1] = -0.02
        a_ts_out = np.zeros(6)
        a_ts_out = optfunc.opt_scatter_corr(a_ts, wlngth, c_ts, wlngth)
        np.testing.assert_allclose(a_ts_out, a_ts, rtol=1e-8, atol=1e-8)
        # subcase: user selected reference wavelength
        a_ts = np.copy(a_ts_save)
        c_ts = np.copy(c_ts_save)
        a_ts[-2] = -0.01
        c_ts[-2] = -0.02
        a_ts_out = np.zeros(6)
        a_ts_out = optfunc.opt_scatter_corr(a_ts, wlngth, c_ts, wlngth, ref_wave)
        np.testing.assert_allclose(a_ts_out, a_ts, rtol=1e-8, atol=1e-8)

        # case: when c_ref - a_ref = 0, do not apply the scatter correction to abs.
        # subcase: default reference wavelength
        a_ts = np.copy(a_ts_save)
        c_ts = np.copy(c_ts_save)
        a_ts[-1] = 0.01
        c_ts[-1] = 0.01
        a_ts_out = np.zeros(6)
        a_ts_out = optfunc.opt_scatter_corr(a_ts, wlngth, c_ts, wlngth)
        np.testing.assert_allclose(a_ts_out, a_ts, rtol=1e-8, atol=1e-8)
        # subcase: user selected reference wavelength
        a_ts = np.copy(a_ts_save)
        c_ts = np.copy(c_ts_save)
        a_ts[-2] = 0.01
        c_ts[-2] = 0.01
        a_ts_out = np.zeros(6)
        a_ts_out = optfunc.opt_scatter_corr(a_ts, wlngth, c_ts, wlngth, ref_wave)
        np.testing.assert_allclose(a_ts_out, a_ts, rtol=1e-8, atol=1e-8)

    def test_opt_functions_OPTAA_wrapper_functions(self):
        """
        Test the OPTAA wrapper functions in the opt_functions.py module.
        Use realistically shaped ac-s data arrays in the calling arguments:
            offsets are dimensioned at the number of wavelengths;
            for internal temperature, L1 temperature, and salinity,
                one value per a-c dataset.

        Values based on test data in DPSs and available on Alfresco.

        OOI (2013). Data Product Specification for Optical Beam Attenuation
            Coefficient. Document Control Number 1341-00690.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI
            >> Controlled >> 1000 System Level >>
            1341-00690_Data_Product_SPEC_OPTATTN_OOI.pdf)

        OOI (2013). Data Product Specification for Optical Absorption
            Coefficient. Document Control Number 1341-00700.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI
            >> Controlled >> 1000 System Level >>
            1341-00700_Data_Product_SPEC_OPTABSN_OOI.pdf)

        OOI (2014). OPTAA Unit Test. 1341-00700_OPTABSN Artifact.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI >>
            >> REFERENCE >> Data Product Specification Artifacts >> 1341-00700_OPTABSN >>
            OPTAA_unit_test.xlsx)

        Implemented by Russell Desiderio, 21-Feb-2014. Additional unit test data
            constructed as documented in the OPTAA Unit Test document artifact
            reference above.
        2015-04-17: Russell Desiderio. Changed signal and reference count input from float to fix.
        """

        # test inputs:
        tbins = np.array([14.5036, 15.5200, 16.4706, 17.4833, 18.4831, 19.5196, 20.5565])
        # nrows of tarr = length(wlngth); ncols of tarr = length(tbins)
        tarr = np.array([
            [0.0, -0.004929, -0.004611, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.004611, -0.004418, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.004418, -0.004355, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.004355, -0.004131, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.004131, -0.003422, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.003422, -0.002442, 0.0, 0.0, 0.0, 0.0]
        ])

        # a_off and c_off must have same dimensions as wlngth.
        wlngth = np.array([500., 550., 600., 650., 700., 715.])
        c_sig = np.array([150, 225, 200, 350, 450, 495])
        c_ref = np.array([550, 540, 530, 520, 510, 500])
        c_off = np.array([1.35, 1.30, 1.25, 1.20, 1.15, 1.10])
        a_sig = np.array([250, 300, 210, 430, 470, 495])
        a_ref = np.array([450, 460, 470, 480, 490, 500])
        a_off = np.array([0.35, 0.30, 0.25, 0.20, 0.15, 0.10])

        # traw, T, and PS must be scalar (before vectorization).
        traw = 48355
        Tcal = 20.0
        T = 12.0
        PS = 35.0

        # also test case of user selected reference wavelength for abs scatter correction
        ref_wave = 695  # nearest wavelength (700nm) should become reference wavelength

        # expected outputs
        cpd_ts = np.array([6.553646, 4.807949, 5.161655, 2.788220, 1.666362, 1.184530], ndmin=2)
        apd_ts_s_ref715 = np.array([1.999989, 1.501766, 3.177048, 0.249944, 0.086438, 0.000000],
                                   ndmin=2)
        apd_ts_s_ref700 = np.array([1.750859, 1.320885, 3.068470, 0.111075, 0.000000, -0.064806],
                                   ndmin=2)

        # beam attenuation (beam c) wrapper test
        c_ts = np.zeros(6)
        c_ts = optfunc.opt_beam_attenuation(c_ref, c_sig, traw, wlngth, c_off, Tcal,
                                            tbins, tarr, T, PS)
        np.testing.assert_allclose(c_ts, cpd_ts, rtol=1e-6, atol=1e-6)

        # absorption wrapper test
        # case: default reference wavelength for scatter correction
        a_ts_sdef = np.zeros(6)
        a_ts_sdef = optfunc.opt_optical_absorption(a_ref, a_sig, traw, wlngth, a_off, Tcal,
                                                   tbins, tarr, cpd_ts, wlngth, T, PS)
        np.testing.assert_allclose(a_ts_sdef, apd_ts_s_ref715, rtol=1e-6, atol=1e-6)
        # case: user selected reference wavelength for scatter correction
        a_ts_s700 = np.zeros(6)
        a_ts_s700 = optfunc.opt_optical_absorption(a_ref, a_sig, traw, wlngth, a_off, Tcal,
                                                   tbins, tarr, cpd_ts, wlngth, T, PS, ref_wave)
        np.testing.assert_allclose(a_ts_s700, apd_ts_s_ref700, rtol=1e-6, atol=1e-6)

    def test_opt_functions_OPTAA_a_and_c_wavelength_sets(self):
        """
        Test the OPTAA wrapper functions in the opt_functions.py module.
            Test a-c dataset with different wavelength values for the absorption versus
            beam c optical channels. Included to test that the c optical values are correctly
            interpolated onto the abs wavelengths in the scatter correction algorithm.

        OOI (2013). Data Product Specification for Optical Beam Attenuation
            Coefficient. Document Control Number 1341-00690.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI
            >> Controlled >> 1000 System Level >>
            1341-00690_Data_Product_SPEC_OPTATTN_OOI.pdf)

        OOI (2013). Data Product Specification for Optical Absorption
            Coefficient. Document Control Number 1341-00700.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI
            >> Controlled >> 1000 System Level >>
            1341-00700_Data_Product_SPEC_OPTABSN_OOI.pdf)

        OOI (2014). OPTAA Unit Test. 1341-00700_OPTABSN Artifact.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI >>
            >> REFERENCE >> Data Product Specification Artifacts >> 1341-00700_OPTABSN >>
            OPTAA_unit_test.xlsx)

        Implemented by Russell Desiderio, 26-Feb-2014. Additional unit test data
            constructed as documented in the OPTAA Unit Test document artifact
            reference above.
        2015-04-17: Russell Desiderio. Changed signal and reference count input from float to fix.
        """

        # test inputs:
        tbins = np.array([14.5036, 15.5200, 16.4706, 17.4833, 18.4831, 19.5196, 20.5565])
        # nrows of tarr = length(wlngth); ncols of tarr = length(tbins)
        tarr = np.array([
            [0.0, -0.004929, -0.004611, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.004611, -0.004418, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.004418, -0.004355, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.004355, -0.004131, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.004131, -0.003422, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.003422, -0.002442, 0.0, 0.0, 0.0, 0.0]
        ])

        # a_off and c_off must have same dimensions as wlngth.
        cwlngth = np.array([510., 540., 580., 630., 670., 710.])
        awlngth = np.array([500., 550., 600., 650., 700., 715.])
        c_sig = np.array([150, 225, 200, 350, 450, 495])
        c_ref = np.array([550, 540, 530, 520, 510, 500])
        c_off = np.array([1.35, 1.30, 1.25, 1.20, 1.15, 1.10])
        a_sig = np.array([250, 300, 210, 430, 470, 495])
        a_ref = np.array([450, 460, 470, 480, 490, 500])
        a_off = np.array([0.35, 0.30, 0.25, 0.20, 0.15, 0.10])

        # traw, T, and PS must be scalar (before vectorization).
        traw = 48355
        Tcal = 20.0
        T = 12.0
        PS = 35.0

        # also test case of user selected reference wavelength for abs scatter correction
        ref_wave = 695  # nearest abs wavelength (700nm) should become reference wavelength

        # expected outputs
        cpd_ts = np.array([6.553771, 4.807914, 5.156010, 2.788715, 1.655607, 1.171965], ndmin=2)
        apd_ts_s_ref715 = np.array([1.990992, 1.479089, 3.350109, 0.350108, 0.152712, 0.000000],
                                   ndmin=2)
        apd_ts_s_ref700 = np.array([1.379858, 1.021574, 3.235057, 0.099367, 0.000000, -0.156972],
                                   ndmin=2)

        # beam attenuation (beam c) wrapper test
        c_ts = np.zeros(6)
        c_ts = optfunc.opt_beam_attenuation(c_ref, c_sig, traw, cwlngth, c_off, Tcal,
                                            tbins, tarr, T, PS)
        np.testing.assert_allclose(c_ts, cpd_ts, rtol=1e-6, atol=1e-6)

        # absorption wrapper test
        # case: default reference wavelength for scatter correction
        a_ts_sdef = np.zeros(6)
        a_ts_sdef = optfunc.opt_optical_absorption(a_ref, a_sig, traw, awlngth, a_off, Tcal,
                                                   tbins, tarr, cpd_ts, cwlngth, T, PS)
        np.testing.assert_allclose(a_ts_sdef, apd_ts_s_ref715, rtol=1e-6, atol=1e-6)
        # case: user selected reference wavelength for scatter correction
        a_ts_s700 = np.zeros(6)
        a_ts_s700 = optfunc.opt_optical_absorption(a_ref, a_sig, traw, awlngth, a_off, Tcal,
                                                   tbins, tarr, cpd_ts, cwlngth, T, PS, ref_wave)
        np.testing.assert_allclose(a_ts_s700, apd_ts_s_ref700, rtol=1e-6, atol=1e-6)

    def test_opt_functions_OPTAA_vectorization(self):
        """
        Test the vectorization of the OPTAA wrapper functions in the opt_functions.py module.

        OOI (2013). Data Product Specification for Optical Beam Attenuation
            Coefficient. Document Control Number 1341-00690.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI
            >> Controlled >> 1000 System Level >>
            1341-00690_Data_Product_SPEC_OPTATTN_OOI.pdf)

        OOI (2013). Data Product Specification for Optical Absorption
            Coefficient. Document Control Number 1341-00700.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI
            >> Controlled >> 1000 System Level >>
            1341-00700_Data_Product_SPEC_OPTABSN_OOI.pdf)

        OOI (2014). OPTAA Unit Test. 1341-00700_OPTABSN Artifact.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI >>
            >> REFERENCE >> Data Product Specification Artifacts >> 1341-00700_OPTABSN >>
            OPTAA_unit_test.xlsx)

        Implemented by Russell Desiderio, 05-Mar-2014. Additional unit test data
            constructed as documented in the OPTAA Unit Test document artifact
            reference above.
        2015-04-17: Russell Desiderio. Changed signal and reference count input from float to fix.
                                       Replaced 'exec' statements.
        2015-04-21: Russell Desiderio. Corrected input data shapes to current CI implementation.
                                       Added scalar time record case.
        """
        ### set test inputs: scalar (that is, 1) time record case
        # native scalars will come into the DPA with shape (1,)
        traw = np.array([48355])
        Tcal = np.array([20.0])
        T = np.array([12.0])
        PS = np.array([35.0])
        # native 1D arrays will come in as 2D row vectors
        cwlngth = np.array([[510., 540., 580., 630., 670., 710.]])
        awlngth = np.array([[500., 550., 600., 650., 700., 715.]])
        c_sig = np.array([[150, 225, 200, 350, 450, 495]])
        c_ref = np.array([[550, 540, 530, 520, 510, 500]])
        c_off = np.array([[1.35, 1.30, 1.25, 1.20, 1.15, 1.10]])
        a_sig = np.array([[250, 300, 210, 430, 470, 495]])
        a_ref = np.array([[450, 460, 470, 480, 490, 500]])
        a_off = np.array([[0.35, 0.30, 0.25, 0.20, 0.15, 0.10]])
        tbins = np.array([[14.5036, 15.5200, 16.4706, 17.4833, 18.4831, 19.5196, 20.5565]])
        # native 2D arrays will come in as 3D arrays, for this single time record case;
        # a singleton dimension will be prepended to the native 2D array.
        # native dimensions of tarr are (nwavelengths, ntempbins),
        # so that tarr.shape = (1,6,7)
        tarr = np.array([[
            [0.0, -0.004929, -0.004611, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.004611, -0.004418, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.004418, -0.004355, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.004355, -0.004131, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.004131, -0.003422, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.003422, -0.002442, 0.0, 0.0, 0.0, 0.0]
        ]])

        ## expected outputs to scalar time record test: 2D row vectors
        # beam attenuation (beam c) test
        # expected:
        cpd_ts = np.array([[6.553771, 4.807914, 5.156010, 2.788715, 1.655607, 1.171965]])
        # calculated:
        c_ts = optfunc.opt_beam_attenuation(c_ref, c_sig, traw, cwlngth, c_off, Tcal,
                                            tbins, tarr, T, PS)
        np.testing.assert_allclose(c_ts, cpd_ts, rtol=1e-6, atol=1e-6)

        # absorption test
        # case: default reference wavelength for scatter correction
        apd_ts_s_ref715 = np.array([[1.990992, 1.479089, 3.350109, 0.350108, 0.152712, 0.000000]])
        a_ts_sdef = optfunc.opt_optical_absorption(a_ref, a_sig, traw, awlngth, a_off, Tcal,
                                                   tbins, tarr, cpd_ts, cwlngth, T, PS)
        np.testing.assert_allclose(a_ts_sdef, apd_ts_s_ref715, rtol=1e-6, atol=1e-6)

        # absorption test
        # case: user selected reference wavelength for scatter correction to abs values
        ref_wave = 695  # nearest abs wavelength (700nm) should become reference wavelength

        apd_ts_s_ref700 = np.array([[1.379858, 1.021574, 3.235057, 0.099367, 0.000000, -0.156972]])
        a_ts_s700 = optfunc.opt_optical_absorption(a_ref, a_sig, traw, awlngth, a_off, Tcal,
                                                   tbins, tarr, cpd_ts, cwlngth, T, PS, ref_wave)
        np.testing.assert_allclose(a_ts_s700, apd_ts_s_ref700, rtol=1e-6, atol=1e-6)

        ### multiple time records case
        # replicate the inputs to represent 3 data packets
        npackets = 3
        [traw, Tcal, T, PS] = [np.tile(xarray, npackets) for xarray in [traw, Tcal, T, PS]]
        [cwlngth, awlngth, c_sig, c_ref, c_off, a_sig, a_ref, a_off, tbins] = [
            np.tile(xarray, (npackets, 1)) for xarray in [
                cwlngth, awlngth, c_sig, c_ref, c_off, a_sig, a_ref, a_off, tbins]]
        tarr = np.tile(tarr, (npackets, 1, 1))

        # replicate the expected output data products
        cpd_ts = np.tile(cpd_ts, (npackets, 1))
        apd_ts_s_ref715 = np.tile(apd_ts_s_ref715, (npackets, 1))
        apd_ts_s_ref700 = np.tile(apd_ts_s_ref700, (npackets, 1))

        # beam attenuation (beam c) test
        c_ts = optfunc.opt_beam_attenuation(c_ref, c_sig, traw, cwlngth, c_off, Tcal,
                                            tbins, tarr, T, PS)
        np.testing.assert_allclose(c_ts, cpd_ts, rtol=1e-6, atol=1e-6)

        # absorption test
        # case: default reference wavelength for scatter correction
        a_ts_sdef = optfunc.opt_optical_absorption(a_ref, a_sig, traw, awlngth, a_off, Tcal,
                                                   tbins, tarr, cpd_ts, cwlngth, T, PS)
        np.testing.assert_allclose(a_ts_sdef, apd_ts_s_ref715, rtol=1e-6, atol=1e-6)

        # case: user selected reference wavelength for scatter correction
        a_ts_s700 = optfunc.opt_optical_absorption(a_ref, a_sig, traw, awlngth, a_off, Tcal,
                                                   tbins, tarr, cpd_ts, cwlngth, T, PS, ref_wave)
        np.testing.assert_allclose(a_ts_s700, apd_ts_s_ref700, rtol=1e-6, atol=1e-6)

    def test_opt_functions_OPTAA_vectorization_with_tscor_nans(self):
        """
        Test the vectorized OPTAA wrapper functions in the opt_functions.py module;
        include "a" and "c" test wavelengths outside of the (original) range of the
        wavelength keys in the tscor.py file, which contains the dictionary of the
        temperature and salinity correction coefficients. The dictionary keys have
        been extended from [400.0 755.0] to [380.0 775.0] with entries of np.nan.

        OOI (2013). Data Product Specification for Optical Beam Attenuation
            Coefficient. Document Control Number 1341-00690.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI
            >> Controlled >> 1000 System Level >>
            1341-00690_Data_Product_SPEC_OPTATTN_OOI.pdf)

        OOI (2013). Data Product Specification for Optical Absorption
            Coefficient. Document Control Number 1341-00700.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI
            >> Controlled >> 1000 System Level >>
            1341-00700_Data_Product_SPEC_OPTABSN_OOI.pdf)

        OOI (2014). OPTAA Unit Test. 1341-00700_OPTABSN Artifact.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI >>
            >> REFERENCE >> Data Product Specification Artifacts >> 1341-00700_OPTABSN >>
            OPTAA_unit_test.xlsx)

        Initial code by Russell Desiderio, 29-Apr-2014. Added tests for when OPTAA
            wavelengths are outside the range of the empirically derived temperature
            and salinity correction wavelength keys in the original tscor.py dictionary
            by modifying def test_opt_functions_OPTAA_vectorization test data.
        2015-04-17: Russell Desiderio. Use np.nan instead of fill_value.
                                       Changed signal and reference count input from float to fix.
                                       Replaced 'exec' statements.
        2015-04-21: Russell Desiderio. Cosmetic statement re-ordering and shaping of input arrays.
        """
        # test inputs:
        # replicate the inputs to represent 3 data packets
        npackets = 3

        ## natively scalar inputs
        traw = np.array([48355])
        Tcal = np.array([20.0])
        T = np.array([12.0])
        PS = np.array([35.0])
        # time-vectorize
        [traw, Tcal, T, PS] = [np.tile(xarray, npackets) for xarray in [traw, Tcal, T, PS]]

        ## natively 1D inputs
        # valid tscor dictionary entries (those without nan values) are
        # keyed at [400.0 755.0] nm. The dictionary has been extended to
        # to [380.0 775.0] using np.nan as fill values. test:
        cwlngth = np.array([[398.5, 540., 580., 630., 670., 710.]])
        awlngth = np.array([[500., 550., 600., 650., 700., 761.2]])
        c_sig = np.array([[150, 225, 200, 350, 450, 495]])
        c_ref = np.array([[550, 540, 530, 520, 510, 500]])
        c_off = np.array([[1.35, 1.30, 1.25, 1.20, 1.15, 1.10]])
        a_sig = np.array([[250, 300, 210, 430, 470, 495]])
        a_ref = np.array([[450, 460, 470, 480, 490, 500]])
        a_off = np.array([[0.35, 0.30, 0.25, 0.20, 0.15, 0.10]])
        tbins = np.array([[14.5036, 15.5200, 16.4706, 17.4833, 18.4831, 19.5196, 20.5565]])
        # time-vectorize
        [cwlngth, awlngth, c_sig, c_ref, c_off, a_sig, a_ref, a_off, tbins] = [
            np.tile(xarray, (npackets, 1)) for xarray in [
                cwlngth, awlngth, c_sig, c_ref, c_off, a_sig, a_ref, a_off, tbins]]

        # nrows of tarr = length(nwlngth); ncols of tarr = length(tbins)
        tarr = np.array([[
            [0.0, -0.004929, -0.004611, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.004611, -0.004418, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.004418, -0.004355, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.004355, -0.004131, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.004131, -0.003422, 0.0, 0.0, 0.0, 0.0],
            [0.0, -0.003422, -0.002442, 0.0, 0.0, 0.0, 0.0]
        ]])
        # time-vectorize
        # dimensions will be (npackets, nwavelengths, ntempbins)
        tarr = np.tile(tarr, (npackets, 1, 1))

        # expected outputs:
        # 1st c wavelength is out-of-range of the empirical TS corrections, so the DPA
        # should return nan.
        cpd_ts = np.array([[np.nan, 4.807914, 5.156010, 2.788715, 1.655607, 1.171965]])
        # to scatter-correct the 1st abs wavelength value, the c values are interpolated;
        # because the 1st c value is out-of-TScor-range, the 1st abs value will also be
        # a nan value. and, the last abs value will also be a nan value because 761.2>755.0.
        apd_ts_s_ref700 = np.array([[np.nan, 1.021574, 3.235057, 0.099367, 0.000000, np.nan]])

        cpd_ts = np.tile(cpd_ts, (npackets, 1))
        apd_ts_s_ref700 = np.tile(apd_ts_s_ref700, (npackets, 1))

        # beam attenuation (beam c)  test
        c_ts = optfunc.opt_beam_attenuation(c_ref, c_sig, traw, cwlngth, c_off, Tcal,
                                            tbins, tarr, T, PS)
        np.testing.assert_allclose(c_ts, cpd_ts, rtol=1e-6, atol=1e-6)

        # absorption test
        a_ts_s700 = optfunc.opt_optical_absorption(a_ref, a_sig, traw, awlngth, a_off, Tcal,
                                                   tbins, tarr, cpd_ts, cwlngth, T, PS)
        np.testing.assert_allclose(a_ts_s700, apd_ts_s_ref700, rtol=1e-6, atol=1e-6)

    def test_opt_par_satlantic(self):
        """
        Test the opt_par_satlantic function.

        Values based on that described in DPS as available on Alfresco:

        OOI (2012). Data Product Specification for PHOTOSYNTHETICALLY
        ACTIVE RADIATION (PAR) FROM SATLANTIC INSTRUMENT ON RSN SHALLOW
        PROFILER Document Control Number 1341-00720.
        https://alfresco.oceanobservatories.org/
        (See: Company Home >> OOI >> Controlled >> 1000 System Level >>
        1341-00720_Data_Product_SPEC_OPTPARW_Satl_OOI.pdf)

        Implemented by Craig Risien, February 2014
        Modified by Russell Desiderio, April 09, 2015.
            check function operation without iterating on the input data using for loop.
            check both cal coeff cases: (1) calcoeffs are scalars
                                        (2) calcoeffs are 1d arrays in time.
        Modified by Russell Desiderio, April 21, 2015.
            removed test using 1D data input array with scalar calcoeffs.
            added scalar data with scalar calcoeff test.
        """

    # test inputs
    Im = np.array([1.3589])
    a0 = np.array([2156849800.8])
    a1 = np.array([2.586852835e-006])

    par_count = np.array([2159403328, 2159400384, 2159396992, 2159400384, 2159407488, 2159399296, 2159400384,
                          2159404800, 2159403904, 2159402240, 2159403200, 2159409728, 2159409792, 2159408320,
                          2159407808, 2159402304, 2159402688, 2159407552, 2159404160, 2159403776, 2159402048,
                          2159404544])

    # expected outputs
    par_expected = np.array([8.976348585, 8.965999618, 8.954075807, 8.965999618, 8.990972126, 8.962174999,
                             8.965999618, 8.981523069, 8.978373383, 8.972523967, 8.97589863, 8.998846341,
                             8.999071318, 8.993896835, 8.992097014, 8.972748944, 8.97409881, 8.991197104,
                             8.979273293, 8.977923428, 8.971849034, 8.980623159])

    # compute par value for 'scalar' in time case.
    par_calc = optfunc.opt_par_satlantic(par_count[0], a0, a1, Im)
    # compare calculated results to expected
    np.testing.assert_allclose(par_calc, par_expected[0], rtol=0.000001, atol=0.000001)

    # compute par values: time-vectorized cal coeff case
    # all inputs and output are 1D arrays of shape (tval,)
    tval = par_count.shape[0]
    a0_tiled = np.tile(a0, tval)
    a1_tiled = np.tile(a1, tval)
    Im_tiled = np.tile(Im, tval)
    par_calc = optfunc.opt_par_satlantic(par_count, a0_tiled, a1_tiled, Im_tiled)
    # compare calculated results to expected
    np.testing.assert_allclose(par_calc, par_expected, rtol=0.000001, atol=0.000001)

    def test_opt_par_wetlabs(self):
        """
        Test the opt_par_wetlabs function.

        Values based on that described in DPS as available on Alfresco:

        OOI (2014). Data Product Specification for DATA PRODUCT SPECIFICATION
        FOR PHOTOSYNTHETICALLY ACTIVE RADIATION (PAR) FROM WET LABS INSTRUMENT
        ON COASTAL SURFACE PIERCING PROFILER Document Control Number 1341-00722.
        https://alfresco.oceanobservatories.org/
        (See: Company Home >> OOI >> Controlled >> 1000 System Level >>
        1341-00722_Data_Product_SPEC_OPTPARW_WETLabs_OOI.pdf)

        Implemented by Craig Risien, December 2014
        Modified by Russell Desiderio, April 09, 2015.
            check function operation without iterating on the input data using for loop.
            check both cal coeff cases: (1) calcoeffs are scalars
                                        (2) calcoeffs are 1d arrays in time.
        Modified by Russell Desiderio, April 21, 2015.
            removed test using 1D data input array with scalar calcoeffs.
            added scalar data with scalar calcoeff test.
        """

    # test inputs
    Im = np.array([1.3589])
    a0 = np.array([4381])
    a1 = np.array([2904])

    par_count = np.array([4975, 4999, 5077, 5230, 5413, 5483, 5427, 5344, 5285, 5284,
                          5344, 5470, 5649, 5835, 6007, 11824, 4975, 4999, 5077, 5230,
                          5413, 11824])

    # expected outputs
    par_expected = np.array([2.176371115, 2.218183222, 2.359700693, 2.664053061, 3.08006158,
                             3.255847708, 3.114442636, 2.916077547, 2.782801729, 2.780596116,
                             2.916077547, 3.222459729, 3.713869311, 4.304039059, 4.932928879,
                             496.8256708, 2.176371115, 2.218183222, 2.359700693, 2.664053061,
                             3.08006158, 496.8256708])

    # compute par value for 'scalar' in time case.
    par_calc = optfunc.opt_par_wetlabs(par_count[0], a0, a1, Im)
    # compare calculated results to expected
    np.testing.assert_allclose(par_calc, par_expected[0], rtol=0.000001, atol=0.000001)

    # compute par values: time-vectorized cal coeff case
    tval = par_count.shape[0]
    a0_tiled = np.tile(a0, tval)
    a1_tiled = np.tile(a1, tval)
    Im_tiled = np.tile(Im, tval)
    par_calc = optfunc.opt_par_wetlabs(par_count, a0_tiled, a1_tiled, Im_tiled)
    # compare calculated results to expected
    np.testing.assert_allclose(par_calc, par_expected, rtol=0.000001, atol=0.000001)

    def opt_par_biospherical_mobile(self):
                """
        Test the opt_par_biospherical_mobile function.

        Values based on that described in DPS as available on Alfresco:

        OOI (2012). Data Product Specification for PHOTOSYNTHETICALLY
        ACTIVE RADIATION (PAR) FROM BIOSPHERICAL INSTRUMENT QSP 2100
        ON CGSN MOBILE ASSETS Document Control Number 1341-00721.
        https://alfresco.oceanobservatories.org/
        (See: Company Home >> OOI >> Controlled >> 1000 System Level >>
        1341-00721_Data_Product_SPEC_OPTPARW_Bios_OOI.pdf)

        Implemented by Craig Risien, February 2014
        Modified by Russell Desiderio, April 09, 2015.
            check function operation without iterating on the input data using for loop.
            check both cal coeff cases: (1) calcoeffs are scalars
                                        (2) calcoeffs are 1d arrays in time.
        Modified by Russell Desiderio, April 21, 2015.
            removed test using 1D data input array with scalar calcoeffs.
            added scalar data with scalar calcoeff test.
        """

    # test inputs
    offset = np.array([0.0101])
    scale = np.array([5.897e-04])

    par_volts = np.array([1.016793, 0.599800, 0.452400, 0.305000, 0.187000, 0.178900, 0.069100, 0.039600, 0.010100])

    par_expected = np.array([1707.13, 1000.00, 750.00, 500.00, 300.00, 286.25, 100.00, 50.00, 0.00])

    # scalar time case
    par_calc = optfunc.opt_par_biospherical_mobile(par_volts[0], offset, scale)
    # compare calculated results to expected
    np.testing.assert_allclose(par_calc, par_expected[0], rtol=0.01, atol=0.01)

    # compute par values: time-vectorized cal coeff case
    tval = par_volts.shape[0]
    offset_tiled = np.tile(offset, tval)
    scale_tiled = np.tile(scale, tval)
    par_calc = optfunc.opt_par_biospherical_mobile(par_volts, offset_tiled, scale_tiled)
    # compare calculated results to expected
    np.testing.assert_allclose(par_calc, par_expected, rtol=0.01, atol=0.01)

    def opt_par_biospherical_wfp(self):
                """
        Test the opt_par_biospherical_wfp function.

        Values based on that described in DPS as available on Alfresco:

        OOI (2012). Data Product Specification for PHOTOSYNTHETICALLY
        ACTIVE RADIATION (PAR) FROM BIOSPHERICAL INSTRUMENT QSP 2200
        ON CGSN PROFILERS Document Control Number 1341-00721.
        https://alfresco.oceanobservatories.org/
        (See: Company Home >> OOI >> Controlled >> 1000 System Level >>
        1341-00721_Data_Product_SPEC_OPTPARW_Bios_OOI.pdf)

        Implemented by Craig Risien, March 2014
        Modified by Russell Desiderio, April 09, 2015.
            check function operation without iterating on the input data using for loop.
            check both cal coeff cases: (1) calcoeffs are scalars
                                        (2) calcoeffs are 1d arrays in time.
        Modified by Russell Desiderio, April 21, 2015.
            removed test using 1D data input array with scalar calcoeffs.
            added scalar data with scalar calcoeff test.
        """

    # test inputs
    #offset units = mvolts
    offset = np.array([10.1])
    scale = np.array([5.897e-04 / (6.02 * 10**13)])

    par_mvolts = np.array([1016.793, 599.800, 452.400, 305.000, 187.000, 178.900, 69.100, 39.600, 10.100])

    par_expected = np.array([1707.13, 1000.00, 750.00, 500.00, 300.00, 286.25, 100.00, 50.00, 0.00])

    # scalar case
    par_calc = optfunc.opt_par_biospherical_wfp(par_mvolts[0], offset, scale)
    # compare calculated results to expected
    np.testing.assert_allclose(par_calc, par_expected[0], rtol=0.01, atol=0.01)

    # compute par values: time-vectorized cal coeff case
    tval = par_mvolts.shape[0]
    offset_tiled = np.tile(offset, tval)
    scale_tiled = np.tile(scale, tval)
    par_calc = optfunc.opt_par_biospherical_wfp(par_mvolts, offset_tiled, scale_tiled)
    # compare calculated results to expected
    np.testing.assert_allclose(par_calc, par_expected, rtol=0.01, atol=0.01)

    def test_opt_ocr507_irradiance(self):
        """
        Test opt_ocr507_irradiance function.

        Values calculated using data product algorithm in the DPS in an excel spreadsheet
        as available on Alfresco:

        OOI (2014). Data Product Specification for Downwelling Spectral Irradiance.
            Document Control Number 1341-00730. https://alfresco.oceanobservatories.org/
            (See: Company Home >> OOI >> Controlled >> 1000 System Level >>
            1341-00730__???.pdf)

        OOI (2014). SPECTIR Unit Test. 1341-00730_SPECTIR Artifact.
            https://alfresco.oceanobservatories.org/ (See: Company Home >> OOI >>
            >> REFERENCE >> Data Product Specification Artifacts >> 1341-00730_SPECTIR >>
            SPKIR_SPECTIR_unit_test.xlsx)

        Implemented by Russell Desiderio, March 14, 2014.
        Modified by Russell Desiderio, March 25, 2014.
            Transposed test array.
            Set known output to be 2D row vector.
        Modified by Russell Desiderio, April 09, 2015.
            Added unit tests for multiple data packets and exception checking
            Set known output for 1 data packet case to be 1D array.
        Modified by Russell Desiderio, April 21, 2015.
            Conditioned inputs for 1 data packet case to be 2D row vectors.
            Set known output for 1 data packet case back to 2D row vector.
        """
#              counts         offset          scale        mrsn     Ed
        test_array = np.transpose(np.array([
            [2148370944, 2148377867.8, 2.09023117662E-07, 1.368, -0.00198],
            [2200000000, 2148218092.4, 2.06543624674E-07, 1.410, 15.080],
            [2300000000, 2147607229.7, 2.12484770952E-07, 1.365, 44.200],
            [2400000000, 2147789959.1, 2.07241106309E-07, 1.354, 70.771],
            [2500000000, 2148047456.7, 1.99358530187E-07, 1.372, 96.266],
            [2600000000, 2147335412.8, 2.06033896796E-07, 1.404, 130.943],
            [2700000000, 2146998228.4, 2.14806273478E-07, 1.347, 160.008]
        ]))

        ## one data packet case
        # set inputs
        counts = test_array[[0], :]
        offset = test_array[[1], :]
        scale = test_array[[2], :]
        immersion_factor = test_array[[3], :]

        # set known output to be a 2D row vector
        Ed = test_array[[-1], :]
        #print Ed.shape

        # calculate the downwelling irradiance
        Ed_out = optfunc.opt_ocr507_irradiance(counts, offset, scale, immersion_factor)

        ###########################################################################
        # The DPS specifies a precision of 0.25 uW/cm^2/nm
        np.testing.assert_allclose(Ed_out, Ed, rtol=0.0, atol=0.1)
        ###########################################################################

        ## multiple data packets case
        ## cal coeffs will also be "time-vectorized" by CI
        # create 10 data packets
        z_value = 10
        zcounts = np.tile(counts, (z_value, 1))
        zoffset = np.tile(offset, (z_value, 1))
        zscale = np.tile(scale, (z_value, 1))
        zimmersion_factor = np.tile(immersion_factor, (z_value, 1))

        # set known output
        zEd = np.tile(Ed, (z_value, 1))
        #print Ed.shape

        # calculate the downwelling irradiance
        zEd_out = optfunc.opt_ocr507_irradiance(zcounts, zoffset, zscale, zimmersion_factor)

        ###########################################################################
        # The DPS specifies a precision of 0.25 uW/cm^2/nm
        np.testing.assert_allclose(zEd_out, zEd, rtol=0.0, atol=0.1)
        ###########################################################################

        ## test error-checking
        # data array counts does not have 7 elements
        counts_wrongshape = test_array[0, 0:-1]
        np.testing.assert_raises(ValueError, optfunc.opt_ocr507_irradiance, counts_wrongshape,
                                 offset, scale, immersion_factor)

        # cal coeff scale does not have 7 elements
        scale_wrongshape = test_array[2, 0:-1]
        np.testing.assert_raises(ValueError, optfunc.opt_ocr507_irradiance, counts, offset,
                                 scale_wrongshape, immersion_factor)

max_stars_repo_name: steinermg/ion-functions

max_issues_repo_name: steinermg/ion-functions

max_forks_repo_name: steinermg/ion-functions


text:
from numpy import *

def neighboringelements(t2t, elem):

    t2te = t2t[:,elem];
    nbelem = sort(unique(t2te.flatten('F')));
    if nbelem[0] == -1:
        nbelem = nbelem[1:];

    return nbelem

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim

reason: from numpy


text:
# An example of Markov decision process
from mdp import MDP
import numpy as np
from mdp import FrozenLakeEnv
import matplotlib.pyplot as plt

def get_action_value(mdp, state_values, state, action, gamma):
    Q = 0
    state_primes = mdp.get_next_states(state, action)
    for state_prime, prob in state_primes.items():
        reward = mdp.get_reward(state, action, state_prime)
        Q += prob * (reward + gamma * state_values[state_prime])
    return Q

def get_new_state_value(mdp, state_values, state, gamma):
    if mdp.is_terminal(state):
        return 0
    actions = mdp.get_possible_actions(state)
    new_state_value = max([get_action_value(mdp, state_values, state, action, gamma) for action in actions])
    return new_state_value

def get_optimal_action(mdp, state_values, state, gamma=0.9):
    """ Finds optimal action using formula above. """
    if mdp.is_terminal(state):
        return None
    actions = mdp.get_possible_actions(state)
    final_action_value = {action: get_action_value(mdp, state_values, state, action, gamma) for action in actions}
    optimal_action = sorted(final_action_value, key=lambda x:final_action_value[x], reverse=True)[0]
    return optimal_action

def value_iteration(mdp, state_values=None, gamma=0.9, num_iter=1000, min_difference=1e-5):
    """ performs num_iter value iteration steps starting from state_values. """
    state_values = state_values or {s: 0 for s in mdp.get_all_states()}
    for i in range(num_iter):

        # Compute new state values using the functions you defined above.
        # It must be a dict {state : new_V(state)}
        new_state_values = {state: get_new_state_value(mdp, state_values, state, gamma) for state, values in state_values.items()}

        # Compute difference
        diff = max(abs(new_state_values[s] - state_values[s]) for s in mdp.get_all_states())
        print("iter %4i   |   diff: %6.5f   |   V(start): %.3f " % (i, diff, new_state_values[mdp._initial_state]))

        state_values = new_state_values
        if diff < min_difference:
            print("Iteration Terminated")
            break

    return state_values

def draw_policy(mdp, state_values):
    plt.figure(figsize=(3, 3))
    h, w = mdp.desc.shape
    states = sorted(mdp.get_all_states())
    V = np.array([state_values[s] for s in states])
    Pi = {s: get_optimal_action(mdp, state_values, s, gamma) for s in states}
    plt.imshow(V.reshape(w, h), cmap='gray', interpolation='none', clim=(0, 1))
    ax = plt.gca()
    ax.set_xticks(np.arange(h) - .5)
    ax.set_yticks(np.arange(w) - .5)
    ax.set_xticklabels([])
    ax.set_yticklabels([])
    Y, X = np.mgrid[0:4, 0:4]
    a2uv = {'left': (-1, 0), 'down': (0, -1), 'right': (1, 0), 'up': (-1, 0)}
    for y in range(h):
        for x in range(w):
            plt.text(x, y, str(mdp.desc[y, x].item()),
                     color='g', size=12, verticalalignment='center',
                     horizontalalignment='center', fontweight='bold')
            a = Pi[y, x]
            if a is None:
                continue
            u, v = a2uv[a]
            plt.arrow(x, y, u * .3, -v * .3, color='m', head_width=0.1, head_length=0.1)
    plt.grid(color='b', lw=2, ls='-')
    plt.show()

transition_probs = {
    's0': {
        'a0': {'s0': 0.5, 's2': 0.5},
        'a1': {'s2': 1}
    },
    's1': {
        'a0': {'s0': 0.7, 's1': 0.1, 's2': 0.2},
        'a1': {'s1': 0.95, 's2': 0.05}
    },
    's2': {
        'a0': {'s0': 0.4, 's1': 0.6},
        'a1': {'s0': 0.3, 's1': 0.3, 's2': 0.4}
    }
}
rewards = {
    's1': {'a0': {'s0': +5}},
    's2': {'a1': {'s0': -1}}
}

mdp = MDP(transition_probs, rewards, initial_state='s0')

# Start from state s0
print('initial state =', mdp.reset())
next_state, reward, done, info = mdp.step('a1')
print('next_state = %s, reward = %s, done = %s' % (next_state, reward, done))

# Display all available functions in the MDP class
print("mdp.get_all_states =", mdp.get_all_states())
print("mdp.get_possible_actions('s1') = ", mdp.get_possible_actions('s1'))
print("mdp.get_next_states('s1', 'a0') = ", mdp.get_next_states('s1', 'a0'))
print("mdp.get_reward('s1', 'a0', 's0') = ", mdp.get_reward('s1', 'a0', 's0'))
print("mdp.get_transition_prob('s1', 'a0', 's0') = ", mdp.get_transition_prob('s1', 'a0', 's0'))

# Value Iteration
# parameters
gamma = 0.9  # discount for MDP
num_iter = 100  # maximum iterations, excluding initialization
min_difference = 0.001  # stop VI if new values are this close to old values (or closer)

# initialize V(s)
state_values = {s: 0 for s in mdp.get_all_states()}

for i in range(num_iter):
    # Compute new state values using the functions defined above
    new_state_values = {state: get_new_state_value(mdp, state_values, state, gamma) for state, values in state_values.items()}

    # Compute difference of value
    diff = max(abs(new_state_values[s] - state_values[s]) for s in mdp.get_all_states())
    print("iter %4i   |   diff: %6.5f   |   " % (i, diff), end="")
    print('   '.join("V(%s) = %.3f" % (s, v) for s, v in state_values.items()), end='\n\n')
    state_values = new_state_values

    if diff < min_difference:
        print("Iteration Terminated")
        break

# Measure agent's average reward
s = mdp.reset()
rewards = []
for _ in range(10000):
    s, r, done, _ = mdp.step(get_optimal_action(mdp, state_values, s, gamma))
    rewards.append(r)
print("average reward: ", np.mean(rewards))

# Frozen lake
# initialize the environment
mdp = FrozenLakeEnv(slip_chance=0)
mdp.render()

state_values = value_iteration(mdp)

# Visualize the environment
print("\n Visualize in the environment: ")
s = mdp.reset()
mdp.render()
for t in range(100):
    a = get_optimal_action(mdp, state_values, s, gamma)
    print(a, end='\n\n')
    s, r, done, _ = mdp.step(a)
    mdp.render()
    if done:
        break

# Visualization
print("\nThe state action in each iteration")
state_values = {s: 0 for s in mdp.get_all_states()}

for i in range(10):
    print("after iteration %i" % i)
    state_values = value_iteration(mdp, state_values, num_iter=1)
    draw_policy(mdp, state_values)

print("\n30 Iteration over 8 * 8 grid")

mdp = FrozenLakeEnv(map_name='8x8', slip_chance=0.1)
state_values = {s: 0 for s in mdp.get_all_states()}

for i in range(30):
    state_values = value_iteration(mdp, state_values, num_iter=1)
    if i % 5 == 0:
        print("after iteration %i" % i)
        draw_policy(mdp, state_values)

print("\nMassive test")
mdp = FrozenLakeEnv(slip_chance=0)
state_values = value_iteration(mdp)

total_rewards = []
for game_i in range(1000):
    s = mdp.reset()
    rewards = []
    for t in range(100):
        s, r, done, _ = mdp.step(get_optimal_action(mdp, state_values, s, gamma))
        rewards.append(r)
        if done:
            break
    total_rewards.append(np.sum(rewards))

print("average reward: ", np.mean(total_rewards))
print("Iteration terminated")

print("\n Average reward of the agent")
mdp = FrozenLakeEnv(slip_chance=0.25)
state_values = value_iteration(mdp)

total_rewards = []
for game_i in range(1000):
    s = mdp.reset()
    rewards = []
    for t in range(100):
        s, r, done, _ = mdp.step(get_optimal_action(mdp, state_values, s, gamma))
        rewards.append(r)
        if done:
            break
    total_rewards.append(np.sum(rewards))

print("average reward: ", np.mean(total_rewards))
print("Iteration terminated")

max_stars_repo_name: gesuwen/Machine-Learning

max_issues_repo_name: gesuwen/Machine-Learning

max_forks_repo_name: gesuwen/Machine-Learning


text:
import pygslib
import numpy as np

def gcos(table, vardhole, varmodel, cutoff, model, cbb, decwt = None):
    """
    gcos(table, mask, vardhole, varmodel, cutoff, model, cbb, decwt = None)

    Produce global change of support validation (GCOS)

    Inputs
    ------
    table: pandas dataframe
        drillhole database
    vardhole: str
        variable on drillhole
    varmodel:
        variable in bmodel
    cutoff:
        cutoff list
    model: pandas dataframe
        block model
    cbb:
        block variance
    decwt: str default None
        declustering weight variable name

    """

    # firt the discrete gausian model
    zin = table[vardhole]

    if decwt is None:
        w = np.ones(table[vardhole].shape[0])
    else:
        w = table[decwt]

    zmin = table[vardhole].min(),
    zmax = table[vardhole].max(),

    PCI, H, raw, zana, gauss, z, P, raw_var, PCI_var, fig1 = pygslib.nonlinear.anamor(
       zin, w, ltail=1, utail=4, ltpar=1, utpar=1.5, K=30,
       ymin=-5, ymax=5, ndisc = 1000)

    r = pygslib.nonlinear.get_r(Var_Zv = cbb, PCI = PCI)

    ZV, PV, fig2 = pygslib.nonlinear.anamor_blk( PCI, H, r = r, gauss = gauss, Z = z,
                  ltail=1, utail=1, ltpar=1, utpar=1,
                  raw=raw, zana=zana)

    # calculate GTC
    tt = []
    gg = []
    label = []

    # calculate GTC from gaussian in block support

    cutoff[cutoff>ZV.max()] = ZV.max() - 0.000001 # there is an error of not like this
    cutoff[cutoff<ZV.min()] = ZV.min() + 0.000001
    t,ga,gb = pygslib.nonlinear.gtcurve (cutoff = cutoff, z=ZV, p=PV)

    tt.append(t)
    gg.append(ga)
    label.append(f'DGM {vardhole}')

    fig = pygslib.nonlinear.plotgt(cutoff = cutoff, t = tt, g = gg, label = label)

    # to compare global resources with the one estimated we calculate the CDF of the blocks

    # cdf of kriging estimate
    for v in varmodel:
        parameters_probplt = {
                'iwt'  : 0,
                'va'   : model[v][model[v].notnull()].values,
                'wt'   : np.ones(model[v][model[v].notnull()].shape[0])}

        binval_ok,cl_ok,xpt025,xlqt,xmed,xuqt,xpt975,xmin,xmax, \
                xcvr,xmen,xvar,error = pygslib.gslib.__plot.probplt(**parameters_probplt)

        # calculate GTC ok
        cutoff[cutoff>cl_ok.max()] = cl_ok.max() - 0.000001 # there is an error of not like this
        cutoff[cutoff<cl_ok.min()] = cl_ok.min() + 0.000001
        t,ga,gb = pygslib.nonlinear.gtcurve (cutoff = cutoff, z=cl_ok,
                       p=binval_ok, varred = 1, ivtyp = 2, zmin = 0, zmax = None,
                      ltail = 1, ltpar = 1, middle = 1, mpar = 1, utail = 1,
                      utpar = 1,maxdis = 1000)

        tt.append(t)
        gg.append(ga)
        label.append(f'Estimate {v}')

    return pygslib.nonlinear.plotgt(cutoff = cutoff, t = tt, g = gg, label = label)

max_stars_repo_name: cmrajan/pygslib

max_issues_repo_name: kaufmanno/pygslib

max_forks_repo_name: kaufmanno/pygslib


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

# [Eq. 11.1 & 11.2]
def _theta_c(Q, r, h):
    """
    :param Q: [W] Rate of heat release from the fire.
    :param r: [m] Radial distance from the centre of the fire plume impingement.
    :param h: [m]
    :return theta_c: [K]
    """
    Q = copy.copy(Q) / 1000.0  # Unit: [W] to [kW]

    # Condition for r/h
    if r / h > 0.18:
        theta_c = 5.38 * (Q / r) ** (2 / 3) / h
    else:
        theta_c = 16.9 * Q ** (2 / 3) / h ** (5 / 3)

    theta_c += 273.15  # Unit: [C] to [K]
    return theta_c

def _U(Q, r, h):
    """[Eq. 11.3 & 11.4]
    :param Q: [W] Rate of heat release from the fire.
    :param r: [m] Radial distance from the centre of the fire plume impingement.
    :param h: [m]
    :return U: [m/s]
    """
    Q = copy.copy(Q) / 1000.0  # Unit: [W] to [kW]

    if r / h > 0.15:
        U = 0.195 * Q ** (1 / 3) * h ** (1 / 2) / r ** (5 / 6)
    else:
        U = 0.96 * Q / h ** (1 / 3)

    return U

def _dT_d_dt(U, T_g, T_d, RTI):
    """[Eq. 11.5]
    :param U: [m/s]
    :param T_g: [K]
    :param T_d: [K]
    :param RTI: [-]
    :return dT_d_dt: [K/s]
    """

    dT_d_dt = U ** (1 / 2) * (T_g - T_d) / RTI

    return dT_d_dt

def _Q(alpha, t):
    """[Eq. 6.1]
    :param alpha:
    :param t:
    :return:
    """
    Q = alpha * t ** 2
    return Q

if __name__ == "__main__":
    import numpy as np

    # INPUTS
    time_start = 0
    time_step = 0.5
    time_end = 5 * 60

    alpha = 0.0117e3  # [W/s2]
    r = 2.75  # Estimation
    h = 2.6
    RTI = 115
    T_d_activation = 273.15 + 68  # [K]

    # CONTAINERS
    time = np.arange(time_start, time_end, time_step)
    Q = time * 0
    T_g = time * 0
    T_d = time * 0

    # INITIAL CONDITIONS
    T_g[0] = 273.15
    T_d[0] = 273.15

    # CALCULATION
    # calculate heat release rate
    Q = _Q(alpha, time)
    # calculate jet speed near sprinkler
    U = _U(Q, r, h)

    # calculate temperature near sprinkler
    T_g = _theta_c(Q, r, h)

    # calculate sprinkler temperature
    iter_time = enumerate(time)
    next(iter_time)
    for i, t in iter_time:
        T_d[i] = T_d[i - 1] + _dT_d_dt(U[i], T_g[i], T_d[i - 1], RTI)

    # RE-EVALUATE FIRE FOR SUPPRESSION ACTIVATION
    Q[T_d >= T_d_activation] = -1
    T_g[T_d >= T_d_activation] = -1
    T_d[T_d >= T_d_activation] = -1
    Q[Q == -1] = np.max(Q)
    T_d[T_d == -1] = np.max(T_d)
    T_g[T_g == -1] = np.max(T_g)

    print(
        "{:25}: {:7.2f} [min]".format(
            "Sprinkler activated at", np.min(time[T_d == np.max(T_d)]) / 60.0
        )
    )
    print("{:25}: {:7.2f} [kW]".format("The maximum HRR is", np.max(Q) / 1.0e3))

    import matplotlib.pyplot as plt

    plt.figure(num=1)
    plt.subplot(211)
    plt.plot(time / 60, Q / 1e6, label="HRR")
    plt.ylabel("Heat Release Rate [MW]")
    plt.subplot(212)
    plt.plot(time / 60, T_g - 273.15, label="Temperature (near field)")
    plt.plot(time / 60, T_d - 273.15, label="Temperature (sprinkler)")
    plt.xlabel("Time [min]")
    plt.ylabel("Temperature [$\degree C$]")
    plt.legend()
    plt.show()

max_stars_repo_name: fsepy/sfeprapy

max_issues_repo_name: fsepy/sfeprapy

max_forks_repo_name: fsepy/sfeprapy


text:
#!/usr/bin/env python

import sys, getopt
import glob, os
import numpy as np

help_message = 'usage example: python kmer_cluster_merge.py -r 1 -i /project/home/cluster_vectors/ -o /project/home/cluster_vectors/'
if __name__ == "__main__":
    try:
        opts, args = getopt.getopt(sys.argv[1:],'hr:i:o:',["filerank=","inputdir=","outputdir="])
    except:
        print help_message
        sys.exit(2)
    for opt, arg in opts:
        if opt in ('-h','--help'):
            print help_message
            sys.exit()
        elif opt in ('-i','--inputdir'):
            inputdir = arg
            if inputdir[-1] != '/':
                inputdir += '/'
        elif opt in ('-o','--outputdir'):
            outputdir = arg
            if outputdir[-1] != '/':
                outputdir += '/'
        elif opt in ('-r','--filerank'):
            fr = int(arg) - 1
    FP = glob.glob(os.path.join(inputdir+str(fr),'*.npy'))
    if len(FP) > 0:
        C = np.empty((10**9,),dtype=np.uint64)
        ci = 0
        for fp in FP:
            c = np.load(fp)
            C[ci:ci+c.shape[0]] = c
            ci += c.shape[0]
        np.save(outputdir+str(fr)+'.cluster.npy',C[:ci])
        os.system('rm -r %s%d/' % (outputdir,fr))
    else:
        print 'NO FILES:',fr

max_stars_repo_name: TaylorDennisLee/LatentStrainAnalysis

max_issues_repo_name: TaylorDennisLee/LatentStrainAnalysis

max_forks_repo_name: TaylorDennisLee/LatentStrainAnalysis


text:
import tensorflow as tf
k = tf.keras
import numpy as np
import cv2
from scipy.special import softmax, expit
import sys
import os
sys.path.insert(0, os.getcwd())
from retinaface_camera import get_anchors, nms_oneclass, decode_bbox, decode_landm

def detect_face(retinaface_model: k.Model,
                pfld_model: k.Model,
                anchors: np.ndarray,
                draw_img: np.ndarray,
                obj_thresh=0.7,
                nms_threshold=0.4,
                variances=[0.1, 0.2]):
  """ resize """
  img = cv2.cvtColor(draw_img, cv2.COLOR_BGR2RGB)
  """ normlize """
  det_img = ((img / 255. - 0.5) / 1)[None, ...]
  """ infer """
  predictions = retinaface_model.predict(det_img)
  """ parser """
  bbox, landm, clses = np.split(predictions[0], [4, -2], 1)
  """ softmax class"""
  clses = softmax(clses, -1)
  score = clses[:, 1]
  """ decode """
  bbox = decode_bbox(bbox, anchors, variances)
  bbox = bbox * np.tile([640, 640], [2])
  """ filter low score """
  inds = np.where(score > obj_thresh)[0]
  bbox = bbox[inds]
  score = score[inds]
  """ keep top-k before NMS """
  order = np.argsort(score)[::-1]
  bbox = bbox[order]
  score = score[order]
  """ do nms """
  keep = nms_oneclass(bbox, score, nms_threshold)

  for b, s in zip(bbox[keep].astype(int), score[keep]):
    cv2.rectangle(draw_img, tuple(b[:2]), tuple(b[2:]), (255, 0, 0), 2)
    cx, cy = (b[:2] + b[2:]) // 2
    halfw = np.max(b[2:] - b[:2]) // 2
    croped_img: np.ndarray = img[cy - halfw:cy + halfw, cx - halfw:cx + halfw]
    croped_wh = croped_img.shape[1::-1]
    if croped_wh[0] == croped_wh[1] and min(croped_wh) > 10:
      croped_img = cv2.resize(croped_img, (112, 112))
      croped_img = ((croped_img / 255. - 0.5) / 1)[None, ...]
      landmarks = pfld_model.predict(croped_img)
      s_point = np.array([cx - halfw, cy - halfw])
      for landm in np.reshape(expit(landmarks), (-1, 2)) * croped_wh:
        cv2.circle(draw_img, tuple((s_point + landm).astype(int)), 1, (0, 255, 0))

  return draw_img

if __name__ == "__main__":
  physical_devices = tf.config.experimental.list_physical_devices('GPU')
  assert len(physical_devices) > 0, "Not enough GPU hardware devices available"
  tf.config.experimental.set_memory_growth(physical_devices[0], True)
  retinaface_model: k.Model = k.models.load_model('asset/retinaface_train.h5')
  anchors = get_anchors([640, 640], [[0.025, 0.05], [0.1, 0.2], [0.4, 0.8]],
                        [8, 16, 32])
  pfld_model: k.Model = k.models.load_model('asset/pfld_infer.h5')

  capture = cv2.VideoCapture(0)
  capture.set(cv2.CAP_PROP_FRAME_WIDTH, 640)
  capture.set(cv2.CAP_PROP_FRAME_HEIGHT, 480)
  while (True):
    ret, img = capture.read()
    # img = cv2.flip(img, 1)
    img = cv2.copyMakeBorder(img, 80, 80, 0, 0, cv2.BORDER_CONSTANT, value=0)
    draw_img = detect_face(retinaface_model, pfld_model, anchors, img)
    cv2.imshow('frame', draw_img)
    if cv2.waitKey(1) == ord('q'):
      break

reason: import numpy,from scipy


text:
import argparse
import dill
import jax
import jax.numpy as np1
import numpy as np

import pomdp
pomdp.horizon = 25

parser = argparse.ArgumentParser()
parser.add_argument('--silent', action='store_true')
parser.add_argument('-i', type=int, default=0)
parser.add_argument('-n', type=int, default=5)
parser.add_argument('--cont', action='store_true')
args = parser.parse_args()

key = jax.random.PRNGKey(0)

###
S = 3
A = 2
Z = 12
with open('data/adni.obj', 'rb') as f:
    data = dill.load(f)
    data = data[:len(data)*args.i//args.n] + data[len(data)*(args.i+1)//args.n:]

n = len(data)
tau = max([d['tau'] for d in data])
data_a = -1 * np1.ones((n,tau), 'int')
data_z = -1 * np1.ones((n,tau), 'int')
for i, traj in zip(range(n), data):
    data_a = data_a.at[i,:traj['tau']].set(np1.array(traj['a']))
    data_z = data_z.at[i,:traj['tau']].set(np1.array(traj['z']))

def log_pi(alp, bet, b):
    res = np.zeros(A)
    for a in range(A):
        if alp[a].size == 0:
            res[a] = -1e6
        else:
            res[a] = bet * (alp[a] @ b).max()
    return res - np.log(np.sum(np.exp(res)))

def likelihood(b0, T, O, alp):
    res = 0
    for traj in data:
        b = b0
        for a, z in zip(traj['a'], traj['z']):
            res += log_pi(alp, 1, b)[a]
            b = O[a,:,z] * (T[:,a,:].T @ b)
            b /= b.sum()
    return res

def _messages_alp(alp_T_O, a_z):
    alp, T, O = alp_T_O
    a, z = a_z
    alp1 = jax.lax.select(a >= 0, O[a,:,z] * (T[:,a,:].T @ alp), alp)
    alp1 = alp1 / alp1.sum()
    return (alp1, T, O), alp1

def _messages_bet(bet_T_O, a_z):
    bet, T, O = bet_T_O
    a, z = a_z
    bet1 = jax.lax.select(a >= 0, T[:,a,:] @ (O[a,:,z] * bet), bet)
    bet1 = bet1 / bet1.sum()
    return (bet1, T, O), bet1

def messages(b0, T, O, traj_a, traj_z):
    _, alps = jax.lax.scan(_messages_alp, (b0, T, O), (traj_a, traj_z))
    alps = np1.concatenate((b0[None,...], alps))
    _, bets = jax.lax.scan(_messages_bet, (np1.ones(S), T, O), (traj_a, traj_z), reverse=True)
    bets = np1.concatenate((bets, np1.ones(S)[None,...]))
    gmms = alps * bets
    gmms = gmms / gmms.sum(axis=-1, keepdims=True)
    return gmms
# messages = jax.vmap(messages, in_axes=(None,None,None,0,0), out_axes=0)
# messages = jax.jit(messages)

def _posterior(s_alpT_alpO, a_z_gmm_key):
    s, alpT, alpO = s_alpT_alpO
    a, z, gmm, key = a_z_gmm_key
    s1 = jax.random.choice(key, range(S), p=gmm)
    alpT = jax.lax.select(a >= 0, alpT.at[s,a,s1].add(1), alpT)
    alpO = jax.lax.select(a >= 0, alpO.at[a,s1,z].add(1), alpO)
    return (s1, alpT, alpO), None

def posterior(b0, T, O, traj_a, traj_z, key):
    gmms = messages(b0, T, O, traj_a, traj_z)
    key, subkey = jax.random.split(key)
    s0 = jax.random.choice(subkey, range(S), p=gmms[0])
    alpb0 = np1.zeros(S)
    alpb0 = alpb0.at[s0].add(1)
    alpT = np1.zeros((S,A,S))
    alpO = np1.zeros((A,S,Z))
    subkeys = jax.random.split(key, tau)
    (_, alpT, alpO), _ = jax.lax.scan(_posterior, (s0, alpT, alpO), (traj_a, traj_z, gmms[1:], subkeys))
    return alpb0, alpT, alpO
posterior = jax.vmap(posterior, in_axes=(None,None,None,0,0,0), out_axes=(0,0,0))
posterior = jax.jit(posterior)

if args.cont:

    with open('res/adni/res{}-offpoirl.obj'.format(args.i), 'rb') as f:
        res = dill.load(f)
        key = res['key']
        b0, T, O, R = res['out'][-1]

else:

    res = dict()
    res['out'] = list()

    key, *subkey = jax.random.split(key, 4)
    b0 = np.array(jax.random.dirichlet(subkey[0], np.ones(S)))
    T = np.array(jax.random.dirichlet(subkey[1], np.ones((S,A,S)), shape=(S,A)))
    O = np.array(jax.random.dirichlet(subkey[1], np.ones((A,S,Z)), shape=(A,S)))

    ###
    O[0,:,(1,2,3,5,6,7,9,10,11)] = np.zeros((S,9)).T
    O[1,:,(0,4,8)] = np.zeros((S,3)).T
    O /= O.sum(axis=-1, keepdims=True)

    key, subkey = jax.random.split(key)
    R = np.zeros((S,A))
    R += .001 * np.array(jax.random.normal(subkey, shape=(S,A)))

alp = pomdp.solve(S, A, Z, b0, T, O, R)
like = likelihood(b0, T, O, alp)

rtio = 0
rtio_n = 0
for i in range(len(res['out']), 1000):

    _b0, _T, _O, _R = b0, T, O, R

    key, subkey = jax.random.split(key)
    if jax.random.choice(subkey, [True, False]):

        key, subkey = jax.random.split(key)
        subkeys = jax.random.split(subkey, n)
        alpb0, alpT, alpO = posterior(b0, T, O, data_a, data_z, subkeys)

        alpb0 = 1 + alpb0.sum(axis=0)
        alpT = 1 + alpT.sum(axis=0)
        alpO = 1 + alpO.sum(axis=0)

        key, subkey = jax.random.split(key)
        _b0 = np.array(jax.random.dirichlet(subkey, alpb0))
        _T = np.zeros((S,A,S))
        for s in range(S):
            for a in range(A):
                key, subkey = jax.random.split(key)
                _T[s,a,:] = np.array(jax.random.dirichlet(subkey, alpT[s,a,:]))
        _O = np.zeros((A,S,Z))
        for a in range(A):
            for s in range(S):
                key, subkey = jax.random.split(key)
                _O[a,s,:] = np.array(jax.random.dirichlet(subkey, alpO[a,s,:]))

    else:

        key, subkey = jax.random.split(key)
        _R = R + .001 * np.array(jax.random.normal(subkey, shape=(S,A)))

    _alp = pomdp.solve(S, A, Z, _b0, _T, _O, _R)
    _like = likelihood(_b0, _T, _O, _alp)

    key, subkey = jax.random.split(key)
    unif = jax.random.uniform(subkey)
    if np.log(unif) < _like - like:
        b0, T, O, R = _b0, _T, _O, _R
        like = _like

    rtio += 1 if like == _like else 0
    rtio_n += 1
    if not args.silent:
        print('i = {}, like = {}, {} ({})'.format(i, like, '*' if like == _like else '-', rtio / rtio_n))

    res['key'] = key
    res['out'].append((b0, T, O, R))
    if (i+1) % 100 == 0:
        with open('res/adni/res{}-offpoirl.obj'.format(args.i), 'wb') as f:
            dill.dump(res, f)

with open('res/adni/res{}-offpoirl.obj'.format(args.i), 'wb') as f:
    dill.dump(res, f)

max_stars_repo_name: DaraOrange/mlforhealthlabpub

max_issues_repo_name: DaraOrange/mlforhealthlabpub

max_forks_repo_name: DaraOrange/mlforhealthlabpub

reason: import numpy,import jax