text:
"""
Library for doing sequence design that can be expressed as linear algebra
operations for rapid processing by numpy (e.g., generating all DNA sequences
of a certain length and calculating all their full duplex binding energies
in the nearest neighbor model and filtering those outside a given range).

Based on the DNA single-stranded tile (SST) sequence designer used in the following publication.

"Diverse and robust molecular algorithms using reprogrammable DNA self-assembly"
Woods\*, Doty\*, Myhrvold, Hui, Zhou, Yin, Winfree. (\*Joint first co-authors)
"""  # noqa

from __future__ import annotations

from typing import Tuple, List, Collection, Optional, Union, Sequence, Dict, Iterable
from dataclasses import dataclass
import math
import itertools as it
from functools import lru_cache

import numpy as np

default_rng: np.random.Generator = np.random.default_rng()  # noqa

bits2base = ['A', 'C', 'G', 'T']
base2bits = {'A': 0b00, 'C': 0b01, 'G': 0b10, 'T': 0b11,
             'a': 0b00, 'c': 0b01, 'g': 0b10, 't': 0b11}

def idx2seq(idx: int, length: int) -> str:
    """Return the lexicographic idx'th DNA sequence of given length."""
    seq = ['x'] * length
    for i in range(length - 1, -1, -1):
        seq[i] = bits2base[idx & 0b11]
        idx >>= 2
    return ''.join(seq)

def seq2arr(seq: str, base2bits_local: Optional[Dict[str, int]] = None) -> np.ndarray:
    """Convert seq (string with DNA alphabet) to numpy array with integers 0,1,2,3."""
    if base2bits_local is None:
        base2bits_local = base2bits
    return np.array([base2bits_local[base] for base in seq], dtype=np.ubyte)

def seqs2arr(seqs: Sequence[str]) -> np.ndarray:
    """Return numpy 2D array converting the given DNA sequences to integers."""
    if len(seqs) == 0:
        return np.empty((0, 0), dtype=np.ubyte)
    seq_len = len(seqs[0])
    for seq in seqs:
        if len(seq) != seq_len:
            raise ValueError('All sequences in seqs must be equal length')
    num_seqs = len(seqs)
    arr = np.empty((num_seqs, seq_len), dtype=np.ubyte)
    for i in range(num_seqs):
        arr[i] = [base2bits[base] for base in seqs[i]]
    return arr

def arr2seq(arr: np.ndarray) -> str:
    bases_ch = [bits2base[base] for base in arr]
    return ''.join(bases_ch)

def make_array_with_all_sequences(length: int, digits: Sequence[int]) -> np.ndarray:
    num_digits = len(digits)
    num_seqs = num_digits ** length
    powers_num_digits = [num_digits ** k for k in range(length)]
    digits = np.array(digits, dtype=np.ubyte)

    arr = np.zeros((num_seqs, length), dtype=np.ubyte)
    for i, j, c in zip(reversed(powers_num_digits), powers_num_digits, list(range(length))):
        arr[:, c] = np.tile(np.repeat(digits, i), j)

    return arr

def make_array_with_all_dna_seqs(length: int, bases: Collection[str] = ('A', 'C', 'G', 'T')) -> np.ndarray:
    """Return 2D numpy array with all DNA sequences of given length in
    lexicographic order. Bases contains bases to be used: ('A','C','G','T') by
    default, but can be set to a subset of these.

    Uses the encoding described in the documentation for DNASeqList. The result
    is a 2D array, where each row represents a DNA sequence, and that row
    has one byte per base."""

    if len(bases) == 0:
        raise ValueError('bases cannot be empty')
    if not set(bases) <= {'A', 'C', 'G', 'T'}:
        raise ValueError(f"bases must be a subset of {'A', 'C', 'G', 'T'}; cannot be {bases}")

    base_bits = [base2bits[base] for base in bases]
    digits = np.array(base_bits, dtype=np.ubyte)

    return make_array_with_all_sequences(length, digits)
    # num_bases = len(bases)
    # num_seqs = num_bases ** length
    #
    # # the former code took up too much memory (using int32 or int64)
    # # the following code makes sure it's 1 byte per base
    # powers_num_bases = [num_bases ** k for k in range(length)]
    #
    # list_of_arrays = False
    # if list_of_arrays:
    #     # this one seems to be faster but takes more memory, probably because just before the last command
    #     # there are two copies of the array in memory at once
    #     columns = []
    #     for i, j, c in zip(reversed(powers_num_bases), powers_num_bases, list(range(length))):
    #         columns.append(np.tile(np.repeat(bases, i), j))
    #     arr = np.vstack(columns).transpose()
    # else:
    #     # this seems to be slightly slower but takes less memory, since it
    #     # allocates only the final array, plus one extra column of that
    #     # array at a time
    #     arr = np.empty((num_seqs, length), dtype=np.ubyte)
    #     for i, j, c in zip(reversed(powers_num_bases), powers_num_bases, list(range(length))):
    #         arr[:, c] = np.tile(np.repeat(bases, i), j)
    #
    # return arr

def make_array_with_random_subset_of_dna_seqs(
        length: int, num_random_seqs: int, rng: np.random.Generator = default_rng,
        bases: Collection[str] = ('A', 'C', 'G', 'T')) -> np.ndarray:
    """
    Return 2D numpy array with random subset of size `num_seqs` of DNA sequences of given length.
    Bases contains bases to be used: ('A','C','G','T') by default, but can be set to a subset of these.

    Uses the encoding described in the documentation for DNASeqList. The result is a 2D array,
    where each row represents a DNA sequence, and that row has one byte per base.

    Sequences returned will be unique (i.e., sampled without replacement) and in a random order

    :param length:
        length of each row
    :param num_random_seqs:
        number of rows
    :param bases:
        DNA bases to use
    :param rng:
        numpy random number generator (type returned by numpy.random.default_rng())
    :return:
        2D numpy array with random subset of size `num_seqs` of DNA sequences of given length
    """
    if length < 0:
        raise ValueError(f'length = {num_random_seqs} must be nonnegative')
    elif length == 0:
        return np.array([[]], dtype=np.ubyte)
    if num_random_seqs <= 0:
        raise ValueError(f'num_seqs = {num_random_seqs} must be positive')
    if not set(bases) <= {'A', 'C', 'G', 'T'}:
        raise ValueError(f"bases must be a subset of {'A', 'C', 'G', 'T'}; cannot be {bases}")
    if len(bases) == 0:
        raise ValueError('bases cannot be empty')
    elif len(bases) == 1:
        raise ValueError('bases must have at least two elements')

    max_possible = len(bases) ** length
    if num_random_seqs > max_possible:
        raise ValueError(f'''\
num_random_seqs = {num_random_seqs} is greater than the total number {max_possible}
of sequences of length {length} using alphabet {bases}, so we cannot guarantee
that many unique sequences. Please set num_random_seqs <= {max_possible}.''')

    # If we want sufficiently many sequences, then it's simpler to just generate all sequences
    # of that length and choose a random subset of size num_seqs.
    if num_random_seqs >= max_possible / 4:
        all_seqs = make_array_with_all_dna_seqs(length=length, bases=bases)
        # https://stackoverflow.com/a/27815343/5339430
        idxs = rng.choice(all_seqs.shape[0], num_random_seqs, replace=False)
        sampled_seqs = all_seqs[idxs]
        return sampled_seqs

    # This comment justifies why we sample 2*num_random_seqs sequences randomly (with replacement) in order
    # to get our goal of at least num_random_seqs *unique* sequences. We sample with replacement since that's
    # the only option using numpy's random number generation routines: we are generating a 2D array
    # of numbers and want each *row* to be unique, but numpy can only let us say each *number* is unique,
    # which doesn't describe what we want.
    #
    # Define
    #
    #   m = num_seqs_to_sample = number of sequences we sample with replacement
    #   n = max_possible       = total number of sequences of that length
    #   k = num_random_seqs    = desired number of unique sequences
    #
    # If we sample m sequences with replacement, out of n total, this is tossing m balls into n bins.
    # We want m sufficiently large that at least k bins are non-empty. We throw m=2*k balls.
    #
    # Above we check if k >= n/4 and generate all sequences if so. Thus, if we need the random selection
    # below, then k < n/4. Thus at least 3/4 of bins are empty the entire time, so each time we throw a ball,
    # there is a < 1/4 chance for it to land in a non-empty bin. We are throwing m=2*k balls,
    # so the number of these that land in non-empty bins is at most X=Binomial(2*k, 1/4), with E[X] = k/2.
    # Using Chernoff bounds for X=Binomial(n,p) with E[X]=np and delta=2 giving
    #   Pr[we generate fewer than k unique sequences]
    #      = Pr[more than k out of 2*k balls land in non-empty bins]
    #      = Pr[X > (1+delta)E[X]]
    #      < e^{-delta^2 E[X] / 3}
    #      = e^{-2*k/3},
    # i.e., for large k, a very small probability. Even for k=1, this is about 1/2,
    # so we'll only need to execute expected 2 iterations of the while loop below.

    base_bits = np.array([base2bits[base] for base in bases], dtype=np.ubyte)
    num_seqs_to_sample = 2 * num_random_seqs  # c*k in analysis above
    unique_sorted_arr = None

    # odds are low to have a collision, so for simplicity we just repeat the whole process if needed
    while unique_sorted_arr is None or len(unique_sorted_arr) < num_random_seqs:
        arr = rng.choice(a=base_bits, size=(num_seqs_to_sample, length))
        unique_sorted_arr = np.unique(arr, axis=0)
        if len(unique_sorted_arr) < num_random_seqs:
            print(f'WARNING: did not find {num_random_seqs} unique sequences. If you are seeing this warning '
                  f'repeatedly, check the parameters to make_array_with_random_subset_of_dna_seqs.')

    # We probably have too many, so pick a random subset of sequences to return.
    # We need a random subset, rather than just taking the first num_random_seqs elements,
    # since numpy.unique above sorts the elements.
    idxs = rng.choice(unique_sorted_arr.shape[0], num_random_seqs, replace=False)
    sampled_seqs = unique_sorted_arr[idxs]
    return sampled_seqs

# https://stackoverflow.com/questions/4941753/is-there-a-math-ncr-function-in-python
# In Python 3.8 there's math.comb, but this is about 3x faster somehow.
import operator as op
from functools import reduce

def comb(n: int, k: int) -> int:
    # n choose k = n! / (k! * (n-k)!)
    k = min(k, n - k)
    numer = reduce(op.mul, range(n, n - k, -1), 1)
    denom = reduce(op.mul, range(1, k + 1), 1)
    return numer // denom

def make_array_with_all_dna_seqs_hamming_distance(
        dist: int, seq: str, bases: Collection[str] = ('A', 'C', 'G', 'T')) -> np.ndarray:
    """
    Return 2D numpy array with all DNA sequences of given length in lexicographic order. Bases contains
    bases to be used: ('A','C','G','T') by default, but can be set to a subset of these.

    Uses the encoding described in the documentation for DNASeqList. The result is a 2D array, where each
    row represents a DNA sequence, and that row has one byte per base.
    """
    length = len(seq)
    assert 1 <= dist <= length

    num_bases = len(bases)

    if num_bases == 0:
        raise ValueError('bases cannot be empty')
    if not set(bases) <= {'A', 'C', 'G', 'T'}:
        raise ValueError(f"bases must be a subset of {'A', 'C', 'G', 'T'}; cannot be {bases}")

    num_ways_to_choose_subsequence_indices = comb(length, dist)
    num_different_bases = len(bases) - 1
    num_subsequences = num_different_bases ** dist
    num_seqs = num_ways_to_choose_subsequence_indices * num_subsequences

    # for simplicity of modular arithmetic, we use integers 0,...,len(bases)-1 to represent the bases,
    # then map these back to the correct subset of 0,1,2,3 when we are done
    offsets = range(1, num_different_bases + 1)
    subseq_offsets = make_array_with_all_sequences(length=dist, digits=offsets)
    assert len(subseq_offsets) == num_subsequences
    subseq_offsets_repeats = \
        np.tile(subseq_offsets.flatten(), num_ways_to_choose_subsequence_indices).reshape(num_seqs, dist)

    # all (length choose dist) indices where we could change the bases
    idxs = combnr_idxs(length, dist)
    assert len(idxs) == num_ways_to_choose_subsequence_indices
    idxs_repeat = np.tile(idxs, num_subsequences).reshape(num_seqs, length)
    assert len(idxs_repeat) == num_seqs

    # map subset of bases used to *prefix* of 0,1,2,3
    base2bits_local = {base: digit for base, digit in zip(bases, range(4))}
    seq_as_arr = seq2arr(seq, base2bits_local=base2bits_local)
    new_arr = np.tile(seq_as_arr, num_seqs).reshape(num_seqs, length)

    new_arr[idxs_repeat] += subseq_offsets_repeats.flatten()
    new_arr %= num_bases

    # now map back to correct subset of 0,1,2,3 to represent bases
    for base, digit in zip(['A', 'C', 'G', 'T'], range(4)):
        if base not in bases:
            idxs_to_inc = new_arr >= digit
            new_arr[idxs_to_inc] += 1

    return new_arr

