text:
import re
import collections
from string import Template
from itertools import product

import numpy as np
from joblib import Parallel, delayed
from pyparsing import (
    Word,
    alphanums,
    Suppress,
    Optional,
    CharsNotIn,
    Group,
    nums,
    ZeroOrMore,
    OneOrMore,
    cppStyleComment,
    printables,
)

from pgmpy.models import BayesianModel
from pgmpy.factors.discrete import TabularCPD

class BLOGWriter(object):

    """
    Base class for writing BLOG network file format
    """

    def __init__(self, model):
        """
        Initialise a BLOGWriter Object

        Parameters
        ----------
        model: BayesianModel Instance

        Examples
        ---------
        >>> from pgmpy.readwrite import BLOGWriter
        >>> from pgmpy.utils import get_example_model
        >>> asia = get_example_model('asia')
        >>> writer = BLOGWriter(asia)
        >>> writer
        <writer_BLOG.BLOGWriter at 0x7f05e5ea27b8>
        >>> writer.write_BLOG('asia.BLOG')
        """
        if not isinstance(model, BayesianModel):
            raise TypeError("model must be an instance of BayesianModel")
        self.model = model
        if not self.model.name:
            self.network_name = "unknown"
        else:
            self.network_name = self.model.name
        self.variable_states = self.get_states()
        self.property_tag = self.get_properties()
        self.variable_parents = self.get_parents()
        self.tables = self.get_cpds()

    def BLOG_templates(self):
        """
        Create template for writing in BLOG format
        """
        network_template = Template("network $name {\n}\n")
        # property tag may or may not be present in model,and since no of properties
        # can be more than one , will replace them accoriding to format otherwise null
        variable_template = Template(
            """type Type$name;\ndistinct Type$name $states;\n\n"""
        )
        property_template = Template("    property $prop ;\n")
        # $variable_ here is name of variable, used underscore for clarity
        probability_template = Template(
            """random Type$variable_ $variable_ ~ $parents\n\t$values;\n\n"""
        )
        return (
            variable_template,
            probability_template,
        )

    def __str__(self):
        """
        Returns the BLOG format as string
        """
        (
            variable_template,
            probability_template,
        ) = self.BLOG_templates()
        variables = self.model.nodes()

        network = ""

        for var in sorted(variables):
            no_of_states = str(len(self.variable_states[var]))
            variable_states = [var + s for s in self.variable_states[var]]
            states = ", ".join(variable_states)
            network += variable_template.substitute(
                name=var,
                no_of_states=no_of_states,
                states=states,
            )

        for var in sorted(variables):
            if not self.variable_parents[var]:
                parents = ""
                cpd = "Categorical({"+self.tables[var]+"})\n"
            else:
                if len(self.variable_parents[var]) == 1:
                    parents = "case "+", ".join(self.variable_parents[var])+" in"
                else:
                    parents = "case ["+", ".join(self.variable_parents[var])+"] in"
                cpd = "{"+self.tables[var]+"}"
            network += probability_template.substitute(
                variable_=var, parents=parents, values=cpd
            )

        network += "query "+var+";\n"

        return network

    def get_variables(self):
        """
        Add variables to BLOG

        Returns
        -------
        list: a list containing names of variable

        Example
        -------
        >>> from pgmpy.readwrite import BLOGReader, BLOGWriter
        >>> model = BLOGReader('dog-problem.BLOG').get_model()
        >>> writer = BLOGWriter(model)
        >>> writer.get_variables()
        ['bowel-problem', 'family-out', 'hear-bark', 'light-on', 'dog-out']
        """
        variables = self.model.nodes()
        return variables

    def get_states(self):
        """
        Add states to variable of BLOG

        Returns
        -------
        dict: dict of type {variable: a list of states}

        Example
        -------
        >>> from pgmpy.readwrite import BLOGReader, BLOGWriter
        >>> model = BLOGReader('dog-problem.BLOG').get_model()
        >>> writer = BLOGWriter(model)
        >>> writer.get_states()
        {'bowel-problem': ['bowel-problem_0', 'bowel-problem_1'],
         'dog-out': ['dog-out_0', 'dog-out_1'],
         'family-out': ['family-out_0', 'family-out_1'],
         'hear-bark': ['hear-bark_0', 'hear-bark_1'],
         'light-on': ['light-on_0', 'light-on_1']}
        """
        variable_states = {}
        cpds = self.model.get_cpds()
        for cpd in cpds:
            variable = cpd.variable
            variable_states[variable] = []
            for state in cpd.state_names[variable]:
                variable_states[variable].append(str(state))
        return variable_states

    def get_properties(self):
        """
        Add property to variables in BLOG

        Returns
        -------
        dict: dict of type {variable: list of properties }

        Example
        -------
        >>> from pgmpy.readwrite import BLOGReader, BLOGWriter
        >>> model = BLOGReader('dog-problem.BLOG').get_model()
        >>> writer = BLOGWriter(model)
        >>> writer.get_properties()
        {'bowel-problem': ['position = (335, 99)'],
         'dog-out': ['position = (300, 195)'],
         'family-out': ['position = (257, 99)'],
         'hear-bark': ['position = (296, 268)'],
         'light-on': ['position = (218, 195)']}
        """
        variables = self.model.nodes()
        property_tag = {}
        for variable in sorted(variables):
            properties = self.model.nodes[variable]
            properties = collections.OrderedDict(sorted(properties.items()))
            property_tag[variable] = []
            for prop, val in properties.items():
                property_tag[variable].append(str(prop) + " = " + str(val))
        return property_tag

    def get_parents(self):
        """
        Add the parents to BLOG

        Returns
        -------
        dict: dict of type {variable: a list of parents}

        Example
        -------
        >>> from pgmpy.readwrite import BLOGReader, BLOGWriter
        >>> model = BLOGReader('dog-problem.BLOG').get_model()
        >>> writer = BLOGWriter(model)
        >>> writer.get_parents()
        {'bowel-problem': [],
         'dog-out': ['bowel-problem', 'family-out'],
         'family-out': [],
         'hear-bark': ['dog-out'],
         'light-on': ['family-out']}
        """
        cpds = self.model.get_cpds()
        variable_parents = {}
        for cpd in cpds:
            variable_parents[cpd.variable] = cpd.variables[1:]
        return variable_parents

    def get_cpds(self):
        """
        Adds tables to BLOG

        Returns
        -------
        dict: dict of type {variable: array}

        Example
        -------
        >>> from pgmpy.readwrite import BLOGReader, BLOGWriter
        >>> model = BLOGReader('dog-problem.BLOG').get_model()
        >>> writer = BLOGWriter(model)
        >>> writer.get_cpds()
        {'bowel-problem': array([ 0.01,  0.99]),
         'dog-out': array([ 0.99,  0.97,  0.9 ,  0.3 ,  0.01,  0.03,  0.1 ,  0.7 ]),
         'family-out': array([ 0.15,  0.85]),
         'hear-bark': array([ 0.7 ,  0.01,  0.3 ,  0.99]),
         'light-on': array([ 0.6 ,  0.05,  0.4 ,  0.95])}
        """
        cpds = self.model.get_cpds()
        tables = {}
        for cpd in cpds:
            print(cpd)
            evidence = cpd.variables[1:]
            evidence_card = cpd.cardinality[1:]
            o_value = []
            print(evidence)
            if evidence:
                col_indexes = np.array(list(product(*[range(i) for i in evidence_card])))
                for m in range(np.shape(col_indexes)[0]):
                    col_name = []
                    var_name = []
                    cat_value = []
                    for k in range(np.size(col_indexes[m])):
                        var_name.append("{var}".format(var=evidence[k]))
                        col_name.append("{var}{state}".format(state=cpd.state_names[evidence[k]][col_indexes[m][k]], var=evidence[k]))
                    for j in range(cpd.variable_card):
                        row_name = "{var}{state}".format(
                            var=cpd.variable, state=cpd.state_names[cpd.variable][j]
                            )
                        index = np.shape(col_indexes)[0] * j + m
                        row_index = row_name +" -> " + str(cpd.values.ravel()[index])
                        cat_value.append(row_index)

                    var_name = ', '.join(var_name)
                    cat_value= ', '.join(cat_value)
                    if len(col_name) == 1:
                        o_value_temp = ''.join([ col_name[0], " -> " , "Categorical({", cat_value, "})"])
                    else:
                        col_name = ', '.join(col_name)
                        o_value_temp = ''.join(["[ ", col_name, " ]"," -> " , "Categorical({", cat_value, "})"])
                    o_value.append(o_value_temp)
                tables[cpd.variable] = ',\n'.join(o_value)
            else:
                for j in range(cpd.variable_card):
                    row_name = "{var}({state})".format(
                            var=cpd.variable, state=cpd.state_names[cpd.variable][j]
                            )
                    index = j
                    print(row_name)
                    col_name="{var}{state}".format(state=cpd.state_names[cpd.variable][j], var=cpd.variable)
                    o_value_temp = "".join([col_name, " -> ", str(cpd.values.ravel()[index])])
                    o_value.append(o_value_temp)
                tables[cpd.variable] = ','.join(o_value)

        return tables

    def write_BLOG(self, filename):
        """
        Writes the BLOG data into a file

        Parameters
        ----------
        filename : Name of the file

        Example
        -------
        >>> from pgmpy.utils import get_example_model
        >>> from pgmpy.readwrite import BLOGReader, BLOGWriter
        >>> asia = get_example_model('asia')
        >>> writer = BLOGWriter(asia)
        >>> writer.write_BLOG(filename='asia.BLOG')
        """
        writer = self.__str__()
        with open(filename, "w") as fout:
            fout.write(writer)

from pgmpy.readwrite import BIFReader, BIFWriter
import os
from bif2blog import BLOGWriter
if __name__ ==  "__main__":
    file = "./survey.bif"
    if file.endswith(".bif"):
        print("convert BIF file:", file)
        model = BIFReader(file).get_model()
        writer = BLOGWriter(model)
        new_name = os.path.splitext(file)[0]+".blog"
        writer.write_BLOG(filename=new_name)

ext: py

max_stars_repo_name: shiruizhao/swift

max_issues_repo_name: shiruizhao/swift

max_forks_repo_name: shiruizhao/swift

reason: import numpy


text:
import csv
import numpy as np
import pandas as pd
import itertools as it
import venn
import matplotlib.pyplot as plt

default_infile = 'data/tracks_tags.csv'
default_outfile = 'results/stats.csv'
accepted_types = {
    # 'genre',
    # 'audiodescriptor',
    # 'uncategorized',
    # 'instrument',
    'vartags'
}
accepted_sources = {
    'team',
    'artist',
    # 'bmat',
    # 'mturk'
    # 'auto1',
    # 'auto2',
}

def read_data(filename=default_infile):
    tags_map = {}
    tag_types = set()
    tag_sources = set()
    track_ids = set()

    with open(filename) as fp:
        reader = csv.reader(fp)
        next(reader, None)  # skip header
        for track_id, tag_type, value, source in reader:
            tag_types.add(tag_type)
            tag_sources.add(source)
            if tag_type in accepted_types and source in accepted_sources:
                track_ids.add(track_id)
                if value not in tags_map:
                    tags_map[value] = set()
                tags_map[value].add(track_id)

    return tags_map, len(track_ids)

def get_top_tags(tags_map):
    tags = np.array(list(tags_map.keys()))
    counts = np.array([len(vals) for vals in tags_map.values()])
    map = counts.argsort()[::-1]
    return tags[map], counts[map]

def export_stats(tags, counts, filename=default_outfile):
    with open(filename, 'w') as fp:
        writer = csv.writer(fp)
        writer.writerow(['Tag', 'Count'])
        for tag, count in zip(tags, counts):
            writer.writerow([tag, count])

