text:
#!/usr/bin/python
#
# Project Saturn
# _____________________________________________________________________________
#
#                                                                         _.oo.
# August 2019                                    _.u[[/;:,.         .odMMMMMM'
#                                             .o888UU[[[/;:-.  .o@P^    MMM^
# afterburner.py                             oN88888UU[[[/;::-.        dP^
# Post-hoc data-merging                     dNMMNN888UU[[[/;:--.   .o@P^
#                                          ,MMMMMMN888UU[[/;::-. o@^
#                                          NNMMMNN888UU[[[/~.o@P^
# Markus Ernst                             888888888UU[[[/o@^-..
#                                         oI8888UU[[[/o@P^:--..
#                                      .@^  YUU[[[/o@^;::---..
#                                    oMP     ^/o@P^;:::---..
#                                 .dMMM    .o@^ ^;::---...
#                                dMMMMMMM@^`       `^^^^
#                               YMMMUP^
#                                ^^
# _____________________________________________________________________________
#
#
# Copyright 2019 Markus Ernst
#
# 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 libraries
# ----------------

# standard libraries
# -----
#import tensorflow as tf
import tensorflow.compat.v1 as tf
tf.disable_eager_execution()
import numpy as np
import pickle
import os
import sys
import matplotlib.pyplot as plt
import seaborn as sns
import matplotlib.style as mplstyle
from matplotlib.ticker import MaxNLocator
from matplotlib import collections
from matplotlib.patches import Rectangle
from scipy import stats
from collections import namedtuple
from cycler import cycler

# custom functions
# -----
import utilities.visualizer as visualizer
import utilities.tfevent_handler as tfevent_handler

# commandline arguments
# -----
tf.flags.DEFINE_string('path_to_experiment', '/Users/markus/Research/' +
                           'Code/saturn/experiments/001_noname_experiment/',
                           'path to the experiment, i.e. the data_essences')

FLAGS = tf.flags.FLAGS

class DataEssence(object):
    """docstring for DataEssence."""

    def __init__(self):
        super(DataEssence, self).__init__()

    def write_to_file(self, filename):
        with open(filename, 'wb') as output:  # Overwrites any existing file.
            pickle.dump(self.essence, output, pickle.HIGHEST_PROTOCOL)

    def read_from_file(self, filename):
        with open(filename, 'rb') as input:
            self.essence = pickle.load(input)

    def distill(self, path, evaluation_data, embedding_data=None):
        self.essence = self._read_tfevents(path)

        needless_keys = [l for l in list(self.essence.keys()) if '/' in l]
        for key in needless_keys:
            self.essence.pop(key, None)

        self.essence['evaluation'] = evaluation_data
        if embedding_data:
            self.essence['embedding'] = embedding_data

    def _read_tfevents(self, path):
        em = tfevent_handler.read_multiple_runs(path)
        df = tfevent_handler.convert_em_to_df(em)
        return df

    def plot_essentials(self, savefile):
        # start figure
        fig, axes = plt.subplots(4, 3, sharex='all', figsize=[7, 9])
        # plot images onto figure
        self._plot_traintest_lcurve(axes)
        self._plot_timebased_lcurve(axes)
        self._plot_parameter_lcurve(fig)

        # save figure
        fig.suptitle(savefile.rsplit('/')[-1])
        fig.savefig(savefile)
        pass

    def _plot_traintest_lcurve(self, axes):
        training_acc = \
            self.essence['training']['accuracy_and_error/accuracy'].tolist()
        training_partacc = \
            self.essence['training']['accuracy_and_error/partial_accuracy']\
            .tolist()
        training_loss = \
            self.essence['training']['accuracy_and_error/loss'].tolist()
        training_steps = \
            self.essence['training']['step'].tolist()

        testing_acc = \
            self.essence['testing']['testtime/accuracy'].tolist()
        testing_partacc = \
            self.essence['testing']['testtime/partial_accuracy'].tolist()
        testing_loss = \
            self.essence['testing']['testtime/loss'].tolist()
        testing_steps = \
            self.essence['testing']['step'].tolist()

        axes[0, 0].plot(training_steps, training_acc)
        axes[0, 0].plot(testing_steps, testing_acc)
        axes[0, 0].set_title('accuracy')

        axes[0, 2].plot(training_steps, training_partacc)
        axes[0, 2].plot(testing_steps, testing_partacc)
        axes[0, 2].set_title('partial accuracy')

        axes[0, 1].plot(training_steps, training_loss)
        axes[0, 1].plot(testing_steps, testing_loss)
        axes[0, 1].set_title('loss')

        pass

    def _plot_timebased_lcurve(self, axes):
        sorted_list_of_keys = self.essence['testing'].keys().tolist()
        sorted_list_of_keys.sort()
        timedepth = int(sorted_list_of_keys[-1].rsplit('_')[-1])

        sorted_list_of_keys = \
            [key for key in sorted_list_of_keys if 'testtime_beyond' in key]
        sorted_list_of_keys = \
            [key.rsplit('_', 1)[0] for key in sorted_list_of_keys]
        sorted_list_of_keys = list(set(sorted_list_of_keys))
        sorted_list_of_keys.sort()
        for i in range(len(sorted_list_of_keys)):
            for j in range(timedepth + 1):
                axes[1, i].plot(self.essence['testing']['step'].tolist(),
                                self.essence['testing']['{}_{}'.format(
                                    sorted_list_of_keys[i], j)].tolist(),
                                label=sorted_list_of_keys[i].split('/')[-1])
            axes[1, i].set_title(sorted_list_of_keys[i].split('/')[-1])
        pass

    def _plot_parameter_lcurve(self, fig):
        sorted_list_of_keys = self.essence['training'].keys().tolist()
        sorted_list_of_keys.sort()
        sorted_list_of_keys = \
            [key for key in sorted_list_of_keys if 'weights_and_biases' in key]
        sorted_list_of_keys = \
            [key.rsplit('/', 1)[0] for key in sorted_list_of_keys]
        sorted_list_of_keys = list(set(sorted_list_of_keys))
        sorted_list_of_keys = \
            [key for key in sorted_list_of_keys if '_weights' in key]
        sorted_list_of_keys.sort()
        # for i in range(len(sorted_list_of_keys)):
        for i, ax in enumerate(fig.axes[6:]):
            if i < len(sorted_list_of_keys):
                ax.plot(
                    self.essence['training']['step'].tolist(),
                    self.essence['training']['{}/mean_1'.format(
                        sorted_list_of_keys[i])].tolist(),
                    label=sorted_list_of_keys[i].split('/')[1],
                    color='tab:blue')
                ax.plot(
                    self.essence['training']['step'].tolist(),
                    self.essence['training']['{}/min_1'.format(
                        sorted_list_of_keys[i])].tolist(),
                    label=sorted_list_of_keys[i].split('/')[1],
                    color='tab:blue')
                ax.plot(
                    self.essence['training']['step'].tolist(),
                    self.essence['training']['{}/max_1'.format(
                        sorted_list_of_keys[i])].tolist(),
                    label=sorted_list_of_keys[i].split('/')[1],
                    color='tab:blue')
                ax.fill_between(
                    self.essence['training']['step'],
                    self.essence['training']['{}/mean_1'.format(
                        sorted_list_of_keys[i])] +
                    self.essence['training']['{}/stddev_1'.format(
                        sorted_list_of_keys[i])],
                    self.essence['training']['{}/mean_1'.format(
                        sorted_list_of_keys[i])] -
                    self.essence['training']['{}/stddev_1'.format(
                        sorted_list_of_keys[i])],
                    color='tab:blue',
                    alpha=0.3
                )
                ax.set_title(
                    sorted_list_of_keys[i].split('/')[1].split('_')[1])
        pass

    # TODO: write these visualizer functions
    def visualize_softmax_output(self):
        pass

    def vizualize_tsne_analysis(self):
        pass

