text:
# -*- coding: utf-8 -*-
"""pointer_gen_modules.ipynb

Automatically generated by Colaboratory.

"""

import numpy as np
import random
import tensorflow as tf
import tensorflow.nn as nn

from utils import Linear
from utils import apply_mask_normalize
from utils import _mask_and_avg
from utils import _calc_final_dist

class Encoder():
  """ A simple encoder class to encode the input via Bidirectional LSTM
      Args:
          hpm : hyperparameters
          rand_unif_init : Initializer Object (random uniform) to initialize LSTMs parameters
          rand_norm_init : Initializer object (truncate normal) to initialize weights and biases for linear transf.
  """
  def __init__(self, hpm, rand_unif_init, rand_norm_init):
    self.hpm= hpm
    self.rand_unif_init = rand_unif_init
    self.rand_norm_init = rand_norm_init

    with tf.variable_scope('encoder'):
      self.lstm_cell_fw = tf.contrib.rnn.LSTMCell(self.hpm["hidden_size"],
                                           state_is_tuple= True, initializer=self.rand_unif_init) # forward lstm cell
      self.lstm_cell_bw = tf.contrib.rnn.LSTMCell(self.hpm["hidden_size"],
                                           state_is_tuple= True, initializer=self.rand_unif_init) # backward lstm cell

      self.w_c = Linear(self.hpm['hidden_size'], True, "reduce_c", self.rand_norm_init) # Parameters for the concatenated state linear transf.
      self.w_h = Linear(self.hpm['hidden_size'], True, 'reduce_h', self.rand_norm_init) # Parameters for the concatenated hidden output linear transf.

  def __call__(self, encoder_inputs, seq_lens):
    """ Call method for the encoding feedforward
        Args:
            encoder_inpputs : 3D tensor, shape : [batch_size, max_enc_len, embed_size]
            seq_lens : 1D tensor, lengths of the sequences (without padding) in the batch, shape : [batch_size]

        Returns:
            encoder_outputs : 3D tensor, output of the bidirectional dynamic rnn, shape : [batch_size, None, 2*hidden_size] (None because the max seq len vary with the batch)
            new state : tuple object made of two tensors : c => state, h=> last hidden output, shape : [2,batch_size, hidden_size]
    """
    with tf.variable_scope('encoder', reuse = tf.AUTO_REUSE):
      (encoder_outputs, (fw_st,bw_st)) = tf.nn.bidirectional_dynamic_rnn(
                                      self.lstm_cell_fw, self.lstm_cell_bw, encoder_inputs,
                                      dtype=tf.float32, swap_memory=True,
                                    sequence_length = seq_lens)

      encoder_outputs=tf.concat(encoder_outputs, axis= 2)

      old_c= tf.concat(values=[fw_st.c,bw_st.c], axis= 1) # we concatenate the forward and backward state, shape: [batch_size, 2*hidden_size]
      old_h= tf.concat(values=[fw_st.h,bw_st.h], axis= 1) # we concatenate the forwarrd and backward last hidden output, shape : [batch_size, 2*hidden_size]
      new_c= tf.nn.relu(self.w_c(old_c)) # linear transformation + relu activation, shape : [batch_size, hidden_size]
      new_h= tf.nn.relu(self.w_h(old_h)) # same as above

    return encoder_outputs, tf.contrib.rnn.LSTMStateTuple(new_c,new_h)

class Decoder():
  """
    A simple decoder class made of a unidirectional LSTM cell which decodes the next word given a previous one, a context vector and a previous state
    Args :
        hpm : hyperparameters
        rand_unif_init : Initializer Object (random uniform) to initialize LSTM parameters
  """
  def __init__(self,hpm,rand_unif_init):
    self.hpm= hpm
    self.rand_unif_init = rand_unif_init

    with tf.variable_scope('decoder'):
      self.lstm_cell= tf.contrib.rnn.LSTMCell(self.hpm["hidden_size"],
                                           state_is_tuple= True, initializer=self.rand_unif_init) # unidirectional lstm cell

  def __call__(self, dec_inputs, prev_state):
    """ Feedforward method for the simple decoder

        Args:
            dec_inputs : 2D tensor, list of words time step t for each sequence in the batch, shape = [batch_size, embed_size]
            prev_state : tuple object made of two vectors : c => state, h => last hidden output, shape : [2, batch_size, hidden_size]

        Returns:
            decoder_outputs : 2D tensor, shape = [batch_size, hidden_size]
            curr_st : current state of the decoder, shape : [2, batch_size, hidden_size]
    """
    with tf.variable_scope('decoder', reuse = tf.AUTO_REUSE):
    	decoder_outputs, curr_st= tf.nn.dynamic_rnn(self.lstm_cell, dec_inputs,
                                           dtype= tf.float32, initial_state=prev_state, swap_memory= True, time_major=True)
    return decoder_outputs, curr_st

class Attention_decoder():
  """
      An attentional based encoder-decoder model (bhadanau attention, additive style)
      Args:
          hpm : hyperparameters
          rand_unif_init : Initializer Object (random uniform) to initialize LSTMs parameters
          rand_norm_init : Initializer object (truncate normal) to initialize weights and biases for linear transf.

  """
  def __init__(self,hpm, rand_unif_init, rand_norm_init ):
    self.rand_unif_init = rand_unif_init
    self.rand_norm_init = rand_norm_init
    self.hpm=hpm

    with tf.variable_scope('attention_decoder', reuse = tf.AUTO_REUSE):
      self.decoder= Decoder(self.hpm, self.rand_unif_init) # simple decoder object (unidirecitional lstm)

      # Almost all the parameters (weights and biases) for the linear transformations (see below in the call method)

      self.w_h = Linear(self.hpm['attn_hidden_size'], True, "h")
      self.w_s = Linear(self.hpm['attn_hidden_size'], True, "s" )
      self.v = Linear(1, False, 'V')

      self.w_dec = Linear(self.hpm['emb_size'],True, "dec_inp")
      self.w_out = Linear(self.hpm['vocab_size'], True, 'out')

      if self.hpm['pointer_gen']:
        self.w_c_reduce = Linear(1, True, 'c_reduce')
        self.w_s_reduce = Linear(1, True, 's_reduce')
        self.w_i_reduce = Linear(1, True, 'i_reduce')

  def __call__(self, enc_outputs, enc_mask, enc_state, decoder_inputs,batch_max_oov_len = None, encoder_input_with_oov = None, cov_vec=None):
    """
        Attentional feedforward graph .
        We call this method once during training for each batch, and max_dec_len times for decode mode.

        Args:
            enc_outputs : 3D tensor, encoder outputs, shape : [batch_size, batch_max_enc_len, 2*hidden_size]
            enc_mask : 2D tensor, encoder sequence mask, shape : [batch_size, batch_max_enc_len]
            decoder_inputs: 3D tensor, decoder inputs, shape : [batch_size, max_dec_len, embed_size]
            batch_max_oov_len : Integer, Maximum number of oov for the current batch, (None if pointer_gen = False)
            encoder_input_with_oov : 2D tensor, encoder input with oovs ids, shape : [batch_size, batch_max_enc_len]

            !!! NB : batch_max_enc_len is None when we build graph, and vary during the feedforward with the current batch treated,
                      it is the maximum length of sequences of the current batch

        Returns : A dictionary
            output : list max_dec_en of 2D tensors of shape [batch_size, vocab_size + batch_max_oov_len (if pointer_gen)]
            last_context_vector : 2D tensor, shape : [batch_size, 2*hidden_size], this will be useful in the decode mode
            dec_state : 2D tensor, decoder last state, shape : [2, batch_size, hidden_size]
            p_gen : max_dec_len-many list of 1D tensors of length[batch_size] (only if pointer_gen is true)
            attention_vec : max_dec_len-many list of 2D tensors of shape [batch_size, batch_max_enc_len] (only if coverage is true)
    """

    if(self.hpm["pointer_gen"]):
      p_gens=[] # if pointer gen, we add an array to store the probability of each word in the sequences to be generated or pointed on

    attn_dists = [] # array to store the attention distributions over the enc seq
    dec_state = enc_state # we init the decoder state with the encoder last state
    outputs=[] # array to store the final probability distributions (decoded sequence)
    dec_inp = tf.unstack(decoder_inputs) # we unstack the decoder input to be able to enumerate over this tensor

    if  self.hpm['decode_using_prev']:
      argmax_arr = []
      samples_arr = []
      argmax_logprob_arr = []
      samples_logprob_arr = []

    # nested function
    def attention(dec_state, cov_vec=None):
      """
          Attention mechanism

          Args:
              dec_state : previous state of the decoder. shape : [2, batch_size, hidden_size]. For the first step, it corresponds to the encoder last state
              cov_vec : only if coverage is True (default None).  shape : [batch_size, <batch_max_enc_len>]. The previous coverage vector.

          Returns:
              attn_vec : 2D tensor, the attention vector at time step t. shape : [batch_size, <batch_max_enc_len>]
              context_vector : 2D tensor, shape: [batch_size, 2*hidden_size]
              cov_vec : 2D tensor, shape : [batch_size, <batch_max_enc_len>], the current coverage vector
      """
      if(self.hpm["coverage"]):
        with tf.variable_scope('coverage', reuse = tf.AUTO_REUSE ):
          w_c = tf.get_variable("w_c", [1,1,1,self.hpm['attn_hidden_size']]) # we add additional parameters for the coverage vector linear transf.

        cov_features = tf.expand_dims(tf.expand_dims(cov_vec, axis=2),axis=2) # given that the encoder max length is unknown and variable, we cannot just apply a
        cov_features = tf.nn.conv2d(cov_features, w_c, [1,1,1,1], "SAME")     # linear transformation as above. To avoid this issue, we can apply a convolution layer
                                                                              # which will transform the cov vector as a simple linear transf. would.

      # e = V*tanh(w_h*h + w_s*s + w_c*c ) (the last term, only is coverage = True)
      # attention weights all over the encoder input sequence
      # shape : [batch_size, <batch_max_enc_len>, 1]
        e=tf.nn.tanh(self.w_h(enc_outputs) +
                   tf.expand_dims(self.w_s(dec_state.c), axis=1) +
                   tf.squeeze(cov_features, [2]))
      else:
        e=tf.nn.tanh(self.w_h(enc_outputs) +
                   tf.expand_dims(self.w_s(dec_state.c), axis=1))
      e = self.v(e)

      # we take off the last dimension which equals 1
      e =  tf.reshape(e, [ e.get_shape().as_list()[0], -1]) # shape : [batch_size, <batch_max_enc_len>]

      attn_vec = tf.nn.softmax(e, axis=-1) # we apply a softmax on the attention weights to normalize them and obtain the attention vector.
      attn_vec = apply_mask_normalize(attn_vec, enc_mask) # Given that the input is padded with <PAD> token, the attentions weights over those tokens
                                                          # are not relevant, we apply the encoder input masks on the attention vectors to drop those 'irrelevant' attention weights
                                                          # and finally we re-normalize the attention weights to obtain probability distributions

      # context vector computation
      # we multiply the encoder outputs by the attention vector weigths (a weight for each output vector, when we consider only one sequence for the example)
      weighted_enc_outputs = tf.multiply(enc_outputs, tf.expand_dims(attn_vec, axis=-1)) # context vector at time step t, shape : [batch_size, ]
      context_vec = tf.reduce_sum(weighted_enc_outputs, axis=1)

      if self.hpm['coverage']:
          cov_vec = cov_vec + attn_vec # we update the coverage

      return attn_vec, context_vec, cov_vec
      # end of nested function

    with tf.variable_scope('attention_decoder', reuse = tf.AUTO_REUSE):
      # we compute the initial context vector
      _ , context_vec, _  = attention( dec_state, cov_vec)
      timesteps = self.hpm['max_dec_len']
      decoder_input = dec_inp[0]
      a=0
      if not self.hpm['decode_using_prev']:
        a = 1
      for i in range (a,timesteps):
        # for each item in the decoder inputs (this loops only once for decode mode)

        # concatenation of input (previous word) and context vector at timestep t
        new_dec_inp = tf.concat([decoder_input, context_vec], axis = -1) # shape : [batch_size, embed_size+2*hidden_size]
        new_dec_inp = self.w_dec(new_dec_inp) #shape : [batch_size, embed_size]

        # We apply the LSTM decoder on the new input
        dec_output, dec_state = self.decoder(tf.expand_dims(new_dec_inp, axis=0), dec_state) # dec_output shape : [1, batch_size, hidden_size]
                                                                                # dec_state shape : [2, batch_size, hidden_size] (2 for the state c and the last hidden output h)
        # attention vector of the current step, context vector for the next step
        # we update the coverage vector
        attn_vec, context_vec, cov_vec  = attention( dec_state, cov_vec)
        attn_dists.append(attn_vec)

        dec_output = tf.reshape(dec_output, [-1, dec_output.get_shape().as_list()[-1]]) # shape : [batch_size, hidden_size]
        dec_output = self.w_out(dec_output) # shape : [batch_size, vocab_size]
        vocab_dist = dec_output

        if not self.hpm['pointer_gen']:
          outputs.append(vocab_dist) # we do not apply yet the softmax function because this function is integrated in some futures ops like the loss function
        else:
          # if pointer_gen=True, we need to compute the softmax function because of the scatter op with the attention distribution
          outputs.append(tf.nn.softmax(dec_output, axis=-1))
          state = tf.concat([dec_state.c, dec_state.h], axis=1)

          #p_gen computation with the current concatenated state, context vector and the decoder input
          p_gen = tf.nn.sigmoid(self.w_c_reduce(context_vec)+
                             self.w_s_reduce(state )+
                             self.w_i_reduce(new_dec_inp)) # shape : [batch_size, 1]
          p_gens.append(p_gen)

      if self.hpm['pointer_gen']:
        # we apply the scatter op between the output distibutions (over the vocabulary) with the attention distributions
        outputs = _calc_final_dist(encoder_input_with_oov, outputs, attn_dists, p_gens, batch_max_oov_len, self.hpm)

      if not self.hpm['decode_using_prev']:
        decoder_input = dec_inp[i]
      else:

        batch_nums = tf.range(0, limit=self.hpm['batch_size'], dtype=tf.int64)
        argmax_seqs = []
        argmax_seqs_log_probs = []
        for i , x in enumerate(outputs):
          max_ids = tf.argmax(x, axis=-1)
          indices = tf.stack((batch_nums, max_ids), axis = -1)
          log_probs = tf.gather_nd(x, indices)
          argmax_seqs.append(max_ids)
          argmax_seqs_log_probs.append(log_probs)

        soft_outputs = tf.stack(outputs)
        if not self.hpm['pointer_gen']:
          soft_outputs = tf.softmax(soft_outputs)

        argmax_seqs = tf.stack(argmax_seqs)
        argmax_seqs_log_probs = tf.stack(argmax_seqs_log_probs)

        sampler = tf.distributions.Categorical(logits=soft_outputs)
        samples = sampler.sample()
        samples_log_probs = sampler.log_prob(samples)
        samples_log_probs = tf.identity(samples_log_probs)

        argmax_arr.append(argmax_seqs)
        argmax_logprob_arr.append(argmax_seqs_log_probs)
        samples_arr.append(samples)
        samples_logprob_arr.append(samples_log_probs)

        decoder_input = samples

    if  self.hpm['decode_using_prev']:
      argmax_arr = tf.stack(argmax_arr)
      argmax_logprob_arr = tf.stack(argmax_logprob_arr)
      samples_arr = tf.stack(samples_arr)
      samples_logprob_arr = tf.stack(samples_logprob_arr)

    dic = { 'output':outputs, 'last_context_vector':context_vec, 'dec_state':dec_state, 'attention_vec':attn_dists}
    if(self.hpm['pointer_gen']):
      dic['p_gen'] = p_gens
    if(self.hpm['coverage']):
      dic['coverage'] = cov_vec

    if  self.hpm['decode_using_prev']:
      dic.update({
        "argmax_seqs" : argmax_arr,
        "argmax_log_probs" : argmax_logprob_arr,
        "samples_seqs" : samples_arr,
        "samples_log_probs" : samples_logprob_arr
        })

    return dic