def print_tags(tags_map):
    for tag_value, track_ids in tags_map.items():
        print('{}: {}'.format(tag_value, len(track_ids)))

if __name__ == '__main__':
    data, total = read_data()
    df = pd.read_csv('results/stats_intersected.csv')
    quads = {1: 'I', 2: 'II', 3: 'III', 4: 'IV'}
    tracks = {quad_name: set() for quad_name in quads.values()}
    for quad_i, quad_name in quads.items():
        for tag in df[df.Quad == quad_i].Tags:
            tracks[quad_name] |= data[tag]
        print(quad_name, len(tracks[quad_name]))

    venn.venn(tracks)
    plt.savefig('results/quads.png', bbox_inches='tight')
    plt.title('Tracks ')
    plt.show()

    # print('Total: {}'.format(total))
    # tags, counts = get_top_tags(data)
    # export_stats(tags, counts)

max_stars_repo_name: philtgun/medieval-esr

max_issues_repo_name: philtgun/medieval-esr

max_forks_repo_name: philtgun/medieval-esr


text:
import random
from os import listdir
from os.path import join, isfile

import numpy as np
from skimage.io import imread

class SemanticWorms:
    def __init__(self, path, balance=False):
        self.path = path

        files = [f for f in listdir(self.path) if isfile(join(self.path, f))]
        self.all_names = [n[:-len('.tif')] for n in files if n.endswith('.tif')]

        self.consts = np.load('data/consts.npy').tolist()

        # split data
        self.all_names = np.random.permutation(self.all_names)
        number = len(self.all_names)
        self.train_names = self.all_names[:-number / 10 * 2]
        self.test_names = self.all_names[-number / 10 * 2:-number / 10]
        self.valid_names = self.all_names[-number / 10:]

        self.cache = {}

        self.train_probs = [1. / float(len(self.train_names))] * len(self.train_names)
        self.balance = balance
        if balance:
            self.train_probs = self.build_stats()

    def build_stats(self):

        leaf_stats = []
        for j in range(len(self.train_names)):
            x, y, m = self.get_train(j)
            leafs = []
            for i in range(1, y.shape[0]):
                leaf = np.sum(y[i])
                if leaf > 0:
                    leafs.append(leaf)
            leaf_stats.append(np.mean(leafs))

        hist = np.histogram(np.array(leaf_stats), bins=5)
        leaf_prob = np.array([1.] * len(leaf_stats))
        for i in range(len(hist[0])):
            pos = np.where(np.logical_and(leaf_stats >= hist[1][i], leaf_stats <= hist[1][i + 1]))
            leaf_prob[pos] = 1. - hist[0][i] / float(len(leaf_stats))
        return leaf_prob / np.sum(leaf_prob)

    def make_all_train(self):
        self.train_names = self.all_names
        self.train_probs = [1. / float(len(self.train_names))] * len(self.train_names)
        if self.balance:
            self.train_probs = self.build_stats()

    def __get_image__(self, basename):
        rgb = imread(join(self.path, basename + '.tif')).astype(np.float32)
        label = imread(join(self.path, basename + '_labels.png'), 'pil')
        # label[label>0]=1
        rgb = (rgb - self.consts['min']) / (self.consts['max'] - self.consts['min']) - 0.5

        return rgb, label

    def __get_labeled__(self, names, n=None):
        if n is None:
            if 1 == len(names):
                n = 0
            else:
                n = np.random.randint(0, len(names) - 1)

        name = names[n]
        if name in self.cache:
            return self.cache[name]

        # print(name)
        rgb, binary = self.__get_image__(name)

        self.cache[name] = (rgb, binary)
        return self.cache[name]

    def get_test(self, n=None):
        x, y = self.__get_image__(self.test_names[n])
        return x, y, None

    def get_train(self, n=None):
        if n is None:
            n = np.random.choice(range(len(self.train_names)), 1, p=self.train_probs)[0]
        return self.__get_labeled__(self.train_names, n)

    def get_valid(self, n=None):
        return self.__get_labeled__(self.valid_names, n)

    def __get_batch__(self, get_func, batch_size=3, transforms=[]):
        images = []
        labels = []

        for t in transforms:
            if hasattr(t, 'prepare'):
                t.prepare()

        for index in range(batch_size):
            i, l = get_func(None)

            i = i
            l = l

            for t in transforms:
                seed = random.randint(0, 2 ** 32 - 1)
                np.random.seed(seed)
                i = t(i, True)
                np.random.seed(seed)
                l = t(l)

            images.append(np.expand_dims(np.expand_dims(i, 0),0))
            labels.append(np.expand_dims(l, 0))

        labels = np.clip(np.vstack(labels).astype(np.int), 0, 3)

        return np.vstack(images), labels

    def get_train_batch(self, batch_size=3, transforms=[]):
        return self.__get_batch__(self.get_train, batch_size=batch_size, transforms=transforms)

    def get_valid_batch(self, batch_size=3, transforms=[]):
        return self.__get_batch__(self.get_valid, batch_size=batch_size, transforms=transforms)

max_stars_repo_name: kulikovv/DeepWorm

max_issues_repo_name: kulikovv/DeepWorm

max_forks_repo_name: kulikovv/DeepWorm


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-

import os
import h5py
import numpy as np
import pandas as pd
from tqdm import tqdm
from multiprocessing import Pool

from acannon.data import load_light_curves
from acannon.spectrum import get_power_spectrum

OUTPUT_DIR = "spectra"
os.makedirs(OUTPUT_DIR, exist_ok=True)

def run(row):
    kicid = int(row.kicid)
    t, f, fe = load_light_curves(kicid)
    freq, power = get_power_spectrum(t, f, fe)
    data = np.array(list(zip(freq, power)),
                    dtype=[("freq_uHz", float), ("power", float)])
    with h5py.File(os.path.join(OUTPUT_DIR, "{0}.h5".format(kicid)), "w") as f:
        f.attrs["kicid"] = kicid
        f.attrs["nu_max"] = float(row.nu_max)
        f.attrs["delta_nu"] = float(row.delta_nu)
        f.create_dataset("power", data=data)

df = pd.read_csv("stello_2013.dat", delim_whitespace=True,
                 names=["kicid", "nu_max", "delta_nu"])
pool = Pool()
list(tqdm(pool.imap_unordered(run, (
    row for _, row in df.iterrows()
)), total=len(df)))

max_stars_repo_name: dfm/AsteroCannon

max_issues_repo_name: dfm/AsteroCannon

max_forks_repo_name: dfm/AsteroCannon


text:
from __future__ import division
import os
from pyomo.environ import *

def define_arguments(argparser):
    argparser.add_argument('--psip-force', action='store_true', default=True,
        help="Force following of PSIP plans (retiring AES and building certain technologies).")
    argparser.add_argument('--psip-relax', dest='psip_force', action='store_false',
        help="Relax PSIP plans, to find a more optimal strategy.")
    argparser.add_argument('--psip-minimal-renewables', action='store_true', default=False,
        help="Use only the amount of renewables shown in PSIP plans, and no more (should be combined with --psip-relax).")
    argparser.add_argument('--force-build', nargs=3, default=None,
        help="Force construction of at least a certain quantity of a particular technology during certain years. Space-separated list of year, technology and quantity.")