class EssenceCollection(object):
    """docstring for EssenceCollection."""

    def __init__(self, remove_files=False):
        super(EssenceCollection, self).__init__()
        path_to_file = os.path.realpath(__file__)
        self.path_to_experiment = path_to_file.rsplit('/', 3)[0]
        self.collection = self.collect_data_essences(
            self.path_to_experiment, remove_files)

    def collect_data_essences(self, path_to_experiment, remove_files):
        # gather and read all files in files/essence/
        collection = {}
        essence = DataEssence()
        path_to_data = path_to_experiment + "/data/"
        for file in os.listdir(path_to_data):
            if file.endswith(".pkl") and file.startswith("conf"):
                config_name = file.split('.')[0].rsplit('i', 1)[0]
                iteration_number = file.split('.')[0].rsplit('i', 1)[-1]
                essence.read_from_file(os.path.join(path_to_data, file))
                if config_name in collection.keys():
                    collection[config_name][iteration_number] = essence.essence
                else:
                    collection[config_name] = {}
                    collection[config_name][iteration_number] = essence.essence
                # delete file
                if remove_files:
                    os.remove(os.path.join(path_to_data, file))
                    print("[INFO] File '{}' deleted.".format(file))
                else:
                    pass
        return collection

    def write_to_file(self, filename):
        with open(filename, 'wb') as output:  # Overwrites any existing file.
            pickle.dump(self.collection, output, pickle.HIGHEST_PROTOCOL)

    def read_from_file(self, filename):
        with open(filename, 'rb') as input:
            self.collection = pickle.load(input)

    def _plot_barplot_comparison(self, filename):

        configurations_to_plot = "all"

        # mplstyle.use('default')

        rects = []
        n_groups = len(self.collection.keys())
        n_iterations = 100  # high number so that true value is smaller
        for configuration in self.collection.keys():
            n_iterations = min(n_iterations,
                               len(self.collection[configuration].keys()))

        n_bars = 1
        data = np.zeros([n_iterations, n_bars, n_groups])
        means = np.zeros([n_bars, n_groups])
        stderror = np.zeros([n_bars, n_groups])
        mcnemar_data = {}
        ev_it = 1

        list_of_configurations = list(self.collection.keys())
        list_of_configurations.sort()

        for j, configuration in enumerate(list_of_configurations):
            list_of_iterations = list(
                self.collection['config{}'.format(j)].keys())
            list_of_iterations.sort()

            for it, iteration in enumerate(list_of_iterations):
                data[it, 0, j] = 1 - \
                    self.collection[configuration][iteration]['testing']['testtime/partial_accuracy'].values[-1]
                mcnemar_data[it, 0, j] = \
                    self.collection[configuration][iteration]['evaluation']['boolean_classification']

        chancelevel = 1 - 1. / 10.

        means = np.mean(data, axis=0)
        stderror = np.std(data, axis=0) / np.sqrt(n_iterations)

        # printout
        # -----
        print((means).round(3))
        print((stderror).round(3))

        fig, ax = plt.subplots(figsize=[8.0 / 6 * n_groups, 6.0])

        index = np.arange(n_groups)
        bar_width = 0.8

        opacity = 1
        error_config = {'ecolor': '0.3'}

        current_palette = sns.dark_palette(
            sns.color_palette('colorblind')[0], 5, reverse=True)
        current_palette = sns.dark_palette(
            sns.color_palette('Set3')[0], 5, reverse=True)
        current_palette = sns.dark_palette(
            'LightGray', 5, reverse=True)
        current_palette = sns.dark_palette(
            'Black', 5, reverse=False)

        for k in range(n_bars):
            rects.append(ax.bar(index + k * bar_width, means[k], bar_width,
                                alpha=opacity,
                                yerr=stderror[k], error_kw=error_config,
                                label='{}'.format('configuration'),
                                linewidth=0, color=current_palette[k]))
            # , edgecolor=current_palette[(k+1)],
            # hatch=hatches[k]))

        ax.set_xlabel('Network configuration')
        ax.set_ylabel('Error rate')
        ax.set_title('experiment: {}'.format(
            self.path_to_experiment.rsplit('/')[-1]), fontsize=10)
        ax.set_ylim([0, ax.get_ylim()[1] * 1])  # *1.25])
        ax.set_xticks(index + bar_width * (n_bars / 2. - 0.5))
        # ax.set_xticklabels([arch[:-1] for arch in ARCHARRAY])
        ax.legend(frameon=True, facecolor='white',
                  edgecolor='white', framealpha=1.0)
        ax.grid(axis='y', zorder=0, alpha=0.5)
        ax.set_axisbelow(True)
        ax.spines['top'].set_visible(False)
        ax.spines['right'].set_visible(False)
        ax.tick_params(
            axis='x',          # changes apply to the x-axis
            which='both',      # both major and minor ticks are affected
            bottom=False,      # ticks along the bottom edge are off
            top=False,         # ticks along the top edge are off
            labelbottom=True)  # labels along the bottom edge are off
        ax.spines['left'].set_visible(False)

        ax.axhline(y=chancelevel, xmin=0, xmax=6, color='black')
        ax.text(0.2, chancelevel + 0.005, r'chancelevel')

        # Shrink current axis by 20%
        box = ax.get_position()
        ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
        # Put a legend to the right of the current axis
        ax.legend(frameon=True, facecolor='white', edgecolor='white',
                  framealpha=1.0, bbox_to_anchor=(1, .9), loc='center left')

        # mcnemar's table
        for i in range(n_bars):
            qstar = 0.05
            pval = np.ones([n_groups * n_bars, n_groups * n_bars])
            chi_table = np.ones([n_groups * n_bars, n_groups * n_bars])
            significance_table = np.zeros(
                [n_groups * n_bars, n_groups * n_bars])
            vertexlist = []
            for k in range(n_groups * n_bars):
                for j in range(n_groups * n_bars):
                    if k != j and k > j:
                        mcnemar_table = 2 * \
                            mcnemar_data[(ev_it - 1), i, k] - \
                            mcnemar_data[(ev_it - 1), i, j]
                        # i.e. nhot encoding
                        if (len(mcnemar_data[(ev_it - 1), i, k].shape) > 1):
                            ak = np.sum(
                                np.array(
                                    [e for e in mcnemar_table == 1.]), axis=1)
                            bk = np.sum(
                                np.array(
                                    [e for e in mcnemar_table == 2.]), axis=1)
                            ck = np.sum(
                                np.array(
                                    [e for e in mcnemar_table == -1.]), axis=1)
                            dk = np.sum(
                                np.array(
                                    [e for e in mcnemar_table == 0]), axis=1)
                            chi_sq_denom = np.sum(
                                ((bk - ck) / mcnemar_table.shape[-1]))**2
                            chi_sq_count = np.sum(
                                ((bk - ck) / mcnemar_table.shape[-1])**2)
                            chi_sq = chi_sq_denom / chi_sq_count
                        else:
                            a = mcnemar_table[mcnemar_table == 1.].shape[0]
                            b = mcnemar_table[mcnemar_table == 2.].shape[0]
                            c = mcnemar_table[mcnemar_table == -1.].shape[0]
                            d = mcnemar_table[mcnemar_table == 0].shape[0]
                            chi_sq = ((b - c)**2) / (b + c)
                        # stats.chi2.pdf(chi_sq , 1) #
                        pval[k, j] = 1 - stats.chi2.cdf(chi_sq, 1)
                        chi_table[k, j] = chi_sq

            print(np.round(pval, 4))
            print(np.round(chi_table, 4))
            sorted_pvals = np.sort(pval[pval < 1])
            bjq = np.arange(1, len(sorted_pvals) + 1) / \
                len(sorted_pvals) * qstar

            for k in range(n_groups * n_bars):
                for j in range(n_groups * n_bars):
                    if k != j and k > j:
                        if pval[k, j] in sorted_pvals[sorted_pvals - bjq <= 0]:
                            significance_table[k, j] = 1

                            if j < 3 and k > 2 and (means[i][j] < means[i][k]):
                                vertexlist.append(
                                    [[j + 0.25, k - 0.5],
                                     [j + 0.5, k - 0.25],
                                     [j + 0.5, k - 0.5]])
                        else:
                            significance_table[k, j] = 0

            ax = plt.axes([0.8, 0.3, .2, .2])
            ax.matshow(significance_table, cmap='Greys')
            # ax.set_xticklabels(, fontsize=8)#fontsize=65)
            ax.tick_params(labelbottom='on', labeltop='off',
                           top='off', right='off')
            # ax.set_yticklabels(, fontsize=8)#fontsize=65)
            ax.set_xlim([-0.5, n_groups * n_bars - 1 - 0.5])
            ax.set_ylim([n_groups * n_bars - 0.5, -0.5 + 1])
            ax.spines['top'].set_visible(False)
            ax.spines['right'].set_visible(False)

            ax.axhline(y=2.5, xmin=-0.5, xmax=2.5, color='white', linewidth=1)
            ax.axvline(x=2.5, ymin=-0.5, ymax=2.5, color='white', linewidth=1)
            for vertices in vertexlist:
                pc = collections.PolyCollection(
                    (vertices,), color='white', edgecolor="none")
                ax.add_collection(pc)

            ax.add_patch(Rectangle((-2.5, 6.9), 1, 1, fill='black',
                                   color='black', alpha=1, clip_on=False))
            ax.text(-2.5, 9.4, 'Significant difference \n(two-sided McNemar \
                test, \nexpected FDR=0.05)',
                    fontsize=8, horizontalalignment='left',
                    verticalalignment='center')
            # plt.show()
            plt.savefig(filename)

            fig, ax = plt.subplots(figsize=(12, 12))
            ax.matshow(significance_table, cmap='Greys')
            # ax.set_xticklabels(
            #     [''] + [arch[:-1] for arch in ARCHARRAY], fontsize=65)
            ax.tick_params(labelbottom='on', labeltop='off',
                           top='off', right='off')
            # ax.set_yticklabels(
            #     [arch[:-1] for arch in ARCHARRAY], fontsize=65)
            ax.set_xlim([-0.5, n_groups * n_bars - 1 - 0.5])
            ax.set_ylim([n_groups * n_bars - 0.5, -0.5 + 1])
            ax.spines['top'].set_visible(False)
            ax.spines['right'].set_visible(False)
            ax.axhline(y=2.5, xmin=-0.5, xmax=2.5, color='white')
            ax.axvline(x=2.5, ymin=-0.5, ymax=2.5, color='white')
            for vertices in vertexlist:
                pc = collections.PolyCollection(
                    (vertices,), color='white', edgecolor="none")
                ax.add_collection(pc)
            # ax.set_ylabel('Network Architecture')
        plt.tight_layout()
        plt.savefig(filename.rsplit('.', 1)[0] + '_sigtable.pdf')

# TODO: Think about how to implement this, low priority
def freeze_model(path_to_model, path_to_checkpoint):
    pass

def unfreeze_model(path_to_model, path_to_checkpoint):
    pass

evaluation_filenames = []
evaluate_ckpt = False

if __name__ == '__main__':
    print('[INFO] afterburner running, collecting data')
    ess_coll = EssenceCollection(remove_files=True)

    ess_coll._plot_barplot_comparison(
        ess_coll.path_to_experiment +
        '/visualization/{}_comparision.pdf'.format(
            ess_coll.path_to_experiment.rsplit('/')[-1]))
    ess_coll.write_to_file(
        ess_coll.path_to_experiment +
        '/data/{}.pkl'.format(
            ess_coll.path_to_experiment.rsplit('/')[-1]))

# # save output of boolean comparison
# boolfile = open(WRITER_DIRECTORY + 'checkpoints/' +
#                 'evaluation/' + "bool_classification.txt", "a")
#
# if CONFIG['label_type'] == "onehot":
#     for i in list(bc):
#         boolfile.write(str(int(i)) + '\n')
# else:
#     for i in range(len(bc[0])):
#         for el in bc[0][i]:
#             boolfile.write(
#                 str(int(el in set(bc[1][i]))) + '\t')
#         boolfile.write('\n')
# boolfile.close()

# np.savez_compressed(WRITER_DIRECTORY + 'checkpoints/' + 'evaluation/' +
#                     'softmax_output.npz',
#                     np.array(list_of_output_values))

# _____________________________________________________________________________
# Description:
#
# This program is supposed to be invoked after the experiment is over,
# it includes functions to create a numpy pickle file with all the important
# data for plotting.
# If invoked as a main program
# the afterburner tries to combine the experiments into useful comparison plots
# _____________________________________________________________________________


text:
"""This module implements the CircuitGate class."""
from __future__ import annotations

from typing import TYPE_CHECKING

import numpy as np
import numpy.typing as npt

from bqskit.qis.unitary.unitary import RealVector
from bqskit.qis.unitary.unitarymatrix import UnitaryMatrix

if TYPE_CHECKING:
    from bqskit.ir.circuit import Circuit

from bqskit.ir.gate import Gate

class CircuitGate(Gate):
    """
    The CircuitGate class.

    A CircuitGate is a immutable circuit represented as a gate.
    """

    def __init__(self, circuit: Circuit, move: bool = False) -> None:
        """
        CircuitGate Constructor.

        Args:
            circuit (Circuit): The circuit to copy into gate format.

            move (bool): If true, the constructor will not copy the circuit.
                This should only be used when you are sure `circuit` will no
                longer be used on caller side. If unsure use the default.
                (Default: False)
        """

        self._circuit = circuit if move else circuit.copy()
        self._num_qudits = self._circuit.num_qudits
        self._radixes = self._circuit.radixes
        self._num_params = self._circuit.num_params
        self._name = 'CircuitGate(%s)' % str(self._circuit)

    def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
        """Return the unitary for this gate, see :class:`Unitary` for more."""
        return self._circuit.get_unitary(params)

    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
        """
        Return the gradient for this gate.

        See :class:`DifferentiableUnitary` for more info.
        """
        return self._circuit.get_grad(params)

    def get_unitary_and_grad(
        self,
        params: RealVector = [],
    ) -> tuple[UnitaryMatrix, npt.NDArray[np.complex128]]:
        """
        Return the unitary and gradient for this gate.

        See :class:`DifferentiableUnitary` for more info.
        """
        return self._circuit.get_unitary_and_grad(params)

    def is_differentiable(self) -> bool:
        """Return true if the circuit is differentiable."""
        return self._circuit.is_differentiable()

    def __hash__(self) -> int:
        hashes: list[int] = [hash(self.name)]
        for op in self._circuit:
            hashes.append(hash(op))

            # Don't let the hash list grow too large.
            if len(hashes) >= 100:
                hashes = [hash(tuple(hashes))]

        hash_val = hash(tuple(hashes)) if len(hashes) > 1 else hashes[0]
        return hash_val

    def __eq__(self, other: object) -> bool:
        if not isinstance(other, CircuitGate):
            return NotImplemented

        if self._circuit.num_qudits != other._circuit.num_qudits:
            return False

        if self._circuit.radixes != other._circuit.radixes:
            return False

        return all(
            op1.gate == op2.gate and op1.location == op2.location
            for op1, op2 in zip(self._circuit, other._circuit)
        )

ext: py

reason: import numpy


text:
import cv2
import numpy as np
import matplotlib.pyplot as plt
import iris.utils
import iris.utils.contours
from iris.utils import ImageDataReader, ImageViewer
import settings

reader = ImageDataReader()
viewer = ImageViewer()

data = reader.read(settings.DATA_FOLDER_PATH, gray=True)

# choose a sample image
img = data["14"]["left"][4]
viewer.addImage("1-Left", img)

i = iris.utils.extractPupil(img)
viewer.addImage("EyeLid", i)

num_labels, labels_im = cv2.connectedComponents(i)

def imshow_components(labels):
    # Map component labels to hue val
    label_hue = np.uint8(179*labels/np.max(labels))
    blank_ch = 255*np.ones_like(label_hue)
    labeled_img = cv2.merge([label_hue, blank_ch, blank_ch])

    # cvt to BGR for display
    labeled_img = cv2.cvtColor(labeled_img, cv2.COLOR_HSV2BGR)

    # set bg label to black
    labeled_img[label_hue==0] = 0

    viewer.addImage('labeled.png', labeled_img)

imshow_components(labels_im)

viewer.show()
plt.show()