max_stars_repo_name: AIKevin/PGN

max_issues_repo_name: AIKevin/PGN

max_forks_repo_name: AIKevin/PGN


text:
# -*- coding: utf-8 -*-
"""
Compare with the CIFAR-10 Papers With Code pareto frontier here
https://paperswithcode.com/sota/image-classification-on-cifar-10?dimension=PARAMS
"""
import functools
import click
import numpy as np
import utils
import tqdm
import time
import torch
import torch.nn.functional as F
import torch.optim as optim
from torch.optim.lr_scheduler import OneCycleLR
from models import OriginalModel, ResidualNetwork

from torch.utils.tensorboard import SummaryWriter
from dataloaders import get_cifar10_dataloader

torch.backends.cudnn.benchmark = True

global_step = 0

def train_epoch(model, loader, optimizer, scheduler, tb_writer: SummaryWriter):
    global global_step
    epoch_start_time = time.time()

    epoch_training_num_correct = 0
    epoch_duration = time.time() - epoch_start_time
    n_images = 0
    for i, (images, labels) in enumerate(loader):
        n_images += images.shape[0]
        images = utils.to_cuda(images)
        labels = utils.to_cuda(labels)
        preds = model(images)
        loss = F.cross_entropy(preds, labels)

        optimizer.zero_grad()
        loss.backward()
        optimizer.step()
        scheduler.step()
        if i % 4 == 0: # Only log every 4th iteration (reduce I/O)
            tb_writer.add_scalar("loss", loss.detach().cpu().item(), global_step)
            tb_writer.add_scalar("learning_rate", scheduler.get_last_lr()[0], global_step)
        epoch_training_num_correct += (
            preds.detach().argmax(dim=1).eq(labels).sum().cpu().item())
        global_step += 1
    accuracy = epoch_training_num_correct / n_images
    tb_writer.add_scalar("accuracy", accuracy, global_step)
    print("train/accuracy: ", f"{accuracy:.4f}\t", end="", sep="")
    epoch_duration = time.time() - epoch_start_time
    return epoch_duration

@torch.no_grad()
def evaluate(model, dataloader, tb_writer):
    test_preds = torch.tensor([])
    n_images = 0
    epoch_testing_num_correct = 0
    epoch_testing_loss = 0

    for images, labels in dataloader:
        n_images += images.shape[0]
        images = utils.to_cuda(images)
        labels = utils.to_cuda(labels)
        preds = model(images)
        loss = F.cross_entropy(preds, labels)
        test_preds = torch.cat((test_preds, preds.cpu()), dim = 0)

        epoch_testing_loss += loss.item() * images.shape[0]
        epoch_testing_num_correct += (
            preds.argmax(dim=1).eq(labels).sum().cpu().item())

    testing_loss = epoch_testing_loss / n_images
    testing_accuracy = epoch_testing_num_correct / n_images
    print("test/loss: ", f"{testing_loss:.4f}\t ", "test/accuracy: ", f"{testing_accuracy:.4f}\t", end="", sep="")
    tb_writer.add_scalar("accuracy", testing_accuracy, global_step)
    tb_writer.add_scalar("loss", testing_loss, global_step)

model_zoo = dict(
    original=OriginalModel,
    revised=functools.partial(ResidualNetwork, start_ch=32, num_blocks_per_level=1, use_residual=False),
    revised_residual=functools.partial(ResidualNetwork, start_ch=32, num_blocks_per_level=1, use_residual=True),
    deeper=functools.partial(ResidualNetwork, start_ch=32, num_blocks_per_level=3, use_residual=True)
)
@click.command()
@click.argument("run_name", default="original")
@click.option("--batch-size", default=128, type=int)
@click.option("--n-epochs", default=100, type=int)
@click.option("--max-lr", default=0.1, type=float)
@click.option("--model", default="original", type=click.Choice(model_zoo.keys()))
def main(run_name: str, batch_size: int, n_epochs: int, max_lr: int, model: str):
    training_loader, testing_loader = get_cifar10_dataloader(batch_size)
    network = utils.to_cuda(model_zoo[model]())
    example_inputs = utils.to_cuda(torch.randn((batch_size, 3, 32, 32)))
    # Got a 30% runtime boost with scripting the module (on a NVIDIA 1060 GPU)
    print(network)
    network = torch.jit.script(network, example_inputs=(example_inputs,))
    num_params = sum([np.prod(p.shape) for p in network.parameters()])

    print("Number of parameters: ", num_params/10**6, "M",  sep="")

    optimizer = optim.Adam(network.parameters(), lr=max_lr)
    scheduler = OneCycleLR(
        optimizer,
        max_lr=max_lr,
        steps_per_epoch=len(training_loader),
        epochs=n_epochs)

    tb_train_writer = SummaryWriter(log_dir=f"outputs/{run_name}/train")
    tb_train_writer.add_scalar("num_params", num_params)
    tb_test_writer = SummaryWriter(log_dir=f"outputs/{run_name}/test")
    print("Starting train.")
    for i_epoch in range(n_epochs):
        print(f"epoch: {i_epoch}\t", end="")
        epoch_duration = train_epoch(network, training_loader, optimizer, scheduler, tb_train_writer)
        evaluate(network, testing_loader, tb_test_writer)
        print(f"epoch duration: {epoch_duration:.4f}")

main()

ext: py

reason: import numpy


text:
# coding = utf-8

import os
import numpy as np

from mykit.core.log import Verbose
from mykit.core.utils import (conv_string, get_filename_wo_ext, trim_after,
                              trim_both_sides)
from mykit.wien2k.constants import (IN1_UNIT_READER_nmr, IN1_UNIT_READER_v142,
                                    IN1_UNIT_READER_v171)
from mykit.wien2k.utils import find_complex_file, get_casename

class InputError(Exception):
    pass