def define_components(m):
    ###################
    # resource rules to match HECO's 2016-04-01 PSIP
    ##################

    # decide whether to enforce the PSIP preferred plan
    # if an environment variable is set, that takes precedence
    # (e.g., on a cluster to override options.txt)
    psip_env_var = os.environ.get('USE_PSIP_PLAN')
    if psip_env_var is None:
        # no environment variable; use the --psip-relax flag
        psip = m.options.psip_force
    elif psip_env_var.lower() in ["1", "true", "y", "yes", "on"]:
        psip = True
    elif psip_env_var.lower() in ["0", "false", "n", "no", "off"]:
        psip = False
    else:
        raise ValueError('Unrecognized value for environment variable USE_PSIP_PLAN={} (should be 0 or 1)'.format(psip_env_var))

    if psip:
        print "Using PSIP construction plan."
    else:
        print "Relaxing PSIP construction plan."

    # force conversion to LNG in 2021
    # force use of containerized LNG
    # don't allow addition of anything other than those specified here
    # force retirement of AES in 2021

    # targets for individual generation technologies
    # (year, technology, MW added)
    # TODO: allow either CentralFixedPV or CentralTrackingPV for utility-scale solar
    # (not urgent now, since CentralFixedPV is not currently in the model)

    # technologies that are definitely being built (we assume near-term
    # are underway and military projects are being built for their own
    # reasons)
    technology_targets_definite = [
        (2016, 'CentralTrackingPV', 27.6),  # Waianae Solar by Eurus Energy America
        (2018, 'IC_Schofield', 54.0),
        (2018, 'IC_Barge', 100.0),         # JBPHH plant
        (2021, 'IC_MCBH', 27.0),
        # Distributed PV from Figure J-19
        (2016, 'DistPV',  443.993468266547 - 210), # net of 210 MW of pre-existing DistPV
        (2017, 'DistPV',  92.751756737742),
        (2018, 'DistPV',  27.278236032368),
        (2019, 'DistPV',  26.188129564885),
    ]
    # technologies proposed in PSIP but which may not be built if a
    # better plan is found
    technology_targets_psip = [
        (2018, 'OnshoreWind', 24),      # NPM wind
        (2018, 'CentralTrackingPV', 109.6),  # replacement for canceled SunEdison projects
        (2018, 'OnshoreWind', 10),      # CBRE wind
        (2018, 'CentralTrackingPV', 15),  # CBRE PV
        (2020, 'OnshoreWind', 30),
        (2020, 'CentralTrackingPV', 60),
        (2021, 'CC_383', 383.0),
        (2030, 'CentralTrackingPV', 100),
        (2030, 'OffshoreWind', 200),
        (2040, 'CentralTrackingPV', 200),
        (2040, 'OffshoreWind', 200),
        (2045, 'CentralTrackingPV', 300),
        (2045, 'OffshoreWind', 400),
        (2020, 'DistPV',  21.8245069017911),
        (2021, 'DistPV',  15.27427771741),
        (2022, 'DistPV',  12.0039583149589),
        (2023, 'DistPV',  10.910655054315),
        (2024, 'DistPV',  10.913851847475),
        (2025, 'DistPV',  10.910655054316),
        (2026, 'DistPV',  9.82054858683205),
        (2027, 'DistPV',  10.910655054316),
        (2028, 'DistPV',  10.910655054315),
        (2029, 'DistPV',  14.1873680430859),
        (2030, 'DistPV',  9.82054858683205),
        (2031, 'DistPV',  10.913851847475),
        (2032, 'DistPV',  9.82054858683193),
        (2033, 'DistPV',  14.1841712499261),
        (2034, 'DistPV',  7.64033565186492),
        (2035, 'DistPV',  13.094064782442),
        (2036, 'DistPV',  9.82054858683205),
        (2037, 'DistPV',  10.9202454337949),
        (2038, 'DistPV',  9.66989970917803),
        (2039, 'DistPV',  12.1514103994531),
        (2040, 'DistPV',  12.2397218104919),
        (2041, 'DistPV',  11.7673956211361),
        (2042, 'DistPV',  10.9106550543149),
        (2043, 'DistPV',  9.82054858683205),
        (2044, 'DistPV',  15.27747451057),
        (2045, 'DistPV',  10.291675978754),
    ]

    if m.options.force_build is not None:
        b = list(m.options.force_build)
        b[0] = int(b[0])    # year
        b[2] = float(b[2])  # quantity
        b = tuple(b)
        print "Forcing build: {}".format(b)
        technology_targets_definite.append(b)

    if psip:
        technology_targets = technology_targets_definite + technology_targets_psip
    else:
        technology_targets = technology_targets_definite

    def technology_target_init(m, per, tech):
        """Find the amount of each technology that is targeted to be built by the end of each period."""
        start = 2000 if per == m.PERIODS.first() else per
        end = per + m.period_length_years[per]
        target = sum(
            mw for (tyear, ttech, mw) in technology_targets
                if ttech == tech and start <= tyear and tyear < end
        )
        return target
    m.technology_target = Param(m.PERIODS, m.GENERATION_TECHNOLOGIES, initialize=technology_target_init)

    # with PSIP: BuildGen is zero except for technology_targets
    #     (sum during each period or before first period)
    # without PSIP: BuildGen is >= definite targets
    def Enforce_Technology_Target_rule(m, per, tech):
        """Enforce targets for each technology; exact target for PSIP cases, minimum target for non-PSIP."""

        def adjust_psip_credit(g, target):
            if g in m.DISCRETELY_SIZED_GENS and target > 0.0:
                # Rescale so that the n integral units that come closest
                # to the target gets counted as the n.n fractional units
                # needed to exactly meet the target.
                # This is needed because some of the targets are based on
                # nominal unit sizes rather than actual max output.
                return (target / m.gen_unit_size[g]) / round(target / m.gen_unit_size[g])
            else:
                return 1.0

        target = m.technology_target[per, tech]
        build = sum(
            m.BuildGen[g, per] * adjust_psip_credit(g, target)
            for g in m.GENERATION_PROJECTS
            if m.gen_tech[g] == tech and (g, per) in m.GEN_BLD_YRS
        )

        if type(build) is int and build == 0:
            # no matching projects found
            if target == 0:
                return Constraint.Skip
            else:
                print(
                    "WARNING: target was set for {} in {}, but no matching projects are available. "
                    "Model will be infeasible.".format(tech, per)
                )
                return Constraint.Infeasible
        elif psip:
            return (build == target)
        elif m.options.psip_minimal_renewables and any(txt in tech for txt in ["PV", "Wind", "Solar"]):
            # only build the specified amount of renewables, no more
            return (build == target)
        else:
            # treat the target as a lower bound
            return (build >= target)
    m.Enforce_Technology_Target = Constraint(
        m.PERIODS, m.GENERATION_TECHNOLOGIES, rule=Enforce_Technology_Target_rule
    )

    aes_g = 'Oahu_AES'
    aes_size = 180
    aes_bld_year = 1992
    m.AES_OPERABLE_PERIODS = Set(initialize = lambda m:
        m.PERIODS_FOR_GEN_BLD_YR[aes_g, aes_bld_year]
    )
    m.OperateAES = Var(m.AES_OPERABLE_PERIODS, within=Binary)
    m.Enforce_AES_Deactivate = Constraint(m.TIMEPOINTS, rule=lambda m, tp:
        Constraint.Skip if (aes_g, tp) not in m.GEN_TPS
        else (m.DispatchGen[aes_g, tp] <= m.OperateAES[m.tp_period[tp]] * aes_size)
    )
    m.AESDeactivateFixedCost = Expression(m.PERIODS, rule=lambda m, per:
        0.0 if per not in m.AES_OPERABLE_PERIODS
        else - m.OperateAES[per] * aes_size * m.gen_fixed_om[aes_g, aes_bld_year]
    )
    m.Cost_Components_Per_Period.append('AESDeactivateFixedCost')

    if psip:
        # keep AES active until 2022 or just before; deactivate after that
        m.PSIP_Retire_AES = Constraint(m.AES_OPERABLE_PERIODS, rule=lambda m, per:
            (m.OperateAES[per] == 1) if per + m.period_length_years[per] <= 2022
            else (m.OperateAES[per] == 0)
        )

        # before 2040: no biodiesel, and only 100-300 GWh of non-LNG fossil fuels
        # period including 2040-2045: <= 300 GWh of oil; unlimited biodiesel or LNG

        # no biodiesel before 2040 (then phased in fast enough to meet the RPS)
        m.EARLY_BIODIESEL_MARKETS = Set(dimen=2, initialize=lambda m: [
            (rfm, per)
                for per in m.PERIODS if per + m.period_length_years[per] <= 2040
                    for rfm in m.REGIONAL_FUEL_MARKETS if m.rfm_fuel == 'Biodiesel'
        ])
        m.NoEarlyBiodiesel = Constraint(m.EARLY_BIODIESEL_MARKETS, rule=lambda m, rfm, per:
            m.FuelConsumptionInMarket[rfm, per] == 0
        )

        # # 100-300 GWh of non-LNG fuels in 2021-2040 (based on 2016-04 PSIP fig. 5-5)
        # # Note: this is needed because we assume HECO plans to burn LNG in the future
        # # even in scenarios where it costs more than oil.
        # m.PSIP_HIGH_LNG_PERIODS = Set(initialize=lambda m:
        #     [per for per in m.PERIODS if per + m.period_length_years[per] > 2021 and per < 2045]
        # )
        # m.OilProductionGWhPerYear = Expression(m.PERIODS, rule=lambda m, per:
        #     sum(
        #         m.DispatchGenByFuel[g, tp, f] * m.tp_weight_in_year[tp] * 0.001 # convert from MWh to GWh
        #             for f in ['Diesel', 'LSFO', 'LSFO-Diesel-Blend']
        #                 for g in m.GENS_BY_FUEL[f]
        #                     for tp in m.TPS_IN_PERIOD[per] if (g, tp) in m.GEN_TPS
        #     )
        # )
        # m.Upper_Limit_Oil_Power = Constraint(m.PERIODS, rule=lambda m, per:
        #     (m.OilProductionGWhPerYear[per] <= 300)
        #         if per + 0.5 * m.period_length_years[per] >= 2021
        #     else
        #         Constraint.Skip
        # )
        # # lower limit is in place to roughly reflect HECO's plan
        # m.Lower_Limit_Oil_Power = Constraint(m.PERIODS, rule=lambda m, per:
        #     (m.OilProductionGWhPerYear[per] >= 100)
        #         if per + m.period_length_years[per] < 2040  # relax constraint if period ends after 2040
        #     else
        #         Constraint.Skip
        # )

        # force LNG conversion in 2021 (modeled on similar constraint in lng_conversion.py)
        # This could have extra code to skip the constraint if there are no periods after 2021,
        # but it is unlikely ever to be run that way.
        # Note: this is not needed if some plants are forced to run on LNG
        # NOTE: this is no longer used; use '--force-lng-tier container' instead
        # m.PSIP_Force_LNG_Conversion = Constraint(m.LOAD_ZONES, rule=lambda m, z:
        #         m.ConvertToLNG[
        #             z,
        #             min(per for per in m.PERIODS if per + m.period_length_years[per] > 2021)
        #         ] == 1
        #     )

        # # Kahe 5, Kahe 6, Kalaeloa and CC_383 only burn LNG after 2021
        # # This is not used because it creates a weird situation where HECO runs less-efficient non-LNG
        # # plants instead of more efficient LNG-capable plants on oil.
        # # there may be a faster way to build this, but it's not clear what
        # m.PSIP_Force_LNG_Use = Constraint(m.GEN_TP_FUELS, rule=lambda m, g, tp, fuel:
        #     (m.GenFuelUseRate[g, tp, fuel] == 0)
        #         if g in m.LNG_CONVERTED_PLANTS
        #             and fuel != 'LNG'
        #             and m.tp_period[tp] + m.period_length_years[m.tp_period[tp]] > 2021
        #     else
        #         Constraint.Skip
        # )

        # don't allow construction of any advanced technologies (e.g., batteries, pumped hydro, fuel cells)
        advanced_tech_vars = [
            "BuildBattery",
            "BuildPumpedHydroMW", "BuildAnyPumpedHydro",
            "BuildElectrolyzerMW", "BuildLiquifierKgPerHour", "BuildLiquidHydrogenTankKg",
            "BuildFuelCellMW",
        ]
        def no_advanced_tech_rule_factory(v):
            return lambda m, *k: (getattr(m, v)[k] == 0)
        for v in advanced_tech_vars:
            try:
                var = getattr(m, v)
                setattr(m, "PSIP_No_"+v, Constraint(var._index, rule=no_advanced_tech_rule_factory(v)))
            except AttributeError:
                pass    # model doesn't have this var

        # # don't allow any changes to the fuel market, including bulk LNG
        # # not used now; use "--force-lng-tier container" instead
        # m.PSIP_Deactivate_Limited_RFM_Supply_Tiers = Constraint(m.RFM_SUPPLY_TIERS,
        #     rule=lambda m, r, p, st:
        #         Constraint.Skip if (m.rfm_supply_tier_limit[r, p, st] == float('inf'))
        #         else (m.RFMSupplyTierActivate[r, p, st] == 0)
        # )

reason: from pyomo


text:
import cv2
import numpy as np
from copy_paste import CopyPaste
from coco import CocoDetectionCP
from visualize import display_instances
import albumentations as A
import random
from matplotlib import pyplot as plt
#%%

def xywh2xyxy(x):
    # Convert nx4 boxes from [x, y, w, h] to [x1, y1, x2, y2] where xy1=top-left, xy2=bottom-right
    y =  np.zeros_like(x)
    y[:, 0] = x[:, 0]  # top left x
    y[:, 1] = x[:, 1]  # top left y
    y[:, 2] = x[:, 0] + x[:, 2]  # bottom right x
    y[:, 3] = x[:, 1] + x[:, 3]  # bottom right y
    y[:, 4] = x[:,4]
    return y
transform = A.Compose([
    A.RandomScale(scale_limit=(-0.9, 1), p=1), #LargeScaleJitter from scale of 0.1 to 2
    A.PadIfNeeded(512, 512, border_mode=0), #pads with image in the center, not the top left like the paper
    A.RandomCrop(512, 512),
    CopyPaste(blend=True, sigma=1, pct_objects_paste=0.8, p=1.) #pct_objects_paste is a guess
], bbox_params=A.BboxParams(format="coco", min_visibility=0.05)
)

data = CocoDetectionCP(
    '../data/foof_full_coco/',
    '../data/foof_full_coco/annotations.json',
    transform
)

#%%
path='../copy-paste-aug/fake_image/'
for i in range(100):
    #f, ax = plt.subplots(1, 1, figsize=(16, 16))
    index = random.randint(0, len(data)-1)
    img_data = data[index]
    image = img_data['image']
    masks = img_data['masks']
    bboxes = img_data['bboxes']
    bbox=xywh2xyxy(np.array(bboxes)[:,:5]).astype(int)

    cv2.imwrite(path+'{}.jpg'.format(i),cv2.cvtColor(image, cv2.COLOR_BGR2RGB))
    with open('lable.txt','a',encoding='utf-8') as f:
        f.write(path+'{}.jpg'.format(i)+'')
        for i in bbox:
            f.write(" " + ",".join([str(a) for a in i[:5]]))
        f.write('\n')

    # empty = np.array([])
    # #display_instances(image, empty, empty, empty, empty, show_mask=False, show_bbox=False, ax=ax)
    # f1, ax1 = plt.subplots(1, 1, figsize=(16, 16))
    # if len(bboxes) > 0:
    #     boxes = np.stack([b[:4] for b in bboxes], axis=0)
    #     box_classes = np.array([b[-2] for b in bboxes])
    #     mask_indices = np.array([b[-1] for b in bboxes])
    #     show_masks = np.stack(masks, axis=-1)[..., mask_indices]
    #     class_names = {k: data.coco.cats[k]['name'] for k in data.coco.cats.keys()}
        #display_instances(image, boxes, show_masks, box_classes, class_names, show_bbox=True, ax=ax1)
    #else:
        #display_instances(image, empty, empty, empty, empty, show_mask=False, show_bbox=False, ax=ax1)


text:
# coding: utf-8

import numpy
import scipy
from scipy import ndimage
import numpy as np
from scipy import misc

