text:
import asyncio
import os
import time
from typing import List

import numpy as np
import pytest

from aqueduct.exceptions import (
    FlowError,
    NotRunningError,
)
from aqueduct.flow import (
    Flow,
    FlowState,
    FlowStep,
)
from aqueduct.handler import BaseTaskHandler
from aqueduct.metrics import MAIN_PROCESS
from aqueduct.task import BaseTask
from tests.unit.conftest import (
    Task,
    run_flow,
    terminate_worker,
)
from tests.unit.test_shm import ArrayFieldTask

TASKS_BATCH_SIZE = 16

class HandlerWithError(BaseTaskHandler):
    def handle(self, *tasks: BaseTask):
        raise ValueError

@pytest.fixture
async def flow_with_failed_handler(loop):
    async with run_flow(Flow(HandlerWithError())) as flow:
        yield flow

class CatDetector:
    """GPU model emulator that predicts the presence of the cat in the image."""
    IMAGE_PROCESS_TIME = 0.01
    BATCH_REDUCTION_FACTOR = 0.7
    OVERHEAD_TIME = 0.02
    BATCH_PROCESS_TIME = IMAGE_PROCESS_TIME * TASKS_BATCH_SIZE * BATCH_REDUCTION_FACTOR + OVERHEAD_TIME

    def predict(self, images: np.array) -> np.array:
        """Always says that there is a cat in the image.

        The image is represented by a one-dimensional array.
        The model spends less time for processing batch of images due to GPU optimizations. It's emulated
        with BATCH_REDUCTION_FACTOR coefficient.
        """
        batch_size = images.shape[0]
        if batch_size == 1:
            time.sleep(self.IMAGE_PROCESS_TIME)
        else:
            time.sleep(self.IMAGE_PROCESS_TIME * batch_size * self.BATCH_REDUCTION_FACTOR)
        return np.ones(batch_size, dtype=bool)

class CatDetectorHandler(BaseTaskHandler):
    def handle(self, *tasks: ArrayFieldTask):
        images = np.array([task.array for task in tasks])
        predicts = CatDetector().predict(images)
        for task, predict in zip(tasks, predicts):
            task.result = predict

@pytest.fixture
async def flow_with_not_batch_handler(loop) -> Flow:
    """Actually the handler has batching, but batch size is one."""
    async with run_flow(Flow(CatDetectorHandler())) as flow:
        yield flow

@pytest.fixture
async def flow_with_batch_handler(loop) -> Flow:
    async with run_flow(Flow(
            FlowStep(
                CatDetectorHandler(),
                batch_size=TASKS_BATCH_SIZE,
            ))) as flow:
        yield flow

async def process_tasks(flow: Flow, tasks: List[ArrayFieldTask]):
    await asyncio.gather(*(flow.process(task) for task in tasks))

@pytest.fixture
def tasks_batch(array):
    return [ArrayFieldTask(array) for _ in range(TASKS_BATCH_SIZE)]

class ShareArrayHandler(BaseTaskHandler):
    def __init__(self, array: np.array):
        super().__init__()
        self._array = array

    def handle(self, *tasks: BaseTask):
        for task in tasks:
            task.array = self._array
            task.share_value('array')

class UpdateSharedArrayFieldHandler(BaseTaskHandler):
    def handle(self, *tasks: BaseTask):
        for task in tasks:
            task.array[0] = task.array[-1]  # noqa

@pytest.fixture
async def flow_for_shared_array(loop, array) -> Flow:
    async with run_flow(Flow(ShareArrayHandler(array), UpdateSharedArrayFieldHandler())) as flow:
        yield flow

class MultiProcBatchHandler(BaseTaskHandler):
    def handle(self, *tasks: Task):
        for task in tasks:
            task.result = os.getpid()

@pytest.fixture
async def flow_with_multiproc_bathcer(loop):
    async with run_flow(
            Flow(FlowStep(MultiProcBatchHandler(), nprocs=2, batch_size=TASKS_BATCH_SIZE))
    ) as flow:
        yield flow

class SlowStartingHandler(BaseTaskHandler):
    start_time: int

    def on_start(self):
        time.sleep(self.start_time)

    def handle(self, task: Task):
        task.result = os.getpid()

class TestFlow:
    async def test_start_waiting_for_subproccess_to_start(self, task: Task):
        """flow.start() should block until all subprocesses fully initialized."""
        handler = SlowStartingHandler()
        handler.start_time = 1
        slowest_handler = SlowStartingHandler()
        slowest_handler.start_time = 3
        flow = Flow(FlowStep(handler, nprocs=3), FlowStep(slowest_handler, nprocs=3))

        t0 = time.time()
        flow.start()
        t1 = time.time()

        assert t1 - t0 == pytest.approx(slowest_handler.start_time, 0.5)

    async def test_start_waiting_timeout(self, task: Task):
        """flow.start(timeout) should raise FlowError when starting_timeout expired."""
        handler = SlowStartingHandler()
        handler.start_time = 10
        flow = Flow(FlowStep(handler, nprocs=3))

        timeout = 3
        t0 = time.time()
        with pytest.raises(TimeoutError):
            flow.start(timeout=timeout)
        t1 = time.time()

        assert t1 - t0 == pytest.approx(timeout, 0.5)

    async def test_process_set_result(self, simple_flow: Flow, task: Task):
        assert task.result is None
        result = await simple_flow.process(task)
        assert result
        assert task.result == 'test'

    def test_get_queues_info(self, simple_flow: Flow):
        info = simple_flow._get_queues_info()
        queues = simple_flow._queues
        assert len(info) == len(queues)
        assert MAIN_PROCESS in info[queues[0]]
        assert MAIN_PROCESS in info[queues[-1]]

    async def test_process_expired_task(self, log_file, aqueduct_logger,
                                        slow_sleep_handlers, slow_simple_flow, task):
        timeouts = [handler._handler_sec for handler in slow_sleep_handlers]

        with pytest.raises(FlowError, match='timeout'):
            await slow_simple_flow.process(task, timeout_sec=timeouts[0])
        await asyncio.sleep(sum(timeouts))
        log_str = log_file.read().decode()

        # task is processed by first handler
        assert f'[{slow_sleep_handlers[0].__class__.__name__}] Have message' in log_str
        # task arrives in the second handler, but it is not processed because it has expired
        assert f'[{slow_sleep_handlers[1].__class__.__name__}] Task expired. Skip' in log_str
        # doesn't arrive in the third handler
        assert f'[{slow_sleep_handlers[2].__class__.__name__}] Have message' not in log_str

    async def test_process_not_running_flow(self, simple_flow, task):
        await simple_flow.stop(graceful=False)
        with pytest.raises(NotRunningError):
            await simple_flow.process(task)

    async def test_process_huge_count_of_tasks(self, simple_flow):
        """Count of tasks much more than queue size."""
        # queue size = 20
        n_tasks = 90

        coros = [simple_flow.process(Task()) for _ in range(n_tasks)]
        res = await asyncio.gather(*coros)
        assert len(res) == n_tasks

    async def test_stopped_flow_due_to_terminated_step_process(self, simple_flow: Flow, caplog):
        """Checks that flow and event loop were stopped when child Step process was terminated."""
        assert simple_flow.state == FlowState.RUNNING

        process, handler = await terminate_worker(simple_flow)

        assert simple_flow.state == FlowState.STOPPED
        assert f'The process {process.pid} for {handler.__class__.__name__} handler is dead' in caplog.text

    async def test_stopped_flow_due_to_failed_step_process(self, flow_with_failed_handler, task):
        """Checks that flow and event loop were stopped when child Step process finished with error."""
        assert flow_with_failed_handler.state == FlowState.RUNNING

        with pytest.raises(FlowError):
            await flow_with_failed_handler.process(task)
        await asyncio.sleep(1)

        assert flow_with_failed_handler.state == FlowState.STOPPED

    async def test_stopped_flow_due_to_stop_method(self, simple_flow):
        await simple_flow.stop(graceful=False)
        await asyncio.sleep(0.5)

        assert simple_flow.state == FlowState.STOPPED

    async def test_stop_finish_task_on_graceful(self, sleep_handlers, simple_flow, task):
        """Checks that task will be processed to the end after Flow stop command with graceful mode."""
        t = asyncio.ensure_future(simple_flow.process(task))
        await asyncio.sleep(sleep_handlers[0]._handler_sec)
        assert t.done() is False
        await simple_flow.stop()
        await asyncio.sleep(0.1)

        assert t.result()
        assert simple_flow.state == FlowState.STOPPED

    async def test_stop_abort_task(self, sleep_handlers, simple_flow, task):
        """Checks that task will be aborted after Flow stop command without graceful mode."""
        t = asyncio.ensure_future(simple_flow.process(task))
        await asyncio.sleep(sleep_handlers[0]._handler_sec)
        await simple_flow.stop(graceful=False)
        await asyncio.sleep(0.1)

        assert isinstance(t.exception(), FlowError)
        assert simple_flow.state == FlowState.STOPPED

    async def test_process_performance_without_batching(self, flow_with_not_batch_handler, tasks_batch):
        """Checks that will be received TimeoutError with not batch handler."""
        with pytest.raises(asyncio.TimeoutError):
            await asyncio.wait_for(
                process_tasks(flow_with_not_batch_handler, tasks_batch),
                timeout=CatDetector.BATCH_PROCESS_TIME,
            )

    async def test_process_performance_with_batching(self, flow_with_batch_handler, tasks_batch):
        """Checks that all tasks will be processed on time with batch handler."""
        await asyncio.wait_for(
            process_tasks(flow_with_batch_handler, tasks_batch),
            timeout=CatDetector.BATCH_PROCESS_TIME,
        )
        assert all(task.result for task in tasks_batch)

    async def test_process_seq_batch_filling(self, tasks_batch, flow_with_multiproc_bathcer):
        """Checks that batches are being filled sequentially."""
        await asyncio.gather(*[flow_with_multiproc_bathcer.process(task) for task in tasks_batch])
        assert len(set(task.result for task in tasks_batch)) == 1

async def test_segfault_due_to_update_shared_fields(array, array_sh_data, simple_flow):
    t1, t2 = BaseTask(), BaseTask()
    t1.sh_array = array_sh_data
    t2.sh_array = array_sh_data

    await simple_flow.process(t1)
    await simple_flow.process(t2)

    # segfault during access to t2.sh_array due to updated `_shared_fields` field
    assert all(t2.sh_array == array)
    # segfault during access to t1.sh_array due to updated `sh_array` field
    assert all(t1.sh_array == array)

async def test_segfault_due_to_share_field_in_handler(array, flow_for_shared_array):
    t = BaseTask()
    await flow_for_shared_array.process(t)
    assert t.array[0] == array[-1]  # noqa

max_stars_repo_name: artemcpp/aqueduct

max_issues_repo_name: artemcpp/aqueduct

max_forks_repo_name: artemcpp/aqueduct


text:
import json
import numpy as np
from pathlib import Path

# Converts `val` into a bitstring that is `bw` characters wide
def bitstring(val, bw):
    # first truncate val by `bw` bits
    val &= 2 ** bw - 1
    return "{:x}".format(val)

def parse_dat(path):
    with path.open("r") as f:
        lines = f.readlines()
        arr = np.array(list(map(lambda v: int(v.strip(), 16), lines)))
        return arr

# go through the json data and create a file for each key,
# flattening the data and then converting it to bitstrings
def convert2dat(output_dir, data, extension):
    output_dir = Path(output_dir)
    shape = {}
    for k, item in data.items():
        path = output_dir / f"{k}.{extension}"
        path.touch()
        arr = np.array(item["data"])
        shape[k] = {"shape": list(arr.shape), "bitwidth": item["bitwidth"]}
        with path.open("w") as f:
            for v in arr.flatten():
                f.write(bitstring(v, item["bitwidth"]) + "\n")

    # commit shape.json file
    shape_json_file = output_dir / "shape.json"
    with shape_json_file.open("w") as f:
        json.dump(shape, f, indent=2)

# converts a directory of *.dat files back into a json file
def convert2json(input_dir, extension):
    input_dir = Path(input_dir)
    data = {}
    shape_json_path = input_dir / "shape.json"
    shape_json = None
    if shape_json_path.exists():
        shape_json = json.load(shape_json_path.open("r"))

    # TODO: change to use shape json
    for f in input_dir.glob(f"*.{extension}"):
        arr = parse_dat(f)
        if (
            shape_json is not None
            and shape_json.get(f.stem) is not None
            and shape_json[f.stem]["shape"] != [0]
        ):
            try:
                arr = arr.reshape(tuple(shape_json[f.stem]["shape"]))
            except Exception:
                raise Exception(f.stem)
            name = f.stem
            data[name] = arr.tolist()
    return data

