text:
# MIT License
#
# Copyright (C) The Adversarial Robustness Toolbox (ART) Authors 2020
#
# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated
# documentation files (the "Software"), to deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit
# persons to whom the Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all copies or substantial portions of the
# Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
# WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
# TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
from __future__ import absolute_import, division, print_function, unicode_literals

import logging
import pytest

import numpy as np
import torch.nn as nn
import torch.optim as optim

from art.attacks.inference import (
    AttributeInferenceWhiteBoxDecisionTree,
    AttributeInferenceBlackBox,
    AttributeInferenceWhiteBoxLifestyleDecisionTree,
)
from art.estimators.classification.pytorch import PyTorchClassifier
from art.estimators.estimator import BaseEstimator
from art.estimators.classification.scikitlearn import ScikitlearnDecisionTreeClassifier

from tests.attacks.utils import backend_test_classifier_type_check_fail

logger = logging.getLogger(__name__)

def test_black_box(get_tabular_classifier_list, get_iris_dataset):
    classifier_list = get_tabular_classifier_list(AttributeInferenceBlackBox)
    if not classifier_list:
        logging.warning("Couldn't perform  this test because no classifier is defined")
        return

    attack_feature = 2  # petal length

    # need to transform attacked feature into categorical
    def transform_feature(x):
        x[x > 0.5] = 2.0
        x[(x > 0.2) & (x <= 0.5)] = 1.0
        x[x <= 0.2] = 0.0

    values = [0.0, 1.0, 2.0]

    (x_train_iris, y_train_iris), (x_test_iris, y_test_iris) = get_iris_dataset
    # training data without attacked feature
    x_train_for_attack = np.delete(x_train_iris, attack_feature, 1)
    # only attacked feature
    x_train_feature = x_train_iris[:, attack_feature].copy().reshape(-1, 1)
    transform_feature(x_train_feature)
    # training data with attacked feature (after transformation)
    x_train = np.concatenate((x_train_for_attack[:, :attack_feature], x_train_feature), axis=1)
    x_train = np.concatenate((x_train, x_train_for_attack[:, attack_feature:]), axis=1)

    # test data without attacked feature
    x_test_for_attack = np.delete(x_test_iris, attack_feature, 1)
    # only attacked feature
    x_test_feature = x_test_iris[:, attack_feature].copy().reshape(-1, 1)
    transform_feature(x_test_feature)

    for classifier in classifier_list:
        # print(type(classifier).__name__)
        if type(classifier).__name__ == "ScikitlearnDecisionTreeClassifier":
            attack = AttributeInferenceBlackBox(classifier, attack_feature=attack_feature)
            # get original model's predictions
            x_train_predictions = np.array([np.argmax(arr) for arr in classifier.predict(x_train_iris)]).reshape(-1, 1)
            x_test_predictions = np.array([np.argmax(arr) for arr in classifier.predict(x_test_iris)]).reshape(-1, 1)
            # train attack model
            attack.fit(x_train)
            # infer attacked feature
            inferred_train = attack.infer(x_train_for_attack, x_train_predictions, values=values)
            inferred_test = attack.infer(x_test_for_attack, x_test_predictions, values=values)
            # check accuracy
            train_acc = np.sum(inferred_train == x_train_feature.reshape(1, -1)) / len(inferred_train)
            test_acc = np.sum(inferred_test == x_test_feature.reshape(1, -1)) / len(inferred_test)
            assert train_acc == pytest.approx(0.8285, abs=0.03)
            assert test_acc == pytest.approx(0.8888, abs=0.03)

def test_black_box_with_model(get_tabular_classifier_list, get_iris_dataset):
    classifier_list = get_tabular_classifier_list(AttributeInferenceBlackBox)
    if not classifier_list:
        logging.warning("Couldn't perform  this test because no classifier is defined")
        return

    attack_feature = 2  # petal length

    # need to transform attacked feature into categorical
    def transform_feature(x):
        x[x > 0.5] = 2.0
        x[(x > 0.2) & (x <= 0.5)] = 1.0
        x[x <= 0.2] = 0.0

    values = [0.0, 1.0, 2.0]

    (x_train_iris, y_train_iris), (x_test_iris, y_test_iris) = get_iris_dataset
    # training data without attacked feature
    x_train_for_attack = np.delete(x_train_iris, attack_feature, 1)
    # only attacked feature
    x_train_feature = x_train_iris[:, attack_feature].copy().reshape(-1, 1)
    transform_feature(x_train_feature)
    # training data with attacked feature (after transformation)
    x_train = np.concatenate((x_train_for_attack[:, :attack_feature], x_train_feature), axis=1)
    x_train = np.concatenate((x_train, x_train_for_attack[:, attack_feature:]), axis=1)

    # test data without attacked feature
    x_test_for_attack = np.delete(x_test_iris, attack_feature, 1)
    # only attacked feature
    x_test_feature = x_test_iris[:, attack_feature].copy().reshape(-1, 1)
    transform_feature(x_test_feature)

    model = nn.Linear(4, 3)

    # Define a loss function and optimizer
    loss_fn = nn.CrossEntropyLoss()
    optimizer = optim.Adam(model.parameters(), lr=0.01)
    attack_model = PyTorchClassifier(
        model=model, clip_values=(0, 1), loss=loss_fn, optimizer=optimizer, input_shape=(4,), nb_classes=3
    )

    for classifier in classifier_list:
        if type(classifier).__name__ == "ScikitlearnDecisionTreeClassifier":
            attack = AttributeInferenceBlackBox(classifier, attack_model=attack_model, attack_feature=attack_feature)
            # get original model's predictions
            x_train_predictions = np.array([np.argmax(arr) for arr in classifier.predict(x_train_iris)]).reshape(-1, 1)
            x_test_predictions = np.array([np.argmax(arr) for arr in classifier.predict(x_test_iris)]).reshape(-1, 1)
            # train attack model
            attack.fit(x_train)
            # infer attacked feature
            inferred_train = attack.infer(x_train_for_attack, x_train_predictions, values=values)
            inferred_test = attack.infer(x_test_for_attack, x_test_predictions, values=values)
            # check accuracy
            train_acc = np.sum(inferred_train == x_train_feature.reshape(1, -1)) / len(inferred_train)
            test_acc = np.sum(inferred_test == x_test_feature.reshape(1, -1)) / len(inferred_test)
            # assert train_acc == pytest.approx(0.5523, abs=0.03)
            # assert test_acc == pytest.approx(0.5777, abs=0.03)

def test_white_box(get_tabular_classifier_list, get_iris_dataset):
    classifier_list = get_tabular_classifier_list(AttributeInferenceWhiteBoxDecisionTree)
    if not classifier_list:
        logging.warning("Couldn't perform  this test because no classifier is defined")
        return

    attack_feature = 2  # petal length
    values = [0.14, 0.42, 0.71]  # rounded down
    priors = [50 / 150, 54 / 150, 46 / 150]

    (x_train_iris, y_train_iris), (x_test_iris, y_test_iris) = get_iris_dataset
    x_train_for_attack = np.delete(x_train_iris, attack_feature, 1)
    x_train_feature = x_train_iris[:, attack_feature]
    x_test_for_attack = np.delete(x_test_iris, attack_feature, 1)
    x_test_feature = x_test_iris[:, attack_feature]

    for classifier in classifier_list:
        if type(classifier).__name__ == "ScikitlearnDecisionTreeClassifier":
            attack = AttributeInferenceWhiteBoxDecisionTree(classifier, attack_feature=attack_feature)
            x_train_predictions = np.array([np.argmax(arr) for arr in classifier.predict(x_train_iris)]).reshape(-1, 1)
            x_test_predictions = np.array([np.argmax(arr) for arr in classifier.predict(x_test_iris)]).reshape(-1, 1)
            inferred_train = attack.infer(x_train_for_attack, x_train_predictions, values=values, priors=priors)
            inferred_test = attack.infer(x_test_for_attack, x_test_predictions, values=values, priors=priors)
            train_diff = np.abs(inferred_train - x_train_feature.reshape(1, -1))
            test_diff = np.abs(inferred_test - x_test_feature.reshape(1, -1))
            assert np.sum(train_diff) / len(inferred_train) == pytest.approx(0.2108, abs=0.03)
            assert np.sum(test_diff) / len(inferred_test) == pytest.approx(0.1988, abs=0.03)

