text:
from __future__ import annotations

from types import MappingProxyType
from typing import Any, Dict

import meshio
import numpy as np
import pyvista as pv

from .._doc import DocFormatterMeta, doc
from ..region_markers import RegionMarkerList

@doc(prefix='Generic mesh class', dim_points='n', dim_cells='j')
class Mesh(object, metaclass=DocFormatterMeta):
    """{prefix}.

    Depending on the number of dimensions of the cells, the appropriate
    subclass will be chosen if possible.

    Parameters
    ----------
    points : (m, {dim_points}) numpy.ndarray[float]
        Array with points.
    cells : (i, {dim_cells}) numpy.ndarray[int]
        Index array describing the cells of the mesh.
    fields : Dict[str, int]:
        Mapping from field names to labels
    region_markers : RegionMarkerList, optional
        List of region markers used for assigning labels to regions.
        Defaults to an empty list.
    **cell_data
        Additional cell data. Argument must be a 1D numpy array
        matching the number of cells defined by `i`.
    """
    _registry: Dict[int, Any] = {}
    cell_type: str = 'generic'

    def __init_subclass__(cls, cell_dim: int, **kwargs):
        super().__init_subclass__(**kwargs)
        cls._registry[cell_dim] = cls

    def __new__(cls, points: np.ndarray, cells: np.ndarray, *args, **kwargs):
        cell_dim = cells.shape[1]
        subclass = cls._registry.get(cell_dim, cls)
        return super().__new__(subclass)

    def __init__(self,
                 points: np.ndarray,
                 cells: np.ndarray,
                 *,
                 fields: Dict[str, int] = None,
                 region_markers: RegionMarkerList = None,
                 **cell_data):
        default_key = 'physical'
        if (not cell_data) or (default_key in cell_data):
            self.default_key = default_key
        else:
            self.default_key = list(cell_data.keys())[0]

        self.fields = dict(fields) if fields else {}

        self.points = points
        self.cells = cells
        self.field_to_number = MappingProxyType(self.fields)
        self.region_markers = RegionMarkerList()
        if region_markers:
            self.region_markers.extend(region_markers)
        self.cell_data = cell_data

    def __repr__(self, indent: int = 0):
        """Canonical string representation."""
        region_markers = set(m.name if m.name else m.label
                             for m in self.region_markers)
        s = (
            f'{self.__class__.__name__}(',
            f'    points = {self.points.shape},',
            f'    cells = {self.cells.shape},',
            f'    fields = {tuple(self.field_to_number.keys())},',
            f'    region_markers = {tuple(region_markers)},',
            f'    cell_data = {tuple(self.cell_data.keys())},',
            ')',
        )

        prefix = ' ' * indent
        return f'\n{prefix}'.join(s)

    @property
    def number_to_field(self):
        """Mapping from numbers to fields, proxy to
        :attr:`{classname}.field_to_number`."""
        return MappingProxyType(
            {v: k
             for k, v in self.field_to_number.items()})

    def to_meshio(self) -> 'meshio.Mesh':
        """Return instance of :func:`meshio.Mesh`."""
        cells = [
            (self.cell_type, self.cells),
        ]

        mesh = meshio.Mesh(self.points, cells)

        for key, value in self.cell_data.items():
            mesh.cell_data[key] = [value]

        return mesh

    @classmethod
    def from_meshio(cls, mesh: 'meshio.Mesh'):
        """Return :class:`{classname}` from meshio object."""
        points = mesh.points
        cells = mesh.cells[0].data
        cell_data = {}

        for key, value in mesh.cell_data.items():
            # PyVista chokes on ':ref' in cell_data
            key = key.replace(':ref', '-ref')
            cell_data[key] = value[0]

        return Mesh(points=points, cells=cells, **cell_data)

    def write(self, *args, **kwargs):
        """Simple wrapper around :func:`meshio.write`."""
        self.to_meshio().write(*args, **kwargs)

    @classmethod
    def read(cls, filename, **kwargs):
        """Simple wrapper around :func:`meshio.read`."""
        mesh = meshio.read(filename, **kwargs)
        return cls.from_meshio(mesh)

    def to_pyvista_unstructured_grid(self) -> 'pv.PolyData':
        """Return instance of :class:`pyvista.UnstructuredGrid`.

        References
        ----------
        https://docs.pyvista.org/core/point-grids.html#pv-unstructured-grid-class-methods
        """
        return pv.from_meshio(self.to_meshio())

    def plot(self, **kwargs):
        raise NotImplementedError(
            f'Not implemented for {self.__class__.__name__}, '
            'use one of the specific subclasses.')

    def plot_mpl(self, **kwargs):
        raise NotImplementedError(
            f'Not implemented for {self.__class__.__name__}, '
            'use one of the specific subclasses.')

    def plot_itk(self, **kwargs):
        """Wrapper for :func:`pyvista.plot_itk`.

        Parameters
        ----------
        **kwargs
            These parameters are passed to :func:`pyvista.plot_itk`
        """
        return pv.plot_itk(self.to_meshio(), **kwargs)

    def plot_pyvista(self, **kwargs):
        """Wrapper for :func:`pyvista.plot`.

        Parameters
        ----------
        **kwargs
            These parameters are passed to :func:`pyvista.plot`
        """
        return pv.plot(self.to_meshio(), **kwargs)

    @property
    def cell_centers(self):
        """Return centers of cells (mean of corner points)."""
        return self.points[self.cells].mean(axis=1)

    def get_cell_data(self, key: str, default_value: float = 0) -> np.ndarray:
        """Get cell data with optional default value.

        Parameters
        ----------
        key : str
            Key of the cell data to retrieve.
        default_value : float, optional
            Optional default value (if cell data does not exist)

        Returns
        -------
        numpy.ndarray
        """
        if not key:
            key = self.default_key

        try:
            return self.cell_data[key]
        except KeyError:
            return np.ones(len(self.cells), dtype=int) * default_value

    @property
    def zero_labels(self) -> np.ndarray:
        """Return zero labels as fallback."""
        return np.zeros(len(self.cells), dtype=int)

    @property
    def labels(self) -> np.ndarray:
        """Shortcut for cell labels."""
        try:
            return self.cell_data[self.default_key]
        except KeyError:
            return self.zero_labels

    @labels.setter
    def labels(self, data: np.ndarray):
        """Shortcut for setting cell labels.

        Updates :attr:`{classname}.cell_data`.
        """
        self.cell_data[self.default_key] = data

    @property
    def dimensions(self) -> int:
        """Return number of dimensions for point data."""
        return self.points.shape[1]

    def reverse_cell_order(self):
        """Reverse order of cells and cell data.

        Updates :attr:`{classname}.cell_data`.
        """
        self.cells = self.cells[::-1]
        for key, data in self.cell_data.items():
            self.cell_data[key] = data[::-1]

    def remove_cells(self, *, label: int, key: str = None):
        """Remove cells with cell data matching given label.

        Parameters
        ----------
        label : int
            All cells with this label will be removed.
        key : str, optional
            The key of the cell data to use.
            Uses :attr:`Mesh.default_key` if None.
        """
        if not key:
            key = self.default_key
        idx = self.cell_data[key] != label

        for k, v in self.cell_data.items():
            self.cell_data[k] = v[idx]
            pass

        self.cells = self.cells[idx]
        self.remove_loose_points()

    def remove_loose_points(self):
        """Remove points that do not belong to any cells."""
        cell_indices = np.unique(self.cells)
        self.points = np.take(self.points, cell_indices, axis=0)
        self._regenerate_cell_indices(cell_indices)

    def _regenerate_cell_indices(self, indices: np.ndarray = None):
        """Re-generate cell indices to remove gaps, i.e., 1,3,5 -> 1,2,3.

        Parameters
        ----------
        indices : np.ndarray
            Current list of cell indices.
        """
        if indices is None:
            indices = np.unique(self.cells)
        indices = indices.ravel()

        mapping = np.vstack([indices, np.arange(len(indices))])

        shape = self.cells.shape
        new_cells = self.cells.ravel()

        mask = np.in1d(new_cells, mapping[0, :])
        new_cells[mask] = mapping[1,
                                  np.searchsorted(mapping[
                                      0, :], new_cells[mask])]
        self.cells = new_cells.reshape(shape)

    def crop(
        self,
        xmin: float = -np.inf,
        xmax: float = np.inf,
        ymin: float = -np.inf,
        ymax: float = np.inf,
        zmin: float = -np.inf,
        zmax: float = np.inf,
        include_partial: bool = False,
    ):
        """Crop mesh to given region.

        Parameters
        ----------
        xmin : float, optional
            Minimum x value.
        xmax : float, optional
            Maximum x value.
        ymin : float, optional
            Minimum y value.
        ymax : float, optional
            Maximum y value.
        zmin : float, optional
            Minimum z value (3D point data only).
        zmax : float, optional
            Maximum z value (3D point data only).
        include_partial : bool, optional
            If True, include cells that are partially inside the
            given crop region, i.e. one of its points is inside.

        Returns
        -------
        cropped_mesh : :class:`{classname}`
            Cropped mesh
        """
        points = self.points
        cells = self.cells
        cell_data = self.cell_data

        cell_coords = points[cells]
        dim = points.shape[1]

        if dim not in (2, 3):
            raise NotImplementedError('Cropping not supported on '
                                      f'{dim}d data ({self.points.shape=})')

        if dim >= 2:
            x_idx = (xmin <= cell_coords[:, :, 0]) & (cell_coords[:, :, 0] <=
                                                      xmax)
            y_idx = (ymin <= cell_coords[:, :, 1]) & (cell_coords[:, :, 1] <=
                                                      ymax)
            idx = x_idx & y_idx

        if dim == 3:
            z_idx = (zmin <= cell_coords[:, :, 2]) & (cell_coords[:, :, 2] <=
                                                      zmax)
            idx = idx & z_idx

        f = np.any if include_partial else np.all
        cells_to_keep = f(idx, axis=1)

        new_cells = cells[cells_to_keep]
        new_cell_data = {}

        for name, data in cell_data.items():
            new_cell_data[name] = data[cells_to_keep]

        new_mesh = self.__class__(points=points,
                                  cells=new_cells,
                                  **new_cell_data)
        new_mesh.remove_loose_points()
        return new_mesh

max_stars_repo_name: hpgem/nanomesher

max_issues_repo_name: hpgem/nanomesher

max_forks_repo_name: hpgem/nanomesher


text:
# encoding:utf/8
import sys
from mmdet.apis import inference_detector, init_detector
import json
import os
import numpy as np
import argparse
from tqdm import tqdm

class MyEncoder(json.JSONEncoder):
    def default(self, obj):
        if isinstance(obj, np.integer):
            return int(obj)
        elif isinstance(obj, np.floating):
            return float(obj)
        elif isinstance(obj, np.ndarray):
            return obj.tolist()
        else:
            return super(MyEncoder, self).default(obj)

        # generate result