text:
import sys
import os
sys.path.append('/Users/adam2392/Documents/tvb/')
# sys.path.append('/Users/adam2392/Documents/tvb/_tvbdata/')
# sys.path.append('/Users/adam2392/Documents/tvb/_tvblibrary/')
# from tvb.simulator.lab import *
import numpy as np
import pandas as pd
import scipy
from sklearn import cluster
from sklearn.preprocessing import MinMaxScaler
import tvbsim

'''
Module for preprocessing and organizing data computed by FFT on iEEG data
into a huge compressed data structure

Data will be stored as npz compressed file

matrix will be H x W x F x T (height x width x frequency band x time window),
T x H x W x F
where height and width define a grid where power is projected.

It will depend on the number of parcellated regions, so for 84, it will be a
12x7 image x4 frequency bands (alpha, beta, gamma, high freq)
'''

# get factors of the number of regions

def get_factors(x):
   # This function takes a number and prints the factors

    factors = []
    for i in range(1, x + 1):
        if x % i == 0:
            factors.append(i)
    return factors

class PreProcess(object):
    def __init__(self, datafiles=[], freqbands=None):
        if not freqbands:
            # establish frequency bands
            freqbands = {
                'dalpha': [0, 15],
                'beta': [15, 30],
                'gamma': [30, 90],
                'high': [90, 200],
                'hfo': [200, 500]
            }
        self.datafiles = datafiles
        self.freqbands = freqbands

    def _computefreqindices(self, freqs, freqbands):
        freqbandindices = {}
        for band in freqbands:
            lowerband = freqbands[band][0]
            upperband = freqbands[band][1]

            # get indices where the freq bands are put in
            freqbandindices[band] = np.where(
                (freqs >= lowerband) & (freqs < upperband))
            freqbandindices[band] = [freqbandindices[band]
                                     [0][0], freqbandindices[band][0][-1]]
        return freqbandindices

    def binFrequencyValues(self, power, freqs):
        # List of physiological frequency ranges
        freqbands = self.freqbands

        # Create an empty array
        power_binned = np.zeros(shape=(power.shape[0],
                                       len(freqbands),
                                       power.shape[2]))

        # compute the freq indices for each band
        freqbandindices = self._computefreqindices(freqs, freqbands)

        # compress data using frequency bands
        for idx, band in enumerate(freqbandindices):
            indices = freqbandindices[band]
            # average between these two indices
            power_binned[:, idx, :] = np.mean(
                np.abs(power[:, indices[0]:indices[1]+1, :]), axis=1)  # [np.newaxis,:,:]

        return power_binned

    def compresspowermat(self, datapath):
        print(os.path.join(datapath))
        freqbands = self.freqbands

        # load data from this datafile and run
        powerbands = {}
        data = np.load(os.path.join(datapath), encoding='bytes')

        # extract data from the numpy file
        power = data['power']
        freqs = data['freqs']
        timepoints = data['timepoints']
        metadata = data['metadata'].item()
        metadata[b'freqbands'] = freqbands

        # compute the freq indices for each band
        freqbandindices = self._computefreqindices(freqs, freqbands)

        # compress data using frequency bands
        for idx, band in enumerate(freqbandindices):
            indices = freqbandindices[band]
            # average between these two indices
            powerbands[band] = np.mean(
                power[:, indices[0]:indices[1]+1, :], axis=1)  # [np.newaxis,:,:]

        return powerbands, timepoints

    def getseiztimes(self, datafile):
        '''
        Create an accompanying list of datafiles to save, so that it is a metadata
        list of all the data that is compressed into the final data structure
        '''
        data = np.load(datafile, encoding='bytes')
        metadata = data['metadata'].item()
        metadata = {k.decode("utf-8"): (v.decode("utf-8")
                                        if isinstance(v, bytes) else v) for k, v in metadata.items()}

        onsettimes = metadata['onsettimes']
        offsettimes = metadata['offsettimes']

        if len(onsettimes) > 0:
            if onsettimes[0] > offsettimes[0]:
                onsettimes = metadata['offsettimes']
                offsettimes = metadata['onsettimes']
        return onsettimes, offsettimes

    def loadmetadata(self, datafile):
        data = np.load(datafile, encoding='bytes')
        # extract data from the numpy file
        metadata = data['metadata'].item()
        metadata = {k.decode("utf-8"): (v.decode("utf-8")
                                        if isinstance(v, bytes) else v) for k, v in metadata.items()}

        return metadata

    def getimdims(self, metadata):
        regions = metadata['regions']

        # reshape the regions of 84 into a parcellated rectangular "image"
        # height = np.ce
        factors = get_factors(len(regions))
        height = factors[int(len(factors)/2)]
        width = int(len(regions) / height)

        return height, width

    def projectpower_gain(self, powerbands, metadata, verts, areas, regmap):
        seeg_xyz = metadata['seeg_xyz']
        seeg_labels = metadata['seeg_contacts']

        # extract the seeg_xyz coords and the region centers
        region_centers = metadata['region_centers']
        regions = metadata['regions']

        height, width = self.getimdims(metadata)

        # check seeg_xyz moved correctly - In early simulation data results, was not correct
        buff = seeg_xyz - region_centers[:, np.newaxis]
        buff = np.sqrt((buff**2).sum(axis=-1))
        test = np.where(buff == 0)
        indice = test[1]

        modgain = tvbsim.util.gain_matrix_inv_square(verts, areas,
                                                     regmap, len(regions), seeg_xyz)
        modgain = modgain.T

        # map seeg activity -> epileptor source and create data structure
        for idx, band in enumerate(powerbands):
            mapped_power_band = np.tensordot(
                modgain, powerbands[band], axes=([1], [0]))

            if idx == 0:
                mapped_power_bands = mapped_power_band.reshape(height, width, mapped_power_band.shape[1],
                                                               order='C')[np.newaxis, :, :, :]
            else:
                mapped_power_bands = np.append(mapped_power_bands,
                                               mapped_power_band.reshape(height,
                                                                         width,
                                                                         mapped_power_band.shape[1],
                                                                         order='C')[np.newaxis, :, :, :], axis=0)

        # new condensed data structure is H x W x F x T, to concatenate more, add to T dimension
        mapped_power_bands = mapped_power_bands.swapaxes(0, 3).swapaxes(1, 2)

        return mapped_power_bands

    def projectpower_knn(self, powerbands, metadata):
        # extract the seeg_xyz coords and the region centers
        region_centers = metadata['region_centers']
        regions = metadata['regions']
        seeg_xyz = metadata['seeg_xyz']
        seeg_labels = metadata['seeg_contacts']

        numfreqbands = len(self.freqbands)

        height, width = self.getimdims(metadata)

        # map seeg_xyz to 3 closest region_centers
        tree = scipy.spatial.KDTree(region_centers)
        seeg_near_indices = []

        seeg_counter = np.zeros(len(regions))

        for ichan in range(0, len(seeg_labels)):
            near_regions = tree.query(seeg_xyz[ichan, :].squeeze(), k=3)
            near_indices = near_regions[1]

            # go through each frequency band and map activity onto those near indices
            for idx, band in enumerate(powerbands):
                # initialize the final result data structure
                if ichan == 0 and idx == 0:
                    mapped_power_bands = np.zeros((len(regions),
                                                   numfreqbands,
                                                   powerbands[band].shape[1]), dtype='float64')

                chanpower = powerbands[band][ichan, :]
                mapped_power_bands[near_indices, idx,
                                   :] += chanpower.astype('float64')

            seeg_counter[near_indices] += 1
            # seeg_near_indices.append(near_indices)

        # get the average based on how many contributions of the seeg power was to this region
        mapped_power_bands = np.divide(mapped_power_bands, seeg_counter[:, np.newaxis, np.newaxis],
                                       out=np.zeros_like(mapped_power_bands), where=seeg_counter[:, np.newaxis, np.newaxis] != 0)

        # reshape for the correct output
        mapped_power_bands = mapped_power_bands.reshape(height,
                                                        width,
                                                        numfreqbands,
                                                        powerbands[band].shape[1],
                                                        order='C')
        mapped_power_bands = mapped_power_bands.swapaxes(0, 3).swapaxes(2, 3)

        return mapped_power_bands

    def projectpower_invsquare(self, powerbands, metadata):
        # extract the seeg_xyz coords and the region centers
        region_centers = metadata['region_centers']
        regions = metadata['regions']
        seeg_xyz = metadata['seeg_xyz']
        seeg_labels = metadata['seeg_contacts']

        height, width = self.getimdims(metadata)

        # map seeg_xyz to the rest of the regions from a factor of
        # computes distance along each axis
        dr = region_centers - seeg_xyz[:, np.newaxis]
        ndr = np.sqrt((dr**2).sum(axis=-1))  # computes euclidean distance
        Vr = 1/(ndr**2)  # fall off as a function of r^2

        inf_indices = np.where(np.isinf(Vr))
        small_indices = np.where(ndr <= 1)

        # can either set to 1, or the max that there currently is + some small epsilon
        # the problem with setting to 1 is that the signal might drown out everything else
        Vr[small_indices] = np.nanmax(
            np.ma.masked_invalid(Vr[:])) + np.nanmin(Vr[:])

        # normalize Vr with minmax
        scaler = MinMaxScaler(feature_range=(0, 1))
        Vr = scaler.fit_transform(Vr).T

        # map seeg activity -> epileptor source and create data structure
        for idx, band in enumerate(powerbands):
            mapped_power_band = np.tensordot(
                Vr, powerbands[band], axes=([1], [0]))

            # store the formatted power bands
            if idx == 0:
                mapped_power_bands = mapped_power_band.reshape(height,
                                                               width,
                                                               mapped_power_band.shape[1],
                                                               order='C')[np.newaxis, :, :, :]
            else:
                mapped_power_bands = np.append(mapped_power_bands,
                                               mapped_power_band.reshape(height,
                                                                         width,
                                                                         mapped_power_band.shape[1],
                                                                         order='C')[np.newaxis, :, :, :], axis=0)

        # new condensed data structure is H x W x F x T, to concatenate more, add to T dimension
        mapped_power_bands = mapped_power_bands.swapaxes(0, 3).swapaxes(1, 2)

        return mapped_power_bands

reason: import numpy,import scipy


text:
'''Example to illustrate Quantile Regression

Author: Josef Perktold

polynomial regression with systematic deviations above

'''

import numpy as np
from statsmodels.compat.python import zip
from scipy import stats
import statsmodels.api as sm

from statsmodels.regression.quantile_regression import QuantReg

sige = 0.1
nobs, k_vars = 500, 3
x = np.random.uniform(-1, 1, size=nobs)
x.sort()
exog = np.vander(x, k_vars+1)[:,::-1]
mix = 0.1 * stats.norm.pdf(x[:,None], loc=np.linspace(-0.5, 0.75, 4), scale=0.01).sum(1)
y = exog.sum(1) + mix + sige * (np.random.randn(nobs)/2 + 1)**3

p = 0.5
res_qr = QuantReg(y, exog).fit(p)
res_qr2 = QuantReg(y, exog).fit(0.1)
res_qr3 = QuantReg(y, exog).fit(0.75)
res_ols = sm.OLS(y, exog).fit()

params = [res_ols.params, res_qr2.params, res_qr.params, res_qr3.params]
labels = ['ols', 'qr 0.1', 'qr 0.5', 'qr 0.75']

import matplotlib.pyplot as plt

plt.figure()
plt.plot(x, y, '.', alpha=0.5)
for lab, beta in zip(['ols', 'qr 0.1', 'qr 0.5', 'qr 0.75'], params):
    print('%-8s' % lab, np.round(beta, 4))
    fitted = np.dot(exog, beta)
    lw = 2
    plt.plot(x, fitted, lw=lw, label=lab)
plt.legend()
plt.title('Quantile Regression')

plt.show()

max_stars_repo_name: haribharadwaj/statsmodels

max_issues_repo_name: haribharadwaj/statsmodels

max_forks_repo_name: haribharadwaj/statsmodels

reason: import numpy,from scipy,import statsmodels,from statsmodels