from PIL import Image
from pyinterp2.interp2 import interp2linear

# In[36]:

import numpy
import scipy
#noinspection PyPep8Naming
import matplotlib.pyplot as plt

# In[38]:

def rgb2gray(rgb):
    #return np.dot(rgb[...,:3], [0.299, 0.587, 0.144])
    return np.dot(rgb[..., :3], [0.2989, 0.5870, 0.1140]).round()  #matlab's version

class MyStruct:
    def __init__(self):
        self.ololo = 1

#TODO: check for several sizes
#noinspection PyPep8Naming,PyPep8Naming,PyPep8Naming
def extend(I, shape):
    """
    Extends image I to shape
    """
    r = [0, 0]
    Inew = I
    for i in [0, 1]:
        r = numpy.ones(Inew.shape[i], dtype='int32')
        r[-1] = shape[i] - Inew.shape[i] + 1
        Inew = np.repeat(Inew, r, axis=i)
    return Inew

#TODO: add antialiasing to work like matlab
#noinspection PyPep8Naming
def my_resize(I, factor):
    if factor != 0.5:
        raise Exception("Unsupported resize factor!")
    newshape = numpy.round((numpy.array(I.shape) ) / 2.).astype('int32')
    #print "new shape = ", newshape
    #print I[0::2,0::2].shape
    #print I[0::2,1::2].shape
    #print I[1::2,0::2].shape
    #print I[1::2,1::2].shape
    res = (
              extend(I[0::2, 0::2], newshape) +
              extend(I[0::2, 1::2], newshape) +
              extend(I[1::2, 0::2], newshape) +
              extend(I[1::2, 1::2], newshape)) / 4
    return res

import scipy.signal

def median_filter_symmetric(x, r):
    border = r
    padded = numpy.pad(x, pad_width=border, mode='symmetric')
    filt = scipy.signal.medfilt2d(padded, r)
    cropped = filt[border:-border, border:-border]
    return cropped

# In[60]:

from scipy.sparse import coo_matrix

def makeweights_helper(abs_diff, params, ):
    epsilon = 1e-5
    vd_min = numpy.min(abs_diff)
    vd_max = numpy.max(abs_diff)
    d = vd_max - vd_min
    if d == 0:
        d = 1.

    t = (abs_diff - vd_min) / d
    weights_new = numpy.exp(-(params.beta * t)) + epsilon
    weights1 = params.nu + (1 - params.nu) * weights_new

    #return numpy.ones(abs_diff.shape)
    return weights1

def makeweights(I0_, params):
    args = dict(order='F')
    M, N = I0_.shape
    MN = M * N
    X, Y = np.meshgrid(np.arange(N), np.arange(M))
    #print X, Y
    abs_diff = numpy.abs(I0_[:, :-1] - I0_[:, 1:])
    weights = makeweights_helper(abs_diff, params)
    coords1 = numpy.vstack([Y[:, :-1].ravel(**args), X[:, :-1].ravel(**args)])
    coords2 = numpy.vstack([Y[:, 1:].ravel(**args), X[:, 1:].ravel(**args)])
    #print coords1
    edges1 = np.ravel_multi_index(coords1, I0_.shape, **args)
    edges2 = np.ravel_multi_index(coords2, I0_.shape, **args)

    D1 = coo_matrix(
        (
            numpy.hstack([-weights.ravel(**args), weights.ravel(**args)]),
            (numpy.hstack([edges1, edges1]), numpy.hstack([edges1, edges2])),
        ),
        shape=(MN, MN)
    )

    abs_diff = numpy.abs(I0_[:-1, :] - I0_[1:, :])
    weights = makeweights_helper(abs_diff, params)
    coords1 = numpy.vstack([Y[:-1, :].ravel(**args), X[:-1, :].ravel(**args)])
    coords2 = numpy.vstack([Y[1:, :].ravel(**args), X[1:, :].ravel(**args)])
    edges1 = np.ravel_multi_index(coords1, I0_.shape, **args)
    edges2 = np.ravel_multi_index(coords2, I0_.shape, **args)

    D2 = coo_matrix(
        (
            numpy.hstack([-weights.ravel(**args), weights.ravel(**args)]),
            (numpy.hstack([edges1, edges1]), numpy.hstack([edges1, edges2])),
        ),
        shape=(MN, MN)
    )

    return D1, D2,

#%% Computes Lipschitz constant
def compute_lipschitz_constant(normAtA, alpha, mu_data, mu_tv):
    L = 8 * alpha / mu_tv + normAtA / mu_data
    return L

#%% Computes the derivative of the TV term
def df_tv(x, x0, mu, Dt, D1, D2):
    #%% see Eq. 17 at Ayvaci, Raptis, Soatto, NIPS'10
    x = x + x0

    D1x = D1.dot(x)
    D2x = D2.dot(x)
    tvx = (D1x * D1x + D2x * D2x) ** .5
    #print "tvx", tvx.T
    #print "mu", mu
    w = numpy.maximum(tvx, mu)
    #print "w", w.T
    u1 = D1x / w
    u2 = D2x / w

    df = Dt.dot(numpy.vstack([u1, u2]))

    f = tvx.sum()

    return df, f

#%% Computes the derivative of the Huber-L1 norm on x
def huber_l1(x, mu):
    inds = x <= mu
    #print ~inds
    fx = numpy.zeros(x.shape)
    fx[inds] = x[inds] ** 2 / (2 * mu)
    fx[~inds] = numpy.abs(x[~inds]) - (mu / 2)
    #print "fx", fx
    f = numpy.sum(fx)

    return f

# Computes the derivative of the Huber-L1 norm on the data term |Ax + b|
def df_huber_l1_Axplusb(x, A, b, mu):
    Axplusb = A.dot(x) + b
    #print "Axplusb", Axplusb
    Axbplus_abs = numpy.abs(Axplusb)
    #print "Axbplus_abs", Axbplus_abs
    #print mu
    max_ = numpy.maximum(Axbplus_abs, mu)
    #print max_
    #print "Axplusb / max_", Axplusb / max_,
    df = A.T.dot(Axplusb / max_)
    f = huber_l1(Axplusb, mu)
    return df, f

# Computes the derivative
def compute_df(x, A, b, u0, v0, alpha, mu_data, mu_tv, Dt, D1, D2):
    MN = x.shape[0] / 2
    MN2 = 2 * MN

    df1, f1 = df_huber_l1_Axplusb(x, A, b, mu_data)
    #print "df1", df1
    df2, f2 = df_tv(x[0: MN], u0, mu_tv, Dt, D1, D2)
    #print "df2", df2, df2.shape, f2.shape
    df3, f3 = df_tv(x[MN:MN2], v0, mu_tv, Dt, D1, D2)
    #print "df3", df3
    df = df1 + alpha * numpy.vstack([df2, df3])
    f = f1 + alpha * (f2 + f3)
    #print "df2", df2.T

    return df, f

def visualize(k, x, stats, A, b, u0, v0, M, N):
    #print x.shape, u0.shape, v0.shape
    MN = M * N
    u = (x[0:MN] + u0 ).reshape([N, M]).T
    v = (x[MN:] + v0 ).reshape([N, M]).T
    plt.figure(100)
    #plt.quiver( u, v, scale_units='xy', scale=0.1, angles='xy')
    plt.quiver(u, v, scale_units='xy', angles='xy')
    plt.xlim(-1, M + 1)
    plt.ylim(-1, N + 1)

# pyramid generation and flow estimation
def Huber_L1_wTV_nesterov_core(A, b, u0_, v0_, D1, D2, M, N, params):
    MN = M * N
    MN2 = MN * 2

    At = A.T
    Atb = At.dot(b)
    AtA = At.dot(A)
    AAt = A.dot(At)

    normAtA = max(AAt.diagonal())

    D1 = D1
    D2 = D2
    Dt = scipy.sparse.vstack([D1, D2]).T

    # %% Parameters
    alpha = params.alpha
    mu_data = params.mu_data
    mu_tv = params.mu_tv

    # %% Initialize
    x0 = numpy.zeros([MN2, 1])
    xk = x0
    xold = xk

    L = compute_lipschitz_constant(normAtA, alpha, mu_data, mu_tv)

    # initialize statistics storage
    stats = MyStruct()
    stats.f = numpy.zeros([1, params.maxiters])
    stats.energy = numpy.zeros([params.maxiters])
    stats.conver = numpy.zeros([params.maxiters])

    k = 0
    iteration = 1
    stop = False
    wdf = 0

    while not stop and iteration < params.maxiters:
    #if 0:
    #    print "iter", iter

        # step (1) compute the derivative
        df, f = compute_df(xk, A, b, u0_, v0_, alpha, mu_data, mu_tv, Dt, D1, D2)
        # step (2) update yk
        yk = xk - (1 / L) * df

        # step (3) update zk
        alphak = (k + 1) / 2.
        wdf += alphak * df
        zk = x0 - (1 / L) * wdf

        # step (4) blend yk and zk
        tauk = 2. / (k + 3)
        xkp = tauk * zk + (1 - tauk) * yk
        xk = xkp

        # save statistics
        stats.energy[iteration] = f

        if iteration > 10:
            iterm10 = iteration - 10
            fbar = numpy.mean(stats.energy[iterm10:iteration])
            convergence = abs(f - fbar) / fbar
            #print "convergence", convergence
            if convergence < 1e-4:
                stop = True

        # visualize
        if params.display and (((k + 1) % 100) == 1 or stop or ((iteration + 1) == params.maxiters)):
        #if 1:
            visualize(k, xk, stats, A, b, u0_, v0_, M, N)

        xold = xk
        k += 1
        iteration += 1

    xk = xk + numpy.vstack([u0_, v0_])
    return xk

# In[116]:

#np.set_printoptions(precision=5, linewidth = 250)
##At.todense()
##print L, normAtA, alpha, params.alpha, df
#df, f = compute_df(xk, A, b, u0_, v0_, alpha, mu_data, mu_tv, Dt, D1, D2);

def Huber_L1_wTV_nesterov(I0, I1warped, u0, v0, Ix, Iy, It, params):
    M, N = I0.shape
    b = It.T.reshape([-1, 1])
    A = scipy.sparse.hstack([
        scipy.sparse.diags([Ix.T.ravel()], [0]),
        scipy.sparse.diags([Iy.T.ravel()], [0]),
        ])
    D1, D2 = makeweights(I0, params)
    x = Huber_L1_wTV_nesterov_core(A, b, u0.reshape(-1, 1, order='F'), v0.reshape(-1, 1, order='F'), D1, D2, M, N,
                                   params)
    x = x.reshape([2, M * N]).T
    u = x[:, 0].reshape([N, M]).T
    v = x[:, 1].reshape([N, M]).T
    return u, v

import sys

def construct_image_pyramid(I, pyrlevels, pyrfactor):
    #print >> sys.stderr, "Pyramid is constructed to [0,pyrlevels-1] range!!"
    factor = 2. ** .5
    smooth_sigma = (1. / pyrfactor) ** .5 / factor

    pyr = []
    tmp = I
    pyr.append(I)
    for m in range(pyrlevels - 1):
        filt1 = scipy.ndimage.filters.gaussian_filter(tmp, smooth_sigma, order=0, output=None, mode='reflect', cval=0.0,
                                                      truncate=2.0)
        #filt1r  = scipy.ndimage.zoom(filt1, pyrfactor, order=1)
        filt1r = my_resize(filt1, pyrfactor)
        tmp = filt1r
        pyr.append(filt1r)
    return pyr