def result_from_dir():
    index = {1: 1, 2: 9, 3: 5, 4: 5, 5: 4, 6: 2, 7: 8, 8: 6, 9: 10, 10: 7}
    # build the model from a config file and a checkpoint file
    model = init_detector(config2make_json, model2make_json, device='cuda:0')
    pics = os.listdir(pic_path)
    meta = {}
    images = []
    annotations = []
    num = 0
    for im in tqdm(pics):
        num += 1
        img = os.path.join(pic_path, im)
        result_ = inference_detector(model, img)
        images_anno = {}
        images_anno['file_name'] = im
        images_anno['id'] = int(num)
        images.append(images_anno)
        for i, boxes in enumerate(result_, 1):
            if len(boxes):
                defect_label = index[i]
                for box in boxes:
                    anno = {}
                    anno['image_id'] = int(num)
                    anno['category_id'] = defect_label
                    anno['bbox'] = [round(float(i), 2) for i in box[0:4]]
                    anno['bbox'][2] = anno['bbox'][2] - anno['bbox'][0]
                    anno['bbox'][3] = anno['bbox'][3] - anno['bbox'][1]
                    anno['score'] = float(box[4])
                    annotations.append(anno)
    meta['images'] = images
    meta['annotations'] = annotations
    with open(json_out_path, 'w') as fp:
        json.dump(meta, fp, cls=MyEncoder, indent=4, separators=(',', ': '))
    res = check_json_with_error_msg(json_out_path)
    print(res)

def check_json_with_error_msg(pred_json, num_classes=10):
    '''
    Args:
        pred_json (str): Json path
        num_classes (int): number of foreground categories
    Returns:
        Message (str)
    Example:
        msg = check_json_with_error_msg('./submittion.json')
        print(msg)
    '''
    if not pred_json.endswith('.json'):
        return "the prediction file should ends with .json"
    with open(pred_json) as f:
        data = json.load(f)
    if not isinstance(data, dict):
        return "the prediction data should be a dict"
    if not 'images' in data:
        return "missing key \"images\""
    if not 'annotations' in data:
        return "missing key \"annotations\""
    images = data['images']
    annotations = data['annotations']
    if not isinstance(images, (list, tuple)):
        return "\"images\" format error"
    if not isinstance(annotations, (list, tuple)):
        return "\"annotations\" format error"
    for image in images:
        if not 'file_name' in image:
            return "missing key \"file_name\" in \"images\""
        if not 'id' in image:
            return "missing key \"id\" in \"images\""
    for annotation in annotations:
        if not 'image_id' in annotation:
            return "missing key \"image_id\" in \"annotations\""
        if not 'category_id' in annotation:
            return "missing key \"category_id\" in \"annotations\""
        if not 'bbox' in annotation:
            return "missing key \"bbox\" in \"annotations\""
        if not 'score' in annotation:
            return "missing key \"score\" in \"annotations\""
        if not isinstance(annotation['bbox'], (tuple, list)):
            return "bbox format error"
        if len(annotation['bbox'])==0:
            return "empty bbox"
        if annotation['category_id'] > num_classes or annotation['category_id'] < 0:
            return "category_id out of range"
    return "OK"

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="Generate result")
    parser.add_argument("-m", "--model", help="Model path", type=str, )
    parser.add_argument("-c", "--config", help="Config path", type=str, )
    parser.add_argument("-im", "--im_dir", help="Image path", type=str,
                        default='/workspace/meijinpeng/dataset/tianchi1/chongqing1_round1_testA_20191223/images')
    parser.add_argument('-o', "--out", help="Save path", type=str,
                        default='submit.json')
    args = parser.parse_args()
    model2make_json = args.model
    config2make_json = args.config
    json_out_path = args.out
    pic_path = args.im_dir
    result_from_dir()

ext: py

reason: import numpy


text:
#Requiered Libraries
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button,CheckButtons

class fourier:
    #PARAMETERS OF WAVE FUNCTIONS
    def __init__(self,dim,Amplitude,Frequency,init_time,end_time,dt):
        self.dim = dim
        self.Amplitude = Amplitude
        self.Frequency = Frequency
        self.init_time = init_time
        self.end_time = end_time
        self.dt = dt

    #Sin/Cos Wave Function
    def wave(self):
        time_array = np.arange(self.init_time,self.end_time,self.dt)
            #[time_array,wave_points]
        return [time_array,(self.Amplitude)*np.sin(2*np.pi*(self.Frequency)*time_array)]

    #WRAP FUNCTION

#######################################3
#GENERO MIS ONDAS
w_1 = fourier(8,1,10,0,2,0.001)
w_2 = fourier(3,1,25,0,2,0.001)

t = w_1.wave()[0]
y1 = w_1.wave()[1]
y2 = w_2.wave()[1]
A1 = w_1.Amplitude
A2 = w_2.Amplitude

###############################3
class wrap_plot:
    def __init__(self,w_array,Amplitude_array,j):
        #PARAMETERS
        prec = 1/(1.1e3)
        freq_min = 1
        freq_max = 100
        freq_loop = np.linspace(freq_min,freq_max,5000)
        factor = 1/prec
        freq_wrap = 3
        #CENTER OF MASS
        wrap = np.sum(w_array,axis=0)*np.exp(2j*np.pi*t*freq_wrap)
        ms_sum = np.sum(wrap)/factor
        ms_abs = np.abs(ms_sum)
        ms_x = []
        ms_y = []
        #line plot of center mass
        for freq_wrap in freq_loop:
            arg = np.sum(w_array,axis=0)*np.exp(-2j*np.pi*t*freq_wrap)
            ms_x.append(freq_wrap)
            ms_y.append(np.abs(np.sum(arg.real)/factor+np.sum(arg.imag)/factor))

        # Create the figure and the line that we will manipulate
        plt.style.use('dark_background')
        fig, axs = plt.subplots(1,2)

        plt.subplots_adjust(bottom=.25)

        if j == 0:
            p1, = axs[0].plot(t,w_array[0],color = '#ae00ff',label='WAVE 1') # COMMA CONNECT WITH SLIDER
            p2, = axs[0].plot(t,w_array[1],color = '#ff8c00',label='WAVE 2') # COMMA CONNECT WITH SLIDER
            p, = axs[0].plot(t,np.sum(w_array,axis=0),color = 'white',label = 'SUM OF WAVES')
            w, = axs[1].plot(wrap.real,wrap.imag,color = 'white')
            ms, = axs[1].plot(ms_sum.real,ms_sum.imag,'ro')
            #PLOT 0:
            ax_slide1 = plt.axes([0.25, 0.15, 0.65, 0.03]) # CREATE THE SLIDER BAR 1
            ax_slide2 = plt.axes([0.25, 0.1, 0.65, 0.03]) # CREATE THE SLIDER BAR 2
            ax_slide_wrap = plt.axes([0.25, 0.05, 0.65, 0.03]) # CREATE THE SLIDER BAR 3
            freq1_factor = Slider(ax_slide1,"FREQ1",valmin=0.2,valmax=60.0,valinit=w_1.Frequency,color ='#ae00ff')
            freq2_factor = Slider(ax_slide2,"FREQ2",valmin=0.2,valmax=60.0,valinit=w_2.Frequency,color ='#ff8c00')
            freq_wrap_factor = Slider(ax_slide_wrap,"FREQ_WRAP",valmin=0.2,valmax=60.0,valinit=1,color ='white')
            axs[0].set_title('Harmonics Waves',fontsize = 16)
            axs[0].set_xlabel('Time Domain',fontsize = 13)
            axs[0].set_ylabel('Amplitude',fontsize = 13)
            axs[1].set_title('Wrapping the Unknown Wave',fontsize = 16)
            axs[1].set_xlabel('Real Axis',fontsize = 13)
            axs[1].set_ylabel('Imaginary Axis',fontsize = 13)
            axs[1].set_xlim(-2,2)
            axs[1].set_ylim(-2,2)

            ########################################
            #            BUTTONS
            # Make checkbuttons with all plotted lines with correct visibility
            lines = [p, p1, p2]
            col = (0.4,0.3,0.75,0.4)
            rax = plt.axes([0, 0, 0.15, 0.15],facecolor=col)
            labels = [str(line.get_label()) for line in lines]
            visibility = [line.get_visible() for line in lines]
            check = CheckButtons(rax, labels, visibility)

            def func(label):
                index = labels.index(label)
                lines[index].set_visible(not lines[index].get_visible())
                plt.draw()

            check.on_clicked(func)
            def update(val):
                f1 = freq1_factor.val
                f2 = freq2_factor.val
                freq_wrap = freq_wrap_factor.val
                p1.set_ydata(1*np.sin(2*np.pi*f1*t))
                p2.set_ydata(1*np.sin(2*np.pi*f2*t))
                p.set_ydata(1*np.sin(2*np.pi*f1*t) + 1*np.sin(2*np.pi*f2*t))
                y1 = Amplitude_array[0]*np.sin(2*np.pi*f1*t)
                y2 = Amplitude_array[1]*np.sin(2*np.pi*f2*t)
                arg = (y1+y2)*np.exp(-2j*np.pi*t*freq_wrap)
                w.set_xdata(arg.real)
                w.set_ydata(arg.imag)
                ms.set_xdata(np.sum(arg.real)/factor)
                ms.set_ydata(np.sum(arg.imag)/factor)

            freq1_factor.on_changed(update)
            freq2_factor.on_changed(update)
            freq_wrap_factor.on_changed(update)
        else:
            w, = axs[0].plot(wrap.real,wrap.imag,color = 'white')
            ms, = axs[0].plot(ms_sum.real,ms_sum.imag,'ro')
            FT, = axs[1].plot(ms_sum.real,ms_sum.imag,'ro')
            msxy = axs[1].plot(ms_x,ms_y/np.max(ms_y),'#8c5eff')
            ax_slide_wrap = plt.axes([0.25, 0.05, 0.65, 0.03]) # CREATE THE SLIDER BAR 2
            freq_wrap_factor = Slider(ax_slide_wrap,"FREQ_WRAP",valmin=0.2,valmax=60.0,valinit=1,color ='white')
            axs[0].set_title('Wrapping the Unknown Wave',fontsize = 16)
            axs[0].set_xlabel('Real Axis',fontsize = 13)
            axs[0].set_ylabel('Imaginary Axis',fontsize = 13)
            axs[0].set_xlim(-2,2)
            axs[0].set_ylim(-2,2)
            axs[1].set_title('Fourier Transform',fontsize = 16)
            axs[1].set_xlabel('Frequency Domain',fontsize = 13)
            axs[1].set_ylabel('Weights',fontsize = 13)
            axs[1].set_xlim(0,100)
            axs[1].set_ylim(0,1)
            ######################################################
            def update(val):
                freq_wrap = freq_wrap_factor.val
                arg = np.sum(w_array,axis=0)*np.exp(-2j*np.pi*t*freq_wrap)
                FT.set_xdata(freq_wrap)
                FT.set_ydata(np.abs(np.sum(arg.real)/factor+np.sum(arg.imag)/factor))
                w.set_xdata(arg.real)
                w.set_ydata(arg.imag)
                ms.set_xdata(np.sum(arg.real)/factor)
                ms.set_ydata(np.sum(arg.imag)/factor)

            freq_wrap_factor.on_changed(update)
        plt.show()