text:
from ledfx.effects.audio import AudioReactiveEffect, FREQUENCY_RANGES_SIMPLE
from ledfx.effects.gradient import GradientEffect
import voluptuous as vol
import numpy as np

class ScrollAudioEffect(AudioReactiveEffect):

    NAME = "Scroll"

    CONFIG_SCHEMA = vol.Schema({
        vol.Optional('blur', description='Amount to blur the effect', default = 3.0): vol.Coerce(float),
        vol.Optional('mirror', description='Mirror the effect', default = True): bool,
        vol.Optional('speed', description='Speed of the effect', default = 4):  vol.All(vol.Coerce(int), vol.Range(min=1)),
        vol.Optional('decay', description='Decay rate of the scroll', default = 0.97):  vol.All(vol.Coerce(float), vol.Range(min=0.2, max=1.0)),
    })

    def config_updated(self, config):

        # TODO: Determine how buffers based on the pixels should be
        # allocated. Technically there is no guarantee that the effect
        # is bound to a device while the config gets updated. Might need
        # to move to a model where effects are created for a device and
        # must be destroyed and recreated to be moved to another device.
        self.output = None

    def audio_data_updated(self, data):

        if self.output is None:
            self.output = self.pixels

        # Divide the melbank into lows, mids and highs
        lows_max = np.clip(np.max(data.melbank_lows() ** 2), 0, 1)
        mids_max = np.clip(np.max(data.melbank_mids() ** 2), 0, 1)
        highs_max = np.clip(np.max(data.melbank_highs() ** 2), 0, 1)

        if lows_max < 0.2:
            lows_max = 0
        if mids_max < 0.2:
            mids_max = 0
        if highs_max < 0.2:
            highs_max = 0

        # Compute the value for each range based on the max
        lows_val = (np.array((255,0,0)) * lows_max)
        mids_val = (np.array((0,255,0)) * mids_max)
        high_val = (np.array((0,0,255)) * highs_max)

        # Roll the effect and apply the decay
        speed = self.config['speed']
        self.output[speed:,:] = self.output[:-speed,:]
        self.output = (self.output * self.config['decay'])

        # Add in the new color from the signal maxes
        self.output[:speed, 0] = lows_val[0] + mids_val[0] + high_val[0]
        self.output[:speed, 1] = lows_val[1] + mids_val[1] + high_val[1]
        self.output[:speed, 2] = lows_val[2] + mids_val[2] + high_val[2]

        # Set the pixels
        self.pixels = self.output

max_issues_repo_name: simon-wh/LedFx

max_forks_repo_name: simon-wh/LedFx


text:
# -*- coding: utf-8 -*-
"""This module defines MSA analysis functions."""

__author__ = 'Anindita Dutta, Ahmet Bakan, Wenzhi Mao'

from numbers import Integral
import os

from numpy import dtype, zeros, empty, ones, where, ceil, shape, eye
from numpy import indices, tril_indices, array, ndarray, isscalar, unique

from prody import LOGGER
from prody.utilities import which, MATCH_SCORE, MISMATCH_SCORE
from prody.utilities import GAP_PENALTY, GAP_EXT_PENALTY, ALIGNMENT_METHOD

from prody.sequence.msa import MSA, refineMSA
from prody.sequence.msafile import parseMSA, writeMSA
from prody.sequence.sequence import Sequence
from prody.atomic import Atomic
from prody.measure import calcDistance

from Bio import pairwise2
import sys

__all__ = ['calcShannonEntropy', 'buildMutinfoMatrix', 'calcMSAOccupancy',
           'applyMutinfoCorr', 'applyMutinfoNorm', 'calcRankorder', 'filterRankedPairs',
           'buildSeqidMatrix', 'uniqueSequences', 'buildOMESMatrix',
           'buildSCAMatrix', 'buildDirectInfoMatrix', 'calcMeff',
           'buildPCMatrix', 'buildMSA', 'showAlignment', 'alignTwoSequencesWithBiopython',
           'alignSequenceToMSA', 'calcPercentIdentities', 'alignSequencesByChain',
           'trimAtomsUsingMSA']

doc_turbo = """

    By default, *turbo* mode, which uses memory as large as the MSA array
    itself but runs four to five times faster, will be used.  If memory
    allocation fails, the implementation will fall back to slower and
    memory efficient mode."""

def calcPercentIdentities(msa):
    percent_ids = []
    aas = ['A','C','D','E','F','G','H','I','J','K','L', \
           'M','N','P','Q','R','S','T','V','W','Y','-']

    for i in range(len(msa)):
        col_list = list(msa.getArray()[:,i])
        max_count = 0
        for aa in aas:
            if col_list.count(aa) > max_count:
                max_count = col_list.count(aa)
        percent_ids.append(float(max_count)/float(len(col_list))*100)

    return percent_ids

def getMSA(msa):
    """Returns MSA character array."""

    try:
        msa = msa._getArray()
    except AttributeError:
        pass

    try:
        dtype_, ndim, shape = msa.dtype, msa.ndim, msa.shape
    except AttributeError:
        raise TypeError('msa must be an MSA instance or a 2D character array')

    if dtype_ != dtype('|S1') or ndim != 2:
        raise TypeError('msa must be an MSA instance or a 2D character array')

    return msa

def calcShannonEntropy(msa, ambiguity=True, omitgaps=True, **kwargs):
    """Returns Shannon entropy array calculated for *msa*, which may be
    an :class:`.MSA` instance or a 2D Numpy character array.  Implementation
    is case insensitive and handles ambiguous amino acids as follows:

      * **B** (Asx) count is allocated to *D* (Asp) and *N* (Asn)
      * **Z** (Glx) count is allocated to *E* (Glu) and *Q* (Gln)
      * **J** (Xle) count is allocated to *I* (Ile) and *L* (Leu)
      * **X** (Xaa) count is allocated to the twenty standard amino acids

    Selenocysteine (**U**, Sec) and pyrrolysine (**O**, Pyl) are considered
    as distinct amino acids.  When *ambiguity* is set **False**, all alphabet
    characters as considered as distinct types.

    All non-alphabet characters are considered as gaps, and they are handled
    in two ways:

      * non-existent, the probability of observing amino acids in a given
        column is adjusted, by default
      * as a distinct character with its own probability, when *omitgaps* is
        **False**"""

    msa = getMSA(msa)
    length = msa.shape[1]
    entropy = empty(length, float)
    from .msatools import msaentropy
    return msaentropy(msa, entropy,
                      ambiguity=bool(ambiguity), omitgaps=bool(omitgaps))

def buildMutinfoMatrix(msa, ambiguity=True, turbo=True, **kwargs):
    """Returns mutual information matrix calculated for *msa*, which may be an
    :class:`.MSA` instance or a 2D Numpy character array.  Implementation
    is case insensitive and handles ambiguous amino acids as follows:

      * **B** (Asx) count is allocated to *D* (Asp) and *N* (Asn)
      * **Z** (Glx) count is allocated to *E* (Glu) and *Q* (Gln)
      * **J** (Xle) count is allocated to *I* (Ile) and *L* (Leu)
      * **X** (Xaa) count is allocated to the twenty standard amino acids
      * Joint probability of observing a pair of ambiguous amino acids is
        allocated to all potential combinations, e.g. probability of **XX**
        is allocated to 400 combinations of standard amino acids, similarly
        probability of **XB** is allocated to 40 combinations of *D* and *N*
        with the standard amino acids.

    Selenocysteine (**U**, Sec) and pyrrolysine (**O**, Pyl) are considered
    as distinct amino acids.  When *ambiguity* is set **False**, all alphabet
    characters as considered as distinct types.  All non-alphabet characters
    are considered as gaps.

    Mutual information matrix can be normalized or corrected using
    :func:`applyMINormalization` and :func:`applyMICorrection` methods,
    respectively.  Normalization by joint entropy can performed using this
    function with *norm* option set **True**."""

    msa = getMSA(msa)

    from .msatools import msamutinfo
    LOGGER.timeit('_mutinfo')
    length = msa.shape[1]
    mutinfo = empty((length, length), float)
    mutinfo = msamutinfo(msa, mutinfo,
                         ambiguity=bool(ambiguity), turbo=bool(turbo),
                         norm=bool(kwargs.get('norm', False)),
                         debug=bool(kwargs.get('debug', False)))
    LOGGER.report('Mutual information matrix was calculated in %.2fs.',
                  '_mutinfo')

    return mutinfo

buildMutinfoMatrix.__doc__ += doc_turbo

def calcMSAOccupancy(msa, occ='res', count=False):
    """Returns occupancy array calculated for residue positions (default,
    ``'res'`` or ``'col'`` for *occ*) or sequences (``'seq'`` or ``'row'``
    for *occ*) of *msa*, which may be an :class:`.MSA` instance or a 2D
    NumPy character array.  By default, occupancy [0-1] will be calculated.
    If *count* is **True**, count of non-gap characters will be returned.
    Implementation is case insensitive."""

    from .msatools import msaocc

    msa = getMSA(msa)

    try:
        dim = occ.startswith('res') or occ.startswith('col')
    except AttributeError:
        raise TypeError('occ must be a string')
    occ = zeros(msa.shape[int(dim)], float)
    return msaocc(msa, occ, dim, count=bool(count))

def applyMutinfoNorm(mutinfo, entropy, norm='sument'):
    """Apply one of the normalizations discussed in [MLC05]_ to *mutinfo*
    matrix.  *norm* can be one of the following:

      * ``'sument'``: :math:`H(X) + H(Y)`, sum of entropy of columns
      * ``'minent'``: :math:`min\{H(X), H(Y)\}`, minimum entropy
      * ``'maxent'``: :math:`max\{H(X), H(Y)\}`, maximum entropy
      * ``'mincon'``: :math:`min\{H(X|Y), H(Y|X)\}`, minimum conditional
         entropy
      * ``'maxcon'``: :math:`max\{H(X|Y), H(Y|X)\}`, maximum conditional
         entropy

    where :math:`H(X)` is the entropy of a column, and
    :math:`H(X|Y) = H(X) - MI(X, Y)`.  Normalization with joint entropy, i.e.
    :math:`H(X, Y)`, can be done using :func:`.buildMutinfoMatrix` *norm*
    argument.

    .. [MLC05] Martin LC, Gloor GB, Dunn SD, Wahl LM. Using information theory
       to search for co-evolving residues in proteins. *Bioinformatics*
       **2005** 21(22):4116-4124."""

    try:
        ndim, shape = mutinfo.ndim, mutinfo.shape
    except AttributeError:
        raise TypeError('mutinfo must be a 2D square array')

    if ndim != 2 or shape[0] != shape[1]:
        raise ValueError('mutinfo must be a 2D square array')

    try:
        ndim, shapent = entropy.ndim, entropy.shape
    except AttributeError:
        raise TypeError('entropy must be a numpy array')

    if ndim != 1:
        raise ValueError('entropy must be a 1D array')

    if shapent[0] != shape[0]:
        raise ValueError('shape of mutinfo and entropy does not match')

    try:
        sw = norm.startswith
    except AttributeError:
        raise TypeError('norm must be a string')

    if sw('sument'):
        norm = lambda i_val, j_val, val: i_val + j_val

    elif sw('minent'):
        norm = lambda i_val, j_val, val: min(i_val, j_val)

    elif sw('maxent'):
        norm = lambda i_val, j_val, val: max(i_val, j_val)

    elif sw('mincon'):
        norm = lambda i_val, j_val, val: min(i_val - val, j_val - val)

    elif sw('maxcon'):
        norm = lambda i_val, j_val, val: max(i_val - val, j_val - val)

    elif sw('joint'):
        raise ValueError('for joint entropy normalization, use '
                         'buildMutinfoMatrix function')
    else:
        raise ValueError('norm={0} is not a valid normalization type'
                         .format(norm))

    mi = mutinfo.copy()
    for i, i_val in enumerate(entropy):
        for j, j_val in enumerate(entropy):
            val = mi[i, j]
            div = norm(i_val, j_val, val)
            if div == 0:
                mi[i, j] = 0
            else:
                mi[i, j] /= div

    return mi