class In1(Verbose):
    """class for manipulating WIEN2k in1 file

    TODO:
        Read data in and after KPOINT UNIT line

    Args:
        casename (str)
        switch (str)
        efermi (float) : Fermi energy
        rkmax (float) : RKmax
        lmax (int) : the maximum angular quantum of APW basis in MT region
        lnsmax (int)
        cmplx (bool)
        elparams (dict): the various information of exceptions of each atom
    """

    def __init__(self, casename, switch, efermi, rkmax, lmax, lnsmax, cmplx,
                 unit, emin, emax, nbands, *elparams, **kwargs):
        self.casename = casename
        self.switch = switch
        self.efermi = efermi
        self.rkmax = rkmax
        self.lmax = lmax
        self.lnsmax = lnsmax
        self.cmplx = cmplx
        self.unit = unit
        self.emin = emin
        self.emax = emax
        self.nbands = nbands
        self.elparams = elparams

    def add_exception(self, atomId, l, e, search=0.000, cont=True, apw=1):
        """Add one exception for a particular atom and l channel

        Args:
            atomId (int): index of inquivalent atom
            l (int): the angular momentum
            e (float): the linearization energy to set
            search (float): step size for energy search
            cont (bool)
            apw (0, 1)
        """
        # add to existing item only
        try:
            elparam = self.elparams[atomId]
        except IndexError:
            raise IndexError("atom index %d not available" % atomId)
        assert isinstance(cont, bool)
        assert apw in [0, 1]
        contStr = {True: "CONT", False: "STOP"}[cont]
        elparam["Ndiff"] += 1
        if not l in elparam["exceptions"]:
            elparam["exceptions"][l] = []
        elparam["exceptions"][l].append([e, search, contStr, apw])

    def remove_last_exception(self, atomId, l):
        """Remove the last linearization energy of atomId in channel l

        Args:
            atomId (int): index of inquivalent atom
            l (int): the angular momentum
        """
        try:
            elparam = self.elparams[atomId]
        except IndexError:
            raise IndexError("atom index %d not available" % atomId)
        if l in elparam["exceptions"]:
            elparam["Ndiff"] -= 1
            del elparam["exceptions"][l][-1]

    @property
    def natoms(self):
        """int. number of inequivalent atoms
        """
        return len(self.elparams)

    def get_exceptions(self, atomId, l=None, eonly=False):
        """Get exceptions of a particular atom

        Args:
            atomId (int)
            l (int): if specified, the exceptions in particular l channel is returned
            eonly (bool): get the linearization energy only

        Returns:
            dict if l is None. Otherwise list.
        """
        assert isinstance(atomId, int)
        try:
            elparam = self.elparams[atomId]["exceptions"]
            if eonly:
                es = {}
                for l in elparam.keys():
                    e = [x[0] for x in elparam[l]]
                    es[l] = e
                elparam = es
        except IndexError:
            return None
        if l is not None:
            elparam = elparam.get(l, None)
        return elparam

    def __str__(self):
        s = ["%5s  EF=%12.11f" % (self.switch, self.efermi), ]
        s += ["%4.1f %8d %4d" % (self.rkmax, self.lmax, self.lnsmax)]
        for _ia, elparams in enumerate(self.elparams):
            s.append(_write_el_block(elparams))
        last = "K-VECTORS FROM UNIT:%1d %6.1f" % (self.unit, self.emin)
        if self.emax > 100:
            last += "%10.1E" % self.emax
        else:
            last += "%10.1f" % self.emax
        if self.nbands is not None:
            last += "%6d" % self.nbands
        s.append(last)
        return '\n'.join(s)

    def write(self, pathIn=None, backup=False, suffix="_bak"):
        """Write to in1 file
        """
        path = pathIn
        if path is None:
            path = self.casename + 'in1' + 'c' * int(self.cmplx)
        if os.path.isdir(path):
            raise IOError("Trying to write in1 as an existing directory")
        if os.path.isfile(path) and backup:
            bakpath = path + suffix.strip()
            os.rename(path, bakpath)
        with open(path, 'w') as f:
            print(self.__str__(), file=f)

    @classmethod
    def read_from_file(cls, pathIn1=None):
        """Return In1 instance by reading an exisiting

        Args:
            filePath (str): the file name
        """
        cmplx = False
        if pathIn1 is None:
            casename = get_casename()
            path, cmplx = find_complex_file(casename, "in1")
        else:
            path = pathIn1
            if path.endswith("c"):
                cmplx = True
            casename = get_filename_wo_ext(path)

        with open(path, "r") as h:
            w2klines = h.readlines()
        switch = w2klines[0][:5]
        efermi = float(trim_both_sides(w2klines[0], r"=", r"\("))
        _params = trim_after(w2klines[1], r"\(").split()
        rkmax = float(_params[0])
        lmax, lnsmax = tuple(map(int, _params[1:3]))

        elparams = []
        i = 2
        _flagEnergy = True
        _flagInExcept = False
        while i < len(w2klines):
            line = trim_after(w2klines[i], r"\(")
            if line.startswith("K-VECTORS FROM UNIT"):
                _flagEnergy = False
                try:
                    # wien2k 17.1
                    data = IN1_UNIT_READER_v171.read(line)
                    unit, emin, emax, nbands = data
                    de = emax - emin
                except ValueError:
                    try:
                        # wien2k 14.2
                        data = IN1_UNIT_READER_v142.read(line)
                        unit, emin, de, nbands = data
                        emax = efermi + de
                    except ValueError:
                        # meet NMR in1
                        data = IN1_UNIT_READER_nmr.read(line)
                        unit, emin, emax = data
                        nbands = None
            else:
                words = line.split()
                if _flagEnergy and len(words) == 3:
                    ndiff = int(words[1])
                    atomEl = _read_el_block(w2klines[i : i + ndiff + 1])
                    elparams.append(atomEl)
                    i += ndiff
            i += 1
        return cls(casename, switch, efermi, rkmax, lmax, lnsmax, cmplx,
                   unit, emin, emax, nbands, *elparams)

def _read_el_block(elBlock):
    """
    Args:
        elBlock (list or tuple):  strings containing el information

    Returns
        dict
    """
    try:
        assert isinstance(elBlock, (list, tuple))
        for i in elBlock:
            assert isinstance(i, str)
    except AssertionError:
        raise InputError("elBlock should be a list of strings")
    elParams = {}
    globParams = elBlock[0].split()
    elParams["Etrial"] = float(globParams[0])
    ndiff = int(globParams[1])
    if len(elBlock) != ndiff + 1:
        raise InputError(
            "Inconsistent El block: need {}, parsed {}".format(ndiff, len(elBlock) - 1)
        )
    elParams["Ndiff"] = ndiff
    elParams["Napw"] = int(globParams[2])
    exceptions = {}
    for j in range(ndiff):
        # l, e, eIncre, cont, apw = EL_READER.read(elBlock[1 + j])
        s = elBlock[1+j]
        l, apw = conv_string(s, int, 0, -1)
        e, eIncre = conv_string(s, float, 1, 2)
        cont = s.split()[-2]
        if not l in exceptions:
            exceptions[l] = []
        exceptions[l].append([e, eIncre, cont, apw])
    elParams["exceptions"] = exceptions
    return elParams

def _write_el_block(elParams):
    """
    Args:
        elParams (dict)

    Returns
        str
    """
    ep = elParams
    ret = ["%6.2f %4d %2d" % (ep["Etrial"], ep["Ndiff"], ep["Napw"])]
    # sort by l
    exceptions = sorted(list(ep["exceptions"].items()), key=lambda x: x[0])
    for l, els in exceptions:
        for el in els:
            ret.append("%2d %9.5f %9.5f %4s %1d" % (l, *el))
    return "\n".join(ret)

max_stars_repo_name: minyez/mykit

max_issues_repo_name: minyez/mykit

max_forks_repo_name: minyez/mykit


text:
import numpy as np
import copy

from am_sim.model_parameters import high_en_exp_cutoff

detpop_color = 'C0'

def plot_exp_normalized_pop(ax, det_pf):
    '''
    This function plots the binding energy distribution of the population
    function on the specified axis. The distribution is normalized only for on
    the part below the high-energy experimental detection limit, so that it can
    be compared with experimental measurements.

    Args:
    - ax (matplotlib ax object): ax on which the binding energy distribution
        is plotted.
    - det_pf (det_pop object): det_pop whose binding energy distribution must
        be plotted.
    '''
    # domain and discretization step
    x, dx = det_pf.x, det_pf.dx
    # distribution
    y = det_pf.varphi
    # mask to select the range of energies lower than the experimental cutoff
    norm_mask = x <= high_en_exp_cutoff
    # normalization factor for the part below-cutoff
    norm = np.sum(y[norm_mask]) * dx
    # normalized distribution
    y_norm = y / norm
    # plotting function
    ax.plot(x, y_norm, color=detpop_color, label='theory')


text:
from typing import Tuple

import numpy as np
from numpy.typing import NDArray

def split_training_test(patterns: NDArray, targets: NDArray, ratio: float = 0.8) -> Tuple[
    NDArray, NDArray, NDArray, NDArray]:
    """
    Splits the entire dataset into a training set and a testing set
    :param ratio: How much of the entire set should be training
    :param patterns:
    :param targets:
    :return: training_patterns, training_targets, testing_patterns, testing_targets
    """
    mask = np.random.rand(len(patterns)) <= ratio
    return patterns[mask], targets[mask], patterns[~mask], targets[~mask]

max_stars_repo_name: theswampire/matura

max_issues_repo_name: theswampire/matura

max_forks_repo_name: theswampire/matura

reason: import numpy,from numpy


text:
#!/usr/bin/env python3

import os
import time
import datetime
from shutil import copyfile
import errno
import subprocess
import numpy as np
import numbers

from . import firefly

__all__ = ['PESScan', 'make_absolute_path', 'make_sure_path_exists']

def make_absolute_path(path):
    return os.path.normpath(os.path.join(os.getcwd(), path))

# creates path if not exists
def make_sure_path_exists(path):
    try:
        os.makedirs(path)
    except OSError as exception:
        if exception.errno != errno.EEXIST:
            raise
    return path

class PESScan(object):
    def __init__(self, molecule):
        self.molecule = molecule
        self.dry_run = False
        self.inputfile = self.punchfile = self.modes = self.coords = self.working_directory = None

    def modes_to_string(self, Qi=None):
        if Qi is None:
            return 'Q' + '_Q'.join([str(m) for m in self.modes])
        else:
            return 'Q' + '_Q'.join([str(m) + '_' + str('%.3f'%qi) for m, qi in zip(self.modes, Qi)])

    def perform_scan(self):
        if any(e is None for e in (self.inputfile, self.punchfile, self.modes, self.coords, self.working_directory)):
            raise ValueError('not all required parameters are set')

        os.chdir(self.working_directory)

        last_punchfile = self.punchfile
        for i, Q in enumerate(self.coords):
            # set molecule to coordinates for this step
            if isinstance(Q, (numbers.Integral, numbers.Real)):
                self.molecule.set_atom_displacement(self.modes, [Q]*len(self.modes))
                # create filenames
                this_step_inpfile = self.modes_to_string() + '_' + '%2.3f'%Q + '.inp'
                this_step_punchfile = self.modes_to_string() + '_' + '%2.3f'%Q + '.dat'
                this_step_outfile = self.modes_to_string() + '_' + '%2.3f'%Q + '.out'
            elif isinstance(Q, (list, np.ndarray)):
                Qlist = Q.tolist()
                self.molecule.set_atom_displacement(self.modes, Qlist)
                this_step_inpfile = self.modes_to_string(Qlist) + '.inp'
                this_step_punchfile = self.modes_to_string(Qlist) + '.dat'
                this_step_outfile = self.modes_to_string(Qlist) + '.out'
            else:
                raise ValueError('unsupported type of Q')

            # print('### Q = ', '%.3f'%Q, ' ###')
            print('###  Q =', Q, ' (' + str(i+1) + '/' + str(len(self.coords)) + ')   ###')
            print('\tInput file:', this_step_inpfile)

            copyfile(self.inputfile, this_step_inpfile)
            with open(this_step_inpfile, 'a') as f:
                f.write(self.molecule.firefly_coords())
                f.write(' $end\n')
                # dry run: no vec block to copy, skip
                if not self.dry_run:
                    vec = firefly.get_vec_block(last_punchfile)
                    f.writelines(vec)

            # dry run: skip rest of the loop, continue with next Q
            if self.dry_run:
                continue

            # include temporary files from previous steps, if they exist
            # tempfolders = ''
            # if os.path.exists('fftemp.0'):
            #     tempfolders = ' fftemp.*'

            t_start = datetime.datetime.now()
            print('\tStart:', t_start)
            # os.system('ffsub ' + this_step_inpfile) # + tempfolders)
            cp = subprocess.run(['ffsub', this_step_inpfile], stdout=subprocess.PIPE, universal_newlines=True)
            print('\tPBS_ID:', cp.stdout.strip('\n'))

            # wait for calculation to finish
            while not os.path.exists(this_step_punchfile):
                time.sleep(1)

            t_end = datetime.datetime.now()
            print('\tEnd:', t_end)
            print('\tUsed Walltime:', t_end-t_start)
            last_punchfile = this_step_punchfile

            # safety delay to wait for file transfer from node
            time.sleep(10)

    def perform_scan_with_basis(self):
        if any(e is None for e in (self.inputfile, self.punchfile, self.modes, self.coords, self.working_directory)):
            raise ValueError('not all required parameters are set')

        os.chdir(self.working_directory)

        last_punchfile = self.punchfile
        for Q in self.coords:
            # set molecule to coordinates for this step
            self.molecule.set_atom_displacement(self.modes, [Q]*len(self.modes))

            # create filenames
            this_step_inpfile = self.modes_to_string() + '_' + '%2.3f'%Q + '.inp'
            this_step_punchfile = self.modes_to_string() + '_' + '%2.3f'%Q + '.dat'
            this_step_outfile = self.modes_to_string() + '_' + '%2.3f'%Q + '.out'

            print('### Q = ', '%.3f'%Q, ' ###')
            print('\tInput file:', this_step_inpfile)

            copyfile(self.inputfile, this_step_inpfile)
            with open(this_step_inpfile, 'a') as f:
                f.write(self.molecule.firefly_coords_with_basis())
                f.write(' $end\n')
                # dry run: no vec block to copy, skip
                if not self.dry_run:
                    vec = firefly.get_vec_block(last_punchfile)
                    f.writelines(vec)

            # dry run: skip rest of the loop, continue with next Q
            if self.dry_run:
                continue

            # include temporary files from previous steps, if they exist
            # tempfolders = ''
            # if os.path.exists('fftemp.0'):
            #     tempfolders = ' fftemp.*'

            print('\tStart:', datetime.datetime.now())
            # os.system('ffsub ' + this_step_inpfile) # + tempfolders)
            cp = subprocess.run(['ffsub', this_step_inpfile], stdout=subprocess.PIPE, universal_newlines=True)
            print('\tPBS_ID:', cp.stdout.strip('\n'))

            # wait for calculation to finish
            while not os.path.exists(this_step_punchfile):
                time.sleep(1)

            print('\tEnd:', datetime.datetime.now())
            last_punchfile = this_step_punchfile

            # safety delay to wait for file transfer from node
            time.sleep(10)