def test_white_box_lifestyle(get_tabular_classifier_list, get_iris_dataset):
    classifier_list = get_tabular_classifier_list(AttributeInferenceWhiteBoxLifestyleDecisionTree)
    if not classifier_list:
        logging.warning("Couldn't perform  this test because no classifier is defined")
        return

    attack_feature = 2  # petal length
    values = [0.14, 0.42, 0.71]  # rounded down
    priors = [50 / 150, 54 / 150, 46 / 150]

    (x_train_iris, y_train_iris), (x_test_iris, y_test_iris) = get_iris_dataset
    x_train_for_attack = np.delete(x_train_iris, attack_feature, 1)
    x_train_feature = x_train_iris[:, attack_feature]
    x_test_for_attack = np.delete(x_test_iris, attack_feature, 1)
    x_test_feature = x_test_iris[:, attack_feature]

    for classifier in classifier_list:
        if type(classifier).__name__ == "ScikitlearnDecisionTreeClassifier":
            attack = AttributeInferenceWhiteBoxLifestyleDecisionTree(classifier, attack_feature=attack_feature)
            x_train_predictions = np.array([np.argmax(arr) for arr in classifier.predict(x_train_iris)]).reshape(-1, 1)
            x_test_predictions = np.array([np.argmax(arr) for arr in classifier.predict(x_test_iris)]).reshape(-1, 1)
            inferred_train = attack.infer(x_train_for_attack, x_train_predictions, values=values, priors=priors)
            inferred_test = attack.infer(x_test_for_attack, x_test_predictions, values=values, priors=priors)
            train_diff = np.abs(inferred_train - x_train_feature.reshape(1, -1))
            test_diff = np.abs(inferred_test - x_test_feature.reshape(1, -1))
            assert np.sum(train_diff) / len(inferred_train) == pytest.approx(0.3357, abs=0.03)
            assert np.sum(test_diff) / len(inferred_test) == pytest.approx(0.3149, abs=0.03)
            # assert np.sum(train_diff) / len(inferred_train) < np.sum(test_diff) / len(inferred_test)

def test_classifier_type_check_fail():
    backend_test_classifier_type_check_fail(AttributeInferenceBlackBox, [BaseEstimator])
    backend_test_classifier_type_check_fail(
        AttributeInferenceWhiteBoxLifestyleDecisionTree, [BaseEstimator, ScikitlearnDecisionTreeClassifier]
    )
    backend_test_classifier_type_check_fail(
        AttributeInferenceWhiteBoxDecisionTree, [BaseEstimator, ScikitlearnDecisionTreeClassifier]
    )

if __name__ == "__main__":
    pytest.cmdline.main("-q -s {} --mlFramework=scikitlearn --durations=0".format(__file__).split(" "))


text:
# ===============================================================================#
import re, sys
import os
import csv
import numpy as np
from datetime import datetime
import itertools
from collections import namedtuple, defaultdict, OrderedDict, deque
import six

if six.PY3:
    from configparser import ConfigParser
elif six.PY2:
    from ConfigParser import ConfigParser
import logging
import json
from collections import defaultdict

try:
    from PIL import Image, ImageFont, ImageDraw

    # imagetitle = True
except ImportError:
    pass
# imagetitle = False

RefseqInfo = namedtuple(
    "RefseqInfo", "taxonid, geneid, homologeneid,proteingi,genefraglen"
)

def byte_formatter(b):
    conv = b / (2 ** 10)
    if conv < 1000:
        return "{:.4f} KB".format(conv)
    elif conv > 1000:
        conv = conv / (2 ** 10)
        if conv < 1000:
            return "{:.4f} MB".format(conv)
        elif conv < 1000:
            conv = conv / (2 ** 10)
            return "{:.4f} GB".format(conv)

def bufcount(filename):
    """fast way to count lines in a file"""
    f = open(filename, mode="rt", encoding="latin1")
    lines = 0
    buf_size = 2048
    read_f = f.read  # loop optimization

    buf = read_f(buf_size)
    while buf:
        lines += buf.count("\n")
        buf = read_f(buf_size)

    return lines + 1  # plus 1 since I use csv_reader, which ignores the header

def rolling_window(seq, length_of_window):
    it = iter(seq)
    window = deque(maxlen=length_of_window)
    for _ in range(length_of_window):
        try:
            window.append(next(it))
        except StopIteration:
            return
    yield tuple(window)
    while it:
        try:
            _next = next(it)
            # window.append(next(it)) # this is wrong!
            window.append(_next)  # this is correct
            yield tuple(window)
        # except (StopIteration, RuntimeError):
        except StopIteration:
            return

from pyteomics import parser

def protease(seq, miscuts=2, rule="[KR]", min_pept_len=7, semi=False, **kwargs):
    nomiscuts = parser.cleave(
        seq, rule=rule, semi=semi, missed_cleavages=0, min_length=min_pept_len
    )
    all_seqs = parser.cleave(
        seq, rule=rule, semi=semi, missed_cleavages=miscuts, min_length=min_pept_len
    )

    # without M
    nomiscuts_without_m = parser.cleave(
        seq[1:], rule=rule, missed_cleavages=0, min_length=min_pept_len, semi=semi
    )
    all_seqs_without_m = parser.cleave(
        seq[1:], rule=rule, missed_cleavages=miscuts, min_length=min_pept_len, semi=semi
    )

    nomiscuts = nomiscuts | nomiscuts_without_m
    all_seqs = all_seqs | all_seqs_without_m
    return (all_seqs, len(nomiscuts))

def protease_x(
    seq,
    minlen=0,
    cutsites=tuple(),
    exceptions=None,
    miscuts=2,
    semi_tryptic=False,
    semi_tryptic_iter=6,
):
    frags = []
    chop = ""
    if exceptions is None:
        exceptions = tuple()
    import ipdb

    ipdb.set_trace()
    while seq:
        cuts = [seq.find(x) for x in cutsites]
        if len(cuts) > 0:
            if min(cuts) >= 0:
                cut = min(cuts) + 1
            else:
                cut = max(cuts) + 1
        else:
            cut = 0
        chop += seq[:cut]
        seq = seq[cut:]
        if chop == "LTEAKPAKPAPAK":
            import ipdb

            ipdb.set_trace()
        if cut == 0 or len(seq) == 0:
            if cut == 0 and len(frags) == 0:
                frags.append(chop + seq)
            elif cut == 0 and chop and seq:
                frags.append(chop + seq)
            # print(chop)
            elif len(seq) != 0:
                frags.append(seq)
            elif len(seq) == 0 and chop:
                frags.append(chop)  # special case for
            # if last amino acid is one of the cut sites
            break  # no more cuts
        if seq[0] not in exceptions:
            frags.append(chop)
            chop = ""
    merged_list = list()
    no_met = list()
    for k in range(0, miscuts):
        for ix, chunk in enumerate(rolling_window(frags, k + 2)):  # +2 to join adjacent
            merged_list.append("".join(chunk))
            if ix == 0:
                no_met = list()
                for ix, c in enumerate(chunk):
                    if ix == 0:
                        no_met.append(c[1:])
                    else:
                        no_met.append(c)
        merged_list.append("".join(no_met))

    if semi_tryptic:  #
        # MAX_ITER = 6 # let's not go beyond this, too many peptides
        MAX_ITER = semi_tryptic_iter
        # maybe have option later
        for frag in frags:
            for i in range(1, len(frag)):
                if i > MAX_ITER:
                    break
                semiN = frag[i:]
                semiC = frag[: -1 * i]
                if len(semiN) < minlen and len(semiC) < minlen:
                    break
                if len(semiN) >= minlen and semiN not in merged_list:
                    merged_list.append(semiN)
                if len(semiC) >= minlen:
                    merged_list.append(semiC)

    if len(frags) > 1:
        frags.append(frags[0][1:])  # chop off methionine
    nomiscuts_len = len([x for x in frags if len(x) >= minlen])

    return [x for x in frags + merged_list if len(x) >= minlen], nomiscuts_len

def fancyprint(ShowText, string_size=12):
    try:
        font = ImageFont.truetype("arialbd.ttf", string_size)  # load the font
        font_import = True
    except IOError:
        font_import = False
        print(ShowText)
    if font_import:
        size = font.getsize(ShowText)  # calc the size of text in pixels
        image = Image.new("1", size, 1)  # create a b/w image
        draw = ImageDraw.Draw(image)
        draw.text((0, 0), ShowText, font=font)  # render the text to the bitmap
        for rownum in range(size[1]):
            # scan the bitmap:
            # print ' ' for black pixel and
            # print '#' for white one
            line = []
            for colnum in range(size[0]):
                if image.getpixel((colnum, rownum)):
                    line.append(" "),
                else:
                    line.append("#"),
            print("".join(line))

label_pos = re.compile(
    """
    (?<!Label\:)  # not preceded by
    (\d+|N-Term)      # Match numbers,
    (?!       # only if it's not followed by..."
     [^(]*    #   any number of characters except opening parens" +
     \\)      #   followed by a closing parens" +
    )         # End of lookahead""",
    re.VERBOSE,
)

inside_paren = re.compile(r"\((\S+)\)")
amino_acids = list("ACDEFGHIKLMNPQRSTVWY")