wrap_plot([y1,y2],[A1,A2],0)

max_stars_repo_name: LopezBanos/PyFourier

max_issues_repo_name: LopezBanos/PyFourier

max_forks_repo_name: LopezBanos/PyFourier


text:
from __future__ import print_function

from nose.tools import eq_, ok_

import nose, functools
import os
import numpy as np

import load
from utils import dna_io

def expected_failure(test):
    @functools.wraps(test)
    def inner(*args, **kwargs):
        try:
            test(*args, **kwargs)
        except Exception:
            raise nose.SkipTest
        else:
            raise AssertionError('Failure expected')
    return inner

### DNase data fixtures
train_size = 3201397
valid_size = 70000
output_size = 164
chunk_size = 4096
batch_size = 128
num_chunks_train = train_size // chunk_size

full_train_shape = (train_size, 4, 1, 600)
chunk_train_shape = (chunk_size,4,1,600)
chunk_out_shape = (chunk_size,output_size)
batch_train_shape = (batch_size,4,1,600)

path = os.path.expanduser("~/projects/SeqDemote/data/DNase/encode_roadmap.h5")

### Kmer fixtures

aaa_3mer_string = 'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAA'
aaa_3mer_position = 0

ttt_3mer_string = 'TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT'
ttt_3mer_position = 63

aaa_4mer_string = 'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA'
ttt_4mer_string = 'TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT'

rando_3mer_string = 'ATGGGGTAGAGAATGGGGTAGAGACCAGGT'

alphabet_size_3 = 64
alphabet_size_4 = 256

threemer_chunk_shape = (chunk_size,alphabet_size_3, 1, 598)
fourmer_chunk_shape = (chunk_size, alphabet_size_4, 1, 597)

def test_build_data_loader():
    """ Can I build a data loader for the DNase data """

    data_loader = load.StandardDataLoader(chunk_size=chunk_size, batch_size=batch_size, num_chunks_train=num_chunks_train)
    data_loader.load_train()
    eq_(data_loader.train_in.shape, full_train_shape)

def test_build_data_loader_kwargs():
    """ Can I build a data loader for the DNase data specifying the data path """

    data_loader = load.StandardDataLoader(data_path=path,chunk_size=chunk_size, batch_size=batch_size, num_chunks_train=num_chunks_train)
    data_loader.load_train()
    eq_(data_loader.train_in.shape, full_train_shape)
    eq_(data_loader.chunk_size, chunk_size)
    eq_(data_loader.num_chunks_train, num_chunks_train)

def test_dnase_data_shape():
    """ Is my DNase data the right size and shape """

    data_loader = load.StandardDataLoader(data_path=path)
    data_loader.load_train()
    num_chunks = range(num_chunks_train)
    for e, (x_chunk, y_chunk) in zip(num_chunks,data_loader.create_batch_gen()):
        print("Chunk: ", str(e + 1), " of ", num_chunks_train)
        eq_(x_chunk.shape, chunk_train_shape)
        eq_(y_chunk.shape, chunk_out_shape)

def test_exhaust_data():
    """ If I iterate through all the chunks, how many data points do I see? """
    seen_pts = 0

    data_loader = load.StandardDataLoader(data_path=path)
    data_loader.load_train()
    num_chunks = range(num_chunks_train)
    for e, (x_chunk, y_chunk) in zip(num_chunks,data_loader.create_batch_gen()):
        seen_pts = seen_pts + x_chunk.shape[0]
    print("Saw ", str(seen_pts), " points total")
    ok_(seen_pts <= train_size)

### Semantic tests: are all examples properly one-hot encoded?

def test_training_batch_encoding_sum():
    """ If I sum all elements of a chunk of training data, do I get the number of expected ones? """

    data_loader = load.StandardDataLoader(data_path=path)
    data_loader.load_train()
    num_chunks = range(num_chunks_train)
    expected_chunk_sum = 600 * chunk_size
    ### each chunk of should have |seq_length| * |batch_size| * |chunk_size| number of 1s
    for e, (x_chunk, y_chunk) in zip(num_chunks,data_loader.create_batch_gen()):
        ones_per_chunk = np.sum(x_chunk)
        #print("Chunk  ", str(e + 1), " has ", ones_per_chunk, " bits turned on")
        eq_(ones_per_chunk, expected_chunk_sum)

def test_validation_batch_encoding_sum():
    """ If I sum all elements of a chunk of validation data, do I get the number of expected ones? """

    data_loader = load.StandardDataLoader(data_path=path)
    data_loader.load_validation()
    num_chunks = range(num_chunks_train)
    expected_chunk_sum = 600 * chunk_size
    ### each chunk of should have |seq_length| * |batch_size| * |chunk_size| number of 1s
    for e, (x_chunk, y_chunk) in zip(num_chunks,data_loader.create_valid_gen()):
        ones_per_chunk = np.sum(x_chunk)
        #print("Chunk  ", str(e + 1), " has ", ones_per_chunk, " bits turned on")
        eq_(ones_per_chunk, expected_chunk_sum)

def test_any_always_negatives_training():
    """ Are there any always-negative examples in the training set? """

    data_loader = load.StandardDataLoader(data_path=path)
    data_loader.load_train()
    for e, (x_chunk, y_chunk) in zip(range(num_chunks_train),data_loader.create_batch_gen()):
        for label_vector in y_chunk:
            total_peaks_on = np.sum(label_vector)
            ok_(total_peaks_on > 0)

def test_any_always_negatives_validation():
    """ Are there any always-negative examples in the validation set? """

    data_loader = load.StandardDataLoader(data_path=path)
    data_loader.load_validation()
    for e, (x_chunk, y_chunk) in zip(range(num_chunks_valid),data_loader.create_valid_gen()):
        for label_vector in y_chunk:
            total_peaks_on = np.sum(label_vector)
            ok_(total_peaks_on > 0)

def test_encoding_kmerizing_threemers():

    # encode the fixtures as a one-hot encoded matrix
    my_seqs = [aaa_3mer_string,ttt_3mer_string,aaa_3mer_string,ttt_3mer_string]
    seqs = np.vstack([dna_io.dna_one_hot(seq) for seq in my_seqs])

    # reshape into expected shape
    seqs = seqs.reshape((seqs.shape[0],4,1,seqs.shape[1]/4))

    # re-encode the one-hot examples as positional 3-mers
    seqs_3mers = dna_io.one_hot_to_kmerized(seqs, 3)
    seqs_3mers_reshaped = seqs_3mers.reshape((seqs_3mers.shape[0],alphabet_size_3,1,len(aaa_3mer_string)-2))

    # make up the positional 3-mer fixture
    aaa_vec = dna_io.dna_one_hot_kmer(aaa_3mer_string, 3)
    aaa_mat = aaa_vec.reshape((1,alphabet_size_3,1,len(aaa_3mer_string) - 2))
    ttt_vec = dna_io.dna_one_hot_kmer(ttt_3mer_string, 3)
    ttt_mat = ttt_vec.reshape((1,alphabet_size_3,1,len(ttt_3mer_string) - 2))

    fixture = np.vstack((aaa_mat,ttt_mat,aaa_mat,ttt_mat))

    for fix, elem in zip(fixture,seqs_3mers_reshaped):
        ok_(np.allclose(fix, elem))

def test_decoding_kmerizing_threemers():
    # encode the fixtures as a one-hot encoded matrix
    my_seqs = [aaa_3mer_string,ttt_3mer_string,aaa_3mer_string,ttt_3mer_string]
    seqs = np.vstack([dna_io.dna_one_hot(seq) for seq in my_seqs])

    # reshape into expected shape
    seqs = seqs.reshape((seqs.shape[0],4,1,seqs.shape[1]/4))

    # re-encode the one-hot examples as positional 3-mers
    seqs_3mers = dna_io.one_hot_to_kmerized(seqs, 3)
    seqs_3mers_reshaped = seqs_3mers.reshape((seqs_3mers.shape[0],alphabet_size_3,1,len(aaa_3mer_string)-2))

    # de-code the encoded 3-mers
    decoded_seqs = dna_io.kmer_vecs_to_dna(seqs_3mers_reshaped, 3)
    for seq, fix in zip(decoded_seqs,my_seqs):
        eq_(seq,fix)

def test_encoding_kmerizing_fourmers():

    # encode the fixtures as one-hot encoded matrix
    my_seqs = [aaa_4mer_string,ttt_4mer_string, aaa_4mer_string, ttt_4mer_string]
    seqs = np.vstack([dna_io.dna_one_hot(seq) for seq in my_seqs])

    # reshape into expected shape
    seqs = seqs.reshape((seqs.shape[0], 4, 1, seqs.shape[1]/4))

    # re-encode the one-hot examles as positional 4-mers
    seqs_4mers = dna_io.one_hot_to_kmerized(seqs, 4)
    seqs_4mers_reshaped = seqs_4mers.reshape((seqs_4mers.shape[0],alphabet_size_4,1,len(aaa_4mer_string)-3))

    # make up the positional 4-mer fixture
    aaa_vec = dna_io.dna_one_hot_kmer(aaa_4mer_string,4)
    aaa_mat = aaa_vec.reshape((1,alphabet_size_4,1,len(aaa_4mer_string)-3))
    ttt_vec = dna_io.dna_one_hot_kmer(ttt_4mer_string, 4)
    ttt_mat = ttt_vec.reshape((1,alphabet_size_4,1,len(ttt_4mer_string)-3))

    fixture = np.vstack((aaa_mat,ttt_mat,aaa_mat,ttt_mat))

    for fix, elem in zip(fixture, seqs_4mers_reshaped):
        ok_(np.allclose(fix,elem))

def test_decoding_kmerizing_fourmers():

    # encode the fixtures as one-hot encoded matrix
    my_seqs = [aaa_4mer_string, ttt_4mer_string, aaa_4mer_string, ttt_4mer_string]
    seqs = np.vstack([dna_io.dna_one_hot(seq) for seq in my_seqs])

    # reshape into expected shape
    seqs = seqs.reshape((seqs.shape[0], 4, 1, seqs.shape[1]/4))

    # re-encode the one-hot examples as positional 4-mers
    seqs_4mers = dna_io.one_hot_to_kmerized(seqs,4)
    seqs_4mers_reshaped = seqs_4mers.reshape((seqs_4mers.shape[0], alphabet_size_4, 1, len(aaa_4mer_string)-3))

    # decode the encoded 4-mers
    decoded_seqs = dna_io.kmer_vecs_to_dna(seqs_4mers_reshaped, 4)
    for seq, fix in zip(decoded_seqs,my_seqs):
        eq_(seq,fix)