def combnr_idxs(length: int, number: int) -> np.ndarray:
    # Gives 2D Boolean numpy array, with `length` columns and (`length` choose `number`) rows,
    # representing all ways to set exactly `number` elements of the row True and the others to False.
    # Useful for indexing into a same-shape numpy array, changing exactly `number` elements in each row.
    #
    # :param length:
    #     number of columns
    # :param number:
    #     number of True values in each row
    # :return:
    #     numpy array, with `length` columns and (`length` choose `number`) rows,
    #     representing all ways to set exactly `number` elements of the row True and the others to False.
    if number == 0:
        return np.array([[False] * length], dtype=bool)
    elif number > length // 2:
        vals = combnr_idxs(length, length - number)
        return np.logical_not(vals)
    x = np.array(np.meshgrid(*([np.arange(0, length)] * number))).T.reshape(-1, number)
    z = np.sum(np.identity(length)[x], 1, dtype=bool).astype(int)
    return np.unique(z[np.sum(z, axis=1) == number], axis=0).astype(bool)

def make_array_with_random_subset_of_dna_seqs_hamming_distance(
        num_seqs: int, dist: int, seq: str, rng: np.random.Generator = default_rng,
        bases: Collection[str] = ('A', 'C', 'G', 'T')) -> np.ndarray:
    """
    Return 2D numpy array with random subset of size `num_seqs` of DNA sequences of given length.
    Bases contains bases to be used: ('A','C','G','T') by default, but can be set to a subset of these.

    Uses the encoding described in the documentation for DNASeqList. The result is a 2D array,
    where each row represents a DNA sequence, and that row has one byte per base.

    Sampled *with* replacement, so the same row may appear twice in the returned array

    :param num_seqs:
        number of sequences to generate
    :param dist:
        Hamming distance to be from `seq`
    :param seq:
        sequence to generate other sequences close to
    :param bases:
        DNA bases to use
    :param rng:
        numpy random number generator (type returned by numpy.random.default_rng())
    :return:
        2D numpy array with random subset of size `num_seqs` of DNA sequences of given length
    """
    if not set(bases) <= {'A', 'C', 'G', 'T'}:
        raise ValueError(f"bases must be a subset of {'A', 'C', 'G', 'T'}; cannot be {bases}")
    if len(bases) == 0:
        raise ValueError('bases cannot be empty')
    elif len(bases) == 1:
        raise ValueError('bases must have at least two elements')

    length = len(seq)

    if dist < 1:
        raise ValueError(f'dist must be positive, but dist = {dist}')
    if dist > length:
        raise ValueError(f'should have dist <= len("{seq}") = {length}, but dist = {dist}')

    num_different_bases = len(bases) - 1

    # map subset of bases used to *prefix* of 0,1,2,3
    base2bits_local = {base: digit for base, digit in zip(bases, range(4))}
    seq_as_arr = seq2arr(seq, base2bits_local=base2bits_local)
    seqs = np.tile(seq_as_arr, num_seqs).reshape(num_seqs, length)

    # for simplicity of modular arithmetic, we use integers 0,...,len(bases)-1 to represent the bases,
    # then map these back to the correct subset of 0,1,2,3 when we are done
    subseq_offsets = rng.integers(low=1, high=num_different_bases + 1, size=(num_seqs, dist))
    assert len(subseq_offsets) == num_seqs

    # all (length choose dist) indices where we could change the bases
    idxs = random_choice_noreplace(np.arange(length), dist, num_seqs, rng)
    assert len(idxs) == num_seqs
    changes = rng.integers(1, num_different_bases + 1, size=idxs.shape, dtype=np.uint8)
    seqs[np.arange(num_seqs)[:, None], idxs] += changes
    seqs = np.mod(seqs, num_different_bases + 1)

    # now map back to correct subset of 0,1,2,3 to represent bases
    for base, digit in zip(['A', 'C', 'G', 'T'], range(4)):
        if base not in bases:
            idxs_to_inc = seqs >= digit
            seqs[idxs_to_inc] += 1

    # use the next two lines to return only unique rows
    # seqs = np.unique(seqs, axis=0)
    # rng.shuffle(seqs)

    return seqs

# def random_hamming(sequence: Union[List[int], np.ndarray], distance: int, number: int,
#                    rng: np.random.Generator) -> np.ndarray:
#     sequence = np.array(sequence, dtype=np.uint8)
#     length = len(sequence)
#     seqrepeats = np.tile(sequence, number).reshape((number, length))
#     places = random_choice_noreplace(np.arange(length), distance, number, rng)
#     changes = rng.integers(1, 4, size=places.shape, dtype=np.uint8)
#     seqrepeats[np.arange(number)[:, None], places] += changes
#     seqrepeats = np.mod(seqrepeats, 4)
#     return seqrepeats

# https://stackoverflow.com/a/59328647/5339430
def random_choice_noreplace(l: np.ndarray, n_sample: int, num_draw: int,
                            rng: np.random.Generator) -> np.ndarray:
    '''
    l: 1-D array or list
    n_sample: sample size for each draw
    num_draw: number of draws

    Intuition: Randomly generate numbers, get the index of the smallest n_sample number for each row.
    '''
    l = np.array(l)
    random_array_floats = rng.random((num_draw, len(l)))
    return l[np.argpartition(random_array_floats, n_sample - 1, axis=-1)[:, :n_sample]]

# @lru_cache(maxsize=10000000)
def longest_common_substring(a1: np.ndarray, a2: np.ndarray, vectorized: bool = True) -> Tuple[int, int, int]:
    """Return start and end indices (a1start, a2start, length) of longest common
    substring (subarray) of 1D arrays a1 and a2."""
    assert len(a1.shape) == 1
    assert len(a2.shape) == 1
    counter = np.zeros(shape=(len(a1) + 1, len(a2) + 1), dtype=np.int)
    a1idx_longest = a2idx_longest = -1
    len_longest = 0

    if vectorized:
        for i1 in range(len(a1)):
            idx = (a2 == a1[i1])
            idx_shifted = np.hstack([[False], idx])
            counter[i1 + 1, idx_shifted] = counter[i1, idx] + 1
        idx_longest = np.unravel_index(np.argmax(counter), counter.shape)
        if idx_longest[0] > 0:
            len_longest = counter[idx_longest]
            a1idx_longest = int(idx_longest[0] - len_longest)
            a2idx_longest = int(idx_longest[1] - len_longest)
    else:
        for i1 in range(len(a1)):
            for i2 in range(len(a2)):
                if a1[i1] == a2[i2]:
                    c = counter[i1, i2] + 1
                    counter[i1 + 1, i2 + 1] = c
                    if c > len_longest:
                        len_longest = c
                        a1idx_longest = i1 + 1 - c
                        a2idx_longest = i2 + 1 - c
    return a1idx_longest, a2idx_longest, len_longest

# @lru_cache(maxsize=10000000)
def longest_common_substrings_singlea1(a1: np.ndarray, a2s: np.ndarray) \
        -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Return start and end indices (a1starts, a2starts, lengths) of longest common
    substring (subarray) of 1D array a1 and rows of 2D array a2s.

    If length[i]=0, then a1starts[i]=a2starts[i]=0 (not -1), so be sure to check
    length[i] to see if any substrings actually matched."""
    assert len(a1.shape) == 1
    assert len(a2s.shape) == 2
    numa2s = a2s.shape[0]
    len_a1 = len(a1)
    len_a2 = a2s.shape[1]
    counter = np.zeros(shape=(len_a1 + 1, numa2s, len_a2 + 1), dtype=np.int)

    for i1 in range(len(a1)):
        idx = (a2s == a1[i1])
        idx_shifted = np.insert(idx, 0, np.zeros(numa2s, dtype=np.bool), axis=1)
        counter[i1 + 1, idx_shifted] = counter[i1, idx] + 1

    counter = np.swapaxes(counter, 0, 1)

    counter_flat = counter.reshape(numa2s, (len_a1 + 1) * (len_a2 + 1))
    idx_longest_raveled = np.argmax(counter_flat, axis=1)
    len_longest = counter_flat[np.arange(counter_flat.shape[0]), idx_longest_raveled]

    idx_longest = np.unravel_index(idx_longest_raveled, dims=(len_a1 + 1, len_a2 + 1))
    a1idx_longest = idx_longest[0] - len_longest
    a2idx_longest = idx_longest[1] - len_longest

    return a1idx_longest, a2idx_longest, len_longest

# @lru_cache(maxsize=10000000)
def longest_common_substrings_product(a1s: np.ndarray, a2s: np.ndarray) \
        -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Return start and end indices (a1starts, a2starts, lengths) of longest common
    substring (subarray) of each pair in the cross product of rows of a1s and a2s.

    If length[i]=0, then a1starts[i]=a2starts[i]=0 (not -1), so be sure to check
    length[i] to see if any substrings actually matched."""
    numa1s = a1s.shape[0]
    numa2s = a2s.shape[0]
    a1s_cp = np.repeat(a1s, numa2s, axis=0)
    a2s_cp = np.tile(a2s, (numa1s, 1))

    a1idx_longest, a2idx_longest, len_longest = _longest_common_substrings_pairs(a1s_cp, a2s_cp)

    a1idx_longest = a1idx_longest.reshape(numa1s, numa2s)
    a2idx_longest = a2idx_longest.reshape(numa1s, numa2s)
    len_longest = len_longest.reshape(numa1s, numa2s)

    return a1idx_longest, a2idx_longest, len_longest

def pair_index(n: int) -> np.ndarray:
    index = np.fromiter(it.chain.from_iterable(it.combinations(range(n), 2)), int, count=n * (n - 1))
    return index.reshape(-1, 2)

def _longest_common_substrings_pairs(a1s: np.ndarray, a2s: np.ndarray) \
        -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    assert len(a1s.shape) == 2
    assert len(a2s.shape) == 2
    assert a1s.shape[0] == a2s.shape[0]

    numpairs = a1s.shape[0]

    len_a1 = a1s.shape[1]
    len_a2 = a2s.shape[1]

    counter = np.zeros(shape=(len_a1 + 1, numpairs, len_a2 + 1), dtype=np.int)

    for i1 in range(len_a1):
        a1s_cp_col = a1s[:, i1].reshape(numpairs, 1)
        a1s_cp_col_rp = np.repeat(a1s_cp_col, len_a2, axis=1)

        idx = (a2s == a1s_cp_col_rp)
        idx_shifted = np.hstack([np.zeros(shape=(numpairs, 1), dtype=np.bool), idx])
        counter[i1 + 1, idx_shifted] = counter[i1, idx] + 1

    counter = np.swapaxes(counter, 0, 1)

    counter_flat = counter.reshape(numpairs, (len_a1 + 1) * (len_a2 + 1))
    idx_longest_raveled = np.argmax(counter_flat, axis=1)
    len_longest = counter_flat[np.arange(counter_flat.shape[0]), idx_longest_raveled]

    idx_longest = np.unravel_index(idx_longest_raveled, dims=(len_a1 + 1, len_a2 + 1))
    a1idx_longest = idx_longest[0] - len_longest
    a2idx_longest = idx_longest[1] - len_longest

    return a1idx_longest, a2idx_longest, len_longest

def longest_common_substrings_all_pairs_strings(seqs1: Sequence[str], seqs2: Sequence[str]) \
        -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """For Python strings"""
    a1s = seqs2arr(seqs1)
    a2s = seqs2arr(seqs2)
    return _longest_common_substrings_pairs(a1s, a2s)

