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

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

from absl.testing import parameterized
import numpy as np
import tensorflow as tf

from tensorflow_model_optimization.python.core.quantization.keras import quantize
from tensorflow_model_optimization.python.core.quantization.keras import quantize_aware_activation
from tensorflow_model_optimization.python.core.quantization.keras import quantize_layer
from tensorflow_model_optimization.python.core.quantization.keras import quantizers
from tensorflow_model_optimization.python.core.quantization.keras.default_8bit import default_8bit_quantize_configs
from tensorflow_model_optimization.python.core.quantization.keras.default_8bit import default_8bit_transforms
from tensorflow_model_optimization.python.core.quantization.keras.graph_transformations import model_transformer
from tensorflow_model_optimization.python.core.quantization.keras.layers import conv_batchnorm_test_utils

ModelTransformer = model_transformer.ModelTransformer

Conv2DModel = conv_batchnorm_test_utils.Conv2DModel
DepthwiseConv2DModel = conv_batchnorm_test_utils.DepthwiseConv2DModel

keras = tf.keras

Conv2DBatchNormActivationQuantize = default_8bit_transforms.Conv2DBatchNormActivationQuantize
Conv2DBatchNormReLUQuantize = default_8bit_transforms.Conv2DBatchNormReLUQuantize

# TODO(alanchiao): reduce redundancy by parameterizing on Depthwise vs Conv.
class DefaultTransformsTest(tf.test.TestCase, parameterized.TestCase):

  def testTransformsConvBNReLUPattern(self):
    model = Conv2DModel.get_nonfolded_batchnorm_model(
        post_bn_activation=keras.layers.ReLU(6.0), model_type='functional')
    folded_model = Conv2DModel.get_folded_batchnorm_model(
        post_bn_activation=keras.layers.ReLU(6.0), is_quantized=True)

    transformed_model, _ = ModelTransformer(
        model,
        [default_8bit_transforms.Conv2DBatchNormReLU6Fold()]).transform()

    inputs = np.random.standard_normal(Conv2DModel.get_batched_input_shape())
    self.assertAllClose(
        transformed_model.predict(inputs), folded_model.predict(inputs))

  def testTransformsConvBNReLUPatternPreservesWeights(self):
    # random_init to prevent non-random initialization in resulting
    # in same weights between transformed and non-transformed models.
    model = Conv2DModel.get_nonfolded_batchnorm_model(
        post_bn_activation=keras.layers.ReLU(6.0),
        model_type='functional',
        random_init=True)

    transformed_model, _ = ModelTransformer(
        model,
        [default_8bit_transforms.Conv2DBatchNormReLU6Fold()]).transform()

    transformed_weights = transformed_model.get_weights()
    # Remove quantization related weights.
    del transformed_weights[3:8]

    self.assertEqual(len(transformed_weights), len(model.get_weights()))
    for i in range(len(transformed_weights)):
      self.assertAllEqual(transformed_weights[i], model.get_weights()[i])

  def testTransformsConvBNPattern(self):
    model = Conv2DModel.get_nonfolded_batchnorm_model(
        model_type='functional')
    folded_model = Conv2DModel.get_folded_batchnorm_model(
        is_quantized=True)

    with quantize.quantize_scope():
      transformed_model, _ = ModelTransformer(
          model, [default_8bit_transforms.Conv2DBatchNormFold()]).transform()

    inputs = np.random.standard_normal(Conv2DModel.get_batched_input_shape())
    self.assertAllClose(
        transformed_model.predict(inputs), folded_model.predict(inputs))

  def testTransformsConvBNPatternPreservesWeights(self):
    # random_init to prevent non-random initialization in resulting
    # in same weights between transformed and non-transformed models.
    model = Conv2DModel.get_nonfolded_batchnorm_model(
        model_type='functional',
        random_init=True)

    transformed_model, _ = ModelTransformer(
        model, [default_8bit_transforms.Conv2DBatchNormFold()]).transform()

    transformed_weights = transformed_model.get_weights()
    # Remove quantization related weights.
    del transformed_weights[3:8]

    self.assertEqual(len(transformed_weights), len(model.get_weights()))
    for i in range(len(transformed_weights)):
      self.assertAllEqual(transformed_weights[i], model.get_weights()[i])

  def testTransformsDepthwiseConvBNReLUPattern(self):
    model = DepthwiseConv2DModel.get_nonfolded_batchnorm_model(
        post_bn_activation=keras.layers.ReLU(6.0), model_type='functional')
    folded_model = DepthwiseConv2DModel.get_folded_batchnorm_model(
        post_bn_activation=keras.layers.ReLU(6.0), is_quantized=True)

    transformed_model, _ = ModelTransformer(
        model, [default_8bit_transforms.DepthwiseConv2DBatchNormReLU6Fold()
               ]).transform()

    inputs = np.random.standard_normal(
        DepthwiseConv2DModel.get_batched_input_shape())
    self.assertAllClose(
        transformed_model.predict(inputs), folded_model.predict(inputs))

  def testTransformsDepthwiseConvBNReLUPatternPreservesWeights(self):
    # random_init to prevent non-random initialization in resulting
    # in same weights between transformed and non-transformed models.
    model = DepthwiseConv2DModel.get_nonfolded_batchnorm_model(
        post_bn_activation=keras.layers.ReLU(6.0),
        model_type='functional',
        random_init=True)

    transformed_model, _ = ModelTransformer(
        model, [default_8bit_transforms.DepthwiseConv2DBatchNormReLU6Fold()
               ]).transform()

    transformed_weights = transformed_model.get_weights()
    # Remove quantization related weights.
    del transformed_weights[3:8]

    self.assertEqual(len(transformed_weights), len(model.get_weights()))
    for i in range(len(transformed_weights)):
      self.assertAllEqual(transformed_weights[i], model.get_weights()[i])

  @staticmethod
  def _get_model(layer_type, activation_type):
    activation = None
    if activation_type == 'relu':
      activation = keras.layers.ReLU(6.0)
    elif activation_type == 'act_relu':
      activation = keras.layers.Activation('relu')

    if layer_type == 'Conv2D':
      return Conv2DModel.get_nonfolded_batchnorm_model(
          model_type='functional', post_bn_activation=activation)
    elif layer_type == 'DepthwiseConv2D':
      return DepthwiseConv2DModel.get_nonfolded_batchnorm_model(
          model_type='functional', post_bn_activation=activation)

  @staticmethod
  def _get_input_shape(layer_type):
    if layer_type == 'Conv2D':
      return Conv2DModel.get_batched_input_shape()
    elif layer_type == 'DepthwiseConv2D':
      return DepthwiseConv2DModel.get_batched_input_shape()

  @parameterized.parameters('Conv2D', 'DepthwiseConv2D')
  def testConv2DBatchNormQuantize(self, layer_type):
    model = self._get_model(layer_type, False)
    input_shape = self._get_input_shape(layer_type)

    transformed_model, updated_metadata = ModelTransformer(
        model,
        [default_8bit_transforms.Conv2DBatchNormQuantize()],
    ).transform()

    conv_layer = transformed_model.layers[1]
    bn_layer = transformed_model.layers[2]

    self.assertIsInstance(
        conv_layer.activation, quantize_aware_activation.NoOpActivation)
    self.assertIsInstance(
        updated_metadata.get(bn_layer.name).get('quantize_config'),
        default_8bit_quantize_configs.Default8BitOutputQuantizeConfig)

    inputs = np.random.standard_normal(input_shape)
    self.assertAllClose(
        transformed_model.predict(inputs), model.predict(inputs))

  @parameterized.parameters(
      ('Conv2D', 'relu', Conv2DBatchNormReLUQuantize),
      ('Conv2D', 'act_relu', Conv2DBatchNormActivationQuantize),
      ('DepthwiseConv2D', 'relu', Conv2DBatchNormReLUQuantize),
      ('DepthwiseConv2D', 'act_relu',
       Conv2DBatchNormActivationQuantize),
  )
  def testConv2DBatchNormReLUQuantize(
      self, layer_type, activation_type, transform_type):
    model = self._get_model(layer_type, activation_type)
    input_shape = self._get_input_shape(layer_type)

    transformed_model, updated_metadata = ModelTransformer(
        model,
        [transform_type()],
    ).transform()

    conv_layer = transformed_model.layers[1]
    bn_layer = transformed_model.layers[2]

    self.assertIsInstance(
        conv_layer.activation, quantize_aware_activation.NoOpActivation)
    self.assertIsInstance(
        updated_metadata.get(bn_layer.name).get('quantize_config'),
        default_8bit_quantize_configs.NoOpQuantizeConfig)

    inputs = np.random.standard_normal(input_shape)
    self.assertAllClose(
        transformed_model.predict(inputs), model.predict(inputs))

  @parameterized.named_parameters(
      ('padding_valid', {'padding': 'valid'}),
      ('padding_same', {'padding': 'same'}),
      ('padding_same_dilation_2', {'padding': 'same', 'dilation_rate': 2}),
      ('strides', {'strides': 2}),
      ('dilation_rate', {'dilation_rate': 2}),
      ('depth_multiplier', {'depth_multiplier': 2}),
      ('regularizer', {
          'depthwise_regularizer': 'l2',
          'pointwise_regularizer': 'l2',
          'bias_regularizer': 'l2',
          'activity_regularizer': 'l2'}),
      ('constraint', {
          'depthwise_constraint': tf.keras.constraints.max_norm(2.),
          'pointwise_constraint': tf.keras.constraints.min_max_norm(0., 2.),
          'bias_constraint': tf.keras.constraints.unit_norm()}),
      ('activation_relu', {'activation': 'relu'}),
      # TODO(pulkitb): Temporarily disabling due to numerical errors resulting
      # from caching of activation logits in TF code.
      # ('activation_softmax', {'activation': 'softmax'}),
  )
  def testSeparableConv1DQuantize_(self, kwargs):
    kwargs['filters'] = 2
    kwargs['kernel_size'] = 3
    num_samples = 2
    stack_size = 3
    num_row = 7

    sepconv_model = tf.keras.Sequential([
        tf.keras.Input(
            shape=(num_row, stack_size), batch_size=num_samples),
        tf.keras.layers.SeparableConv1D(**kwargs)])

    transformed_model, updated_metadata = ModelTransformer(
        sepconv_model,
        [default_8bit_transforms.SeparableConv1DQuantize()],
    ).transform()

    self.assertContainsSubset(
        {'sepconv1d_expand_1', 'separable_conv1d_QAT_SepConv2D',
         'sepconv1d_squeeze_1'},
        updated_metadata.keys())
    self.assertEqual(sepconv_model.output_shape, transformed_model.output_shape)

    x = np.random.rand(*sepconv_model.input_shape)
    y = np.random.rand(*sepconv_model.output_shape)

    # Ensure model is equivalent, and forward pass results are the same.
    self.assertAllClose(sepconv_model.predict(x), transformed_model.predict(x))

    # Ensure model is equivalent, and training results are the same.
    sepconv_model.compile(loss='categorical_crossentropy', optimizer='sgd')
    sepconv_model.fit(x, y, epochs=100)
    transformed_model.compile(loss='categorical_crossentropy', optimizer='sgd')
    transformed_model.fit(x, y, epochs=100)

    # Over a long training cycle with constraints and regularizers, the model
    # can build very minute differences. Hence reducing tol to 1e-5.
    self.assertAllClose(sepconv_model.predict(x), transformed_model.predict(x),
                        atol=1e-5, rtol=1e-5)

  @parameterized.named_parameters(
      ('padding_valid', {'padding': 'valid'}),
      ('padding_same', {'padding': 'same'}),
      ('padding_same_dilation_2', {'padding': 'same', 'dilation_rate': 2}),
      ('strides', {'strides': 2}),
      ('dilation_rate', {'dilation_rate': 2}),
      ('depth_multiplier', {'depth_multiplier': 2}),
      ('regularizer', {
          'depthwise_regularizer': 'l2',
          'pointwise_regularizer': 'l2',
          'bias_regularizer': 'l2',
          'activity_regularizer': 'l2'}),
      ('constraint', {
          'depthwise_constraint': tf.keras.constraints.max_norm(2.),
          'pointwise_constraint': tf.keras.constraints.min_max_norm(0., 2.),
          'bias_constraint': tf.keras.constraints.unit_norm()})
  )
  def testSeparableConvQuantize_(self, kwargs):
    kwargs['filters'] = 2
    kwargs['kernel_size'] = 3
    num_samples = 2
    stack_size = 3
    num_row = 7
    num_col = 6

    sepconv_model = tf.keras.Sequential([
        tf.keras.Input(
            shape=(num_row, num_col, stack_size), batch_size=num_samples),
        tf.keras.layers.SeparableConv2D(**kwargs)])

    transformed_model, updated_metadata = ModelTransformer(
        sepconv_model,
        [default_8bit_transforms.SeparableConvQuantize()],
    ).transform()

    self.assertContainsSubset(
        updated_metadata.keys(), {'depthwise_conv2d', 'conv2d'})
    # Transformed model should have the same output shape
    self.assertEqual(sepconv_model.output_shape, transformed_model.output_shape)

    x = np.random.rand(*sepconv_model.input_shape)
    y = np.random.rand(*sepconv_model.output_shape)

    # Ensure model is equivalent, and forward pass results are the same.
    self.assertAllClose(sepconv_model.predict(x), transformed_model.predict(x))

    # Ensure model is equivalent, and training results are the same.
    sepconv_model.compile(loss='categorical_crossentropy', optimizer='sgd')
    sepconv_model.fit(x, y, epochs=100)
    transformed_model.compile(loss='categorical_crossentropy', optimizer='sgd')
    transformed_model.fit(x, y, epochs=100)

    # Over a long training cycle with constraints and regularizers, the model
    # can build very minute differences. Hence reducing tol to 1e-5.
    self.assertAllClose(sepconv_model.predict(x), transformed_model.predict(x),
                        atol=1e-5, rtol=1e-5)

  # TODO(pulkitb): Add individual tests for the following transforms.
  # Conv2DReshapeBatchNormQuantize, Conv2DReshapeBatchNormReLUQuantize
  # Conv2DReshapeBatchNormActivationQuantize

  @parameterized.parameters(
      ('relu', default_8bit_transforms.LayerReLUQuantize),
      ('act_relu', default_8bit_transforms.LayerReluActivationQuantize),
  )
  def testAddReLUQuantize(self, activation_type, transform_type):
    add = keras.layers.Add()
    if activation_type == 'relu':
      activation = keras.layers.ReLU(6.0)
    elif activation_type == 'act_relu':
      activation = keras.layers.Activation('relu')

    inp1 = keras.layers.Input((3,))
    inp2 = keras.layers.Input((3,))
    x = activation(add([inp1, inp2]))
    model = keras.Model([inp1, inp2], x)

    transformed_model, updated_metadata = ModelTransformer(
        model,
        [transform_type()],
    ).transform()

    add_layer = transformed_model.layers[2]

    self.assertIsInstance(
        updated_metadata.get(add_layer.name).get('quantize_config'),
        default_8bit_quantize_configs.NoOpQuantizeConfig)

  @parameterized.parameters(
      ('relu', default_8bit_transforms.LayerReLUQuantize),
      ('act_relu', default_8bit_transforms.LayerReluActivationQuantize))
  def testLayerReLUQuantize(self, activation_type, transform_type):
    # TODO(tfmot): Add tests for DepthConv and Dense
    input_shape = (1, 3, 3, 3)
    conv_layer = tf.keras.layers.Conv2D(5, 2, input_shape=input_shape)
    if activation_type == 'relu':
      act_layer = keras.layers.ReLU(6.0)
    elif activation_type == 'act_relu':
      act_layer = keras.layers.Activation('relu')

    model = tf.keras.Sequential([conv_layer, act_layer])

    transformed_model, updated_metadata = ModelTransformer(
        model,
        [transform_type()],
    ).transform()

    self.assertIsInstance(
        updated_metadata.get(model.layers[0].name).get('quantize_config'),
        default_8bit_quantize_configs.NoOpQuantizeConfig)

    inputs = np.random.standard_normal(input_shape)
    self.assertAllClose(
        transformed_model.predict(inputs), model.predict(inputs))

  def testAddsQuantizeLayerAfterInputLayer(self):
    inp1 = keras.layers.Input((3,))
    inp2 = keras.layers.Input((3,))
    x = keras.layers.Concatenate()([inp1, inp2])
    model = keras.Model([inp1, inp2], x)

    transformed_model, _ = ModelTransformer(
        model,
        [default_8bit_transforms.InputLayerQuantize()]).transform()

    for input_layer in transformed_model._input_layers:
      layer_after_input = input_layer._outbound_nodes[0].outbound_layer
      self.assertIsInstance(
          layer_after_input,
          quantize_layer.QuantizeLayer)
      self.assertIsInstance(
          layer_after_input.quantizer, quantizers.AllValuesQuantizer)

  def testConcatTransform(self):
    r"""Tests the Concat Transform.

               Input
              /     \
         Dense       Dense
             \      /
              Concat

      One Dense layer has a pre-specified QuantizeConfig, whereas the other does
      not. The Transform should ensure both the output FakeQuants are disabled,
      and only a FakeQuant after Concat is present.
    """
    dense_1 = keras.layers.Dense(3)
    dense_2 = keras.layers.Dense(3)
    concat = keras.layers.Concatenate()

    inp = keras.layers.Input((2,))
    x1 = dense_1(inp)
    x2 = dense_2(inp)
    x = concat([x1, x2])
    model = keras.Model(inp, x)

    layer_metadata = {
        # dense_1 has an existing quantize_config.
        dense_1.name: {
            'quantize_config':
                default_8bit_quantize_configs.Default8BitOutputQuantizeConfig()
        }
    }
    _, updated_metadata = ModelTransformer(
        model, [default_8bit_transforms.ConcatTransform()],
        layer_metadata=layer_metadata).transform()

    concat_quantize_config = updated_metadata.get(
        concat.name).get('quantize_config')
    # Concat should quantize the output.
    self.assertIsInstance(
        concat_quantize_config,
        default_8bit_quantize_configs.Default8BitOutputQuantizeConfig)
    self.assertNotEmpty(concat_quantize_config.get_output_quantizers(None))

    dense_1_quantize_config = updated_metadata.get(
        dense_1.name).get('quantize_config')
    # The existing quantize_config should do nothing for outputs.
    self.assertIsInstance(
        dense_1_quantize_config,
        default_8bit_quantize_configs.Default8BitOutputQuantizeConfig)
    self.assertEmpty(dense_1_quantize_config.get_output_quantizers(None))

    dense_2_quantize_config = updated_metadata.get(
        dense_2.name).get('quantize_config')
    # The quantize_config from registry should do nothing at output.
    self.assertEqual('Default8BitQuantizeConfig',
                     dense_2_quantize_config.__class__.__name__)
    self.assertEmpty(dense_2_quantize_config.get_output_quantizers(None))

  def testConcatMultipleLevels(self):
    r"""Tests case when concats applied to concats.

            Input --------------.
           /      \      |      |
         Dense   Dense   |      |
            \    /       |      |
             Concat    Dense   Dense
                 \     /        |
                  Concat        |
                        \      /
                         Concat

    The last Concat layer should be quantized but the rest
    of the outputs should just feed into it.
    """
    inp = keras.layers.Input((3,))
    x1 = keras.layers.Dense(3)(inp)
    x2 = keras.layers.Dense(3)(inp)
    x3 = keras.layers.Dense(3)(inp)
    x4 = keras.layers.Dense(3)(inp)
    c1 = keras.layers.Concatenate()([x1, x2])
    c2 = keras.layers.Concatenate()([c1, x3])
    c3 = keras.layers.Concatenate()([c2, x4])
    model = keras.Model(inp, c3)
    model.summary()

    _, layer_metadata = ModelTransformer(
        model,
        [default_8bit_transforms.ConcatTransform()]).transform()

    for layer in model.layers[1:-1]:
      quantize_config = layer_metadata[layer.name].get('quantize_config')
      self.assertEmpty(quantize_config.get_output_quantizers(None))

    c3_layer = model.layers[-1]
    quantize_config = layer_metadata[c3_layer.name].get('quantize_config')
    self.assertIsInstance(
        quantize_config,
        default_8bit_quantize_configs.Default8BitOutputQuantizeConfig)
    self.assertNotEmpty(quantize_config.get_output_quantizers(None))