max_stars_repo_name: oh-fv/vchamtools

max_issues_repo_name: oh-fv/vchamtools

max_forks_repo_name: oh-fv/vchamtools


text:
#!/usr/bin/env python -u

"""staticbyte.py - encode a stream of bytes as static noise

Usage:
    python -u staticbyte.py (-e | -d) <FILE> [-s <N>]
    python -u staticbyte.py -h

Options:
    -s  Set the chip size (default: 4096)
    -h  Print this help
"""

from __future__ import print_function

import numpy
import os
import sys
import wave

__author__ = 'Mansour Moufid'
__email__ = 'mansourmoufid@gmail.com'
__license__ = 'ISC'
__status__ = 'Development'
__version__ = '0.1'

PCM_DTYPE = numpy.int16

PCM_MAX = 2.0 ** (16 - 1) - 1.0

def readbytes(f):
    while True:
        char = f.read(1)
        if char == '':
            break
        byte = ord(char)
        yield byte

def xcor(x, y):
    x = numpy.copy(x)
    y = numpy.copy(y[::-1])
    n = x.size + y.size
    x.resize(n)
    y.resize(n)
    X = numpy.fft.rfft(x)
    Y = numpy.fft.rfft(y)
    Z = X * Y
    z = numpy.fft.irfft(Z, n=n)
    return z