def applyMutinfoCorr(mutinfo, corr='prod'):
    """Returns a copy of *mutinfo* array after average product correction
    (default) or average sum correction is applied.  See [DSD08]_ for details.

    .. [DSD08] Dunn SD, Wahl LM, Gloor GB. Mutual information without the
       influence of phylogeny or entropy dramatically improves residue
       contact prediction. *Bioinformatics* **2008** 24(3):333-340."""

    try:
        ndim, shape = mutinfo.ndim, mutinfo.shape
    except AttributeError:
        raise TypeError('mutinfo must be a 2D square array')

    if ndim != 2 or shape[0] != shape[1]:
        raise ValueError('mutinfo must be a 2D square array')

    try:
        sw = corr.startswith
    except AttributeError:
        raise TypeError('correction must be a string')

    avg_mipos = mutinfo.sum(1) / (shape[0] - 1)
    avg_mi = avg_mipos.mean()

    mi = mutinfo.copy()
    if sw('prod') or sw('apc'):
        for i, i_avg in enumerate(avg_mipos):
            for j, j_avg in enumerate(avg_mipos):
                mi[i, j] -= (i_avg * j_avg) / avg_mi
    elif sw('sum') or sw('asc'):
        for i, i_avg in enumerate(avg_mipos):
            for j, j_avg in enumerate(avg_mipos):
                mi[i, j] -= i_avg + j_avg - avg_mi
    else:
        raise ValueError('correction must be prod or sum, not ' + corr)

    return mi

def filterRankedPairs(pdb, indices, msa_indices, rank_row, rank_col, zscore_sort, \
                      num_of_pairs=20, seqDistance=5, resi_range=None, \
                      pdbDistance=8, chain1='A', chain2='A'):
    '''
    indices and msa_indices are lists output from alignSequenceToMSA

    rank_row, rank_col and zscore_sort are the outputs from calcRankorder

    :arg num_of_pairs: The number of pairs to be output, if no value is given
        then all pairs are output. Default is 20
    :type num_of_pairs: int

    :arg seqDistance: Remove pairs that are closer than this in the reference sequence
        Default is 5
    :type seqDistance: int

    :arg pdbDistance: Remove pairs with Calpha atoms further apart than this in the PDB
        Default is 8
    :type pdbDistance: int

    :arg chain1: The chain used for the residue specified by rank_row when measuring distances
    :type chain1: str

    :arg chain2: The chain used for the residue specified by rank_col when measuring distances
    :type chain2: str
    '''

    if isscalar(indices):
        raise TypeError('Please provide a valid indices list')

    if isscalar(msa_indices):
        raise TypeError('Please provide valid msa_indices, which should be a list')

    if isscalar(rank_row):
        raise TypeError('Please provide ranked row from calcRankorder')

    if isscalar(rank_col):
        raise ValueError('Please provide ranked col from calcRankorder')

    if isscalar(zscore_sort):
        raise ValueError('Please provide sorted Z scores from calcRankorder')

    if num_of_pairs is None:
        num_of_pairs = len(rank_row)

    pairList = []
    i = -1
    j = 0
    while j < num_of_pairs:

        i += 1

        row_idx = indices[where(msa_indices == rank_row[i])[0][0]]
        col_idx = indices[where(msa_indices == rank_col[i])[0][0]]
        if not isinstance(row_idx, Integral) or not isinstance(col_idx, Integral):
            continue

        if row_idx - col_idx < seqDistance:
            continue

        distance = calcDistance(pdb.select('chain %s and resid %s' % (chain1, \
                                           row_idx)).copy(), \
                                pdb.select('chain %s and resid %s' % (chain2, \
                                           row_idx)).copy())
        if distance > pdbDistance:
            continue

        if resi_range is not None:
            if not row_idx in resi_range and not col_idx in resi_range:
                continue

        pairList.append('%3d:\t%3d\t%3d\t%5.1f\t%5.1f\n'%(i, row_idx, col_idx, zscore_sort[i], distance))

        j += 1

    return pairList

def buildSeqidMatrix(msa, turbo=True):
    """Returns sequence identity matrix for *msa*."""

    msa = getMSA(msa)

    LOGGER.timeit('_seqid')
    from .seqtools import msaeye

    dim = msa.shape[0]
    seqid = msaeye(msa, ones((dim, dim), float), turbo=bool(turbo))

    LOGGER.report('Sequence identity matrix was calculated in %.2fs.',
                  '_seqid')
    return seqid

buildSeqidMatrix.__doc__ += doc_turbo

def uniqueSequences(msa, seqid=0.98, turbo=True):
    """Returns a boolean array marking unique sequences in *msa*.  A sequence
    sharing sequence identity of *seqid* or more with another sequence coming
    before itself in *msa* will have a **True** value in the array."""

    msa = getMSA(msa)

    from .seqtools import msaeye

    if not (0 < seqid <= 1):
        raise ValueError('seqid must satisfy 0 < seqid <= 1')

    return msaeye(msa, zeros(msa.shape[0], bool),
                  unique=seqid, turbo=bool(turbo))

uniqueSequences.__doc__ += doc_turbo

def calcRankorder(matrix, zscore=False, **kwargs):
    """Returns indices of elements and corresponding values sorted in
    descending order, if *descend* is **True** (default). Can apply a zscore
    normalization; by default along *axis* - 0 such that each column has
    ``mean=0`` and ``std=1``.  If *zcore* analysis is used, return value contains the
    zscores. If matrix is symmetric only lower triangle indices will be
    returned, with diagonal elements if *diag* is **True** (default)."""

    try:
        ndim, shape = matrix.ndim, matrix.shape
    except AttributeError:
        raise TypeError('matrix must be a 2D array')

    if ndim != 2:
        raise ValueError('matrix must be a 2D array')

    threshold = kwargs.get('thredhold', 0.0001)
    try:
        symm = abs((matrix.transpose() - matrix).max()) < threshold
    except:
        symm = False

    if zscore:
        axis = int(bool(kwargs.get('axis', 0)))
        matrix = (matrix - matrix.mean(axis)) / matrix.std(axis)
        LOGGER.info('Zscore normalization has been applied.')

    descend = kwargs.get('descend', True)
    if not symm:
        if descend:
            sorted_index = matrix.argsort(axis=None)[::-1]
        else:
            sorted_index = matrix.argsort(axis=None)
        row = indices(shape)[0].flatten()[sorted_index]
        column = indices(shape)[1].flatten()[sorted_index]
    else:
        LOGGER.info('Matrix is symmetric, only lower triangle indices '
                    'will be returned.')
        if kwargs.get('diag', True):
            k = 0
        else:
            k = -1
        ind_row, ind_column = tril_indices(shape[0], k=k)
        matrix_lt = matrix[ind_row, ind_column]
        if descend:
            sorted_index = matrix_lt.argsort(axis=None)[::-1]
        else:
            sorted_index = matrix_lt.argsort(axis=None)
        row = ind_row[sorted_index]
        column = ind_column[sorted_index]

    return (row, column, matrix[row, column])

def buildOMESMatrix(msa, ambiguity=True, turbo=True, **kwargs):
    """Returns OMES (Observed Minus Expected Squared) covariance matrix
    calculated for *msa*, which may be an :class:`.MSA` instance or a 2D
    NumPy character array. OMES is defined as::

                        (N_OBS - N_EX)^2              (f_i,j - f_i * f_j)^2
      OMES_(i,j) = sum(------------------) = N * sum(-----------------------)
                             N_EX                           f_i * f_j

    Implementation is case insensitive and handles ambiguous amino acids
    as follows:

      * **B** (Asx) count is allocated to *D* (Asp) and *N* (Asn)
      * **Z** (Glx) count is allocated to *E* (Glu) and *Q* (Gln)
      * **J** (Xle) count is allocated to *I* (Ile) and *L* (Leu)
      * **X** (Xaa) count is allocated to the twenty standard amino acids
      * Joint probability of observing a pair of ambiguous amino acids is
        allocated to all potential combinations, e.g. probability of **XX**
        is allocated to 400 combinations of standard amino acids, similarly
        probability of **XB** is allocated to 40 combinations of *D* and *N*
        with the standard amino acids.

    Selenocysteine (**U**, Sec) and pyrrolysine (**O**, Pyl) are considered
    as distinct amino acids.  When *ambiguity* is set **False**, all alphabet
    characters as considered as distinct types.  All non-alphabet characters
    are considered as gaps."""

    msa = getMSA(msa)

    from .msatools import msaomes
    LOGGER.timeit('_omes')
    length = msa.shape[1]
    omes = empty((length, length), float)
    omes = msaomes(msa, omes, ambiguity=bool(ambiguity), turbo=bool(turbo),
                   debug=bool(kwargs.get('debug', False)))
    LOGGER.report('OMES matrix was calculated in %.2fs.',
                  '_omes')

    return omes

buildOMESMatrix.__doc__ += doc_turbo

def buildSCAMatrix(msa, turbo=True, **kwargs):
    """Returns SCA matrix calculated for *msa*, which may be an :class:`.MSA`
    instance or a 2D Numpy character array.

    Implementation is case insensitive and handles ambiguous amino acids
    as follows:

      * **B** (Asx) count is allocated to *D* (Asp) and *N* (Asn)
      * **Z** (Glx) count is allocated to *E* (Glu) and *Q* (Gln)
      * **J** (Xle) count is allocated to *I* (Ile) and *L* (Leu)
      * **X** (Xaa) count is allocated to the twenty standard amino acids
      * Joint probability of observing a pair of ambiguous amino acids is
        allocated to all potential combinations, e.g. probability of **XX**
        is allocated to 400 combinations of standard amino acids, similarly
        probability of **XB** is allocated to 40 combinations of *D* and *N*
        with the standard amino acids.

    Selenocysteine (**U**, Sec) and pyrrolysine (**O**, Pyl) are considered
    as distinct amino acids.  When *ambiguity* is set **False**, all alphabet
    characters as considered as distinct types.  All non-alphabet characters
    are considered as gaps."""

    msa = getMSA(msa)
    if msa.shape[0]<100:
        LOGGER.warning('SCA performs the best with higher number of sequences, and '
                       'minimal number of sequences is recommended as 100.')

    from .msatools import msasca
    LOGGER.timeit('_sca')
    length = msa.shape[1]
    sca = zeros((length, length), float)
    sca = msasca(msa, sca, turbo=bool(turbo))
    LOGGER.report('SCA matrix was calculated in %.2fs.', '_sca')
    return sca

buildSCAMatrix.__doc__ += doc_turbo

def buildPCMatrix(msa, turbo=False, **kwargs):
    """Returns PC matrix calculated for *msa*, which may be an :class:`.MSA`
    instance or a 2D Numpy character array.

    Implementation is case insensitive and handles ambiguous amino acids
    as follows:

      * **B** (Asx) count is allocated to *D* (Asp) and *N* (Asn)
      * **Z** (Glx) count is allocated to *E* (Glu) and *Q* (Gln)
      * **J** (Xle) count is allocated to *I* (Ile) and *L* (Leu)
      * **X** (Xaa) count is allocated to the twenty standard amino acids
      * Joint probability of observing a pair of ambiguous amino acids is
        allocated to all potential combinations, e.g. probability of **XX**
        is allocated to 400 combinations of standard amino acids, similarly
        probability of **XB** is allocated to 40 combinations of *D* and *N*
        with the standard amino acids.

    Selenocysteine (**U**, Sec) and pyrrolysine (**O**, Pyl) are considered
    as distinct amino acids.  When *ambiguity* is set **False**, all alphabet
    characters as considered as distinct types.  All non-alphabet characters
    are considered as gaps.
    """

    msa = getMSA(msa)
    from .msatools import msapsicov
    LOGGER.timeit('_psicov')
    length = msa.shape[1]
    pc = zeros((length, length), float)
    pc = msapsicov(msa, pc, turbo=bool(turbo))
    LOGGER.report('PC matrix was calculated in %.2fs.', '_psicov')
    return pc