ext: py

max_stars_repo_name: pbonh/calyx

max_issues_repo_name: pbonh/calyx

max_forks_repo_name: pbonh/calyx

reason: import numpy


text:
import os
from pathlib import Path
import pandas as pd
import numpy as np
from skimage import io
from matplotlib import pyplot as plt

data_root_dir = './data'
label_table_name = 'label_color'
image_dir = 'CamSeq01'
image_dir_path = Path(data_root_dir, image_dir)

label_table = pd.read_csv(Path(data_root_dir, label_table_name))
label_image_file_names =\
    [f for f in os.listdir(image_dir_path) if '_L.png' in f]

for image_file_name in label_image_file_names:
    file_path = Path(image_dir_path, image_file_name)
    print('processing {}'.format(file_path))
    img = io.imread(file_path)
    # plt.figure()
    # io.imshow(img)
    # plt.show()

    class_map = np.empty((img.shape[0], img.shape[1]), dtype='i')
    for i_row in range(len(label_table)):
        class_color = label_table.iloc[i_row]
        idx = np.where(
            (img[:, :, 0] == class_color[0])  # r
            & (img[:, :, 1] == class_color[1])  # g
            & (img[:, :, 2] == class_color[2]))  # b
        class_map[idx] = i_row
    label_file_name = Path(image_dir_path,
                           os.path.splitext(image_file_name)[0]+'.npz')
    np.savez(label_file_name, data=class_map)

max_issues_repo_name: universebang/pytorch-EncNet

max_forks_repo_name: universebang/pytorch-EncNet


text:
"""
Some utility functions for the map2parameters project.
"""
import numpy as np, healpy as hp
from pixell import curvedsky

def symmetrize(mat):
    """Make a upper diagonal or lower diagonal matrix symmetric

    Parameters
    ----------
    mat : 2d array
      the matrix we want symmetric
    """
    return mat + mat.T - np.diag(mat.diagonal())

def get_noise_matrix_spin0and2(noise_data_dir, exp, freqs, lmax, nsplits, lcut=0):

    """This function uses the noise power spectra computed by 'maps_to_params_prepare_sim_data'
    and generate a three dimensional array of noise power spectra [nfreqs,nfreqs,lmax] for temperature
    and polarisation.
    The different entries ([i,j,:]) of the arrays contain the noise power spectra
    for the different frequency channel pairs.
    for example nl_array_t[0,0,:] =>  nl^{TT}_{f_{0},f_{0}),  nl_array_t[0,1,:] =>  nl^{TT}_{f_{0},f_{1})
    this allows to have correlated noise between different frequency channels.

    Parameters
    ----------
    noise_data_dir : string
      the folder containing the noise power spectra
    exp : string
      the experiment to consider ('LAT' or 'Planck')
    freqs: 1d array of string
      the frequencies we consider
    lmax: integer
      the maximum multipole for the noise power spectra
    n_splits: integer
      the number of data splits we want to simulate
      nl_per_split= nl * n_{splits}
    lcut: integer
      the noise of SO being very red we will cut out the lowest modes to avoid
      leakage
    """

    nfreqs = len(freqs)
    nl_array_t = np.zeros((nfreqs, nfreqs, lmax))
    nl_array_pol = np.zeros((nfreqs, nfreqs, lmax))

    for c1, freq1 in enumerate(freqs):
        for c2, freq2 in enumerate(freqs):
            if c1>c2 : continue

            l, nl_t = np.loadtxt("%s/noise_t_%s_%sx%s_%s.dat"%(noise_data_dir, exp, freq1, exp, freq2), unpack=True)
            l, nl_pol = np.loadtxt("%s/noise_pol_%s_%sx%s_%s.dat"%(noise_data_dir, exp, freq1, exp, freq2), unpack=True)

            nl_array_t[c1, c2, lcut:lmax] = nl_t[lcut - 2:lmax - 2] * nsplits
            nl_array_pol[c1, c2, lcut:lmax] = nl_pol[lcut - 2:lmax - 2] * nsplits

    for i in range(lmax):
        nl_array_t[:,:,i] = symmetrize(nl_array_t[:,:,i])
        nl_array_pol[:,:,i] = symmetrize(nl_array_pol[:,:,i])

    return l, nl_array_t, nl_array_pol

def get_foreground_matrix(foreground_data_dir, fg_components, all_freqs, lmax):

    """This function uses the foreground power spectra generated by 'maps_to_params_prepare_sim_data'
    and generate a three dimensional array of foregroung power spectra [nfreqs,nfreqs,lmax].
    The different entries ([i,j,:]) of the array contains the fg power spectra for the different
    frequency channel pairs.
    for example fl_array_T[0,0,:] =>  fl_{f_{0},f_{0}),  fl_array_T[0,1,:] =>  fl_{f_{0},f_{1})
    this allows to have correlated fg between different frequency channels.
    (Not that for now, no fg are including in pol)

    Parameters
    ----------
    foreground_data_dir : string
      the folder containing the foreground power spectra
    fg_components: list of string
      the list of the foreground components we wish to include
    all_freqs: 1d array of string
      the frequencies we consider
    lmax: integer
      the maximum multipole for the noise power spectra
    """

    nfreqs = len(all_freqs)
    fl_array = np.zeros((2*nfreqs, 2*nfreqs, lmax))

    for c1, freq1 in enumerate(all_freqs):
        for c2, freq2 in enumerate(all_freqs):
            if c1 > c2 : continue

            fl_all_dict = {}
            for mode in ["tt", "te", "ee"]:

                fl_all = 0
                for fg_comp in fg_components[mode]:
                    l, fl = np.loadtxt("%s/%s_%s_%sx%s.dat"%(foreground_data_dir, mode, fg_comp, freq1, freq2), unpack=True)
                    fl_all += fl * 2 * np.pi / (l * (l + 1))

                fl_all_dict[mode] = fl_all
                if mode == "te" and freq1 != freq2:
                    fl_all = 0
                    for fg_comp in fg_components[mode]:
                        l, fl = np.loadtxt("%s/%s_%s_%sx%s.dat"%(foreground_data_dir, mode, fg_comp, freq2, freq1), unpack=True)
                        fl_all += fl * 2 * np.pi / (l * (l + 1))
                    fl_all_dict["et"] = fl_all

            fl_array[2*c1, 2*c2, 2:lmax] = fl_all_dict["tt"][:lmax-2]
            fl_array[2*c1+1, 2*c2+1, 2:lmax] = fl_all_dict["ee"][:lmax-2]
            fl_array[2*c1, 2*c2+1, 2:lmax] = fl_all_dict["te"][:lmax-2]
            if "et" in fl_all_dict:
                fl_array[2*c1+1, 2*c2, 2:lmax] = fl_all_dict["et"][:lmax-2]

    for i in range(lmax):
        fl_array[:,:,i] = symmetrize(fl_array[:,:,i])

    return l, fl_array

def multiply_alms(alms, bl, ncomp):

    """This routine mutliply the alms by a function bl

    Parameters
    ----------
    alms : 1d array
      the alms to be multiplied
    bl : 1d array
      the function to multiply the alms
    ncomp: interger
      the number of components
      ncomp = 3 if T,Q,U
      ncomp = 1 if T only
    """

    alms_mult = alms.copy()
    if ncomp == 1:
        alms_mult = hp.sphtfunc.almxfl(alms_mult, bl)
    else:
        for i in range(ncomp):
            alms_mult[i] = hp.sphtfunc.almxfl(alms_mult[i], bl)
    return alms_mult

def generate_noise_alms(nl_array_t, lmax, n_splits, ncomp, nl_array_pol=None):

    """This function generates the alms corresponding to the noise power spectra matrices
    nl_array_t, nl_array_pol. The function returns a dictionnary nlms["T", i].
    The entry of the dictionnary are for example nlms["T", i] where i is the index of the split.
    note that nlms["T", i] is a (nfreqs, size(alm)) array, it is the harmonic transform of
    the noise realisation for the different frequencies.

    Parameters
    ----------
    nl_array_t : 3d array [nfreq, nfreq, lmax]
      noise power spectra matrix for temperature data

    lmax : integer
      the maximum multipole for the noise power spectra
    n_splits: integer
      the number of data splits we want to simulate
    ncomp: interger
      the number of components
      ncomp = 3 if T,Q,U
      ncomp = 1 if T only
    nl_array_pol : 3d array [nfreq, nfreq, lmax]
      noise power spectra matrix for polarisation data
      (in use if ncomp==3)
    """

    nlms = {}
    if ncomp == 1:
        for k in range(n_splits):
            nlms[k] = curvedsky.rand_alm(nl_array_t,lmax=lmax)
    else:
        for k in range(n_splits):
            nlms["T", k] = curvedsky.rand_alm(nl_array_t, lmax=lmax)
            nlms["E", k] = curvedsky.rand_alm(nl_array_pol, lmax=lmax)
            nlms["B", k] = curvedsky.rand_alm(nl_array_pol, lmax=lmax)

    return nlms

def remove_mean(so_map, window, ncomp):

    """This function removes the mean value of the map after having applied the
    window function
    Parameters
    ----------
    so_map : so_map
      the map we want to subtract the mean from
    window : so_map or so_map tuple
      the window function, if ncomp=3 expect
      (win_t,win_pol)
    ncomp : integer
      the number of components
      ncomp = 3 if T,Q,U
      ncomp = 1 if T only

     """

    if ncomp == 1:
        so_map.data -= np.mean(so_map.data * window.data)
    else:
        so_map.data[0] -= np.mean(so_map.data[0] * window[0].data)
        so_map.data[1] -= np.mean(so_map.data[1] * window[1].data)
        so_map.data[2] -= np.mean(so_map.data[2] * window[1].data)

    return so_map

def get_effective_noise(nl_file_t, bl1, bl2, lmax,nl_file_pol=None, lcut=0):
    """This function returns the effective noise power spectrum
        which is defined as the ratio of the noise power spectrum and the
        beam window function.
        Parameters
        ----------
        nl_file_t: string
          name of the temperature noise power spectrum file, expect l,nl
        nl_file_p: string
          name of the polarisation noise power spectrum file, expect l,nl
        bl1: 1 d array
          beam harmonic transform
        bl2: 1 d array
          beam harmonic transform
        lmax : integer
          the maximum multipole for the noise power spectra
        lcut : integer
          the noise of SO being very red we will cut out the lowest modes to avoid
          leakage
    """

    bl1 = bl1[lcut:lmax]
    bl2 = bl2[lcut:lmax]
    if nl_file_pol == None:
        l, noise = np.loadtxt(nl_file_t, unpack=True)
        noise[lcut:lmax] /= (bl1 * bl2)
    else:
        noise = {}
        spectra = ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"]
        for spec in spectra:
            noise[spec] = np.zeros(lmax)

        l, noise_t = np.loadtxt(nl_file_t,unpack=True)
        l, noise_pol = np.loadtxt(nl_file_pol,unpack=True)

        noise["TT"][lcut:lmax] = noise_t[lcut:lmax] / (bl1 * bl2)
        noise["EE"][lcut:lmax] = noise_pol[lcut:lmax] / (bl1 * bl2)
        noise["BB"][lcut:lmax] = noise_pol[lcut:lmax] / (bl1 * bl2)

    return noise

def is_symmetric(mat, tol=1e-8):
    return np.all(np.abs(mat-mat.T) < tol)

def is_pos_def(mat):
    return np.all(np.linalg.eigvals(mat) > 0)

max_stars_repo_name: adrien-laposta/PSpipe

max_issues_repo_name: adrien-laposta/PSpipe

max_forks_repo_name: adrien-laposta/PSpipe


text:
#!/usr/bin/env python
# Copyright (c) Facebook, Inc. and its affiliates.
"""
Detectron2 training script with a plain training loop.

This script reads a given config file and runs the training or evaluation.
It is an entry point that is able to train standard models in detectron2.

In order to let one script support training of many models,
this script contains logic that are specific to these built-in models and therefore
may not be suitable for your own project.
For example, your research project perhaps only needs a single "evaluator".

Therefore, we recommend you to use detectron2 as a library and take
this file as an example of how to use the library.
You may want to write your own script with your datasets and other customizations.

Compared to "train_net.py", this script supports fewer default features.
It also includes fewer abstraction, therefore is easier to add custom logic.
"""