def test_threemer_data_loader_shape():
    """ Can I build a data loader for the kmerized DNase data specifying the data path? """
    data_loader = load.KmerDataLoader(data_path=path, kmer_length=3)
    data_loader.load_train()
    num_chunks = range(num_chunks_train)
    for e, (x_chunk, y_chunk) in zip(num_chunks,data_loader.create_batch_gen()):
        print("Chunk: ", str(e + 1), " of ", num_chunks_train)
        eq_(x_chunk.shape, threemer_chunk_shape)
        eq_(y_chunk.shape, chunk_out_shape)

def test_fourmer_data_loader_shape():
    """ Can I build a data loader for the kmerized DNase data specifying the data path? """
    data_loader = load.KmerDataLoader(data_path=path, kmer_length=4)
    data_loader.load_train()
    num_chunks = range(num_chunks_train)
    for e, (x_chunk, y_chunk) in zip(num_chunks,data_loader.create_batch_gen()):
        print("Chunk: ", str(e + 1), " of ", num_chunks_train)
        eq_(x_chunk.shape, fourmer_chunk_shape)
        eq_(y_chunk.shape, chunk_out_shape)

def test_exhaust_threemer_data_loader():
    """ Is my kmerized data loader generator producing the right shape chunks? """
    seen_pts = 0
    data_loader = load.KmerDataLoader(data_path=path, kmer_length=3)
    data_loader.load_train()
    num_chunks = range(num_chunks_train)
    for e, (x_chunk, y_chunk) in zip(num_chunks,data_loader.create_batch_gen()):
        seen_pts = seen_pts + x_chunk.shape[0]
    print("Saw ", str(seen_pts), " points total")
    ok_(seen_pts <= train_size)

if __name__ == "__main__":
    #test_build_data_loader()
    #test_build_data_loader_kwargs()
    #test_dnase_data_shape()
    #test_exhaust_data()
    #test_training_batch_encoding_sum()
    #test_validation_batch_encoding_sum()
    #test_encoding_kmerizing_threemers()
    #test_decoding_kmerizing_threemers()
    #test_encoding_kmerizing_fourmers()
    #test_decoding_kmerizing_fourmers()
    test_threemer_data_loader_shape()

max_stars_repo_name: lzamparo/SeqDemote

max_issues_repo_name: lzamparo/SeqDemote

max_forks_repo_name: lzamparo/SeqDemote


text:
# Pablo Gainza 2019 - Assign the raw data from Intpred's webserver to each PDB's b factor column.
# Also, verify the data.
import numpy
from Bio.PDB import *
from IPython.core.debugger import set_trace

def parse_intpred_raw(pdbid, filename):
    ipfile = open(filename)
    lines = ipfile.readlines()
    line_ix = 0
    # First line to verify pdbid - to make sure I made no mistakes.
    verify = lines[0]
    verify = lines[0].split()
    verify_pdb = (''.join(verify[5][0:4])).upper()
    verify_chain = ''.join(verify[5][4])
    verify = '{}_{}'.format(verify_pdb,verify_chain)
    assert(verify == pdbid)
    print(verify)
    all_scores = {}

    # Ignore lines until we reach 'Residue Score'
    while "Residue Score" not in lines[line_ix]:
        line_ix += 1
    line_ix += 1
    for line in lines[line_ix:]:
        fields = line.split()
        if len(fields) > 0:
            chain_id = fields[0][0]
            residue = int(fields[0][1:])
            score = float(fields[1])
            all_scores[(chain_id, residue)] = score

    return all_scores

parser = PDBParser()
for line in open('intpred_list.txt').readlines():
    pdbid = line.rstrip()
    # Open the intpred prediction
    pred_fn = 'raw/'+pdbid
    try:
        intpred_scores = parse_intpred_raw(pdbid, pred_fn)
    except:
        continue

    # Open the pdb file.
    mystruct = parser.get_structure(pdbid, 'orig_pdbs/'+pdbid+'.pdb')
    for model in mystruct:
        for chain in model:
            for residue in chain:
                for atom in residue:
                    key = (chain.get_id(), int(residue.get_id()[1]))
                    if key in intpred_scores:
                        atom.set_bfactor(intpred_scores[key]*100)
                        #atom.set_bfactor(99)
                    else:
                        atom.set_bfactor(-99)
    outio = PDBIO()
    outio.set_structure(mystruct)
    outio.save('intpred_pdbs/'+pdbid+'.pdb')
    print(line)

max_stars_repo_name: NBDsoftware/masif

max_issues_repo_name: NBDsoftware/masif

max_forks_repo_name: NBDsoftware/masif


text:
"""
Main module containing the napari widget for reading from a TissueMAPS interface.
"""
import os
from xml.dom import minidom
import json
from typing import List
from io import BytesIO
import dask
from napari_plugin_engine import napari_hook_implementation
import numpy as np
import requests
import PIL
from napari.types import LayerDataTuple
from qtpy.QtWidgets import QLineEdit  # pylint: disable=E0611
from magicgui.widgets import FunctionGui

from .lazy_array import LazyArray  # pylint: disable=E0401

TILE_SIZE = 256  # hard-coded in TissueMAPS

def authenticate(url, username, password):
    """
    Helper function that returns an authentication token

    :param url: The url of the authentication service
    :param username: The username to be used for authentication
    :param password: The password to be used for authentication
    :return: The access token
    """
    response = requests.post(
        url + '/auth',
        data=json.dumps({'username': username, 'password': password}),
        headers={'content-type': 'application/json'},
    )
    response.raise_for_status()
    data = response.json()
    return data['access_token']

def http_get(url, api_uri, token, **params):
    """
    Helper function to perform an http get, with optional parameters

    :param url: The url of the endpoint
    :param api_url: The URI of the specific api
    :param token: The access token
    :param params: kwargs for optional parameter to the http get
    :return: The full requests response
    """
    response = requests.get(
        url + '/' + api_uri, params=params,
        headers={'Authorization': 'JWT ' + token},
    )
    response.raise_for_status()
    return response

def get_data(url, api_uri, token, **params):
    """
    Helper function to perform an http get to a json enpoint which contains a 'data' field
    Supports optinal parameters for get with optional parameters

    :param url: The url of the endpoint
    :param api_url: The URI of the specific api
    :param token: The access token
    :param params: kwargs for optional parameter to the http get
    :return: The content of the 'data' field in the json response
    """
    response = http_get(url, api_uri, token, **params)
    data = response.json()
    return data['data']

def tissuemaps_interface(url, token, experiment_id, channel_layer_id):
    # pylint: disable=R0914
    """
    Function which reads a channel layer from the TissueMAPs API.
    Given NAPARI_OCTREE==1, it returns a multi-scale (pyramidal) image as a delayed Dask
    array. Otherwise, it returns a high-resolution image as numpy array.
    Note: The later will download the full image into memory!

    :param url: The base url of the TissueMAPs server
    :param token: The authentication token obtained by the TissueMAPs server
    :param experiment_id: The TissueMAPs experiment id
    :param channel_layer_id: The TissueMAPs channel layer id

    :return: A tiled multi-scale image if NAPARI_OCTREE==1,
             otherwise a numpy array with the high-resolution image
    """
    class LazyTiledTMArray(LazyArray):
        # pylint: disable=R0903
        """
        A numpy-like array which lazily loads tiles from a TissueMAPS server.
        """
        def __init__(self, shape, dtype, tile_size, zoom):
            super().__init__(shape, dtype, tile_size)
            self.zoom = zoom

        @dask.delayed
        def read_tile(self, y_tile, x_tile):
            '''
            Reads a tile from a TissuMAPS server
            :param y_tile: the y coordinate of the tile
            :param x_tile: the x coordinate of the tile
            :return: numpy array with the the cooresponding tile updated
            '''
            api_url = (
                'api/experiments/' + str(experiment_id) + '/channel_layers/' +
                str(channel_layer_id) + '/tiles'
            )
            tiles_resp = http_get(url, api_url, token, x=x_tile, y=y_tile, z=self.zoom)
            img = PIL.Image.open(BytesIO(tiles_resp.content))
            data = np.zeros((self.tile_size, self.tile_size))
            data[:img.size[1], :img.size[0]] = np.asarray(img)
            return data

    channel_layers = get_data(
        url,
        'api/experiments/' + str(experiment_id) + '/channel_layers',
        token
    )

    channel_layer = next(
        item for item in channel_layers if item["id"] == channel_layer_id
    )

    image_data = {}
    image_data['image_height'] = channel_layer['image_size']['height']
    image_data['image_width'] = channel_layer['image_size']['width']
    image_data['max_zoom'] = channel_layer['max_zoom']

    pyramid = []
    for zoom in reversed(range(1, image_data['max_zoom'] + 1)):
        if image_data['image_width'] < TILE_SIZE:
            break

        array = LazyTiledTMArray(
            shape=(image_data['image_height'], image_data['image_width']),
            dtype=np.uint8,
            tile_size=TILE_SIZE,
            zoom=zoom
        )

        if 'NAPARI_OCTREE' not in os.environ or os.environ['NAPARI_OCTREE'] != '1':
            # given we don't have a spatial index (e.g. an octree), we can
            # directly return the highest resolution.
            # NOTE: this will download the full image into memory!
            array = np.asarray(array)
            return array

        pyramid.append(array)
        image_data['image_height'] //= 2
        image_data['image_width'] //= 2

    return pyramid

def tissuemaps_connector(path):
    """
    Function which reads an XML specifying TissueMAPS credential and experiment/channel_layer data
    and return a multi-scale (pyramidal) JPEG from TissueMAPs api as delayed Dask
    array.

    :param path: XML file path
    :return: List of LayerData tuple
    """
    xmldoc = minidom.parse(path)
    auth_data = {}
    auth_data['url'] = xmldoc.getElementsByTagName('url')[0].attributes['url'].value
    auth_data['user'] = \
        xmldoc.getElementsByTagName('user')[0].attributes['name'].value
    auth_data['password'] = \
        xmldoc.getElementsByTagName('user')[0].attributes['password'].value
    auth_data['token'] = \
        authenticate(auth_data['url'], auth_data['user'], auth_data['password'])

    query_data = {}
    query_data['experiment_id'] = \
        xmldoc.getElementsByTagName('layerdata')[0].attributes['experiment_id'].value

    query_data['channel_layer_id'] = \
        xmldoc.getElementsByTagName('layerdata')[0].attributes['channel_layer_id'].value

    pyramid = tissuemaps_interface(
        auth_data['url'],
        auth_data['token'],
        query_data['experiment_id'],
        query_data['channel_layer_id']
    )

    return [(pyramid, {})]

# Widget code