def Huber_L1_wTV_nesterov_pyramid(I0, I1):
    params = MyStruct()
    pyrfactor = .5
    warps = 5
    pyrlevels = 2
    # ALPHA is the coefficient of the regularizer. When the option do_varying_alpha is selected,
    # for each pyramid level, its value varies between the values alpha0 and alphamax with the
    # multiplier alphamult at each warping step.
    params.do_varying_alpha = True
    if params.do_varying_alpha:
        params.alpha0 = 0.006
        params.alphamult = 5
        params.alphamax = 0.8
    else:
        params.alpha = 0.2

    # Thresholds for Huber-L1 norm for data term and regularizer
    params.mu_tv = 0.01
    params.mu_data = 0.01

    # Parameters of weights for gradients: w(x) = NU - (1-NU) exp(-BETA |\nabla I(x)|^2_2)
    params.beta = 30
    params.nu = 0.01

    params.maxiters = 500  # max number of iterations for each optimization loop
    params.display = True     # display results

    iscolor = False

    I0pyr = construct_image_pyramid(I0, pyrlevels, pyrfactor)
    I1pyr = construct_image_pyramid(I1, pyrlevels, pyrfactor)

    import scipy.ndimage

    #level = pyrlevels - 1
    for level in range(pyrlevels - 1, -1, -1):
        M, N = I0pyr[level].shape
        if level == pyrlevels - 1:
            u = numpy.zeros_like(I0pyr[level])
            v = numpy.zeros_like(I0pyr[level])
        else:
            # rescale motion vector field and length of the vectors
            rescale_v, rescale_u = numpy.array([M, N], dtype='float32') / I0pyr[level + 1].shape
            #rescale_u = float(N) / I0pyr[level+1].shape[1]
            #rescale_v = float(M) / I0pyr[level+1].shape[0]

            u = scipy.ndimage.zoom(u, [rescale_v, rescale_u], order=3) * rescale_u
            v = scipy.ndimage.zoom(v, [rescale_v, rescale_u], order=3) * rescale_v
            if u.shape != I0pyr[level].shape or v.shape != I0pyr[level].shape:
                raise Exception("Resize failed during transition to higher levels. Need better resize implementation.")

        u0 = u
        v0 = v
        I0_ = I0pyr[level]
        I1_ = I1pyr[level]

        idx, idy = np.meshgrid(np.arange(N), np.arange(M))   # indexing from 0 !!

        # Compute the spatial derivatives
        mask = numpy.array([1, -8, 0, 8, -1], ndmin=2) / 12.0
        Ix = ndimage.correlate(I0_, mask, mode='nearest') #Ix = imfilter(I0_, mask, 'replicate');
        Iy = ndimage.correlate(I0_, mask.T, mode='nearest') #Iy = imfilter(I0_, mask','replicate');

        if params.do_varying_alpha:
            params.alpha = params.alpha0

        for i in range(warps):
            if params.display:
                print 'Pyramid level %d, Warp %d' % (level, i)

            # Median filtering of the motion field
            u0 = median_filter_symmetric(u, 5)  #medfilt2(u, [5 5], 'symmetric');
            v0 = median_filter_symmetric(v, 5)  # medfilt2(v, [5 5], 'symmetric');
            #e0 = numpy.zeros([M,N]);

            idxx = idx + u0
            idyy = idy + v0

            I1warped = interp2linear(I1_, idxx, idyy)
            It = I1warped - I0_

            # boundary handling
            m = (idxx > N - 1) | (idxx < 0) | (idyy > M - 1) | (idyy < 0)
            Ix[m] = 0.0
            Iy[m] = 0.0
            It[m] = 0.0

            # Estimate the motion from I0 to warped I
            u, v = Huber_L1_wTV_nesterov(I0_, I1warped, u0, v0, Ix, Iy, It, params)

            if params.do_varying_alpha:
                params.alpha = min(params.alphamult * params.alpha, params.alphamax)

    return u, v

def main(argv):
    I0color = misc.imread('../data/j/j_frm_00000.png')
    I1color = misc.imread('../data/j/j_frm_00001.png')

    I0 = rgb2gray(I0color)
    I1 = rgb2gray(I1color)

    #I0 = numpy.array([ 0, 1, 0, 0, 2, 0, 0, 1, 0]).reshape(3,3) * 255 / 2.0;
    #I1 = numpy.array([ 1, 0, 0, 2, 0, 0, 1, 0, 0]).reshape(3,3) * 255 / 2.0;

    #I0 = numpy.array([ 0, 1, 0, 0, 2, 0, 0, 1, 0, 2, 0, 0]).reshape(4,3) * 255 / 2.0;
    #I1 = numpy.array([ 1, 0, 0, 2, 0, 0, 1, 0, 0, 0, 2, 0]).reshape(4,3) * 255 / 2.0;

    I0 /= 255.
    I1 /= 255.0

    u, v = Huber_L1_wTV_nesterov_pyramid(I0, I1)

    import scipy.io

    scipy.io.savemat('python_weighted__.mat', {'u': u, 'v': v})
    print "done"

if __name__ == "__main__":
    main(sys.argv)

reason: import numpy,import scipy,from scipy


text:
# Copyright (C) 2002, Thomas Hamelryck (thamelry@binf.ku.dk)
#
# This file is part of the Biopython distribution and governed by your
# choice of the "Biopython License Agreement" or the "BSD 3-Clause License".
# Please see the LICENSE file that should have been included as part of this
# package.
"""Base class for Residue, Chain, Model and Structure classes.

It is a simple container class, with list and dictionary like properties.
"""

from copy import copy

from Bio.PDB.PDBExceptions import PDBConstructionException

class Entity:
    """Basic container object for PDB heirachy.

    Structure, Model, Chain and Residue are subclasses of Entity.
    It deals with storage and lookup.
    """

    def __init__(self, id):
        """Initialize the class."""
        self._id = id
        self.full_id = None
        self.parent = None
        self.child_list = []
        self.child_dict = {}
        # Dictionary that keeps additional properties
        self.xtra = {}

    # Special methods

    def __len__(self):
        """Return the number of children."""
        return len(self.child_list)

    def __getitem__(self, id):
        """Return the child with given id."""
        return self.child_dict[id]

    def __delitem__(self, id):
        """Remove a child."""
        return self.detach_child(id)

    def __contains__(self, id):
        """Check if there is a child element with the given id."""
        return id in self.child_dict

    def __iter__(self):
        """Iterate over children."""
        yield from self.child_list

    # Generic id-based comparison methods considers all parents as well as children
    # Works for all Entities - Atoms have comparable custom operators
    def __eq__(self, other):
        """Test for equality. This compares full_id including the IDs of all parents."""
        if isinstance(other, type(self)):
            if self.parent is None:
                return self.id == other.id
            else:
                return self.full_id[1:] == other.full_id[1:]
        else:
            return NotImplemented

    def __ne__(self, other):
        """Test for inequality."""
        if isinstance(other, type(self)):
            if self.parent is None:
                return self.id != other.id
            else:
                return self.full_id[1:] != other.full_id[1:]
        else:
            return NotImplemented

    def __gt__(self, other):
        """Test greater than."""
        if isinstance(other, type(self)):
            if self.parent is None:
                return self.id > other.id
            else:
                return self.full_id[1:] > other.full_id[1:]
        else:
            return NotImplemented

    def __ge__(self, other):
        """Test greater or equal."""
        if isinstance(other, type(self)):
            if self.parent is None:
                return self.id >= other.id
            else:
                return self.full_id[1:] >= other.full_id[1:]
        else:
            return NotImplemented

    def __lt__(self, other):
        """Test less than."""
        if isinstance(other, type(self)):
            if self.parent is None:
                return self.id < other.id
            else:
                return self.full_id[1:] < other.full_id[1:]
        else:
            return NotImplemented

    def __le__(self, other):
        """Test less or equal."""
        if isinstance(other, type(self)):
            if self.parent is None:
                return self.id <= other.id
            else:
                return self.full_id[1:] <= other.full_id[1:]
        else:
            return NotImplemented

    def __hash__(self):
        """Hash method to allow uniqueness (set)."""
        return hash(self.full_id)

    # Private methods

    def _reset_full_id(self):
        """Reset the full_id (PRIVATE).

        Resets the full_id of this entity and
        recursively of all its children based on their ID.
        """
        for child in self:
            try:
                child._reset_full_id()
            except AttributeError:
                pass  # Atoms do not cache their full ids.
        self.full_id = self._generate_full_id()

    def _generate_full_id(self):
        """Generate full_id (PRIVATE).

        Generate the full_id of the Entity based on its
        Id and the IDs of the parents.
        """
        entity_id = self.get_id()
        parts = [entity_id]
        parent = self.get_parent()
        while parent is not None:
            entity_id = parent.get_id()
            parts.append(entity_id)
            parent = parent.get_parent()
        parts.reverse()
        return tuple(parts)

    # Public methods

    @property
    def id(self):
        """Return identifier."""
        return self._id

    @id.setter
    def id(self, value):
        """Change the id of this entity.

        This will update the child_dict of this entity's parent
        and invalidate all cached full ids involving this entity.

        @raises: ValueError
        """
        if value == self._id:
            return
        if self.parent:
            if value in self.parent.child_dict:
                raise ValueError(
                    "Cannot change id from `{}` to `{}`. The id `{}` is already used "
                    "for a sibling of this entity.".format(self._id, value, value)
                )
            del self.parent.child_dict[self._id]
            self.parent.child_dict[value] = self

        self._id = value
        self._reset_full_id()

    def get_level(self):
        """Return level in hierarchy.

        A - atom
        R - residue
        C - chain
        M - model
        S - structure
        """
        return self.level

    def set_parent(self, entity):
        """Set the parent Entity object."""
        self.parent = entity
        self._reset_full_id()

    def detach_parent(self):
        """Detach the parent."""
        self.parent = None

    def detach_child(self, id):
        """Remove a child."""
        child = self.child_dict[id]
        child.detach_parent()
        del self.child_dict[id]
        self.child_list.remove(child)

    def add(self, entity):
        """Add a child to the Entity."""
        entity_id = entity.get_id()
        if self.has_id(entity_id):
            raise PDBConstructionException("%s defined twice" % entity_id)
        entity.set_parent(self)
        self.child_list.append(entity)
        self.child_dict[entity_id] = entity

    def insert(self, pos, entity):
        """Add a child to the Entity at a specified position."""
        entity_id = entity.get_id()
        if self.has_id(entity_id):
            raise PDBConstructionException("%s defined twice" % entity_id)
        entity.set_parent(self)
        self.child_list[pos:pos] = [entity]
        self.child_dict[entity_id] = entity

    def get_iterator(self):
        """Return iterator over children."""
        yield from self.child_list

    def get_list(self):
        """Return a copy of the list of children."""
        return copy(self.child_list)

    def has_id(self, id):
        """Check if a child with given id exists."""
        return id in self.child_dict

    def get_parent(self):
        """Return the parent Entity object."""
        return self.parent

    def get_id(self):
        """Return the id."""
        return self.id

    def get_full_id(self):
        """Return the full id.

        The full id is a tuple containing all id's starting from
        the top object (Structure) down to the current object. A full id for
        a Residue object e.g. is something like:

        ("1abc", 0, "A", (" ", 10, "A"))

        This corresponds to:

        Structure with id "1abc"
        Model with id 0
        Chain with id "A"
        Residue with id (" ", 10, "A")

        The Residue id indicates that the residue is not a hetero-residue
        (or a water) because it has a blank hetero field, that its sequence
        identifier is 10 and its insertion code "A".
        """
        if self.full_id is None:
            self._reset_full_id()
        return self.full_id

    def transform(self, rot, tran):
        """Apply rotation and translation to the atomic coordinates.

        :param rot: A right multiplying rotation matrix
        :type rot: 3x3 Numeric array

        :param tran: the translation vector
        :type tran: size 3 Numeric array

        Examples
        --------
        This is an incomplete but illustrative example::

            from numpy import pi, array
            from Bio.PDB.vectors import Vector, rotmat
            rotation = rotmat(pi, Vector(1, 0, 0))
            translation = array((0, 0, 1), 'f')
            entity.transform(rotation, translation)

        """
        for o in self.get_list():
            o.transform(rot, tran)

    def copy(self):
        """Copy entity recursively."""
        shallow = copy(self)

        shallow.child_list = []
        shallow.child_dict = {}
        shallow.xtra = copy(self.xtra)

        shallow.detach_parent()

        for child in self.child_list:
            shallow.add(child.copy())
        return shallow