def seq_modi(sequence, modifications, no_count=None):
    """
    function to output a modified string of sequence that includes all modifications at the appropriate place.
    Modified amino acids are indicated in the input 'sequence' by a lowercase letter.
    A dictionary of potential modification keys and their corresponding substitutions in the primary sequence is provided;
    all other modifications should be ignored.
    Sometimes N-terminal modifications are present, this should be considered the same as a modification at the first amino acid.

    modifications is a 1-based index of amino acid and modification delimited by semicolon (; ):
    M1(Oxidation) is first amino acid methionine Oxidized

    """
    if no_count is None:
        no_count = tuple()
    # TODO : have separate collections for labeled and non-labeled modifications
    seqmodi = ""
    seqlist = list(sequence)
    modi_len = 0  # default length is zero, can change
    label = 0  # default label = 0, can change
    label_search = re.search("(Label\S+)(?=\))", modifications)  # for SILAC
    if label_search:
        label = label_search.group(0)

    # if not any(c for c in sequence if (c.islower() or c == 'X')): # check if any modifications to deal with
    if (
        not any(c == "X" for c in sequence) and not modifications
    ):  # check if any modifications to deal with
        return sequence, sequence, 0, label
    modkeys = inside_paren.findall(modifications)
    modi_len = len([x for x in modkeys if not any(y in x for y in no_count)])
    # any(mod in x for mod in to_count)])
    # modkeys = re.findall(r'(\([^\)]+\))',modifications)  #get all modifications
    # modpos = re.findall(r'[A-Z]([0-9]+)',modifications)  # get all modification positions
    modpos = label_pos.findall(modifications)
    modpos = [
        0 if x.lower() == "n-term" else int(x) - 1 for x in modpos
    ]  # n terminal is at first position
    mod_dict = defaultdict(list)
    for (key, value) in zip(modpos, modkeys):
        mod_dict[key].append(value.lower())
    for key, values in mod_dict.items():
        mod_dict[key] = sorted(values)
    for ix, s in enumerate(sequence.upper()):
        # print(ix, seqmodi)
        if ix in mod_dict:
            if s == "X":  # deal with the X first
                to_replace = [x for x in mod_dict[ix] if x in amino_acids]
                if len(to_replace) == 1:
                    if seqlist[ix].islower():
                        replaced_aa = to_replace[0].lower()
                    else:
                        replaced_aa = to_replace[0]
                    seqmodi += replaced_aa.upper()  # all seqmodi AAs are upper case
                    seqlist[ix] = replaced_aa
                    mod_dict[ix].remove(replaced_aa.upper())
                elif len(to_replace) == 0:
                    pass  # not an amino acid (listed above at least)
                    # Probably an unidentified mass of the form X10(110.0)
                else:
                    print(
                        "Error parsing sequence {} with modifications {}".format(
                            sequence, modifications
                        )
                    )
                    seqmodi += s
            else:
                seqmodi += s
        else:
            seqmodi += s
        if ix not in mod_dict:
            continue
        for modi in mod_dict[ix]:
            if sequence[ix].islower() and modi not in amino_acids:
                modi_ = modi.lower()
            else:
                modi_ = modi
            # seqmodi += modi
            seqmodi += "({})".format(modi_)
    sequence = "".join(seqlist)
    return sequence, seqmodi, modi_len, label
    # return sequence, seqmodi, modi_len

def _seq_modi_old(sequence, modifications, keeplog=True):

    """
    function to output a modified string of sequence that includes all
    modifications at the appropriate place.
    Modified amino acids are indicated in the input 'sequence' by a lowercase
    letter.
    A dictionary of potential modification keys and their corresponding
    substitutions in the primary sequence is provided;
    all other modifications should be ignored.
    Sometimes N-terminal modifications are present, this should be considered
    the same as a modification at the first amino acid.

    also modifies sequence to fill in Xs with the predicted modification
    example : AAAx(P)R --> AAAPR
    """
    amino_acids = list("ACDEFGHIKLMNPQRSTVWY")
    amino_acids = ["(" + x + ")" for x in amino_acids]
    modtext = {
        "(DeStreak)": "(des)",
        "(Deamidated)": "(dam)",
        "(Carbamidomethyl)": "(car)",
        "(Oxidation)": "(oxi)",
        "(Phospho)": "(pho)",
        "(Prot)(Acetyl)": "(ace)",
        "(Acetyl)": "(ace)",
        "(GlyGly)": "(gg)",
        "(Label:13C(6))": "(lab)",
        "(Label:13C(6)+GlyGly)": "(labgg)",
    }
    seqmodi = ""
    seqlist = list(sequence)
    modi_len = 0  # default length is zero, can change
    label = 0  # default label = 0, can change
    if any(c for c in sequence if (c.islower() or c == "X")):
        # check if any modifications to deal with
        if "Label" in modifications:
            # routine for grabbing everything enclosed in parenthesis, including
            # nested. Easier to write than regex when dealing with nested
            label = 1  # for indication of label on peptide
            openb, modi, modkeys = 0, "", []
            for e in modifications:
                if e == "(":
                    openb += 1
                elif e == ")":
                    openb += -1
                if openb:
                    modi += e
                elif not openb:
                    if modi:
                        modi += e  # add the last parenthesis
                        modkeys.append(modi)
                    modi = ""
        else:
            modkeys = re.findall(r"(\([^\)]+\))", modifications)
            # get all modifications

        modpos = re.findall(r"[A-Z]([0-9]+)", modifications)
        # get all modification positions
        modpos = [int(d) - 1 for d in modpos]  # convert all to integers
        if "N-Term" in modifications:
            modpos.insert(0, 0)
            # first of the modkeys will be the N-Term modification
        modi_len = len(modpos)
        mod_dict = defaultdict(list)
        for (key, value) in zip(modpos, modkeys):
            mod_dict[key].append(value)
        for key in mod_dict:
            mod_dict[key] = sorted(mod_dict[key])

        for ix, s in enumerate(sequence.upper()):
            seqmodi += s
            # print(ix, seqmodi)
            if ix in mod_dict:
                if s == "X":  # deal with the X first
                    to_replace = [x for x in mod_dict[ix] if x in amino_acids]

                    if len(to_replace) == 1:
                        # print(seqlist[ix],to_replace)
                        if seqlist[ix].islower():
                            seqlist[ix] = to_replace[0][1].lower()
                        else:
                            seqlist[ix] = to_replace[0][1]
                        modi_len += -1  # don't count X == amino acid as a
                        # modification since technically it is not a PTM
                        # mod_dict[ix].remove(to_replace)  # not sure if we want
                        # to remove this or not, but we certainly can if we want
                    elif len(to_replace) == 0:
                        pass  # not an amino acid (listed above at least)
                    # Probably an unidentified mass of the form X10(110.0)

                    else:
                        print(
                            "Error parsing sequence {} with "
                            "modifications {}".format(sequence, modifications)
                        )
                for modi in mod_dict[ix]:
                    if modi in modtext:
                        seqmodi += modtext[modi]
                    elif modi in amino_acids:
                        seqmodi += modi  # for the X amino acids
                    else:
                        if keeplog:
                            logging.warning(
                                "New modification {} that"
                                " is not found in sequence {}".format(modi, sequence)
                            )

    sequence = "".join(seqlist)
    if not seqmodi:
        seqmodi = sequence
    return sequence, seqmodi, modi_len, label

def _count_modis_maxquant(modi, labeltype):
    if modi == "Unmodified":
        return 0
    count = modi.count(",") + 1  # multiple modis are separated by a comma
    if labeltype == "TMT":
        count -= modi.lower().count("tmt")
    elif labeltype == "iTRAQ":
        count -= modi.lower().count("itraq")
    return count

def count_modis_maxquant(df, labeltype):
    return df.apply(
        lambda x: _count_modis_maxquant(x["Modifications"], labeltype), axis=1
    )

def count_modis_seqmodi(df, labeltype):
    if labeltype == "TMT":
        return df.SequenceModi.str.count("\(") - df.SequenceModi.str.count(
            "\([TMT|tmt]"
        )
    elif labeltype == "iTRAQ":
        return df.SequenceModi.str.count("\(") - df.SequenceModi.str.count(
            "\([iTRAQ|itraq]"
        )
    else:
        return df.SequenceModi.str.count("\(")

def calculate_miscuts(seq, targets=None, exceptions=None):
    """Calculates number of miscuts for a given sequence
    using amino acids given in targets"""
    if exceptions is None:
        exceptions = tuple()
    not_miscut = ["".join(x) for x in itertools.product(targets, exceptions)]
    miscuts = sum(seq.count(x) for x in targets)
    miscuts_not = sum(seq.count(x) for x in not_miscut)
    miscuts -= miscuts_not
    if not any(seq[-1] == x for x in targets):  # then at the C terminal
        return miscuts
    return miscuts - 1

import hashlib

def md5sum(filename, blocksize=65536):
    hash = hashlib.md5()
    with open(filename, "rb") as f:
        for block in iter(lambda: f.read(blocksize), b""):
            hash.update(block)
    return hash.hexdigest()

def write_md5(filename, chsum):
    with open(filename, "w") as f:
        f.write(chsum)