class TissueMAPSGetTokenWidget(FunctionGui):
    # pylint: disable=R0901
    # Disabled check for number of ancestors, since this class will have 17 ancestors,
    # way more than the 5 suggested by pylint
    """
    Inner widget to handle connection to a given TissueMAPS server.
    This widget stores username, password and access token
    """
    def __init__(self, value=None, name="tm_connector", **kwargs):
        # pylint: disable=W0613
        if value is None:
            value = ("", "", "")
        url, username, password = value
        super().__init__(
            TissueMAPSGetTokenWidget.apply,
            call_button=False,
            layout='vertical',
            param_options={
                "url": {"widget_type": "LineEdit"},
                "username": {"widget_type": "LineEdit"},
                "password": {"widget_type": "LineEdit"},
                "add_button": {
                    "widget_type": "PushButton", "text": "Connect",
                }
            },
            name=name
        )

        self.password.native.setEchoMode(QLineEdit.Password)
        self.url.value = url
        self.username.value = username
        self.password.value = password
        self.token = ""

        self.native.layout().setContentsMargins(0, 0, 0, 0)

        @self.add_button.changed.connect
        def on_press_import_button(event):
            # pylint: disable=W0613
            self.token = authenticate(self.url.value, self.username.value, self.password.value)

    def __setitem__(self, key, value):
        """Prevent assignment by index."""
        raise NotImplementedError("magicgui.Container does not support item setting.")

    @staticmethod
    def apply(url="", username="", password="", add_button=True):
        # pylint: disable=W0613
        """
        Dummy function to respect the FunctionGui logic. Not used since
        call_button is False in this widget
        """

    @property
    def value(self):
        """
        Associates the value field of the TissueMAPSGetTokenWidget to the access token
        """
        return self.token

class TissueMAPSConnectionWidget(FunctionGui):
    # pylint: disable=R0901,R0903
    # Disabled check for number of ancestors, since this class will have 17 ancestors,
    # way more than the 5 suggested by pylint
    """
    Main widget to manage TissueMAPS data. It contains the TissueMAPSGetTokenWidget
    """
    def __init__(self, value=None, **kwargs):
        # pylint: disable=W0613
        super().__init__(
            self.apply,
            call_button="Load Data",
            layout="vertical",
            param_options={
                "token": {"widget_type": TissueMAPSGetTokenWidget, "name": "tm_connector"},
                "experiment_name": {"choices": [""]},
                "channel_name": {"choices": [""]},
            },
        )

        def get_experiment_names(*args):
            # pylint: disable=W0613
            if len(self.experiments) > 0:
                return [experiment['name'] for experiment in self.experiments]
            return []

        def get_channel_names(*args):
            # pylint: disable=W0613
            channel_names = []
            if len(self.channels) > 0:
                channel_names = [channel['name'] for channel in self.channels]
            if len(channel_names) > 1:
                channel_names.insert(0, "-- All --")
            return channel_names

        @self.tm_connector.changed.connect
        def update_experiments(event):
            # pylint: disable=W0613
            if self.tm_connector.token != "":
                resp = get_data(
                    self.tm_connector.url.value,
                    "/api/experiments",
                    self.tm_connector.token
                )
                self.experiments = resp
                self.experiment_name.choices = []
                self.experiment_name.reset_choices()

        @self.experiment_name.changed.connect
        def update_channels(event):
            # pylint: disable=W0613
            exp_id = [
                exp["id"] for exp in self.experiments if exp["name"] == self.experiment_name.value
            ][0]
            resp = get_data(
                self.tm_connector.url.value,
                'api/experiments/' + str(exp_id) + '/channels',
                self.tm_connector.token
            )
            self.channels = resp
            self.channel_name.choices = []
            self.channel_name.reset_choices()

        self.experiments = []
        self.channels = []
        self.experiment_name._default_choices = get_experiment_names
        self.channel_name._default_choices = get_channel_names

        self.native.layout().addStretch()

    def __setitem__(self, key, value):
        """Prevent assignment by index."""
        raise NotImplementedError("magicgui.Container does not support item setting.")
        # pylint: disable=C0301

    def apply(self, token=("", "", ""),
              experiment_name="",
              channel_name="") -> List[LayerDataTuple]:
        # pylint: disable=W0613
        """
        Function executed when the "Load Data" button is pressed.
        It calls the tissuemaps tissuemaps_interface and returns a napari
        Image layer

        :param token: The access token for querying TissueMAPS
        :param experiment_name: The name of a TissueMAPS experiment
        :param channel_name: The channel name of a TissueMAPS experiment
        :return: napari_layers.Image object, with access_token stored as metadata
        """
        exp_id = [exp["id"] for exp in self.experiments if exp["name"] == experiment_name][0]
        if channel_name == '-- All --':
            multi_layer = []
            for channel in self.channels:
                assert len(channel["layers"]) == 1
                for layer in channel["layers"]:
                    pyramid = tissuemaps_interface(
                        self.tm_connector.url.value,
                        token,
                        exp_id,
                        str(layer["id"])
                    )
                    res = (
                        pyramid, {
                            'name': channel['name'],
                            'metadata': {'token': self.tm_connector.token},
                            'opacity': 1.0 / len(self.channels),
                            'blending': 'additive'
                        },
                        'image'
                    )
                    multi_layer.append(res)
            return list(
                reversed(
                    sorted(multi_layer, key=lambda layer_data: layer_data[1]["name"])
                )
            )

        channels = [ch for ch in self.channels if ch["name"] == channel_name]
        assert len(channels) == 1
        layer_ids = [layer["id"] for layer in channels[0]["layers"]]
        assert len(layer_ids) == 1
        pyramid = tissuemaps_interface(
            self.tm_connector.url.value,
            token,
            exp_id,
            layer_ids[0]
        )
        return [(
            pyramid,
            {
                'name': channel_name,
                'metadata': {'token': self.tm_connector.token},
                'blending': 'additive'
            },
            'image'
        )]

@napari_hook_implementation
def napari_get_reader(path):
    """
    Napari plugin that returns a reader interface for TissueMAPs .

    .. note::
       This hook does not support a list of paths

    :param path:  The path of the image
    :return: The tissuemaps_interface function or None
    """
    if isinstance(path, str) and path.endswith(".xmld"):
        return tissuemaps_connector
    return None

@napari_hook_implementation
def napari_experimental_provide_dock_widget():
    """
    Napari plugin that returns a Magicui widget
    :return: The TissueMAPS connection widget
    """
    return TissueMAPSConnectionWidget


text:
import numpy as np
from scipy.ndimage import imread
from scipy.misc import imsave
import sys

img_fn = sys.argv[1]
out_fn = img_fn + '_alpha.png'
print('Input: %s\nOutput: %s' % (img_fn, out_fn))
img = imread(img_fn)
img[:, :, 3] = 255
imsave(out_fn, img)

max_stars_repo_name: chili-epfl/qml-ar

max_issues_repo_name: chili-epfl/qml-ar

max_forks_repo_name: chili-epfl/qml-ar

reason: import numpy,from scipy


text:
#!/usr/bin/env python
"""
create Latex table with comparison of polarization integrals for
cases 1.1, 2.1 and 2.2 with the
original analytical (wrong), the corrected analytical and the numerical integration methods.
"""
import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)

import numpy as np
import numpy.linalg as la

# fast C++ implementation
from polarization_integrals import PolarizationIntegral

# numerical integral using Becke's multicenter grids
from polarization_ints_numerical import polarization_integral as polarization_integral_numerical
# increase resolution of integration grids
from becke import settings
settings.radial_grid_factor = 3  # increase number of radial points by factor 3
settings.lebedev_order = 23      # angular Lebedev grid of order 23

from polarization_ints_reference import polarization_integral as polarization_integral_reference

# exponents
beta_i = 1.0
beta_j = 1.0

# cutoff function
alpha = 50.0
q = 2

# centers of GTOs
xi,yi = 0,0
xj,yj = 0,0
# The z-coordiantes of GTOs are varied to get different cases (b > 0 and b == 0)

# powers of Cartesian coordinates
nxi,nyi,nzi = 0,0,0
nxj,nyj,nzj = 0,0,0

# powers in operator
mx,my = 0,0
# mz is varied to get different cases (1, 2.1 or 2.2)

def write_row(k=3, mz=0, zi=0.5, zj=0.5, nzj=0):
    li = nxi+nyi+nzi
    lj = nxj+nyj+nzj

    # prepare for integrals between shells with angular momenta li and lj
    pol = PolarizationIntegral(xi,yi,zi, li, beta_i,
                               xj,yj,zj, lj, beta_j,
                               k, mx,my,mz,
                               alpha, q)
    # fast C++ implementation
    pint = pol.compute_pair(nxi,nyi,nzi,
                            nxj,nyj,nzj)

    # numerical integrals
    # The numerical integration fails if the centers coincide exactly, so we add a small number
    eps = 1.0e-6
    if (abs(zi) < 10*eps) and (abs(zj) < 10*eps):
        zi -= eps
        zj += eps
        zi_ = zi - eps
        zj_ = zj + eps
    else:
        zi_ = zi
        zj_ = zj

    pint_numerical = polarization_integral_numerical(xi,yi,zi_, nxi,nyi,nzi, beta_i,
                                                     xj,yj,zj_, nxj,nyj,nzj, beta_j,
                                                     k, mx,my,mz,
                                                     alpha, q)

    try:
        # Xiao's reference integrals using the original definition of H(1,a) in Schwerdtfeger's article
        pint_reference = polarization_integral_reference(xi,yi,zi, nxi,nyi,nzi, beta_i,
                                                         xj,yj,zj, nxj,nyj,nzj, beta_j,
                                                         k, mx,my,mz,
                                                         alpha, q,
                                                         original=True)
    except ZeroDivisionError:
        pint_reference = np.nan

    #### determine which sub case we are dealing with

    # 1, 2.1 or 2.2 ?
    if (k % 2 == 0):
        case = "     1  "
    else:
        j = k//2
        case = "     2."
        s_max = li+lj+mx+my+mz
        if (s_max - j >= 0):
            case += "1"
        else:
            case += "2"

    # b > 0 or b = 0 ?

    ri = np.array([xi,yi,zi])
    rj = np.array([xj,yj,zj])
    b = la.norm(beta_i * ri + beta_j * rj)
    if (abs(b) < 1.0e-10):
        case += " (b = 0) "
    else:
        case += "         "

    print(f"{case} &  {k} & {mz}     & ${zi}$  & ${zj}$  &  ${pint_reference:+12.7f}$ &  ${pint:+12.7f}$   & ${pint_numerical:+12.7f}$     \\\\")

# write header
print(f"Case              &  k & $m_3$ & $z_1$  & $z_2$  &       original  &       corrected   &      numerical     \\\\ \\midrule")

# Case 1  (b > 0)
write_row(k=4, mz=0, zi=-1.0, zj=0.5)
## Case 1  (b == 0)
#write_row(k=4, mz=0, zi=-1.0, zj=1.0)
# Case 1  (b == 0, all centers coincide)
#write_row(k=4, mz=0, zi=0.0, zj=0.0)