from dataclasses import dataclass
import logging
import os,sys,time
import numpy as np
import random
from pyexpat import features
from collections import OrderedDict
import torch, argparse
import torch.nn as nn
from torch.nn.parallel import DistributedDataParallel

import detectron2.utils.comm as comm
from detectron2.checkpoint import DetectionCheckpointer, PeriodicCheckpointer
from detectron2.config import get_cfg
from detectron2.data import (
    MetadataCatalog,
    build_detection_test_loader,
    build_detection_train_loader,
)
from detectron2.engine import default_argument_parser, default_setup, default_writers, launch
from detectron2.evaluation import (
    CityscapesInstanceEvaluator,
    CityscapesSemSegEvaluator,
    COCOEvaluator,
    COCOPanopticEvaluator,
    DatasetEvaluators,
    LVISEvaluator,
    PascalVOCDetectionEvaluator,
    SemSegEvaluator,
    inference_on_dataset,
    print_csv_format,
)
from detectron2.modeling import build_model
from detectron2.solver import build_lr_scheduler, build_optimizer
from detectron2.utils.events import EventStorage

from detectron2.data.dataset_mapper_pairs import DatasetMapperChangePairs
logger = logging.getLogger("detectron2")

def get_evaluator(cfg, dataset_name, output_folder=None):
    """
    Create evaluator(s) for a given dataset.
    This uses the special metadata "evaluator_type" associated with each builtin dataset.
    For your own dataset, you can simply create an evaluator manually in your
    script and do not have to worry about the hacky if-else logic here.
    """
    if output_folder is None:
        output_folder = os.path.join(cfg.OUTPUT_DIR, "inference")
    evaluator_list = []
    evaluator_type = MetadataCatalog.get(dataset_name).evaluator_type
    if evaluator_type in ["sem_seg", "coco_panoptic_seg"]:
        evaluator_list.append(
            SemSegEvaluator(
                dataset_name,
                distributed=True,
                output_dir=output_folder,
            )
        )
    if evaluator_type in ["coco", "coco_panoptic_seg"]:
        evaluator_list.append(COCOEvaluator(dataset_name, output_dir=output_folder))
    if evaluator_type == "coco_panoptic_seg":
        evaluator_list.append(COCOPanopticEvaluator(dataset_name, output_folder))
    if evaluator_type == "cityscapes_instance":
        return CityscapesInstanceEvaluator(dataset_name)
    if evaluator_type == "cityscapes_sem_seg":
        return CityscapesSemSegEvaluator(dataset_name)
    if evaluator_type == "pascal_voc":
        return PascalVOCDetectionEvaluator(dataset_name)
    if evaluator_type == "lvis":
        return LVISEvaluator(dataset_name, cfg, True, output_folder)
    if len(evaluator_list) == 0:
        raise NotImplementedError(
            "no Evaluator for the dataset {} with the type {}".format(dataset_name, evaluator_type)
        )
    if len(evaluator_list) == 1:
        return evaluator_list[0]
    return DatasetEvaluators(evaluator_list)

def do_test(cfg, model):
    results = OrderedDict()
    for dataset_name in cfg.DATASETS.TEST:
        data_loader = build_detection_test_loader(cfg, dataset_name)
        evaluator = get_evaluator(
            cfg, dataset_name, os.path.join(cfg.OUTPUT_DIR, "inference", dataset_name)
        )
        results_i = inference_on_dataset(model, data_loader, evaluator)
        results[dataset_name] = results_i
        if comm.is_main_process():
            logger.info("Evaluation results for {} in csv format:".format(dataset_name))
            print_csv_format(results_i)
    if len(results) == 1:
        results = list(results.values())[0]
    return results

def do_update_g(model_det, model_g, cfg_det, data):
        distributed = comm.get_world_size() > 1
        if "instances" in data[0]:
            gt_instances = [x["instances"].to(model_det.device) for x in data]
        else:
            gt_instances = None
        if distributed:
            images = model_det.module.preprocess_image(data)
            with torch.no_grad():
                features_orig = model_det.module.backbone(images.tensor)
                proposals, _ = model_det.module.proposal_generator(images, features_orig, None)
                proposals = model_det.module.roi_heads.label_and_sample_proposals(proposals, gt_instances)
                features = [features_orig[f] for f in model_det.module.roi_heads.box_in_features]
                box_features = model_det.module.roi_heads.box_pooler(features, [x.proposal_boxes for x in proposals])
        else:
            images = model_det.preprocess_image(data)
            with torch.no_grad():
                features_orig = model_det.backbone(images.tensor)
                proposals, _ = model_det.proposal_generator(images, features_orig, None)
                proposals = model_det.roi_heads.label_and_sample_proposals(proposals, gt_instances)
                features = [features_orig[f] for f in model_det.roi_heads.box_in_features]
                box_features = model_det.roi_heads.box_pooler(features, [x.proposal_boxes for x in proposals])

        masks_slide = []
        win_width=3
        win_overlap=1

        for i in range(3):
            for j in range(3):
                mask_ = torch.ones_like(box_features[0,0,:,:])
                x_start = i*win_width-i*win_overlap
                x_end = x_start+win_width
                y_start = i*win_width-i*win_overlap
                y_end = y_start+win_width

                mask_[x_start:x_end,y_start:y_end]=0
                masks_slide.append(mask_)

        #random perturbation on mask
        max_scores=[]
        for mask_ in masks_slide:

            with torch.no_grad():
                if cfg_det.NET_G.UPDATE_MODE == 'icassp' or cfg_det.NET_G.UPDATE_MODE == 'loss-cls':
                    model_det.train()
                    if distributed:
                        _, detector_loss_g = model_det.module.roi_heads(images,features_orig,proposals,gt_instances,model_g,mask_)
                    else:
                        _, detector_loss_g = model_det.roi_heads(images,features_orig,proposals,gt_instances,model_g,mask_)
                    if cfg_det.NET_G.UPDATE_MODE == 'loss-cls':
                        scores=detector_loss_g['loss_cls'].cpu()
                    else:
                        scores=sum(detector_loss_g.values()).cpu()
                    model_det.eval()
                else:
                    if cfg_det.NET_G.MASK_TYPE == 'icassp':
                        box_features_ = box_features * mask_
                    elif cfg_det.NET_G.MASK_TYPE == 'residual':
                        box_features_ = box_features + mask_
                    else:
                        raise NotImplementedError

                    if distributed:
                        box_features_ = model_det.module.roi_heads.box_head(box_features_)
                        predictions = model_det.module.roi_heads.box_predictor(box_features_)
                        pred_instances = model_det.module.roi_heads.box_predictor.inference(predictions,proposals)
                    else:
                        box_features_ = model_det.roi_heads.box_head(box_features_)
                        predictions = model_det.roi_heads.box_predictor(box_features_)
                        pred_instances = model_det.roi_heads.box_predictor.inference(predictions,proposals)

                    scores=[]

                    for it in pred_instances[0]:
                        if len(it.scores) == 0: #TODO: no such field
                            scores.append(torch.tensor(0,dtype=box_features_.dtype))
                        elif cfg_det.NET_G.UPDATE_MODE=='confuse':
                            idx = torch.argmin(torch.abs(it.scores-0.5))
                            scores.append(it.scores[idx].cpu())
                        else:
                            scores.append(it.scores.max().cpu()) #batchsize * 1

            max_scores.append(scores)
            #masks_pert.append(mask_)

        masks_pert = masks_slide

        #TODO: select target mask and update netG
        def select_mask(max_scores):
            '''
                max_scores: list(np.array(batchsize,1))
            '''
            if cfg_det.NET_G.UPDATE_MODE == 'icassp':
                return np.argmax(max_scores)
            elif cfg_det.NET_G.UPDATE_MODE == 'minmax':
                max_scores = [np.max(it) for it in max_scores]
                return np.argmin(max_scores)
            elif cfg_det.NET_G.UPDATE_MODE == 'minmin':
                min_score = [np.min(it) for it in max_scores]
                return np.argmin(min_score)
            elif cfg_det.NET_G.UPDATE_MODE == 'minmean':
                mean_score = [np.mean(it) for it in max_scores]
                return np.argmin(mean_score)
            elif cfg_det.NET_G.UPDATE_MODE == 'confuse':
                confuse = [np.abs(np.array(it)-0.5).mean() for it in max_scores]
                return np.argmin(confuse)
            else:
                raise NotImplementedError

        idx = select_mask(max_scores)
        mask_tar = masks_pert[idx]

        del masks_pert, max_scores

        # 'Updating net G ...'
        model_g.train(True)

        with torch.enable_grad():
            if distributed:
                mask_pre = model_g.module.netG(box_features)
            else:
                mask_pre = model_g.netG(box_features)
            if cfg_det.NET_G.G_MODE == 'spatial':
                mask_pre = torch.mean(mask_pre,dim=1,keepdim=True)
            if cfg_det.NET_G.G_MODE == 'channel':
                mask_pre = torch.mean(mask_pre,dim=(2,3),keepdim=True)

        loss = nn.L1Loss()(mask_pre,mask_tar.expand(mask_pre.shape))

        return {'loss_g':loss}

def do_train(cfg_g, model_g, cfg_det, model_det, resume=False):
    model_det.train()
    model_g.eval()

    optimizer = build_optimizer(cfg_det, model_det)
    scheduler = build_lr_scheduler(cfg_det, optimizer)

    optimizer_g = build_optimizer(cfg_g, model_g)
    scheduler_g = torch.optim.lr_scheduler.StepLR(
        optimizer_g,
        step_size= int(((cfg_det.SOLVER.MAX_ITER - cfg_det.NET_G.UPDATE_START)/cfg_det.NET_G.UPDATE_INTERVAL)*cfg_det.NET_G.G_STEP),
        gamma=1.0,
        )

    checkpointer_g = DetectionCheckpointer(
        model_g, cfg_g.OUTPUT_DIR,
    )
    # checkpointer_g.load(os.path.join(cfg_g.OUTPUT_DIR,'model_final.pth'))

    checkpointer = DetectionCheckpointer(
        model_det, cfg_det.OUTPUT_DIR, optimizer=optimizer, scheduler=scheduler
    )
    start_iter = (
        checkpointer.resume_or_load(cfg_det.MODEL.WEIGHTS, resume=resume).get("iteration", -1) + 1
    )
    max_iter = cfg_det.SOLVER.MAX_ITER

    periodic_checkpointer = PeriodicCheckpointer(
        checkpointer, cfg_det.SOLVER.CHECKPOINT_PERIOD, max_iter=max_iter
    )

    writers = default_writers(cfg_det.OUTPUT_DIR, max_iter) if comm.is_main_process() else []

    # compared to "train_net.py", we do not support accurate timing and
    # precise BN here, because they are not trivial to implement in a small training loop
    data_loader = build_detection_train_loader(cfg_det)
    logger.info("Starting training from iteration {}".format(start_iter))
    with EventStorage(start_iter) as storage:
        for data, iteration in zip(data_loader, range(start_iter, max_iter)):
            storage.iter = iteration

            #pretrain
            if iteration == cfg_det.NET_G.PRETRAIN_START:
                model_det.eval()
                model_g.train()
            if iteration  == cfg_det.NET_G.UPDATE_START+1:
                model_g.eval()
                model_det.train()

            if cfg_det.NET_G.PRETRAIN_START <= iteration <= cfg_det.NET_G.UPDATE_START:
                tic = time.time()
                model_det.eval()
                model_g.train()
                loss_dict_g = do_update_g(model_det,model_g,cfg_det,data)
                losses = sum(loss_dict_g.values())

                assert torch.isfinite(losses).all(), loss_dict_g

                loss_dict_reduced = {k: v.item() for k, v in comm.reduce_dict(loss_dict_g).items()}
                losses_reduced = sum(loss for loss in loss_dict_reduced.values())
                if comm.is_main_process():
                    storage.put_scalars(total_loss=losses_reduced, **loss_dict_reduced)

                optimizer_g.zero_grad()
                losses.backward()
                optimizer_g.step()
                scheduler_g.step()

                logger.info('Net-G Updated at iter {}   loss_g: {:.4f}   sec: {:.5f}'.format(iteration,  loss_dict_reduced['loss_g'], time.time()-tic))
                continue

            if (
                iteration > cfg_det.NET_G.UPDATE_START and
                (iteration + 1)%cfg_det.NET_G.UPDATE_INTERVAL==0 and
                cfg_det.NET_G.UPDATE_MODE != 'no'
            ):
                for i in range(cfg_det.NET_G.UPDATE_TIMES):
                #TODO: update netG
                    tic = time.time()
                    model_det.eval()
                    model_g.train()
                    loss_dict_g = do_update_g(model_det,model_g,cfg_det,data)
                    losses = sum(loss_dict_g.values())

                    assert torch.isfinite(losses).all(), loss_dict_g

                    loss_dict_reduced = {k: v.item() for k, v in comm.reduce_dict(loss_dict_g).items()}
                    losses_reduced = sum(loss for loss in loss_dict_reduced.values())
                    if comm.is_main_process():
                        storage.put_scalars(total_loss=losses_reduced, **loss_dict_reduced)

                    optimizer_g.zero_grad()
                    losses.backward()
                    optimizer_g.step()
                    scheduler_g.step()

                    model_g.eval()
                    model_det.train()
                    logger.info('Net-G Updated at iter {}  #{}  loss_g: {:.4f}   sec: {:.5f}'.format(iteration, i, loss_dict_reduced['loss_g'], time.time()-tic))
            if (
                iteration > cfg_det.NET_G.START_ITER and
                (iteration + 1)%cfg_det.NET_G.INTERVAL == 0
            ):
                loss_dict = model_det(data, model_g)
            else:
                loss_dict = model_det(data)

            losses = sum(loss_dict.values())
            assert torch.isfinite(losses).all(), loss_dict

            loss_dict_reduced = {k: v.item() for k, v in comm.reduce_dict(loss_dict).items()}
            losses_reduced = sum(loss for loss in loss_dict_reduced.values())
            if comm.is_main_process():
                storage.put_scalars(total_loss=losses_reduced, **loss_dict_reduced)

            optimizer.zero_grad()
            losses.backward()
            optimizer.step()
            storage.put_scalar("lr", optimizer.param_groups[0]["lr"], smoothing_hint=False)
            scheduler.step()

            if (
                cfg_det.TEST.EVAL_PERIOD > 0
                and (iteration + 1) % cfg_det.TEST.EVAL_PERIOD == 0
                and iteration != max_iter - 1
            ):
                do_test(cfg_det, model_det)
                # Compared to "train_net.py", the test results are not dumped to EventStorage
                comm.synchronize()

            if iteration - start_iter > 5 and (
                (iteration + 1) % 20 == 0 or iteration == max_iter - 1
            ):
                for writer in writers:
                    writer.write()
            periodic_checkpointer.step(iteration)