class DisorderedEntityWrapper:
    """Wrapper class to group equivalent Entities.

    This class is a simple wrapper class that groups a number of equivalent
    Entities and forwards all method calls to one of them (the currently selected
    object). DisorderedResidue and DisorderedAtom are subclasses of this class.

    E.g.: A DisorderedAtom object contains a number of Atom objects,
    where each Atom object represents a specific position of a disordered
    atom in the structure.
    """

    def __init__(self, id):
        """Initialize the class."""
        self.id = id
        self.child_dict = {}
        self.selected_child = None
        self.parent = None

    # Special methods

    def __getattr__(self, method):
        """Forward the method call to the selected child."""
        if method == "__setstate__":
            # Avoid issues with recursion when attempting deepcopy
            raise AttributeError
        if not hasattr(self, "selected_child"):
            # Avoid problems with pickling
            # Unpickling goes into infinite loop!
            raise AttributeError
        return getattr(self.selected_child, method)

    def __getitem__(self, id):
        """Return the child with the given id."""
        return self.selected_child[id]

    # XXX Why doesn't this forward to selected_child?
    # (NB: setitem was here before getitem, iter, len, sub)
    def __setitem__(self, id, child):
        """Add a child, associated with a certain id."""
        self.child_dict[id] = child

    def __contains__(self, id):
        """Check if the child has the given id."""
        return id in self.selected_child

    def __iter__(self):
        """Return the number of children."""
        return iter(self.selected_child)

    def __len__(self):
        """Return the number of children."""
        return len(self.selected_child)

    def __sub__(self, other):
        """Subtraction with another object."""
        return self.selected_child - other

    # Sorting
    # Directly compare the selected child
    def __gt__(self, other):
        """Return if child is greater than other."""
        return self.selected_child > other

    def __ge__(self, other):
        """Return if child is greater or equal than other."""
        return self.selected_child >= other

    def __lt__(self, other):
        """Return if child is less than other."""
        return self.selected_child < other

    def __le__(self, other):
        """Return if child is less or equal than other."""
        return self.selected_child <= other

    # Public methods

    def get_id(self):
        """Return the id."""
        return self.id

    def disordered_has_id(self, id):
        """Check if there is an object present associated with this id."""
        return id in self.child_dict

    def detach_parent(self):
        """Detach the parent."""
        self.parent = None
        for child in self.disordered_get_list():
            child.detach_parent()

    def get_parent(self):
        """Return parent."""
        return self.parent

    def set_parent(self, parent):
        """Set the parent for the object and its children."""
        self.parent = parent
        for child in self.disordered_get_list():
            child.set_parent(parent)

    def disordered_select(self, id):
        """Select the object with given id as the currently active object.

        Uncaught method calls are forwarded to the selected child object.
        """
        self.selected_child = self.child_dict[id]

    def disordered_add(self, child):
        """Add disordered entry.

        This is implemented by DisorderedAtom and DisorderedResidue.
        """
        raise NotImplementedError

    def is_disordered(self):
        """Return 2, indicating that this Entity is a collection of Entities."""
        return 2

    def disordered_get_id_list(self):
        """Return a list of id's."""
        # sort id list alphabetically
        return sorted(self.child_dict)

    def disordered_get(self, id=None):
        """Get the child object associated with id.

        If id is None, the currently selected child is returned.
        """
        if id is None:
            return self.selected_child
        return self.child_dict[id]

    def disordered_get_list(self):
        """Return list of children."""
        return list(self.child_dict.values())

max_stars_repo_name: fredricj/biopython

max_issues_repo_name: fredricj/biopython

max_forks_repo_name: fredricj/biopython

reason: from numpy


text:
import numpy
from math import pi
from laplace_helper import L2_rel_error

def poisson_2d(p, b, dx, dy, l2_target):
    '''Performs Jacobi relaxation

    Parameters:
    ----------
    p : 2D array of floats
        Initial guess
    b : 2D array of floats
        Source term
    dx: float
        Mesh spacing in x direction
    dy: float
        Mesh spacing in y direction
    l2_target: float
        Error target

    Returns:
    -------
    p: 2D array of float
        Distribution after relaxation
    '''

    l2_norm = 1
    iterations = 0
    l2_conv = []

    while l2_norm > l2_target:

        pd = p.copy()

        p[1:-1,1:-1] = 1/(2*(dx**2 + dy**2)) * \
                        ((pd[1:-1,2:]+pd[1:-1,:-2])*dy**2 +\
                        (pd[2:,1:-1] + pd[:-2,1:-1])*dx**2 -\
                         b[1:-1,1:-1]*dx**2*dy**2)

        # BCs are automatically enforced

        l2_norm = L2_rel_error(pd,p)
        iterations += 1
        l2_conv.append(l2_norm)

    print('Number of Jacobi iterations: {0:d}'.format(iterations))
    return p, l2_conv

def p_analytical(X,Y,L):
    ''' Computes analytical solution to Poisson problem
    Parameters:
    ----------
    X: 2D array of float
        Mesh with x components
    Y: 2D array of float
        Mesh with y components
    L: float
        Size of domain
    Returns:
    -------
    Analytical solution
    '''
    return numpy.sin(X*pi/L)*numpy.cos(Y*pi/L)

max_stars_repo_name: sergiommr/numerical-mooc

max_issues_repo_name: Anivarth/numerical-mooc

max_forks_repo_name: Anivarth/numerical-mooc


text:
from __future__ import print_function
from __future__ import division
import sys
import math
import os
import numpy

# ***** readme *****
# This code finds the distribution of the chimeric read in each cell for each gene.

def ReverseComplement(str):
    return str[::-1].replace('A', 't').replace('T', 'a').replace('G', 'c').replace('C', 'g').upper()

def chr2num(str):
    if str.startswith('chr'):
        str = str[3:]
    if str == 'X':
        return 23
    if str == 'Y':
        return 24
    return int(str)

def SolveClip(str, readlen):
    num = []
    alp = []
    res = []
    start = 0
    scount = 0
    Ndel = 0
    length = len(str)
    lastalp = -1
    for i in range(length):
        if str[i].isalpha():
            alp.append(str[i])
            num.append(int(str[lastalp + 1:i]))
            lastalp = i
    i = 0
    while i < len(alp):
        if alp[i] == 'N':
            for j in range(i):
                if alp[j] == 'S':
                    scount = 1
                    break
            del alp[i]
            Ndel += num[i]
            del num[i]
            continue
        i += 1
    lastm = -1
    i = 0
    while i < len(alp):
        if alp[i] == 'M':
            if lastm == -1:
                lastm = i
                continue
            if lastm == i - 1:
                num[i - 1] += num[i]
                del alp[i]
                del num[i]
                continue
            lastm = i
        i += 1
    numsum = [num[0]]
    for i in range(1, len(num)):
        numsum.append(numsum[-1] + num[i])
    for i in range(len(alp)):
        if alp[i] == 'M':
            if i == 0:
                res = [numsum[0]]
            elif i == len(alp) - 1:
                res.append(numsum[i - 1])
            else:
                res.append(numsum[i - 1])
                res.append(numsum[i])
    for i in range(len(alp)):
        if alp[i] == 'M':
            if i > 0:
                start = numsum[i - 1]
            break
    return [res, alp, start, 0, Ndel, scount]

def AddCount(genename, q=1.0):
    genename = genename.split(';')
    for gene in genename:
        gene = gene.split('.')[0]
        if gene in ExprDic:
            if ExprDic[gene][1] == i:
                ExprDic[gene][0][-1] += q
            else:
                for k1 in range(ExprDic[gene][1]+1, i):
                    ExprDic[gene][0].append(0)
                ExprDic[gene][0].append(q)
                ExprDic[gene][1] = i
        else:
            ExprDic[gene] = [[], i]
            for k1 in range(start, i):
                ExprDic[gene][0].append(0)
            ExprDic[gene][0].append(q)

def AddCount2(genename, q=1.0):
    genename = genename.split('.')[0]
    if genename in ChimericExprDic:
        ChimericExprDic[genename] += q
    else:
        ChimericExprDic[genename] = q

argnum = len(sys.argv)
filedir = sys.argv[1]
start = int(sys.argv[2])
last = int(sys.argv[3])
ExprDir = sys.argv[4]
HomologyResultFile = open(sys.argv[5])
FusionMatrix = {}
EncompassingMatrix = {}
FusionPos = {}
Clusters = {}
ClusterSize = {}
ClusterCount = {}
ExprDic = {}
ChimericExprDic = {}
GeneDictionary = {}
HomologyDic = {}
good_dist = 1
goodcount = 0
prefix = ''
if argnum == 7:
    prefix = sys.argv[6]
if prefix == '.':
    prefix = ''
readfile = open(filedir + '../ChiDist/' + prefix + 'FusionRead.txt', 'w')
# for finding reads near the brkpnts
CandidateList = []
CandidateVector = {}

for line in HomologyResultFile.readlines():
    if line[0] == '#':
        continue
    info = line.split('\t')
    try:
        fusionscore = float(info[2])
        chr1 = info[3]
        pos1 = info[5]
        chr2 = info[4].replace(' ', '')
        pos2 = info[6]
        homoscore = info[11]
        gccontent = [info[12], info[13]]
        cellcount = info[8]
        read1 = info[16]
        read2 = info[17].replace('\n', '').replace('\r', '')
        HomologyDic[chr1 + '\t' + chr2 + '\t' + pos1 + '\t' + pos2] = [homoscore, gccontent, cellcount, read1, read2]
    except:
        pass