# Case 2.1  (b > 0)
write_row(k=3, mz=1, zi=-1.0, zj=0.5)
## Case 2.1  (b == 0)
#write_row(k=3, mz=1, zi=-1.0, zj=1.0)
# Case 2.1  (b == 0, all centers coincide)
#write_row(k=3, mz=1, zi=0.0, zj=0.0)

# Case 2.2  (b > 0)
write_row(k=3, mz=0, zi=-1.0, zj=0.5)
## Case 2.2  (b == 0)
#write_row(k=3, mz=0, zi=-1.0, zj=1.0)
# Case 2.2  (b == 0, all centers coincide)
write_row(k=3, mz=0, zi=0.0, zj=0.0)


text:
import pytest
from astropy.io import fits as pyfits
from astropy.utils.data import get_pkg_data_filename
from ..kepbls import kepbls
from ..kepio import delete

# fake light curve with transits of period = 2.02 days
# and duration 0.18 days
fake_lc = get_pkg_data_filename("data/golden-lc.fits")

def test_kepbls():
    kepbls(fake_lc, "kepbls.fits", datacol='SAP_FLUX', errcol='SAP_FLUX_ERR',
           minper=1, maxper=3, mindur=0, maxdur=.5, nsearch=100, nbins=100,
           overwrite=True)
    f = pyfits.open("kepbls.fits")
    assert abs(f[3].header['PERIOD'] - 2.02) < 0.1
    assert abs(f[3].header['TRANSDUR'] - 0.18) < 0.01
    f.close()
    delete("kepbls.fits", "log_kepextract.txt", False)

max_stars_repo_name: ecalifornica/pyke

max_issues_repo_name: ecalifornica/pyke

max_forks_repo_name: mirca/PyKE

reason: from astropy


text:
import unittest
import amorf.metrics as metrics
import numpy as np
import torch as torch

class TestError(unittest.TestCase):

    def test_average_correlation_coefficient(self):
        labels = np.array([[2, -6], [3, 5], [4, 4], [5, 3], [-6, 2]])
        predicted_labels = np.array([[3, 7], [2, -5], [4, 4], [-5, 2], [7, 3]])
        self.assertAlmostEqual(metrics.average_correlation_coefficient(
            predicted_labels, labels), -0.725679, places=5)
    def test_tensor_average_correlation_coefficient(self):
        labels = torch.from_numpy(
            np.array([[2, -6], [3, 5], [4, 4], [5, 3], [-6, 2]])).float()
        predicted_labels = torch.from_numpy(
            np.array([[3, 7], [2, -5], [4, 4], [-5, 2], [7, 3]])).float()
        self.assertAlmostEqual(metrics.average_correlation_coefficient(
            predicted_labels, labels).item(), -0.725679, places=5)

    def test_average_relative_error(self):
        labels = np.array([[2, -6], [3, 5], [4, 4], [5, 3], [-6, 2]])
        predicted_labels = np.array([[3, 7], [2, -5], [4, 4], [-5, 2], [7, 3]])
        self.assertAlmostEqual(metrics.average_relative_error(
            predicted_labels, labels), 0.133333, places=5)

    def test_tensor_average_relative_error(self):
        labels = torch.from_numpy(
            np.array([[2, -6], [3, 5], [4, 4], [5, 3], [-6, 2]])).float()
        predicted_labels = torch.from_numpy(
            np.array([[3, 7], [2, -5], [4, 4], [-5, 2], [7, 3]])).float()
        self.assertAlmostEqual(metrics.average_relative_error(
            predicted_labels, labels).item(), 0.133333, places=5)

    def test_average_relative_root_mean_squared_error(self):
        labels = np.array(
            [[2, -6], [3, 5], [4, 4], [5, 3], [-6, 2]]).astype(np.float)
        predicted_labels = np.array(
            [[3, 7], [2, -5], [4, 4], [-5, 2], [7, 3]]).astype(np.float)
        self.assertAlmostEqual(metrics.average_relative_root_mean_squared_error(
            predicted_labels, labels), 1.8735961929670222, places=5)

    def test_tensor_average_relative_root_mean_squared_error(self):
        labels = torch.from_numpy(
            np.array([[2, -6], [3, 5], [4, 4], [5, 3], [-6, 2]])).float()
        predicted_labels = torch.from_numpy(
            np.array([[3, 7], [2, -5], [4, 4], [-5, 2], [7, 3]])).float()
        self.assertAlmostEqual(metrics.average_relative_root_mean_squared_error(
            predicted_labels, labels), 1.8735961929670222, places=5)

    def test_tensor_mean_squared_error(self):
        labels = torch.from_numpy(
            np.array([[2, -6], [3, 5], [4, 4], [5, 3], [-6, 2]])).float()
        predicted_labels = torch.from_numpy(
            np.array([[3, 7], [2, -5], [4, 4], [-5, 2], [7, 3]])).float()
        self.assertEqual(metrics.mean_squared_error(
            predicted_labels, labels), 108.4)

    def test_mean_squared_error(self):
        labels = np.array([[2, -6], [3, 5], [4, 4], [5, 3], [-6, 2]])
        predicted_labels = np.array([[3, 7], [2, -5], [4, 4], [-5, 2], [7, 3]])
        self.assertEqual(metrics.mean_squared_error(
            predicted_labels, labels), 108.4)

    def test_tensor_average_root_mean_squared_error(self):
        labels = torch.from_numpy(
            np.array([[2, -6], [3, 5], [4, 4], [5, 3], [-6, 2]])).float()
        predicted_labels = torch.from_numpy(
            np.array([[3, 7], [2, -5], [4, 4], [-5, 2], [7, 3]])).float()
        self.assertAlmostEqual(metrics.average_root_mean_squared_error(
            predicted_labels, labels).item(), 7.36206, places=3)

    def test_average_root_mean_squared_error(self):
        labels = np.array([[2, -6], [3, 5], [4, 4], [5, 3], [-6, 2]])
        predicted_labels = np.array([[3, 7], [2, -5], [4, 4], [-5, 2], [7, 3]])
        self.assertAlmostEqual(metrics.average_root_mean_squared_error(
            predicted_labels, labels), 7.36206, places=3)

if __name__ == '__main__':
    unittest.main()

max_stars_repo_name: DSAAR/amorf

max_issues_repo_name: DSAAR/amorf

max_forks_repo_name: DSAAR/amorf


text:
# floor
# The tool floor returns the floor of the input element-wise.
# The floor of x is the largest integer i where i <= x.

# import numpy
# my_array = numpy.array([1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7, 8.8, 9.9])
# print numpy.floor(my_array)         #[ 1.  2.  3.  4.  5.  6.  7.  8.  9.]

# ceil
# The tool ceil returns the ceiling of the input element-wise.
# The ceiling of x is the smallest integer i where i >= x.

# import numpy
# my_array = numpy.array([1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7, 8.8, 9.9])
# print numpy.ceil(my_array)          #[  2.   3.   4.   5.   6.   7.   8.   9.  10.]

# rint
# The rint tool rounds to the nearest integer of input element-wise.

# import numpy
# my_array = numpy.array([1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7, 8.8, 9.9])
# print numpy.rint(my_array)          #[  1.   2.   3.   4.   6.   7.   8.   9.  10.]

# Task
# You are given a 1-D array, A.
# Your task is to print the floor, ceil and rint of all the elements of A.

# Input Format
# A single line of input containing the space separated elements of array A.

# Output Format
# On the first line, print the floor of A.
# On the second line, print the ceil of A.
# On the third line, print the rint of A.

# Sample Input
# 1.1 2.2 3.3 4.4 5.5 6.6 7.7 8.8 9.9

# Sample Output
# [ 1.  2.  3.  4.  5.  6.  7.  8.  9.]
# [  2.   3.   4.   5.   6.   7.   8.   9.  10.]
# [  1.   2.   3.   4.   6.   7.   8.   9.  10.]

import numpy as np

# Important part of the solution
np.set_printoptions(sign = ' ') # Turns out numpy has print-formatting settings

l = np.array(list(map(float, input().rstrip().split())))
print(np.floor(l), np.ceil(l), np.rint(l), sep = '\n')


text:
"""
Implementation of the class `IBAMRSimulation`.
"""

import os

import numpy

from ..simulation import Simulation
from ..force import Force

class IBAMRSimulation(Simulation):
  """
  Contains info about a IBAMR simulation.
  Inherits from class Simulation.
  """

  def __init__(self, description=None, directory=os.getcwd(), **kwargs):
    """
    Initializes by calling the parent constructor.

    Parameters
    ----------
    description: string, optional
      Description of the simulation;
      default: None.
    directory: string, optional
      Directory of the simulation;
      default: <current working directory>.
    """
    super(IBAMRSimulation, self).__init__(software='ibamr',
                                          description=description,
                                          directory=directory,
                                          **kwargs)

  def read_forces(self, file_path=None, labels=None):
    """
    Reads forces from files.

    Parameters
    ----------
    file_path: string, optional
      Path of the file containing the forces;
      default: None.
    labels: list of strings, optional
      Label of each force to read;
      default: None.
    """
    if not file_path:
      file_path = os.path.join(self.directory,
                               'dataIB',
                               'ib_Drag_force_struct_no_0')
    print('[info] reading forces from {} ...'.format(file_path)),
    with open(file_path, 'r') as infile:
      times, force_x, force_y = numpy.loadtxt(infile,
                                              dtype=float,
                                              usecols=(0, 4, 5),
                                              unpack=True)
    self.forces = []
    self.forces.append(Force(times, force_x, '$F_x$'))
    self.forces.append(Force(times, force_y, '$F_y$'))
    print('done')

  def write_visit_summary_files(self, time_steps):
    """
    Writes summary files for Visit with list of sub-directories to look into.

    Parameters
    ----------
    time_steps: 3-tuple of integers
      Staring and and ending time_steps followed by the time-step increment.
    """
    print('[info] writing summary files for VisIt ...'),
    time_steps = numpy.arange(time_steps[0], time_steps[1] + 1, time_steps[2])
    # list of SAMRAI files to VisIt
    dumps_visit_list = [os.path.join('visit_dump.{0:05}'.format(time_step),
                                     'summary.samrai')
                        for time_step in time_steps]
    dumps_visit = numpy.array(dumps_visit_list)
    # list of SILO file to VisIt
    lag_visit_list = [os.path.join('lag_data.cycle_{:06d}'.format(time_step),
                                   'lag_data.cycle_{:06d}.summary.silo'
                                   ''.format(time_step))
                      for time_step in time_steps]
    lag_data_visit = numpy.array(lag_visit_list)
    # write files
    with open(os.path.join(self.directory, 'dumps.visit'), 'w') as outfile:
      numpy.savetxt(outfile, dumps_visit, fmt='%s')
    with open(os.path.join(self.directory, 'lag_data.visit'), 'w') as outfile:
      numpy.savetxt(outfile, lag_data_visit, fmt='%s')
    print('done')

  def plot_field_contours_visit(self, field_name,
                                field_range,
                                body=None,
                                solution_folder='numericalSolution',
                                states=(0, 20000, 1),
                                view=(-2.0, -2.0, 2.0, 2.0),
                                width=800):
    """
    Plots the contour of a given field using VisIt.

    Parameters
    ----------
    field_name: string
      Name of field to plot;
      choices: vorticity, pressure, velocity-magnitude, x-velocity, y-velocity.
    field_range: 2-tuple of floats
      Range of the field to plot (min, max).
    body: string, optional
      Name of the immersed body;
      default: None.
    solution_folder: string, optional
      Relative path of the folder containing the numerical solution;
      default: 'numericalSolution'.
    states: 3-tuple of integers, optional
      Limits of index of the states to plot followed by the increment;
      default: (0, 20000, 1).
    view: 4-tuple of floats, optional
      Bottom-left and top-right coordinates of the view to display;
      default: (-2.0, -2.0, 2.0, 2.0).
    width: integer, optional
      Width (in pixels) of the figure;
      default: 800.
    """
    args = {}
    args['--directory'] = self.directory
    args['--field'] = field_name
    args['--range'] = '{} {}'.format(*field_range)
    if body:
      args['--body'] = body
    args['--solution-folder'] = solution_folder
    args['--states'] = '{} {} {}'.format(*states)
    args['--view'] = '{} {} {} {}'.format(*view)
    args['--width'] = str(width)
    script = os.path.join(os.environ['SNAKE'],
                          'snake',
                          'ibamr',
                          'plotField2dVisIt.py')
    arguments = ' '.join([key + ' ' + value
                          for key, value in args.items()])
    os.system('visit -nowin -cli -s {} {}'.format(script, arguments))

  def compute_mean_number_cells_visit(self,
                                      solution_folder='numericalSolution',
                                      states=(0, 20000, 1),
                                      time_limits=(0.0, float('inf'))):
    """
    Computes the number of a cells, on average, in the AMR mesh at requested
    states.

    Parameters
    ----------
    solution_folder: string, optional
      Relative path of the folder containing the numerical solution;
      default: 'numericalSolution'.
    states: 3-tuple of integers, optional
      Starting and ending states followed by the increment;
      default: (0, 20000, 1).
    time_limits: 2-tuple of floats, optional
      Time-limits within which the mean value is calculated;
      default: (0.0, inf).
    """
    args = {}
    args['--directory'] = self.directory
    args['--solution-folder'] = solution_folder
    args['--states'] = '{} {} {}'.format(*states)
    args['--time-limits'] = '{} {}'.format(*time_limits)
    script = os.path.join(os.environ['SNAKE'],
                          'snake',
                          'ibamr',
                          'getNumberCellsVisIt.py')
    arguments = ' '.join([key + ' ' + value
                          for key, value in args.items()])
    os.system('visit -nowin -cli -s {} {}'.format(script, arguments))