if __name__ == '__main__':
  tf.test.main()

max_stars_repo_name: markostam/model-optimization

max_issues_repo_name: markostam/model-optimization

max_forks_repo_name: markostam/model-optimization


text:
# Author: Eric Larson <larson.eric.d@gmail.com>
#
# License: BSD-3-Clause

import os.path as op

import numpy as np
from numpy.testing import (assert_allclose, assert_array_less,
                           assert_array_equal)
from scipy.interpolate import interp1d
from scipy.spatial.distance import cdist
import pytest

from mne import pick_types, pick_info
from mne.forward._compute_forward import _MAG_FACTOR
from mne.io import (read_raw_fif, read_raw_artemis123, read_raw_ctf, read_info,
                    RawArray, read_raw_kit)
from mne.io.constants import FIFF
from mne.chpi import (compute_chpi_amplitudes, compute_chpi_locs,
                      compute_chpi_snr, compute_head_pos, _setup_ext_proj,
                      _chpi_locs_to_times_dig, _compute_good_distances,
                      extract_chpi_locs_ctf, head_pos_to_trans_rot_t,
                      read_head_pos, write_head_pos, filter_chpi,
                      get_chpi_info, _get_hpi_initial_fit,
                      extract_chpi_locs_kit)
from mne.datasets import testing
from mne.simulation import add_chpi
from mne.transforms import rot_to_quat, _angle_between_quats
from mne.utils import catch_logging, assert_meg_snr, verbose
from mne.viz import plot_head_positions

base_dir = op.join(op.dirname(__file__), '..', 'io', 'tests', 'data')
ctf_fname = op.join(base_dir, 'test_ctf_raw.fif')
hp_fif_fname = op.join(base_dir, 'test_chpi_raw_sss.fif')
hp_fname = op.join(base_dir, 'test_chpi_raw_hp.txt')
raw_fname = op.join(base_dir, 'test_raw.fif')

data_path = testing.data_path(download=False)
sample_fname = op.join(
    data_path, 'MEG', 'sample', 'sample_audvis_trunc_raw.fif')
chpi_fif_fname = op.join(data_path, 'SSS', 'test_move_anon_raw.fif')
pos_fname = op.join(data_path, 'SSS', 'test_move_anon_raw.pos')
sss_fif_fname = op.join(data_path, 'SSS', 'test_move_anon_raw_sss.fif')
sss_hpisubt_fname = op.join(data_path, 'SSS', 'test_move_anon_hpisubt_raw.fif')
chpi5_fif_fname = op.join(data_path, 'SSS', 'chpi5_raw.fif')
chpi5_pos_fname = op.join(data_path, 'SSS', 'chpi5_raw_mc.pos')
ctf_chpi_fname = op.join(data_path, 'CTF', 'testdata_ctf_mc.ds')
ctf_chpi_pos_fname = op.join(data_path, 'CTF', 'testdata_ctf_mc.pos')

art_fname = op.join(data_path, 'ARTEMIS123', 'Artemis_Data_2017-04-04' +
                    '-15h-44m-22s_Motion_Translation-z.bin')
art_mc_fname = op.join(data_path, 'ARTEMIS123', 'Artemis_Data_2017-04-04' +
                       '-15h-44m-22s_Motion_Translation-z_mc.pos')

con_fname = op.join(data_path, 'KIT', 'MQKIT_125_2sec.con')
mrk_fname = op.join(data_path, 'KIT', 'MQKIT_125.mrk')
elp_fname = op.join(data_path, 'KIT', 'MQKIT_125.elp')
hsp_fname = op.join(data_path, 'KIT', 'MQKIT_125.hsp')
berlin_fname = op.join(data_path, 'KIT', 'data_berlin.con')

@testing.requires_testing_data
def test_chpi_adjust():
    """Test cHPI logging and adjustment."""
    raw = read_raw_fif(chpi_fif_fname, allow_maxshield='yes')
    with catch_logging() as log:
        _get_hpi_initial_fit(raw.info, adjust=True, verbose='debug')
        get_chpi_info(raw.info, on_missing='raise', verbose='debug')
    # Ran MaxFilter (with -list, -v, -movecomp, etc.), and got:
    msg = ['HPIFIT: 5 coils digitized in order 5 1 4 3 2',
           'HPIFIT: 3 coils accepted: 1 2 4',
           'Hpi coil moments (3 5):',
           '2.08542e-15 -1.52486e-15 -1.53484e-15',
           '2.14516e-15 2.09608e-15 7.30303e-16',
           '-3.2318e-16 -4.25666e-16 2.69997e-15',
           '5.21717e-16 1.28406e-15 1.95335e-15',
           '1.21199e-15 -1.25801e-19 1.18321e-15',
           'HPIFIT errors:  0.3, 0.3, 5.3, 0.4, 3.2 mm.',
           'HPI consistency of isotrak and hpifit is OK.',
           'HP fitting limits: err = 5.0 mm, gval = 0.980.',
           'Using 5 HPI coils: 83 143 203 263 323 Hz',  # actually came earlier
           ]

    log = log.getvalue().splitlines()
    assert set(log) == set(msg), '\n' + '\n'.join(set(msg) - set(log))

    # Then took the raw file, did this:
    raw.info['dig'][5]['r'][2] += 1.
    # And checked the result in MaxFilter, which changed the logging as:
    msg = msg[:8] + [
        'HPIFIT errors:  0.3, 0.3, 5.3, 999.7, 3.2 mm.',
        'Note: HPI coil 3 isotrak is adjusted by 5.3 mm!',
        'Note: HPI coil 5 isotrak is adjusted by 3.2 mm!'] + msg[-2:]
    with catch_logging() as log:
        _get_hpi_initial_fit(raw.info, adjust=True, verbose='debug')
        get_chpi_info(raw.info, on_missing='raise', verbose='debug')
    log = log.getvalue().splitlines()
    assert set(log) == set(msg), '\n' + '\n'.join(set(msg) - set(log))