HomologyResultFile.close()
for i in range(int(start), int(last) + 1):
    try:
        infile = open(filedir + '/' + str(i) + '_FusionSupport.txt')
        exprfile = open(ExprDir + '/' + str(i) + '.rpkm.txt')
        samfile = open(filedir + '/' + str(i) + '_geneanno.sam')
    except:
        continue
    foundcount = 0
    for key in CandidateVector:
        if CandidateVector[key]:
            foundcount += 1
    sys.stderr.write('Starting: ' + str(i) + '\tCandidate Size: ' + str(len(CandidateVector)) + '\tFound Size: ' + str(foundcount) + '\n')
    lastgene = ''
    totalrc = 0
    for line in exprfile.readlines():
        info = line.split('\t')
        gene = info[0]
        rc = int(info[1])
        if gene == lastgene:
            totalrc += rc
            continue
        if lastgene != '':
            AddCount(lastgene, totalrc)
        lastgene = gene
        totalrc = rc
    AddCount(lastgene, totalrc)
    lastname = ''
    lines = []
    badreadname = ''
    allsamfilelines = samfile.readlines()
    readlength = 1000
    for line in allsamfilelines:
        if line[0] == '#':
            continue
        info = line.split('\t')
        thisname = info[1]
        if thisname == badreadname:
            continue
        if thisname == lastname:
            lines.append(line)
        else:
            if len(lines) == 2:
                info1 = lines[0].split('\t')
                info2 = lines[1].split('\t')
                if info1[0] != info2[0]:
                    gene1 = info1[0]
                    gene2 = info2[0]
                    AddCount(gene1)
                    AddCount(gene2)
                    AddCount2(gene1)
                    AddCount2(gene2)
            elif len(lines) == 3:
                info1 = lines[0].split('\t')
                info2 = lines[1].split('\t')
                info3 = lines[2].split('\t')
                gene1 = info1[0].split('.')[0]
                gene2 = info2[0].split('.')[0]
                gene3 = info3[0].split('.')[0]
                if gene1 != gene2 or gene2 != gene3:
                    readlength = len(info1[10])
                    if info1[10] == info2[10] or info1[10] == ReverseComplement(info2[10]):
                        AddCount(gene1, 0.5)
                        AddCount(gene2, 0.5)
                        AddCount(gene3)
                        AddCount2(gene1, 0.5)
                        AddCount2(gene2, 0.5)
                        AddCount2(gene3)
                        readinfo1 = info1
                        readinfo2 = info2
                    elif info1[10] == info3[10] or info1[10] == ReverseComplement(info3[10]):
                        AddCount(gene1, 0.5)
                        AddCount(gene2)
                        AddCount(gene3, 0.5)
                        AddCount2(gene1, 0.5)
                        AddCount2(gene2)
                        AddCount2(gene3, 0.5)
                        readinfo1 = info1
                        readinfo2 = info3
                    elif info2[10] == info3[10] or info2[10] == ReverseComplement(info3[10]):
                        AddCount(gene1)
                        AddCount(gene2, 0.5)
                        AddCount(gene3, 0.5)
                        AddCount2(gene1)
                        AddCount2(gene2, 0.5)
                        AddCount2(gene3, 0.5)
                        readinfo1 = info3
                        readinfo2 = info2
                    else:
                        #sys.stderr.write('!!!!!' + info1[1])
                        lines = [line]
                        lastname = thisname
                        readinfo2 = [1]
                        readinfo1 = [1]  # don't run codes below, just want to read lines
            lines = [line]
            lastname = thisname

    CandidateList = []
    for line in infile.readlines():
        info = line.split('\t')
        gene1 = info[0]
        gene2 = info[1]
        chromo1 = info[4]
        chromo2 = info[5]
        encompass = int(info[2])
        if gene1 + '\t' + gene2 in EncompassingMatrix:
            if encompass > 0:
                EncompassingMatrix[gene1 + '\t' + gene2].append(encompass)
        elif gene2 + '\t' + gene1 in EncompassingMatrix:
            if encompass > 0:
                EncompassingMatrix[gene2 + '\t' + gene1].append(encompass)
        else:
            if encompass > 0:
                EncompassingMatrix[gene1 + '\t' + gene2] = [encompass]
            else:
                EncompassingMatrix[gene1 + '\t' + gene2] = []
        splitcount = int(info[3])
        if splitcount == 0:
            continue
        splitreadinfo = info[6].split(';')
        splitreadinfo = list(set(splitreadinfo))        # remove duplication
        subscore = []
        Pos = {}
        subgrp = []
        for item in splitreadinfo:
            if len(item) > 3:
                iteminfo = item.split('+')
                if not 6 <= int(iteminfo[1]) <= readlength - 6:
                    continue
                pos = iteminfo[0].split(',')
                if pos[0] + ',' + pos[1] in Pos:
                    Pos[pos[0] + ',' + pos[1]] += 1
                else:
                    Pos[pos[0] + ',' + pos[1]] = 1
                badmap = iteminfo[2].split(',')
                # editscore = max(50 - abs(50 - int(iteminfo[1])) - int(badmap[0]) - int(badmap[1]), 0)
                ssscore = 1
                if len(subgrp) == 0:
                    subgrp = [{pos[0] + ',' + pos[1]: 1}]
                    subscore.append(ssscore)
                else:
                    fin = False
                    for j in range(len(subgrp)):
                        for key in subgrp[j]:
                            grpos0 = key.split(',')[0]
                            grpos1 = key.split(',')[1]
                            if abs(int(pos[0]) - int(grpos0)) <= 3 * good_dist and abs(int(pos[1]) - int(grpos1)) <= 3 * good_dist:
                                fin = True
                                if pos[0] + ',' + pos[1] in subgrp[j]:
                                    subgrp[j][pos[0] + ',' + pos[1]] += 1
                                else:
                                    subgrp[j][pos[0] + ',' + pos[1]] = 1
                                subscore[j] += ssscore
                                break
                        if fin:
                            break
                    if not fin:
                        subgrp.append({pos[0] + ',' + pos[1]: 1})
                        subscore.append(ssscore)
        if len(subgrp) == 0:
            continue
        rawscore = []
        for j in subscore:
            if j > 0.5:
                rawscore.append(j)
            else:
                rawscore.append(0)
        rev = False
        AvePos = []
        for mdict in subgrp:  # calculate average pos for each cluster
            left = 0
            right = 0
            cc = 0
            for key in mdict:
                pos = key.split(',')
                left += int(pos[0]) * mdict[key]
                right += int(pos[1]) * mdict[key]
                cc += mdict[key]
            left = int(left / cc)
            right = int(right / cc)
            AvePos.append([left, right, cc])
        if gene1 + '\t' + gene2 in FusionMatrix:
            added = list(numpy.zeros(len(FusionMatrix[gene1 + '\t' + gene2])))
            for j in range(len(subgrp)):
                dist = []
                for k in FusionPos[gene1 + '\t' + gene2][1]:
                    d1 = abs(k[1][0] - AvePos[j][0])
                    d2 = abs(k[1][1] - AvePos[j][1])
                    if d1 <= good_dist and d2 <= good_dist:
                        dist.append(d1 + d2)
                    else:
                        dist.append(10 * good_dist)
                if min(dist) >= 3 * good_dist:
                    FusionPos[gene1 + '\t' + gene2][1].append([subgrp[j], AvePos[j], [i]])
                    FusionMatrix[gene1 + '\t' + gene2].append([rawscore[j]])
                    added.append(0)
                else:
                    smallindex = dist.index(min(dist))
                    if added[smallindex] == 0:
                        FusionMatrix[gene1 + '\t' + gene2][smallindex].append(rawscore[j])
                        added[smallindex] = 1
                    else:
                        FusionMatrix[gene1 + '\t' + gene2][smallindex][-1] += (rawscore[j])
                    old = FusionPos[gene1 + '\t' + gene2][1][smallindex][1]
                    newcc = AvePos[j][2] + old[2]
                    newleft = int((old[2] * old[0] + AvePos[j][2] * AvePos[j][0]) / newcc)
                    newright = int((old[2] * old[1] + AvePos[j][2] * AvePos[j][1]) / newcc)
                    FusionPos[gene1 + '\t' + gene2][1][smallindex][1] = [newleft, newright, newcc]
                    FusionPos[gene1 + '\t' + gene2][1][smallindex][2].append(i)
                if not [gene1, gene2, chromo1, chromo2, AvePos[j][0], AvePos[j][1]] in CandidateList and not [gene2, gene1, chromo2, chromo1, AvePos[j][1], AvePos[j][0]] in CandidateList:
                    CandidateList.append([gene1, gene2, chromo1, chromo2, AvePos[j][0], AvePos[j][1]])
            gene = gene1 + '\t' + gene2
        elif gene2 + '\t' + gene1 in FusionMatrix:
            added = list(numpy.zeros(len(FusionMatrix[gene2 + '\t' + gene1])))
            for j in range(len(subgrp)):
                dist = []
                for k in FusionPos[gene2 + '\t' + gene1][1]:
                    d1 = abs(k[1][0] - AvePos[j][1])
                    d2 = abs(k[1][1] - AvePos[j][0])
                    if d1 <= good_dist and d2 <= good_dist:
                        dist.append(d1 + d2)
                    else:
                        dist.append(10 * good_dist)
                if min(dist) >= 3 * good_dist:
                    FusionPos[gene2 + '\t' + gene1][1].append(
                        [subgrp[j], [AvePos[j][1], AvePos[j][0], AvePos[j][2]], [i]])
                    FusionMatrix[gene2 + '\t' + gene1].append([rawscore[j]])
                    added.append(0)
                else:
                    smallindex = dist.index(min(dist))
                    if added[smallindex] == 0:
                        FusionMatrix[gene2 + '\t' + gene1][smallindex].append(rawscore[j])
                        added[smallindex] = 1
                    else:
                        FusionMatrix[gene2 + '\t' + gene1][smallindex][-1] += (rawscore[j])
                    old = FusionPos[gene2 + '\t' + gene1][1][smallindex][1]
                    newcc = AvePos[j][2] + old[2]
                    newleft = int((old[2] * old[0] + AvePos[j][2] * AvePos[j][1]) / newcc)
                    newright = int((old[2] * old[1] + AvePos[j][2] * AvePos[j][0]) / newcc)
                    FusionPos[gene2 + '\t' + gene1][1][smallindex][1] = [newleft, newright, newcc]
                    FusionPos[gene2 + '\t' + gene1][1][smallindex][2].append(i)
                if not [gene1, gene2, chromo1, chromo2, AvePos[j][0], AvePos[j][1]] in CandidateList and not [gene2, gene1, chromo2, chromo1, AvePos[j][1], AvePos[j][0]] in CandidateList:
                    CandidateList.append([gene1, gene2, chromo1, chromo2, AvePos[j][0], AvePos[j][1]])
            gene = gene2 + '\t' + gene1
            rev = True
        else:
            FusionMatrix[gene1 + '\t' + gene2] = []
            for item in rawscore:
                FusionMatrix[gene1 + '\t' + gene2].append([item])
            FusionPos[gene1 + '\t' + gene2] = [[chromo1, chromo2], []]
            for k in range(len(subgrp)):
                FusionPos[gene1 + '\t' + gene2][1].append([subgrp[k], AvePos[k], [i]])
                if not [gene1, gene2, chromo1, chromo2, AvePos[k][0], AvePos[k][1]] in CandidateList and not [gene2, gene1, chromo2, chromo1, AvePos[k][1], AvePos[k][0]] in CandidateList:
                    CandidateList.append([gene1, gene2, chromo1, chromo2, AvePos[k][0], AvePos[k][1]])
            gene = gene1 + '\t' + gene2

    for l in range(len(CandidateList)):
        try:
            aa = chr2num(CandidateList[l][2])
            aa = chr2num(CandidateList[l][3])
            thischr1 = chr2num(CandidateList[l][2])
            thischr2 = chr2num(CandidateList[l][3])
        except:
            pass
        pos1 = int(CandidateList[l][4])
        pos2 = int(CandidateList[l][5])
        if thischr1 > thischr2 or (thischr1 == thischr2 and pos2 > pos1):
            thischr1, thischr2 = thischr2, thischr1
            pos1, pos2 = pos2, pos1
        if not (thischr1, thischr2, pos1, pos2) in CandidateVector:
            CandidateVector[(thischr1, thischr2, pos1, pos2)] = []

    for line in allsamfilelines:
        if line[0] == '#':
            continue
        info = line.split('\t')
        thisname = info[1]
        if thisname == badreadname:
            continue
        if thisname == lastname:
            lines.append(line)
        else:
            if len(lines) == 2:
                info1 = lines[0].split('\t')
                info2 = lines[1].split('\t')
            elif len(lines) == 3:
                info1 = lines[0].split('\t')
                info2 = lines[1].split('\t')
                info3 = lines[2].split('\t')
                gene1 = info1[0].split('.')[0]
                gene2 = info2[0].split('.')[0]
                gene3 = info3[0].split('.')[0]
                if gene1 != gene2 or gene2 != gene3:
                    if info1[10] == info2[10] or info1[10] == ReverseComplement(info2[10]):
                        readinfo1 = info1
                        readinfo2 = info2
                    elif info1[10] == info3[10] or info1[10] == ReverseComplement(info3[10]):
                        readinfo1 = info1
                        readinfo2 = info3
                    elif info2[10] == info3[10] or info2[10] == ReverseComplement(info3[10]):
                        readinfo1 = info3
                        readinfo2 = info2
                    else:
                        # sys.stderr.write('!!!!!' + info1[1])
                        lines = [line]
                        lastname = thisname
                        readinfo2 = [1]
                        readinfo1 = [1]  # don't run codes below, just want to read lines
                    if len(readinfo1) > 1 and len(gene1) > 1 and len(gene2) > 1 and len(gene3) > 1:  # normal conditions
                        gene1 = readinfo1[0]
                        gene2 = readinfo2[0]
                        chromo1 = readinfo1[3]
                        chromo2 = readinfo2[3]
                        if not chromo1.startswith('chr'):
                            chromo1 = 'chr' + chromo1
                        if not chromo2.startswith('chr'):
                            chromo2 = 'chr' + chromo2
                        clip1 = readinfo1[6]
                        clip2 = readinfo2[6]
                        readlen = len(readinfo1[10])
                        clipsplit1 = SolveClip(clip1, readlen)
                        clipsplit2 = SolveClip(clip2, readlen)
                        cc = False
                        splitpnt1 = -1
                        splitpnt2 = -1
                        for i in range(len(clipsplit1[0])):
                            for j in range(len(clipsplit2[0])):
                                if clipsplit2[0][j] == clipsplit1[0][i] or clipsplit2[0][j] + clipsplit1[0][i] == readlen:
                                    splitpnt1 = clipsplit1[0][i]
                                    splitpnt2 = clipsplit2[0][j]
                                    cc = True
                                    break
                            if cc:
                                break
                        if readlen - 0.5 > splitpnt1 > -0.5 and readlen - 0.5 > splitpnt2 > -0.5:
                            if clipsplit1[1][i] == 'S':
                                brkpnt1 = int(readinfo1[4])
                                direct1 = '+'
                            else:
                                brkpnt1 = splitpnt1 - clipsplit1[2] + int(readinfo1[4]) + clipsplit1[4] - 1
                                direct1 = '-'
                            if clipsplit2[1][j] == 'S':
                                brkpnt2 = int(readinfo2[4])
                                direct2 = '+'
                            else:
                                brkpnt2 = splitpnt2 - clipsplit2[2] + int(readinfo2[4]) + clipsplit2[4] - 1
                                direct2 = '-'
                            try:
                                aa = chr2num(chromo1)
                                aa = chr2num(chromo2)
                                thisposlist = [chr2num(chromo1), chr2num(chromo2), brkpnt1, brkpnt2]
                            except:
                                thisposlist = [chromo1, chromo2, brkpnt1, brkpnt2]
                            exchange = False
                            if thisposlist[0] > thisposlist[1] or (thisposlist[0] == thisposlist[1] and thisposlist[3] > thisposlist[2]):
                                thisposlist = [thisposlist[1], thisposlist[0], thisposlist[3], thisposlist[2]]
                                exchange = True
                            found = False
                            for k1 in range(-3, 4):
                                for k2 in range(-3, 4):
                                    if (thisposlist[0], thisposlist[1], thisposlist[2] + k1, thisposlist[3] + k2) in CandidateVector:
                                        found = True
                                        thisposindex = (thisposlist[0], thisposlist[1], thisposlist[2] + k1, thisposlist[3] + k2)
                            if found:
                                extractedRead = ''
                                if CandidateVector[thisposindex]:
                                    if CandidateVector[thisposindex][1] != 30:
                                        dist = min(CandidateVector[thisposindex][1], 60 - CandidateVector[thisposindex][1])
                                        if 30 > splitpnt1 > dist:
                                            extractedRead = readinfo1[10][:60]
                                            inbrkpnt = splitpnt1
                                        if 30 <= splitpnt1 <= readlen - 30:
                                            extractedRead = readinfo1[10][splitpnt1 - 30:splitpnt1 + 30]
                                            inbrkpnt = 30
                                        if readlen - dist > splitpnt1 > readlen - 30:
                                            extractedRead = readinfo1[10][readlen - 60:]
                                            inbrkpnt = splitpnt1 - readlen + 60
                                        if extractedRead != '':
                                            if direct1 == '+':
                                                extractedRead = ReverseComplement(extractedRead)
                                                inbrkpnt = 60 - inbrkpnt
                                            if exchange:
                                                extractedRead = ReverseComplement(extractedRead)
                                                inbrkpnt = 60 - inbrkpnt
                                            CandidateVector[thisposindex][0] = extractedRead
                                            CandidateVector[thisposindex][1] = inbrkpnt
                                else:
                                    if splitpnt1 < 30:
                                        extractedRead = readinfo1[10][:60]
                                        inbrkpnt = splitpnt1
                                    if 30 <= splitpnt1 <= readlen - 30:
                                        extractedRead = readinfo1[10][splitpnt1 - 30:splitpnt1 + 30]
                                        inbrkpnt = 30
                                    if splitpnt1 > readlen - 30:
                                        extractedRead = readinfo1[10][readlen - 60:]
                                        inbrkpnt = splitpnt1 - readlen + 60
                                    if direct1 == '+':
                                        extractedRead = ReverseComplement(extractedRead)
                                        inbrkpnt = 60 - inbrkpnt
                                    if exchange:
                                        extractedRead = ReverseComplement(extractedRead)
                                        inbrkpnt = 60 - inbrkpnt
                                        direct1, direct2 = direct2, direct1
                                    CandidateVector[thisposindex] = [extractedRead, inbrkpnt, direct1, direct2]
            lines = [line]
            lastname = thisname
    samfile.close()
    infile.close()
    exprfile.close()