def setup(args):
    """
    Create configs and perform basic setups.
    """
    cfg_det = get_cfg()
    cfg_det.merge_from_file(args.config_det)
    cfg_det.merge_from_list(args.opts)
    cfg_det.freeze()

    cfg_g = get_cfg()
    cfg_g.merge_from_file(args.config_file)
    #cfg_g.merge_from_list(args.opts)
    cfg_g.freeze()

    default_setup(
        cfg_det, args
    )  # if you don't like any of the default setup, write your own setup code
    return cfg_det, cfg_g

def main(args):
    cfg_det,cfg_g = setup(args)

    model_det = build_model(cfg_det)
    logger.info("Detector Model:\n{}".format(model_det))

    model_g = build_model(cfg_g)
    logger.info("Net-G Model:\n{}".format(model_g))

    distributed = comm.get_world_size() > 1
    if distributed:
        model_det = DistributedDataParallel(
            model_det, device_ids=[comm.get_local_rank()], broadcast_buffers=False
        )
        model_g = DistributedDataParallel(
            model_g, device_ids=[comm.get_local_rank()], broadcast_buffers=False
        )

    do_train(cfg_g, model_g, cfg_det, model_det, resume=args.resume)
    return do_test(cfg_det,model_det)

def my_argument_parser(epilog=None):
    """
    Create a parser with some common arguments used by detectron2 users.

    Args:
        epilog (str): epilog passed to ArgumentParser describing the usage.

    Returns:
        argparse.ArgumentParser:
    """
    parser = argparse.ArgumentParser(
        epilog=epilog
        or f"""
Examples:

Run on single machine:
    $ {sys.argv[0]} --num-gpus 8 --config-file cfg.yaml

Change some config options:
    $ {sys.argv[0]} --config-file cfg.yaml MODEL.WEIGHTS /path/to/weight.pth SOLVER.BASE_LR 0.001

Run on multiple machines:
    (machine0)$ {sys.argv[0]} --machine-rank 0 --num-machines 2 --dist-url <URL> [--other-flags]
    (machine1)$ {sys.argv[0]} --machine-rank 1 --num-machines 2 --dist-url <URL> [--other-flags]
""",
        formatter_class=argparse.RawDescriptionHelpFormatter,
    )
    parser.add_argument("--config-file", default="", metavar="FILE", help="path to config file")
    parser.add_argument("--config-det", default="", metavar="FILE", help="path to config file of detector")
    parser.add_argument(
        "--resume",
        action="store_true",
        help="Whether to attempt to resume from the checkpoint directory. "
        "See documentation of `DefaultTrainer.resume_or_load()` for what it means.",
    )
    parser.add_argument("--eval-only", action="store_true", help="perform evaluation only")
    parser.add_argument("--num-gpus", type=int, default=1, help="number of gpus *per machine*")
    parser.add_argument("--num-machines", type=int, default=1, help="total number of machines")
    parser.add_argument(
        "--machine-rank", type=int, default=0, help="the rank of this machine (unique per machine)"
    )

    # PyTorch still may leave orphan processes in multi-gpu training.
    # Therefore we use a deterministic way to obtain port,
    # so that users are aware of orphan processes by seeing the port occupied.
    port = 2 ** 15 + 2 ** 14 + hash(os.getuid() if sys.platform != "win32" else 1) % 2 ** 14
    parser.add_argument(
        "--dist-url",
        default="tcp://127.0.0.1:{}".format(port),
        help="initialization URL for pytorch distributed backend. See "
        "https://pytorch.org/docs/stable/distributed.html for details.",
    )
    parser.add_argument(
        "opts",
        help="""
Modify config options at the end of the command. For Yacs configs, use
space-separated "PATH.KEY VALUE" pairs.
For python-based LazyConfig, use "path.key=value".
        """.strip(),
        default=None,
        nargs=argparse.REMAINDER,
    )
    return parser

if __name__ == "__main__":
    args = my_argument_parser().parse_args()
    print("Command Line Args:", args)
    launch(
        main,
        args.num_gpus,
        num_machines=args.num_machines,
        machine_rank=args.machine_rank,
        dist_url=args.dist_url,
        args=(args,),
    )


text:
# |------------------------------------------------------------------
# | # Flu Prediction  - Time Series Analysis TS2
# |------------------------------------------------------------------
# |
# | ## 1. Introduction
# |
# | This is a notebook to practice the routine procedures
# | commonly used in the time sequence analysis.
# | This notebook is based on the Kaggle [Time Series Analysis](https://www.kaggle.com/learn/time-series)
# | offered by Ryan Holbrook.

# | A time series consists of __trend__, __seasonality__, __cycles__,
# | and __peculiarity__. For each features, we have a procedure to deal with,
# |
# | - For __trend__ : Analytical fitting of the baselines (linear, polynomial, etc)
# | - For __seasonality__ : Fourier decomposition.
# | - For __cycle__ : Lags.
# | - For __peculiarity__ : Categorical features.
# |
# | In this notebook we will learn how to use lagged features to make predictions.

# | We also learn some terminology.
# | - __Forecast horizon__ : a part of the time-sequence where we do forecast.
# | - __Forecast origin__ : a point of the time-sequence where the training data ends.
# | - __Lead time__ : a part of the time-sequence after the forecast origin, but before
# | the forecast horizon starts.

# | When we have a forecast horizon longer than one unit time, the prediction requires,
# | - __Multioutput model__.

# | We have a couple of strategies how to create a multiple output.
# | - __Direct strategy__ :  Create one model for each day in the horizon,
# | and perform the prediction directly.
# | One needs so many models as the forecasting points in the forecast horizon.
# | - __Recursive strategy__ :  First, train a model to predict the first
# | day in the horizon. Only the given training data is used for the training.
# | Use that same model to predict the second day in the horizon, but
# | now we have one new input from the day before (=the forecast on the first
# | day in the horizon).
# | - __DirRec strategy__ : Combination of the  above two. Create a
# | model to forecast on the first day in the horizon. Use that new information
# | as a ground truth, and create the second model to forecast on the second day.
# | One needs so many models as the forecasting points in the forecast horizon.

# | ## 2. Task
# | From the historical record of the visits to the doctor's office in the past,
# | we will forecast the numbers of  such visits  in the future.

# | ## 3. Data
# | 1. The historical record of the visits to the doctor's office
# |     over a week, starting from 2009 and ending in 2016.
# |
# | 2. The data above comes with the Google search records related to
# |    a flu. The keywords and the number of searches are tabulated for
# |    each week.

# | ## 4. Notebook
# -------------------------------------------------------
# | Import packages.

from pathlib import Path
import os
import pandas as pd
import numpy as np

from scipy.signal import periodogram
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
from sklearn.multioutput import MultiOutputRegressor, RegressorChain

from statsmodels.tsa.deterministic import CalendarFourier, DeterministicProcess
from statsmodels.graphics.tsaplots import plot_pacf
from statsmodels.tsa.stattools import pacf
import statsmodels.api as sm
from xgboost import XGBRegressor

from IPython.display import display
import plotly.graph_objs as go
from plotly.subplots import make_subplots
import kaleido
from kaggle_tsa.ktsa import *

# -------------------------------------------------------
# | Set up directories.

CWD = Path('/Users/meg/git7/flu/')
DATA_DIR = Path('../input/ts-course-data/')
KAGGLE_DIR = Path('ryanholbrook/ts-course-data/')
IMAGE_DIR = Path('./images')
HTML_DIR = Path('./html')

os.chdir(CWD)
set_cwd(CWD)

# -------------------------------------------------------
# | If the data is not downloaded yet, do so now.

set_data_dir(KAGGLE_DIR, CWD)
show_whole_dataframe(True)

# -------------------------------------------------------
# | Read the data, first as it is.

flu = pd.read_csv(DATA_DIR/'flu-trends.csv')
print(flu.info())
display(flu.head(3))
# -------------------------------------------------------
# | First, let us only deal with `FluVisits`.
# | We would laso like to parse `Week`. We tried `parse_dates=["Week"]`,
# | and `infer_datetime_format=True`, but neither worked out.
# | There might be a simpler way, but this is all what I can
# | think of.

flu = pd.read_csv(DATA_DIR/'flu-trends.csv',
                  dtype={'FluVisits': 'float64'},
                  usecols=['Week', 'FluVisits'])

flu = split_week(flu, append=False)
is_index_continuous(flu, freq='W-MON')

# -------------------------------------------------------
# | Write a function to check the continuity of the index.
# | Raise flag if there are missing/skipped date in the sequence.
# | => Done.
# -------------------------------------------------------
# | Let us take a look at data again.

trace = go.Scatter(x=flu.index,
                   y=flu['FluVisits'])

data = [trace]
layout = go.Layout(height=512,
                   font=dict(size=16),
                   showlegend=False)

fig = go.Figure(data=data, layout=layout)
fig_wrap(fig, IMAGE_DIR/'fig1.png')

# -------------------------------------------------------
# | Create lag plots.

y = flu['FluVisits'].copy()
n_lag = 12
fig, corr = create_lag_plot(y, n_lag=n_lag, n_cols=3)
fig_wrap(fig, IMAGE_DIR/'fig2.png')

_ = [print(f'Lag {i:-2}: {y.autocorr(lag=i):5.3f}')
     for i in range(1, n_lag+1)]

# -------------------------------------------------------
# | Create PACF (Partial Autocorrelation Function) plot.

fig = create_pacf_plot(y, n_lag=n_lag)
fig_wrap(fig, IMAGE_DIR/'fig3.png')

# -------------------------------------------------------
# | Lag 1, 2, 3, and 4 have significant correlation
# | with the target.

# -------------------------------------------------------
# | Construct the time dummy.
# | __Note this trick!__ One can use dictionary `{}`
# | in `pd.concat` to specify the column names and the values
# | in the same time.
y = flu['FluVisits'].copy()
n_lag = 4
X = pd.concat({f'y_lag_{i}': y.shift(i) for i in range(1, n_lag+1)},
              axis=1).fillna(0.0)