def _fasta_dict_from_file(file_object, header_search="specific"):
    """
    Reads a file of FASTA entries and returns a dict for each entry.
    It parsers the headers such that `>gi|12345|ref|NP_XXXXX| Description`
    returns as `{gi: 12345, ref: 'NP_XXXXX', description : 'Description'}`
    The sequence has all whitespace removed.

    :header_search: One of `specific` or `generic`
    if `specific` tries to parse the header
    if `general` just returns each whitespace separated header
    """

    current_id = dict()
    current_seq = ""
    current_header = None
    pat = re.compile(">(\S+)\s*(.*)")
    # header_pat = re.compile(r'(\w+)\|(\w+\.?\w*)?')
    header_pat = re.compile(r"(\w+)\|(\S+)?")
    headers = ["gi", "ref", "geneid", "homologene", "taxon", "symbol"]
    header_pat = re.compile(
        r"gi\|(\S*)\|ref\|(\S*)\|geneid\|(\S*)\|homologene\|(\S*)\|taxon\|(\S*)\|symbol\|(\S*)"
    )

    def parse_header(header, pairs=True):
        header_dict = dict()
        for h in headers:
            try:
                head_value = re.search(h + "\|(\w+\.?\d?)\|", header).group(1)
            except AttributeError:
                head_value = "" if h != "geneid" else header
            header_dict[h] = head_value

        # try:
        #     keys = header_pat.findall(header)[0]
        # except IndexError:
        #     # warnings.warn('Error parsing {}, using full as gene id'.format(header))
        #     return dict(gi='', ref='', geneid=header, homologene='', taxon='', symbol='')
        # header_data = dict()
        # for h, key in zip(headers, keys):
        #     header_data[h] = key
        # for key in keys:
        #     header_data[key[0]] = key[1]
        # gi -> ProteinGI #, ref -> NP_XXXX
        return header_dict

    for line in file_object:
        line = line.rstrip()
        m = pat.search(line)
        if m:
            ## new residue line matched, purge the existing one, if not the first
            if current_id:
                ## remove all whitespace and save
                current_seq = "".join(current_seq.split())
                current_id["sequence"] = current_seq
                yield current_id
                # current_id.clear()  # this is actually bad for list comprehensions
                # as it returns empty dictionaries

            current_seq = ""
            header = m.group(1)
            if header_search == "specific":
                current_id = parse_header(header)
            elif header_search == "generic":
                current_id = dict(header=header)
            current_id["description"] = m.group(2)

        else:
            ## python 2.6+ makes string concatenation amortized O(n)
            ##  http://stackoverflow.com/a/4435752/1368079
            current_seq += str(line)

    ## don't forget the last one
    current_seq = "".join(current_seq.split())
    current_id["sequence"] = current_seq
    yield current_id

def fasta_dict_from_file(fasta_file, header_search="specific"):
    with open(fasta_file, "r") as f:
        # yield from _fasta_dict_from_file(f, header_search=header_search)
        for v in _fasta_dict_from_file(f, header_search=header_search):
            yield v

fasta_dict_from_file.__doc__ = _fasta_dict_from_file.__doc__

def convert_tab_to_fasta(
    tabfile,
    geneid=None,
    ref=None,
    gi=None,
    homologene=None,
    taxon=None,
    description=None,
    sequence=None,
    symbol=None,
):
    HEADERS = {
        "geneid": geneid,
        "ref": ref,
        "gi": gi,
        "homologene": homologene,
        "taxon": taxon,
        "description": description,
        "sequence": sequence,
        "symbol": symbol,
    }
    CUTOFF = 0.35
    df = pd.read_table(tabfile, dtype=str)
    choices = click.Choice([x for y in [df.columns, ["SKIP"]] for x in y])
    col_names = dict()
    for h, v in HEADERS.items():
        if v is not None or v == "SKIP":
            col_names[h] = v
            continue

        closest_match = difflib.get_close_matches(h, df.columns, n=1, cutoff=CUTOFF)
        if closest_match and h != "homologene":
            col_names[h] = closest_match[0]
        else:
            print("Can not find header match for :", h)
            print("Choose from :", " ".join(choices.choices))
            choice = click.prompt(
                "Selected the correct name or SKIP",
                type=choices,
                show_default=True,
                err=True,
                default="SKIP",
            )
            col_names[h] = choice
            print()
    for _, series in df.iterrows():
        row = series.to_dict()
        gid = row.get(col_names["geneid"], "")
        ref = row.get(col_names["ref"], "")
        hid = row.get(col_names["homologene"], "")
        gi = row.get(col_names["gi"], "")
        taxon = row.get(col_names["taxon"], "")
        desc = row.get(col_names["description"], " ")
        seq = "\n".join(textwrap.wrap(row.get(col_names["sequence"], ""), width=70))
        symbol = row.get(col_names["symbol"], "")

        hid = int(hid) if hid and not np.isnan(hid) else ""
        try:
            gid = int(gid)
        except ValueError:
            pass

        r = dict(
            gid=gid,
            ref=ref,
            taxons=taxon,
            gis=gi,
            homologene=hid,
            description=desc,
            seq=seq,
            symbols=symbol,
        )
        yield (FASTA_FMT.format(**r))

def sniff_fasta(fasta):
    nrecs = 1
    fasta = fasta_dict_from_file(fasta)
    counter = 0
    REQUIRED = ("ref", "sequence")
    while counter < nrecs:
        for rec in fasta:
            if any(x not in rec for x in REQUIRED):
                raise ValueError("Invalid input FASTA")
        counter += 1
    return 0  # all good

ext: py

max_stars_repo_name: asalt/gpgrouper

max_issues_repo_name: asalt/gpgrouper

max_forks_repo_name: asalt/gpgrouper

reason: import numpy


text:
"""Load breakthrough profiles.

Load breakthrough profiles determine what parts of load bind to the
column.

See Also
--------
:class:`bio_rtd.core.ChromatographyLoadBreakthrough`

"""

__all__ = ['ConstantPatternSolution']
__version__ = '0.7.1'
__author__ = 'Jure Sencar'

import numpy as _np

from bio_rtd import core as _core
from bio_rtd.chromatography import bt_curve as _bt_curve

class ConstantPatternSolution(_core.ChromatographyLoadBreakthrough):
    """Breakthrough profile - Constant Pattern Solution.

    r = 1 / (1 + exp(`k` * (`dbc_100` - `m_load` / `cv`)))

    Parameters
    ----------
    dt
        Time step duration.
    dbc_100
        Dynamic binding capacity if the load would last indefinitely.
    k
        Steepness of the breakthrough profile.
    bt_profile_id
         Unique identifier of the PDF instance. Used for logs.

    Examples
    --------
    >>> t = _np.linspace(0, 120, 1001)  # min
    >>> dt = t[1]
    >>> btc = ConstantPatternSolution(dt, dbc_100=240, k=0.05)
    >>> btc.update_btc_parameters(cv=8.0)  # column volume [mL]
    >>> btc.get_total_bc()  # 8.0 [mL] * 240 [mg/mL] = 1920.0 [mg]
    1920.0
    >>> # Load.
    >>> f_in = _np.ones_like(t)  # mL / min
    >>> c_in = _np.ones([1, t.size]) * 15  # mg / mL
    >>> c_captured = btc.calc_c_bound(f_in, c_in)
    >>> c_captured.shape == c_in.shape
    True
    >>> c_captured
    array([[14.99990679, 14.99990574, 14.99990467, ..., 10.22437591,
            10.18768049, 10.15083683]])
    >>> round((c_captured * f_in * dt).sum(), 1)  # < 1920.0
    1739.0

    """
    POSSIBLE_KEY_GROUPS = [['cv']]
    OPTIONAL_KEYS = []

    def __init__(self, dt, dbc_100: float, k: float,
                 bt_profile_id: str = "ConstantPatternSolution"):
        super().__init__(dt, bt_profile_id)
        self.dbc_100 = dbc_100
        """Dynamic binding capacity at 100 % load breakthrough."""
        self.k = k
        """Steepness of the breakthrough profile."""
        # placeholder for column volume
        self._cv = -1

    def _update_btc_parameters(self, kw_pars: dict) -> None:

        self._cv = kw_pars["cv"]

    def _calc_unbound_to_load_ratio(
            self,
            m_cumulative_load: _np.ndarray) -> _np.ndarray:

        assert self._cv > 0, f"CV must be defined by now." \
                             f" Make sure `update_btc_parameters`" \
                             f" function was called before."

        ratio = _bt_curve.btc_constant_pattern_solution(m_cumulative_load,
                                                        self.dbc_100,
                                                        self.k,
                                                        self._cv,
                                                        self.log)

        return ratio

    def get_total_bc(self) -> float:
        assert self._cv > 0
        return self.dbc_100 * self._cv

max_stars_repo_name: open-biotech/bio-rtd

max_issues_repo_name: open-biotech/bio-rtd

max_forks_repo_name: open-biotech/bio-rtd


text:
#!/usr/bin/env python
"""
Test one data file
"""

from numpy import *
import sys

import kernel
import distance
import classifier
import clustering
import distribution
import regression
import preprocessor
from sg import sg

SUPPORTED=['kernel', 'distance', 'classifier', 'clustering', 'distribution',
	'regression', 'preprocessor']

def _get_name_fun (fnam):
	module=None

	for supported in SUPPORTED:
		if fnam.find(supported)>-1:
			module=supported
			break

	if module is None:
		print 'Module required for %s not supported yet!' % fnam
		return None

	return module+'.test'

def _test_mfile (fnam):
	try:
		mfile=open(fnam, mode='r')
	except IOError, e:
		print e
		return False

	indata={}

	name_fun=_get_name_fun(fnam)
	if name_fun is None:
		return False

	for line in mfile:
		line=line.strip(" \t\n;")
		param = line.split('=')[0].strip()

		if param=='name':
			name=line.split('=')[1].strip().split("'")[1]
			indata[param]=name
		elif param=='kernel_symdata' or param=='kernel_data':
			indata[param]=_read_matrix(line)
		elif param.startswith('kernel_matrix') or \
			param.startswith('distance_matrix'):
			indata[param]=_read_matrix(line)
		elif param.find('data_train')>-1 or param.find('data_test')>-1:
			# data_{train,test} might also be prepended by *subkernel*
			indata[param]=_read_matrix(line)
		elif param=='clustering_centers' or param=='clustering_pairs':
			indata[param]=_read_matrix(line)
		else:
			if (line.find("'")==-1):
				indata[param]=eval(line.split('=')[1])
			else:
				indata[param]=line.split('=')[1].strip().split("'")[1]

	mfile.close()
	fun=eval(name_fun)

	# seed random to constant value used at data file's creation
	sg('init_random', indata['init_random'])
	random.seed(indata['init_random'])

	return fun(indata)