@testing.requires_testing_data
def test_read_write_head_pos(tmpdir):
    """Test reading and writing head position quaternion parameters."""
    temp_name = op.join(str(tmpdir), 'temp.pos')
    # This isn't a 100% valid quat matrix but it should be okay for tests
    head_pos_rand = np.random.RandomState(0).randn(20, 10)
    # This one is valid
    head_pos_read = read_head_pos(pos_fname)
    for head_pos_orig in (head_pos_rand, head_pos_read):
        write_head_pos(temp_name, head_pos_orig)
        head_pos = read_head_pos(temp_name)
        assert_allclose(head_pos_orig, head_pos, atol=1e-3)
    # Degenerate cases
    pytest.raises(TypeError, write_head_pos, 0, head_pos_read)  # not filename
    pytest.raises(ValueError, write_head_pos, temp_name, 'foo')  # not array
    pytest.raises(ValueError, write_head_pos, temp_name, head_pos_read[:, :9])
    pytest.raises(TypeError, read_head_pos, 0)
    pytest.raises(IOError, read_head_pos, temp_name + 'foo')

@testing.requires_testing_data
def test_hpi_info(tmpdir):
    """Test getting HPI info."""
    temp_name = op.join(str(tmpdir), 'temp_raw.fif')
    for fname in (chpi_fif_fname, sss_fif_fname):
        raw = read_raw_fif(fname, allow_maxshield='yes').crop(0, 0.1)
        assert len(raw.info['hpi_subsystem']) > 0
        raw.save(temp_name, overwrite=True)
        info = read_info(temp_name)
        assert len(info['hpi_subsystem']) == len(raw.info['hpi_subsystem'])

    # test get_chpi_info()
    info = read_info(chpi_fif_fname)
    hpi_freqs, stim_ch_idx, hpi_on_codes = get_chpi_info(info)

    assert_allclose(hpi_freqs, np.array([83., 143., 203., 263., 323.]))
    assert stim_ch_idx == 378
    assert_allclose(hpi_on_codes, np.array([256, 512, 1024, 2048, 4096]))

    # test get_chpi_info() if no proper cHPI info is available
    info['hpi_subsystem'] = None
    info['hpi_meas'] = []
    info['hpi_results'] = []

    with pytest.raises(ValueError, match='No appropriate cHPI information'):
        get_chpi_info(info)

    with pytest.warns(RuntimeWarning, match='No appropriate cHPI information'):
        get_chpi_info(info, on_missing='warn')

    hpi_freqs, stim_ch_idx, hpi_on_codes = get_chpi_info(info,
                                                         on_missing='ignore')
    assert_array_equal([], hpi_freqs)
    assert stim_ch_idx is None
    assert_array_equal([], hpi_on_codes)

def _assert_quats(actual, desired, dist_tol=0.003, angle_tol=5., err_rtol=0.5,
                  gof_rtol=0.001, vel_atol=2e-3):  # 2 mm/s
    """Compare estimated cHPI positions."""
    __tracebackhide__ = True
    trans_est, rot_est, t_est = head_pos_to_trans_rot_t(actual)
    trans, rot, t = head_pos_to_trans_rot_t(desired)
    quats_est = rot_to_quat(rot_est)
    gofs, errs, vels = desired[:, 7:].T
    gofs_est, errs_est, vels_est = actual[:, 7:].T
    del actual, desired

    # maxfilter produces some times that are implausibly large (weird)
    if not np.isclose(t[0], t_est[0], atol=1e-1):  # within 100 ms
        raise AssertionError('Start times not within 100 ms: %0.3f != %0.3f'
                             % (t[0], t_est[0]))
    use_mask = (t >= t_est[0]) & (t <= t_est[-1])
    t = t[use_mask]
    trans = trans[use_mask]
    quats = rot_to_quat(rot)
    quats = quats[use_mask]
    gofs, errs, vels = gofs[use_mask], errs[use_mask], vels[use_mask]

    # double-check our angle function
    for q in (quats, quats_est):
        angles = _angle_between_quats(q, q)
        assert_allclose(angles, 0., atol=1e-5)

    # limit translation difference between MF and our estimation
    trans_est_interp = interp1d(t_est, trans_est, axis=0)(t)
    distances = np.sqrt(np.sum((trans - trans_est_interp) ** 2, axis=1))
    assert np.isfinite(distances).all()
    arg_worst = np.argmax(distances)
    assert distances[arg_worst] <= dist_tol, (
        '@ %0.3f seconds: %0.3f > %0.3f mm'
        % (t[arg_worst], 1000 * distances[arg_worst], 1000 * dist_tol))

    # limit rotation difference between MF and our estimation
    # (note that the interpolation will make this slightly worse)
    quats_est_interp = interp1d(t_est, quats_est, axis=0)(t)
    angles = 180 * _angle_between_quats(quats_est_interp, quats) / np.pi
    arg_worst = np.argmax(angles)
    assert angles[arg_worst] <= angle_tol, (
        '@ %0.3f seconds: %0.3f > %0.3f deg'
        % (t[arg_worst], angles[arg_worst], angle_tol))

    # error calculation difference
    errs_est_interp = interp1d(t_est, errs_est)(t)
    assert_allclose(errs_est_interp, errs, rtol=err_rtol, atol=1e-3,
                    err_msg='err')  # 1 mm

    # gof calculation difference
    gof_est_interp = interp1d(t_est, gofs_est)(t)
    assert_allclose(gof_est_interp, gofs, rtol=gof_rtol, atol=1e-7,
                    err_msg='gof')

    # velocity calculation difference
    vel_est_interp = interp1d(t_est, vels_est)(t)
    assert_allclose(vel_est_interp, vels, atol=vel_atol,
                    err_msg='velocity')