# -------------------------------------------------------
# | Start the machine learning part.
# | Note __`shuffle=False`__, so that the test data comes
# | after the training data in the time sequence.

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, shuffle=False)

model = LinearRegression()
model.fit(X_train, y_train)
y_fit = model.predict(X_train)
y_pred = model.predict(X_test)

# -------------------------------------------------------
# | Let us see the results.

fig = show_training_results(X, y, X_train, y_fit, X_test, y_pred,
                            titles=('[Year]', '[Cases]',
                                    'Flu-visit predictions in 2015-2016 (Lags only)'))

fig_wrap(fig, IMAGE_DIR/'fig4.png')

# -------------------------------------------------------
# | Error evaluation.

train_rmse = mean_squared_error(y_train, y_fit, squared=False)
print(f'Train RMSE : \033[96m{train_rmse:6.2f}\033[0m')

test_rmse = mean_squared_error(y_test, y_pred, squared=False)
print(f'Test RMSE : \033[96m{test_rmse:6.2f}\033[0m')

# -------------------------------------------------------
# | The test error being smaller than the training error indicates
# | the model is still under fitting. There is a room to improve the
# | forecast.
# |
# | We will see if we can make the error smaller if we include
# | the Google search keywords. First take a look at the data.
# |

search = pd.read_csv(DATA_DIR/'flu-trends.csv')

# | Which words are most often Google-searched?

search.mean(axis=0).sort_values()

# -------------------------------------------------------
# | Let us take the keywords that contains a string 'flu'.

s_cols_week = search.columns[search.columns.str.contains('flu|Week')]
s_cols = search.columns[search.columns.str.contains('flu')]

search = pd.read_csv(DATA_DIR/'flu-trends.csv',
                     dtype={s: 'float64' for s in s_cols},
                     usecols=s_cols_week)

search = split_week(search, append=False)
is_index_continuous(search, freq='W-MON')

# -------------------------------------------------------
# | Create lagged time features.

y = flu['FluVisits'].copy()
n_lag = 4
X_lag = pd.concat({f'y_lag_{i}': y.shift(i) for i in range(1, n_lag+1)},
                  axis=1).fillna(0.0)
X_lag.isna().sum().sum()

# -------------------------------------------------------
# | Create lagged search-words features.

y_search = search[s_cols]
y_search.isna().sum()

X_search = pd.concat({f'y_lag_{i}': y_search.shift(i)
                      for i in range(1, n_lag+1)}, axis=1).fillna(0.0)

X_search.isna().sum().sum()

X = pd.concat([X_lag, X_search], axis=1)
X.isna().sum().sum()

# -------------------------------------------------------
# | Train, predict, and show.

y, X = y.align(X, join='inner', axis=0)

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, shuffle=False)

model = LinearRegression()
model.fit(X_train, y_train)
y_fit = model.predict(X_train)
y_pred = model.predict(X_test)

fig = show_training_results(X, y, X_train, y_fit, X_test, y_pred,
                            titles=('[Year]', '[Cases]',
                                    'Flu-visit predictions in 2015-2016 (Google Trends)'))

fig_wrap(fig, IMAGE_DIR/'fig5.png')
# -------------------------------------------------------
# | Error evaluation.

train_rmse = mean_squared_error(y_train, y_fit, squared=False)
print(f'Train RMSE : \033[96m{train_rmse:6.2f}\033[0m')

test_rmse = mean_squared_error(y_test, y_pred, squared=False)
print(f'Test RMSE : \033[96m{test_rmse:6.2f}\033[0m')

#  -------------------------------------------------------
# | Now the model is not under-fitting, but the forecast is not
# | as good as before. Need trials and errors to select
# | the search words.
# |
# | Overall, the reproduction of the `FluVisits` is excellent.
# |
# | * This is in a sense very much expected though,
# | from the high correlation between `FluVisits` and the data  on the one day before.
# | One can just tell some similar value with yesterday, and
# | the forecast will be quite good.

# -------------------------------------------------------
# | Now we will perform forecasting under more realistic conditions,
# | with a forecast horizon longer than one unit.
# |
# | We first try a simple __Multioutput model__.
# | To output the forecast on the multiple days, what we have to do is only one thing.
# |
# | 1. Make target (=`y`) concatenated for multiple steps in column-wise,
# | ```
# | y_step_1  y_step_2 y_step_3 ...
# | ```

# | Here we perform forecast on a 8-weeks horizon.
# |

y = flu['FluVisits'].copy()
n_step = 8
y = pd.concat({f'y_step_{i}': y.shift(-i) for i in range(0, n_step)},
              axis=1).fillna(0.0)
y.head(3)

# -------------------------------------------------------
# | Do not forget this.

y, X = y.align(X, join='inner', axis=0)
# -------------------------------------------------------
# | Do the training.

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, shuffle=False)

model = LinearRegression()
model.fit(X_train, y_train)
y_fit = model.predict(X_train)
y_pred = model.predict(X_test)

fig = show_training_results(X, y, X_train, y_fit, X_test, y_pred,
                            titles=('[Year]', '[Cases]',
                                    'Flu-visit predictions in 2015-2016 (Multi-Output)'))

fig_wrap(fig, IMAGE_DIR/'fig6.png')

# -------------------------------------------------------
# | Use `XGBoost` to perform Direct strategy.

model = MultiOutputRegressor(XGBRegressor())
model.fit(X_train, y_train)
y_fit = model.predict(X_train)
y_pred = model.predict(X_test)

fig = show_training_results(X, y, X_train, y_fit, X_test, y_pred,
                            titles=('[Year]', '[Cases]',
                                    'Flu-visit predictions in 2015-2016 (XGBRegressor)'))

fig_wrap(fig, IMAGE_DIR/'fig7.png')

# -------------------------------------------------------
# | Use `XGBoost` to perform DirRec strategy

model = RegressorChain(XGBRegressor())
model.fit(X_train, y_train)
y_fit = model.predict(X_train)
y_pred = model.predict(X_test)

fig = show_training_results(X, y, X_train, y_fit, X_test, y_pred,
                            titles=('[Year]', '[Cases]',
                                    'Flu-visit predictions in 2015-2016 (XGBRegressor)'))

fig_wrap(fig, IMAGE_DIR/'fig8.png')

# -------------------------------------------------------
# | END

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


text:
#!/usr/bin/env python3
import pandas as pd
import numpy as np
import gzip
import pickle
import sys
import argparse
import os

def choose_variant_annotation(csq_string, variant_annotation_rank_dict, gene_ix, conseq_ix):
    """Parse lowest ranking consequence and gene ID from CSQ string"""
    minrank = len(variant_annotation_rank_dict)
    gene_id = 'NA'
    conseq = 'NA'

    for tr in csq_string.split(','):
        annots = tr.split('|')
        for v in annots[conseq_ix].split('&'):
            if v in variant_annotation_rank_dict:
                r = variant_annotation_rank_dict[v]
                if r<minrank:
                    minrank = r
                    gene_id = annots[gene_ix]
                    conseq = v
    return gene_id, conseq

def get_vep_format(vep_vcf):
    """Format of the CSQ string"""
    fmt = None
    with gzip.open(vep_vcf) as f:
        for line in f:
            line = line.decode().strip()
            if line[0]!='#':
                break
            elif line[:6]=='##INFO' and 'CSQ' in line:
                fmt = line.split('Format: ')[1].replace('">','').split('|')
                break
    if fmt is None:
        raise ValueError('CSQ format not found in VCF header')
    return fmt

# from http://useast.ensembl.org/info/genome/variation/predicted_data.html
variant_annotation_rank_list = [
    'transcript_ablation',
    'splice_acceptor_variant',
    'splice_donor_variant',
    'stop_gained',
    'frameshift_variant',
    'stop_lost',
    'start_lost',
    'transcript_amplification',
    'inframe_insertion',
    'inframe_deletion',
    'missense_variant',
    'protein_altering_variant',
    'splice_region_variant',
    'incomplete_terminal_codon_variant',
    'stop_retained_variant',
    'synonymous_variant',
    'coding_sequence_variant',
    'mature_miRNA_variant',
    '5_prime_UTR_variant',
    '3_prime_UTR_variant',
    'non_coding_transcript_exon_variant',
    'intron_variant']

variant_annotation_rank_dict = {j:i for i,j in enumerate(variant_annotation_rank_list)}

if __name__=='__main__':

    parser = argparse.ArgumentParser(description='Extract single effect from VEP VCF')
    parser.add_argument('vcf', help='VCF')
    parser.add_argument('vep_vcf', help='Variant Effect Predictor VCF')
    parser.add_argument('-o', '--output_dir', default='.')
    args = parser.parse_args()

    #------------------------------------------------
    # Pre-processing: extract variant set from VCF
    #------------------------------------------------
    prefix = os.path.basename(args.vcf).split('.vcf')[0]
    variant_id_file = os.path.join(args.output_dir, prefix+'.variant_ids.txt.gz')
    if not os.path.exists(variant_id_file):
        print('Extracting variant IDs to {}'.format(variant_id_file))
        subprocess.check_call('zcat '+args.vcf+' | grep -v "#" | cut -f3 | gzip -c -1 > '+variant_id_file)

    print('Loading variant IDs')
    with gzip.open(variant_id_file) as f:
        v = f.read()
    variant_set = set(v.decode().strip().split('\n'))

    # get position of Gene and Consequence fields in CSQ string (changes with VEP versions)
    fmt = get_vep_format(args.vep_vcf)
    gene_ix = np.where(np.array(fmt)=='Gene')[0][0]
    conseq_ix = np.where(np.array(fmt)=='Consequence')[0][0]

    # parse VEP VCF and write variant_id, gene, consequence to file
    # use pandas parser for speed/process by chunks to reduce memory usage
    print('Parsing VEP VCF')
    vep_file = os.path.join(args.output_dir, prefix+'.vep.txt.gz')
    with gzip.open(vep_file, 'wt') as f:
        f.write('variant_id\tgene_id\tvep\n')
        for k,chunk in enumerate(pd.read_csv(args.vep_vcf, sep='\t', comment='#', header=None,
            names=['chr', 'pos', 'id', 'ref', 'alt', 'qual', 'filter', 'info'], dtype=bytes, usecols=['id', 'info'],
            iterator=True, chunksize=100000)):

            print('\rProcessing chunk {0:d}'.format(k+1), end='')
            for i,c in zip(chunk['id'], chunk['info']):
                if len(c.split('|')) > 10:
                    cs = [i for i in c.split(';') if i.startswith('CSQ')][0]
                    gene_id, conseq = choose_variant_annotation(cs, variant_annotation_rank_dict, gene_ix, conseq_ix)
                    if gene_id != 'NA':
                        for j in i.split(';'):
                            if j in variant_set:
                                f.write(j+'\t'+gene_id+'\t'+conseq+'\n')
            print()

    print('Converting to dictionary')
    vep_df = pd.read_csv(vep_file, sep='\t', index_col=0)
    vep_df['combined'] = [(i,j) for i,j in zip(vep_df['gene_id'], vep_df['vep'])]
    vep_dict = vep_df['combined'].to_dict()

    print('Saving as pickle')
    with open(os.path.join(args.output_dir, prefix+'.vep_dict.pickle'), 'wb') as f:
        pickle.dump(vep_dict, f, pickle.HIGHEST_PROTOCOL)

max_stars_repo_name: richardslab/gtex-pipeline

max_issues_repo_name: richardslab/gtex-pipeline

max_forks_repo_name: richardslab/gtex-pipeline


text:
import numpy as np
import cv2

def _rect_to_tuple(rect):
    left = rect.left()
    right = rect.right()
    top = rect.top()
    bottom = rect.bottom()
    return left, top, right, bottom

def _extract_eye(shape, eye_indices):
    points = map(lambda i: shape.part(i), eye_indices)
    return list(points)

def _extract_eye_center(shape, eye_indices):
    points = _extract_eye(shape, eye_indices)
    xs = map(lambda p: p.x, points)
    ys = map(lambda p: p.y, points)
    return sum(xs) // 6, sum(ys) // 6

def extract_left_eye_center(shape):
    LEFT_EYE_INDICES = [36, 37, 38, 39, 40, 41]
    return _extract_eye_center(shape, LEFT_EYE_INDICES)

def extract_right_eye_center(shape):
    RIGHT_EYE_INDICES = [42, 43, 44, 45, 46, 47]
    return _extract_eye_center(shape, RIGHT_EYE_INDICES)