def _strongest_common_substrings_all_pairs_return_energies_and_counter(
        a1s: np.ndarray, a2s: np.ndarray, temperature: float) \
        -> Tuple[np.ndarray, np.ndarray]:
    assert len(a1s.shape) == 2
    assert len(a2s.shape) == 2
    assert a1s.shape[0] == a2s.shape[0]

    numpairs = a1s.shape[0]
    len_a1 = a1s.shape[1]
    len_a2 = a2s.shape[1]
    counter = np.zeros(shape=(len_a1 + 1, numpairs, len_a2 + 1), dtype=np.int)
    energies = np.zeros(shape=(len_a1 + 1, numpairs, len_a2 + 1), dtype=np.float)

    #     if not loop_energies:
    loop_energies = calculate_loop_energies(temperature)

    prev_match_shifted_idxs = None

    for i1 in range(len_a1):
        a1s_col = a1s[:, i1].reshape(numpairs, 1)
        a1s_col_rp = np.repeat(a1s_col, len_a2, axis=1)

        # find matching chars and extend length of substring
        match_idxs = (a2s == a1s_col_rp)
        match_shifted_idxs = np.hstack([np.zeros(shape=(numpairs, 1), dtype=np.bool), match_idxs])
        counter[i1 + 1, match_shifted_idxs] = counter[i1, match_idxs] + 1

        if i1 > 0:
            # calculate energy if matching substring has length > 1
            prev_bases = a1s[:, i1 - 1]
            cur_bases = a1s[:, i1]
            loops = (prev_bases << 2) + cur_bases
            latest_energies = loop_energies[loops].reshape(numpairs, 1)
            latest_energies_rp = np.repeat(latest_energies, len_a2, axis=1)
            match_idxs_false_at_end = np.hstack([match_idxs, np.zeros(shape=(numpairs, 1), dtype=np.bool)])
            both_match_idxs = match_idxs_false_at_end & prev_match_shifted_idxs
            prev_match_shifted_shifted_idxs = np.hstack(
                [np.zeros(shape=(numpairs, 1), dtype=np.bool), prev_match_shifted_idxs])[:, :-1]
            both_match_shifted_idxs = match_shifted_idxs & prev_match_shifted_shifted_idxs
            energies[i1 + 1, both_match_shifted_idxs] = energies[i1, both_match_idxs] + latest_energies_rp[
                both_match_idxs]

        #         prev_match_idxs = match_idxs
        prev_match_shifted_idxs = match_shifted_idxs

    counter = counter.swapaxes(0, 1)
    energies = energies.swapaxes(0, 1)

    return counter, energies

def internal_loop_penalty(n: int, temperature: float) -> float:
    return 1.5 + (2.5 * 0.002 * temperature * math.log(1 + n))

def _strongest_common_substrings_all_pairs(a1s: np.ndarray, a2s: np.ndarray, temperature: float) \
        -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    numpairs = a1s.shape[0]
    len_a1 = a1s.shape[1]
    len_a2 = a2s.shape[1]

    counter, energies = _strongest_common_substrings_all_pairs_return_energies_and_counter(a1s, a2s,
                                                                                temperature)

    counter_flat = counter.reshape(numpairs, (len_a1 + 1) * (len_a2 + 1))
    energies_flat = energies.reshape(numpairs, (len_a1 + 1) * (len_a2 + 1))

    idx_strongest_raveled = np.argmax(energies_flat, axis=1)
    len_strongest = counter_flat[np.arange(counter_flat.shape[0]), idx_strongest_raveled]
    energy_strongest = energies_flat[np.arange(counter_flat.shape[0]), idx_strongest_raveled]

    idx_strongest = np.unravel_index(idx_strongest_raveled, dims=(len_a1 + 1, len_a2 + 1))
    a1idx_strongest = idx_strongest[0] - len_strongest
    a2idx_strongest = idx_strongest[1] - len_strongest

    return a1idx_strongest, a2idx_strongest, len_strongest, energy_strongest

def strongest_common_substrings_all_pairs_string(seqs1: Sequence[str], seqs2: Sequence[str],
                                                 temperature: float) \
        -> Tuple[List[float], List[float], List[float], List[float]]:
    """For Python strings representing DNA; checks for reverse complement matches
    rather than direct matches, and evaluates nearest neighbor energy, returning
    indices lengths, and energies of strongest complementary substrings."""
    a1s = seqs2arr(seqs1)
    a2s = seqs2arr(seqs2)
    a1idx_strongest, a2idx_strongest, len_strongest, energy_strongest = _strongest_common_substrings_all_pairs(
        a1s,
        wc_arr(a2s),
        temperature)
    return list(a1idx_strongest), list(a2idx_strongest), list(len_strongest), list(energy_strongest)

def energies_strongest_common_substrings(seqs1: Sequence[str], seqs2: Sequence[str], temperature: float) \
        -> List[float]:
    a1s = seqs2arr(seqs1)
    a2s = seqs2arr(seqs2)
    a1idx_strongest, a2idx_strongest, len_strongest, energy_strongest = \
        _strongest_common_substrings_all_pairs(a1s, wc_arr(a2s), temperature)
    return list(energy_strongest)