max_stars_repo_name: mesnardo/cuIBM

max_issues_repo_name: mesnardo/cuIBM

max_forks_repo_name: mesnardo/cuIBM


text:
import torch
import numpy as np
import argparse
import time
import Utils.util as util
from engine import trainer
import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt
import pickle
import os
import ipdb

# python train.py --gcn_bool --adjtype doubletransition --addaptadj  --randomadj --num_nodes 207 --seq_length 12 --save ./garage/metr
# python train.py --gcn_bool --adjtype doubletransition --addaptadj  --randomadj --num_nodes 80 --data syn --blocks 5
# python train.py --gcn_bool --adjtype doubletransition --addaptadj  --randomadj --data CRASH

parser = argparse.ArgumentParser()
parser.add_argument('--device',type=str,default='cuda:0',help='')
parser.add_argument('--data',type=str,default='data/METR-LA',help='data path')
parser.add_argument('--adjdata',type=str,default='data/sensor_graph/adj_mx.pkl',help='adj data path')
parser.add_argument('--adjtype',type=str,default='doubletransition',help='adj type')
parser.add_argument('--gcn_bool',action='store_true',help='whether to add graph convolution layer')
parser.add_argument('--aptonly',action='store_true',help='whether only adaptive adj')
parser.add_argument('--addaptadj',action='store_true',help='whether add adaptive adj')
parser.add_argument('--randomadj',action='store_true',help='whether random initialize adaptive adj')
parser.add_argument('--seq_length',type=int,default=48,help='')
parser.add_argument('--nhid',type=int,default=32,help='')
parser.add_argument('--in_dim',type=int,default=2,help='inputs dimension')
parser.add_argument('--num_nodes',type=int,default=80,help='number of nodes')
parser.add_argument('--layers',type=int,default=2,help='number of layers per gwnet block')
parser.add_argument('--blocks',type=int,default=4,help='number of blocks in gwnet model')

parser.add_argument('--batch_size',type=int,default=32,help='batch size')
parser.add_argument('--learning_rate',type=float,default=0.001,help='learning rate')
parser.add_argument('--dropout',type=float,default=0.3,help='dropout rate')
parser.add_argument('--weight_decay',type=float,default=0.0001,help='weight decay rate')
parser.add_argument('--epochs',type=int,default=100,help='')
parser.add_argument('--print_every',type=int,default=50,help='')
# parser.add_argument('--seed',type=int,default=0,help='random seed')
parser.add_argument('--save',type=str,default='./garage/syn',help='save path')
parser.add_argument('--expid',type=int,default=1,help='experiment id')

args = parser.parse_args()

np.random.seed(0)
torch.manual_seed(999)
if torch.cuda.is_available():
    torch.cuda.manual_seed_all(999)
    torch.cuda.empty_cache()