def _angle_between_2_points(p1, p2):
    x1, y1 = p1
    x2, y2 = p2
    if x1 != x2:
        tan = (y2 - y1) / (x2 - x1)
        return np.degrees(np.arctan(tan))
    elif y2 > y1:
        return np.degrees(np.pi / 2)

    else:
        return np.degrees(-np.pi / 2)

def get_rotation_matrix(p1, p2):
    angle = _angle_between_2_points(p1, p2)
    x1, y1 = p1
    x2, y2 = p2
    xc = (x1 + x2) // 2
    yc = (y1 + y2) // 2
    M = cv2.getRotationMatrix2D((xc, yc), angle, 1)
    return M

def crop_image(image, det):
    left, top, right, bottom = _rect_to_tuple(det)
    return image[top:bottom, left:right]

max_stars_repo_name: ZohraRezgui/inaFaceGender

max_issues_repo_name: ZohraRezgui/inaFaceGender

max_forks_repo_name: ZohraRezgui/inaFaceGender


text:
import numpy as np
import torch

from ...image.spatial import rfft2d_freqs

def fft_scale(h: int, w: int, decay_power: float):
    freqs = rfft2d_freqs(h, w)
    scale = np.sqrt(h * w) / np.maximum(freqs, 1. / max(h, w)) ** decay_power
    scale = torch.from_numpy(scale).float()
    scale = scale.view(1, 1, *scale.shape, 1)
    return scale

max_stars_repo_name: HealthML/lucid-torch

max_issues_repo_name: HealthML/lucid-torch

max_forks_repo_name: HealthML/lucid-torch


text:
"""
Here are some classifiers.
"""
import numpy
import theano
import theano.tensor as T
from sklearn.metrics import confusion_matrix
import matplotlib
# matplotlib.use('Agg')
import matplotlib.pyplot as plt
from model import Model
from layer import SigmoidLayer

import pdb

class Classifier(Model):
    def __init__(self, n_in, n_out, varin=None, vartruth=None):
        super(Classifier, self).__init__(n_in, n_out, varin=varin)
        if not vartruth:
            vartruth = T.lvector('truth')
        assert isinstance(vartruth, T.TensorVariable)
        self.vartruth = vartruth

    def output(self):
        """
        output() is defined as the activity of the highest layer, not the
        prediction results. Prediction results are generated in the predict()
        method. The input and output should always be theano variables.
        """
        raise NotImplementedError("Must be implemented by subclass.")

    def predict(self):
        raise NotImplementedError("Must be implemented by subclass.")

    # Following are for analysis ----------------------------------------------

    def analyze_performance(self, data, truth, verbose=True):
        """
        TODO: WRITEME

        data : numpy.ndarray
        truth : numpy.ndarray
        verbose : bool
        """
        assert data.shape[0] == truth.shape[0], "Data and truth shape mismatch."
        if not hasattr(self, '_get_predict'):
            self._get_predict = theano.function([self.varin], self.predict())

        cm = confusion_matrix(truth, self.get_predict(data))
        pr_a = cm.trace()*1.0 / test_truth.size
        pr_e = ((cm.sum(axis=0)*1.0 / test_truth.size) * \
            (cm.sum(axis=1)*1.0 / test_truth.size)).sum()
        k = (pr_a - pr_e) / (1 - pr_e)
        print "OA: %f, kappa index of agreement: %f" % (pr_a, k)
        if verbose: # Show confusion matrix
            if not hasattr(self, '_fig_confusion'):
                self._fig_confusion = plt.gcf()
                self.ax = self._fig_confusion.add_subplot(111)
                self.confmtx = self.ax.matshow(cm)
                plt.title("Confusion Matrix")
                self._fig_confusion.canvas.show()
            else:
                self.confmtx.set_data(cm)
                self._fig_confusion.canvas.draw()
            print "confusion matrix:"
            print cm

class LogisticRegression(Classifier):
    def __init__(self, n_in, n_out, varin=None, vartruth=None,
                 init_w=None, init_b=None, npy_rng=None):
        if not npy_rng:
            npy_rng = numpy.random.RandomState(123)
        assert isinstance(npy_rng, numpy.random.RandomState)
        self.npy_rng = npy_rng

        self.n_out = n_out
        if not init_b:
            init_b = theano.shared(
                value=numpy.zeros(n_out, dtype=theano.config.floatX),
                name=self.__class__.__name__ + '_b', borrow=True)
        else:
            assert init_b.get_value().shape == (n_out,)
        self.b = init_b

        self.vartruth = vartruth
        self.init_w = init_w
        self.n_in, self.varin = n_in, varin
        if self.n_in != None:
            self._init_complete()

    def _init_complete(self):
        assert self.n_in != None, "Need to have n_in attribute to execute."
        super(LogisticRegression, self).__init__(
            self.n_in, self.n_out, varin=self.varin, vartruth=self.vartruth
        )

        if not self.init_w:
            w = numpy.asarray(self.npy_rng.uniform(
                low = -4 * numpy.sqrt(6. / (self.n_in + self.n_out)),
                high = 4 * numpy.sqrt(6. / (self.n_in + self.n_out)),
                size=(self.n_in, self.n_out)), dtype=theano.config.floatX)
            self.init_w = theano.shared(value=w,
                name=self.__class__.__name__ + '_w', borrow=True)
        self.w = self.init_w

        self.params = [self.w, self.b]

    def fanin(self):
        self.varfanin = T.dot(self.varin, self.w) + self.b
        return self.varfanin

    def output(self, fanin=None):
        """The output of a logistic regressor is p_y_given_x."""
        if fanin == None:   fanin = self.fanin()
        return T.nnet.softmax(fanin)

    def activ_prime(self):
        """for the special relationship between softmax and sigmoid, we can
        define this method in this classifier."""
        return SigmoidLayer(
            self.n_in, self.n_out, varin=self.varin,
            init_w=self.w, init_b=self.b
        ).activ_prime()

    def cost(self):
        return -T.mean(
            T.log(self.output())[
                T.arange(self.vartruth.shape[0]), self.vartruth]
        )

    def weightdecay(self, weightdecay=1e-3):
        return weightdecay * (self.w**2).sum()

    def predict(self):
        return T.argmax(self.output(), axis=1)

class BinaryLogisticRegression(LogisticRegression):
    def __init__(self, n_in, n_out, mode='stochastic', varin=None,
                 vartruth=None, init_w=None, init_b=None, npy_rng=None):
        super(BinaryLogisticRegression, self).__init__(
            n_in=n_in, n_out=n_out, varin=varin, vartruth=vartruth,
            init_w=init_w, init_b=init_b, npy_rng=npy_rng
        )
        self.srng = theano.sandbox.rng_mrg.MRG_RandomStreams(
            self.npy_rng.randint(999999))
        self.mode = mode

    def fanin(self):
        self.varfanin = T.dot(self.varin, self.binarized_weight()) + self.b
        return self.varfanin

    def binarized_weight(self):
        self.w0 = (
            4 * numpy.sqrt(6. / (self.n_in + self.n_out)) / 2
        ).astype(theano.config.floatX)
        if self.mode == 'deterministic':
            self.wb = T.switch(T.ge(self.w, 0), self.w0, -self.w0)

        elif self.mode == 'stochastic':
            # probability=hard_sigmoid(w/w0)
            p = T.clip(((self.w / self.w0) + 1) / 2, 0, 1)
            p_mask = T.cast(self.srng.binomial(n=1, p=p, size=T.shape(self.w)),
                            theano.config.floatX)

            # [0,1] -> -W0 or W0
            self.wb = T.switch(p_mask, self.w0, -self.w0)

        else:
            raise ValueError("Parameter 'self.mode' has to be either "
                             "'deterministic' or 'stochastic'")

        return self.wb

    def quantized_bprop(self, cost):
        index_low = T.switch(self.varin > 0.,
            T.floor(T.log2(self.varin)), T.floor(T.log2(-self.varin))
        )
        index_low = T.clip(index_low, -4, 3)
        sign = T.switch(self.varin > 0., 1., -1.)
        # the upper 2**(integer power) though not used explicitly.
        # index_up = index_low + 1
        # percentage of upper index.
        p_up = sign * self.varin / 2**(index_low) - 1
        index_random = index_low + self.srng.binomial(
            n=1, p=p_up, size=T.shape(self.varin), dtype=theano.config.floatX)
        quantized_rep = sign * 2**index_random

        error = T.grad(cost=cost, wrt=self.varfanin)

        self.dEdW = T.dot(quantized_rep.T, error)

class LinearSVM(LogisticRegression):
    def output(self, fanin=None):
        """The output of a logistic regressor is p_y_given_x."""
        if fanin == None:   fanin = self.fanin()
        return fanin

    def activ_prime(self):
        raise NotImplementedError("Must be implemented by subclass.")

    def cost(self):
        """Squared hinge loss"""
        return T.mean(
            T.sqr(T.maximum(
                0.,
                1. - (2 * T.extra_ops.to_one_hot(
                    self.vartruth, self.n_out) - 1) * self.output()
            ))
        )

class BinaryLinearSVM(BinaryLogisticRegression):
    def output(self, fanin=None):
        """The output of a logistic regressor is p_y_given_x."""
        if fanin == None:   fanin = self.fanin()
        return fanin

    def activ_prime(self):
        raise NotImplementedError("Must be implemented by subclass.")

    def cost(self):
        """Squared hinge loss"""
        return T.mean(
            T.sqr(T.maximum(
                0.,
                1. - (2 * T.extra_ops.to_one_hot(
                    self.vartruth, self.n_out) - 1) * self.output()
            ))
        )

class Perceptron(Classifier):
    def __init__(self):
        #
        # TODO:
        #
        pass

class LinearRegression(Classifier):
    def __init__(self, n_in, n_out, varin=None, vartruth=None,
                 init_w=None, init_b=None, npy_rng=None):
        super(LinearRegression, self).__init__(
            n_in, n_out, varin=varin, vartruth=vartruth
        )

        if not npy_rng:
            npy_rng = numpy.random.RandomState(123)
        assert isinstance(npy_rng, numpy.random.RandomState)
        self.npy_rng = npy_rng

        if not init_w:
            w = numpy.asarray(npy_rng.uniform(
                low = -4 * numpy.sqrt(6. / (n_in + n_out)),
                high = 4 * numpy.sqrt(6. / (n_in + n_out)),
                size=(n_in, n_out)), dtype=theano.config.floatX)
            init_w = theano.shared(value=w, name='w_sigmoid', borrow=True)
        # else:
        #     TODO. The following assetion is complaining about an attribute
        #     error while passing w.T to init_w. Considering using a more
        #     robust way of assertion in the future.
        #     assert init_w.get_value().shape == (n_in, n_out)
        self.w = init_w

        if not init_b:
            init_b = theano.shared(value=numpy.zeros(
                                       n_out,
                                       dtype=theano.config.floatX),
                                   name='b_sigmoid', borrow=True)
        else:
            assert init_b.get_value().shape == (n_out,)
        self.b = init_b

        self.params = [self.w, self.b]

    def fanin(self):
        return T.dot(self.varin, self.w) + self.b

    def output(self):
        return self.fanin()

    def activ_prime(self):
        """for the special relationship between LinearRegression and LinearLayer, we can
        define this method in this classifier."""
        return LinearLayer(
            self.n_in, self.n_out, varin=self.varin,
            init_w=self.w, init_b=self.b
        ).activ_prime()

    def cost(self):
        # TODO not known
        raise NotImplementedError("ERROR")
        return str("ERROR")

    def weightdecay(self, weightdecay=1e-3):
        return weightdecay * (self.w**2).sum()

    def predict(self):
        return T.argmax(self.output(), axis=1)

max_stars_repo_name: hantek/NeuroBricks

max_issues_repo_name: hantek/NeuroBricks

max_forks_repo_name: hantek/NeuroBricks

reason: import numpy,import theano


text:
from unittest import TestCase