def _read_matrix (line):
	try:
		str_line=(line.split('[')[1]).split(']')[0]
	except IndexError:
		str_line=(line.split('{')[1]).split('}')[0]

	lines=str_line.split(';')
	lis2d=list()

	for x in lines:
		lis=list()
		for y in x.split(','):
			y=y.replace("'","").strip()
			if(y.isalpha()):
				lis.append(y)
			else:
				if y.find('.')!=-1:
					lis.append(float(y))
				else:
					try:
						lis.append(int(y))
					except ValueError: # not int, RAWDNA?
						lis.append(y)
		lis2d.append(lis)

	return array(lis2d)

for filename in sys.argv:
	if (filename.endswith('.m')):
		res=_test_mfile(filename)
		if res:
			sys.exit(0)
		else:
			sys.exit(1)

max_stars_repo_name: srgnuclear/shogun

max_stars_repo_licenses: Ruby MIT

max_issues_repo_name: waderly/shogun

max_issues_repo_licenses: Ruby MIT

max_forks_repo_name: waderly/shogun

max_forks_repo_licenses: Ruby MIT

reason: from numpy


text:
import sys
import torch
import argparse
import numpy as np
import torch.nn as nn
import torch.nn.functional as F
import shutil

from torch.autograd import Variable
from torch.utils import data
import os

from dataset import IC19TestLoader
from util.misc import get_mean_and_std

#data_loader = IC19Loader(is_transform=False, img_size=640, kernel_num=7, min_scale=0.4)
data_loader = IC19TestLoader(data_dirs="/home/azhar/kuzushiji/data_root/train_images/")
mean, std = get_mean_and_std(data_loader)
print(mean,std)

max_stars_repo_name: prabhatrmishra/IDCardInfoExtr

max_issues_repo_name: prabhatrmishra/IDCardInfoExtr

max_forks_repo_name: prabhatrmishra/IDCardInfoExtr


text:
"""
--- Day 1: The Tyranny of the Rocket Equation ---
https://adventofcode.com/2019/day/1
"""
from aocd import data
import numpy as np