@dataclass
class DNASeqList:
    """
    Represents a list of DNA sequences of identical length. The sequences are stored as a 2D numpy array
    of bytes :py:data:`DNASeqList.seqarr`. Each byte represents a single DNA base (so it is not a compact
    representation; the most significant 6 bits of the byte will always be 0).
    """

    seqarr: np.ndarray
    """
    Uses a (noncompact) internal representation using 8 bits (1 byte, dtype = np.ubyte) per base,
    stored in a numpy 2D array of bytes.
    Each row (axis 0) is a DNA sequence, and each column (axis 1) is a base in a sequence.

    The code used is :math:`A \\to 0, C \\to 1, G \\to 2, T \\to 3`.
    """

    numseqs: int
    """Number of DNA sequences (number of rows, axis 0, in :py:data:`DNASeqList.seqarr`)"""

    seqlen: int
    """Length of each DNA sequence (number of columns, axis 1, in :py:data:`DNASeqList.seqarr`)"""

    rng: np.random.Generator
    """Random number generator to use."""

    def __init__(self,
                 length: Optional[int] = None,
                 num_random_seqs: Optional[int] = None,
                 shuffle: bool = False,
                 alphabet: Collection[str] = ('A', 'C', 'G', 'T'),
                 seqs: Optional[Sequence[str]] = None,
                 seqarr: np.ndarray = None,
                 filename: Optional[str] = None,
                 rng: np.random.Generator = default_rng,
                 hamming_distance_from_sequence: Optional[Tuple[int, str]] = None):
        """
        Creates a set of DNA sequences, all of the same length.

        Create either all sequences of a given length if seqs is not specified,
        or all sequences in seqs if seqs is specified. If neither is specified
        then all sequences of length 3 are created.

        *Exactly one* of the following should be specified:

        - `length` (possibly along with `alphabet` and `num_random_seqs`)

        - `seqs`

        - `seqarr`

        - `filename`

        - `hamming_distance_from_sequence` (possibly along with `alphabet` and `num_random_seqs`)

        :param length:
            length of sequences; `num_seqs` and `alphabet` can also be specified along with it
        :param hamming_distance_from_sequence:
            if specified and equal to `(dist, seq)` of type (int, str),
            then only sequences at Hamming distance `dist` from `seq` will be generated.
            Raises error if `length`, `seqs`, `seqarr`, or `filename` is specified.
        :param num_random_seqs:
            number of sequences to generate; if not specified, then all sequences
            of length `length` using bases from `alphabet` are generated.
            Sequences are sampled *with* replacement, so the same sequence may appear twice.
        :param shuffle:
            whether to shuffle sequences
        :param alphabet:
            a subset of {'A', 'C', 'G', 'T'}
        :param seqs:
            sequence (e.g., list or tuple) of strings, all of the same length
        :param seqarr:
            2D NumPy array, with axis 0 moving between sequences,
            and axis 1 moving between consecutive DNA bases in a sequence
        :param filename:
            name of file containing a :any:`DNASeqList`
            as written by :py:meth:`DNASeqList.write_to_file`
        :param rng:
            numpy random number generator (type returned by numpy.random.default_rng())
        """
        for v1, v2 in it.combinations([length, seqs, seqarr, filename, hamming_distance_from_sequence], 2):
            if v1 is not None and v2 is not None:
                raise ValueError('exactly one of length, seqs, seqarr, filename, or '
                                 'hamming_distance_from_sequence must be non-None')
        self.rng = rng
        if seqarr is not None:
            self.seqarr = seqarr
            self._update_size()

        elif seqs is not None:
            if len(seqs) == 0:
                raise ValueError('seqs must have positive length')
            self.seqlen = len(seqs[0])
            for seq in seqs:
                if len(seq) != self.seqlen:
                    raise ValueError('All sequences in seqs must be equal length')
            self.numseqs = len(seqs)
            self.seqarr = seqs2arr(seqs)

        elif filename is not None:
            self._read_from_file(filename)

        elif length is not None:
            if num_random_seqs is None:
                self.seqarr = make_array_with_all_dna_seqs(length=length, bases=alphabet)
            else:
                self.seqarr = make_array_with_random_subset_of_dna_seqs(
                    length=length, num_random_seqs=num_random_seqs, rng=self.rng, bases=alphabet)
            self.seqlen = length
            self.numseqs = len(self.seqarr) if self.seqlen > 0 else 1

        elif hamming_distance_from_sequence is not None:
            dist, seq = hamming_distance_from_sequence
            if num_random_seqs is None:
                self.seqarr = make_array_with_all_dna_seqs_hamming_distance(dist=dist, seq=seq,
                                                                            bases=alphabet)
            else:
                self.seqarr = make_array_with_random_subset_of_dna_seqs_hamming_distance(
                    num_seqs=num_random_seqs, dist=dist, seq=seq, rng=self.rng, bases=alphabet)
            self.seqlen = len(seq)
            self.numseqs = len(self.seqarr) if self.seqlen > 0 else 1

        else:
            raise ValueError('at least one of length, seqs, seqarr, filename, or '
                             'hamming_distance_from_sequence must be specified')

        self.shift = np.arange(2 * (self.seqlen - 1), -1, -2)

        if shuffle:
            self.shuffle()

    def random_choice(self, num: int, rng: np.random.Generator = default_rng,
                      replace: bool = False) -> List[str]:
        """
        Returns random choice of `num` DNA sequence(s) (represented as list of Python strings).

        :param num:
            number of sequences to sample
        :param replace:
            whether to sample with replacement
        :return:
            sampled sequences
        """
        idxs = rng.choice(np.arange(self.numseqs), num, replace=replace)
        seqs = [self[int(idx)] for idx in idxs]
        return seqs

    def random_sequence(self, rng: np.random.Generator = default_rng) -> str:
        """
        Returns random DNA sequence (represented as Python string).

        :return:
            sampled sequence
        """
        idx = int(rng.integers(0, self.numseqs))
        return self[idx]

    def _update_size(self) -> None:
        # updates numseqs and seqlen based on shape of seqarr
        self.numseqs, self.seqlen = self.seqarr.shape

    def __len__(self) -> int:
        return self.numseqs

    def __contains__(self, seq: str) -> bool:
        if len(seq) != self.seqlen:
            return False
        arr = seq2arr(seq)
        return np.any(~np.any(self.seqarr - arr, 1))

    def _read_from_file(self, filename: str) -> None:
        """Reads from fileName in the format defined in writeToFile.
        Only meant to be called from constructor."""
        with open(filename, 'r+') as f:
            first_line = f.readline()
            num_seqs_str, seq_len_str, temperature = first_line.split()
            self.numseqs = int(num_seqs_str)
            self.seqlen = int(seq_len_str)
            self.seqarr = np.empty((self.numseqs, self.seqlen), dtype=np.ubyte)
            for i in range(self.numseqs):
                line = f.readline()
                seq = line.strip()
                self.seqarr[i] = [base2bits[base] for base in seq]

    def write_to_file(self, filename: str) -> None:
        """Writes text file describing DNA sequence list, in format

        numseqs seqlen
        seq1
        seq2
        seq3
        ...

        where numseqs, seqlen are integers, and seq1,
        ... are strings from {A,C,G,T}"""
        with open(filename, 'w+') as f:
            f.write(str(self.numseqs) + ' ' + str(self.seqlen) + '\n')
            for i in range(self.numseqs):
                f.write(self.get_seq_str(i) + '\n')

    def wcenergy(self, idx: int, temperature: float) -> float:
        """Return energy of idx'th sequence binding to its complement."""
        return wcenergy(self.seqarr[idx], temperature)

    def __repr__(self) -> str:
        return 'DNASeqSet(seqs={})'.format(str([self[i] for i in range(self.numseqs)]))

    def __str__(self) -> str:
        if self.numseqs <= 256:
            ret = [self.get_seq_str(i) for i in range(self.numseqs)]
            return ','.join(ret)
        else:
            ret = [self.get_seq_str(i) for i in range(3)] + ['...'] + \
                  [self.get_seq_str(i) for i in range(self.numseqs - 3, self.numseqs)]
            return ','.join(ret)

    def shuffle(self) -> None:
        self.rng.shuffle(self.seqarr)

    def to_list(self) -> List[str]:
        """Return list of strings representing the sequences, e.g. ['ACG','TAA']"""
        return [self.get_seq_str(idx) for idx in range(self.numseqs)]

    def get_seq_str(self, idx: int) -> str:
        """Return idx'th DNA sequence as a string."""
        return arr2seq(self.seqarr[idx])

    def get_seqs_str_list(self, slice_: slice) -> List[str]:
        """Return a list of strings specified by slice."""
        bases_lst = self.seqarr[slice_]
        ret = []
        for bases in bases_lst:
            bases_ch = [bits2base[base] for base in bases]
            ret.append(''.join(bases_ch))
        return ret

    def keep_seqs_at_indices(self, indices: Iterable[int]) -> None:
        """Keeps only sequences at the given indices."""
        if not isinstance(indices, list):
            indices = list(indices)
        self.seqarr = self.seqarr[indices]
        self._update_size()

    def __getitem__(self, slice_: Union[int, slice]) -> Union[str, List[str]]:
        if isinstance(slice_, int):
            return self.get_seq_str(slice_)
        elif isinstance(slice_, slice):
            return self.get_seqs_str_list(slice_)
        else:
            raise ValueError('slice_ must be int or slice')

    def __setitem__(self, idx: int, seq: str) -> None:
        # cannot set on slice
        self.seqarr[idx] = seq2arr(seq)

    def pop(self) -> str:
        """Remove and return last seq, as a string."""
        seq_str = self.get_seq_str(-1)
        self.seqarr = np.delete(self.seqarr, -1, 0)
        self.numseqs -= 1
        return seq_str

    def pop_array(self) -> np.ndarray:
        """Remove and return last seq, as a string."""
        arr = self.seqarr[-1]
        self.seqarr = np.delete(self.seqarr, -1, 0)
        self.numseqs -= 1
        return arr

    def append_seq(self, newseq: str) -> None:
        self.append_arr(seq2arr(newseq))

    def append_arr(self, newarr: np.ndarray) -> None:
        self.seqarr = np.vstack([self.seqarr, newarr])
        self.numseqs += 1

    def sequences_at_hamming_distance(self, sequence: str, distance: int) -> DNASeqList:
        sequence_1d_array = seq2arr(sequence)
        distances = np.sum(np.bitwise_xor(self.seqarr, sequence_1d_array) != 0, axis=1)
        indices_at_distance = distances == distance
        arr = self.seqarr[indices_at_distance]
        return DNASeqList(seqarr=arr)

    def hamming_map(self, sequence: str) -> Dict[int, DNASeqList]:
        """Return dict mapping each length `d` to a :any:`DNASeqList` of sequences that are
        Hamming distance `d` from `seq`."""
        # import time
        # times = []
        # before = time.perf_counter_ns()
        sequence_1d_array = seq2arr(sequence)
        distances = np.sum(np.bitwise_xor(self.seqarr, sequence_1d_array) != 0, axis=1)
        distance_map = {}
        # after_it_prev = time.perf_counter_ns()
        # print(f'time to calculate Hamming distances: {(after_it_prev - before) / 1e6:.1f} ms')
        for distance in range(self.seqlen + 1):
            indices_at_distance = distances == distance
            arr = self.seqarr[indices_at_distance]
            if arr.shape[0] > 0:  # don't bother putting empty array into map
                distance_map[distance] = DNASeqList(seqarr=arr)
            # after_it_next = time.perf_counter_ns()
            # times.append((after_it_next - after_it_prev) / 1e6)
            # after_it_prev = after_it_next
        # after = time.perf_counter_ns()
        # print(f'time spent finding neighbors: {(after - before) / 1e6:.1f} ms')
        # print(f'times in each iteration: {times}')
        return distance_map

    def sublist(self, start: int, end: Optional[int] = None) -> DNASeqList:
        """Return sublist of DNASeqList from `start`, inclusive, to `end`, exclusive.

        If `end` is not specified, goes until the end of the list."""
        if end is None:
            end = self.numseqs
        arr = self.seqarr[start:end]
        return DNASeqList(seqarr=arr)

    # def filter_hamming(self, threshold: int) -> None:
    #     seq = self.pop_array()
    #     arr_keep = np.array([seq])
    #     self.shuffle()
    #     while self.seqarr.shape[0] > 0:
    #         seq = self.pop_array()
    #         while self.seqarr.shape[0] > 0:
    #             hamming_min = np.min(np.sum(np.bitwise_xor(arr_keep, seq) != 0, axis=1))
    #             too_close = (hamming_min < threshold)
    #             if not too_close:
    #                 break
    #             seq = self.pop_array()
    #         arr_keep = np.vstack([arr_keep, seq])
    #     self.seqarr = arr_keep
    #     self.numseqs = self.seqarr.shape[0]
    #
    # def hamming_min(self, arr: np.ndarray) -> int:
    #     """Returns minimum Hamming distance between arr and any sequence in
    #     this DNASeqList."""
    #     distances = np.sum(np.bitwise_xor(self.seqarr, arr) != 0, axis=1)
    #     return np.min(distances)

    def filter_energy(self, low: float, high: float, temperature: float) -> DNASeqList:
        """Return new DNASeqList with seqs whose wc complement energy is within
        the given range."""
        wcenergies = calculate_wc_energies(self.seqarr, temperature)
        within_range = (low <= wcenergies) & (wcenergies <= high)
        new_seqarr = self.seqarr[within_range]
        return DNASeqList(seqarr=new_seqarr)

    def energies(self, temperature: float) -> np.ndarray:
        """
        :param temperature:
            temperature in Celsius
        :return:
            nearest-neighbor energies of each sequence with its perfect Watson-Crick complement
        """
        wcenergies = calculate_wc_energies(self.seqarr, temperature)
        return wcenergies

    def filter_end_gc(self) -> DNASeqList:
        """Remove any sequence with A or T on the end. Also remove domains that
        do not have an A or T either next to that base, or one away. Otherwise
        we could get a domain ending in {C,G}^3, which, placed next to any
        domain ending in C or G, will create a substring in {C,G}^4 and be
        rejected if we are filtering those."""
        left = self.seqarr[:, 0]
        right = self.seqarr[:, -1]
        left_p1 = self.seqarr[:, 1]
        left_p2 = self.seqarr[:, 2]
        right_m1 = self.seqarr[:, -2]
        right_m2 = self.seqarr[:, -3]
        abits = base2bits['A']
        cbits = base2bits['C']
        gbits = base2bits['G']
        tbits = base2bits['T']
        good = (((left == cbits) | (left == gbits)) & ((right == cbits) | (right == gbits)) &
                ((left_p1 == abits) | (left_p1 == tbits) | (left_p2 == abits) | (left_p2 == tbits)) &
                ((right_m1 == abits) | (right_m1 == tbits) | (right_m2 == abits) | (right_m2 == tbits)))
        seqarrpass = self.seqarr[good]
        return DNASeqList(seqarr=seqarrpass)

    def filter_end_at(self, gc_near_end: bool = False) -> DNASeqList:
        """Remove any sequence with C or G on the end. Also, if gc_near_end is True,
        remove domains that do not have an C or G either next to that base,
        or one away, to prevent breathing."""
        left = self.seqarr[:, 0]
        right = self.seqarr[:, -1]
        abits = base2bits['A']
        tbits = base2bits['T']
        good = ((left == abits) | (left == tbits)) & ((right == abits) | (right == tbits))
        if gc_near_end:
            cbits = base2bits['C']
            gbits = base2bits['G']
            left_p1 = self.seqarr[:, 1]
            left_p2 = self.seqarr[:, 2]
            right_m1 = self.seqarr[:, -2]
            right_m2 = self.seqarr[:, -3]
            good = (good &
                    ((left_p1 == cbits) | (left_p1 == gbits) | (left_p2 == cbits) | (left_p2 == gbits)) &
                    ((right_m1 == cbits) | (right_m1 == gbits) | (right_m2 == cbits) | (right_m2 == gbits)))
        seqarrpass = self.seqarr[good]
        return DNASeqList(seqarr=seqarrpass)

    def filter_base_nowhere(self, base: str) -> DNASeqList:
        """Remove any sequence that has given base anywhere."""
        good = (self.seqarr != base2bits[base]).all(axis=1)
        seqarrpass = self.seqarr[good]
        return DNASeqList(seqarr=seqarrpass)

    def filter_base_count(self, base: str, low: int, high: int) -> DNASeqList:
        """Remove any sequence not satisfying low <= #base <= high."""
        sumarr = np.sum(self.seqarr == base2bits[base], axis=1)
        good = (low <= sumarr) & (sumarr <= high)
        seqarrpass = self.seqarr[good]
        return DNASeqList(seqarr=seqarrpass)

    def filter_base_at_pos(self, pos: int, base: str) -> DNASeqList:
        """Remove any sequence that does not have given base at position pos."""
        mid = self.seqarr[:, pos]
        good = (mid == base2bits[base])
        seqarrpass = self.seqarr[good]
        return DNASeqList(seqarr=seqarrpass)

    def filter_substring(self, subs: Sequence[str]) -> DNASeqList:
        """Remove any sequence with any elements from subs as a substring."""
        if len(set([len(sub) for sub in subs])) != 1:
            raise ValueError('All substrings in subs must be equal length: %s' % subs)
        sublen = len(subs[0])
        subints = [[base2bits[base] for base in sub] for sub in subs]
        powarr = [4 ** k for k in range(sublen)]
        subvals = np.dot(subints, powarr)
        toeplitz = create_toeplitz(self.seqlen, sublen)
        convolution = np.dot(toeplitz, self.seqarr.transpose())
        passall = np.ones(self.numseqs, dtype=np.bool)
        for subval in subvals:
            passsub = np.all(convolution != subval, axis=0)
            passall = passall & passsub
        seqarrpass = self.seqarr[passall]
        return DNASeqList(seqarr=seqarrpass)

    def filter_seqs_by_g_quad(self) -> DNASeqList:
        """Removes any sticky ends with 4 G's in a row (a G-quadruplex)."""
        return self.filter_substring(['GGGG'])

    def filter_seqs_by_g_quad_c_quad(self) -> DNASeqList:
        """Removes any sticky ends with 4 G's or C's in a row (a quadruplex)."""
        return self.filter_substring(['GGGG', 'CCCC'])

    def index(self, sequence: Union[str, np.ndarray]) -> int:
        # finds index of sequence in (rows of) self.seqarr
        # raises IndexError if not present
        # taken from https://stackoverflow.com/questions/40382384/finding-a-matching-row-in-a-numpy-matrix
        if isinstance(sequence, str):
            sequence = seq2arr(sequence)
        matching_condition = (self.seqarr == sequence).all(axis=1)
        all_indices_tuple = np.where(matching_condition)
        all_indices = all_indices_tuple[0]
        first_index = all_indices[0]
        return int(first_index)

def create_toeplitz(seqlen: int, sublen: int, indices: Optional[Sequence[int]] = None) -> np.ndarray:
    """Creates a toeplitz matrix, useful for finding subsequences.

    `seqlen` is length of larger sequence; `sublen` is length of substring we're checking for.
    If `indices` is None, then all rows are created, otherwise only rows for checking those indices
    are created."""
    powarr = [4 ** k for k in range(sublen)]
    if indices is None:
        rows = list(range(seqlen - (sublen - 1)))
    else:
        rows = sorted(list(set(indices)))
        for idx in rows:
            if idx < 0:
                raise ValueError(f'index must be nonnegative, but {idx} is not; all indices = {indices}')
            if idx >= seqlen - (sublen - 1):
                raise ValueError(f'index must be less than {seqlen - (sublen - 1)}, '
                                 f'but {idx} is not; all indices = {indices}')
    num_rows = len(rows)
    num_cols = seqlen
    toeplitz = np.zeros((num_rows, num_cols), dtype=np.int)
    toeplitz[:, 0:sublen] = [powarr] * num_rows
    shift = list(rows)
    for i in range(len(rows)):
        toeplitz[i] = np.roll(toeplitz[i], shift[i])
    return toeplitz

@lru_cache(maxsize=32)
def calculate_loop_energies(temperature: float, negate: bool = False) -> np.ndarray:
    """Get SantaLucia and Hicks nearest-neighbor loop energies for given temperature,
    1 M Na+. """
    energies = (_dH - (temperature + 273.15) * _dS / 1000.0)
    if negate:
        energies = -energies
    return energies
    # SantaLucia & Hicks' values are in cal/mol/K for dS, and kcal/mol for dH.
    # Here we divide dS by 1000 to get the RHS term into units of kcal/mol/K
    # which gives an overall dG in units of kcal/mol.
    # One reason we want dG to be in units of kcal/mol is to
    # give reasonable/readable numbers close to 0 for dG(Assembly).
    # The reason we might want to flip the sign is that, by convention, in the kTAM, G_se
    # (which is computed from the usually negative dG here) is usually positive.