from pyanp.rating import *
from numpy.testing import assert_array_equal, assert_allclose
class TestRating(TestCase):

    def test_crud(self):
        rating = Rating()
        alts = ["a1", "a2", "a3"]
        rating.add_alt(alts)
        self.assertEqual(3, rating.nalts())
        unames = ("Bill", "Lee", "Paul", "Frank")
        users=rating.add_user(unames)
        self.assertEqual(4, rating.nusers())
        rating.vote_column("a1", ["H", "M", "L", "M"])
        rating.vote_column("a2", ["good", "good", "bad", "good"])
        ps = rating.priority()
        assert_allclose([0.5, 0.625, 0.0], ps)

    def test_manual_setting(self):
        rating = Rating()
        alts = ["a1", "a2", "a3"]
        rating.add_alt(alts)
        self.assertEqual(3, rating.nalts())
        unames = ("Bill", "Lee", "Paul", "Frank")
        users=rating.add_user(unames)
        rating.set_word_eval({
            "hi":(1, "high", "h"),
            "med":(0.3, "medium", "m", "md"),
            "low":(0.1, "lo", "l", "bad"),
            "stupid":(0.005, "stupd", "stu", "st", "s")
        })
        rating.vote_column("a1", ["H", "M", "L", "M"])
        rating.vote_column("a2", ["good", "good", "bad", "good"])
        rating.vote_column("a3", ["stupid", "st", "bad", "good"])
        ps = rating.priority()
        print(ps)
        print(rating.vote_values(alt_name="a1"))

max_stars_repo_name: cznewt/pyanp

max_issues_repo_name: ArdalanMarandi/pyanp

max_forks_repo_name: ArdalanMarandi/pyanp

reason: from numpy


text:
import copy
from pylab import *
from numpy import *

if False:
    out = open("Irish.txt", "w")
    times = open("Irish_orig.txt").readlines()
    now = 0
    for l in times:
        if '(' in l:
            l = l.split('(')[0]
        if ('-') in l:
            l = l.split('-')[0]
        l = l.strip()
        print('%02d:%02d' % (now / 60, now % 60), l.lower(), file=out)
        now += 5
    out.close()
    print('wrote', out.name)
times = open("Irish_v5.txt").readlines()
sentences = [l[5:].split() for l in times]
times = [l[:5] for l in times]

words = set()
for l in sentences:
    for w in l:
        words.add(w)
for w in list(words):
    assert len(w) < 16

order = []
for l in sentences:
    words = l
    i = 0
    for word in words: ## see if we need to add word to ordred words
        if word not in order[i:]:
            order.insert(i + 1, word)
            i = i + 1
        else:
            i = order.index(word, i)
        if False:
            print(word)
            print(' '.join(words))
            print(', '.join(order))
            print()

## define the partial ordering on the numbers
## a < b if a occurs before b in a sentence
def getAdj(order):
    n = len(order)
    adj = zeros((n, n), bool)
    for s in numeric:
        for ii, i in enumerate(s[:-1]):
            adj[s[ii], s[ii + 1]] = True
    return adj
words = set(order)

def precedence(adj, max_dist):
    ### a non zero in the i, j, entry is True if word i must preceed word j
    n = len(adj)
    adj = array(adj, copy=True).astype(int)
    A_n = eye(n).astype(int)
    out = eye(n).astype(int)
    for i in range(max_dist):
        A_n = dot(adj, A_n)
        out += A_n
    return (out  - eye(n)).astype(bool).astype(int)
#
# A = array([[0, 1, 0, 1],[0, 0, 1, 1], [0, 0, 0, 1], [0, 0, 0, 0]])
# print A
# print precedence(A, 10)

def maximal_dist(adj, diameter):
    out = array(adj, copy=True).astype(int)
    n = len(adj)
    A_n = eye(n)
    for i in range(1, diameter):
        A_n = dot(adj, A_n).astype(bool)
        not_A_n = True - A_n
        out = out * not_A_n +  A_n * i
    return out

def mat_pow(a, k):
    out = eye(a.shape[0])
    for i in range(k):
        out = dot(out, a)
    return out

def minimal_dist(adj, diameter):
    n = len(adj)
    shape = (diameter, n, n)
    A_ns = zeros(shape, int)
    for i in range(diameter):
        A_ns[i] = mat_pow(adj, i)
    for i in range(diameter):
        A_ns[i] = where(A_ns[i], i, diameter)
    out = amin(A_ns, axis=0)
    return where(equal(out, diameter), 0, out)

V_MAT = zeros((100, 100), bool)
for i in range(100): ## 100:: can be any number larger than the number of words
    V_MAT[i,:i] = 1
def is_poset(adj):
    n = len(adj)
    A = precedence(adj, 10)
    return sum(ravel(A * V_MAT[:n, :n])) == 0

def delete_word(adj, idx):
    n = len(adj)
    out = zeros((n - 1, n - 1), adj.dtype)
    out[:idx, :idx] = adj[:idx, :idx]
    out[:idx,idx:] = adj[:idx,idx+1:]
    out[idx:,:idx] = adj[idx+1:,:idx]
    out[idx:,idx:] = adj[idx+1:,idx+1:]
    return out

def validate(order, sentences):
    for sentence in sentences:
        i = 0
        for w in sentence:
            if w in order[i:]:
                pass
            else:
                return False
            i = order.index(w, i) + 1
    return True
order_cp = order[:]
i = 0
print('pruning redundant words')
print(len(order_cp))

if False:
    rm_count = 0
    for i in range(len(order)):
        w = order[i]
        order_cp = order[:i] + order[i+1:]
        if validate(order_cp, sentences):
            print(w, 'redundent')
order_cp = order[:]
i = 0
while i < len(order_cp):
    w = order_cp[i]
    order_cp = order_cp[:i] + order_cp[i+1:]

    ## see if later word can replace earlier word and still be valid
    if validate(order_cp, sentences):
        # print 'rm', w
        pass
    else:
        order_cp.insert(i, w)
        i += 1
print('before pruning: n_word=%d, n_letter=%d' % (len(order), sum([len(w) for w in order])))
print(' after pruning: n_word=%d, n_letter=%d' % (len(order_cp), sum([len(w) for w in order_cp])))

porder = order_cp
# porder = open("Irish_v5_porder.txt").read().split()
print(porder)
assert validate(porder, sentences)

### now we have a minimal set of words that we can reference by index.  Turn the sentences in to a sequence of induces.
numeric = []

for l in sentences:
    # print ' '.join(l)
    i = 0
    ns = []
    for w in l:
        try:
            i = porder.index(w, i)
        except:
            print(' '.join(l))
            print(w)
            print(porder.index(w), i)
            print(','.join(porder))
            raise
        ns.append(i)
    numeric.append(ns)

adj = getAdj(porder)
min_dist = minimal_dist(adj, 30)
print('Pruned list:', ' '.join(porder))

A = precedence(adj, 20)
assert is_poset(adj)

figure(3)
max_dist = maximal_dist(adj, 100)
min_dist = minimal_dist(adj, 100)
dia = max(max_dist.ravel())
print('max_dist from origin', dia)
subplot(311); pcolormesh(max_dist)
colorbar()
title('Maximal Dist')
subplot(312); pcolormesh(min_dist)
colorbar()
subplot(313); pcolormesh(A)
colorbar()

## choose first word to be the word that preceeds the most other words
first_index = argmax(sum(A, axis=1))
head = porder[first_index] ## first word

# the find maximal distance from first word to each other word

dist_from_head = max_dist[first_index]
tree = []
for i in range(13):
    tree.append([j for j in range(len(porder)) if dist_from_head[j] == i])

n_poss = 1
n  = len(A)
print(tree)
for k in range(len(tree)):
    print(k, len(porder[k]), ' '.join([porder[i] for i in tree[k]]))
    if k + 1 < len(tree):
        for j in tree[k + 1]:
            for i in tree[k]:
                print(int(A[i, j]), end=' ')
            print()
        print()
        print()

def overlap(s1, s2):
    n = min([len(s1), len(s2)])
    out = 0
    for i in range(1, n + 1):
        if s1[-i:] == s2[:i]:
            out = i
    return out
def overlap__test__():
    assert overlap("justin", "in case") == 2
    assert overlap("justin", "boston") == 0
    assert overlap("justin", "justin") == 6
overlap__test__()

def addword(letters, word):
    '''
    add word to letter matrix.  put a space if you don't want to overlap
    '''
    x = overlap(letters, word)
    if len(word) - x <= len(letters) % 16:
        letters = letters + word[x:]
    else:
        letters += ' '  * (len(letters) % 16)
        letters += word
    return letters

def doit(sofar, tree):
    if len(tree) == 0:
        out = sofar
    else:
        tree = [l[:] for l in tree]
        word = porder[tree[0].pop(0)]
        if len(tree[0]) == 0:
            tree.pop(0)
        sofar = addword(sofar, word)
        return doit(sofar, tree)
    return out

figure(1)
A_n = eye(len(porder))
for i in range(dia):
    subplot(5, 3, i + 1)
    A_n = dot(A_n, adj)
    pcolormesh(A_n)

figure(2)
ax = subplot(211)
title("Irish_v3 precedence")
pcolormesh(adj)
subplot(212, sharex=ax, sharey=ax)
pcolormesh(A)
# show()
for w in porder:
    print(w)
show()

max_stars_repo_name: wyolum/ClockIOT

max_issues_repo_name: jasoncoon/ClockIOT

max_forks_repo_name: jasoncoon/ClockIOT


text:
from ._utils_find_tips import *
from ._find_tips import *
from numba import njit
import numpy as np
from ._utils_find_tips import reduce_tips
from .interpolate import *
# from ..operari import count_tips

@njit
def __pbc_1d(x_in,width):
	if np.greater_equal(x_in , width):
		return int(x_in-width)
	elif np.less(x_in, 0):
		return int(x_in+width)
	return int(x_in)

@njit
def bilinear_interpolate(x,y,width,height,txt):
	r0 = np.floor(y)
	r1 = np.ceil(y)
	frac_r = y-r0
	if r0==r1:
		r1+=1
	c0 = np.floor(x)
	c1 = np.ceil(x)
	frac_c = x-c0
	if c0==c1:
		c1+=1
	c0 = __pbc_1d(c0,width)
	c1 = __pbc_1d(c1,width)
	r0 = __pbc_1d(r0,height)
	r1 = __pbc_1d(r1,height)
	a11 = txt[r0,c0]
	a21 = txt[r1,c0]-a11
	a12 = txt[r0,c1]-a11
	a22 = txt[r1,c1]-a11-a12-a21
	states_bilinear = a11 + a21*frac_c + a12*frac_r + a22*frac_c*frac_r
	v,f,s = states_bilinear
	return list(states_bilinear)

@njit
def interpolate_states(x_values,y_values,width,height,txt):
	'''Measure color of txt at zip(x_values,y_values).  This was used to measure EP states.
	Example Usage:
	v_lst, f_lst, s_lst = interpolate_states(x_values,y_values,width,height,txt)
	'''
	v_lst = []; f_lst = []; s_lst = [];
	for x,y in zip(x_values,y_values):
		v,f,s = bilinear_interpolate(x,y,width,height,txt)
		v_lst.append(v)
		f_lst.append(f)
		s_lst.append(s)
	return v_lst, f_lst, s_lst

def measure_system(contours1, contours2, width, height, texture, jump_threshold = 2, size_threshold = 6, pad=1, decimals=11):
	'''returns s1_list, s2_list, x_lst, y_lst, v_lst, f_lst, s_lst
	Nota Bene: The vast majority of the runtime of measure_system results from find_tips due to the list of lists datastructure necessary to encode topological information.
	# Actual acceleration of find_tips using C or LLVM is not straightforward due to memory allocation issues...
	# Speedup is possible, but topological information would be destroyed.'''
	contour1_lst_lst, contour2_lst_lst = preprocess_contours(contours1, contours2, width, height, jump_threshold = jump_threshold, size_threshold = size_threshold)
	s1_list, s2_list, x_lst, y_lst = find_tips(contour1_lst_lst, contour2_lst_lst)
	s1_list, s2_list, x_lst, y_lst = reduce_tips(s1_list, s2_list, x_lst, y_lst, width, height, pad=pad, decimals=decimals)
	v_lst, f_lst, s_lst = interpolate_states(x_lst,y_lst,width,height,texture)
	return s1_list, s2_list, x_lst, y_lst, v_lst, f_lst, s_lst
	# dict_out = {
	# 		't': float(t),
	# 		'n': int(n_tips),
	# 		'x': tuple(x_values),
	# 		'y': tuple(y_values),
	# 		'n1': tuple(s1_list),
	# 		'n2': tuple(s2_list),
	# 		'v':v_lst,
	# 		'f':f_lst,
	# 		's':s_lst,
	# 	}
	# sys.getsizeof(dict_out) one time step fills about 360 [bytes or bits?]
	# return dict_out

max_stars_repo_name: timtyree/bgmc

max_issues_repo_name: timtyree/bgmc