def buildDirectInfoMatrix(msa, seqid=.8, pseudo_weight=.5, refine=False,
                          **kwargs):
    """Returns direct information matrix calculated for *msa*, which may be an
    :class:`.MSA` instance or a 2D Numpy character array.

    Sequences sharing sequence identity of *seqid* or more with another
    sequence are regarded as similar sequences for calculating their weights
    using :func:`.calcMeff`.

    *pseudo_weight* are the weight for pseudo count probability.

    Sequences are not refined by default. When *refine* is set **True**,
    the MSA will be refined by the first sequence and the shape of direct
    information matrix will be smaller.
    """

    msa = getMSA(msa)
    from .msatools import msadipretest, msadirectinfo1, msadirectinfo2
    from numpy import matrix

    LOGGER.timeit('_di')
    if msa.shape[0]<250:
        LOGGER.warning('DI performs the best with higher number of sequences, and '
                       'minimal number of sequences is recommended as 250.')
    refine = 1 if refine else 0
    # msadipretest get some parameter from msa to set matrix size
    length, q = msadipretest(msa, refine=refine)
    c = matrix.dot(matrix(zeros((length*q, 1), float)),
                   matrix(zeros((1, length*q), float)))
    prob = zeros((length, q+1), float)
    # msadirectinfo1 return c to be inversed and prob to be used
    meff, n, length, c, prob = msadirectinfo1(msa, c, prob, theta=1.-seqid,
                                              pseudocount_weight=pseudo_weight,
                                              refine=refine, q=q+1)

    c = c.I

    di = zeros((length, length), float)
    # get final DI
    di = msadirectinfo2(n, length, c, prob, di, q+1)
    del prob, c
    LOGGER.report('DI matrix was calculated in %.2fs.', '_di')
    return di

def calcMeff(msa, seqid=.8, refine=False, weight=False, **kwargs):
    """Returns the Meff for *msa*, which may be an :class:`.MSA`
    instance or a 2D Numpy character array.

    Since similar sequences in an *msa* decreases the diversity of *msa*,
    *Meff* gives a weight for sequences in the *msa*.

    For example: One sequence in MSA has 5 other similar sequences in this
    MSA(itself included). The weight of this sequence is defined as 1/5=0.2.
    Meff is the sum of all sequence weights. In another word, Meff can be
    understood as the effective number of independent sequences.

    Sequences sharing sequence identity of *seqid* or more with another
    sequence are regarded as similar sequences to calculate Meff.

    Sequences are not refined by default. When *refine* is set **True**, the
    MSA will be refined by the first sequence.

    The weight for each sequence are returned when *weight* is **True**."""

    msa = getMSA(msa)
    from .msatools import msameff
    LOGGER.timeit('_meff')
    refine = 1 if refine else 0
    weight = 0 if weight else 1  # A Mark for return weighted array.
    if (not weight):
        w = zeros((msa.shape[0]), float)
        meff = msameff(msa, theta=1.-seqid, meff_only=weight,
                       refine=refine, w=w)
    else:
        meff = msameff(msa, theta=1.-seqid, meff_only=weight, refine=refine)
    LOGGER.report('Meff was calculated in %.2fs.', '_meff')
    return meff

def alignSequencesByChain(PDBs, **kwargs):
    """
    Runs :func:`buildMSA` for each chain and optionally joins the results.
    Returns either a single :class:`MSA` or a dictionary containing an :class:`MSA` for each chain.

    :arg PDBs: a list of :class:`AtomGroup` objects
    :type PDBs: list

    :arg join_chains: whether to join chain alignments
        default is True
    :type join_chains: bool

    :arg join_char: a character for joining chain alignments
        default is '/' as used by PIR format alignments
    :type join_char: str
    """

    if isscalar(PDBs):
        raise TypeError('PDBs should be array-like')

    if not PDBs:
        raise ValueError('PDBs should not be empty')

    pdbs = []
    chains = []
    for i, pdb in enumerate(PDBs):
        if isinstance(pdb, Atomic):
            pdbs.append(pdb)
        else:
            raise TypeError('each entry in PDBs must be a :class:`Atomic` instance')

        chains.append([])
        for chain in list(pdbs[i].getHierView()):
            chains[i].append(chain)

        if i != 0 and len(chains[i]) != len(chains[0]):
            raise ValueError('all pdbs should have the same number of chains')

    labels = []
    for pdb in pdbs:
        chids = ''
        for chain in list(pdb.getHierView()):
            chids += chain.getChid()
        labels.append(pdb.getTitle() + '_' + chids)

    chains = array(chains)
    chain_alignments = []
    alignments = {}
    for j in range(len(chains[0])):
        prefix = 'chain_' + chains[0, j].getChid()
        msa = buildMSA(chains[:, j], title=prefix, labels=labels)
        msa = refineMSA(msa, colocc=1e-9) # remove gap-only cols

        chain_alignments.append(msa)
        alignments[labels[0].split('_')[1][j]] = msa

    join_chains = kwargs.get('join_chains', True)
    join_char = kwargs.get('join_char', '/')

    if len(chains[0]) == 1:
        join_chains = False

    if join_chains:
        joined_msaarr = []
        for i, chain_alignment in enumerate(chain_alignments):
            pdb_seqs = []
            for j, sequence in enumerate(chain_alignment):
                pdb_seqs.append(sequence)
            joined_msaarr.append(join_char.join(pdb_seqs))

        result = MSA(joined_msaarr, title='joined_chains',
                     labels=[label.split('_')[0] for label in labels])

    else:
        result = alignments
        if len(result) == 1:
            result = result[list(result.keys())[0]]

    return result

def buildMSA(sequences, title='Unknown', labels=None, **kwargs):
    """
    Aligns sequences with clustalw or clustalw2 and returns the resulting MSA.

    :arg sequences: a file, MSA object or a list or array containing sequences
       as Atomic objects with :func:`getSequence` or Sequence objects or strings.
       If strings are used then labels must be provided using ``labels``
    :type sequences: :class:`Atomic`, :class:`.MSA`,
        :class:`~numpy.ndarray`, str

    :arg title: the title for the MSA and it will be used as the prefix for output files.
    :type title: str

    :arg labels: a list of labels to go with the sequences
    :type labels: list

    :arg align: whether to align the sequences
        default True
    :type align: bool

    :arg method: alignment method, one of either 'global' (biopython.pairwise2.align.globalms),
        'local' (biopython.pairwise2.align.localms), clustalw(2), or another software in your path.
        Default is 'clustalw'
    :type align: str
    """

    align = kwargs.get('align', True)
    method = kwargs.pop('method', 'clustalw')
    # 1. check if sequences are in a fasta file and if not make one
    if isinstance(sequences, str):
        filename = sequences
    elif not isinstance(sequences, MSA):
        try:
            max_len = 0
            for sequence in sequences:
                if isinstance(sequence, Atomic):
                    if len(sequence.ca.copy()) > max_len:
                        max_len = len(sequence.ca.copy())
                elif isinstance(sequence, MSA):
                    if len(sequence[0]) > max_len:
                        max_len = len(sequence[0])
                else:
                    if len(sequence) > max_len:
                        max_len = len(sequence)

            msa = []
            fetched_labels = []
            for i, sequence in enumerate(sequences):
                if isinstance(sequence, Atomic):
                    strseq = sequence.ca.getSequence()
                    label = sequence.getTitle()
                elif isinstance(sequence, Sequence):
                    strseq = str(sequence)
                    label = sequence.getLabel()
                elif isinstance(sequence, MSA):
                    strseq = str(sequence[0])
                    label = sequence.getLabel(0)
                    LOGGER.warn('Only the first sequence in the MSA at entry {0} is used.'
                                .format(i))
                elif isinstance(sequence, str):
                    strseq = sequence
                    label = str(i + 1)
                else:
                    raise TypeError('sequences should be a list of strings, '
                                    'Atomic, or Sequence instances')
                strseq = strseq + '-'*(max_len - len(strseq))
                msa.append(array(list(strseq)))
                fetched_labels.append(label)
            sequences = array(msa)
        except:
            raise TypeError('sequences should be iterable')

        # "if a list" is a pythonic way to check if a list is empty or not (or none)
        if not labels and fetched_labels:
            labels = fetched_labels

        label = [label.replace(' ','_') for label in labels]
        # labels checkers are removed because they will be properly handled in MSA class initialization
        msa = MSA(msa=sequences, title=title, labels=labels)

        if align and 'clustal' in method:
            filename = writeMSA(title + '.fasta', msa)

    if align:
        # 2. find and run alignment method
        if method in ['biopython', 'local', 'global']:
            if len(sequences) == 2:
                msa, _, _ = alignTwoSequencesWithBiopython(sequences[0], sequences[1], **kwargs)
            else:
                raise ValueError("Provide only two sequences or another method. \
                                  Biopython pairwise alignment can only be used \
                                  to build an MSA with two sequences.")
        elif 'clustalw' in method:
            clustalw = which('clustalw')
            if clustalw is None:
                if which('clustalw2') is not None:
                    clustalw = which('clustalw2')
                else:
                    raise EnvironmentError("The executable for clustalw was not found, \
                                            install clustalw or add it to the path.")

            os.system('"%s" %s -OUTORDER=INPUT'%(clustalw, filename))

            # 3. parse and return the new MSA
            msa = parseMSA(title + '.aln')

        else:
            alignTool = which(method)
            if alignTool is None:
                raise EnvironmentError("The executable for {0} was not found, \
                                        install it or add it to the path.".format(alignTool))

            os.system('"%s" %s -OUTORDER=INPUT'%(clustalw, filename))

            # 3. parse and return the new MSA
            msa = parseMSA(title + '.aln')

    return msa

def showAlignment(alignment, **kwargs):
    """
    Prints out an alignment as sets of short rows with labels.

    :arg alignment: any object with aligned sequences
    :type alignment: :class: `.MSA`, list

    :arg row_size: the size of each row
        default 60
    :type row_size: int

    :arg indices: a set of indices for some or all sequences
        that will be shown above the relevant sequences
    :type indices: :class:`~numpy.ndarray`, list

    :arg index_start: how far along the alignment to start putting indices
        default 0
    :type index_start: int

    :arg index_stop: how far along the alignment to stop putting indices
        default the point when the shortest sequence stops
    :type index_stop: int

    :arg labels: a list of labels
    :type labels: list
    """
    row_size = kwargs.get('row_size', 60)

    labels = kwargs.get('labels', None)
    if labels is not None:
        if isscalar(labels):
            raise TypeError('labels should be array-like')

        for label in labels:
            if not isinstance(label, str):
                raise TypeError('each label should be a string')

        if len(labels) < len(alignment):
            raise ValueError('there should be a label for every sequence shown')
    else:
        labels = []
        for i, sequence in enumerate(alignment):
            if hasattr(sequence, 'getLabel'):
                labels.append(sequence.getLabel())
            else:
                labels.append(str(i+1))

    indices = kwargs.get('indices', None)
    index_start = kwargs.get('index_start', 0)
    index_stop = kwargs.get('index_stop', 0)

    if index_stop == 0 and indices is not None:
        locs = []
        maxes = []
        for index in indices:
            int_index = []
            for i in index:
                if i == '':
                    int_index.append(0)
                else:
                    int_index.append(int(i))
            int_index = array(int_index)
            maxes.append(max(int_index))
            locs.append(where(int_index == max(int_index))[0][0])
        index_stop = locs[where(maxes == min(maxes))[0][0]]

    for i in range(int(ceil(len(alignment[0])/float(row_size)))):
        for j in range(len(alignment)):
            if indices is not None:
                sys.stdout.write('\n' + ' '*15 + '\t')
                for k in range(row_size*i+10,row_size*(i+1)+10,10):
                    try:
                        if k > index_start + 10 and k < index_stop + 10:
                            sys.stdout.write('{:10d}'.format(int(indices[j][k-1])))
                        elif k < index_stop:
                            sys.stdout.write(' '*(k-index_start))
                        else:
                            sys.stdout.write(' '*10)
                    except:
                            sys.stdout.write(' '*10)
                sys.stdout.write('\n')

            sys.stdout.write(labels[j][:15] + \
                             ' ' * (15-len(labels[j][:15])) + \
                             '\t' + str(alignment[j])[60*i:60*(i+1)] + '\n')

        sys.stdout.write('\n')

    return