def main(model_name=None, syn_file='syn_diffG.pkl'): # directly loading trained model/ generated syn data
    #set seed
    #torch.manual_seed(args.seed)
    #np.random.seed(args.seed)
    #load data
    same_G = False
    device = torch.device(args.device)

    if args.data == 'syn':
        if  os.path.isfile(syn_file):
            with open(syn_file, 'rb') as handle:
                pkl_data = pickle.load(handle)
            nTrain, nValid, nTest, num_timestep = pkl_data['nTrain'], pkl_data['nValid'],\
                                                  pkl_data['nTest'], pkl_data['num_timestep']
            dataloader, adj_mx, F_t, G = pkl_data['dataloader'], pkl_data['adj_mx'],\
                                         pkl_data['F_t'], pkl_data['G']
            print('synthetic data loaded')
        else:
            nTrain = 80 # Number of training samples
            nValid = int(0.25 * nTrain) # Number of validation samples
            nTest = int(0.05 * nTrain) # Number of testing samples
            num_timestep = 1000 # 1000
            dataloader, adj_mx, F_t, G = util.load_dataset_syn(args.adjtype, args.num_nodes,
                                                               nTrain, nValid, nTest, num_timestep,
                                                               args.seq_length, args.batch_size,
                                                               args.batch_size, args.batch_size,
                                                               same_G=same_G)
            # pkl_data = {'nTrain': nTrain, 'nValid': nValid, 'nTest': nTest,
            #             'num_timestep': num_timestep, 'dataloader': dataloader,
            #             'adj_mx': adj_mx, 'F_t': F_t, 'G':G}
            # with open(syn_file, 'wb') as handle:
            #     pickle.dump(pkl_data, handle, protocol=pickle.HIGHEST_PROTOCOL)

    elif args.data == 'CRASH':
        util.load_dataset_CRASH(args.adjtype, args.batch_size, args.batch_size, args.batch_size)

    else:
        sensor_ids, sensor_id_to_ind, adj_mx = util.load_adj(args.adjdata,args.adjtype)
        dataloader = util.load_dataset_metr(args.data, args.batch_size, args.batch_size,
                                            args.batch_size)

    if args.data == 'syn' and not same_G: # different graph structure for each sample
        assert len(adj_mx) == nTrain + nValid + nTest

        # separate adj matrices into train-val-test samples
        adj_train = [[],[]]
        for a in adj_mx[:nTrain]:
            adj_train[0].append(a[0])
            adj_train[1].append(a[1])
        adj_train = [np.stack(np.asarray(i)) for i in adj_train]

        adj_val = [[],[]]
        for a in adj_mx[nTrain:-nTest]:
            adj_val[0].append(a[0])
            adj_val[1].append(a[1])
        adj_val = [np.stack(np.asarray(i)) for i in adj_val]

        adj_test = [[],[]]
        for a in adj_mx[-nTest:]:
            adj_test[0].append(a[0])
            adj_test[1].append(a[1])
        adj_test = [np.stack(np.asarray(i)) for i in adj_test]

        scaler = dataloader['scaler']
        print(args)
        supports = {}
        supports['train'] = [torch.tensor(i).to(device) for i in adj_train]
        supports['val'] = [torch.tensor(i).to(device) for i in adj_val]
        supports['test'] = [torch.tensor(i).to(device) for i in adj_test]

        adjinit = {}
        if args.randomadj:
            adjinit['train'] = adjinit['val'] = adjinit['test'] = None
        else:
            adjinit['train'] = supports['train'][0]
            adjinit['val'] = supports['val'][0]
            adjinit['test'] = supports['test'][0]

        if args.aptonly:
            supports['train'] = supports['val'] = supports['test'] = None

        engine = trainer(scaler, args.in_dim, args.seq_length, args.num_nodes, args.nhid, args.dropout,
                         args.learning_rate, args.weight_decay, device, supports, args.gcn_bool,
                         args.addaptadj, adjinit, args.blocks, args.layers)

        if model_name is None:
            print("start training...",flush=True)

            his_loss =[]
            val_time = []
            train_time = []
            for i in range(1,args.epochs+1):
                #if i % 10 == 0:
                    #lr = max(0.000002,args.learning_rate * (0.1 ** (i // 10)))
                    #for g in engine.optimizer.param_groups:
                        #g['lr'] = lr
                train_loss = []
                train_mape = []
                train_rmse = []
                t1 = time.time()
                dataloader['train_loader'].shuffle()
                engine.set_state('train')

                for iter, (x, y, adj_idx) in enumerate(dataloader['train_loader'].get_iterator()):
                    trainx = torch.Tensor(x).to(device) # torch.Size([64, 15, 80, 2])
                    trainx= trainx.transpose(1, 3) # torch.Size([64, 2, 80, 15])
                    trainy = torch.Tensor(y).to(device)
                    trainy = trainy.transpose(1, 3)

                    metrics = engine.train_syn(trainx, trainy, F_t, G['train'], adj_idx)
                    train_loss.append(metrics[0])
                    train_mape.append(metrics[1])
                    train_rmse.append(metrics[2])
                    if iter % args.print_every == 0 :
                        log = 'Iter: {:03d}, Train Loss: {:.4f}, Train MAPE: {:.4f}, Train RMSE: {:.4f}'
                        print(log.format(iter, train_loss[-1], train_mape[-1], train_rmse[-1]),flush=True)

                t2 = time.time()
                train_time.append(t2-t1)

                #validation
                valid_loss = []
                valid_mape = []
                valid_rmse = []

                s1 = time.time()
                engine.set_state('val')
                for iter, (x, y, adj_idx) in enumerate(dataloader['val_loader'].get_iterator()):
                    testx = torch.Tensor(x).to(device)
                    testx = testx.transpose(1, 3)
                    testy = torch.Tensor(y).to(device)
                    testy = testy.transpose(1, 3)
                    # [64, 2, 80, 15]
                    metrics = engine.eval_syn(testx, testy, F_t, G['val'], adj_idx)
                    valid_loss.append(metrics[0])
                    valid_mape.append(metrics[1])
                    valid_rmse.append(metrics[2])
                s2 = time.time()
                log = 'Epoch: {:03d}, Inference Time: {:.4f} secs'
                print(log.format(i,(s2-s1)))
                val_time.append(s2-s1)
                mtrain_loss = np.mean(train_loss)
                mtrain_mape = np.mean(train_mape)
                mtrain_rmse = np.mean(train_rmse)

                mvalid_loss = np.mean(valid_loss)
                mvalid_mape = np.mean(valid_mape)
                mvalid_rmse = np.mean(valid_rmse)
                his_loss.append(mvalid_loss)

                log = 'Epoch: {:03d}, Train Loss: {:.4f}, Train MAPE: {:.4f}, Train RMSE: {:.4f}, Valid Loss: {:.4f}, Valid MAPE: {:.4f}, Valid RMSE: {:.4f}, Training Time: {:.4f}/epoch'
                print(log.format(i, mtrain_loss, mtrain_mape, mtrain_rmse, mvalid_loss, mvalid_mape, mvalid_rmse, (t2 - t1)),flush=True)
                torch.save(engine.model.state_dict(), args.save+"_epoch_"+str(i)+"_"+str(round(mvalid_loss,2))+".pth")
            print("Average Training Time: {:.4f} secs/epoch".format(np.mean(train_time)))
            print("Average Inference Time: {:.4f} secs".format(np.mean(val_time)))

    else:
        scaler = dataloader['scaler']
        supports = [torch.tensor(i).to(device) for i in adj_mx]
        print(args)

        if args.randomadj:
            adjinit = None
        else:
            adjinit = supports[0]

        if args.aptonly:
            supports = None

        engine = trainer(scaler, args.in_dim, args.seq_length, args.num_nodes, args.nhid, args.dropout,
                             args.learning_rate, args.weight_decay, device, supports, args.gcn_bool,
                             args.addaptadj, adjinit, args.blocks, args.layers)

        if model_name is None:
            print("start training...",flush=True)
            his_loss =[]
            val_time = []
            train_time = []
            for i in range(1,args.epochs+1):
                #if i % 10 == 0:
                    #lr = max(0.000002,args.learning_rate * (0.1 ** (i // 10)))
                    #for g in engine.optimizer.param_groups:
                        #g['lr'] = lr
                train_loss = []
                train_mape = []
                train_rmse = []
                t1 = time.time()
                dataloader['train_loader'].shuffle()
                for iter, (x, y) in enumerate(dataloader['train_loader'].get_iterator()):
                    trainx = torch.Tensor(x).to(device) # torch.Size([64, 12, 207, 2])
                    trainx= trainx.transpose(1, 3) # torch.Size([64, 2, 207, 12])
                    trainy = torch.Tensor(y).to(device)
                    trainy = trainy.transpose(1, 3)
                    if args.data == 'syn':
                        metrics = engine.train_syn(trainx, trainy, F_t, G)
                    else:
                        metrics = engine.train(trainx, trainy[:,0,:,:])
                    train_loss.append(metrics[0])
                    train_mape.append(metrics[1])
                    train_rmse.append(metrics[2])
                    if iter % args.print_every == 0 :
                        log = 'Iter: {:03d}, Train Loss: {:.4f}, Train MAPE: {:.4f}, Train RMSE: {:.4f}'
                        print(log.format(iter, train_loss[-1], train_mape[-1], train_rmse[-1]),flush=True)

                t2 = time.time()
                train_time.append(t2-t1)
                #validation
                valid_loss = []
                valid_mape = []
                valid_rmse = []

                s1 = time.time()
                for iter, (x, y) in enumerate(dataloader['val_loader'].get_iterator()):
                    testx = torch.Tensor(x).to(device)
                    testx = testx.transpose(1, 3)
                    testy = torch.Tensor(y).to(device)
                    testy = testy.transpose(1, 3)
                    if args.data == 'syn':
                        metrics = engine.eval_syn(testx, testy, F_t, G)
                    else:
                        metrics = engine.eval(testx, testy[:,0,:,:])
                    valid_loss.append(metrics[0])
                    valid_mape.append(metrics[1])
                    valid_rmse.append(metrics[2])

                s2 = time.time()
                log = 'Epoch: {:03d}, Inference Time: {:.4f} secs'
                print(log.format(i,(s2-s1)))
                val_time.append(s2-s1)
                mtrain_loss = np.mean(train_loss)
                mtrain_mape = np.mean(train_mape)
                mtrain_rmse = np.mean(train_rmse)

                mvalid_loss = np.mean(valid_loss)
                mvalid_mape = np.mean(valid_mape)
                mvalid_rmse = np.mean(valid_rmse)
                his_loss.append(mvalid_loss)

                log = 'Epoch: {:03d}, Train Loss: {:.4f}, Train MAPE: {:.4f}, Train RMSE: {:.4f}, Valid Loss: {:.4f}, Valid MAPE: {:.4f}, Valid RMSE: {:.4f}, Training Time: {:.4f}/epoch'
                print(log.format(i, mtrain_loss, mtrain_mape, mtrain_rmse, mvalid_loss, mvalid_mape, mvalid_rmse, (t2 - t1)),flush=True)
                torch.save(engine.model.state_dict(), args.save+"_epoch_"+str(i)+"_"+str(round(mvalid_loss,2))+".pth")
            print("Average Training Time: {:.4f} secs/epoch".format(np.mean(train_time)))
            print("Average Inference Time: {:.4f} secs".format(np.mean(val_time)))

    ################################ TESTING ################################
    if model_name is None:
        bestid = np.argmin(his_loss)
        print(bestid)

        print("Training finished")
        print("The valid loss on best model is", str(round(his_loss[bestid],4)))

        engine.model.load_state_dict(torch.load(args.save+"_epoch_"+str(bestid+1)+"_"+str(
                                                        round(his_loss[bestid],2))+".pth"))
    else:
        engine.model.load_state_dict(torch.load(model_name))
    amae = []
    amape = []
    armse = []

    if args.data == 'syn':
        if same_G:
            for iter, (x, y) in enumerate(dataloader['test_loader'].get_iterator()):
                testx = torch.Tensor(x).to(device)
                testx = testx.transpose(1, 3)
                testy = torch.Tensor(y).to(device)
                testy = testy.transpose(1, 3)
                # [64, 2, 80, 15]
                metrics = engine.eval_syn(testx, testy, F_t, G)
                amae.append(metrics[0])
                amape.append(metrics[1])
                armse.append(metrics[2])

        else:
            engine.set_state('test')
            reals = []
            pred_Fs = []
            pred_Es = []
            for iter, (x, y, adj_idx) in enumerate(dataloader['test_loader'].get_iterator()):
                testx = torch.Tensor(x).to(device)
                testx = testx.transpose(1, 3)
                testy = torch.Tensor(y).to(device)
                testy = testy.transpose(1, 3)
                # [64, 2, 80, 15]
                metrics = engine.eval_syn(testx, testy, F_t, G['val'], adj_idx)
                amae.append(metrics[0])
                amape.append(metrics[1])
                armse.append(metrics[2])

                reals.append(testy)
                pred_Fs.append(metrics[3])
                pred_Es.append(metrics[4])

            reals = torch.stack(reals).cpu().numpy()
            reals = reals.reshape(-1, *reals.shape[2:])
            pred_Fs = torch.stack(pred_Fs).cpu().numpy()
            pred_Fs = pred_Fs.reshape(-1, *pred_Fs.shape[2:]).squeeze()
            pred_Es = torch.stack(pred_Es).cpu().numpy()
            pred_Es = pred_Es.reshape(-1, *pred_Es.shape[2:]).squeeze()
            # reals shape: (1984, 2, 80, 15); pred_F/Es shape:(1984, 80, 15)

            # reverse slideing window --> results: (num_nodes, total_timesteps)
            ret = util.reverse_sliding_window([reals[:, 0, :, :].squeeze(),
                                         reals[:, 1, :, :].squeeze(),
                                         pred_Fs, pred_Es])

            viz_node_idx = 0
            plt.figure()
            plt.plot(ret[0][viz_node_idx, :], label='real F')
            plt.plot(ret[1][viz_node_idx, :], label='real E')
            plt.plot(ret[2][viz_node_idx, :], label='pred F')
            plt.plot(ret[3][viz_node_idx, :], label='pred E')
            plt.legend()
            plt.show()

        if model_name is None:
            log = 'On average over seq_length horizons, Test MAE: {:.4f}, Test MAPE: {:.4f}, Test RMSE: {:.4f}'
            print(log.format(np.mean(amae),np.mean(amape),np.mean(armse)))
            torch.save(engine.model.state_dict(), args.save+"_exp"+str(args.expid)+"_best_"+str(round(his_loss[bestid],2))+".pth")

    else:
        outputs = []
        realy = torch.Tensor(dataloader['y_test']).to(device)
        realy = realy.transpose(1,3)[:,0,:,:]

        for iter, (x, y) in enumerate(dataloader['test_loader'].get_iterator()):
            testx = torch.Tensor(x).to(device)
            testx = testx.transpose(1,3)
            with torch.no_grad():
                preds = engine.model(testx).transpose(1,3)
            outputs.append(preds.squeeze())

        yhat = torch.cat(outputs,dim=0)
        yhat = yhat[:realy.size(0),...]

        for i in range(args.seq_length):
            pred = scaler.inverse_transform(yhat[:,:,i])
            real = realy[:,:,i]
            metrics = util.metric(pred,real)
            log = 'Evaluate best model on test data for horizon {:d}, Test MAE: {:.4f}, Test MAPE: {:.4f}, Test RMSE: {:.4f}'
            print(log.format(i+1, metrics[0], metrics[1], metrics[2]))
            amae.append(metrics[0])
            amape.append(metrics[1])
            armse.append(metrics[2])

        log = 'On average over seq_length horizons, Test MAE: {:.4f}, Test MAPE: {:.4f}, Test RMSE: {:.4f}'
        print(log.format(np.mean(amae),np.mean(amape),np.mean(armse)))
        torch.save(engine.model.state_dict(), args.save+"_exp"+str(args.expid)+"_best_"+str(round(his_loss[bestid],2))+".pth")

if __name__ == "__main__":
    t1 = time.time()
    main('garage/syn_epoch_36_0.07.pth')
    # main()
    t2 = time.time()
    print("Total time spent: {:.4f}".format(t2-t1))