A = np.fromstring(data, dtype=int, sep="\n")
total = 0
while A.any():
    A = (A // 3 - 2).clip(0)
    if not total:
        print("part a:", A.sum())
    total += A.sum()
print("part b:", total)

max_stars_repo_name: wimglenn/advent-of-code-wim

max_issues_repo_name: wimglenn/advent-of-code-wim

max_forks_repo_name: wimglenn/advent-of-code-wim


text:
from acme import types
import d4rl
import gym
import numpy as np
from tqdm import tqdm

def split_into_trajectories(
    observations, actions, rewards, masks, dones_float, next_observations
):
    trajs = [[]]

    for i in tqdm(range(len(observations))):
        trajs[-1].append(
            (
                observations[i],
                actions[i],
                rewards[i],
                masks[i],
                dones_float[i],
                next_observations[i],
            )
        )
        if dones_float[i] == 1.0 and i + 1 < len(observations):
            trajs.append([])

    return trajs

def merge_trajectories(trajs):
    observations = []
    actions = []
    rewards = []
    masks = []
    dones_float = []
    next_observations = []

    for traj in trajs:
        for (obs, act, rew, mask, done, next_obs) in traj:
            observations.append(obs)
            actions.append(act)
            rewards.append(rew)
            masks.append(mask)
            dones_float.append(done)
            next_observations.append(next_obs)

    return (
        np.stack(observations),
        np.stack(actions),
        np.stack(rewards),
        np.stack(masks),
        np.stack(dones_float),
        np.stack(next_observations),
    )

class Dataset(object):
    def __init__(
        self,
        observations: np.ndarray,
        actions: np.ndarray,
        rewards: np.ndarray,
        masks: np.ndarray,
        dones_float: np.ndarray,
        next_observations: np.ndarray,
        size: int,
    ):
        self.observations = observations
        self.actions = actions
        self.rewards = rewards
        self.masks = masks
        self.dones_float = dones_float
        self.next_observations = next_observations
        self.size = size

    def sample(self, batch_size: int) -> types.Transition:
        indx = np.random.randint(self.size, size=batch_size)
        return types.Transition(
            observation=self.observations[indx],
            action=self.actions[indx],
            reward=self.rewards[indx],
            discount=self.masks[indx],
            next_observation=self.next_observations[indx],
        )

class D4RLDataset(Dataset):
    def __init__(self, env: gym.Env, clip_to_eps: bool = True, eps: float = 1e-5):
        dataset = d4rl.qlearning_dataset(env)

        if clip_to_eps:
            lim = 1 - eps
            dataset["actions"] = np.clip(dataset["actions"], -lim, lim)

        dones_float = np.zeros_like(dataset["rewards"])

        for i in range(len(dones_float) - 1):
            if (
                np.linalg.norm(
                    dataset["observations"][i + 1] - dataset["next_observations"][i]
                )
                > 1e-6
                or dataset["terminals"][i] == 1.0
            ):
                dones_float[i] = 1
            else:
                dones_float[i] = 0

        dones_float[-1] = 1

        super().__init__(
            dataset["observations"].astype(np.float32),
            actions=dataset["actions"].astype(np.float32),
            rewards=dataset["rewards"].astype(np.float32),
            masks=1.0 - dataset["terminals"].astype(np.float32),
            dones_float=dones_float.astype(np.float32),
            next_observations=dataset["next_observations"].astype(np.float32),
            size=len(dataset["observations"]),
        )

max_stars_repo_name: akbir/magi

max_issues_repo_name: akbir/magi

max_forks_repo_name: akbir/magi


text:
import torch as tn
import torchtt as tntt
import matplotlib.pyplot as plt
import tt_iga
import numpy as np
import datetime
import matplotlib.colors
import pandas as pd

tn.set_default_dtype(tn.float64)

def solve(Np, Ns = [40,20,80], deg = 2 ,nl = 8):
    print()
    print('number of parameters ',Np)
    # used to return the results
    dct = {}

    # B-splines
    Ns = np.array([40,20,80])-deg+1
    baza1 = tt_iga.bspline.BSplineBasis(np.linspace(0,1,Ns[0]),deg)
    baza2 = tt_iga.bspline.BSplineBasis(np.linspace(0,1,Ns[1]),deg)
    baza3 = tt_iga.bspline.BSplineBasis(np.linspace(0,1,Ns[2]),deg)
    Basis = [baza1,baza2,baza3]
    N = [baza1.N,baza2.N,baza3.N]

    # Parameter space basis
    var = 0.05
    Basis_param = [tt_iga.lagrange.LagrangeLeg(nl,[-var,var])]*Np

    # B-spline basis for the radius perturbation
    bspl = tt_iga.bspline.BSplineBasis(np.linspace(0,1,Np-2+3),2)
    def interface_func(t1,tp):
        return tn.einsum('ij,ji->j',tn.tensor(bspl(t1)[1:-1,:]),tn.tensor(tp))
    line = lambda t,a,b: t*(b-a)+a
    damp = lambda x: 1 # -4*x*(x-1)

    # parametrization
    w = 1
    h = 0.5
    r = 2
    xparam = lambda t : (2-h+line(t[:,1],0,interface_func(t[:,2],t[:,3:])*damp(t[:,2])+h))*tn.cos(1.5*np.pi+0.25*np.pi*t[:,2])
    yparam = lambda t : (2-h+line(t[:,1],0,interface_func(t[:,2],t[:,3:])*damp(t[:,2])+h))*tn.sin(1.5*np.pi+0.25*np.pi*t[:,2])
    zparam = lambda t : w*t[:,0]

    # instantiate the GeometryMapping object. It is used for intepolating, evaluating and computing the discrete operators corresponding to a parameter dependent geometry
    geom = tt_iga.Geometry(Basis+Basis_param)
    # interpolate the geometry parametrization
    geom.interpolate([xparam, yparam, zparam])

    # compute the mass matrix in TT
    tme = datetime.datetime.now()
    Mass_tt = geom.mass_interp(eps=1e-11)
    tme = datetime.datetime.now() -tme
    print('Time mass matrix ',tme.total_seconds())

    # if tn.cuda.is_available():
    #     tme = datetime.datetime.now()
    #     Stt = geom.stiffness_interp( eps = 1e-9, qtt = True, verb=True, device = tn.device('cuda:0'))
    #     tme = datetime.datetime.now() -tme
    #     print('Time stiffness matrix GPU',tme.total_seconds())
    #     dct['time stiff GPU'] = tme.total_seconds()

    tme = datetime.datetime.now()
    Stt = geom.stiffness_interp( eps = 1e-9, qtt = True, verb=True, device = None)
    tme = datetime.datetime.now() -tme
    print('Time stiffness matrix ',tme.total_seconds())
    dct['time stiff'] = tme.total_seconds()

    # projection operators for enforcing the BCs
    Pin_tt, Pbd_tt = tt_iga.projectors.get_projectors(N,[[1,1],[0,0],[1,1]])
    Pin_tt = Pin_tt ** tntt.eye([nl]*Np)
    Pbd_tt = Pbd_tt ** tntt.eye([nl]*Np)

    # right hand side. Zero since we solve the homogenous equation.
    f_tt = tntt.zeros(Stt.N)

    # interpoalte the excitation and compute the correspinding tensor
    u0 = 1
    extitation_dofs = tt_iga.Function(Basis).interpolate(lambda t: t[:,0]*0+u0)
    tmp = np.zeros(N)
    tmp[:,-1,:] = extitation_dofs[:,-1,:].full()
    g_tt =  Pbd_tt @ (tntt.TT(tmp) ** tntt.ones([nl]*Np))

    # assemble the system matrix
    M_tt = Pin_tt@Stt@Pin_tt + Pbd_tt
    rhs_tt = Pin_tt @ (Mass_tt @ f_tt - Stt @ Pbd_tt @ g_tt) + g_tt
    M_tt = M_tt.round(1e-11)

    # solve the system
    eps_solver = 1e-7
    tme_amen = datetime.datetime.now()
    dofs_tt = tntt.solvers.amen_solve(M_tt, rhs_tt, x0 = tntt.ones(rhs_tt.N), eps = eps_solver, nswp = 50, preconditioner = 'c',  verbose = False)
    tme_amen = (datetime.datetime.now() -tme_amen).total_seconds()
    dct['time solver'] = tme_amen
    print('Time solver', tme_amen)

    if tn.cuda.is_available():
        tme_amen_gpu = datetime.datetime.now()
        dofs_tt = tntt.solvers.amen_solve(M_tt.cuda(), rhs_tt.cuda(), x0 = tntt.ones(rhs_tt.N).cuda(), eps = eps_solver, nswp = 50, preconditioner = 'c', verbose = False).cpu()
        tme_amen_gpu = (datetime.datetime.now() -tme_amen_gpu).total_seconds()
        dct['time solver GPU'] = tme_amen_gpu
        print('Time solver GPU', tme_amen_gpu)

    # save stats in the dictionary
    dct['rank matrix'] = np.mean(M_tt.R)
    dct['rank rhs'] = np.mean(rhs_tt.R)
    dct['rank solution'] = np.mean(dofs_tt.R)
    dct['memory stiff'] = tntt.numel(Stt)*8/1e6
    dct['memory mass'] = tntt.numel(Mass_tt)*8/1e6
    dct['memory system mat'] = tntt.numel(M_tt)*8/1e6
    dct['memory rhs'] = tntt.numel(rhs_tt)*8/1e6
    dct['memory solution'] = tntt.numel(dofs_tt)*8/1e6
    dct['np'] = Np
    # check the error for the case Theta = 0 (cylinder capacitor)
    fspace = tt_iga.Function(Basis+Basis_param)
    fspace.dofs = dofs_tt

    u_val = fspace([tn.linspace(0,1,8),tn.linspace(0,1,128),tn.linspace(0,1,128)]+[tn.tensor([0.0]) for i in range(Np)]).full()
    x,y,z = geom([tn.linspace(0,1,8),tn.linspace(0,1,128),tn.linspace(0,1,128)]+[tn.tensor([0.0]) for i in range(Np)])
    r = tn.sqrt(x.full()**2+y.full()**2)

    a = u0/np.log(2/(2-h))
    b = u0-a*np.log(2)

    u_ref = a*tn.log(r)+b

    err = tn.max(tn.abs(u_val-u_ref))

    dct['max_err'] = err
    print('\nMax err %e\n\n'%(err))
    return dct

if __name__ == '__main__':

    Nps = [2,3,4,5,6,7,8,9,10]

    dct_results = dict()

    for Np in Nps:
        dct = solve(Np)
        dct_results[Np] = dct

    # print header
    df = pd.DataFrame([[el for el in v.values() ] for v in dct_results.values()], columns = [k for k in dct_results[2]])
    print(df)
    df.to_pickle('./data/cshape_gpu.pickle')
    eoc = lambda x,y: np.log(y[1:]/y[:-1])/np.log(x[1:]/x[:-1])

max_stars_repo_name: ion-g-ion/code-paper-tt-iga

max_issues_repo_name: ion-g-ion/code-paper-tt-iga

max_forks_repo_name: ion-g-ion/code-paper-tt-iga


text:
import numpy as np
from collections import defaultdict

class Agent:

    def __init__(self, nA=6):
        """ Initialize agent.

        Params
        ======
        - nA: number of actions available to the agent
        """
        self.nA = nA
        self.Q = defaultdict(lambda: np.zeros(self.nA))
        self.epsilon = 0.005
        self.gamma = 1.0
        self.alpha = 0.3

    def select_action(self, state):
        """ Given the state, select an action.

        Params
        ======
        - state: the current state of the environment

        Returns
        =======
        - action: an integer, compatible with the task's action space
        """
        policy_s = self.epsilon_greedy_policy(self.Q[state], self.nA, self.epsilon)
        return np.random.choice(np.arange(self.nA), p=policy_s)

    def step(self, state, action, reward, next_state, done):
        """ Update the agent's knowledge, using the most recently sampled tuple.

        Params
        ======
        - state: the previous state of the environment
        - action: the agent's previous choice of action
        - reward: last reward received
        - next_state: the current state of the environment
        - done: whether the episode is complete (True or False)
        """
        policy_s = self.epsilon_greedy_policy(self.Q[next_state], self.nA, self.epsilon)
        self.Q[state][action] += self.alpha * (reward + self.gamma * np.dot(self.Q[next_state], policy_s) - self.Q[state][action])

    def epsilon_greedy_policy(self, Q_s, nA, epsilon=1.0):
        policy_s = np.ones(nA) * epsilon / nA
        a = np.argmax(Q_s)
        policy_s[a] = 1 - epsilon + (epsilon / nA)
        return policy_s


text:
import torch
from torch import nn
from reid import models
from reid.trainers import Trainer
from reid.evaluators import extract_features, Evaluator
from reid.dist_metric import DistanceMetric
import numpy as np
from collections import OrderedDict
import os.path as osp
import pickle
from reid.utils.serialization import load_checkpoint
from reid.utils.data import transforms as T
from torch.utils.data import DataLoader
from reid.utils.data.preprocessor import Preprocessor
import random
from sklearn.metrics.pairwise import cosine_similarity
# from run import outf
# import run
import  math
import codecs

import ssl
ssl._create_default_https_context = ssl._create_unverified_context

class EUG():
    def __init__(self, model_name, batch_size, mode, num_classes, data_dir, l_data, u_data, save_path, dropout=0.5, max_frames=900):

        self.model_name = model_name
        self.num_classes = num_classes
        self.mode = mode
        self.data_dir = data_dir
        self.save_path = save_path

        self.l_data = l_data
        self.u_data = u_data
        self.l_label = np.array([label for _,label,_,_ in l_data])
        self.u_label = np.array([label for _,label,_,_ in u_data])

        self.dataloader_params = {}
        self.dataloader_params['height'] = 256
        self.dataloader_params['width'] = 128
        self.dataloader_params['batch_size'] = batch_size
        self.dataloader_params['workers'] = 6

        self.batch_size = batch_size
        self.data_height = 256
        self.data_width = 128
        self.data_workers = 6

        # batch size for eval mode. Default is 1.
        self.eval_bs = 1
        self.dropout = dropout
        self.max_frames = max_frames

    def get_dataloader(self, dataset, training=False) :
        normalizer = T.Normalize(mean=[0.485, 0.456, 0.406],
                                 std=[0.229, 0.224, 0.225])

        if training:
            transformer = T.Compose([
                T.RandomSizedRectCrop(self.data_height, self.data_width),
                T.RandomHorizontalFlip(),
                T.ToTensor(),
                normalizer,
            ])
            batch_size = self.batch_size

        else:
            transformer = T.Compose([
                T.RectScale(self.data_height, self.data_width),
                T.ToTensor(),
                normalizer,
            ])
            batch_size = self.eval_bs

        data_loader = DataLoader(
            Preprocessor(dataset, root=self.data_dir,
                         transform=transformer, is_training=training, max_frames=self.max_frames),
            batch_size=batch_size, num_workers=self.data_workers,
            shuffle=training, pin_memory=True, drop_last=training)

        current_status = "Training" if training else "Test"
        print("create dataloader for {} with batch_size {}".format(current_status, batch_size))
        return data_loader

    def train(self, train_data, step, epochs=10, step_size=55, init_lr=0.1, dropout=0.5):

        """ create model and dataloader """
        model = models.create(self.model_name, dropout=self.dropout, num_classes=self.num_classes, mode=self.mode)
        model = nn.DataParallel(model).cuda()
        # model = nn.DataParallel(model, device_ids=[3,4]).cuda()
        # model.to(device)
        dataloader = self.get_dataloader(train_data, training=True)

        # the base parameters for the backbone (e.g. ResNet50)
        base_param_ids = set(map(id, model.module.CNN.base.parameters()))

        # we fixed the first three blocks to save GPU memory
        base_params_need_for_grad = filter(lambda p: p.requires_grad, model.module.CNN.parameters())

        # params of the new layers
        new_params = [p for p in model.parameters() if id(p) not in base_param_ids]

        # set the learning rate for backbone to be 0.1 times
        param_groups = [
            {'params': base_params_need_for_grad, 'lr_mult': 0.1},
            {'params': new_params, 'lr_mult': 1.0}]

        criterion = nn.CrossEntropyLoss().cuda()  # 标准
        optimizer = torch.optim.SGD(param_groups, lr=init_lr, momentum=0.5, weight_decay = 5e-4, nesterov=True)

        # change the learning rate by step
        def adjust_lr(epoch, step_size):     #学习率的衰减也可以做调整
            lr = init_lr / (10 ** (epoch // step_size))
            for g in optimizer.param_groups:
                g['lr'] = lr * g.get('lr_mult', 1)

            if epoch % step_size == 0:
                print("Epoch {}, current lr {}".format(epoch, lr))

        """ main training process """
        trainer = Trainer(model, criterion)
        for epoch in range(epochs):
            adjust_lr(epoch, step_size)
            trainer.train(epoch, dataloader, optimizer)
            # trainer.train(epoch, dataloader, optimizer, print_freq=len(dataloader)//30 * 10)

        torch.save(model.state_dict(), osp.join(self.save_path,  "{}_step_{}.ckpt".format(self.mode, step)))
        self.model = model

    def get_feature(self, dataset):
        dataloader = self.get_dataloader(dataset, training=False)
        features,_ = extract_features(self.model, dataloader)
        features = np.array([logit.numpy() for logit in features.values()])
        return features

    # def get_Classification_result(self):  #是指的用CNN的分类结果贴标签
    #     logits = self.get_feature(self.u_data)
    #     exp_logits = np.exp(logits)
    #     predict_prob = exp_logits / np.sum(exp_logits,axis=1).reshape((-1,1))
    #     assert len(logits) == len(predict_prob)
    #     assert predict_prob.shape[1] == self.num_classes
    #
    #     pred_label = np.argmax(predict_prob, axis=1)
    #     pred_score = predict_prob.max(axis=1)
    #     print("get_Classification_result", predict_prob.shape)
    #
    #
    #     num_correct_pred = 0
    #     for idx, p_label in enumerate(pred_label):
    #         if self.u_label[idx] == p_label:
    #             num_correct_pred +=1
    #
    #
    #     print("{} predictions on all the unlabeled data: {} of {} is correct, accuracy = {:0.3f}".format(
    #         self.mode, num_correct_pred, pred_label.shape[0], num_correct_pred/pred_label.shape[0]))
    #
    #     return pred_label, pred_score

    def get_Dissimilarity_result(self):

        # extract feature
        u_feas = self.get_feature(self.u_data)
        l_feas = self.get_feature(self.l_data)
        print("u_features", u_feas.shape, "l_features", l_feas.shape)

        scores = np.zeros((u_feas.shape[0]))
        labels = np.zeros((u_feas.shape[0]))
        # 分别用来存 _ufeas的分数和标签

        id_num = {}  #以标签名称作为字典
        # a = 1
        num_correct_pred = 0
        dists = []
        for idx, u_fea in enumerate(u_feas):
            diffs = l_feas - u_fea
            dist = np.linalg.norm(diffs,axis=1)
            index_min = np.argmin(dist)
            dists.append(dist)
            scores[idx] = - dist[index_min]  # "- dist" : more dist means less score
            labels[idx] = self.l_label[index_min] # take the nearest labled neighbor as the prediction label
            # if a:
            #     print("labels :-------------------------------------------", labels[idx])
            #     a = 0
            #     输出的结果是0.0
            # count the correct number of Nearest Neighbor prediction
            if self.u_label[idx] == labels[idx]:
                num_correct_pred +=1
            # 统计各个id的数量
            # if str(labels[idx]) in id_num.keys():
            #     id_num[str(labels[idx])]=id_num[str(labels[idx])]+1 #值加1
            # else:
            #     id_num[str(labels[idx])]=1

        print("{} predictions on all the unlabeled data: {} of {} is correct, accuracy = {:0.3f}".format(
            self.mode, num_correct_pred, u_feas.shape[0], num_correct_pred/u_feas.shape[0]))

        # sorted(id_num.items(),key = lambda item:item[1])
        # print("id_num:--------------------------------------------id_num----------------- ")
        # print(id_num)
        dists = np.vstack(dists)
        return labels, scores,num_correct_pred/u_feas.shape[0],dists

    def get_Dissimilarity_result2(self):
        # l_feas_file = codecs.open("logs/l_feas/test1.txt",'a')
        # extract feature
        u_feas = self.get_feature(self.u_data)
        l_feas = self.get_feature(self.l_data)
        l_mean = l_feas.mean(axis=0)
        l_std = l_feas.std(axis=0)
        l_feas -= l_mean
        l_feas /= l_std
        # np.save("logs/l_feas/test1.npy",l_feas)
        print("u_features", u_feas.shape, "l_features", l_feas.shape)
        scores = np.zeros((u_feas.shape[0]))
        labels = np.zeros((u_feas.shape[0]))
        # 分别用来存 _ufeas的分数和标签
        id_num = {}  #以标签名称作为字典
        num_correct_pred = 0
        for idx, u_fea in enumerate(u_feas):
            # dist = []
            # for l_fea in l_feas:
            #     d = np.dot(l_fea,u_fea)/(np.linalg.norm(l_fea)*(np.linalg.norm(u_fea)))
            #     dist.appen(d)
            # dist = cosine_similarity(l_feas,u_fea) #维度不对
            # dist = cosine_similarity(l_feas,np.array([u_fea])).reshape(-1)
            # print("----------------------------------------------------dist:{}".format(dist))
            u_fea -=l_mean
            u_fea /=l_std
            diffs = l_feas - u_fea
            dist = np.linalg.norm(diffs,axis=1)
            index_min = np.argmin(dist)
            scores[idx] = - dist[index_min]  # "- dist" : more dist means less score
            labels[idx] = self.l_label[index_min] # take the nearest labled neighbor as the prediction label
            # if a:
            #     print("labels :-------------------------------------------", labels[idx])
            #     a = 0
            #     输出的结果是0.0
            # count the correct number of Nearest Neighbor prediction
            if self.u_label[idx] == labels[idx]:
                num_correct_pred +=1
            # 统计各个id的数量
            if str(labels[idx]) in id_num.keys():
                id_num[str(labels[idx])]=id_num[str(labels[idx])]+1 #值加1
            else:
                id_num[str(labels[idx])] =1

        print("{} predictions on all the unlabeled data: {} of {} is correct, accuracy = {:0.3f}".format(
            self.mode, num_correct_pred, u_feas.shape[0], num_correct_pred/u_feas.shape[0]))

        sorted(id_num.items(),key = lambda item:item[1])
        # print("id_num:--------------------------------------------id_num----------------- ")
        # print(id_num)
        return labels, scores,num_correct_pred/u_feas.shape[0],id_num

    def estimate_label(self):

        print("label estimation by {} mode.".format(self.mode))

        if self.mode == "Dissimilarity":
            # predict label by dissimilarity cost
            [pred_label, pred_score,label_pre,id_num] = self.get_Dissimilarity_result()

        elif self.mode == "Classification":
            # predict label by classification
            [pred_label, pred_score] = self.get_Classification_result()
        else:
            raise ValueError

        return pred_label, pred_score,label_pre,id_num

    def select_top_data(self, pred_score, nums_to_select):
        v = np.zeros(len(pred_score))
        index = np.argsort(-pred_score)
        for i in range(nums_to_select):  #排序,求最前面的n个
            v[index[i]] = 1
        return v.astype('bool')

    def select_top_data_nlvm_b1(self,pred_score,dists,new_expend_nums_to_select,new_nums_to_select):
        # pred_score = pred_score.T # if necessary
        # 方案2, 求最近的P%样本的方差
        N_u,N_l = dists.shape
        stds = np.zeros(N_u)
        selection1=selection2= np.zeros(N_u, 'bool')
        index = np.argsort(-pred_score)
        selection1[index[:new_expend_nums_to_select]] = True  # selection1 就是选出 距离在一定范围内的数量
        # 求最近的P%样本的方差
        for i in index[:new_expend_nums_to_select]:
            score = - dists[i]
            # 求k近邻
            # topk = int(N_l * percent_P)
            topk = 2
            topk_idxs = np.argpartition(score, topk)[:topk]
            stds[i] = score[topk_idxs].std()   # 这里要求pre_score 要是二维的才行
        # 根据方差排序
        idxs = np.argsort(-stds)  # 这里该取负还是取正呢?
        # print(stds[idxs[:nums_to_select]])
        selection2[idxs[:new_nums_to_select]] = True
        return selection2

    def select_top_data_nlvm_b5(self,pred_score,dists,new_nums_to_select,percent_vari):
        # pred_score = pred_score.T # if necessary
        # 方案2, 求最近的P%样本的方差
        N_u,N_l = dists.shape
        vari_num = round(N_u * percent_vari)  # 取前面vari_num个
        if vari_num <= new_nums_to_select: # 就不需要在判断方差了
            return self.select_top_data(pred_score,new_nums_to_select)
        stds = np.zeros(N_u)
        selection = np.zeros(N_u,'bool')
        for i in range(N_u):  #对所有的值求方差
            score = -dists[i]
            # 求k近邻
            topk = 2
            topk_idxs = np.argpartition(score, topk)[:topk]
            stds[i] = score[topk_idxs].std()
        idxs = np.argsort(-stds)
        index = np.argsort(-pred_score[idxs[:vari_num]])  # 分数从大到小排序
        selection[index[:new_nums_to_select]] = True
        return selection

    def select_top_data_nlvm_b6(self,pred_score,dists,new_nums_to_select,percent_vari):
        # pred_score = pred_score.T # if necessary
        # 方案2, 求最近的P%样本的方差
        N_u,N_l = dists.shape
        vari_num = round(N_u * percent_vari ) # 取前面vari_num个
        if vari_num <= new_nums_to_select: # 就不需要在判断方差了
            return self.select_top_data(pred_score,new_nums_to_select)
        stds = np.zeros(N_u)
        selection1=selection2 = np.zeros(N_u,'bool')
        for i in range(N_u):  #对所有的值求方差
            score = -dists[i]
            # 求k近邻
            topk = 2
            topk_idxs = np.argpartition(score, topk)[:topk]
            stds[i] = score[topk_idxs].std()
        idxs = np.argsort(-stds)
        selection1[idxs[:vari_num]] = True
        index = np.argsort(-pred_score)  # 分数从大到小排序
        selection2[index[:new_nums_to_select]] = True
        selection = selection1 & selection2
        return selection,np.count_nonzero(selection)

    def select_top_data_NLVM(self, pred_score, nums_to_select, percent_P = 0.1, percent_N = 0.1):
        # pred_score = pred_score.T # if necessary
        N_u,N_l = pred_score.shape
        diam = pred_score.max()
        # 标记距离
        masks = np.zeros_like(pred_score, dtype='int32')
        masks[pred_score < diam * percent_P] = 1
        masks[pred_score > diam * (1-percent_N)] = -1
        stds = np.zeros(N_u)
        selection = np.zeros(N_u,'bool')
        # 计算P样本方差
        for i in range(N_u):
            score = pred_score[i] # 这个是距离
            mask = masks[i] == 1
            # print(score.std(),score[mask].std())
            if sum(mask) > 1:
                stds[i] = score[mask].std()
        # 根据方差排序
        idxs = np.argsort(-stds)
        # print(stds[idxs[:nums_to_select]])
        selection[idxs[:nums_to_select]] = True
        return selection

    def select_top_data_NLVM_2(self, pred_score, nums_to_select, percent_P = 0.1, percent_N = 0.1):
        # pred_score = pred_score.T # if necessary
        # 方案2, 求最近的P%样本的方差
        N_u,N_l = pred_score.shape
        stds = np.zeros(N_u)
        selection = np.zeros(N_u,'bool')
        # 求最近的P%样本的方差
        for i in range(N_u):
            score = pred_score[i]
            # 求k近邻
            topk = int(N_l * percent_P)
            topk_idxs = np.argpartition(score,topk)[:topk]
            stds[i] = score[topk_idxs].std()
        # 根据方差排序
        idxs = np.argsort(-stds)
        # print(stds[idxs[:nums_to_select]])
        selection[idxs[:nums_to_select]] = True
        return selection

    def select_top_data3(self, pred_score, nums_to_select,id_num,pred_y,u_data):
        total_number = 0
        for item in id_num:
            id_num[item] = round(id_num[item] * nums_to_select / len(u_data))  #向下取整/ 四舍五入
            total_number = total_number+id_num[item]

        print("nums_to_select vs total_number = {} vs {}".format(nums_to_select,total_number))
        v = np.zeros(len(pred_score))
        index = np.argsort(-pred_score)
        count = 0
        for i in range(len(pred_score)):
            if count == total_number:
                break
            if round(id_num[str(pred_y[i])]):
                v[index[i]] = 1
                count  = count+1
                id_num[str(pred_y[i])] = id_num[str(pred_y[i])]-1
        return v.astype('bool')

    def generate_new_train_data(self, sel_idx, pred_y):
        """ generate the next training data """

        seletcted_data = []
        correct, total = 0, 0
        for i, flag in enumerate(sel_idx):
            if flag: # if selected
                seletcted_data.append([self.u_data[i][0], int(pred_y[i]), self.u_data[i][2], self.u_data[i][3]])
                total += 1
                if self.u_label[i] == int(pred_y[i]):
                    correct += 1
        acc = correct / total

        new_train_data = self.l_data + seletcted_data
        print("selected pseudo-labeled data: {} of {} is correct, accuracy: {:0.4f}  new train data: {}".format(
                correct, len(seletcted_data), acc, len(new_train_data)))

        return new_train_data,acc

    def resume(self, ckpt_file, step):
        print("continued from step", step)
        model = models.create(self.model_name, dropout=self.dropout, num_classes=self.num_classes, mode=self.mode)
        self.model = nn.DataParallel(model).cuda()
        self.model.load_state_dict(load_checkpoint(ckpt_file))

    def evaluate(self, query, gallery):
        test_loader = self.get_dataloader(list(set(query) | set(gallery)), training = False)
        evaluator = Evaluator(self.model)
        return  evaluator.evaluate(test_loader, query, gallery)

"""
    Get one-shot split for the input dataset.
"""
def get_one_shot_in_cam1(dataset, load_path, seed=0):

    np.random.seed(seed)
    random.seed(seed)

    # if previous split exists, load it and return
    if osp.exists(load_path):
        with open(load_path, "rb") as fp:
            dataset = pickle.load(fp)
            label_dataset = dataset["label set"]
            unlabel_dataset = dataset["unlabel set"]

        print("  labeled  |   N/A | {:8d}".format(len(label_dataset)))
        print("  unlabel  |   N/A | {:8d}".format(len(unlabel_dataset)))
        print("\nLoad one-shot split from", load_path)
        return label_dataset, unlabel_dataset

    #print("random create new one-shot split and save it to", load_path)

    label_dataset = []
    unlabel_dataset = []

    # dataset indexed by [pid, cam]
    dataset_in_pid_cam = [[[] for _ in range(dataset.num_cams)] for _ in range(dataset.num_train_ids) ]
    for index, (images, pid, camid, videoid) in enumerate(dataset.train):
        dataset_in_pid_cam[pid][camid].append([images, pid, camid, videoid])

    # generate the labeled dataset by randomly selecting a tracklet from the first camera for each identity
    for pid, cams_data  in enumerate(dataset_in_pid_cam):
        for camid, videos in enumerate(cams_data):
            if len(videos) != 0:
                selected_video = random.choice(videos)
                break
        label_dataset.append(selected_video)
    assert len(label_dataset) == dataset.num_train_ids
    labeled_videoIDs =[vid for _, (_,_,_, vid) in enumerate(label_dataset)]

    # generate unlabeled set
    for (imgs, pid, camid, videoid) in dataset.train:
        if videoid not in labeled_videoIDs:
            unlabel_dataset.append([imgs, pid, camid, videoid])

    with open(load_path, "wb") as fp:
        pickle.dump({"label set":label_dataset, "unlabel set":unlabel_dataset}, fp)

    print("  labeled    | N/A | {:8d}".format(len(label_dataset)))
    print("  unlabeled  | N/A | {:8d}".format(len(unlabel_dataset)))
    print("\nCreate new one-shot split, and save it to", load_path)
    return label_dataset, unlabel_dataset

def get_one_shot_in_cam2(dataset, load_path, seed=0):

    np.random.seed(seed)
    random.seed(seed)

    # if previous split exists, load it and return
    if osp.exists(load_path):
        with open(load_path, "rb") as fp:
            dataset = pickle.load(fp)
            label_dataset = dataset["label set"]
            unlabel_dataset = dataset["unlabel set"]

        print("  labeled  |   N/A | {:8d}".format(len(label_dataset)))
        print("  unlabel  |   N/A | {:8d}".format(len(unlabel_dataset)))
        print("\nLoad one-shot split from", load_path)
        return label_dataset+unlabel_dataset,[]

    #print("random create new one-shot split and save it to", load_path)

    label_dataset = []
    unlabel_dataset = []

    # dataset indexed by [pid, cam]
    dataset_in_pid_cam = [[[] for _ in range(dataset.num_cams)] for _ in range(dataset.num_train_ids) ]
    for index, (images, pid, camid, videoid) in enumerate(dataset.train):
        dataset_in_pid_cam[pid][camid].append([images, pid, camid, videoid])

    # generate the labeled dataset by randomly selecting a tracklet from the first camera for each identity
    for pid, cams_data  in enumerate(dataset_in_pid_cam):
        for camid, videos in enumerate(cams_data):
            if len(videos) != 0:
                selected_video = random.choice(videos)
                break
        label_dataset.append(selected_video)
    assert len(label_dataset) == dataset.num_train_ids
    labeled_videoIDs =[vid for _, (_,_,_, vid) in enumerate(label_dataset)]

    # generate unlabeled set
    for (imgs, pid, camid, videoid) in dataset.train:
        if videoid not in labeled_videoIDs:
            unlabel_dataset.append([imgs, pid, camid, videoid])

    with open(load_path, "wb") as fp:
        pickle.dump({"label set":label_dataset, "unlabel set":unlabel_dataset}, fp)

    print("  labeled    | N/A | {:8d}".format(len(label_dataset)))
    print("  unlabeled  | N/A | {:8d}".format(len(unlabel_dataset)))
    print("\nCreate new one-shot split, and save it to", load_path)
    return label_dataset+unlabel_dataset,[]