for genename in ExprDic:
    if ExprDic[genename][1] < last:
        for k in range(ExprDic[genename][1]+1, last+1):
            ExprDic[genename][0].append(0)
for gene in FusionMatrix:
    temp = FusionPos[gene][1]
    for i in range(len(temp)):
        keylist = sorted(temp[i][0], key=temp[i][0].__getitem__, reverse=True)
        possum = 0
        countlist = []
        for key in keylist:
            possum += temp[i][0][key]
            countlist.append(int(temp[i][0][key]))
        if possum > 0 and max(countlist) / possum >= 0.4:
            for key in temp[i][0]:
                if key == keylist[0]:
                    pos0 = int(keylist[0].split(',')[0])
                    pos1 = int(keylist[0].split(',')[1])
                    if abs(pos0 - temp[i][1][0]) + abs(pos1 - temp[i][1][1]) < abs(pos1 - temp[i][1][0]) + abs(
                            pos0 - temp[i][1][1]):
                        FusionPos[gene][1][i].append([pos0, pos1])
                    else:
                        FusionPos[gene][1][i].append([pos1, pos0])
                    break
        else:
            FusionPos[gene][1][i].append([-1, -1])
templines = []
for gene in FusionMatrix:
    for i in range(len(FusionPos[gene][1])):
        if FusionPos[gene][1][i][3][0] > -0.5:
            genepart = gene.split('\t')
            usecell = FusionPos[gene][1][i][2]
            usecell = list(set(usecell))
            try:
                aaaaa = []
                thisgenename = genepart[0].split('.')[0]
                for cell in usecell:
                    aaaaa.append(ExprDic[thisgenename][0][cell-start])
            except:
                aaaaa = list(numpy.zeros(last-start+1, int))
                continue
            try:
                bbbbb = []
                thisgenename = genepart[1].split('.')[0]
                for cell in usecell:
                    bbbbb.append(ExprDic[thisgenename][0][cell-start])
            except:
                bbbbb = list(numpy.zeros(last-start+1, int))
                continue
            if FusionPos[gene][0][0] + '\t' + FusionPos[gene][0][1] + '\t' + str(FusionPos[gene][1][i][3][0]) + '\t' \
                    + str(FusionPos[gene][1][i][3][1]) in HomologyDic:
                [homoscore, gccontent, cellcount, read1, read2] = HomologyDic[FusionPos[gene][0][0] + '\t' + FusionPos[gene][0][1]
                                                                + '\t' + str(FusionPos[gene][1][i][3][0]) + '\t' + str(
                    FusionPos[gene][1][i][3][1])]
            elif FusionPos[gene][0][1] + '\t' + FusionPos[gene][0][0] + '\t' + str(FusionPos[gene][1][i][3][1]) + '\t' \
                    + str(FusionPos[gene][1][i][3][0]) in HomologyDic:
                [homoscore, gccontent, cellcount, read1, read2] = HomologyDic[FusionPos[gene][0][1] + '\t' + FusionPos[gene][0][0]
                                                                + '\t' + str(FusionPos[gene][1][i][3][1]) + '\t' + str(
                    FusionPos[gene][1][i][3][0])]
            else:
                continue
            engene = gene
            if gene not in EncompassingMatrix:
                engene = genepart[1] + '\t' + genepart[0]
            templine = genepart[0] + '\t' + genepart[1] + '\t' + \
                       str(len(FusionMatrix[gene][i])) + '\t' + str(FusionMatrix[gene][i]) + '\t' + \
                       str(len(EncompassingMatrix[engene])) + '\t' + str(EncompassingMatrix[engene]) + '\t' + \
                       FusionPos[gene][0][0] + '\t' + FusionPos[gene][0][1] + '\t' + str(FusionPos[gene][1][i][3][0]) + \
                       '\t' + str(FusionPos[gene][1][i][3][1]) + '\t' + str(i) + '\t' + homoscore + '\t' + \
                       gccontent[0] + '\t' + gccontent[1] + '\t' + str(aaaaa) + '\t' + str(bbbbb) + '\t' + read1 + '\t' + read2
            templines.append(templine)
uselines = []
for i in range(len(templines)):
    info = templines[i].split('\t')
    try:
        aa = chr2num(info[6])
        aa = chr2num(info[7])
        thischr1 = chr2num(info[6])
        thischr2 = chr2num(info[7])
    except:
        pass
    thispos1 = int(info[8])
    thispos2 = int(info[9])
    exchange = False
    if thischr1 > thischr2 or (thischr1 == thischr2 and thispos2 > thispos1):
        thischr1, thischr2 = thischr2, thischr1
        thispos1, thispos2 = thispos2, thispos1
        exchange = True
    found = False
    for k1 in range(-3, 4):
        for k2 in range(-3, 4):
            if (thischr1, thischr2, thispos1 + k1, thispos2 + k2) in CandidateVector:
                found = True
                thisposindex = (thischr1, thischr2, thispos1 + k1, thispos2 + k2)
    if found:
        if CandidateVector[thisposindex]:
            if exchange:
                uselines.append(templines[i].rstrip() + '\t' + ReverseComplement(CandidateVector[thisposindex][
                    0]) + '\t' + str(60-CandidateVector[thisposindex][1]) + '\t' + CandidateVector[thisposindex][3] + '\t' +
                                CandidateVector[thisposindex][2])
            else:
                uselines.append(templines[i].rstrip() + '\t' + CandidateVector[thisposindex][
                    0] + '\t' + str(CandidateVector[thisposindex][1]) + '\t' + CandidateVector[thisposindex][2] + '\t' +
                                CandidateVector[thisposindex][3])

linedic = {}
for line in uselines:
    info = line.split('\t')
    score = len(info[3])
    linedic[line] = score
for key in sorted(linedic, key=linedic.__getitem__, reverse=True):
    info = key.split('\t')
    read = info[-4]
    splitpos = info[-3]
    chromo1 = info[6]
    chromo2 = info[7]
    pos1 = info[8]
    pos2 = info[9]
    if read.find('N') > -1:
        continue
    if len(read) == 60 or True:
        readfile.write(read + '\t' + splitpos + '\t' + chromo1 + ':' + pos1 + ':' + info[-2] + '\t' + chromo2 + ':' +
                       pos2 + ':' + info[-1] + '\n')
        print(key)
readfile.close()

max_stars_repo_name: ZijieJin/scFusion

max_issues_repo_name: ZijieJin/scFusion

max_forks_repo_name: ZijieJin/scFusion