def alignSequenceToMSA(seq, msa, **kwargs):
    """
    Align a sequence from a PDB or Sequence to a sequence from an MSA
    and create two sets of indices.

    The sequence from the MSA (*seq*), the alignment and
    the two sets of indices are returned.

    The first set (*indices*) maps the residue numbers in the PDB to
    the reference sequence. The second set (*msa_indices*) indexes the
    reference sequence in the msa and is used for retrieving values
    from the first indices.

    :arg seq: an object with an associated sequence string
         or a sequence string itself
    :type seq: :class:`.Atomic`, :class:`.Sequence`, str

    :arg msa: a multiple sequence alignment
    :type msa: :class:`.MSA`

    :arg label: a label for a sequence in msa or a PDB ID
        ``msa.getIndex(label)`` must return a sequence index
    :type label: str

    :arg chain: which chain from pdb to use for alignment, default is **None**,
        which does no selection on *seq*. This value will be ignored if seq is
        not an :class:`.Atomic` object.
    :type chain: str

    Parameters for Biopython ``pairwise2`` alignments can be provided as
    keyword arguments. Default values are originally from ``proteins.compare``
    module, but now found in ``utilities.seqtools``.

    :arg match: a positive integer, used to reward finding a match
    :type match: int

    :arg mismatch: a negative integer, used to penalise finding a mismatch
    :type mismatch: int

    :arg gap_opening: a negative integer, used to penalise opening a gap
    :type gap_opening: int

    :arg gap_extension: a negative integer, used to penalise extending a gap
    :type gap_extension: int

    :arg method: method for pairwise2 alignment.
        Possible values are ``"local"`` and ``"global"``
    :type method: str
    """
    label = kwargs.get('label', None)
    chain = kwargs.get('chain', None)

    match = kwargs.get('match', MATCH_SCORE)
    mismatch = kwargs.get('mismatch', MISMATCH_SCORE)
    gap_opening = kwargs.get('gap_opening', GAP_PENALTY)
    gap_extension = kwargs.get('gap_extension', GAP_EXT_PENALTY)
    method = kwargs.get('method', ALIGNMENT_METHOD)

    if isinstance(seq, Atomic):
        if isinstance(chain, str):
            ag = seq.select('chain {0}'.format(chain))
        elif chain is None:
            ag = seq

            chids = ag.getChids()
            if len(unique(chids)) > 1:
                LOGGER.warn('%s consists of multiple chains. Please consider selecting one chain'%(seq.getTitle()))
        else:
            raise TypeError('chain should be a string or None')

        if ag is None:
            raise ValueError('seq may be None or chain ID may be invalid')
        sequence = ag.select('ca').getSequence()

    elif isinstance(seq, Sequence):
         sequence = str(seq)
         ag = None
    elif isinstance(seq, str):
        sequence = seq
        ag = None
    else:
        raise TypeError('seq must be an atomic class, sequence class, or str not {0}'
                        .format(type(seq)))

    if not isinstance(msa, MSA):
        raise TypeError('msa must be an MSA instance')

    if label is None:
        if ag:
            label = ag.getTitle().split('_')[0]
        elif isinstance(seq, Sequence):
            label = seq.getLabel()
        else:
            raise ValueError('A label cannot be extracted from seq so please provide one.')

    index = msa.getIndex(label)

    if index is None and (len(label) == 4 or len(label) == 5):
        from prody import parsePDB
        try:
            structure, header = parsePDB(label[:4], header=True)
        except Exception as err:
            raise IOError('failed to parse header for {0} ({1})'
                            .format(label[:4], str(err)))

        chid = chain
        for poly in header['polymers']:
            if chid and poly.chid != chid:
                continue
            for dbref in poly.dbrefs:
                if index is None:
                    index = msa.getIndex(dbref.idcode)
                    if index is not None:
                        LOGGER.info('{0} idcode {1} for {2}{3} '
                                    'is found in {4}.'.format(
                                    dbref.database, dbref.idcode,
                                    label[:4], poly.chid, str(msa)))
                        label = dbref.idcode
                        break
                if index is None:
                    index = msa.getIndex(dbref.accession)
                    if index is not None:
                        LOGGER.info('{0} accession {1} for {2}{3} '
                                    'is found in {4}.'.format(
                                    dbref.database, dbref.accession,
                                    label[:4], poly.chid, str(msa)))
                        label = dbref.accession
                        break
        if index is not None:
            chain = structure[poly.chid]

    if index is None:
        raise ValueError('label is not in msa, or msa is not indexed')
    try:
        len(index)
    except TypeError:
        pass
    else:
        raise ValueError('label {0} maps onto multiple sequences, '
                            'so cannot be used for refinement'.format(label))

    if isinstance(index, int):
        refMsaSeq = str(msa[index]).upper().replace('-','.')
    else:
        raise TypeError('The output from querying that label against msa is not a single sequence.')

    if method == 'local':
        alignment = pairwise2.align.localms(sequence, str(refMsaSeq),
                                            match, mismatch, gap_opening, gap_extension,
                                            one_alignment_only=1)
    elif method == 'global':
        alignment = pairwise2.align.globalms(sequence, str(refMsaSeq),
                                       match, mismatch, gap_opening, gap_extension,
                                       one_alignment_only=1)
    else:
        raise ValueError('method should be local or global')

    seq_indices = [0]
    msa_indices = [0]

    for i in range(len(alignment[0][0])):
        if alignment[0][0][i] != '-':
            seq_indices.append(seq_indices[i]+1)
        else:
            seq_indices.append(seq_indices[i])

        if alignment[0][1][i] != '-':
            msa_indices.append(msa_indices[i]+1)
        else:
            msa_indices.append(msa_indices[i])

    seq_indices.pop(0) # The first element was extra for initialisation
    msa_indices.pop(0) # The first element was extra for initialisation

    seq_indices = array(seq_indices)
    msa_indices = array(msa_indices)

    if ag:
        seq_indices += ag.getResnums()[0] - 1

    alignment = MSA(msa=array([array(list(alignment[0][0])), \
                               array(list(alignment[0][1]))]), \
                    labels=[ag.getTitle(), label])

    return alignment, seq_indices, msa_indices

def alignTwoSequencesWithBiopython(seq1, seq2, **kwargs):
    """Easily align two sequences with Biopython's globalms or localms.
    Returns an MSA and indices for use with :func:`.showAlignment`.

    Alignment parameters can be provided as keyword arguments.
    Default values are as originally set in the proteins.compare
    module, but now found in utilities.seqtools.

    :arg match: a positive integer, used to reward finding a match
    :type match: int

    :arg mismatch: a negative integer, used to penalise finding a mismatch
    :type mismatch: int

    :arg gap_opening: a negative integer, used to penalise opening a gap
    :type gap_opening: int

    :arg gap_extension: a negative integer, used to penalise extending a gap
    :type gap_extension: int

    :arg method: method for pairwise2 alignment.
        Possible values are 'local' and 'global'
    :type method: str
    """

    match = kwargs.get('match', MATCH_SCORE)
    mismatch = kwargs.get('mismatch', MISMATCH_SCORE)
    gap_opening = kwargs.get('gap_opening', GAP_PENALTY)
    gap_extension = kwargs.get('gap_extension', GAP_EXT_PENALTY)
    method = kwargs.get('method', ALIGNMENT_METHOD)

    if method == 'local':
        alignment = pairwise2.align.localms(seq1, seq2, match, mismatch, gap_opening, gap_extension)
    elif method == 'global':
        alignment = pairwise2.align.globalms(seq1, seq2, match, mismatch, gap_opening, gap_extension)
    else:
        raise ValueError('method should be local or global')

    seq_indices = [0]
    msa_indices = [0]

    for i in range(len(alignment[0][0])):
        if alignment[0][0][i] != '-':
            seq_indices.append(seq_indices[i]+1)
        else:
            seq_indices.append(seq_indices[i])

        if alignment[0][1][i] != '-':
            msa_indices.append(msa_indices[i]+1)
        else:
            msa_indices.append(msa_indices[i])

    seq_indices = array(seq_indices)
    msa_indices = array(msa_indices)

    alignment = MSA(msa=array([array(list(alignment[0][0])), \
                               array(list(alignment[0][1]))]))

    return alignment, seq_indices, msa_indices

def trimAtomsUsingMSA(atoms, msa, **kwargs):
    """This function uses :func:`.alignSequenceToMSA` and has the same kwargs.

    :arg atoms: an atomic structure for trimming
    :type atoms: :class:`.Atomic`

    :arg msa: a multiple sequence alignment
    :type msa: :class:`.MSA`
    """
    aln, idx_1, idx_2 = alignSequenceToMSA(atoms, msa, **kwargs)

    u, i = unique(idx_2, return_index=True)

    resnums_str = ' '.join([str(x) for x in idx_1[i]])
    chain = kwargs.get('chain', 'A')

    return atoms.select('chain {0} and resnum {1}'.format(chain, resnums_str))

max_stars_repo_name: kaynakb/ProDy

max_issues_repo_name: kaynakb/ProDy

max_forks_repo_name: kaynakb/ProDy

reason: from numpy


text:
"""
System expansion, reduction, and corrections functions for Sparse systems.

@author: Joseph C. Slater and Sainag Immidisetty
"""
__license__ = "Joseph C. Slater"

__docformat__ = 'reStructuredText'

#import math
import numpy as np
#import scipy.signal as sig
import scipy.linalg as la
import scipy.sparse.linalg as spla
from scipy.sparse import lil_matrix
#from scipy.sparse import csr_matrix
import scipy.sparse as sps

def sos_modal_forsparse(M, K, C=False, damp_diag=0.03, shift=1):
	r"""Eigen analysis of proportionally damped system.

	Optimally find mass normalized mode shapes and natural frequencies
	of a system modelled by :math:`M\ddot{x}+Kx=0`.

	If provided, obtain damping ratios presuming :math:`C` can be decoupled.

	Provides a warning if diagonalization of damping matrix fails worse than
	relative error of `damp_diag`.

	Parameters
	----------
	M, K : float arrays
		Mass and stiffness matrices
	C : float array, optional
		Damping matrix
	damp_diag : float, optional
		Maximum amount of off-diagonal error allowed in assuming C can be
		diagonalized
	shift : float, optional
		Shift used in eigensolution. Should be approximately equal to the first
		non-zero eigenvalue.

	Returns
	-------
	omega : float array (1xN)
		Vector of natural frequencies (rad/sec)
	zeta : float array (1xN)
		Vector of damping ratios
	Psi : float array (NxN)
		Matrix of mass normalized mode shapes by column

	Examples
	--------

	"""

	#Kdiag = K.diagonal().reshape(-1,1)

	#Mdiag = M.diagonal().reshape(-1,1)

	#minvals = np.sort((Kdiag/Mdiag),axis=0)

	#shift = minvals[min(7,len(minvals))]

	#shift = shift[0]

	#K = ((K.tocsr() + (K.T).tocsr()).tolil()) * 0.5 + shift * ((M.tocsr() + (M.T).tocsr()).tolil()) * 0.5

	#M = ((M.tocsr() + (M.T).tocsr()).tolil()) * 0.5

	K = lil_matrix(K)

	M = lil_matrix(M)

	K = K + M * shift

	[lam, Psi] = la.eigh(K.toarray(), M.toarray())

	lam = lam.reshape(-1,1)

	omega = np.sqrt(np.abs(lam - shift))

	omega = omega.reshape(-1,)

	norms = np.diag(1.0 / np.sqrt(np.diag(Psi.T@M@Psi)))

	Psi = Psi @ norms

	zeta = np.zeros_like(omega)

	if C is not False:
		diagonalized_C = Psi.T@C@Psi

		diagonal_C = np.diag(diagonalized_C)

		if min(omega) > 1e-5:
			zeta = diagonal_C / 2 / omega  # error if omega = 0
			max_off_diagonals = np.amax(np.abs(diagonalized_C
											   - np.diag(diagonal_C)), axis=0)
			# error if no damping
			damp_error = np.max(max_off_diagonals / diagonal_C)
		else:
			zeta = np.zeros_like(omega)
			damp_error = np.zeros_like(omega)
			de_diag_C = diagonalized_C - np.diag(diagonal_C)
			for mode_num, omega_i in enumerate(omega):
				if omega[mode_num] > 1e-5:
					zeta[mode_num] = diagonal_C[mode_num] / 2 / omega_i
					damp_error = (np.max(np.abs(de_diag_C[:, mode_num]))
								  / diagonal_C[mode_num])

		if damp_error > damp_diag:
			print('Damping matrix cannot be completely diagonalized.')
			print('Off diagonal error of {:4.0%}.'.format(damp_error))

	return omega, zeta, Psi