# _dH and _dS come from Table 1 in SantaLucia and Hicks, Annu Rev Biophys Biomol Struct. 2004;33:415-40.
#                 AA    AC    AG    AT    CA    CC     CG    CT
_dH = np.array([-7.6, -8.4, -7.8, -7.2, -8.5, -8.0, -10.6, -7.8,
                # GA    GC    GG    GT    TA    TC    TG    TT
                -8.2, -9.8, -8.0, -8.4, -7.2, -8.2, -8.5, -7.6],
               dtype=np.float32)

#                  AA     AC     AG     AT     CA     CC     CG     CT
_dS = np.array([-21.3, -22.4, -21.0, -20.4, -22.7, -19.9, -27.2, -21.0,
                #  GA     GC     GG     GT     TA     TC     TG     TT
                -22.2, -24.4, -19.9, -22.4, -21.3, -22.2, -22.7, -21.3],
               dtype=np.float32)

#  AA  AC  AG  AT  CA  CC  CG  CT  GA  GC  GG  GT  TA  TC  TG  TT
#  00  01  02  03  10  11  12  13  20  21  22  23  30  31  32  34

# nearest-neighbor energies for Watson-Crick complements at 37C
# (Table 1 in SantaLucia and Hicks 2004)
# ordering of array is
# #                      AA    AC    AG    AT    CA    CC    CG    CT
# _nndGwc = np.array([-1.00,-1.44,-1.28,-0.88,-1.45,-1.84,-2.17,-1.28,
#                     #  GA    GC    GG    GT    TA    TC    TG    TT
#                     -1.30,-2.24,-1.84,-1.44,-0.58,-1.30,-1.45,-1.00],
#                    dtype=np.float32)
#                    # AA   AC   AG   AT   CA   CC   CG   CT
# _nndGwc = np.array([1.00,1.44,1.28,0.88,1.45,1.84,2.17,1.28,
#                    # GA   GC   GG   GT   TA   TC   TG   TT
#                    1.30,2.24,1.84,1.44,0.58,1.30,1.45,1.00],
#                   dtype=np.float32)
# _nndGwcStr = {'AA':1.00,'AC':1.44,'AG':1.28,'AT':0.88,'CA':1.45,'CC':1.84,
#              'CG':2.17,'CT':1.28,'GA':1.30,'GC':2.24,'GG':1.84,'GT':1.44,
#              'TA':0.58,'TC':1.30,'TG':1.45,'TT':1.00}

# nearest-neighbor energies for single mismatches (Table 2 in SantaLucia)
# ordering of array is
#                  # GA/CA GA/CG    AG    AT    CA    CC    CG    CT   GA    GC    GG    GT    TA   TC    TG   TT
# _nndGsmm = np.array([0.17,-1.44,-1.28,-0.88,-1.45,-1.84,-2.17,-1.28,-1.3,-2.24,-1.84,-1.44,-0.58,-1.3,-1.45,-1.0], dtype=np.float32)

_all_pairs = [((i << 2) + j, bits2base[i] + bits2base[j])
              for i in range(4) for j in range(4)]

@lru_cache(maxsize=32)
def calculate_loop_energies_dict(temperature: float, negate: bool = False) -> Dict[str, float]:
    loop_energies = calculate_loop_energies(temperature, negate)
    return {pair[1]: loop_energies[pair[0]] for pair in _all_pairs}

@lru_cache(maxsize=100000)
def wcenergy(seq: str, temperature: float, negate: bool = False) -> float:
    """Return the wc energy of seq binding to its complement."""
    loop_energies = calculate_loop_energies_dict(temperature, negate)
    return sum(loop_energies[seq[i:i + 2]] for i in range(len(seq) - 1))

def wcenergies_str(seqs: Sequence[str], temperature: float, negate: bool = False) -> List[float]:
    seqarr = seqs2arr(seqs)
    return list(calculate_wc_energies(seqarr, temperature, negate))

def wcenergy_str(seq: str, temperature: float, negate: bool = False) -> float:
    seqarr = seqs2arr([seq])
    return list(calculate_wc_energies(seqarr, temperature, negate))[0]

def hash_ndarray(arr: np.ndarray) -> int:
    writeable = arr.flags.writeable
    if writeable:
        arr.flags.writeable = False
    h = hash(bytes(arr.data))  # hash(arr.data)
    arr.flags.writeable = writeable
    return h

CACHE_WC = False
_calculate_wc_energies_cache: Optional[np.ndarray] = None
_calculate_wc_energies_cache_hash: int = 0

def calculate_wc_energies(seqarr: np.ndarray, temperature: float, negate: bool = False) -> np.ndarray:
    """Calculate and store in an array all energies of all sequences in seqarr
    with their Watson-Crick complements."""
    global _calculate_wc_energies_cache
    global _calculate_wc_energies_cache_hash
    if CACHE_WC and _calculate_wc_energies_cache is not None:
        if _calculate_wc_energies_cache_hash == hash_ndarray(seqarr):
            return _calculate_wc_energies_cache
    loop_energies = calculate_loop_energies(temperature, negate)
    left_index_bits = seqarr[:, :-1] << 2
    right_index_bits = seqarr[:, 1:]
    pair_indices = left_index_bits + right_index_bits
    pair_energies = loop_energies[pair_indices]
    energies: np.ndarray = np.sum(pair_energies, axis=1)
    if CACHE_WC:
        _calculate_wc_energies_cache = energies
        _calculate_wc_energies_cache_hash = hash_ndarray(_calculate_wc_energies_cache)
    return energies

def wc_arr(seqarr: np.ndarray) -> np.ndarray:
    """Return numpy array of complements of sequences in `seqarr`."""
    return (3 - seqarr)[:, ::-1]

def prefilter_length_10_11(low_dg: float, high_dg: float, temperature: float, end_gc: bool,
                           convert_to_list: bool = True) \
        -> Union[Tuple[List[str], List[str]], Tuple[DNASeqList, DNASeqList]]:
    """Return sequences of length 10 and 11 with wc energies between given values."""
    s10: DNASeqList = DNASeqList(length=10)
    s11: DNASeqList = DNASeqList(length=11)
    s10 = s10.filter_energy(low=low_dg, high=high_dg, temperature=temperature)
    s11 = s11.filter_energy(low=low_dg, high=high_dg, temperature=temperature)
    forbidden_subs = [f'{a}{b}{c}{d}' for a in ['G', 'C']
                      for b in ['G', 'C']
                      for c in ['G', 'C']
                      for d in ['G', 'C']]
    s10 = s10.filter_substring(forbidden_subs)
    s11 = s11.filter_substring(forbidden_subs)
    if end_gc:
        print(
            'Removing any domains that end in either A or T; '
            'also ensuring every domain has an A or T within 2 indexes of the end')
        s10 = s10.filter_end_gc()
        s11 = s11.filter_end_gc()
    for seqs in (s10, s11):
        if len(seqs) == 0:
            raise ValueError(
                f'low_dg {low_dg:.2f} and high_dg {high_dg:.2f} too strict! '
                f'no sequences of length {seqs.seqlen} found')
    return (s10.to_list(), s11.to_list()) if convert_to_list else (s10, s11)

def all_cats(seq: Sequence[int], seqs: Sequence[int]) -> np.ndarray:
    """
    Return all sequences obtained by concatenating seq to either end of a sequence in seqs.

    For example,

    .. code-block:: Python

        all_cats([0,1,2,3], [[3,3,3], [0,0,0]])

    returns the numpy array

    .. code-block:: Python

        [[0,1,2,3,3,3,3],
         [3,3,3,0,1,2,3],
         [0,1,2,3,0,0,0],
         [0,0,0,0,1,2,3]]
    """
    seqarr = np.asarray([seq])
    seqsarr = np.asarray(seqs)
    ar = seqarr.repeat(seqsarr.shape[0], axis=0)
    ret = np.concatenate((seqsarr, ar), axis=1)
    ret2 = np.concatenate((ar, seqsarr), axis=1)
    ret = np.concatenate((ret, ret2))
    return ret

max_stars_repo_name: dsd-test-nupack-in-github-workflow/dsd

max_issues_repo_name: dsd-test-nupack-in-github-workflow/dsd

max_forks_repo_name: dsd-test-nupack-in-github-workflow/dsd


text:
import numpy as np
from tools import get_array, dags

data = dags[0]

find_or_add_index = get_array(data, "find_or_add index")
find_or_add_time = get_array(data, "find_or_add time")
find_or_add_add = get_array(data, "find_or_add add")
find_or_add_level = get_array(data, "find_or_add level")

indices_add = find_or_add_add == 1
indices_noadd = find_or_add_add == 0

print("Num add: ", np.sum(indices_add))
print("Num no add: ", np.sum(indices_noadd))

assert np.sum(indices_add) + np.sum(indices_noadd) == len(find_or_add_index)

for level in range(32):
    level_indices = level == find_or_add_level
    if np.sum(level_indices) == 0:
        continue

    print("level {}: mean no add: {}; mean add: {}; sum no add: {}; sum add: {}".format(
        level,
        np.mean(find_or_add_time[np.logical_and(indices_noadd, level_indices)]),
        np.mean(find_or_add_time[np.logical_and(indices_add, level_indices)]),
        np.sum(find_or_add_time[np.logical_and(indices_noadd, level_indices)]),
        np.sum(find_or_add_time[np.logical_and(indices_add, level_indices)])))

ext: py

max_stars_repo_name: MaximeDup/HashDAG

max_issues_repo_name: MaximeDup/HashDAG

max_forks_repo_name: MaximeDup/HashDAG

reason: import numpy


text:
#!/usr/bin/env python

import networkx as nx
import sys

with open('soc-sign-bitcoinotc.csv', 'r') as fin:
	with open('bitcoinotc.in', 'a') as fout:
		G = nx.DiGraph()
		fout.write("Relation Edges" + '\n')
		fout.write("From To" + '\n')
		for line in fin:
			n1, n2, rating, time = line.split(',')
			fout.write(str(n1) + ' ' + str(n2) + ' ' + str(10 - int(rating)) + '\n')
			G.add_edge(int(n1), int(n2))
		fout.write("End of Edges" + '\n')

		print "===== BitcoinOTC ===="
		print "Number of nodes = " + str(nx.number_of_nodes(G))
		print "Number of edges = " + str(nx.number_of_edges(G))
		degrees = nx.degree(G)
		degree_list = [d for (_, d) in degrees]
		print "Max degree = " + str(max(degree_list))
		avg_degree = sum(degree_list) * 1.0 / len(degree_list)
		print "Average degree = " + str(avg_degree)

max_stars_repo_name: northeastern-datalab/anyk-code

max_issues_repo_name: northeastern-datalab/anyk-code

max_forks_repo_name: northeastern-datalab/anyk-code

reason: import networkx


text:
import numpy as np

class Util:
    @staticmethod
    # Maakt wat variaties op de wachtwoorden
    def getCapitalPermutations(words):
        return np.concatenate([(word, word.upper(), word.lower(), word.capitalize())
                               for word in words])