max_forks_repo_name: timtyree/bgmc

reason: import numpy,from numba


text:
import numpy as np
import pandas
import random
import math
import matplotlib.pyplot as plt

class perceptron:
	def __init__(self,input_size,output_size,learning_rate):
		self.input_size = input_size
		self.output_size = output_size
		self.learning_rate = learning_rate
		self.weights = self.init_weights()
		self.bias = np.random.randn()
		print(self.weights)
		print(self.bias)

	def init_weights(self):
		w = []
		for x in range(self.input_size):
			w.append(np.random.randn())
		return w

	def train(self,data_x,data_y):
		for i in range(data_x.shape[0]):
			x = data_x[i].tolist()
			y = data_x[i].tolist()
			x1 = 0.0
			for j in range(self.input_size):
				x1 += self.weights[j]*x[j]
			x1 += self.bias
			x2 = (1+np.exp(-x1)) ** -1
			x3 = (y[0]-x2) ** 2
			print("loss = ",x3)
			x4 = x2*(1-x2)
			x5 = 2*(y[0]-x2)
			self.bias -= self.learning_rate * (x4*x5)
			for j in range(self.input_size):
				self.weights[j] -= self.learning_rate*(x4*x5*x[j])

	@classmethod

	def standardization(self,x):
		mean = np.mean(x,axis = 0)
		sd = np.var(x,axis = 0) ** .5
		for i in range(x.shape[1]):
			x[:,i] = x[:,i] - mean[i]
			x[:,i] = x[:,i] / sd[i]
		return x

data = pandas.read_csv("log_data.csv")
x = np.array(data.iloc[:, [2, 3]].values,dtype = np.float64)
y = data.iloc[:, [4]].values
x = perceptron.standardization(x)
new_p = perceptron(2,1,0.001)
new_p.train(x,y)


text:
import pandas as pd
import numpy as np

def clean_static_df(static_df):
    static_df_clean = static_df
    static_df_clean = pd.get_dummies(data=static_df_clean, columns=[
                                     'sex']).drop('sex_MALE', axis=1)

    static_df_clean.drop('discharge_destination', axis=1, inplace=True)

    static_df_clean.loc[:, 'datetime_min'] = pd.to_datetime(static_df_clean['datetime_min'],
        dayfirst=True)

    first_lab_time_col = pd.to_datetime(
        static_df_clean['datetime_min'], dayfirst=True)
    static_df_clean.loc[:, 'icu_in_year'] = first_lab_time_col.dt.year
    static_df_clean.loc[:, 'icu_in_month'] = first_lab_time_col.dt.month
    static_df_clean.loc[:, 'icu_in_day'] = first_lab_time_col.dt.day
    static_df_clean.drop('datetime_min', axis=1, inplace=True)
    static_df_clean.drop('datetime_max', axis=1, inplace=True)
    # TODO: why is there no hourly info?

    static_df_clean.set_index('patient_id')

    return static_df_clean

def clean_dynamic_df(dynamic_df, static_df):
    dynamic_df_times = dynamic_df.merge(static_df[['patient_id','datetime_min']],
        how='inner', on='patient_id')

    dynamic_df_times.loc[:, 'datetime_min'] = pd.to_datetime(dynamic_df_times['datetime_min'],
        dayfirst=True)
    dynamic_df_times.loc[:, 'datetime'] = pd.to_datetime(dynamic_df_times['datetime'],
        dayfirst=True)

    dynamic_df_times.loc[:, 'hours_in'] = \
         (dynamic_df_times['datetime'] - dynamic_df_times['datetime_min'])

    dynamic_df_times = dynamic_df_times.drop(['datetime','datetime_min'], axis=1)

    # upsample to hourly
    dynamic_df_hourly = dynamic_df_times.set_index('hours_in').groupby('patient_id').resample('H').mean()
    dynamic_df_hourly = dynamic_df_hourly.drop('patient_id', axis=1).reset_index()
#    dynamic_df_hourly = dynamic_df_hourly.set_index('patient_id')

    dynamic_df_clean = dynamic_df_hourly.set_index(['patient_id','hours_in']).dropna(how='all')

    return dynamic_df_clean

def clean_treat_df(treat_df, static_df):
    treat_df_times = treat_df.merge(static_df[['patient_id','datetime_min']],
        how='inner', on='patient_id')

    treat_df_times.loc[:, 'datetime_min'] = pd.to_datetime(treat_df_times['datetime_min'],
        dayfirst=True)
    treat_df_times.loc[:, 'datetime'] = pd.to_datetime(treat_df_times['date'],
        dayfirst=True)

    treat_df_times.loc[:, 'hours_in'] = \
         (treat_df_times['datetime'] - treat_df_times['datetime_min'])

    treat_df_times = treat_df_times.drop(['datetime','datetime_min'], axis=1)

    treat_df_hourly = treat_df_times.set_index('hours_in').groupby('patient_id').resample('H').mean()
    treat_df_hourly = treat_df_hourly.drop('patient_id', axis=1).reset_index()

    treat_df_clean = treat_df_hourly.set_index(['patient_id','hours_in']).dropna(how='all')

    return treat_df_clean

def clean_outcome_df(static_df):
    outcome_df = static_df[['patient_id','discharge_destination','datetime_min','datetime_max']]
    outcome_df.loc[:,'datetime_min'] = pd.to_datetime(outcome_df['datetime_min'])
    outcome_df.loc[:,'datetime_max'] = pd.to_datetime(outcome_df.loc[:,'datetime_max'])

    outcome_df['hours_in'] = outcome_df['datetime_max'] - outcome_df['datetime_min']

    outcome_df.loc[:, 'death'] = (outcome_df['discharge_destination'] == 'Fallecimiento') * 1.0

    outcome_df = outcome_df.drop(['discharge_destination','datetime_min','datetime_max'], axis=1)
    outcome_df = outcome_df.set_index('patient_id')

    return outcome_df

def get_CXYT_for_modeling(static_df, dynamic_df, treat_df, outcome_df, output_filepath):
    static_df = static_df.reset_index()
    dynamic_df = dynamic_df.reset_index()
    treat_df = treat_df.reset_index()
    outcome_df = outcome_df.reset_index()

    # upsample treatemnt to hourly before merging
    treat_df.loc[:, 'hours_in'] = pd.to_timedelta(treat_df['hours_in'])
    dynamic_df.loc[:, 'hours_in'] = pd.to_timedelta(dynamic_df['hours_in'])
    outcome_df.loc[:, 'hours_in'] = pd.to_timedelta(outcome_df['hours_in'])
    treat_df.set_index('hours_in').groupby('patient_id').resample('H').ffill()
    treat_df_upsampled = treat_df.\
        set_index('hours_in').groupby('patient_id').resample('H').\
        ffill().drop('patient_id', axis=1)

    dynamic_treat_outcome_df = dynamic_df.merge(
        treat_df_upsampled, how='outer',on=['patient_id', 'hours_in']).merge(
        outcome_df, how='outer',on=['patient_id', 'hours_in'])

    treat_cols = dynamic_treat_outcome_df.columns.str.contains('HIBOR')
    outcome_cols = dynamic_treat_outcome_df.columns == 'death'

    #fill in not treated timepoints and not dead timepoints with 0 before padding time
    dynamic_treat_outcome_df.loc[:,treat_cols | outcome_cols] = \
        dynamic_treat_outcome_df.loc[:,treat_cols | outcome_cols].fillna(0)

    date_first = dynamic_treat_outcome_df['hours_in'].min()
    date_first_trim = pd.Timedelta('-2 days')
    date_last = dynamic_treat_outcome_df['hours_in'].max()
    date_last_trim = pd.Timedelta('30 days 0 hours')
    print(date_first, date_last)

    mux = pd.MultiIndex.from_product([list(set(dynamic_treat_outcome_df['patient_id'])),
                                    pd.timedelta_range(date_first_trim, date_last_trim, freq='H')],
                                    names=['patient_id', 'hours_in_filled'])

    num_hours = len(pd.timedelta_range(date_first_trim, date_last_trim, freq='H'))
    #print(num_hours)
    print('size of new time-padded data ' + str(len(mux)))

    temp = dynamic_treat_outcome_df.set_index(['patient_id', 'hours_in'])
    result = temp.reindex(mux, fill_value=np.nan).reset_index()

    #fill in not treated timepoints and not dead timepoints with 0, after padding time too
    result.loc[:,treat_cols | outcome_cols] = result.loc[:,treat_cols | outcome_cols].fillna(0)

    # need to drop patients whose outcomes was after the max date cutoff
    death_times_dist = dynamic_treat_outcome_df.loc[dynamic_treat_outcome_df['death'] == 1,:][['patient_id','hours_in','death']]
    died_too_late = death_times_dist.loc[death_times_dist['hours_in'] > date_last_trim, :]['patient_id']
    print('omitting these patients' + str(died_too_late.values) + ' because their outcome times was after the date cutoff')

    XYT_df = result.loc[~result['patient_id'].isin(died_too_late),:]
    num_patients = len(set(XYT_df['patient_id']))

    CXYT_df = static_df.merge(XYT_df, on='patient_id', how='right')

    return CXYT_df


text:
import re
import numpy as np
import scipy.misc
import os
from scipy.misc import imresize

from PIL import Image, ImageDraw, ImageFont

def mkdir_p(path):
    try:
        os.makedirs(path)
    except OSError as exc:  # Python >2.5
        if exc.errno == errno.EEXIST and os.path.isdir(path):
            pass

def drawCaption(img, caption, max_len):
    img_txt = Image.fromarray(img)
    # get a font
    fnt = ImageFont.truetype('/eai/project/.fonts/FreeMono.ttf', 30)
    # get a drawing context
    d = ImageDraw.Draw(img_txt)

    d.text((10, 256), 'Stage-I', font=fnt, fill=(0, 0, 0, 255))
    d.text((10, 512), 'Stage-II', font=fnt, fill=(0, 0, 0, 255))
    d.text((10, 768), 'Stage-III', font=fnt, fill=(0, 0, 0, 255))

    caption = caption.split(' ')

    cap1 = ' '.join(caption[:max_len])
    cap2 = ' '.join(caption[max_len + 1:])
    d.text((256, 10), cap1, font=fnt, fill=(0, 0, 0, 255))
    d.text((256, 60), cap2, font=fnt, fill=(127, 127, 127, 255))

    return img_txt

def save_images_with_text(
        lr_sample_batchs, hr_sample_batchs, sr_sample_batchs,
        reals_batch, texts_batch, batch_size, max_len,
        startID, save_dir=None):

    if save_dir and not os.path.isdir(save_dir):
        print('Make a new folder: ', save_dir)
        mkdir_p(save_dir)

    # Save up to 16 samples for each text embedding/sentence
    img_shape = sr_sample_batchs[0][0].shape
    super_images = []
    for i in range(batch_size):
        if not re.search('[a-zA-Z]+', texts_batch[i]):
            continue

        padding = 255 + np.zeros(img_shape)
        row1 = [padding]
        row2 = [padding]
        row3 = [padding]

        for j in range(lr_sample_batchs[0].shape[0]):
            lr_img = lr_sample_batchs[i][j]
            hr_img = hr_sample_batchs[i][j]
            sr_img = sr_sample_batchs[i][j]
            if j == 0:
                row1.append(imresize(reals_batch[0][i], sr_img.shape[:2]))
                row2.append(imresize(reals_batch[1][i], sr_img.shape[:2]))
                row3.append(imresize(reals_batch[2][i], sr_img.shape[:2]))

            lr_re_sample = imresize(lr_img, sr_img.shape[:2])
            hr_re_sample = imresize(hr_img, sr_img.shape[:2])
            row1.append(lr_re_sample)
            row2.append(hr_re_sample)
            row3.append(sr_img)

        row1 = np.concatenate(row1, axis=1)
        row2 = np.concatenate(row2, axis=1)
        row3 = np.concatenate(row3, axis=1)
        superimage = np.concatenate([row1, row2, row3], axis=0)

        top_padding = 255 + np.zeros((128, superimage.shape[1], 3))
        superimage = np.concatenate([top_padding, superimage], axis=0)

        fullpath = '%s/sentence_%04d.jpg' % (save_dir, startID + i)
        superimage = drawCaption(np.uint8(superimage), texts_batch[i], max_len)

        if save_dir:
            scipy.misc.imsave(fullpath, superimage)
        super_images.append(superimage)

    return super_images

reason: import numpy,import scipy,from scipy