def guyan_forsparse(M, K, master=None, fraction=None):
	r"""Guyan reduced model.

	Applies Guyan Reductions to second order system of equations of the form

	.. math:: M \ddot{x} + K x = 0

	which are reduced to the form

	.. math:: M_r \ddot{x}_m + K_r x_m = 0

	where

	.. math::

		x = T x_m

		M_r= T^T M T

		K_r= T^T K T

	Parameters
	----------
	M, K : float arrays
		Mass and Stiffness matrices
	master : float array or list, optional
		List of retained degrees of freedom (0 indexing)
	fraction : float, optional
		Fraction of degrees of freedom (0< `fraction` <1.0) to retain in model.
		If both master and fraction are neglected, fraction is set to 0.25.

	Returns
	-------
	Mred, Kred, T : float arrays
		Reduced Mass matric, Reduced Stiffness matrix, Transformation matrix
	master_dofs : int list
		List of master degrees of freedom (0 indexing)
	truncated_dofs : int list
		List of truncated degrees of freedom (0 indexing)

	Examples
	--------

	Notes
	-----
	Reduced coordinate system forces can be obtained by `Fr = T.T @ F`.

	Reduced damping matrix can be obtained using `Cr = T.T @ C @ T`.

	If mode shapes are obtained for the reduced system, full system mode shapes
	are `phi = T @ phi_r`.

	Code is not as efficient as possible. Using block submatrices would be
	more efficient.

	"""

	if master is None:
		if fraction is None:
			fraction = 0.25

		ratios = np.diag(M) / np.diag(K)
		ranked = [i[0] for i in sorted(enumerate(ratios), key=lambda x: x[1])]
		thresh = int(fraction * ratios.size)
		if (thresh >= ratios.size) or thresh == 0:
			print("Can't keep", thresh, 'DOFs.')
			print("Fraction of", fraction, "is too low or too high.")
			return 0, 0, 0, 0, 0

		master = ranked[-thresh:]

	master = np.array(master)

	ncoord = M.shape[0]

	i = np.arange(0, ncoord)

	i = i.reshape(1,-1)

	i = i + np.ones((1,i.shape[1]),int)

	lmaster = master.shape[1]

	i[0,master-1] = np.transpose(np.zeros((lmaster,1)))

	i = np.sort((i), axis =1)

	slave = i[0,lmaster + 0:ncoord]

	K= lil_matrix(K)

	slave = slave.reshape(1,-1)

	master = master-np.ones((1,master.shape[0]),int)

	master = master.ravel()

	slave = slave - np.ones((1,slave.shape[0]),int)

	slave = slave.ravel()

	kss = slice_forSparse(K, slave, slave)

	ksm = slice_forSparse(K, slave, master)

	T= np.zeros((len(master)+len(slave), len(master)))

	T= lil_matrix(T)

	T[master,:lmaster] = sps.eye(lmaster,lmaster)

	T[slave,0:lmaster]=spla.spsolve(-kss,ksm)

	Mred = T.T * M * T

	Kred = T.T * K * T

	return Mred, Kred, master

def mode_expansion_from_model_forsparse(Psi, omega, M, K, measured):
    r"""Deflection extrapolation to full FEM model coordinates, matrix method.

    Provided an equation  of the form:

    :math:`\begin{pmatrix}-\begin{bmatrix}M_{mm}&M_{mu}\\ M_{um}&M_{uu}
    \end{bmatrix} \omega_i^2
    +\begin{bmatrix}K_{mm}&K_{mu}\\ K_{um}&K_{uu}\end{bmatrix}\end{pmatrix}`
    :math:`\begin{bmatrix}\Psi_{i_m}\\ \Psi_{i_u}\end{bmatrix}= 0`

    Where:

    - :math:`M` and :math:`K` are the mass and stiffness matrices, likely from
      a finite element model
    - :math:`\Psi_i` and :math:`\omega_i` represent a mode/frequency pair
    - subscripts :math:`m` and :math:`u` represent measure and unmeasured
      of the mode

    Determines the unknown portion of the mode (or operational deflection)
    shape, :math:`\Psi_{i_u}` by
    direct algebraic solution, aka

    :math:`\Psi_{i_u} = - (K_{uu}- M_{ss} \omega_i^2) ^{-1}
    (K_{um}-M_{um}\omega_i^2)\Psi_{i_m}`

    Parameters
    ----------
    Psi : float array
        mode shape or shapes, 2-D array columns of which are mode shapes
    omega : float or 1-D float array
        natural (or driving) frequencies
    M, K : float arrays
        Mass and Stiffness matrices
    measured : float or integer array or list
        List of measured degrees of freedom (0 indexed)

    Returns
    -------
    Psi_full : float array
        Complete mode shape

    Examples
    --------
    >>> import vibrationtesting as vt
    >>> M = np.array([[4, 0, 0],
    ...               [0, 4, 0],
    ...               [0, 0, 4]])
    >>> K = np.array([[8, -4, 0],
    ...               [-4, 8, -4],
    ...               [0, -4, 4]])
    >>> measured = np.array([[0, 2]])
    >>> omega, zeta, Psi = vt.sos_modal(M, K)
    >>> Psi_measured = np.array([[-0.15], [-0.37]])
    >>> Psi_full = vt.mode_expansion_from_model(Psi_measured, omega[0], M, K,
    ... measured)
    >>> print(np.hstack((Psi[:,0].reshape(-1,1), Psi_full)))
    [[-0.164  -0.15  ]
     [-0.2955  0.2886]
     [-0.3685 -0.37  ]]

    Notes
    -----
    .. seealso:: incomplete multi-mode update. Would require each at a
      different frequency.

    """

    measured = measured.reshape(-1)  # retained dofs
    num_measured = len(measured)
    ndof = int(M.shape[0])  # length(M);
    unmeasured_dofs = list(set(np.arange(ndof)) - set(measured))
    num_unmeasured = len(unmeasured_dofs)

    M= lil_matrix(M)

    K= lil_matrix(K)

    Muu = slice_forSparse(M, unmeasured_dofs, unmeasured_dofs)

    Kuu = slice_forSparse(K, unmeasured_dofs, unmeasured_dofs)

    Mum = slice_forSparse(M, unmeasured_dofs, measured)

    Kum = slice_forSparse(K, unmeasured_dofs, measured)

    if isinstance(omega, float):
        omega = np.array(omega).reshape(1)

    Psi_full = np.zeros((num_measured + num_unmeasured, Psi.shape[1]))
    Psi_full[measured] = Psi

    for i, omega_n in enumerate(omega):
        Psi_i = Psi[:, i].reshape(-1, 1)
        Psi_unmeasured = la.solve((Kuu - Muu * omega_n**2),
                                  (Kum - Mum * omega_n**2)@Psi_i)
        Psi_unmeasured = Psi_unmeasured.reshape(-1, )
        Psi_full[unmeasured_dofs, i] = Psi_unmeasured
        # Psi_full = Psi_full.reshape(-1, 1)
    return Psi_full

def slice_forSparse(Matrix, a, b):
	"""

	Parameters
	----------
	Matrix : float array
		Arbitrary array
	a, b : int lists or arrays
		list of rows and columns to be selected from `Matrix`

	Returns
	-------
	Matrix : float array

	"""
	Maa = Matrix[a,:].toarray()

	Maa = Maa[:, a]

	Mab = Matrix[a,:].toarray()

	Mab = Mab[:, b]

	if(np.array_equiv(a,b)):
		return Maa

	else:
		return Mab

max_stars_repo_name: Vibration-Testing/vibrationtesting

max_issues_repo_name: bagustris/vibrationtesting

max_forks_repo_name: bagustris/vibrationtesting

reason: import numpy,import scipy,from scipy


text:
import numpy as np
import cv2
from matplotlib import pyplot as plt

img = cv2.imread('9.jpg', 0)

dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)

magnitude_spectrum = 20 * np.log(cv2.magnitude(dft_shift[:, :, 0], dft_shift[:, :, 1]))

plt.figure(figsize=(10, 10))
plt.subplot(221), plt.imshow(img, cmap='gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(222), plt.imshow(magnitude_spectrum, cmap='gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])

rows, cols = img.shape
crow, ccol = rows // 2, cols // 2
# create a mask first, center square is 1, remaining all zeros
mask = np.zeros((rows, cols, 2), np.uint8)
mask[crow - 30:crow + 30, ccol - 30:ccol + 30] = 1
fshift = dft_shift * mask
# apply mask and inverse DFT

f_ishift = np.fft.ifftshift(fshift)
img_back = cv2.idft(f_ishift)
img_back = cv2.magnitude(img_back[:, :, 0], img_back[:, :, 1])

plt.subplot(223), plt.imshow(img_back, cmap='gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(224), plt.imshow(img_back)
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])

plt.show()

max_stars_repo_name: swpucwf/Deeplearning

max_issues_repo_name: swpucwf/deeplearning

max_forks_repo_name: swpucwf/deeplearning


text:
#!/usr/bin/python
"""
Plot parameter 1D marginal constraints as a fn. of FG subtraction efficiency.
"""
import numpy as np
import pylab as P
from rfwrapper import rf

nsurveys = 2 #5
nbins = 14

name = ['SD', 'Interferom.']
#name = ['A', 'B', 'C', 'D', 'E']
cols = ['b', 'g', 'c', 'r', 'y']

P.subplot(111)
for j in range(nsurveys):
    for i in range(1): #range(nbins): # FIXME: Just z=0.1
        #if j == 0: continue
        fname_smooth = "tmp/pk_surveys_smooth_powspec_"+str(j)+"_"+str(i)+".npy"
        fname_constr = "tmp/pk_surveys_constraints_"+str(j)+"_"+str(i)+".npy"

        k, pk, fbao, pksmooth = np.load(fname_smooth).T
        kc, pkc, pkerr, fbao_c, pkcsmooth = np.load(fname_constr).T

        pkerr[np.where(np.isinf(pkerr))] = 1e9
        #pkerr[np.where(np.isnan(pkerr))] = 1e9
        P.plot(kc, pkerr, color=cols[j], label=name[j], marker='.', lw=1.5)

        print pkerr

        # FIXME: h^-1 units?
        #P.plot(k, pk)
        #yup, ylow = rf.fix_log_plot(pkc, pkerr*pkc)
        #P.errorbar(kc, pkc, yerr=[ylow, yup], marker='.', ls='none', color='r')

P.xscale('log')
P.yscale('log')

P.xlim((4e-3, 2.))
#P.ylim((3e-3, 1e1))

# Display options
P.legend(loc='upper left', prop={'size':'x-large'}, ncol=2)
P.ylabel("$\Delta P(k) / P(k) \; (z=0)$", fontdict={'fontsize':'22'})
P.xlabel("$k \,[Mpc^{-1}]$", fontdict={'fontsize':'20'})

fontsize = 18.
for tick in P.gca().yaxis.get_major_ticks():
  tick.label1.set_fontsize(fontsize)
for tick in P.gca().xaxis.get_major_ticks():
  tick.label1.set_fontsize(fontsize)
P.tight_layout()
P.show()

max_stars_repo_name: sjforeman/RadioFisher

max_issues_repo_name: sjforeman/RadioFisher

max_forks_repo_name: sjforeman/RadioFisher