def _decimate_chpi(raw, decim=4):
    """Decimate raw data (with aliasing) in cHPI-fitting compatible way."""
    raw_dec = RawArray(
        raw._data[:, ::decim], raw.info, first_samp=raw.first_samp // decim)
    raw_dec.info['sfreq'] /= decim
    for coil in raw_dec.info['hpi_meas'][0]['hpi_coils']:
        if coil['coil_freq'] > raw_dec.info['sfreq']:
            coil['coil_freq'] = np.mod(coil['coil_freq'],
                                       raw_dec.info['sfreq'])
            if coil['coil_freq'] > raw_dec.info['sfreq'] / 2.:
                coil['coil_freq'] = raw_dec.info['sfreq'] - coil['coil_freq']
    return raw_dec

# A shortcut method for testing that does both steps
@verbose
def _calculate_chpi_positions(raw, t_step_min=0.01, t_step_max=1.,
                              t_window='auto', too_close='raise',
                              dist_limit=0.005, gof_limit=0.98,
                              ext_order=1, verbose=None):
    chpi_amplitudes = compute_chpi_amplitudes(
        raw, t_step_min=t_step_min, t_window=t_window,
        ext_order=ext_order, verbose=verbose)
    chpi_locs = compute_chpi_locs(
        raw.info, chpi_amplitudes, t_step_max=t_step_max,
        too_close=too_close, verbose=verbose)
    head_pos = compute_head_pos(
        raw.info, chpi_locs, dist_limit=dist_limit, gof_limit=gof_limit,
        verbose=verbose)
    return head_pos

@pytest.mark.slowtest
@testing.requires_testing_data
def test_calculate_chpi_positions_vv():
    """Test calculation of cHPI positions."""
    # Check to make sure our fits match MF decently
    mf_quats = read_head_pos(pos_fname)
    raw = read_raw_fif(chpi_fif_fname, allow_maxshield='yes')
    raw.crop(0, 5).load_data()
    # check "auto" t_window estimation at full sampling rate
    with catch_logging() as log:
        compute_chpi_amplitudes(raw, t_step_min=0.1, t_window='auto',
                                tmin=0, tmax=2, verbose=True)
    assert '83.3 ms' in log.getvalue()
    # This is a little hack (aliasing while decimating) to make it much faster
    # for testing purposes only. We can relax this later if we find it breaks
    # something.
    raw_dec = _decimate_chpi(raw, 15)
    with catch_logging() as log:
        with pytest.warns(RuntimeWarning, match='cannot determine'):
            py_quats = _calculate_chpi_positions(raw_dec, t_window=0.2,
                                                 verbose='debug')
    log = log.getvalue()
    assert '\nHPIFIT' in log
    assert 'Computing 4385 HPI location guesses' in log
    _assert_quats(py_quats, mf_quats, dist_tol=0.001, angle_tol=0.7)
    # degenerate conditions
    raw_no_chpi = read_raw_fif(sample_fname)
    with pytest.raises(ValueError, match='No appropriate cHPI information'):
        _calculate_chpi_positions(raw_no_chpi)
    raw_bad = raw.copy()
    del raw_bad.info['hpi_meas'][0]['hpi_coils'][0]['coil_freq']
    with pytest.raises(ValueError, match='No appropriate cHPI information'):
        _calculate_chpi_positions(raw_bad)
    raw_bad = raw.copy()
    for d in raw_bad.info['dig']:
        if d['kind'] == FIFF.FIFFV_POINT_HPI:
            d['coord_frame'] = FIFF.FIFFV_COORD_UNKNOWN
            break
    with pytest.raises(RuntimeError, match='coordinate frame incorrect'):
        _calculate_chpi_positions(raw_bad)
    for d in raw_bad.info['dig']:
        if d['kind'] == FIFF.FIFFV_POINT_HPI:
            d['coord_frame'] = FIFF.FIFFV_COORD_HEAD
            d['r'] = np.ones(3)
    raw_bad.crop(0, 1.)
    picks = np.concatenate([np.arange(306, len(raw_bad.ch_names)),
                            pick_types(raw_bad.info, meg=True)[::16]])
    raw_bad.pick_channels([raw_bad.ch_names[pick] for pick in picks])
    with pytest.warns(RuntimeWarning, match='Discrepancy'):
        with catch_logging() as log_file:
            _calculate_chpi_positions(raw_bad, t_step_min=1., verbose=True)
    # ignore HPI info header and [done] footer
    assert '0/5 good HPI fits' in log_file.getvalue()

    # half the rate cuts off cHPI coils
    raw.info['lowpass'] /= 2.
    with pytest.raises(RuntimeError, match='above the'):
        _calculate_chpi_positions(raw)

@pytest.mark.slowtest
def test_calculate_chpi_snr():
    """Test cHPI SNR calculation."""
    raw = read_raw_fif(chpi_fif_fname, allow_maxshield='yes')
    result = compute_chpi_snr(raw)
    # make sure all the entries are there
    keys = {f'{ch_type}_{key}' for ch_type in ('mag', 'grad') for key in
            ('snr', 'power', 'resid')}
    assert set(result) == keys.union({'times', 'freqs'})
    # make sure the values are plausible, given the sample data file
    assert result['mag_snr'].min() > 1
    assert result['mag_snr'].max() < 40
    assert result['grad_snr'].min() > 1
    assert result['grad_snr'].max() < 40

@testing.requires_testing_data
@pytest.mark.slowtest
def test_calculate_chpi_positions_artemis():
    """Test on 5k artemis data."""
    raw = read_raw_artemis123(art_fname, preload=True)
    mf_quats = read_head_pos(art_mc_fname)
    mf_quats[:, 8:] /= 100  # old code errantly had this factor
    py_quats = _calculate_chpi_positions(raw, t_step_min=2., verbose='debug')
    _assert_quats(
        py_quats, mf_quats,
        dist_tol=0.001, angle_tol=1., err_rtol=0.7, vel_atol=1e-2)

@testing.requires_testing_data
def test_initial_fit_redo():
    """Test that initial fits can be redone based on moments."""
    raw = read_raw_fif(chpi_fif_fname, allow_maxshield='yes')
    slopes = np.array(
        [[c['slopes'] for c in raw.info['hpi_meas'][0]['hpi_coils']]])
    amps = np.linalg.norm(slopes, axis=-1)
    amps /= slopes.shape[-1]
    assert_array_less(amps, 5e-11)
    assert_array_less(1e-12, amps)
    proj, _, _ = _setup_ext_proj(raw.info, ext_order=1)
    chpi_amplitudes = dict(times=np.zeros(1), slopes=slopes, proj=proj)
    chpi_locs = compute_chpi_locs(raw.info, chpi_amplitudes)

    # check GOF
    coil_gof = raw.info['hpi_results'][0]['goodness']
    assert_allclose(chpi_locs['gofs'][0], coil_gof, atol=0.3)  # XXX not good

    # check moment
    # XXX our forward and theirs differ by an extra mult by _MAG_FACTOR
    coil_moment = raw.info['hpi_results'][0]['moments'] / _MAG_FACTOR
    py_moment = chpi_locs['moments'][0]
    coil_amp = np.linalg.norm(coil_moment, axis=-1, keepdims=True)
    py_amp = np.linalg.norm(py_moment, axis=-1, keepdims=True)
    assert_allclose(coil_amp, py_amp, rtol=0.2)
    coil_ori = coil_moment / coil_amp
    py_ori = py_moment / py_amp
    angles = np.rad2deg(np.arccos(np.abs(np.sum(coil_ori * py_ori, axis=1))))
    assert_array_less(angles, 20)

    # check resulting dev_head_t
    head_pos = compute_head_pos(raw.info, chpi_locs)
    assert head_pos.shape == (1, 10)
    nm_pos = raw.info['dev_head_t']['trans']
    dist = 1000 * np.linalg.norm(nm_pos[:3, 3] - head_pos[0, 4:7])
    assert 0.1 < dist < 2
    angle = np.rad2deg(_angle_between_quats(
        rot_to_quat(nm_pos[:3, :3]), head_pos[0, 1:4]))
    assert 0.1 < angle < 2
    gof = head_pos[0, 7]
    assert_allclose(gof, 0.9999, atol=1e-4)

@testing.requires_testing_data
def test_calculate_head_pos_chpi_on_chpi5_in_one_second_steps():
    """Comparing estimated cHPI positions with MF results (one second)."""
    # Check to make sure our fits match MF decently
    mf_quats = read_head_pos(chpi5_pos_fname)
    raw = read_raw_fif(chpi5_fif_fname, allow_maxshield='yes')
    # the last two seconds contain a maxfilter problem!
    # fiff file timing: 26. to 43. seconds
    # maxfilter estimates a wrong head position for interval 16: 41.-42. sec
    raw = _decimate_chpi(raw.crop(0., 10.).load_data(), decim=8)
    # needs no interpolation, because maxfilter pos files comes with 1 s steps
    py_quats = _calculate_chpi_positions(
        raw, t_step_min=1.0, t_step_max=1.0, t_window=1.0, verbose='debug')
    _assert_quats(py_quats, mf_quats, dist_tol=0.002, angle_tol=1.2,
                  vel_atol=3e-3)  # 3 mm/s

@pytest.mark.slowtest
@testing.requires_testing_data
def test_calculate_head_pos_chpi_on_chpi5_in_shorter_steps():
    """Comparing estimated cHPI positions with MF results (smaller steps)."""
    # Check to make sure our fits match MF decently
    mf_quats = read_head_pos(chpi5_pos_fname)
    raw = read_raw_fif(chpi5_fif_fname, allow_maxshield='yes')
    raw = _decimate_chpi(raw.crop(0., 5.).load_data(), decim=8)
    with pytest.warns(RuntimeWarning, match='cannot determine'):
        py_quats = _calculate_chpi_positions(
            raw, t_step_min=0.1, t_step_max=0.1, t_window=0.1, verbose='debug')
    # needs interpolation, tolerance must be increased
    _assert_quats(py_quats, mf_quats, dist_tol=0.002, angle_tol=1.2,
                  vel_atol=0.02)  # 2 cm/s is not great but probably fine

def test_simulate_calculate_head_pos_chpi():
    """Test calculation of cHPI positions with simulated data."""
    # Read info dict from raw FIF file
    info = read_info(raw_fname)
    # Tune the info structure
    chpi_channel = u'STI201'
    ncoil = len(info['hpi_results'][0]['order'])
    coil_freq = 10 + np.arange(ncoil) * 5
    hpi_subsystem = {'event_channel': chpi_channel,
                     'hpi_coils': [{'event_bits': np.array([256, 0, 256, 256],
                                                           dtype=np.int32)},
                                   {'event_bits': np.array([512, 0, 512, 512],
                                                           dtype=np.int32)},
                                   {'event_bits':
                                       np.array([1024, 0, 1024, 1024],
                                                dtype=np.int32)},
                                   {'event_bits':
                                       np.array([2048, 0, 2048, 2048],
                                                dtype=np.int32)}],
                     'ncoil': ncoil}

    info['hpi_subsystem'] = hpi_subsystem
    for fi, freq in enumerate(coil_freq):
        info['hpi_meas'][0]['hpi_coils'][fi]['coil_freq'] = freq
    picks = pick_types(info, meg=True, stim=True, eeg=False, exclude=[])
    info['sfreq'] = 100.  # this will speed it up a lot
    info = pick_info(info, picks)
    info['chs'][info['ch_names'].index('STI 001')]['ch_name'] = 'STI201'
    info._update_redundant()
    info['projs'] = []

    info_trans = info['dev_head_t']['trans'].copy()

    dev_head_pos_ini = np.concatenate([rot_to_quat(info_trans[:3, :3]),
                                       info_trans[:3, 3]])
    ez = np.array([0, 0, 1])  # Unit vector in z-direction of head coordinates

    # Define some constants
    duration = 10  # Time / s

    # Quotient of head position sampling frequency
    # and raw sampling frequency
    head_pos_sfreq_quotient = 0.01

    # Round number of head positions to the next integer
    S = int(duration * info['sfreq'] * head_pos_sfreq_quotient)
    assert S == 10
    dz = 0.001  # Shift in z-direction is 0.1mm for each step

    dev_head_pos = np.zeros((S, 10))
    dev_head_pos[:, 0] = np.arange(S) * info['sfreq'] * head_pos_sfreq_quotient
    dev_head_pos[:, 1:4] = dev_head_pos_ini[:3]
    dev_head_pos[:, 4:7] = dev_head_pos_ini[3:] + \
        np.outer(np.arange(S) * dz, ez)
    dev_head_pos[:, 7] = 1.0

    # m/s
    dev_head_pos[:, 9] = dz / (info['sfreq'] * head_pos_sfreq_quotient)

    # Round number of samples to the next integer
    raw_data = np.zeros((len(picks), int(duration * info['sfreq'] + 0.5)))
    raw = RawArray(raw_data, info)
    add_chpi(raw, dev_head_pos)
    quats = _calculate_chpi_positions(
        raw, t_step_min=raw.info['sfreq'] * head_pos_sfreq_quotient,
        t_step_max=raw.info['sfreq'] * head_pos_sfreq_quotient, t_window=1.0)
    _assert_quats(quats, dev_head_pos, dist_tol=0.001, angle_tol=1.,
                  vel_atol=4e-3)  # 4 mm/s

def _calculate_chpi_coil_locs(raw, verbose):
    """Wrap to facilitate change diff."""
    chpi_amplitudes = compute_chpi_amplitudes(raw, verbose=verbose)
    chpi_locs = compute_chpi_locs(raw.info, chpi_amplitudes, verbose=verbose)
    return _chpi_locs_to_times_dig(chpi_locs)

def _check_dists(info, cHPI_digs, n_bad=0, bad_low=0.02, bad_high=0.04):
    __tracebackhide__ = True
    orig = _get_hpi_initial_fit(info)
    hpi_coil_distances = cdist(orig, orig)
    new_pos = np.array([d['r'] for d in cHPI_digs])
    mask, distances = _compute_good_distances(hpi_coil_distances, new_pos)
    good_idx = np.where(mask)[0]
    assert len(good_idx) >= 3
    meds = np.empty(len(orig))
    for ii in range(len(orig)):
        idx = np.setdiff1d(good_idx, ii)
        meds[ii] = np.median(distances[ii][idx])
    meds = np.array(meds)
    assert_array_less(meds[good_idx], 0.003)
    bad_idx = np.where(~mask)[0]
    if len(bad_idx):
        bads = meds[bad_idx]
        assert_array_less(bad_low, bads)
        assert_array_less(bads, bad_high)

@pytest.mark.slowtest
@testing.requires_testing_data
def test_calculate_chpi_coil_locs_artemis():
    """Test computing just cHPI locations."""
    raw = read_raw_fif(chpi_fif_fname, allow_maxshield='yes', preload=True)
    # This is a little hack (aliasing while decimating) to make it much faster
    # for testing purposes only. We can relax this later if we find it breaks
    # something.
    raw_dec = _decimate_chpi(raw, 15)
    times, cHPI_digs = _calculate_chpi_coil_locs(raw_dec, verbose='debug')

    # spot check
    assert_allclose(times[0], 9., atol=1e-2)
    assert_allclose(cHPI_digs[0][2]['r'],
                    [-0.01937833, 0.00346804, 0.06331209], atol=1e-3)
    assert_allclose(cHPI_digs[0][2]['gof'], 0.9957, atol=1e-3)

    assert_allclose(cHPI_digs[0][4]['r'],
                    [-0.0655, 0.0755, 0.0004], atol=3e-3)
    assert_allclose(cHPI_digs[0][4]['gof'], 0.9323, atol=1e-3)
    _check_dists(raw.info, cHPI_digs[0], n_bad=1)

    # test on 5k artemis data
    raw = read_raw_artemis123(art_fname, preload=True)
    times, cHPI_digs = _calculate_chpi_coil_locs(raw, verbose='debug')

    assert len(np.setdiff1d(times, raw.times + raw.first_time)) == 0
    assert_allclose(times[5], 1.5, atol=1e-3)
    assert_allclose(cHPI_digs[5][0]['gof'], 0.995, atol=5e-3)
    assert_allclose(cHPI_digs[5][0]['r'],
                    [-0.0157, 0.0655, 0.0018], atol=1e-3)
    _check_dists(raw.info, cHPI_digs[5])
    coil_amplitudes = compute_chpi_amplitudes(raw)
    with pytest.raises(ValueError, match='too_close'):
        compute_chpi_locs(raw.info, coil_amplitudes, too_close='foo')
    # ensure values are in a reasonable range
    amps = np.linalg.norm(coil_amplitudes['slopes'], axis=-1)
    amps /= coil_amplitudes['slopes'].shape[-1]
    assert amps.shape == (len(coil_amplitudes['times']), 3)
    assert_array_less(amps, 1e-11)
    assert_array_less(1e-13, amps)
    # with nan amplitudes (i.e., cHPI off) it should return an empty array,
    # but still one that is 3D
    coil_amplitudes['slopes'].fill(np.nan)
    chpi_locs = compute_chpi_locs(raw.info, coil_amplitudes)
    assert chpi_locs['rrs'].shape == (0, 3, 3)
    pos = compute_head_pos(raw.info, chpi_locs)
    assert pos.shape == (0, 10)

def assert_suppressed(new, old, suppressed, retained):
    """Assert that some frequencies are suppressed and others aren't."""
    __tracebackhide__ = True
    from scipy.signal import welch
    picks = pick_types(new.info, meg='grad')
    sfreq = new.info['sfreq']
    new = new.get_data(picks)
    old = old.get_data(picks)
    f, new = welch(new, sfreq, 'hann', nperseg=1024)
    _, old = welch(old, sfreq, 'hann', nperseg=1024)
    new = np.median(new, axis=0)
    old = np.median(old, axis=0)
    for freqs, lim in ((suppressed, (10, 60)), (retained, (-3, 3))):
        for freq in freqs:
            fidx = np.argmin(np.abs(f - freq))
            this_new = np.median(new[fidx])
            this_old = np.median(old[fidx])
            suppression = -10 * np.log10(this_new / this_old)
            assert lim[0] < suppression < lim[1], freq

@testing.requires_testing_data
def test_chpi_subtraction_filter_chpi():
    """Test subtraction of cHPI signals."""
    raw = read_raw_fif(chpi_fif_fname, allow_maxshield='yes', preload=True)
    raw.info['bads'] = ['MEG0111']
    raw.del_proj()
    raw_orig = raw.copy().crop(0, 16)
    with catch_logging() as log:
        filter_chpi(raw, include_line=False, t_window=0.2, verbose=True)
    log = log.getvalue()
    assert 'No average EEG' not in log
    assert '5 cHPI' in log
    # MaxFilter doesn't do quite as well as our algorithm with the last bit
    raw.crop(0, 16)
    # remove cHPI status chans
    raw_c = read_raw_fif(sss_hpisubt_fname).crop(0, 16).load_data()
    raw_c.pick_types(
        meg=True, eeg=True, eog=True, ecg=True, stim=True, misc=True)
    assert_meg_snr(raw, raw_c, 143, 624)
    # cHPI suppressed but not line freqs (or others)
    assert_suppressed(raw, raw_orig, np.arange(83, 324, 60), [30, 60, 150])
    raw = raw_orig.copy()
    with catch_logging() as log:
        filter_chpi(raw, include_line=True, t_window=0.2, verbose=True)
    log = log.getvalue()
    assert '5 cHPI' in log
    assert '6 line' in log
    # cHPI and line freqs suppressed
    suppressed = np.sort(np.concatenate([
        np.arange(83, 324, 60), np.arange(60, 301, 60),
    ]))
    assert_suppressed(raw, raw_orig, suppressed, [30, 150])

    # No HPI information
    raw = read_raw_fif(sample_fname, preload=True)
    raw_orig = raw.copy()
    assert raw.info['line_freq'] is None
    with pytest.raises(RuntimeError, match='line_freq.*consider setting it'):
        filter_chpi(raw, t_window=0.2)
    raw.info['line_freq'] = 60.
    with pytest.raises(ValueError, match='No appropriate cHPI information'):
        filter_chpi(raw, t_window=0.2)
    # but this is allowed
    with catch_logging() as log:
        filter_chpi(raw, t_window='auto', allow_line_only=True, verbose=True)
    log = log.getvalue()
    assert '0 cHPI' in log
    assert '1 line' in log
    # Our one line freq suppressed but not others
    assert_suppressed(raw, raw_orig, [60], [30, 45, 75])

    # When MaxFliter downsamples, like::
    #     $ maxfilter -nosss -ds 2 -f test_move_anon_raw.fif \
    #           -o test_move_anon_ds2_raw.fif
    # it can strip out some values of info, which we emulate here:
    raw = read_raw_fif(chpi_fif_fname, allow_maxshield='yes')
    raw = raw.crop(0, 1).load_data().resample(600., npad='auto')
    raw.info['lowpass'] = 200.
    del raw.info['maxshield']
    del raw.info['hpi_results'][0]['moments']
    del raw.info['hpi_subsystem']['event_channel']
    with catch_logging() as log:
        filter_chpi(raw, t_window='auto', verbose=True)
    with pytest.raises(ValueError, match='must be > 0'):
        filter_chpi(raw, t_window=-1)
    assert '2 cHPI' in log.getvalue()

@testing.requires_testing_data
def test_calculate_head_pos_ctf():
    """Test extracting of cHPI positions from CTF data."""
    raw = read_raw_ctf(ctf_chpi_fname)
    chpi_locs = extract_chpi_locs_ctf(raw)
    quats = compute_head_pos(raw.info, chpi_locs)
    mc_quats = read_head_pos(ctf_chpi_pos_fname)
    mc_quats[:, 9] /= 10000  # had old factor in there twice somehow...
    _assert_quats(quats, mc_quats, dist_tol=0.004, angle_tol=2.5, err_rtol=1.,
                  vel_atol=7e-3)  # 7 mm/s
    plot_head_positions(quats, info=raw.info)

    raw = read_raw_fif(ctf_fname)
    with pytest.raises(RuntimeError, match='Could not find'):
        extract_chpi_locs_ctf(raw)

@testing.requires_testing_data
def test_calculate_head_pos_kit():
    """Test calculation of head position using KIT data."""
    raw = read_raw_kit(con_fname, mrk_fname, elp_fname, hsp_fname)
    assert len(raw.info['hpi_results']) == 1
    chpi_locs = extract_chpi_locs_kit(raw)
    assert chpi_locs['rrs'].shape == (2, 5, 3)
    assert_array_less(chpi_locs['gofs'], 1.)
    assert_array_less(0.98, chpi_locs['gofs'])
    quats = compute_head_pos(raw.info, chpi_locs)
    assert quats.shape == (2, 10)
    # plotting works
    plot_head_positions(quats, info=raw.info)
    raw_berlin = read_raw_kit(berlin_fname)
    assert_allclose(raw_berlin.info['dev_head_t']['trans'], np.eye(4))
    assert len(raw_berlin.info['hpi_results']) == 0
    with pytest.raises(ValueError, match='Invalid value'):
        extract_chpi_locs_kit(raw_berlin)
    with pytest.raises(RuntimeError, match='not find appropriate'):
        extract_chpi_locs_kit(raw_berlin, 'STI 014')
    with pytest.raises(RuntimeError, match='no initial cHPI'):
        compute_head_pos(raw_berlin.info, chpi_locs)

ext: py

max_stars_repo_name: GanshengT/mne-python

max_issues_repo_name: GanshengT/mne-python

max_forks_repo_name: GanshengT/mne-python

reason: import numpy,from numpy,from scipy


text:
import pytest
import tensorflow as tf
import numpy as np

@pytest.fixture(scope='session', name='graph_tuple')
def simple_graph_and_outputs():
    graph = tf.Graph()
    with graph.as_default():
        weight = tf.constant(np.random.randn(3, 3),
                             dtype=tf.float32,
                             name='weight')
        bias = tf.constant(np.pi,
                           dtype=tf.float32,
                           name='bias')
        x1 = tf.constant(np.random.randn(10, 3),
                         dtype=tf.float32,
                         name='x1')
        x2 = tf.add(tf.matmul(x1, weight), bias, name='x2')
        bias2 = tf.constant(3.69, name='bias2', dtype=tf.float32)
        x3 = tf.multiply(x2, bias2, name='x3')
    return graph.as_graph_def(), [x2.op.name, x3.op.name]

reason: import numpy


text:
import numpy as np
from flask import Flask, request, jsonify, render_template
import pickle

from predict import predict_out

app = Flask(__name__)

@app.route('/')
def home():
    return render_template('index.html')

@app.route('/predict',methods=['POST'])
def predict():
    rate = request.form['rate']
    output = predict_out(rate)[0][0]

    output = np.round(output)

    return render_template('index.html', prediction_text='Predicted Sales should be: $ {}K'.format(output))

if __name__ == "__main__":
    app.run(debug=True)


text:
import numpy as np

from tests.test_mcts_player import FirstChoiceHeuristic
from zero_play.play_controller import PlayController, PlayerResults
from zero_play.game_state import GridGameState
from zero_play.mcts_player import MctsPlayer

class FirstPlayerWinsGame(GridGameState):
    game_name = 'First Player Wins'

    def __init__(self,
                 board_height: int = 1,
                 board_width: int = 1,
                 spaces: np.ndarray = None):
        super().__init__(board_height, board_width, spaces=spaces)

    def is_win(self, player: int) -> bool:
        if (self.board == 0).sum() != 0:
            return False
        x_count = (self.board == self.X_PLAYER).sum()
        o_count = (self.board == self.O_PLAYER).sum()
        if player == self.X_PLAYER:
            return o_count < x_count
        return x_count <= o_count

class SecondPlayerWinsGame(FirstPlayerWinsGame):
    game_name = 'Second Player Wins'

    def __init__(self,
                 board_height: int = 1,
                 board_width: int = 2,
                 spaces: np.ndarray = None):
        super().__init__(board_height, board_width, spaces=spaces)

class NoPlayerWinsGame(FirstPlayerWinsGame):
    game_name = 'No Player Wins'

    def is_win(self, player: int):
        return False

def test_play_one_first_wins():
    start_state = FirstPlayerWinsGame()
    controller = PlayController(start_state=start_state,
                                players=[MctsPlayer(start_state),
                                         MctsPlayer(start_state)])
    expected_p1_wins = 1
    expected_ties = 0
    expected_p2_wins = 0

    p1_wins, ties, p2_wins = controller.play()

    assert expected_p1_wins == p1_wins
    assert expected_ties == ties
    assert expected_p2_wins == p2_wins

def test_play_one_second_wins():
    start_state = SecondPlayerWinsGame()
    controller = PlayController(start_state=start_state,
                                players=[MctsPlayer(start_state),
                                         MctsPlayer(start_state)])
    expected_p1_wins = 0
    expected_ties = 0
    expected_p2_wins = 1

    p1_wins, ties, p2_wins = controller.play()

    assert expected_p1_wins == p1_wins
    assert expected_ties == ties
    assert expected_p2_wins == p2_wins

def test_play_one_tie():
    start_state = NoPlayerWinsGame()
    controller = PlayController(start_state=start_state,
                                players=[MctsPlayer(start_state),
                                         MctsPlayer(start_state)])
    expected_p1_wins = 0
    expected_ties = 1
    expected_p2_wins = 0

    p1_wins, ties, p2_wins = controller.play()

    assert expected_p1_wins == p1_wins
    assert expected_ties == ties
    assert expected_p2_wins == p2_wins

def test_play_two():
    start_state = FirstPlayerWinsGame()
    controller = PlayController(start_state=start_state,
                                players=[MctsPlayer(start_state),
                                         MctsPlayer(start_state)])
    expected_p1_wins = 2
    expected_ties = 0
    expected_p2_wins = 0

    p1_wins, ties, p2_wins = controller.play(games=2)

    assert expected_p1_wins == p1_wins
    assert expected_ties == ties
    assert expected_p2_wins == p2_wins

def test_play_two_flip():
    start_state = FirstPlayerWinsGame()
    controller = PlayController(start_state=start_state,
                                players=[MctsPlayer(start_state),
                                         MctsPlayer(start_state)])
    expected_p1_wins = 1
    expected_ties = 0
    expected_p2_wins = 1

    p1_wins, ties, p2_wins = controller.play(games=2, flip=True)

    assert expected_p1_wins == p1_wins
    assert expected_ties == ties
    assert expected_p2_wins == p2_wins

def test_display(capsys):
    game = SecondPlayerWinsGame()
    heuristic = FirstChoiceHeuristic()
    players = [MctsPlayer(game,
                          game.X_PLAYER,
                          iteration_count=10,
                          heuristic=heuristic),
               MctsPlayer(game,
                          game.O_PLAYER,
                          iteration_count=20,
                          heuristic=heuristic)]
    controller = PlayController(game, players)
    expected_output = """\
  AB
1 ..

  AB
1 X.

  AB
1 XO

  AB
1 ..

  AB
1 X.

  AB
1 XO

('mcts', 'first choice', '10 iterations') - 1 wins,
('mcts', 'first choice', '20 iterations') - 1 wins,
0 ties
"""

    controller.play(games=2, flip=True, display=True)

    out, err = capsys.readouterr()
    expected_lines = expected_output.splitlines()
    lines = out.splitlines()
    assert lines[:-3] == expected_lines[:-3]
    for line, expected_line in zip(lines[-3:], expected_lines[-3:]):
        assert line.startswith(expected_line)

def test_player_results():
    player = MctsPlayer(FirstPlayerWinsGame(), iteration_count=100)
    player_results = PlayerResults(player)
    player_results.move_count = 100
    player_results.total_time = 1246  # seconds
    player_results.win_count = 1

    summary = player_results.get_summary()

    assert summary == "('mcts', 'playout', '100 iterations') - 1 wins, 12.5s/move"

max_stars_repo_name: donkirkby/zero-play

max_issues_repo_name: donkirkby/zero-play

max_forks_repo_name: donkirkby/zero-play


text:
import os
import errno
import itertools
import scipy.misc
import numpy as np

curr_dir = os.getcwd()
proj_dir = os.path.normpath(os.path.join(curr_dir))
image_dir = 'data/train'
input_filepath = os.path.normpath(os.path.join(proj_dir, image_dir))

validation_dir = 'data/processed/validation'
train_dir = 'data/processed/train'

class_dict = {'dog': 0, 'cat': 1}

def put_classes_into_separate_folders(parent_dir, images):
    make_sure_path_exists(parent_dir)

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

train_shuffled_images = []
test_shuffled_images = []
def create_train_validation_shuffle_images():
    test_percentage = 0.2
    for key, _ in class_dict.items():
        class_images = [x for x in os.listdir(input_filepath) if x.startswith(key + ".")]
        k = int(len(class_images) * test_percentage)
        test_shuffled_images = class_images[0:k]
        train_shuffled_images = class_images[k:]
    return train_shuffled_images, test_shuffled_images

train_images = []
validation_images = []
def resize_of_images():
    train_images = [scipy.misc.imresize(scipy.misc.imread(input_filepath + '/' + image), (100,50)) for image in train_shuffled_images]
    validation_images = [scipy.misc.imresize(scipy.misc.imread(input_filepath +
        '/' + image), (100,50)) for image in test_shuffled_images]
    return train_images, validation_images

def put_class_images_in_folders(save_dir, image_files, class_feature):
    counter = 0
    class_counter = 0
    first = True

    make_sure_path_exists(save_dir)

    for image in image_files:
        counter += 1

        if (counter % int(len(image_files) / len(class_feature)) == 0) and (first == False) != (counter == int(len(image_files))):
            class_counter += 1
        first = False

        class_dir = os.path.join(save_dir,class_feature[class_counter])
        make_sure_path_exists(class_dir)

        save_image = os.path.join(class_dir,'{}_{}.jpeg'.format(class_feature[class_counter],"".join((map(str,np.random.randint(0,9,8))))))

        scipy.misc.imsave(save_image, image)

if __name__ == '__main__':
    train_shuffled_images, test_shuffled_images = create_train_validation_shuffle_images()
    train_images, test_images = resize_of_images()
    keys = []
    for k, _ in class_dict.items():
        keys.append(k)
    print(input_filepath)
    put_class_images_in_folders(input_filepath + '/train/', train_images,
            keys)
    put_class_images_in_folders(input_filepath + '/validation/', test_images,
            keys)

max_stars_repo_name: eleijonmarck/image-preprocessing

max_issues_repo_name: eleijonmarck/image-preprocessing

max_forks_repo_name: eleijonmarck/image-preprocessing

reason: import numpy,import scipy


text:
"""Utility functions for measuring performances."""
import numpy as np
from time import time as tst

def timeit(method):
    """Get computing time of a function.

    @timeit
    def my_fcn()
    """
    def timed(*args, **kw):
        ts = tst()
        result = method(*args, **kw)
        te = tst()
        print(f"{method.__name__} : {(te - ts) * 1000}ms")
        return result
    return timed

def id(x):
    """Get the memory block address of an array."""
    return x.__array_interface__['data'][0]

def get_data_base(arr):
    """For a given array, finds the base array that "owns" the actual data."""
    base = arr
    while isinstance(base.base, np.ndarray):
        base = base.base
    return base

def arrays_share_data(x, y):
    """Return if two arrays share an offset."""
    return get_data_base(x) is get_data_base(y)

max_stars_repo_name: StanSStanman/frites

max_issues_repo_name: StanSStanman/frites

max_forks_repo_name: StanSStanman/frites


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
# CODE NAME HERE

# CODE DESCRIPTION HERE

Created on 2020-05-21

@author: cook
"""
import numpy as np
from astropy.io import fits
from  pathlib import Path
from typing import Tuple, Dict

# =============================================================================
# Define variables
# =============================================================================
# define workspace
WORKSPACE = '/spirou/cfht_nights/common/raw/'
WORKSPACE = '/spirou/cfht_nights/common/rawsym/'

# use header
USE_HDR = True
# skip engineering
SKIP_ENGINEERING = True

# define required file types
REQ_FILES = dict()
REQ_FILES['DARK_DARK'] = 'd.fits', 1
REQ_FILES['DARK_FLAT'] = 'f.fits', 3
REQ_FILES['FLAT_DARK'] = 'f.fits', 3
REQ_FILES['FLAT_FLAT'] = 'f.fits', 3
REQ_FILES['FP_FP'] = 'a.fits', 3
REQ_FILES['HC_HC'] = 'c.fits', 1

# a more in depth requirement
# SBCCAS_P, SBCREF_P, EXT, MIN NUMBER
REQ_HDR = dict()
REQ_HDR['DARK_DARK'] = 'pos_pk', 'pos_pk', 'd.fits', 1
REQ_HDR['DARK_FLAT'] = 'pos_pk', 'pos_wl', 'f.fits', 3
REQ_HDR['FLAT_DARK'] = 'pos_wl', 'pos_pk', 'f.fits', 3
REQ_HDR['FLAT_FLAT'] = 'pos_wl', 'pos_wl', 'f.fits', 3
REQ_HDR['FP_FP'] = 'pos_fp', 'pos_fp', 'a.fits', 3
REQ_HDR['HC_HC'] = 'pos_hc1', 'pos_hc1', 'c.fits', 1
# define the header keys to look at (related to values in REQ_HDR)
KEY1 = 'SBCCAS_P'
KEY2 = 'SBCREF_P'

# =============================================================================
# Define functions
# =============================================================================
def valid_night(night_path: Path) -> Tuple[bool, bool, bool, Dict]:
    """

    :param night_path: Path

    :return: passed_exist, passed_number and count dictionary
    """

    # deal with engineering
    if SKIP_ENGINEERING:
        obj_files = list(night_path.glob('*o.fits'))
        if len(obj_files) == 0:
            return True, True, True, dict()

    passed_exist = True
    passed_number = True

    count = dict()

    for req_file in REQ_FILES:
        # get extension and number
        ext, num = REQ_FILES[req_file]
        # get files
        files = np.sort(list(night_path.glob('*' + ext)))
        # if use header filter these files by header values
        if USE_HDR:
            files = filter_by_header(files, req_file)
        # get number of files
        numfiles = len(files)

        # deal with files existing
        if numfiles > 0:
            passed_exist = True
            # add to count dictionary
            count[req_file] = [numfiles, num, True]
        else:
            passed_exist = False
            # add to count dictionary
            count[req_file] = [numfiles, num, False]
        # deal with minimum number of files
        if numfiles < num:
            passed_number = False
            # add to count dictionary
            count[req_file] = [numfiles, num, False]
        else:
            passed_number = True
            # add to count dictionary
            count[req_file] = [numfiles, num, True]
    # return passed criteria and count
    return passed_exist, passed_number, False, count

def filter_by_header(files, key):
    # store filtered files
    filtered_files = []
    # loop through files
    for filename in files:
        # try to load file
        try:
            header = fits.getheader(filename)
        except:
            continue
        # if keys not in header skip this file (it is not valid)
        if KEY1 not in header:
            continue
        if KEY2 not in header:
            continue
        # if keys not equal to the required skip this file
        if not header[KEY1] == REQ_HDR[key][0]:
            continue
        if not header[KEY2] == REQ_HDR[key][1]:
            continue
        # add to filtered files
        filtered_files.append(filename)
    # return filtered files
    return filtered_files

# =============================================================================
# Start of code
# =============================================================================
if __name__ == "__main__":
    # get a list of night names
    files = np.sort(list(Path(WORKSPACE).glob('*')))
    # only keep directories
    directories = []
    for filename in files:
        if filename.is_dir():
            directories.append(filename)
    # sort directories
    directories = np.sort(directories)
    # storage
    list_valid = []
    list_invalid = []
    list_engineering = []
    # for all directories find those that are valid
    for directory in directories:

        # get validity
        p1, p2, p3, cnt = valid_night(directory)

        # deal with engineering nights
        if p3:
            print('Processing directory: {0}'.format(directory))
            print('\tDirectory is engineering night - skipping')
            list_engineering.append(directory)
        # deal with
        elif not p1:
            # print directory
            print('Processing directory: {0}'.format(directory))
            print('\tDirectory not valid for DRS: {0}'.format(directory))
            for key in cnt:
                if not cnt[key][2]:
                    pargs = [key, cnt[key][0], cnt[key][1]]
                    print('\t\t {0}: not found'.format(*pargs))
            list_invalid.append(directory)
        elif not p2:
            print('Processing directory: {0}'.format(directory))
            print('\tMinimum number of files not met')
            for key in cnt:
                if not cnt[key][2]:
                    pargs = [key, cnt[key][0], cnt[key][1]]
                    print('\t\t {0}: found {1} required {2}'.format(*pargs))
            list_invalid.append(directory)
        else:
            list_valid.append(directory)

# =============================================================================
# End of code
# =============================================================================

reason: import numpy,from astropy


text:
from sklearn import linear_model
from sklearn import neural_network
import matplotlib.pyplot as plt
import numpy as np

def plot(model, title, test_data_x):
    test_data_y = model.predict(test_data_x)
    print(test_data_y)

    row_idx = 0
    for features in test_data_x:
        format = 'ro'
        if test_data_y[row_idx] == 1:
            format = 'bo'

        plt.plot(features[0], features[1], format)
        row_idx += 1

    plt.title(title)
    plt.show()

# Can't fix xor data with linear regression directly
# since xor is not linear itself, We need to choose
# a feature map which can be done by neural network?
def regression(train_x, train_y):
    regr = linear_model.LinearRegression()
    regr.fit(train_x, train_y)
    print(regr.coef_)
    plot(regr, 'Linear regression', train_x)

if __name__ == '__main__':
    # xor data
    train_data_x = np.array([
        [0, 0],
        [0, 1],
        [1, 0],
        [1, 1]])

    train_data_y = np.array([0, 1, 1, 0])

    regression(train_data_x, train_data_y)

    mlp = neural_network.MLPRegressor()
    mlp.fit(train_data_x, train_data_y)
    y_pred = mlp.predict(train_data_x)
    print(y_pred)

max_stars_repo_name: tqtifnypmb/ML

max_issues_repo_name: tqtifnypmb/ML

max_forks_repo_name: tqtifnypmb/ML


text:
# -*- coding: utf-8 -*-
"""
Created on Mon Apr  6 17:03:33 2020

@author: Leonardo Saccotelli
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc
from Algoritmi_Quadratura import CompositeTrapezoid

# Funzione test da integrare
def f(x):
    y = 1/x - 1/(x**2)
    return y

# Primitiva della funzione test
def F(x):
    y = np.log(abs(x))+1/x
    return y

# Intervallo di integrazione e integrale esatto
a =2.0; b = 4.0
I = F(b) - F(a)

# Formula trapezi composti al variare di N
N_max = 1000

N_range = range(2,N_max,10)

Err_Trap = np.zeros(len(N_range))

k = 0

for N in N_range:
    TC = CompositeTrapezoid(f,a,b,N)
    Err_Trap[k] = abs(I-TC)
    k = k + 1

# Grafico errori al variare di N
plt.figure(1,figsize=(14,8))
plt.semilogy(N_range,Err_Trap,'k--o',label='Trap. Comp.')
plt.xlabel('N',fontsize=26)
plt.ylabel('Errore',fontsize=26)
plt.legend(prop={'size':26})
plt.title('Errore formule di quadratura composte al variare di N',fontsize=26)
rc('xtick',labelsize=26)
rc('ytick',labelsize=26)


text:
import time
from matplotlib import pyplot as plt
import numpy as np
from sklearn.metrics.pairwise import cosine_similarity
import random
from rnngen.word2vec.embedding_testing import test_embeddings
from rnngen.setupvariables import preprocess_setup, word2vec_setup, test_emb_setup
from rnngen.recurrentnetwork.trainer import load_data

def get_data(data, window_size):
    """
    Creates input training data by looking at words in context.

    If window_size=2 in the sentence "the black cat took all my cookies.
    If we are to find the context words to cat, these will be: the, black, took, all.
    These are xs and the y for all of them is cat.

    :param data: All text to train word2vec on
    :param window_size: How many words to look on beside the chosen one
    :return: training data, one y for every x.
    """
    xs = []
    ys = []
    for sentence in data:
        for ind in range(len(sentence)):
            indexes = list(range(max(ind-window_size, 0), max(ind, 0))) \
                      + list(range(min(ind+1, len(sentence)), min(ind + window_size+1, len(sentence))))
            for y in indexes:
                xs.append(sentence[ind])
                ys.append(sentence[y])
    return np.array(xs), np.array(ys)

def forward_propagation(x_batch, weights):
    # x_batch: batches x emb
    # y_batch: batches x vocab_size
    # weights: emb x vocab_size
    z = np.dot(x_batch, weights)
    softmax = np.exp(z) / (np.sum(np.exp(z)))
    return softmax

def back_propagation(y_batch, softmax, weights, x_batch):
    # Softmax:       64, 4036
    # Weights:       50, 4036
    # Batch Emb:     64, 50
    # d_embeddings:  4036, 50

    d_softmax = softmax - y_batch
    d_weights = np.dot(x_batch.T, d_softmax)
    d_embeddings = np.dot(d_softmax, weights.T)

    return d_weights, d_embeddings

def calculate_loss(softmax, y):
    m = softmax.shape[1]
    cost = -(1 / m) * np.sum(np.sum(y * np.log(softmax)))
    return cost

def one_hot(ys, unique_words):
    # Backprop needs ys in sparse vectors in order to function
    m = len(ys)
    one_hot_y = np.zeros((m, unique_words))
    for row in range(m):
        one_hot_y[row][int(ys[row])] = 1
    return one_hot_y

def get_words_to_check(dicts, test_words, random_words):
    words_to_check = []
    words_to_check += test_words
    vocab_size = len(dicts['word_to_id'])
    for _ in range(random_words):
        try:
            random_words1 = dicts['id_to_word'][random.randint(0, vocab_size-4)]
            random_words2 = dicts['id_to_word'][random.randint(0, vocab_size-4)]
            words_to_check.append((random_words1, random_words2))
        except KeyError:
            pass
    return words_to_check

def print_word_similarity(words_to_check, word_emb, dicts):
    checked_words = []
    for word1, word2 in words_to_check:
        try:
            sim = cosine_similarity(word_emb[dicts['word_to_id'][word1]].reshape(1, -1),
                                    word_emb[dicts['word_to_id'][word2]].reshape(1, -1))[0][0]
            print(f'{word1} | {word2}:   {round(sim, 4)}')
        except KeyError:
            pass
    print('')
    return checked_words

def get_batches(x, y, batches, word_emb, unique_words, iteration):
    x_indexes = x[batches * iteration:batches * iteration + batches]
    x_batch = word_emb[x_indexes.astype(int), :]
    y_batch = one_hot(y[batches * iteration:batches * iteration + batches], unique_words)
    return x_indexes, x_batch, y_batch

def count_print_loss(loss, losses, iters_before_reset, m, batches, epoch, epochs, iteration):
    loss = round(loss / iters_before_reset * 100, 4)
    losses.append(loss)
    print_losses = losses
    if len(losses) > 40:
        print_losses = losses[::3]
    if len(losses) > 200:
        print_losses = losses[::10]
    if len(losses) > 600:
        print_losses = losses[::20]
    if len(losses) > 1200:
        print_losses = losses[::60]
    print('Loss: {}'.format(loss), print_losses)
    print('Iter: {} of {}\nEpoch: {} of {}'.format(iteration, m // batches, epoch + 1, epochs))
    return losses

def word2vec_trainer(text_dir, save_emb_dir, previously_trained_emb=''):
    """
    Trains word embeddings on data with 2 layer neural network.

    :param text_dir: Directory to training data
    :param save_emb_dir: Directory to save word embeddings
    :param previously_trained_emb: If set to a directory, will open this emb and train it
    :return: Returns matrix of word embeddings
    """
    if previously_trained_emb and not previously_trained_emb.endswith('.npy'):
        previously_trained_emb += '.npy'
    data, dicts, sentences_lengths = load_data(text_dir, preprocess_setup)
    print('\nWord2Vec trainer starting...')
    batches = word2vec_setup['BATCHES']
    embedding_size = word2vec_setup['EMBEDDING_SIZE']
    epochs = word2vec_setup['EPOCHS']
    lr = word2vec_setup['LEARNING_RATE']
    short = word2vec_setup['SHORT_MODE']
    iters_before_decrease = word2vec_setup['ITERATIONS_BEFORE_LR_DECREASE']
    lr_decrease = word2vec_setup['LR_DECREASE']
    window_size = word2vec_setup['WINDOW_SIZE']

    use_test_embeddings = test_emb_setup['USE_TEST_EMBEDDINGS']
    verbose_cosine_distance = test_emb_setup['VERBOSE_COSINE']
    test_words = test_emb_setup['TESTING_WORDS']

    unique_words = len(dicts['word_to_id'])
    # Initialize or load word embedding matrix
    if not previously_trained_emb:
        word_emb = np.random.randn(unique_words, embedding_size)
    else:
        word_emb = np.load(previously_trained_emb)
    weights = np.random.randn(embedding_size, unique_words)
    xs, ys = get_data(data, window_size)
    m = len(xs)
    losses = []  # saves accumulated losses
    all_losses = []  # saves all individual losses so it can be printed in matplot
    accumulated_loss = 0  # Accumulates losses before each LR_DECREASE
    time_one_decrease = time.time()

    assert ys.shape == xs.shape
    print('Word2vec is up running. (May take a few minutes before verbose)\n')
    for epoch in range(epochs):  # Iters epochs
        for iteration in range(m//batches-3):
            # Gets batch data
            x_indexes, x_batch, y_batch = get_batches(xs, ys,
                                                      batches, word_emb,
                                                      unique_words, iteration)

            softmax = forward_propagation(x_batch, weights)  # Forward propagates
            loss = calculate_loss(softmax, y_batch)  # Calculates loss
            accumulated_loss += loss
            d_weights, d_word_embeddings = back_propagation(y_batch,  # Back Propagation
                                                            softmax,
                                                            weights,
                                                            x_batch)
            # Every iters_before_decrease iter, losses and cosine similarities will be printed
            if (iteration % iters_before_decrease == 0) and (iteration != 0 and epochs != 0):
                print('Time per output: ', round(time.time() - time_one_decrease), ' Seconds')
                time_one_decrease = time.time()
                losses = count_print_loss(accumulated_loss, losses,
                                          iters_before_decrease, m,
                                          batches, epoch, epochs,
                                          iteration)
                accumulated_loss = 0
                lr = lr * lr_decrease  # Decrease learning rate
                if verbose_cosine_distance:
                    words_to_check = get_words_to_check(dicts, test_words, random_words=4)
                    print_word_similarity(words_to_check, word_emb, dicts)
                if iteration > 30000 or epoch >= 1:
                    np.save(save_emb_dir, word_emb)

            if short and iteration > 100:   # Only for developing, easier checking so things work out.
                np.save(save_emb_dir, word_emb)
                break

            # Standard SGD. Should be upgraded to Adam
            weights = weights - d_weights * lr / len(x_batch) / batches
            word_emb[x_indexes.astype(int), :] = word_emb[x_indexes.astype(int), :] \
                - (d_word_embeddings * lr / len(x_batch)) / batches

            all_losses.append(loss)

    np.save(save_emb_dir, word_emb)
    if use_test_embeddings:
        test_embeddings(word_emb, dicts['id_to_word'], dicts['word_to_id'], test_emb_setup)
    plt.plot(all_losses)
    plt.show()

    print('\nOBS: Word vectors must be separately trained for every new change. '
          '\nThey cannot be used for another set of text'
          ' or with another set of preprocessing parameters.\n')

max_stars_repo_name: gabrielpetersson/rnngen

max_issues_repo_name: gabrielpetersson/rnngen

max_forks_repo_name: gabrielpetersson/rnngen


text:
import pickle
import theano
import theano.tensor as T
from pylearn2.models import mlp
from pylearn2.training_algorithms import sgd
from pylearn2.training_algorithms.learning_rule import AdaGrad
from pylearn2.training_algorithms.learning_rule import Momentum
from pylearn2.training_algorithms.sgd import ExponentialDecay
from pylearn2.training_algorithms.sgd import LinearDecay
from pylearn2.training_algorithms.learning_rule import AdaDelta

from pylearn2.termination_criteria import EpochCounter
from pylearn2.datasets.dense_design_matrix import DenseDesignMatrix
from pylearn2.costs.cost import Cost
from pylearn2.costs.cost import DefaultDataSpecsMixin
import numpy as np
from random import randint
import os
import matplotlib.pyplot as plt
from features import logfbank

class XOR(DenseDesignMatrix):
    def __init__(self):
        self.class_names = ['0', '1']
        X = [[randint(0, 1), randint(0, 1)] for _ in range(1000)]
        y = []
        for a, b in X:
            if a + b == 1:
                y.append([0, 1])
            else:
                y.append([1, 0])
        X = np.array(X)
        y = np.array(y)
        super(XOR, self).__init__(X=X, y=y)

class NegativeLogLikelihoodCost(DefaultDataSpecsMixin, Cost):
    supervised = True

    def expr(self, model, data, **kwargs):
        space, source = self.get_data_specs(model)
        space.validate(data)

        inputs, targets = data

        outputs = model.fprop(inputs)

        loss = -(targets * T.log(outputs)).sum(axis=1)

        return loss.mean()

left_context = 10
right_context = 5
keywords = ["she", "had"]

def getDataForFrames(data, X, y):
    fbanks = data[0]
    segments = data[1]
    feat_cnt = len(fbanks[0])
    for start, end, class_id in segments:

        if randint(0, 1) == 0:
            for i in xrange(start, end):
                features = np.empty(0)
                for j in xrange(-left_context, right_context):
                    if i + j < 0 or i + j >= len(fbanks):
                        features = np.append(features, np.zeros(feat_cnt))
                    else:
                        features = np.append(features, fbanks[i + j])
                result = np.zeros(len(keywords) + 1)
                result[class_id] = 1.0

                X.append(features)
                y.append(result)

        if randint(0, 1) == 0:
            for i in xrange(start, end):
                features = np.empty(0)
                for j in xrange(-left_context, right_context):
                    if i + j < start or i + j >= end:
                        features = np.append(features, np.zeros(feat_cnt))
                    else:
                        features = np.append(features, fbanks[i + j])
                result = np.zeros(len(keywords) + 1)
                result[class_id] = 1.0

                X.append(features)
                y.append(result)

def getWindowedFeats(fbanks):
    X = []
    feat_cnt = len(fbanks[0])
    for i in xrange(len(fbanks)):
        features = np.empty(0)
        for j in xrange(-left_context, right_context):
            if i + j < 0 or i + j >= len(fbanks):
                features = np.append(features, np.zeros(feat_cnt))
            else:
                features = np.append(features, fbanks[i + j])
        X.append(features)
    return X

def getDataFromPath(path, maxFiles = 1e9):
    X = []
    y = []
    for file in os.listdir(path):
        maxFiles -= 1
        if maxFiles < 0:
            break
        data = np.load(path + "/" + file)
        getDataForFrames(data, X, y)

    X = np.array(X)
    y = np.array(y)

    return X, y

class SHEHAD(DenseDesignMatrix):
    feat_cnt = 0

    def __init__(self, path):
        self.class_names = ["she", "had", "filler"]
        X, y = getDataFromPath(path)
        # print X.shape
        self.feat_cnt = len(X[0])
        super(SHEHAD, self).__init__(X=X, y=y)

def test(model):
    X, y = getDataFromPath("/Users/evgeny/data/TEST")

    confusion = np.zeros([3, 3])

    ypred = np.log(model.fprop(theano.shared(np.array(X), name='inputs')).eval())

    cnt = np.zeros(3)

    for a, b in zip(ypred, y):
        pos = np.argmax(b)
        i = np.argmax(a)
        #for i in xrange(3):
        confusion[pos][i] += 1
        # cnt[pos] += 1

    for i in xrange(3):
        for j in xrange(3):
            print "%.0f" % confusion[i][j],
        print   # confusion[i][j] /= cnt[i]
    print confusion

# create hidden layer with 2 nodes, init weights in range -0.1 to 0.1 and add
# a bias with value 1

rng = 0.001

modelName = "2x128relu50epochs-v7-momentum.mdl"
debug = True

if debug or not os.path.exists(modelName):
    ds = SHEHAD("/Users/evgeny/data/TRAIN")
    vds = SHEHAD("/Users/evgeny/data/TEST")
    hidden_layer = mlp.RectifiedLinear(layer_name='hidden', dim=128, irange=0.001, init_bias=0)
    hidden_layer2 = mlp.RectifiedLinear(layer_name='hidden2', dim=128, irange=0.01, init_bias=0)
    hidden_layer3 = mlp.RectifiedLinear(layer_name='hidden3', dim=128, irange=0.01, init_bias=0)
    # create Softmax output layer
    output_layer = mlp.Softmax(3, 'output', irange=.1)
    # create Stochastic Gradient Descent trainer that runs for 400 epochs
    cost = NegativeLogLikelihoodCost()
    rule = Momentum(0.9)
    # rule = Momentum(0.9, True)
    # update_callbacks=ExponentialDecay(1 + 1e-5, 0.001)
    trainer = sgd.SGD(learning_rate=0.01, cost=cost,  batch_size=128, termination_criterion=EpochCounter(1000),
                      monitoring_dataset=vds, learning_rule=rule)
    layers = [hidden_layer, hidden_layer2, output_layer]
    # create neural net that takes two inputs
    ann = mlp.MLP(layers, nvis=ds.feat_cnt)

    trainer.setup(ann, ds)
    print trainer.cost
    # train neural net until the termination criterion is true

    iteration = 0

    while True:
        trainer.train(dataset=ds)
        ann.monitor.report_epoch()
        ann.monitor()
        if iteration % 10 == 0:
            if not debug:
                with open(modelName, 'wb') as f:
                    pickle.dump(ann, f)
        if not trainer.continue_learning(ann):
            break

        iteration += 1

    if not debug:
        with open(modelName, 'wb') as f:
            pickle.dump(ann, f)

else:
    with open(modelName) as f:
        ann = pickle.load(f)

#test(ann)

#exit(0)

window = 0.025
step = 0.01
nfilt = 40
fftsize = 512

def extractLogFBank(rate, sig):
    feats = logfbank(sig, rate, window, step, nfilt, fftsize, 0, None, 0)
    return feats

#sph2pipe = "/Users/evgeny/kaldi3/tools/sph2pipe_v2.5/sph2pipe"

#os.system(sph2pipe + " -f wav " + "SA1.WAV" + " tmp.wav")

def computeFile(model, path):
    # import extractFeats
    import scipy.io.wavfile as wav

    (rate, sig) = wav.read(path)

    fbanks = extractLogFBank(rate, sig)

    X = getWindowedFeats(fbanks)

    ypred = np.log(model.fprop(theano.shared(np.array(X), name='inputs')).eval())

    ypred = np.transpose(ypred)

    ig, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex=True)
    ax1.plot(ypred[0])
    ax2.plot(ypred[1])
    ax3.plot(ypred[2])
    plt.show()

computeFile(ann, "tmp2.wav")

max_stars_repo_name: notimesea/VADnnet

max_issues_repo_name: notimesea/VADnnet

max_forks_repo_name: notimesea/VADnnet

reason: import numpy,import scipy,import theano


text:
###################################################################################################
# Copyright (c) 2021 Jonas Nicodemus
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
#
###################################################################################################
#
# This file incorporates work and modifications to the originally published code
# according to the previous license by the following contributors under the following licenses
#
#   Copyright (c) 2022 Birgit Hillebrecht
#
#   Permission is hereby granted, free of charge, to any person obtaining a copy
#   of this software and associated documentation files (the "Software"), to deal
#   in the Software without restriction, including without limitation the rights
#   to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
#   copies of the Software, and to permit persons to whom the Software is
#   furnished to do so, subject to the following conditions:
#
#   The above copyright notice and this permission notice shall be included in all
#   copies or substantial portions of the Software.
#
#   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#   OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
#   SOFTWARE.
#
###################################################################################################

import logging
import os

import matplotlib.pyplot as plt
import numpy as np

def fig_size(width_pt, fraction=1, ratio=(5 ** .5 - 1) / 2, subplots=(1, 1)):
    """
    Returns the width and heights in inches for a matplotlib figure.

    :param float width_pt: document width in points, in latex can be determined with \showthe\linewidth
    :param float fraction: fraction of the width with which the figure will occupy
    :param float ratio: ratio of the figure, default is the golden ratio
    :param tuple subplots: the shape of subplots
    :return: float fig_width_in: width in inches of the figure, float fig_height_in: height in inches of the figure
    """
    # Width of figure (in pts)
    fig_width_pt = width_pt * fraction
    # Convert from pt to inches
    inches_per_pt = 1 / 72.27

    # Golden ratio to set aesthetic figure height
    golden_ratio = ratio

    # Figure width in inches
    fig_width_in = fig_width_pt * inches_per_pt
    # Figure height in inches
    fig_height_in = fig_width_in * golden_ratio * (subplots[0] / subplots[1])

    return fig_width_in, fig_height_in

def new_fig(width_pt=410, fraction=1, ratio=(5 ** .5 - 1) / 2, subplots=(1, 1)):
    """
    Creates new instance of a `matplotlib.pyplot.figure` fig by using the `fig_size` function.

    :param float width_pt: document width in points, in latex can be determined with \showthe\textwidth
    :param float fraction: fraction of the width with which the figure will occupy
    :param float ratio: ratio of the figure, default is the golden ratio
    :param tuple subplots: the shape of subplots
    :return: matplotlib.pyplot.figure fig: instance of a `matplotlib.pyplot.figure` with desired width and height
    """
    fig = plt.figure(figsize=fig_size(width_pt, fraction, ratio, subplots))
    return fig

def save_fig(fig, name, path, tight_layout=True):
    """
    Saves a `matplotlib.pyplot.figure` as pdf file.

    :param matplotlib.pyplot.figure fig: instance of a `matplotlib.pyplot.figure` to save
    :param str name: filename without extension
    :param str path: path where the figure is saved, if None the figure is saved at the results directory
    :param bool crop: bool if the figure is cropped before saving
    """
    if tight_layout:
        fig.tight_layout()

    if not os.path.exists(path):
        os.makedirs(path)

    fig.savefig(os.path.join(path, f'{name}.pdf'), transparent=True)

def plot_states_1D(T, Z_ref, name, filepath, ylabeli=r'$x(t)$', legend_ref =r'$x(t)$', legend_pred=r'$\hat{x}(t)$', Z_pred=None, Z_mpc=None):
    linewidth = 2
    fig = new_fig()
    ax = fig.add_subplot(1, 1, 1)
    ax.set(xlabel='Time $t$', ylabel=ylabeli)
    ax.set(xlim=[np.min(T), np.max(T)])
    colors = ['tab:blue', 'tab:orange', 'tab:red', 'tab:green']
    ax.plot(T, Z_ref, linewidth=linewidth, label=legend_ref, c=colors[2])
    if Z_pred is not None:
        ax.plot(T, Z_pred, linestyle='--', linewidth=linewidth, label=legend_pred,
                    c=colors[1])

    ax.grid('on')
    ax.legend(loc='best')
    fig.tight_layout()
    if filepath is not None:
        save_fig(fig, name, filepath)
    plt.show()

def plot_trustzone(T, Z, Z_err, Z_ref,name, filepath):
    Z_err = np.reshape(Z_err, Z.shape)
    linewidth = 2
    fig = new_fig()
    ax = fig.add_subplot(1, 1, 1)
    ax.set(xlabel='Time $t$', ylabel=r'x(t)')
    ax.set(xlim=[np.min(T), np.max(T)])
    colors = ['tab:blue', 'tab:orange', 'tab:red', 'tab:green']
    ax.plot(T, Z, linewidth=linewidth, c=colors[2], label=r'$\hat{x}(t)$')
    ax.plot(T, Z_ref, linewidth=linewidth, linestyle='--', c=colors[0], label=r'$x(t)$')
    ax.plot(T, Z-Z_err, linestyle='--', linewidth=linewidth, c=colors[1])
    ax.plot(T, Z+Z_err, linestyle='--', linewidth=linewidth, c=colors[1])

    ax.grid('on')
    ax.legend(loc='best')
    fig.tight_layout()
    save_fig(fig, name, filepath)
    plt.show()

def plot_absolute_errors(T, Z_pred_err, name, filepath, Z_ref=None,  Z_pred=None):
    colors = ['tab:orange', 'tab:red', 'tab:green', 'tab:blue']
    fig = new_fig()
    ax = fig.add_subplot(1, 1, 1)
    ax.set(xlabel='Time $t$ ', ylabel=r'Error')
    ax.set_yscale('log')
    ax.grid(True, which='both')
    ax.set(xlim=[np.min(T), np.max(T)])

    ax.plot(T, Z_pred_err[:,0]+Z_pred_err[:,1], linewidth=2, color= colors[0], label=r'$E_\mathrm{Init} + E_\mathrm{PI}$')
    ax.plot(T, Z_pred_err[:,0], linewidth=2, linestyle="--",  color= colors[3], label=r'$E_\mathrm{Init}$')
    ax.plot(T, Z_pred_err[:,1], linewidth=2, linestyle="--", color= colors[2], label=r'$E_\mathrm{PI}$')

    if Z_pred is not None and Z_ref is not None:
        logging.info(f'{Z_ref - Z_pred}')
        if  Z_ref.ndim>1 :
            abs_errors = (np.sum((Z_ref - Z_pred)**2, axis=1))**0.5
        else :
            abs_errors = np.absolute(Z_ref-Z_pred)
        ax.plot(T, abs_errors, linewidth=2, color= colors[1], label=r'$\|\hat{x}(t)-x(t)|$')

    ax.legend(loc='best')
    save_fig(fig, name, filepath)
    fig.tight_layout()
    plt.show()

def plot_loss(T, filename, filepath, xname, yname, Loss1, Loss2=None, Loss3= None, Loss4= None) :
    t_sorted = np.argsort(T[:,0])
    if Loss1.ndim > 1:
        Loss1 = Loss1[t_sorted, 0]
    else:
        Loss1 = Loss1[t_sorted]
    T_sort = T[t_sorted]

    axxmax = 1
    axymax = 1
    if Loss2 is not None:
        axxmax = 2
    if Loss3 is not None or Loss4 is not None:
        axymax = 2

    colors = ['tab:blue', 'tab:blue', 'tab:blue', 'tab:blue']
    fig = new_fig()
    ax = fig.add_subplot(axxmax, axymax, 1)
    ax.set( ylabel=yname)
    # ax.set_yscale('log')
    ax.grid(True, which='both')
    ax.set(xlim=[np.min(T), np.max(T)])

    ax.plot(T_sort, Loss1, linewidth=2,c=colors[0])

    if Loss2 is not None:
        Loss2 = Loss2[t_sorted]
        ax2= fig.add_subplot(axxmax,axymax,2)
        ax2.set(xlim=[np.min(T), np.max(T)])
        ax2.plot(T_sort, Loss2, linewidth=2, c=colors[1])
        if Loss4 is None:
            ax2.set(xlabel=xname)
    else:
        ax.set(xlabel=xname)

    if Loss3 is not None:
        Loss3 = Loss3[t_sorted]
        ax3= fig.add_subplot(axxmax,axymax,3)
        ax3.set(xlim=[np.min(T), np.max(T)])
        ax3.plot(T_sort, Loss3, linewidth=2, c=colors[2])
        ax3.set(xlabel=xname)

    if Loss4 is not None:
        Loss4 = Loss4[t_sorted]
        ax4= fig.add_subplot(axxmax,axymax,4)
        ax4.set(xlim=[np.min(T), np.max(T)])
        ax4.plot(T_sort, Loss4, linewidth=2, c=colors[3])
        ax4.set(xlabel=xname)

    save_fig(fig, filename, filepath)
    plt.show()

max_stars_repo_name: bhillebrecht/CertifiedML-ODE

max_issues_repo_name: bhillebrecht/CertifiedML-ODE

max_forks_repo_name: bhillebrecht/CertifiedML-ODE


text:
import numpy as np

t = int(input())
for i in range(t):
    rc = input().split()
    r1 = int(rc[0])
    c1 = int(rc[1])
    img1 = []
    for j in range(r1):
        img1.append(list(input()))

    rc = input().split()
    r2 = int(rc[0])
    c2 = int(rc[1])
    img2 = []
    for j in range(r2):
        img2.append(list(input()))

    img1 = np.array(img1)
    img2 = np.array(img2)
    if (r1*c1) < (r2*c2):
        changer = img1
        fixed = img2
        changerC = c1
        changerR = r1
        fixedC = c2
        fixedR = r2
    else:
        changer = img2
        fixed = img1
        changerC = c2
        changerR = r2
        fixedC = c1
        fixedR = r1

    changer2 = np.flip(changer, 0)
    changer11 = np.rot90(changer)
    changer12 = np.rot90(changer11)
    changer13 = np.rot90(changer12)
    changer21 = np.rot90(changer2)
    changer22 = np.rot90(changer21)
    changer23 = np.rot90(changer22)
    changes1 = [changer,changer12,changer2,changer22]
    changes2 = [changer11, changer13, changer21, changer23]

    fixed1 = fixed
    for j in range(changerR-1):
        fixed1 = np.vstack([['.']*fixedC, fixed1])
        fixed1 = np.vstack([fixed1, ['.']*fixedC])
    for j in range(changerC-1):
        fixed1 = np.insert(fixed1, 0, '.', axis=1)
        fixed1 = np.insert(fixed1, fixedC+1+j, '.', axis=1)

    answer = 0
    for j in range(fixedR+changerR-1):
        for k in range(fixedC+changerC-1):
            for z in changes1:
                slidedWindow = fixed1[j:j+changerR,  k:k+changerC]
                temp = np.sum((slidedWindow=='#') & (z=='#'))
                if temp>answer:
                    answer=temp

    fixed2 = fixed
    for j in range(changerC-1):
        fixed2 = np.vstack([['.']*fixedC, fixed2])
        fixed2 = np.vstack([fixed2, ['.']*fixedC])
    for j in range(changerR-1):
        fixed2 = np.insert(fixed2, 0, '.', axis=1)
        fixed2 = np.insert(fixed2, fixedC+1+j, '.', axis=1)
    for j in range(fixedR+changerC-1):
        for k in range(fixedC+changerR-1):
            for z in changes2:
                slidedWindow = fixed2[j:j+changerC,  k:k+changerR]
                temp = np.sum((slidedWindow=='#') & (z=='#'))
                if temp>answer:
                    answer=temp

    print(answer)

max_stars_repo_name: FkrBk/IEEEXtreme

max_issues_repo_name: FkrBk/IEEEXtreme

max_forks_repo_name: FkrBk/IEEEXtreme