class Wordlist:
    def __init__(self,
                 preWords:tuple=None,
                 midWords:tuple=None,
                 endWords:tuple=None,
                 sorting:bool=False,
                 filename:str="wwNL",
                 ):
        # De officiele oplossing voor lijst-type parameter mutability oplossen is 2 lines
        # maar ik vond een manier met 1-line yesss
        preWords = {None: []}.get(preWords, preWords)
        midWords = {None: []}.get(midWords, midWords)
        endWords = {None: []}.get(endWords, endWords)
        self.filename = filename

        # Start woorden (Het grootste stuk van een wachtwoord)
        self._preWords = set(Util.getCapitalPermutations((
            *preWords,
            # Basis dingen
            "welkom",
            "hallo",
            "wachtwoord",
            "beveiliging",
            "geheim",
            "veilig",
            "1234hoedjevanpapier",
            "1234HoedjeVanPapier",
            "vakantie",
            "weekend",
            # Seizoenen
            "herfst",
            "winter",
            "lente",
            "zomer",
            # Maanden
            "januari",
            "februari",
            "maart",
            "april",
            "mei",
            "juni",
            "juli",
            "augustus",
            "september",
            "oktober",
            "november",
            "december",
            # Dag-namen
            "MaanDag",
            "DinsDag",
            "WoensDag",
            "DonderDag",
            "VrijDag",
            "ZaterDag",
            "ZonDag",
            # Voetbal clubs
            "Feyenoord",
            "SC Heerenveen",
            "SCHeerenveen",
            "Ajax",
            "PSV",
            "FC Twente",
            "FCTwente",
            "Vitesse",
            "SC Cambuur",
            "SCCambuur",
            # Bekendheden
            "WillemAlexander",
            "Willem-Alexander",
            "Willem",
            "Alexander",
            "Maxima",
            "Rutte",
            "Geert",
            "GeertWilders",
            "Wilders",
            "Beatrix",
            # Artiesten / Performers
            "Vengaboys",
            "ArminVanBuuren",
            "MartinGarrix",
            "Martin",
            "Garrix",
            "MaxVerstappen",
            "Verstappen",
        )))
        # Achtervoegsel (Dingen zoals 01 of 123 om net te voldoen aan eisen zoals 8-tekens)
        self._midWords = {
            *midWords,
            # Blank, anders kunnen losse preWords niet gebruikt worden
            ""
            # Makkelijk om te typen
            "123",
            # Dit getal staat voor iets ongepasts - https://en.wikipedia.org/wiki/69
            "69",
            # Weed getal - https://www.urbandictionary.com/define.php?term=420
            "420",
            # Duivel getal - https://en.wikipedia.org/wiki/Number_of_the_beast
            "666",
            # Elite getal - https://www.urbandictionary.com/define.php?term=1337
            "1337",
            # Afkomstig van een Dragon-ball Z meme - https://en.wikipedia.org/wiki/It%27s_Over_9000!
            "9000",
            "9001",
            # 0 - 52 want er zijn 52 weken in een jaar
            *map(lambda x: f"{x}".rjust(2, "0"), range(0, 53)),
            *map(lambda x: f"{x}", range(0, 53)),
            # Wachtwoorden met datums staan altijd ergens voor
            # En meestal is het verjaardagen, geboorte datums enz,
            # Dus alles binnen 100 jaar van vandaag is goed
            *map(lambda x: f"{x}", range(1900, 2101)),
        }
        # Extra's (Sommige mensen denken dat 1 symbool toevoegen het wachtwoord goed maakt)
        self._endWords = {
            *endWords,
            # Blank
            "",
            # Meest veelvoorkomende symboolen in wachtwoorden
            "!",
            "@",
            "#",
            # Dit gebeurt als je 123 indrukt met shift
            "!@#",
        }
        # set() wordt hier niet gebruikt voor deduplicatie maar om het een __len__ te geven
        # ⚠ Door set() is de volgorde elke keer willekeurig
        self.wordlist = set((
            # Dit is 2x sneller dan string formatting
            pre + mid + end
            for pre in self._preWords
            for mid in self._midWords
            for end in self._endWords
        ))
        if sorting:
            self.wordlist = sorted(self.wordlist)

    def writeFile(self):
        fn = f"{self.filename}-{int(len(self.wordlist)/1000)}k.txt"
        with open(fn, "w") as f:
            f.write("\n".join(self.wordlist))
        print(f"Bestand {fn} opgeslagen")

    @property
    # Laat wachtwoorden in notatie zien
    def passwordsNotation(self):
        return (f"[{'/'.join(sorted(self._preWords))}]\n"
                f"[{'/'.join(sorted(self._midWords))}]\n"
                f"[{'/'.join(sorted(self._endWords))}]")

max_stars_repo_name: Walter-o/WalterWordlistNL

max_issues_repo_name: Walter-o/WalterWordlistNL

max_forks_repo_name: Walter-o/WalterWordlistNL


text:
from pycalphad import Database, Model, variables as v
from pycalphad.tests.datasets import ROSE_TDB
from pycalphad import equilibrium
import numpy as np

my_phases_rose = ['TEST']
comps = ['H', 'HE', 'LI', 'BE', 'B', 'C', 'N', 'O', 'F']
comps = sorted(comps)
conds = dict({v.T: 1000, v.P: 101325})
for comp in comps[1:]:
    conds[v.X(comp)] = 1.0/float(len(comps))
dbf = Database(ROSE_TDB)
eqx = equilibrium(Database(ROSE_TDB), comps, my_phases_rose, conds, calc_opts={'pdens': 10}, verbose=True)

max_stars_repo_name: richardotis/pycalphad-sandbox

max_issues_repo_name: richardotis/pycalphad-sandbox

max_forks_repo_name: richardotis/pycalphad-sandbox


text:
from os.path import expanduser
home = expanduser('~')
import numpy as np
import random
from data_utils import get_pretrain_folder
from conf_utils import get_conf_base
from keras.optimizers import Adam, SGD, RMSprop

class Conf(object):
    def __init__(self, data_name, param_dict=None):
        self.data_name = data_name
        self.max_epoch = 20
        self.num_negatives = 5
        self.batch_size_p = 64
        self.eval_topk = 50

        self.loss = 'skip-gram'
        #self.loss = 'max-margin'
        #self.loss_margin = 1
        #self.loss = 'log-loss'
        self.neg_loss_weight = 1
        self.use_content_id = False
        self.v_reg = 0
        self.interaction_multiplier = False
        self.interaction_bias = 'item'
        self.emb_normalization = None

        self.learn_rate = 0.01
        self.user_dim = self.item_dim = 50
        self.u_reg = 1e-5
        self.neg_dist = 'uniform'
        self.neg_sampling_power = 1
        self.chop_size = 1
        self.shuffle_st = 'by_item'

        self.evaluation_mode = False

        if param_dict is not None:
            self.__dict__.update(param_dict)
        self._post_init()

    def _post_init(self):
        if self.emb_normalization is None:
            self.emb_normalization = True \
                if self.loss == 'max-margin' or self.loss == 'log-loss' \
                else False

        if self.learn_rate > 0:
            self.optimizer = Adam(self.learn_rate)
        else:
            from optimizer import AdamOptimizer
            #import tensorflow as tf
            #self.optimizer = tf.train.GradientDescentOptimizer(1e4)
            self.optimizer = AdamOptimizer(-self.learn_rate)

        pretrain_folder = get_pretrain_folder(self.data_name, aug=False)
        sentvec_filepath = pretrain_folder + 'sentence_vectors_50d.txt'
        self.pretrain = {'wordvec_filepath': None,
                         'sentvec_filepath': sentvec_filepath,
                         'transform': True,
                         'pretrain_combine_dropout': 0.5,
                         'pretrain_combine_actv': 'relu',
                         'pretrain_combine_mode': None}
        self.pretrain['transform'] = False

        self.contextual_spatial_gated_input = None
        self.contextual_temporal_gated_input = None

def get_conf_default(data_name, param_dict=None):
    conf_seed, conf_var = get_conf_base(param_dict)
    conf = Conf(data_name, param_dict=param_dict)
    if conf_seed:
        conf.conf_seed = conf_seed
    if conf_var:
        conf.conf_var = conf_var
    return conf

def get_conf_evaluation(data_name, param_dict=None):
    # user/word dim is fixed to 50 for comparisons
    conf_seed, conf_var = get_conf_base(param_dict)
    conf = Conf(data_name, param_dict=param_dict)
    if conf_seed:
        conf.conf_seed = conf_seed
    if conf_var:
        conf.conf_var = conf_var
    conf.evaluation_mode = True
    conf.max_epoch = 3
    conf.interaction_bias = None
    conf.emb_normalization = False
    conf.optimizer = SGD(0)
    data_name = data_name.replace('_fold1', '')
    emb_folder = home + '/pbase/x/p607/results-aaai/baselines/CDL/results/{}/fold1/cdl/'.format(data_name)
    if data_name.startswith('citeulike'):
        user_emb_file = emb_folder + '10-U.dat'
        item_emb_file = emb_folder + '10-V.dat'
    else:
        user_emb_file = emb_folder + 'final-U.dat'
        item_emb_file = emb_folder + 'final-V.dat'
    user_emb = np.loadtxt(user_emb_file)
    item_emb = np.loadtxt(item_emb_file)
    conf.pretrain = {'wordvec_filepath': None,
                     'sentvec_filepath': None,
                     'user_emb': user_emb,
                     'item_emb': item_emb,
                     'transform': False,
                     'pretrain_combine_dropout': 0.,
                     'pretrain_combine_actv': 'linear',
                     'pretrain_combine_mode': None}
    return conf

def get_conf_best(data_name, param_dict=None):
    # user/word dim is fixed to 50 for comparisons
    conf_seed, conf_var = get_conf_base(param_dict)
    conf = Conf(data_name, param_dict=param_dict)
    if conf_seed:
        conf.conf_seed = conf_seed
    if conf_var:
        conf.conf_var = conf_var
    conf.c_reg = 0
    conf.max_epoch = 30
    conf.num_negatives = 10
    conf.interaction_bias = None
    conf.pretrain['transform'] = True
    if data_name.startswith('citeulike_title_only'):
        conf.u_reg = 1e-6
        conf.interaction_multiplier = False
        conf.pretrain['pretrain_combine_dropout'] = 0.3
    elif data_name.startswith('citeulike_title_and_abstract'):
        conf.u_reg = 1e-6
        conf.interaction_multiplier = False
        conf.pretrain['pretrain_combine_dropout'] = 0.3
    elif data_name.startswith('news_title_only'):
        conf.u_reg = 1e-5
        conf.interaction_multiplier = False
        conf.pretrain['pretrain_combine_dropout'] = 0.1
    elif data_name.startswith('news_title_and_abstract'):
        conf.u_reg = 1e-5
        conf.interaction_multiplier = False
        conf.pretrain['pretrain_combine_dropout'] = 0.1
    else:
        assert False, 'unknown data_name: %s' % data_name
    try:
        param_dict['reset_after_getconf']
        conf.__dict__.update(param_dict)
    except:
        pass
    return conf

def get_conf_random(data_name, param_dict=None):
    conf_seed, conf_var = get_conf_base(param_dict)
    if conf_seed:
        random.seed(conf_seed)
    conf = Conf(data_name, param_dict=param_dict)
    if conf_seed:
        conf.conf_seed = conf_seed
    if conf_var:
        conf.conf_var = conf_var
    conf.u_reg = random.choice([0., 1e-5, 1e-6, 1e-7])
    conf.pretrain['transform'] = random.choice([True, False])
    if conf.pretrain['transform']:
        conf.pretrain['pretrain_combine_dropout'] = [0.1, 0.3, 0.5]  # choose one?

    return conf

def get_conf(data_name, conf_choice, param_dict=None):
    if conf_choice == 'best':
        conf = get_conf_best(data_name, param_dict=param_dict)
    elif conf_choice == 'default':
        conf = get_conf_default(data_name, param_dict=param_dict)
    elif conf_choice == 'random':
        conf = get_conf_random(data_name, param_dict=param_dict)
    elif conf_choice == 'evaluation':
        conf = get_conf_evaluation(data_name, param_dict=param_dict)
    else:
        assert False, '[Error!] conf_choice {} unknonw'.format(conf_choice)
    return conf

max_stars_repo_name: lastBird/NNCF

max_issues_repo_name: lastBird/NNCF

max_forks_repo_name: lastBird/NNCF


text:
import cv2
import numpy as np

from converter_types import ConverterType

print(ConverterType.GrayScaleCV.value)

def convert(image, type):
    if type == ConverterType.GrayScaleCV.value:
        image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
    elif type == ConverterType.RGB_Red.value:
        image = image[:, :, 2]
    elif type == ConverterType.RGB_Green.value:
        image = image[:, :, 1]
    elif type == ConverterType.RGB_Blue.value:
        image = image[:, :, 0]
    elif type == ConverterType.HSV_Hue.value:
        hsv_image = cv2.cvtColor(image, cv2.COLOR_BGR2HSV)
        image = hsv_image[:, :, 0]
    elif type == ConverterType.GrayScaleMid.value:
        image = (np.int32(image[:, :, 0]) + image[:, :, 1] + image[:, :, 2]) // 3
    else:
        raise Exception("Incorrect converter type")

    return image


text:
# Import modules
import pickle
import pandas as pd
from psychopy import visual, monitors
from psychopy import core, event
import numpy as np
from titta import Titta, helpers_tobii as helpers

#%% Monitor/geometry participant screen
MY_MONITOR                  = 'testMonitor' # needs to exists in PsychoPy monitor center
FULLSCREEN                  = False
SCREEN_RES                  = [1920, 1080]
SCREEN_WIDTH                = 52.7 # cm
VIEWING_DIST                = 63 #  # distance from eye to center of screen (cm)