def modxcor(x, y):
    z = xcor(x, y)
    n = y.size
    return z[(n // 2):-(n // 2)]

class Chip(object):

    def __init__(self, size, seed=None):
        self.size = size
        numpy.random.seed(seed)

    def codes(self, n):
        return numpy.float32(numpy.random.normal(0.0, 0.2, self.size * n))

if __name__ == '__main__':

    if len(sys.argv) == 2 and sys.argv[1] == '-h':
        print(__doc__)
        sys.exit(0)
    try:
        encode = '-e' in sys.argv
        i = sys.argv.index('-e' if encode else '-d')
        filename = sys.argv[i + 1]
        if '-s' in sys.argv:
            i = sys.argv.index('-s')
            size = int(sys.argv[i + 1])
        else:
            size = 2 ** 12
    except:
        print(__doc__)
        sys.exit(os.EX_USAGE)

    tty = os.ctermid()
    prompt = '\xF0\x9F\x94\x91  '
    with open(tty, 'w') as f:
        f.write(prompt)
    with open(tty, 'r') as f:
        line = f.readline()
    key = line.rstrip('\n')
    key = int(key)

    chip = Chip(size, seed=key)

    FS = 44100
    SAMPWIDTH = 2
    NCHANNELS = 1

    if encode:
        wav = wave.open(filename, 'wb')
        wav.setframerate(FS)
        wav.setsampwidth(SAMPWIDTH)
        wav.setnchannels(NCHANNELS)
    else:
        wav = wave.open(filename, 'rb')
        sampwidth = wav.getsampwidth()
        nchannels = wav.getnchannels()
        nframes = wav.getnframes()
        fs = wav.getframerate()
        assert fs == FS
        assert sampwidth == SAMPWIDTH
        assert nchannels == NCHANNELS

    if encode:
        for byte in readbytes(sys.stdin):
            codes = chip.codes(2 ** 8)
            i = byte * chip.size
            j = i + chip.size
            pcm = PCM_DTYPE(codes[i:j] * PCM_MAX)
            frames = pcm.tostring()
            wav.writeframes(frames)
    else:
        while True:
            frames = wav.readframes(chip.size)
            if len(frames) == 0:
                break
            pcm = numpy.fromstring(frames, dtype=PCM_DTYPE)
            data = numpy.float32(pcm)
            data = data / (PCM_MAX + 1.0)
            codes = chip.codes(2 ** 8)
            xc = modxcor(codes, data)
            i = numpy.argmax(xc ** 2)
            byte = i / chip.size
            if sys.stdout.isatty():
                byte &= 127
            sys.stdout.write(chr(byte))

    wav.close()

max_stars_repo_name: eliteraspberries/staticbit

max_issues_repo_name: eliteraspberries/staticbit

max_forks_repo_name: eliteraspberries/staticbit


text:
#!/usr/bin/env python

import argparse
import os
import os.path as osp
import scipy.misc

import rospkg

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument('--overwrite', action='store_true')
    parser.add_argument('-d', '--dir', default='v5')
    args = parser.parse_args()

    # answer candidates
    # default is no
    yes = ['yes', 'y']
    no = ['no', 'n', '']

    rospack = rospkg.RosPack()
    datasetdir = osp.join(
        rospack.get_path('selective_dualarm_stowing'), 'dataset', args.dir)
    for i, dirname in enumerate(os.listdir(datasetdir)):
        datadir = osp.join(datasetdir, dirname, 'after_stow')
        print('count: {}'.format(i+1))
        print(datadir)
        afterdirs = os.listdir(datadir)
        afterdirs = [x for x in afterdirs if osp.isdir(osp.join(datadir, x))]
        assert len(afterdirs) == 1
        afterdir = osp.join(datadir, afterdirs[0])

        if not args.overwrite and osp.exists(osp.join(afterdir, 'label.txt')):
            continue

        # is_bimanual
        with open(osp.join(afterdir, 'is_bimanual.txt'), 'r') as f:
            data = f.read()
            if data == 'True':
                is_bimanual = True
            else:
                is_bimanual = False
        if is_bimanual:
            style = 'dualarm'
        else:
            style = 'singlearm'

        imgpath = osp.join(afterdir, 'default_camera_image.png')
        img = scipy.misc.imread(imgpath)
        scipy.misc.imshow(img)

        # ask dropped
        while True:
            print('Dropped? [y/N]')
            choice = raw_input().lower()
            if choice in yes:
                dropped = True
                break
            elif choice in no:
                dropped = False
                break
            else:
                print('Type yes/no')

        # ask protruded
        while True:
            print('protruded? [y/N]')
            choice = raw_input().lower()
            if choice in yes:
                protruded = True
                break
            elif choice in no:
                protruded = False
                break
            else:
                print('Type yes/no')

        # ask damaged
        while True:
            print('damaged? [y/N]')
            choice = raw_input().lower()
            if choice in yes:
                damaged = True
                break
            elif choice in no:
                damaged = False
                break
            else:
                print('Type yes/no')

        label = []
        if is_bimanual:
            style = 'dualarm'
        else:
            style = 'singlearm'

        if not dropped and not protruded:
            label.append('success')
        else:
            if dropped:
                label.append('drop')
            if protruded:
                label.append('protrude')
            if damaged:
                label.append('damage')

        label = ['{0}_{1}'.format(style, x) for x in label]
        labeltxt = '\n'.join(label)

        labelpath = osp.join(afterdir, 'label.txt')
        with open(labelpath, 'w+') as f:
            f.write(labeltxt)

if __name__ == '__main__':
    main()

reason: import scipy


text:
# -*- coding: utf-8 -*-
import os
import cv2
import imageio
from scipy import misc
from PIL import Image
import numpy as np

# --- alterado -------------------
print(os.system('pwd'))
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
# --------------------------------

dataset_dir = 'datasets' # datasets
dataset = 'data_crop_512_jpg'
train_dir = 'train'
test_dir = 'val_test'
val_dir = 'val'
resize_to = 256
scale = 4

if not os.path.exists(os.path.join(dataset_dir, train_dir)):
    os.mkdir(os.path.join(dataset_dir, train_dir))
if not os.path.exists(os.path.join(dataset_dir, test_dir)):
    os.mkdir(os.path.join(dataset_dir, test_dir))
if not os.path.exists(os.path.join(dataset_dir, val_dir)):
    os.mkdir(os.path.join(dataset_dir, val_dir))

i = 1
fileList = os.listdir(os.path.join(dataset_dir, dataset))
for index, file in enumerate(fileList):
    imgPath = os.path.join(dataset_dir, dataset, file)
    if os.path.isdir(imgPath):
        continue
    if(i > 1000):
        i = i + 1
        raise Exception('Forced finish')
    print("procesing " + file + " " + str(i+1) + '/' + str(len(fileList)))

    #img = cv2.imread(imgPath)
    #print("DBG - img sizes 1", img.shape)
    #imgplot = plt.imshow(img)
    #plt.show()

    img = cv2.cvtColor(cv2.imread(imgPath), cv2.COLOR_BGR2RGB)
    #print("DBG - img sizes 2", img.shape)
    #img = misc.imresize(img, (resize_to, resize_to), interp='bilinear')
    img = np.array(Image.fromarray(img).resize((resize_to, resize_to), Image.BILINEAR))
    #print("DBG - img sizes 3", img.shape)
    size = int(img.shape[0] / scale)
    #print("DBG - size", size)
    #resizeImg = misc.imresize(img, (size, size), interp='bilinear')
    resizeImg = np.array(Image.fromarray(img).resize((size, size), Image.BILINEAR))
    #print("DBG - img sizes 4", resizeImg.shape)
    #finalImg = misc.imresize(resizeImg, (img.shape[0], img.shape[0]), interp='bilinear')
    finalImg = np.array(Image.fromarray(resizeImg).resize((img.shape[0], img.shape[0]), Image.BILINEAR))
    #print("DBG - img sizes 5", finalImg.shape)
    combineImg = Image.new('RGB', (img.shape[0]*2, img.shape[0]))
    combineImg.paste(Image.fromarray(finalImg), (0,0))
    combineImg.paste(Image.fromarray(img), (img.shape[0]+1,0))
    savePath = ""
    # sample ratio of train:test:val is 6:2:2
    if index % 10 < 6:
        savePath = os.path.join(dataset_dir, train_dir, file)
    elif index % 10 < 8:
        savePath = os.path.join(dataset_dir, test_dir, file)
    else:
        savePath = os.path.join(dataset_dir, val_dir, file)
    misc.imsave(savePath, combineImg)
    imageio.imwrite(savePath, combineImg)

max_stars_repo_name: fialhocoelho/projetoic

max_issues_repo_name: fialhocoelho/projetoic

max_forks_repo_name: fialhocoelho/projetoic

reason: import numpy,from scipy


text:
import argparse

import numpy as np
from caffe2.python import core, workspace

def benchmark_concat(num_inputs, input_dim, axis, add_axis, iterations):
    input_names = [f"input{i}" for i in range(num_inputs)]
    for n in input_names:
        workspace.FeedBlob(n, np.random.randn(*input_dim).astype(np.float32))

    net = core.Net("benchmark_net")
    net.Concat(input_names, ["output", "split_info"], axis=axis, add_axis=add_axis)
    workspace.CreateNet(net)

    runtimes = workspace.BenchmarkNet(net.Name(), 1, iterations, True)
    print(f"{num_inputs * np.prod(input_dim) * 4 / runtimes[1] / 1e6} GB/s")

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="minimal benchmark for concat.")
    parser.add_argument("--num_inputs", type=int, default=2)
    parser.add_argument("--input_dim", nargs="+", type=int, required=True)
    parser.add_argument("--axis", type=int, default=-1)
    parser.add_argument("--add_axis", type=int, default=0)
    parser.add_argument("--iterations", type=int, default=64)
    args, extra_args = parser.parse_known_args()
    core.GlobalInit(["python"] + extra_args)
    benchmark_concat(
        args.num_inputs, args.input_dim, args.axis, args.add_axis, args.iterations
    )

max_stars_repo_name: Hacky-DH/pytorch

max_issues_repo_name: Hacky-DH/pytorch

max_forks_repo_name: Hacky-DH/pytorch


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Date    : 2018-11-28 16:27:19
# @Author  : Bolun Wang (bolunwang@cs.ucsb.edu)
# @Link    : http://cs.ucsb.edu/~bolunwang

import os
import sys
import time

import numpy as np
from keras.preprocessing import image

##############################
#        PARAMETERS          #
##############################

RESULT_DIR = 'results'  # directory for storing results
IMG_FILENAME_TEMPLATE = 'gtsrb_visualize_%s_label_%d.png'  # image filename template for visualization results

# input size
IMG_ROWS = 32
IMG_COLS = 32
IMG_COLOR = 3
INPUT_SHAPE = (IMG_ROWS, IMG_COLS, IMG_COLOR)

NUM_CLASSES = 43  # total number of classes in the model

##############################
#      END PARAMETERS        #
##############################

def outlier_detection(l1_norm_list, idx_mapping):

    consistency_constant = 1.4826  # if normal distribution
    median = np.median(l1_norm_list)
    mad = consistency_constant * np.median(np.abs(l1_norm_list - median))
    min_mad = np.abs(np.min(l1_norm_list) - median) / mad

    print('median: %f, MAD: %f' % (median, mad))
    print('anomaly index: %f' % min_mad)

    flag_list = []
    for y_label in idx_mapping:
        if l1_norm_list[idx_mapping[y_label]] > median:
            continue
        if np.abs(l1_norm_list[idx_mapping[y_label]] - median) / mad > 2:
            flag_list.append((y_label, l1_norm_list[idx_mapping[y_label]]))

    if len(flag_list) > 0:
        flag_list = sorted(flag_list, key=lambda x: x[1])

    print('flagged label list: %s' %
          ', '.join(['%d: %2f' % (y_label, l_norm)
                     for y_label, l_norm in flag_list]))

    pass

def analyze_pattern_norm_dist():

    mask_flatten = []
    idx_mapping = {}

    for y_label in range(NUM_CLASSES):
        mask_filename = IMG_FILENAME_TEMPLATE % ('mask', y_label)
        if os.path.isfile('%s/%s' % (RESULT_DIR, mask_filename)):
            img = image.load_img(
                '%s/%s' % (RESULT_DIR, mask_filename),
                color_mode='grayscale',
                target_size=INPUT_SHAPE)
            mask = image.img_to_array(img)
            mask /= 255
            mask = mask[:, :, 0]

            mask_flatten.append(mask.flatten())

            idx_mapping[y_label] = len(mask_flatten) - 1

    l1_norm_list = [np.sum(np.abs(m)) for m in mask_flatten]

    print('%d labels found' % len(l1_norm_list))

    outlier_detection(l1_norm_list, idx_mapping)

    pass

if __name__ == '__main__':

    print('%s start' % sys.argv[0])

    start_time = time.time()
    analyze_pattern_norm_dist()
    elapsed_time = time.time() - start_time
    print('elapsed time %.2f s' % elapsed_time)

max_stars_repo_name: ByzanTine/backdoor

max_issues_repo_name: ByzanTine/backdoor

max_forks_repo_name: ByzanTine/backdoor


text:
import unittest
import subprocess
import os

import pandas as pd
import numpy as np

import sys
test_dir_1 = os.path.dirname(os.path.realpath(__file__))
sys.path.append(test_dir_1)

from base import TestBase

core_stats = \
""",GSM1606107,GSM1606108
total_reads,1000000,1000000
clipped_reads,966503,961142
clipped_reads_%,96.65,96.11
filtered_out,831157,765609
filtered_out_%,86.0,79.66
filter_kept,135346,195533
filter_kept_%,14.0,20.34
transcriptome_aligned_once,103931,150356
transcriptome_aligned_once_%,76.79,76.9
transcriptome_aligned_many,15757,23560
transcriptome_aligned_many_%,11.64,12.05
transcriptome_total_aligned,119688,173916
transcriptome_total_aligned_%,88.43,88.94
transcriptome_unaligned,15658,21617
transcriptome_unaligned_%,11.57,11.06
qpass_aligned_reads,108390,157151
qpass_aligned_reads_%,90.56,90.36
after_dedup,104742,148386
after_dedup_%,96.63,94.42"""

genome_stats = \
""",GSM1606107,GSM1606108
genome_aligned_once,2.0,8.0
genome_aligned_once_%,0.01,0.04
genome_aligned_many,0.0,0.0
genome_aligned_many_%,0.0,0.0
genome_total_aligned,2.0,8.0
genome_total_aligned_%,0.01,0.04
genome_unaligned,15656.0,21609.0
genome_unaligned_%,99.99,99.96"""

post_genome_stats = \
""",GSM1606107,GSM1606108
post_genome_aligned_once,0.0,0.0
post_genome_aligned_once_%,0.0,0.0
post_genome_aligned_many,0.0,0.0
post_genome_aligned_many_%,0.0,0.0
post_genome_total_aligned,0.0,0.0
post_genome_total_aligned_%,0.0,0.0
post_genome_unaligned,15656.0,21609.0
post_genome_unaligned_%,100.0,100.0"""

class TestConcat(TestBase):

    def setUp(self):
        self.core_stats_file        = "core_stats.csv"
        self.genome_stats_file      = "genome_stats.csv"
        self.post_genome_stats_file = "postgenome_stats.csv"
        self.output_file            = "merged_stats.csv"

        self.files = [self.core_stats_file,        self.genome_stats_file,
                      self.post_genome_stats_file, self.output_file]

        with open(self.core_stats_file, "w" ) as output_stream:
            print(core_stats, file = output_stream)

        with open(self.genome_stats_file, "w" ) as output_stream:
            print(genome_stats, file = output_stream)

        with open(self.post_genome_stats_file, "w" ) as output_stream:
            print(post_genome_stats, file = output_stream)

    def test_concat_1(self):
        command = ["rfc", "merge", "concat-csv",
                   "--out", self.output_file,
                   self.core_stats_file,
                   self.genome_stats_file,
                   self.post_genome_stats_file]

        output, error = self.run_command(command)

        if error:
            print("Script Error:\n", error)

        result_df = pd.read_csv(self.output_file, index_col = 0)

        self.assertTrue(np.allclose(result_df.loc["post_genome_unaligned"], [15656.0, 21609.0 ]))

if __name__ == '__main__':

    unittest.main()

max_stars_repo_name: ribosomeprofiling/RFCommands

max_issues_repo_name: ribosomeprofiling/RFCommands

max_forks_repo_name: ribosomeprofiling/RFCommands


text:
import numpy as np
import pandas as pd
import time
from collections import defaultdict
from argparse import ArgumentParser

class GridWorld:
    def __init__(self, path):
        self.path = path  # loaded file path
        self.grid = np.loadtxt(self.path, dtype='str')  # initial grid
        self.n_row, self.n_col = np.shape(self.grid)  # # of row, # of colum
        self.start_state = list()  # start sate
        self.barrier_list = list()  # list of barrier state
        self.end_state = list()  # lsit of end state
        self.possible_state = list()  # list of reachable states
        self.act_n = ['up', 'down', 'right', 'left']  # name list of possible actions
        self.sta_n = 0  # number of possible states
        self._Q = None  # initial Q table
        self.view_dic = defaultdict(int)  # dictionary of number of reach times in each square during training

    def init_func(self):
        for i in range(self.n_row):
            for j in range(self.n_col):
                value = self.grid[i][j]
                if value == 'X':
                    self.barrier_list.append((i, j))
                else:
                    if value == 'S':
                        self.start_state.append((i, j))
                        self.possible_state.append((i, j))
                    elif value == '0':
                        self.possible_state.append((i, j))
                    else:
                        self.end_state.append((i, j))
        self.sta_n = len(self.possible_state)
        Q_table = np.zeros((self.sta_n, len(self.act_n)))
        self._Q = pd.DataFrame(Q_table, index=self.possible_state, columns=self.act_n, dtype='float')
        for i in self.possible_state:
            self.view_dic[i] = 0

class RLearningAgent(GridWorld):
    def __init__(self, path, learning_rate, gamma, suc_p, reward, time, epsilon=0.9):
        super(RLearningAgent, self).__init__(path)
        self._lr = learning_rate  # how fast the algorithm to converge
        self._gamma = gamma  # discounting rate
        self._sp = suc_p  # action success probability
        self.reward = reward  # reward for each action takes
        self._rt = time  # trainning time
        self._ep = epsilon  # epsilon value, probability of choose max Q in action determination
        self.score = 0  # score each trail get
        self.total_score = 0  # total score of training
        self.n_round = 0  # number of training trails

    def get_max_Q_direction(self, state):
        data = self._Q.loc[[state], :]
        max_value = data.max(axis=1)[0]
        result = data.isin([max_value])
        col_list = list(result.any()[result.any() == True].index)
        index = np.random.randint(0, len(col_list))
        direction = col_list[index]
        return direction

    def determineAction(self, state):
        rand_num = np.random.uniform(0, 1)
        if rand_num < self._ep:
            direction = self.get_max_Q_direction(state)  # direction of highest Q(s,a) value
        else:
            index = np.random.randint(0, len(self.act_n))
            direction = self.act_n[index]  # direction of highest Q(s,a) value
        return direction

    def takeAction_p2(self, state, direction):
        rand_num = np.random.uniform(0, 1)
        if rand_num <= self._sp:  # successful, forward 1 squre
            steps = 1
        elif rand_num >= (1 - (1 - self._sp) / 2):  # random action 2, backward 1 squre
            if direction == 'up':
                direction = 'down'
                steps = 1
            elif direction == 'down':
                direction = 'up'
                steps = 1
            elif direction == 'left':
                direction = 'right'
                steps = 1
            else:
                direction = 'left'
                steps = 1
        else:  # random action 1, jump 2 squres
            steps = 2
        row = state[0]
        col = state[1]
        if steps == 1:
            if direction == 'up':
                next_row = row - steps
                next_col = col
                next_state = (next_row, next_col)
                if (next_row, next_col) in self.barrier_list or next_row < 0:  # reach barrier or boundary
                    next_state = (row, col)  # stay
            elif direction == 'down':
                next_row = row + steps
                next_col = col
                next_state = (next_row, next_col)
                if (next_row, next_col) in self.barrier_list or next_row >= self.n_row:  # reach barrier
                    next_state = (row, col)  # stay
            elif direction == 'right':
                next_row = row
                next_col = col + steps
                next_state = (next_row, next_col)
                if (next_row, next_col) in self.barrier_list or next_col >= self.n_col:  # reach barrier
                    next_state = (row, col)  # stay
            else:
                next_row = row
                next_col = col - steps
                next_state = (next_row, next_col)
                if (next_row, next_col) in self.barrier_list or next_col < 0:  # reach barrier
                    next_state = (row, col)  # stay
        else:
            if direction == 'up':
                next_row = row - steps
                next_col = col
                next_state = (next_row, next_col)
                if (next_row + 1, next_col) in self.barrier_list or next_row + 1 < 0:  # reach barrier or boundary
                    next_state = (row, col)  # stay
                elif (next_row, next_col) in self.barrier_list or next_row < 0:  # 1 squre from barrier or boundary
                    next_state = (next_row + 1, col)  # move up 1 square
            elif direction == 'down':
                next_row = row + steps
                next_col = col
                next_state = (next_row, next_col)
                if (
                next_row - 1, next_col) in self.barrier_list or next_row - 1 >= self.n_row:  # reach barrier or boundary
                    next_state = (row, col)  # stay
                elif (
                next_row, next_col) in self.barrier_list or next_row >= self.n_row:  # 1 squre from barrier or boundary
                    next_state = (next_row - 1, next_col)
            elif direction == 'right':
                next_row = row
                next_col = col + steps
                next_state = (next_row, next_col)
                if (
                next_row, next_col - 1) in self.barrier_list or next_col - 1 >= self.n_col:  # reach barrier or boundary
                    next_state = (row, col)  # stay
                elif (
                next_row, next_col) in self.barrier_list or next_col >= self.n_col:  # 1 squre from barrier or boundary
                    next_state = (next_row, next_col - 1)
            else:
                next_row = row
                next_col = col - steps
                next_state = (next_row, next_col)
                if (next_row, next_col + 1) in self.barrier_list or next_col + 1 < 0:  # reach barrier or boundary
                    next_state = (row, col)  # stay
                elif (next_row, next_col) in self.barrier_list or next_col < 0:  # 1 squre from barrier or boundary
                    next_state = (next_row, next_col + 1)
        return next_state

    def notTerminal(self, state):
        if state in self.end_state:
            return False
        else:
            return True

    def QL_update(self, state, direction, next_state):
        predict_Q = self._Q.at[state, direction]
        if self.notTerminal(next_state):
            next_direction = self.get_max_Q_direction(next_state)  # direction of highest Q(s,a) value in next state
            target_Q = self.reward + self._gamma * self._Q.at[next_state, next_direction]
            self.score += self.reward
            self.view_dic[next_state] += 1
        else:
            target_Q = int(self.grid[next_state[0]][next_state[1]])
            self.score += int(self.grid[next_state[0]][next_state[1]])
        self._Q.at[state, direction] += self._lr * (target_Q - predict_Q)

    def QL_learning_p2(self):
        st = time.time()
        timeRemains = 0
        while timeRemains < self._rt:
            self.score = 0
            self.n_round += 1
            index = np.random.randint(0, len(self.start_state))
            state = self.start_state[index]
            self.view_dic[state] += 1
            while self.notTerminal(state):
                direction = self.determineAction(state)
                next_state = self.takeAction_p2(state, direction)
                self.QL_update(state, direction, next_state)
                state = next_state
            # print('the score of {}th round ------ {}'.format(self.n_round, self.score))  # can be deleted
            et = time.time()
            timeRemains = et - st
            self.total_score += self.score
        return self.total_score / self.n_round

    def SARSA_update(self, state, direction, next_state, next_direction):
        predict_Q = self._Q.at[state, direction]
        if self.notTerminal(next_state):
            target_Q = self.reward + self._gamma * self._Q.at[next_state, next_direction]
            self.score += self.reward
            self.view_dic[next_state] += 1
        else:
            target_Q = int(self.grid[next_state[0]][next_state[1]])
            self.score += int(self.grid[next_state[0]][next_state[1]])
        self._Q.at[state, direction] += self._lr * (target_Q - predict_Q)

    def SARSA_learning(self):
        st = time.time()
        timeRemains = 0
        while timeRemains < self._rt:
            self.score = 0
            self.n_round += 1
            index = np.random.randint(0, len(self.start_state))
            state = self.start_state[index]
            self.view_dic[state] += 1
            while self.notTerminal(state):
                direction = self.determineAction(state)
                next_state = self.takeAction_p2(state, direction)
                if self.notTerminal(next_state):
                    next_direction = self.determineAction(next_state)
                else:
                    next_direction = ''
                self.SARSA_update(state, direction, next_state, next_direction)
                state = next_state
                direction = next_direction
            # print('the score of {}th round ------ {}'.format(self.n_round, self.score)) # can be deleted
            et = time.time()
            timeRemains = et - st
            self.total_score += self.score
        return self.total_score / self.n_round

    def draw_heat_grid_Q(self):
        index_list = self._Q.index
        heat_grid = np.zeros((self.n_row, self.n_col))
        d = pd.DataFrame(0, index=index_list, columns=['Q_max'], dtype='float')
        ter_sum = 0
        for i in index_list:
            max_value = self._Q.loc[[i], :].max(axis=1)[0]
            d.loc[[i], 'Q_max'] = abs(max_value)
        max_Q = d.max()[0]
        d = d / max_Q
        total_Q = d.sum()[0]
        for i in self.end_state:
            row = i[0]
            col = i[1]
            value = int(self.grid[row][col])
            ter_sum += value
            heat_grid[row][col] = value
        for i in self.barrier_list:
            heat_grid[i[0]][i[1]] = float('-inf')
        for i in index_list:
            row = i[0]
            col = i[1]
            heat_grid[row][col] = d.loc[[i], 'Q_max'] / total_Q * (100 - ter_sum)
        return heat_grid

    def draw_heat_grid_V(self):
        index_list = list(self.view_dic.keys())
        heat_grid = np.zeros((self.n_row, self.n_col))
        view_pd = pd.DataFrame.from_dict(self.view_dic, orient='index', columns=['num'])
        total_view = view_pd.sum()[0]
        ter_sum = 0
        for i in self.end_state:
            row = i[0]
            col = i[1]
            value = int(self.grid[row][col])
            ter_sum += value
            heat_grid[row][col] = value
        for i in self.barrier_list:
            heat_grid[i[0]][i[1]] = float('-inf')
        for i in index_list:
            row = i[0]
            col = i[1]
            value = self.view_dic[i] / total_view * (100 - ter_sum)
            heat_grid[row][col] = value
        return heat_grid

    def draw_policy_grid(self):
        index_list = self._Q.index
        policy_grid = self.grid[:]
        for i in index_list:
            row = i[0]
            col = i[1]
            direction = self._Q.loc[[i], :].idxmax(axis=1)[0]
            if direction == 'right':
                policy_grid[row][col] = '>'
            elif direction == 'left':
                policy_grid[row][col] = '<'
            elif direction == 'up':
                policy_grid[row][col] = '^'
            else:
                policy_grid[row][col] = 'v'
        return policy_grid

if __name__ == "__main__":
    parser = ArgumentParser()
    parser.add_argument("-f", dest="file_path", help="the file path for the grid with postfix",
                        default='grid_sample.txt')
    parser.add_argument("-lr", dest="learning_rate", help="learning rate of RL", type=float,
                        default=0.1)
    parser.add_argument("-g", dest="gamma", help="discounting rate", type=float,
                        default=0.9)
    parser.add_argument("-s", dest="action_success", help="action success probability", type=float,
                        default=1.0)
    parser.add_argument("-r", dest="reward", help="RL rewards", type=float,
                        default=-0.04)
    parser.add_argument("-t", dest="run_time", help="training time", type=float,
                        default=10.0)
    parser.add_argument("-e", dest="epsilon", help="epsilon greedy exploration rate", type=float,
                        default=0.9)

    parsed_args = parser.parse_args()

    path = parsed_args.file_path # grid file path
    lr = parsed_args.learning_rate # learning rate
    dr = parsed_args.gamma  # discounting rate
    sr = parsed_args.action_success # action success probability
    rw = parsed_args.reward # reward for each action takes
    run_t = parsed_args.run_time # training time
    epsilon = parsed_args.epsilon # epsilon value, probability of choose max Q in action determination

    # temp = GridWorld(path)
    # temp.init_func()
    # test0 = temp._Q
    # print(test0)

    start_time = time.time()

    LearningClass = RLearningAgent(path, lr, dr, sr, rw, run_t, epsilon)
    LearningClass.init_func()
    original_mtx = LearningClass.grid
    print('The original board is')
    print(original_mtx)

    mean_rw = LearningClass.QL_learning_p2()
    updated_mtx = LearningClass._Q

    heat_mtx = LearningClass.draw_heat_grid_V()
    policy_mtx = LearningClass.draw_policy_grid()

    # print('The total run time is %s' %(time.time()-start_time))

    # print('The updated Q table is:')
    # print(updated_mtx)

    print('The heat map after training is:')
    print(heat_mtx)
    print('The policy after training is:')
    print(policy_mtx)

    print("------------------------------------------")
    # print("The total reward is: %s" % (q_temp.total_score))
    print("The number of times for reaching a terminal is: %s" % (LearningClass.n_round))
    print("The mean reward per trail is: %s" % (mean_rw))


text:
"""
Field interpolation example using equivalent surface currents
==============================================================

"""
import numpy as np
from bfieldtools.mesh_conductor import MeshConductor, StreamFunction
from mayavi import mlab
import trimesh

import mne

PLOT = True
SAVE_FIGURES = True
IMPORT_MNE_DATA = True

SAVE_MNE_DATA = True
SAVE_DIR = "./MNE interpolation/"

from pyface.api import GUI

_gui = GUI()

#%%
# First, let's import the MEG data

if IMPORT_MNE_DATA:

    from mne.datasets import sample

    data_path = sample.data_path()
    fname = data_path + "/MEG/sample/sample_audvis-ave.fif"
    # Reading
    condition = "Left Auditory"
    evoked = mne.read_evokeds(fname, condition=condition, baseline=(None, 0), proj=True)
    evoked.pick_types(meg="mag")
    # evoked.plot(exclude=[], time_unit="s")

    i0, i1 = evoked.time_as_index(0.08)[0], evoked.time_as_index(0.09)[0]
    field = evoked.data[:, i0:i1].mean(axis=1)

    # Read BEM for surface geometry and transform to correct coordinate system
    import os.path as op

    subject = "sample"
    subjects_dir = op.join(data_path, "subjects")
    bem_fname = op.join(
        subjects_dir, subject, "bem", subject + "-5120-5120-5120-bem-sol.fif"
    )
    bem = mne.read_bem_solution(bem_fname)

    # Head mesh 0
    # Innerskull mesh 2
    surf_index = 2

    trans_fname = op.join(data_path, "MEG", "sample", "sample_audvis_raw-trans.fif")
    trans0 = mne.read_trans(trans_fname)
    R = trans0["trans"][:3, :3]
    t = trans0["trans"][:3, 3]
    # Surface from MRI to HEAD
    rr = (bem["surfs"][surf_index]["rr"] - t) @ R
    # Surface from HEAD to DEVICE
    trans1 = evoked.info["dev_head_t"]
    R = trans1["trans"][:3, :3]
    t = trans1["trans"][:3, 3]
    rr = (rr - t) @ R

    innerskull = trimesh.Trimesh(rr, bem["surfs"][surf_index]["tris"])

    surf_index = 0

    R = trans0["trans"][:3, :3]
    t = trans0["trans"][:3, 3]
    # Surface from MRI to HEAD
    rr = (bem["surfs"][surf_index]["rr"] - t) @ R
    # Surface from HEAD to DEVICE
    R = trans1["trans"][:3, :3]
    t = trans1["trans"][:3, 3]
    rr = (rr - t) @ R
    head = trimesh.Trimesh(rr, bem["surfs"][surf_index]["tris"])

    mesh = head

    # Sensor locations and directions in DEVICE coordinate system
    p = np.array(
        [
            ch["loc"][:3]
            for ch in evoked.info["chs"]
            if ch["ch_name"][-1] == "1" and ch["ch_name"][:3] == "MEG"
        ]
    )
    n = np.array(
        [
            ch["loc"][-3:]
            for ch in evoked.info["chs"]
            if ch["ch_name"][-1] == "1" and ch["ch_name"][:3] == "MEG"
        ]
    )

    if PLOT:
        # Plot sensor locations and directions
        fig = mlab.figure(bgcolor=(1, 1, 1))
        mlab.triangular_mesh(*innerskull.vertices.T, innerskull.faces)
        mlab.triangular_mesh(
            *head.vertices.T, head.faces, color=(0.5, 0.5, 0.5), opacity=0.5
        )
        mlab.quiver3d(*p.T, *n.T, mode="arrow")

        fig.scene.isometric_view()

        if SAVE_FIGURES:
            mlab.savefig(SAVE_DIR + "MEG_geometry.png", magnification=4, figure=fig)

    if SAVE_MNE_DATA:
        np.savez(
            SAVE_DIR + "mne_data.npz",
            mesh=head,
            p=p,
            n=n,
            vertices=head.vertices,
            faces=head.faces,
        )
        evoked.save(SAVE_DIR + "left_auditory-ave.fif")

else:

    with np.load(SAVE_DIR + "mne_data.npz", allow_pickle=True) as data:
        mesh = data["mesh"]
        p = data["p"]
        n = data["n"]
        mesh = trimesh.Trimesh(vertices=data["vertices"], faces=data["faces"])

    evoked = mne.Evoked(SAVE_DIR + "left_auditory-ave.fif")

#%%
# Fit the surface current for the auditory evoked response

c = MeshConductor(mesh_obj=mesh, basis_name="suh", N_suh=150)
M = c.mass

sensor_coupling = np.einsum("ijk,ij->ik", c.B_coupling(p), n)
# a = np.linalg.pinv(sensor_coupling, rcond=1e-15) @ field
ss = np.linalg.svd(sensor_coupling @ sensor_coupling.T, False, False)

# reg_exps = [0.5, 1, 2, 3, 4, 5, 6, 7, 8]
reg_exps = [1]
rel_errors = []

for reg_exp in reg_exps:
    _lambda = np.max(ss) * (10 ** (-reg_exp))
    # Laplacian in the suh basis is diagonal
    BB = sensor_coupling.T @ sensor_coupling + _lambda * (-c.laplacian) / np.max(
        abs(c.laplacian)
    )
    a = np.linalg.solve(BB, sensor_coupling.T @ field)
    s = StreamFunction(a, c)
    b_filt = sensor_coupling @ s

    rel_error = np.linalg.norm(b_filt - field) / np.linalg.norm(field)
    print("Relative error:", rel_error * 100, "%")
    rel_errors.append(rel_error)

    if PLOT:
        fig = mlab.figure(bgcolor=(1, 1, 1))
        surf = s.plot(False, figure=fig)
        surf.actor.mapper.interpolate_scalars_before_mapping = True
        surf.module_manager.scalar_lut_manager.number_of_colors = 16

        if SAVE_FIGURES:
            mlab.savefig(
                SAVE_DIR + "SUH_scalp_streamfunction.png", magnification=4, figure=fig
            )

#%%
# Interpolate MEG data to the sensor surface

from bfieldtools.utils import load_example_mesh

helmet = load_example_mesh("meg_helmet", process=False)
# Bring the surface roughly to the correct place
helmet.vertices[:, 2] -= 0.05

# Reset coupling by hand
c.B_coupling.reset()
mlab.figure(bgcolor=(1, 1, 1))
B_surf = np.sum(
    c.B_coupling(helmet.vertices) * helmet.vertex_normals[:, :, None], axis=1
)

if PLOT:
    fig = mlab.quiver3d(*p.T, *n.T, mode="arrow")
    scalars = B_surf @ s
    surf = mlab.triangular_mesh(
        *helmet.vertices.T, helmet.faces, scalars=scalars, colormap="seismic"
    )
    surf.actor.mapper.interpolate_scalars_before_mapping = True
    surf.module_manager.scalar_lut_manager.number_of_colors = 15

    if SAVE_FIGURES:
        mlab.savefig(
            SAVE_DIR + "SUH_sensors_streamfunction.png", magnification=4, figure=fig
        )

#%%
# Calculate magnetic field in volumetric grid

Nvol = 30
x = np.linspace(-0.125, 0.125, Nvol)
vol_points = np.array(np.meshgrid(x, x, x, indexing="ij")).reshape(3, -1).T
# mlab.points3d(*vol_points.T)

c.B_coupling.reset()
Bvol_coupling = c.B_coupling(vol_points, Nchunks=100, analytic=True)
s = StreamFunction(a, c)
# s = StreamFunction(a, c)
Bvol = Bvol_coupling @ s

#%%
# Now, plot the computed magnetic field using streamlines

if PLOT:
    from bfieldtools.mesh_calculus import gradient

    fig = mlab.figure(bgcolor=(1, 1, 1))

    surf = s.plot(False, figure=fig)
    surf.actor.mapper.interpolate_scalars_before_mapping = True
    surf.module_manager.scalar_lut_manager.number_of_colors = 16

    vecs = mlab.pipeline.vector_field(
        *vol_points.T.reshape(3, Nvol, Nvol, Nvol), *Bvol.T.reshape(3, Nvol, Nvol, Nvol)
    )
    vecnorm = mlab.pipeline.extract_vector_norm(vecs)

    seed_points = mesh.vertices[mesh.faces].mean(axis=1) - 0.01 * mesh.face_normals

    seed_vals = c.basis @ c.inductance @ s
    seed_vals_grad = np.linalg.norm(gradient(seed_vals, c.mesh), axis=0)

    seed_vals = abs(seed_vals[mesh.faces].mean(axis=1)) ** 2
    seed_vals[seed_vals_grad > seed_vals_grad.max() / 1.8] = 0

    Npoints = 500
    seed_inds = np.random.choice(
        np.arange(len(seed_vals)), Npoints, False, seed_vals / seed_vals.sum()
    )

    seed_points = seed_points[seed_inds]

    streams = []

    for pi in seed_points:
        streamline = mlab.pipeline.streamline(
            vecnorm,
            integration_direction="both",
            colormap="BuGn",
            seed_visible=False,
            seedtype="point",
        )
        streamline.seed.widget.position = pi
        streamline.stream_tracer.terminal_speed = 3e-13
        streamline.stream_tracer.maximum_propagation = 0.1
        streamline.actor.property.render_lines_as_tubes = True
        streamline.actor.property.line_width = 4.0
        streams.append(streamline)

    # Custom colormap with alpha channel
    streamine = streams[0]
    lut = streamline.module_manager.scalar_lut_manager.lut.table.to_array()
    lut[:, -1] = np.linspace(0, 255, 256)
    streamline.module_manager.scalar_lut_manager.lut.table = lut
    streamline.module_manager.scalar_lut_manager.data_range = np.array(
        [1.0e-13, 1.0e-12]
    )

    for streamline in streams:
        streamline.stream_tracer.terminal_speed = 1e-13
        streamline.seed.widget.hot_spot_size = 0.1
        streamline.stream_tracer.initial_integration_step = 0.01
        streamline.stream_tracer.minimum_integration_step = 0.1

    sensors = mlab.quiver3d(*p.T, *n.T, mode="cylinder")
    sensors.glyph.glyph_source.glyph_source.height = 0.1
    sensors.actor.property.color = (0.5, 0.5, 0.5)
    sensors.actor.mapper.scalar_visibility = False
    sensors.glyph.glyph_source.glyph_source.resolution = 32
    sensors.glyph.glyph.scale_factor = 0.03
    # sensors.glyph.glyph_source.glyph_source.shaft_radius = 0.05

    fig.scene.camera.position = [
        0.637392177469018,
        0.07644693029292644,
        -0.07183513804689762,
    ]
    fig.scene.camera.focal_point = [
        -6.413459777832031e-05,
        0.01716560870409012,
        -0.0229007127850005,
    ]
    fig.scene.camera.view_angle = 30.0
    fig.scene.camera.view_up = [
        0.04390624852005244,
        0.3114421192517664,
        0.9492502555685007,
    ]
    fig.scene.camera.clipping_range = [0.3366362817578398, 1.0281065506557443]
    fig.scene.camera.compute_view_plane_normal()

    while fig.scene.light_manager is None:
        _gui.process_events()

    camera_light = fig.scene.light_manager.lights[0]
    camera_light.intensity = 0.7

    if SAVE_FIGURES:
        mlab.savefig(
            SAVE_DIR + "SUH_streamlines_lateral.png", figure=fig, magnification=4
        )

    fig.scene.camera.position = [
        -6.413459777832031e-05,
        0.01716560870409012,
        0.6191735842078244,
    ]
    fig.scene.camera.focal_point = [
        -6.413459777832031e-05,
        0.01716560870409012,
        -0.0229007127850005,
    ]
    fig.scene.camera.view_angle = 30.0
    fig.scene.camera.view_up = [0.0, 1.0, 0.0]
    fig.scene.camera.clipping_range = [0.3381552363433513, 1.0261944997830243]
    fig.scene.camera.compute_view_plane_normal()

    if SAVE_FIGURES:
        mlab.savefig(
            SAVE_DIR + "SUH_streamlines_coronal.png", figure=fig, magnification=4
        )

#%%
# Finally, plot the surface current density itself

fig = mlab.figure(bgcolor=(1, 1, 1))
grad_s = gradient(c.basis @ s, mesh, rotated=True)
q = mlab.quiver3d(
    *(mesh.vertices[mesh.faces].mean(axis=1).T),
    *grad_s,
    colormap="viridis",
    mode="arrow"
)

mlab.triangular_mesh(*head.vertices.T, head.faces, color=(0.8, 0.8, 0.8), opacity=1.0)

fig.scene.camera.position = [
    0.4987072212753703,
    0.06469079487766746,
    -0.0014732384935239248,
]
fig.scene.camera.focal_point = [
    0.0018187984824180603,
    0.012344694641686624,
    -0.04367139294433087,
]
fig.scene.camera.view_angle = 30.0
fig.scene.camera.view_up = [
    -0.10720122151366927,
    0.23975383168819672,
    0.9648968848000314,
]
fig.scene.camera.clipping_range = [0.28329092545021717, 0.7772019991936254]
fig.scene.camera.compute_view_plane_normal()

if SAVE_FIGURES:
    mlab.savefig(
        SAVE_DIR + "SUH_surface_currents_lateral.png", figure=fig, magnification=4
    )

max_stars_repo_name: SAFedorov/bfieldtools

max_issues_repo_name: SAFedorov/bfieldtools

max_forks_repo_name: SAFedorov/bfieldtools


text:
# -*- coding: utf-8 -*-
"""Copy of inspect_balloon_model.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1MLAg2B1HwiLTDBG8xdXnPir66v8KTjft

# Mask R-CNN - Inspect Ballon Trained Model

Code and visualizations to test, debug, and evaluate the Mask R-CNN model.

## Setup
"""

import os
import sys
import skimage.io
import matplotlib.pyplot as plt
import numpy as np
import tensorflow as tf
import cv2

# Root directory of the project
ROOT_DIR = os.path.abspath("../../")
ROOT_DIR = os.path.join(ROOT_DIR,"Mask_RCNN")
#how many frames you want to take(intervall of 1 is every frame)
INTERVALL = 10

# Import Mask RCNN
sys.path.append(ROOT_DIR)  # To find local version of the library
from mrcnn import visualize
from mrcnn.visualize import display_images
import mrcnn.model as modellib
from mrcnn.model import log
import mrcnn.utils as utils

from samples.balloon import balloon

# %matplotlib inline

# Directory to save logs and trained model
MODEL_DIR = os.path.join(ROOT_DIR, "logs")

# Path to Ballon trained weights
# You can download this file from the Releases page
# https://github.com/matterport/Mask_RCNN/releases
# This path will differ depending on where you saved the Weights and how they are named
BALLON_WEIGHTS_PATH = os.path.join(ROOT_DIR, "data","model", "mask_rcnn_tool_disassembly_0030.h5")

"""## Configurations"""

config = balloon.BalloonConfig()
BALLOON_DIR = os.path.join(ROOT_DIR, "data","gear")

# Override the training configurations with a few
# changes for inferencing.
class InferenceConfig(config.__class__):
    # Run detection on one image at a time
    GPU_COUNT = 1
    IMAGES_PER_GPU = 1

config = InferenceConfig()
config.display()

"""## Notebook Preferences"""

# Device to load the neural network on.
# Useful if you're training a model on the same
# machine, in which case use CPU and leave the
# GPU for training.
DEVICE = "/cpu:0"  # /cpu:0 or /gpu:0

# Inspect the model in training or inference modes
# values: 'inference' or 'training'
# TODO: code for 'training' test mode not ready yet
TEST_MODE = "inference"

def get_ax(rows=1, cols=1, size=16):
    """Return a Matplotlib Axes array to be used in
    all visualizations in the notebook. Provide a
    central point to control graph sizes.

    Adjust the size attribute to control how big to render images
    """
    _, ax = plt.subplots(rows, cols, figsize=(size * cols, size * rows))
    return ax

"""## Load Validation Dataset"""

# Load validation dataset
dataset = balloon.BalloonDataset()
dataset.load_balloon(BALLOON_DIR, "val")

# Must call before using the dataset
dataset.prepare()

"""## Load Model"""

# Create model in inference mode
with tf.device(DEVICE):
    model = modellib.MaskRCNN(mode="inference", model_dir=MODEL_DIR,
                              config=config)

# Run this cell to mount your Google Drive.
# click on the link and copy the password

# copy the weights from your drive to the logs directory

# Set path to balloon weights file

# Download file from the Releases page and set its path
# https://github.com/matterport/Mask_RCNN/releases
# weights_path = "/path/to/mask_rcnn_balloon.h5"

# Or, load the last model you trained
weights_path = BALLON_WEIGHTS_PATH
images_path = os.path.join(ROOT_DIR,"data","video")

# Load weights
print("Loading weights ", weights_path)
model.load_weights(weights_path, by_name=True)

"""## Run Detection"""
cam = cv2.VideoCapture(os.path.join(images_path,"video1.mp4"))
result = []
counter = 0
while(True):
    ret, image = cam.read()

    if not ret:
        break
    if (not counter%INTERVALL==0):
        counter+=1
        continue
    # Run object detection
    results = model.detect([image], verbose=1)
    # Display results
    ax = get_ax(1)
    r = results[0]
    visualize.display_instances(image, r['rois'], r['masks'], r['class_ids'],
                               dataset.class_names, r['scores'], ax=ax,
                                title="Predictions")

    for index,roi in enumerate(r['rois']):
        result.append([counter//INTERVALL,np.array(roi),r['class_ids'][index]])

    counter+=1

#comment this line if you do not want to inspect the results
plt.show()

file = open("allAOIs","w")

for counter, aoi in enumerate(result):

    text = str(aoi[0]) + " " + str(aoi[1][0]) + " " + str(aoi[1][1]) + " " + str(aoi[1][2]) + " " + str(aoi[1][3]) + " " + dataset.class_names[aoi[2]] + "\n"
    file.write(text)
file.close()


text:
import numpy as np
import matplotlib.pyplot as plt
import h5py
import scipy.io
import math
import sklearn
import sklearn.datasets

def init_adam(parameters) :

    L = len(parameters) // 2 # number of layers
    v = {}  # the exponentially weighted average of the gradient
    s = {}  # the exponentially weighted average of the squared gradient

    for l in range(L):
        v["dW" + str(l+1)] = np.zeros((parameters["W" + str(l+1)].shape[0], parameters["W" + str(l+1)].shape[1]))
        v["db" + str(l+1)] = np.zeros((parameters["b" + str(l+1)].shape[0], parameters["b" + str(l+1)].shape[1]))
        s["dW" + str(l+1)] = np.zeros((parameters["W" + str(l+1)].shape[0], parameters["W" + str(l+1)].shape[1]))
        s["db" + str(l+1)] = np.zeros((parameters["b" + str(l+1)].shape[0], parameters["b" + str(l+1)].shape[1]))

    return v, s

def update_params_adam(parameters, grads, v, s, t, learning_rate = 0.01,
                                beta1 = 0.9, beta2 = 0.999,  epsilon = 1e-8):

    L = len(parameters) // 2
    v_corrected = {}   # first moment estimate
    s_corrected = {}   # second moment estimate

    for l in range(L):
        # Moving average of the gradients
        v["dW" + str(l+1)] = beta1*v["dW" + str(l+1)] + (1 - beta1)*grads['dW' + str(l+1)]
        v["db" + str(l+1)] = beta1*v["db" + str(l+1)] + (1 - beta1)*grads['db' + str(l+1)]

        # Compute bias-corrected first moment estimate
        v_corrected["dW" + str(l+1)] = v["dW" + str(l+1)]/(1 - beta1**t)
        v_corrected["db" + str(l+1)] = v["db" + str(l+1)]/(1 - beta1**t)

        # Moving average of the squared gradients
        s["dW" + str(l+1)] = beta2*s["dW" + str(l+1)] + (1 - beta2)*np.square(grads['dW' + str(l+1)])
        s["db" + str(l+1)] = beta2*s["db" + str(l+1)] + (1 - beta2)*np.square(grads['db' + str(l+1)])

        # Compute bias-corrected second raw moment estimate
        s_corrected["dW" + str(l+1)] = s["dW" + str(l+1)]/(1 - beta2**t)
        s_corrected["db" + str(l+1)] = s["db" + str(l+1)]/(1 - beta2**t)

        # Update parameters
        parameters["W" + str(l+1)] = parameters["W" + str(l+1)] - \
        learning_rate*v_corrected["dW" + str(l+1)]/(np.sqrt(s_corrected["dW" + str(l+1)])+epsilon)
        parameters["b" + str(l+1)] = parameters["b" + str(l+1)] - \
        learning_rate*v_corrected["db" + str(l+1)]/(np.sqrt(s_corrected["db" + str(l+1)])+epsilon)

    return parameters, v, s

reason: import numpy,import scipy


text:
"""
Handling csr vectors
"""

import numba
from numpy import array, int32, float64
from pickle import loads, dumps
from scipy.sparse import csr_matrix
from .vectors import normalize_vec

class FastCsrMatrix(csr_matrix):
    """a modified csr_matrix class that does not perform checks"""

    def check_format(self, full_check=True):
        pass

def csr_to_bytes(x):
    """Convert a csr vector (one-row matrix) into a bytes-like object

    :param x: csr matrix
    :return: bytes-like object
    """
    return dumps((tuple([float(_) for _ in x.data]),
                  tuple([int(_) for _ in x.indices])))

def bytes_to_csr(x, ncol):
    """convert a bytes object into a csr one-row matrix

    :param x: bytes object
    :param ncol: integer, number of columns in csr matrix
    :return: csr_matrix
    """
    raw = loads(x)
    return FastCsrMatrix((raw[0], raw[1], (0, len(raw[1]))), shape=(1, ncol))

def bytes_to_arrays(x):
    """convert a bytes object into a pair of arrays

    :param x: bytes object
    :return: arrays with data, indices, and a sum of the data
        Note the first two elements are ready to use with csr_matrix.
        The abs max and abs-runner-up are ad-hoc extras.
    """
    raw = loads(x)
    values = array(raw[0], dtype=float64, copy=False)
    indices = array(raw[1], dtype=int32, copy=False)
    return values, indices

def normalize_csr(v):
    """normalize a csr vector

    Note this performs a normalization in-place,
    i.e. the original vector will change.

    :param v: csr_matrix, assumed a single row vector
    :return: csr_matrix of same size as a, with normalized entries
    """

    v.data = normalize_vec(v.data)
    return v

@numba.jit
def threshold_csr_arrays(data, indices, threshold):
    """helper to threshold_csr"""
    result_data, result_indices = [], []
    for i in range(len(data)):
        d = data[i]
        if abs(d) > threshold:
            result_data.append(d)
            result_indices.append(indices[i])
    return result_data, result_indices

def csr_vector(data, indices, n):
    """create a csr_matrix with only one row, i.e. a vector"""
    result = FastCsrMatrix(([], [], (0, 0)), shape=(1, n))
    result.data = array(data, dtype=float64, copy=False)
    result.indices = array(indices, dtype=int32, copy=False)
    result.indptr = array([0, len(indices)], dtype=int32, copy=False)
    return result

def threshold_csr(v, threshold=0.001):
    """set elements in v below a threshold to zero

    :param v: csr vector
    :param threshold: threshold level
    :return: new csr vector with some elements set to zero
    """

    data, indices = threshold_csr_arrays(v.data, v.indices, threshold)
    return csr_vector(data, indices, v.shape[1])

@numba.jit
def _pos_neg_csr_arrays(data, indices):
    """helper to split a data vector into positive and negative components"""
    pos_data, neg_data = [], []
    pos_indices, neg_indices = [], []
    for i in range(len(data)):
        d = data[i]
        if d > 0:
            pos_data.append(d)
            pos_indices.append(indices[i])
        elif d < 0:
            neg_data.append(d)
            neg_indices.append(indices[i])
    return pos_data, pos_indices, neg_data, neg_indices

def dimcollapse_csr(v, indexes=(), normalize=True):
    """dimensional collapse of a vector

    :param v: csr vector
    :param indexes: allowed dimensional indexes
    :param normalize: logical, set True to rescale values to unit norm
    :return: new csr vector, values outside indexes reset to zero
    """

    data = []
    indices = []
    for d, i in zip(v.data, v.indices):
        if i in indexes:
            data.append(d)
            indices.append(i)
    if normalize:
        data = normalize_vec(array(data))
    result = FastCsrMatrix((data, indices, [0, len(indices)]), shape=v.shape)
    return result

@numba.jit
def get_value_csr(data, indices, index):
    """get one value from a sparse vector"""
    for i in range(len(indices)):
        if indices[i] == index:
            return data[i]
    return 0.0

@numba.jit
def diffuse_sparse(arr, concentrations, origin_index, data, indices):
    """add sparse data to a dense array with a gradient add

    :param arr: dense array of floats
    :param concentrations: dense array of concentrations
    :param origin_index: index in concentrations for the source of diffusion
    :param data: array of floats with diffusion magnitudes
    :param indices: array of integers for destination of diffusion
    :return: array consisting of arr+magnitude*(delta concentration)
    """
    c0 = concentrations[origin_index]
    if c0 > 0:
        for i in range(len(indices)):
            index = indices[i]
            diff = c0-concentrations[index]
            if diff > 0:
                arr[index] += data[i] * diff
                arr[index] = min(c0, arr[index])
    else:
        for i in range(len(indices)):
            index = indices[i]
            diff = c0-concentrations[index]
            if diff < 0:
                arr[index] += data[i] * diff
                arr[index] = max(c0, arr[index])

@numba.jit
def add_sparse_skip(arr, data, indices, skip_index=-1):
    """add sparse data to a dense array

    (Similar to add_sparse, but with an option to skip an index)

    :param arr: dense array of floats
    :param data: array of floats
    :param indices: array of integers
    :param skip_index: integer, skip adding this index
    :return: array consisting of arr+data
    """
    for i in range(len(indices)):
        if indices[i] != skip_index:
            arr[indices[i]] += data[i]
    return arr

@numba.jit
def harmonic_multiply_sparse(factors, data, indices, harmonic_factor):
    """multiply sparse data by a harmonic-mean-scaled factor

    :param factors: dense array of floats
    :param data: dense array of floats (from a sparse vector)
    :param indices: array of integers (from a sparse vector)
    :param harmonic_factor: float, factor for use with harmonic
        multiplication.
    :return: array consisting of factors*data in sparse format.
        When a harmonic factor is present, the multiplication is by
        (factors)/(factors+harmonic_factor)
    """

    for i in range(len(indices)):
        fi = factors[indices[i]]
        data[i] *= fi / (fi + harmonic_factor)
    return data

@numba.jit
def max_multiply_sparse(factors, data, indices, max_factor):
    """multiply sparse data by a factor with upper bound

    :param factors: dense array of floats
    :param data: array of floats (from a sparse vector)
    :param indices: array of integers (from a sparse vector)
    :param max_factor: float, maximal factor forfor use with harmonic
        multiplication.
    :return: array consisting of factors*data in sparse format.
        When max_factor is present, the factors are at most
        equal to the max_factor float
    """

    result = array([0.0]*len(data))
    for i in range(len(indices)):
        result[i] = data[i] * min(max_factor, factors[indices[i]])
    return result

@numba.jit
def csr_data_indices(arr):
    """extract data and indices arrays from a dense vector

    (this may seem like it may be done faster with np.array.nonzero
    and similar tools, but this implementation is faster)

    :param arr: dense array
    :return: arrays with nonzero data elements, and corresponding indices
    """

    data, indices = [], []
    for i in range(len(arr)):
        if arr[i] != 0.0:
            data.append(arr[i])
            indices.append(i)
    return data, indices

max_stars_repo_name: tkonopka/crossmap

max_issues_repo_name: tkonopka/crossmap

max_forks_repo_name: tkonopka/crossmap

reason: from numpy,from scipy,import numba