mon = monitors.Monitor(MY_MONITOR)  # Defined in defaults file
mon.setWidth(SCREEN_WIDTH)          # Width of screen (cm)
mon.setDistance(VIEWING_DIST)       # Distance eye / monitor (cm)
mon.setSizePix(SCREEN_RES)

# Monitor/geometry operator screen
MY_MONITOR_OP                  = 'default' # needs to exists in PsychoPy monitor center
FULLSCREEN_OP                  = False
SCREEN_RES_OP                  = [1920, 1080]
SCREEN_WIDTH_OP                = 52.7 # cm
VIEWING_DIST_OP                = 63 #  # distance from eye to center of screen (cm)

mon_op = monitors.Monitor(MY_MONITOR_OP)  # Defined in defaults file
mon_op.setWidth(SCREEN_WIDTH_OP)          # Width of screen (cm)
mon_op.setDistance(VIEWING_DIST_OP)       # Distance eye / monitor (cm)
mon_op.setSizePix(SCREEN_RES_OP)

# Window set-up (this color will be used for calibration)
win = visual.Window(monitor = mon, fullscr = FULLSCREEN,
                    screen=1, size=SCREEN_RES, units = 'deg')

win_op = visual.Window(monitor = mon_op, fullscr = FULLSCREEN_OP,
                    screen=0, size=SCREEN_RES_OP, units = 'norm')

fixation_point = helpers.MyDot2(win)
image = visual.ImageStim(win, image='im1.jpeg', units='norm', size = (2, 2))

#%% ET settings
et_name = 'Tobii Pro Spectrum'
# et_name = 'Tobii4C'

dummy_mode = False
bimonocular_calibration = False

# Change any of the default dettings?e
settings = Titta.get_defaults(et_name)
settings.FILENAME = 'testfile.tsv'

#%% Connect to eye tracker and calibrate
tracker = Titta.Connect(settings)
if dummy_mode:
    tracker.set_dummy_mode()
tracker.init()

# Calibrate
if bimonocular_calibration:
    tracker.calibrate(win, win_operator=win_op, eye='left', calibration_number = 'first')
    tracker.calibrate(win, win_operator=win_op, eye='right', calibration_number = 'second')
else:
    tracker.calibrate(win, win_operator=win_op)

#%% Record some data
tracker.start_recording(gaze_data=True, store_data=True)

# Present fixation dot and wait for one second
fixation_point.draw()
t = win.flip()
tracker.send_message('fixation target onset')
core.wait(1)
tracker.send_message('fixation target offset')

image.draw()
t = win.flip()
tracker.send_message('image onset')
core.wait(3)
tracker.send_message('image offset')

win.flip()

tracker.stop_recording(gaze_data=True)

# Close window and save data
win.close()
win_op.close()
tracker.save_data()

#%% Open pickle and write et-data and messages to tsv-files.
f = open(settings.FILENAME[:-4] + '.pkl', 'rb')
gaze_data = pickle.load(f)
msg_data = pickle.load(f)

# Save data and messages
df = pd.DataFrame(gaze_data, columns=tracker.header)
df.to_csv(settings.FILENAME[:-4] + '.tsv', sep='\t')
df_msg = pd.DataFrame(msg_data,  columns = ['system_time_stamp', 'msg'])
df_msg.to_csv(settings.FILENAME[:-4] + '_msg.tsv', sep='\t')

core.quit()

max_stars_repo_name: mintzer/pupillometry-rf-back

max_issues_repo_name: mintzer/pupillometry-rf-back

max_forks_repo_name: mintzer/pupillometry-rf-back


text:
import logging

from webthing import (Action, Event, Property, MultipleThings, Thing, Value,
                      WebThingServer)

from security.dlm import DLM

import tornado.ioloop
import random
import numpy as np

import uuid

# For debug
# use_static = True
# static_values = {
#     'urn:dev:ops:fake-position-1':1,
#     'urn:dev:ops:fake-position-2':2,
#     'urn:dev:ops:fake-position-3':3,
#     'urn:dev:ops:fake-position-4':3,
#     'urn:dev:ops:fake-position-5':2,
#     'urn:dev:ops:fake-position-6':4,
# }

class FakePulseSensor(Thing):
    """A sensor that detects movment."""
    _id_prefix = 'mock:pulse:' + str(uuid.uuid4())
    _instance_number = 0
    _values = [1,2,3,4]
    _max_value = 170
    _min_value = 40

    def __init__(self, device_id, classification=None, name=None):
        type(self)._instance_number+=1
        self.device_id = device_id #FakePulseSensor._id_prefix + str(type(self)._instance_number) #id
        self.name = name if name else self.device_id
        self.classification = classification if classification else DLM()

        print("Create new device:", self.name,"With classification:",classification)
        # self._position_index = random.randint(0, len(FakePositionDevice._values)-1)
        # self._position = FakePositionDevice._values[self._position_index] # Initial value
        self._value = self.getInitialValue()
        Thing.__init__(
            self,
            self.device_id,
            name, #title
            ['DiscretePostitionDevice'], #type
            'A device reporting its position', #description
            classification = classification
        )

        self.value = Value(self._value)
        self.add_property(
            Property(self,
                     'pulse',
                     self.value,
                     metadata={
                         '@type': 'PulseData',
                         'title': 'Pulse',
                         'type': 'float',
                         'description': 'Real time pulse value',
                         'readOnly': True,
                     }))

        logging.debug('starting the sensor update looping task')
        self.timer = tornado.ioloop.PeriodicCallback(
            self.update_level,
            3000
        )
        self.timer.start()

    def update_level(self):
        new_value = self.read_from_gpio()
        logging.debug('New pulse value: %s', new_value)
        self.value.notify_of_external_update(new_value)

    def cancel_update_level_task(self):
        self.timer.stop()

    # def _decrease_position_index(self):
    #     if (self._position_index==0):
    #         self._position_index = len(FakePositionDevice._values)-1
    #     else:
    #         self._position_index-+1

    # def _increase_position_index(self):
    #     self._position_index = (self._position_index+1)%len(FakePositionDevice._values)

    def getInitialValue(self):
        return 55

    # @staticmethod
    def read_from_gpio(self):
        """Mimic an actual sensor updating its reading every couple seconds."""
        # if (use_static):
        #     if (self.device_id in static_values.keys()):
        #         return static_values[self.device_id]
        #     else:
        #         return 1;

        # if (random.random()<0.2):
        #     if (random.random()<0.5):
        #         self._decrease_position_index()
        #     else:
        #         self._increase_position_index()
                # self._position_index = (self._position_index+1)%len(FakePositionDevice._values)

        # self.position = FakePositionDevice._values[self._position_index]
        mu = self._value
        sig = 1
        self._value = np.random.normal(mu,sig,1)[0]
        if self._value < self._min_value:
            self._value = self._min_value
        if self._value > self._max_value:
            self._value = self._max_value

        return self._value
        return 4
        # return FakePositionDevice._values[self._position_index]

max_stars_repo_name: marbirg/secure-iot-dashboard-prototype

max_issues_repo_name: marbirg/secure-iot-dashboard-prototype

max_forks_repo_name: marbirg/secure-iot-dashboard-prototype


text:
import time

import numpy as np

def pp_time(seconds_remaining, space=5):
    """
    Pretty prints the given time in the largest useful unit
    If you get to years, good luck
    :param seconds_remaining: The number of seconds to be pretty printed
    :param space: The number of characters the number should take. Allows for constant width printing, potentially more
    visually appealing
    """
    time_counters = [
        ("s", 60),
        ("m", 60),
        ("h", 24),
        ("d", 30),
        (" months", 12),
    ]
    time_remaining = seconds_remaining
    for s, t in time_counters:
        if time_remaining > t:
            time_remaining /= t
        else:
            return f"{time_remaining: {space}.2f}{s}"
    return f"{time_remaining: {space}.2f} years"

class ProgressTracker:
    """
    A progress tracker that prints the current status of a job being run.
    Requires knowing how many tasks/items are going to be run and computes estimated time remaining assuming all items
    are the same size.
    Can be used as a context manager where the update function is assigned to the named variable.
    Note that using the context manager prevents reviewing the times after running
    Example:
    n_items = len(work)
    with ProgressTracker(n_items) as pt:
        for i, job in enumerate(work):
            func(job)
            pt(i)
    """

    def __init__(self, n_items, percent_increment=10, print_func=print):
        """
        Creates the progress tracker
        :param n_items: the number of items to be worked on.
        :param percent_increment: How often to print the progress. Will output when the percent complete first hits each
        integer multiple of this value
        :param print_func: The function to call when outputting the progress. Allows for logging functions to be used
        instead of print
        """
        self.n_items = n_items
        self.percent_increment = percent_increment
        self.current_increment = percent_increment
        self.started = False
        self.times = []
        self.print_func = print_func

    def __enter__(self):
        """
        Context manager entry function
        Starts the timer and returns the update function
        """
        self._start()
        return self.update

    def _start(self):
        """
        Start the timer
        Silently exits if called again
        """
        if not self.started:
            self.times.append(time.time())
            self.started = True

    def __exit__(self, exc_type, exc_val, exc_tb):
        """
        Context manager exit method
        Converts the times from a list to an array
        """
        self.update(self.n_items)
        self._stop()

    def _stop(self):
        """
        Converts the list of absolute times to an array of times relative to the start time
        Allows for run time analysis later if required
        Note that times is not available if the tracker is used with the context manager
        """
        self.times = np.asarray(self.times) - self.times[0]

    def update(self, update_index):
        """
        Updates the tracker with the index of the item currently being run
        Adds the time to the list of times and outputs a string with percent completion, time taken since last update,
        time take since start, estimated time remaining and estimated total time
        """
        if not self.started:
            self._start()

        percent_done = 100 * update_index / self.n_items
        if percent_done >= self.current_increment:
            self.times.append(time.time())
            while percent_done >= self.current_increment:
                self.current_increment += self.percent_increment

            elapsed_time = self.times[-1] - self.times[0]
            step_time = self.times[-1] - self.times[-2]
            remaining_time = elapsed_time * self.n_items / update_index - elapsed_time
            self.print_func(
                f"{percent_done:5.1f}% complete. Time taken for this block: {pp_time(step_time)}. "
                f"Time elapsed: {pp_time(elapsed_time)}. Time remaining: {pp_time(remaining_time)}. "
                f"Total time: {pp_time(elapsed_time + remaining_time)}."
            )

max_stars_repo_name: QuakeKorea/qcore

max_issues_repo_name: QuakeKorea/qcore

max_forks_repo_name: QuakeKorea/qcore


text:
from __future__ import print_function

from builtins import str, zip

def interp_to_obs(var, df, lat, lon, radius=12000.):
    """Short summary.

    Parameters
    ----------
    var : type
        Description of parameter `var`.
    df : type
        Description of parameter `df`.
    lat : type
        Description of parameter `lat`.
    lon : type
        Description of parameter `lon`.
    radius : type
        Description of parameter `radius` (the default is 12000.).

    Returns
    -------
    type
        Description of returned object.

    """
    from numpy import NaN, vstack
    from pyresample import geometry, image
    from pandas import to_timedelta, DataFrame
    # define CMAQ pyresample grid (source)
    grid1 = geometry.GridDefinition(lons=lon, lats=lat)
    # get unique sites from df
    dfn = df.drop_duplicates(subset=['Latitude', 'Longitude'])
    # define site grid (target)
    lats = dfn.Latitude.values
    lons = dfn.Longitude.values
    grid2 = geometry.GridDefinition(lons=vstack(lons), lats=vstack(lats))
    # Create image container
    i = image.ImageContainerNearest(var.transpose('y', 'x', 'time').values, grid1, radius_of_influence=radius,
                                    fill_value=NaN)
    # resample
    ii = i.resample(grid2).image_data.squeeze()
    # recombine data
    e = DataFrame(ii, index=dfn.SCS, columns=var.time.values)
    w = e.stack().reset_index().rename(columns={'level_1': 'datetime', 0: 'model'})
    w = w.merge(dfn.drop(['datetime', 'datetime_local', 'Obs'], axis=1), on='SCS', how='left')
    w = w.merge(df[['datetime', 'SCS', 'Obs']], on=['SCS', 'datetime'], how='left')
    # calculate datetime local

    w['datetime_local'] = w.datetime + to_timedelta(w.utcoffset, 'H')

    return w

def find_nearest_latlon_xarray(arr, lat=37.102400, lon=-76.392900, radius=12e3):
    """Short summary.

    Parameters
    ----------
    arr : type
        Description of parameter `arr`.
    lat : type
        Description of parameter `lat` (the default is 37.102400).
    lon : type
        Description of parameter `lon` (the default is -76.392900).
    radius : type
        Description of parameter `radius` (the default is 12e3).

    Returns
    -------
    type
        Description of returned object.

    """
    from pyresample import utils, geometry
    from numpy import array, vstack
    grid1 = geometry.GridDefinition(lons=arr.longitude, lats=arr.latitude)
    grid2 = geometry.GridDefinition(lons=vstack([lon]), lats=vstack([lat]))
    row, col = utils.generate_nearest_neighbour_linesample_arrays(grid1, grid2, radius)
    row = row.flatten()
    col = col.flatten()
    return arr.sel(x=col).sel(y=row).squeeze()

def get_smops_area_def(nx=1440, ny=720):
    """Short summary.

    Parameters
    ----------
    nx : type
        Description of parameter `nx` (the default is 1440).
    ny : type
        Description of parameter `ny` (the default is 720).

    Returns
    -------
    type
        Description of returned object.

    """
    p = Proj(proj='eqc', lat_ts=0., lat_0=0., lon_0=0., x_0=0., y_0=0., a=6378137, b=6378137, units='m')
    proj4_args = p.srs
    area_name = 'Global .25 degree SMOPS Grid'
    area_id = 'smops'
    proj_id = area_id
    aa = p([-180, 180], [-90, 90])
    area_extent = (aa[0][0], aa[1][0], aa[0][1], aa[1][1])
    area_def = utils.get_area_def(area_id, area_name, proj_id, proj4_args, nx, ny, area_extent)
    return area_def

def get_gfs_area_def(nx=1440, ny=721):
    """Short summary.

    Parameters
    ----------
    nx : type
        Description of parameter `nx` (the default is 1440).
    ny : type
        Description of parameter `ny` (the default is 721).

    Returns
    -------
    type
        Description of returned object.

    """
    # proj4_args = '+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m'
    p = Proj(proj='eqc', lat_ts=0., lat_0=0., lon_0=0., x_0=0., y_0=0., a=6378137, b=6378137, units='m')
    proj4_args = p.srs
    area_name = 'Global .25 degree SMOPS Grid'
    area_id = 'smops'
    proj_id = area_id
    aa = p([0, 360 - .25], [-90, 90.])
    area_extent = (aa[0][0], aa[1][0], aa[0][1], aa[1][1])
    area_def = utils.get_area_def(area_id, area_name, proj_id, proj4_args, nx, ny, area_extent)
    return area_def

def geotiff_meta_to_areadef(meta):
    """
    Transform (Rasterio) geotiff meta dictionary to pyresample area definition
    Arguments:
     meta (dictionary) : dictionary containing projection and image geometry
                         information (formed by Rasterio)
    Returns:
         area_def (pyresample.geometry.AreaDefinition) : Area definition object
    """
    import pyresample
    area_id = ""
    name = ""
    proj_id = "Generated from GeoTIFF"
    proj_dict = meta['crs']
    proj_dict_with_string_values = dict(
        list(zip([str(key) for key in list(proj_dict.keys())], [str(value) for value in list(proj_dict.values())])))
    x_size = meta['width']
    x_res = meta['transform'][0]
    y_res = meta['transform'][4] * -1
    y_size = meta['height']
    x_ll = meta['transform'][2]
    y_ur = meta['transform'][5]
    y_ll = y_ur - y_size * y_res
    x_ur = x_ll + x_size * x_res
    area_extent = [x_ll, y_ll, x_ur, y_ur]
    print(area_extent, x_size, y_size, x_res, y_res)

    area_def = pyresample.geometry.AreaDefinition(area_id,
                                                  name,
                                                  proj_id,
                                                  proj_dict_with_string_values,
                                                  x_size,
                                                  y_size,
                                                  area_extent)
    print(area_extent, xsize, ysize)
    return area_def

def geotiff_meta_to_areadef2(meta):
    """
    Transform (Rasterio) geotiff meta dictionary to pyresample area definition
    Arguments:
     meta (dictionary) : dictionary containing projection and image geometry
                         information (formed by Rasterio)
    Returns:
         area_def (pyresample.geometry.AreaDefinition) : Area definition object
    """
    import pyresample
    area_id = ""
    name = ""
    proj_id = "Generated from GeoTIFF"
    proj_dict = meta['crs']
    proj_dict_with_string_values = dict(
        list(zip([str(key) for key in list(proj_dict.keys())], [str(value) for value in list(proj_dict.values())])))
    x_size = meta['width']
    x_res = 50000.
    y_res = 50000.
    y_size = meta['height']
    x_ll = meta['transform'][2]
    y_ur = meta['transform'][5]
    y_ll = y_ur - y_size * y_res
    x_ur = x_ll + x_size * x_res
    area_extent = [x_ll, y_ll, x_ur, y_ur]
    print(area_extent, x_size, y_size, x_res, y_res)

    area_def = pyresample.geometry.AreaDefinition(area_id,
                                                  name,
                                                  proj_id,
                                                  proj_dict_with_string_values,
                                                  x_size,
                                                  y_size,
                                                  area_extent)
    return area_def
    """Short summary.

    Parameters
    ----------
    meta : type
        Description of parameter `meta`.

    Returns
    -------
    type
        Description of returned object.

    """

def get_grid_def(lon, lat):
    """Short summary.

    Parameters
    ----------
    lon : type
        Description of parameter `lon`.
    lat : type
        Description of parameter `lat`.

    Returns
    -------
    type
        Description of returned object.

    """
    return geometry.GridDefinition(lons=lon, lats=lat)

max_stars_repo_name: barronh/MONET

max_issues_repo_name: barronh/MONET

max_forks_repo_name: barronh/MONET

reason: from numpy


text:
#%%
import os
from pkgutil import extend_path
import sys
sys.path.append('/Users/tsukada/git/RAFT/core')
import argparse

import cv2
import glob
import torch
import numpy as np
import xarray as xr
from PIL import Image
import matplotlib.pyplot as plt

from raft import RAFT
from utils import flow_viz
from utils.utils import InputPadder

%load_ext autoreload
%autoreload 2

#%%
DEVICE = 'cpu'

def load_image(imfile):
    img = np.array(Image.open(imfile)).astype(np.uint8)
    if img.shape[-1] >= 4:
        img = img[:,:,:3] # delete alpha cnannel
    img = torch.from_numpy(img).permute(2, 0, 1).float()
    return img[None].to(DEVICE)

def viz(img, flo, fname="flow.png"):
    img = img[0].permute(1,2,0).cpu().numpy()
    flo = flo[0].permute(1,2,0).cpu().numpy()

    # map flow to rgb image
    flo = flow_viz.flow_to_image(flo)
    # img_flo = np.concatenate([img, flo], axis=1)

    import matplotlib.pyplot as plt
    fig, ax = plt.subplots(nrows=1, ncols=2)
    ax[0].imshow(img / 255.0)
    ax[1].imshow(flo / 255.0)
    ax[0].axis("off")
    ax[1].axis("off")
    fig.savefig(fname, bbox_inches="tight", pad_inches=0.1, dpi=144,)
    plt.close()

    # cv2.imshow('image', img_flo[:, :, [2,1,0]]/255.0)
    # cv2.waitKey()

def flow_to_netcdf(flow_uv, fname="flow.nc"):
    flow_uv = flow_uv[0].permute(1,2,0).cpu().numpy()
    u = flow_uv[:,:,0]
    v = -flow_uv[:,:,1] # vは画像下向き(南向き)が正なので反転
    ds = xr.Dataset(
        {"u": (['y','x'], u), "v": (['y','x'], v)},
        coords=({"y":np.r_[0:u.shape[0]],"x":np.r_[0:u.shape[1]]}))
    ds.to_netcdf(fname, mode="w")

def demo(args):
    model = torch.nn.DataParallel(RAFT(args))
    # model.load_state_dict(torch.load(args.model))
    model.load_state_dict(torch.load(args.model, map_location=torch.device('cpu')))

    model = model.module
    model.to(DEVICE)
    model.eval()

    uv_datadir = args.project+f"/results/{os.path.basename(args.model).split('.')[0]}/uv"
    os.makedirs(uv_datadir, exist_ok=True)
    with torch.no_grad():
        images = glob.glob(os.path.join(args.path, '*.png')) + \
                 glob.glob(os.path.join(args.path, '*.jpg'))

        images = sorted(images)
        flow_nc_names = [""]*len(images)
        for i, (imfile1, imfile2) in enumerate(zip(images[:-1], images[1:])):
            image1 = load_image(imfile1)
            image2 = load_image(imfile2)

            padder = InputPadder(image1.shape)
            image1, image2 = padder.pad(image1, image2)

            if args.warm_start and i > 0:
                flow_low, flow_up = model(image1, image2, iters=20, test_mode=True, flow_init=flow_low)
            else:
                flow_low, flow_up = model(image1, image2, iters=20, test_mode=True)
            flow_nc_name = uv_datadir+"/"+os.path.basename(imfile1).replace(".png",".nc").replace(".jpg",".nc")
            flow_nc_names[i] = flow_nc_name
            flow_to_netcdf(flow_up, flow_nc_name)
        return images, flow_nc_names

parser = argparse.ArgumentParser()
parser.add_argument('--model', help="restore checkpoint")
parser.add_argument('--path', help="dataset for evaluation")
parser.add_argument('--project', help="project directory")
parser.add_argument('--small', action='store_true', help='use small model')
parser.add_argument('--mixed_precision', action='store_true', help='use mixed precision')
parser.add_argument('--alternate_corr', action='store_true', help='use efficent correlation implementation')
parser.add_argument('--warm_start', action='store_true', help='use previous flow as next step initial flow')
args = parser.parse_args([
    '--model', 'models/raft-things.pth',
    '--path', '/Users/tsukada/git/RAFT/projects/test_project/frames',
    '--project', '/Users/tsukada/git/RAFT/projects/test_project',
    '--warm_start',
    ])
#%%
model = torch.nn.DataParallel(RAFT(args))
model.load_state_dict(torch.load(args.model, map_location=torch.device('cpu')))
model = model.module
model.to(DEVICE)
model.eval()

# %%
# Print model's state_dict
print("Model's state_dict:")
for param_tensor in model.state_dict():
    print(param_tensor, "\t", model.state_dict()[param_tensor].size())

# %%
net = "cnet"
weights = np.array(model.state_dict()[f'{net}.conv1.weight'])
root = int(np.sqrt(weights.shape[0]))
fig, axs = plt.subplots(root, root, figsize=(8,8), facecolor="w")
for i in range(weights.shape[0]):
    w = weights[i]
    w = np.transpose(w, (1,2,0))
    wmax = w.max()
    wmin = w.min()
    w = 255 * (w-wmin)/(wmax-wmin)

    x, y = i%root, i//root
    ax = axs[y,x]
    ax.imshow(w.astype(np.uint8))
    ax.axis("off")
fig.savefig(f"{net}_weights.png", bbox_inches="tight", pad_inches=0.1)
# %%
from torch import nn
class Conv(nn.Module):
    def __init__(self):
        super(Conv, self).__init__()
        self.conv1 = nn.Conv2d(3, 64, kernel_size=7, stride=2, padding=3) # このレイヤーで1/2

    def forward(self, x):
        x = self.conv1(x)
        return x

conv1 = Conv()
conv1.load_state_dict({"conv1.weight": model.state_dict()[f"{net}.conv1.weight"], "conv1.bias": model.state_dict()[f"{net}.conv1.bias"]})
conv1.eval()

img = np.array(Image.open("/Users/tsukada/git/RAFT/projects/test_project/frames/test_t0.png")).astype(np.uint8)
if img.shape[-1] >= 4:
    img = img[:,:,:3] # delete alpha cnannel
img = torch.from_numpy(img).permute(2, 0, 1).float()
img = img[None]

# img = torch.tensor(img)

padder = InputPadder(img.shape) # 8の倍数に
img = padder.pad(img)[0]

img = 2 * (img / 255.0) - 1.0
img = img.contiguous()

conved = conv1(img).detach().numpy()[0,:,:,:]

weights = np.array(model.state_dict()['cnet.conv1.weight'])
root = int(np.sqrt(weights.shape[0]))
fig, axs = plt.subplots(root, root, figsize=(8,8), facecolor="w")
for i in range(weights.shape[0]):
    c = conved[i]
    x, y = i%root, i//root
    ax = axs[y,x]
    ax.imshow(c, cmap="gray")
    ax.axis("off")
fig.savefig(f"{net}_conved.png", bbox_inches="tight", pad_inches=0.1)
# %%

max_stars_repo_name: tsukada-cs/RAFT

max_issues_repo_name: tsukada-cs/RAFT

max_forks_repo_name: tsukada-cs/RAFT